A Julia package for fast wavelet transforms (1D, 2D, 3D, by filtering or lifting). The package includes discrete wavelet transforms, columnwise discrete wavelet transforms, and wavelet packet transforms.
1st generation wavelets using filter banks (periodic and orthogonal). Filters are included for the following types: Haar, Daubechies, Coiflet, Symmlet, BattleLemarie, Beylkin, Vaidyanathan.
2nd generation wavelets by lifting (periodic and general type including orthogonal and biorthogonal). Included lifting schemes are currently only for Haar and Daubechies (under development). A new lifting scheme can be easily constructed by users. The current implementation of the lifting transforms is 2x faster than the filter transforms.
Thresholding, best basis and denoising functions, e.g. TI denoising by cycle spinning, best basis for WPT, noise estimation, matching pursuit. See example code and image below.
Wavelet utilities e.g. indexing and size calculation, scaling and wavelet functions computation, test functions, up and down sampling, filter mirrors, coefficient counting, inplace circshifts, and more.
Plotting/visualization utilities for 1D and 2D signals.
See license (MIT) in LICENSE.md.
Install via the package manager and load with using
julia> Pkg.add("Wavelets")
julia> using Wavelets
The functions dwt
and wpt
(and their inverses) are linear operators. See wavelet
below for construction of the type wt
.
Discrete Wavelet Transform
# 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
# 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))
The function wavelet
is a type contructor for the transform functions. The transform type t
can be either WT.Filter
or WT.Lifting
.
wavelet(c, t=WT.Filter, boundary=WT.Periodic)
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.
Class Type  Namebase  Supertype  Numbers 

Haar 
haar  OrthoWaveletClass 

Coiflet 
coif  OrthoWaveletClass 
2:2:8 
Daubechies 
db  OrthoWaveletClass 
1:Inf 
Symlet 
sym  OrthoWaveletClass 
4:10 
Battle 
batt  OrthoWaveletClass 
2:2:6 
Beylkin 
beyl  OrthoWaveletClass 

Vaidyanathan 
vaid  OrthoWaveletClass 

CDF 
cdf  BiOrthoWaveletClass 
(9,7) 
Class information
WT.class(::WaveletClass) ::String # class full name
WT.name(::WaveletClass) ::String # type short name
WT.vanishingmoments(::WaveletClass) # vanishing moments of wavelet function
Transform type information
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)
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)
wt = wavelet(WT.Coiflet{4}(), WT.Filter, WT.Periodic)
For a periodic biorthogonal CDF 9/7 lifting scheme:
wt = wavelet(WT.cdf97, WT.Lifting)
Perform a transform of vector x
# 5 level transform
xt = dwt(x, wt, 5)
# inverse tranform
xti = idwt(xt, wt, 5)
# a full transform
xt = dwt(x, wt)
Other examples:
# scaling filters is easy
wt = wavelet(WT.haar)
wt = WT.scale(wt, 1/sqrt(2))
# signals can be transformed inplace with a lowlevel command
# requiring very little memory allocation (especially for L=1 for filters)
dwt!(x, wt, L) # inplace (lifting)
dwt!(xt, x, wt, L) # write to xt (filter)
# denoising with default parameters (VisuShrink hard thresholding)
x0 = testfunction(128, "HeaviSine")
x = x0 + 0.3*randn(128)
y = denoise(x)
# plotting utilities 1d (see images and code in /example)
x = testfunction(128, "Bumps")
y = dwt(x, wavelet(WT.cdf97, WT.Lifting))
d,l = wplotdots(y, 0.1, 128)
A = wplotim(y)
# plotting utilities 2d
img = imread("lena.png")
x = permutedims(img.data, [ndims(img.data):1:1])
L = 2
xts = wplotim(x, L, wavelet(WT.db3))
See Bumps and Lena for plot images.
The Wavelets.Threshold
module includes the following utilities
# threshold types with parameter
threshold!(x::AbstractArray, TH::THType, t::Real)
threshold(x::AbstractArray, TH::THType, t::Real)
# without parameter (PosTH, NegTH)
threshold!(x::AbstractArray, TH::THType)
threshold(x::AbstractArray, TH::THType)
# denoising
denoise(x::AbstractArray,
wt=DEFAULT_WAVELET;
L::Int=min(maxtransformlevels(x),6),
dnt=VisuShrink(size(x,1)),
estnoise::Function=noisest,
TI=false,
nspin=tuple([8 for i=1:ndims(x)]...) )
# entropy
coefentropy(x, et::Entropy, nrm)
# best basis for WPT limited to active inital tree nodes
bestbasistree(y::AbstractVector, wt::DiscreteWavelet,
L::Integer=nscales(y), et::Entropy=ShannonEntropy() )
bestbasistree(y::AbstractVector, wt::DiscreteWavelet,
tree::BitVector, et::Entropy=ShannonEntropy() )
Type  Details 

