Merge pull request #25 from JuliaDSP/gfa/types

Refactor wavelet types and change API

.gitignore

@@ -27,4 +27,5 @@
 /usr
 
 lena*
+testmod.jl
 

README.md

@@ -1,7 +1,5 @@
 <img src="wavelets.png" alt="Wavelets">
 
----------
-
 [![Build Status](https://travis-ci.org/JuliaDSP/Wavelets.jl.svg?branch=master)](https://travis-ci.org/JuliaDSP/Wavelets.jl)
 [![Coverage Status](https://coveralls.io/repos/JuliaDSP/Wavelets.jl/badge.svg?branch=master)](https://coveralls.io/r/JuliaDSP/Wavelets.jl?branch=master)
 
@@ -17,8 +15,6 @@ A [Julia](https://github.com/JuliaLang/julia) package for fast wavelet transform
 
 * Plotting/visualization utilities for 1-D and 2-D signals.
 
-Loosely inspired by [this](https://github.com/tomaskrehlik/Wavelets) and [this](http://statweb.stanford.edu/~wavelab). 
-
 See license (MIT) in LICENSE.md.
 
 
@@ -36,33 +32,38 @@ julia> using Wavelets
 API
 ---------
 
-#### Wavelet transforms and types
+#### Wavelet Transforms
+See `wavelet` below for construction of the type `wt`.
+
+**Discrete Wavelet Transform**
+```julia
+# DWT
+dwt(x, wt, L=maxtransformlevels(x))
+idwt(x, wt, L=maxtransformlevels(x))
+dwt!(y, x, filter, L=maxtransformlevels(x))
+idwt!(y, scheme, L=maxtransformlevels(x))
+```
+
+**Wavelet Packet Transform**
+```julia
+# WPT (tree can also be an integer, equivalent to maketree(length(x), L, :full))
+wpt(x, wt, tree::BitVector=maketree(x, :full))
+iwpt(x, wt, tree::BitVector=maketree(x, :full))
+wpt!(y, x, filter, tree::BitVector=maketree(x, :full))
+iwpt!(y, scheme, tree::BitVector=maketree(y, :full))
+```
+
+#### Wavelet Types
+The function `wavelet` is a type contructor for the transform functions. The transform type `t` can be either `WT.Filter` or `WT.Lifting`.
+
 ```julia
-# Transform Type construction
-wavelet(w::WT.WaveletClass, boundary::WaveletBoundary=Periodic)  # defaults to filter
-waveletfilter(w::WT.WaveletClass, boundary::WaveletBoundary=Periodic)
-waveletls(w::WT.WaveletClass, boundary::WaveletBoundary=Periodic)
-
-# DWT (discrete wavelet transform)
-dwt(x::AbstractArray, wt::DiscreteWavelet, L::Integer=maxtransformlevels(x))
-idwt(x::AbstractArray, wt::DiscreteWavelet, L::Integer=maxtransformlevels(x))
-dwt!(y::AbstractArray, x::AbstractArray, filter::OrthoFilter, L::Integer, fw::Bool)
-dwt!(y::AbstractArray, scheme::GLS, L::Integer, fw::Bool)
-# DWTC (discrete wavelet transform along dimension td (default is last dim.))
-dwtc(x::AbstractArray, wt::DiscreteWavelet, L::Integer=maxtransformlevels(x), td::Integer=ndims(x))
-idwtc(x::AbstractArray, wt::DiscreteWavelet, L::Integer=maxtransformlevels(x), td::Integer=ndims(x))
-# WPT (wavelet packet transform)
-# Ltree can be L::Integer=maxtransformlevels(x) or tree::BitVector=maketree(length(y), L, :full)
-wpt(x::AbstractArray, wt::DiscreteWavelet, Ltree)
-iwpt(x::AbstractArray, wt::DiscreteWavelet, Ltree)
-wpt!(y::AbstractArray, x::AbstractArray, filter::OrthoFilter, Ltree, fw::Bool)
-wpt!(y::AbstractArray, scheme::GLS, Ltree, fw::Bool)
+wavelet(c, t=WT.Filter, boundary=WT.Periodic)
 ```
 
-#### Wavelet classes
+#### Wavelet Classes
 
 The module WT contains the types for wavelet classes. The module defines constants of the form e.g. `WT.coif4` as shortcuts for `WT.Coiflet{4}()`.
-The numbers for orthogonal wavelets indicate the number vanishing moments of the wavelet function. The length of a `wt::OrthoFilter` can be obtained with `length(wt)`.
+The numbers for orthogonal wavelets indicate the number vanishing moments of the wavelet function.
 
 | Class Type | Namebase | Supertype | Numbers |
 |:------- |:------ |:----- |:----- |
@@ -77,65 +78,72 @@ The numbers for orthogonal wavelets indicate the number vanishing moments of the
 
 Class information
 ```julia
-class(::WaveletClass) ::ASCIIString         # class full name
-name(::WaveletClass) ::ASCIIString          # type short name
-vanishingmoments(::WaveletClass)            # vanishing moments of wavelet function
+WT.class(::WaveletClass) ::ASCIIString         # class full name
+WT.name(::WaveletClass) ::ASCIIString          # type short name
+WT.vanishingmoments(::WaveletClass)            # vanishing moments of wavelet function
 ```
 Transform type information
 ```julia
-length(f::OrthoFilter)						# length of filter
-qmf(f::OrthoFilter)							# quadrature mirror filter
-name(f::OrthoFilter)
-makeqmfpair(f::OrthoFilter, fw::Bool=true, T::Type=eltype(qmf(f)))
-makereverseqmfpair(f::OrthoFilter, fw::Bool=true, T::Type=eltype(qmf(f)))
-name(s::GLS)
+WT.name(wt)                                     # type short name
+WT.length(f::OrthoFilter)                       # length of filter
+WT.qmf(f::OrthoFilter)                          # quadrature mirror filter
+WT.makeqmfpair(f::OrthoFilter, fw=true)
+WT.makereverseqmfpair(f::OrthoFilter, fw=true)
 ```
 
 Examples
 ---------
 
+The simplest way to transform a signal x is
 ```julia
-# the simplest way to transform a signal x is
 xt = dwt(x, wavelet(WT.db2))
+```
+
+The transform type can be more explicitly specified (filter, Periodic, Orthogonal, 4 vanishing moments)
+```julia
+wt = wavelet(WT.Coiflet{4}(), WT.Filter, WT.Periodic)
+```
+
+For a periodic biorthogonal CDF 9/7 lifting scheme:
+```julia
+wt = wavelet(WT.cdf97, WT.Lifting)
+```
 
-# the transform type can be more explicitly specified
-# set up wavelet type (filter, Periodic, Orthogonal, 4 vanishing moments)
-wt = wavelet(WT.Coiflet{4}(), Periodic)  # or
-wt = waveletfilter(WT.coif4)
-# which is equivalent to 
-wt = OrthoFilter(WT.coif4)
-# the object wt determines the transform type 
-# wt now contains instructions for a periodic biorthogonal CDF 9/7 lifting scheme
-wt = waveletls(WT.cdf97)
-# xt is a 5 level transform of vector x
+Perform a transform of vector x
+```julia
+# 5 level transform
 xt = dwt(x, wt, 5)
 # inverse tranform
 xti = idwt(xt, wt, 5)
 # a full transform
 xt = dwt(x, wt)
+```
 
+Other examples:
+```julia
 # scaling filters is easy
-wt = scale(wt, 1/sqrt(2))
+wt = wavelet(WT.haar)
+wt = WT.scale(wt, 1/sqrt(2))
 # signals can be transformed inplace with a low-level command
 # requiring very little memory allocation (especially for L=1 for filters)
-dwt!(x, wt, L, true)      # inplace (lifting)
-dwt!(xt, x, wt, L, true)  # write to xt (filter)
+dwt!(x, wt, L)      # inplace (lifting)
+dwt!(xt, x, wt, L)  # write to xt (filter)
 
 # denoising with default parameters (VisuShrink hard thresholding)
-x0 = testfunction(n, "HeaviSine")
-x = x0 + 0.3*randn(n)
+x0 = testfunction(128, "HeaviSine")
+x = x0 + 0.3*randn(128)
 y = denoise(x)
 
 # plotting utilities 1-d (see images and code in /example)
-x = testfunction(n, "Bumps")
-y = dwt(x, waveletls(WT.cdf97))
-d,l = wplotdots(y, 0.1, n)
+x = testfunction(128, "Bumps")
+y = dwt(x, wavelet(WT.cdf97, WT.Lifting))
+d,l = wplotdots(y, 0.1, 128)
 A = wplotim(y)
 # plotting utilities 2-d
 img = imread("lena.png")
 x = permutedims(img.data, [ndims(img.data):-1:1])
 L = 2
-xts = wplotim(x, L, waveletfilter(WT.db3))
+xts = wplotim(x, L, wavelet(WT.db3))
 ```
 
 See [Bumps](/example/transform1d_bumps.png) and [Lena](/example/transform2d_lena.jpg) for plot images.
@@ -209,7 +217,7 @@ threshold!(xtb, BiggestTH(), m)
 threshold!(xt, BiggestTH(), m)
 # compare sparse approximations in ell_2 norm
 vecnorm(x - iwpt(xtb, wt, tree), 2) # best basis wpt
-vecnorm(x - idwt(xt, wt), 2)   		# regular dwt
+vecnorm(x - idwt(xt, wt), 2)        # regular dwt
 ```
 ```
 julia> vecnorm(x - iwpt(xtb, wt, tree), 2)

REQUIRE

@@ -1,3 +1,4 @@
+julia 0.3.7
 Reexport
-julia 0.3-
 Compat
+Docile 0.4.11

example/transform1d.jl

@@ -5,7 +5,7 @@ using PyPlot
 J = 11
 n = 2^J
 x = testfunction(n,"Bumps")
-y = dwt(x, waveletls(WT.cdf97))
+y = dwt(x, wavelet(WT.cdf97, WT.Filter))
 d,l = wplotdots(y, 0.1, n)
 A = wplotim(y)
 

example/transform2d.jl

@@ -5,7 +5,7 @@ using Images  # for imread and imwrite
 img = imread("lena.tiff")
 x = permutedims(img.data, [ndims(img.data):-1:1])
 L = 2
-xts = wplotim(x, L, waveletfilter(WT.db3))
+xts = wplotim(x, L, wavelet(WT.db3, WT.Filter))
 
 imwrite(xts, "lena_2d.png")
 # convert with ImageMagick

perf/bm_dwt_filt_inp.jl

@@ -0,0 +1,10 @@
+
+include("setup_1d.jl")
+
+y = copy(x0)
+
+println("dwt! by filtering (N=",N,"), ", L, " levels")
+f(y, x0, L, wf) = for i = 1:tn; dwt!(y, x0, wf, L); end
+f(y, x0, L, wf);
+@time f(y, x0, L, wf);
+

perf/bm_dwt_ls_inp.jl

@@ -0,0 +1,8 @@
+
+include("setup_1d.jl")
+
+println("dwt! by lifting (N=",N,"), ", L, " levels")
+f(x0, L, wl) = for i = 1:tn; dwt!(x0, wl, L); end
+f(x0, L, wl);
+@time f(x0, L, wl);
+

perf/bm_fft_inp.jl

@@ -0,0 +1,10 @@
+
+include("setup_1d.jl")
+
+x0c = complex(x0)
+
+println("fft! (N=",N,"), (FFTW)")
+f(x0c) = for i = 1:tn; fft!(x0c); end
+f(x0c);
+@time f(x0c);
+

perf/results_1d.txt

@@ -1,7 +1,7 @@
 
 dwt by filtering (N=1048576), 20 levels
-elapsed time: 0.247907616 seconds (125861504 bytes allocated, 8.81% gc time)
+elapsed time: 0.287814759 seconds (125861424 bytes allocated, 14.58% gc time)
 dwt by lifting (N=1048576), 20 levels
-elapsed time: 0.131240966 seconds (104898544 bytes allocated, 17.48% gc time)
+elapsed time: 0.1401407 seconds (104888464 bytes allocated, 21.04% gc time)
 fft (N=1048576), (FFTW)
-elapsed time: 0.487693289 seconds (167805296 bytes allocated, 8.67% gc time)
+elapsed time: 0.428999777 seconds (167805536 bytes allocated, 12.92% gc time)

perf/results_1d_inp.txt

@@ -0,0 +1,7 @@
+
+dwt! by filtering (N=1048576), 20 levels
+elapsed time: 0.2781087 seconds (41974064 bytes allocated, 6.29% gc time)
+dwt! by lifting (N=1048576), 20 levels
+elapsed time: 0.095285706 seconds (21001104 bytes allocated)
+fft! (N=1048576), (FFTW)
+elapsed time: 0.324933232 seconds (32704 bytes allocated)

perf/results_2d.txt

@@ -1,7 +1,7 @@
 
 dwt by filtering (N=1024x1024), 10 levels
-elapsed time: 0.773281141 seconds (85813504 bytes allocated, 2.87% gc time)
+elapsed time: 0.816053876 seconds (85813424 bytes allocated, 1.77% gc time)
 dwt by lifting (N=1024x1024), 10 levels
-elapsed time: 0.317705928 seconds (88765424 bytes allocated, 3.44% gc time)
+elapsed time: 0.41193094 seconds (89424144 bytes allocated, 3.40% gc time)
 fft (N=1024x1024), (FFTW)
-elapsed time: 0.577537263 seconds (167805888 bytes allocated, 5.53% gc time)
+elapsed time: 0.583925703 seconds (167806176 bytes allocated, 9.77% gc time)

perf/run1d_inp.sh

@@ -0,0 +1,20 @@
+#!/bin/bash
+
+usage="Usage: `basename $0` "
+if [ $# -ne 0 ]; then 		# variable supplied?
+	echo $usage 1>&2
+	exit 1 
+fi
+
+results="results_1d_inp.txt"
+echo "" >"$results"
+# warm up
+julia bm_dwt_filt_inp.jl > /dev/null
+
+julia bm_dwt_filt_inp.jl >>"$results"
+julia bm_dwt_ls_inp.jl >>"$results"
+julia bm_fft_inp.jl >>"$results"
+
+exit
+
+

perf/setup_1d.jl

@@ -1,10 +1,10 @@
 
 using Wavelets
 
-wf = waveletfilter(WT.db2)
-wl = waveletls(WT.db2)
+wf = wavelet(WT.db2, WT.Filter)
+wl = wavelet(WT.db2, WT.Lifting)
 N = 1024*1024;
 x0 = rand(N);
-L = nscales(N)  # int(log2(N)-2)
+L = maxtransformlevels(N)  # int(log2(N)-2)
 tn = 10
 

perf/setup_2d.jl

@@ -1,10 +1,10 @@
 
 using Wavelets
 
-wf = waveletfilter(WT.db4)
-wl = waveletls(WT.cdf97)
+wf = wavelet(WT.db4, WT.Filter)
+wl = wavelet(WT.cdf97, WT.Lifting)
 N = 1024;
 x0 = rand(N,N);
-L = nscales(N)  # int(log2(N)-2)
+L = maxtransformlevels(N)  # int(log2(N)-2)
 tn = 10
 

src/Wavelets.jl

@@ -2,14 +2,14 @@ module Wavelets
 
 include("util.jl")
 
-include("wavelettypes.jl")
+include("wt.jl")
 include("transforms.jl")
 
 include("threshold.jl")
 include("plot.jl")
 
 using Reexport
-@reexport using .Util, .WaveletTypes, .Transforms, .Threshold, .Plot
+@reexport using .Util, .WT, .Transforms, .Threshold, .Plot
 
 end
 

src/plot.jl

@@ -1,5 +1,5 @@
 module Plot
-using ..Util, ..WaveletTypes, ..Transforms
+using ..Util, ..WT, ..Transforms
 export wplotdots, wplotim
 
 # PLOTTING UTILITIES