Thresholding  <: THType 
HardTH 
hard thresholding 
SoftTH 
soft threshold 
SemiSoftTH 
semisoft thresholding 
SteinTH 
stein thresholding 
PosTH 
positive thresholding 
NegTH 
negative thresholding 
BiggestTH 
biggest mterm (best mterm) approx. 
Denoising  <: DNFT 
VisuShrink 
VisuShrink denoising 
Entropy  <: Entropy 
ShannonEntropy 
v^2*log(v^2) (CoifmanWickerhauser) 
LogEnergyEntropy 
log(v^2) 
Find best basis tree for wpt
, and compare to dwt
using biggest mterm approximations.
wt = wavelet(WT.db4)
x = sin.(4*range(0, stop=2*pieps(), length=1024))
tree = bestbasistree(x, wt)
xtb = wpt(x, wt, tree)
xt = dwt(x, wt)
# get biggest mterm approximations
m = 50
threshold!(xtb, BiggestTH(), m)
threshold!(xt, BiggestTH(), m)
# compare sparse approximations in ell_2 norm
norm(x  iwpt(xtb, wt, tree), 2) # best basis wpt
norm(x  idwt(xt, wt), 2) # regular dwt
julia> norm(x  iwpt(xtb, wt, tree), 2)
0.008941070750964843
julia> norm(x  idwt(xt, wt), 2)
0.05964431178940861
``` html
Denoising:
n = 2^11; x0 = testfunction(n,"Doppler") x = x0 + 0.05*randn(n) y = denoise(x,TI=true)
![Doppler](/example/denoise_doppler.png)
Benchmarks

Timing of `dwt` (using `db2` filter of length 4) and `fft`. The lifting wavelet transforms are faster and use less memory than `fft` in 1D and 2D. `dwt` by lifting is currently 2x faster than by filtering.
dwt by filtering (N=1048576), 20 levels elapsed time: 0.247907616 seconds (125861504 bytes allocated, 8.81% gc time) dwt by lifting (N=1048576), 20 levels elapsed time: 0.131240966 seconds (104898544 bytes allocated, 17.48% gc time) fft (N=1048576), (FFTW) elapsed time: 0.487693289 seconds (167805296 bytes allocated, 8.67% gc time)
For 2D transforms (using a `db4` filter and CDF 9/7 lifting scheme):
dwt by filtering (N=1024x1024), 10 levels elapsed time: 0.773281141 seconds (85813504 bytes allocated, 2.87% gc time) dwt by lifting (N=1024x1024), 10 levels elapsed time: 0.317705928 seconds (88765424 bytes allocated, 3.44% gc time) fft (N=1024x1024), (FFTW) elapsed time: 0.577537263 seconds (167805888 bytes allocated, 5.53% gc time)
By using the lowlevel function `dwt!` and preallocating temporary arrays, significant gains can be made in terms of memory usage (and some speedup). This is useful when transforming multiple times.
Todo list

* Transforms for nonsquare 2D signals
* Boundary extensions (other than periodic)
* Boundary orthogonal wavelets
* Define more lifting schemes
* WPT in 2D
* Stationary transform
* Continuous wavelets
* Wavelet scalogram
09/26/2014
11 days ago
211 commits