Wavelet Analysis: Applications

Wolfram Research
wolfram-u@wolfram.com

Introduction: Wavelet Analysis

Fourier transforms are applied to signals to obtain frequency information that is not readily available in the time domain.

Frequency is measured in cycles per second or Hertz (Hz)

Rapid change in events is high frequency

Slow changing events are low frequency

The Fourier transform of a signal f(t) is a function ℱ(ω) in the frequency domain. The time variable t is discrete, the frequency variable ω is continuous:

ℱ(ω)=(∑_(t=-∞)^∞)f(t) (e^(-i t ω))

The inverse Fourier transform recovers f(t) :

f(t) = (1/(2π))(∫_(-π)^π)ℱ(ω) (e^(i t ω)) dω

The signal is divided into frequencies split into uniform bands from ω to ω+Δω providing a very good understanding of the frequency information embedded within the signal.

Consider for example a stationary signal (a signal is stationary when frequency information does not change in time).

ToPictureCell[Labeled[Block[{data,st},st[t_]:=Sin[2 π 20 t]+Sin[2 π 40 t];data=Table[st[t],{t,-π,π,((2 π)/2047)}];Row[{Plot[st[t],{t,0,6},ImageSize->450,PlotLabel->Style["Stationary Data",FontFamily->"Verdana",Darker[Brown],Bold,16]],ListLinePlot[Take[fourier[data,((2 π)/2047)],500],ImageSize->450,PlotLabel->Style["Fourier Transform",FontFamily->"Verdana",Darker[Brown],Bold,16],AxesLabel -> {Style["2 frequency",Bold,18], ""}]}]],
Style[TraditionalForm[Sin[2 π 20 𝓉]+Sin[2 π 40 𝓉]],Directive[sty,18]],Top]]

In the above example, the Fourier transform provides good frequency and amplitude information. What it lacks is the time information; i.e., at what time did these frequencies occur? In real life the signals are not stationary and since the Fourier transform models the signal as a periodic signal, when you view a non-stationary signal in the Fourier domain it starts to look distorted.

ToPictureCell[Labeled[Block[{data,nst},
nst[t_]:=({{\[Piecewise], ({{Cos[2 π t], 0≤t<2}, {Cos[10 π t], 2≤t<4}, {Cos[20 π t], 4≤t<8}, {0, True}})}});
data=N[Table[nst[t],{t,0,6,(6/2047)}]];Row[{Plot[nst[t],{t,0,6},ImageSize->450,PlotLabel->Style["Non Stationary Data",FontFamily->"Verdana",Darker[Brown],Bold,16]],ListLinePlot[Take[fourier[data,(6/2047)],100],ImageSize->450,PlotLabel->Style["Fourier Transform",FontFamily->"Verdana",Darker[Brown],Bold,16],AxesLabel -> {Style["2 frequency",Bold,18], ""}]}]],Style[TraditionalForm[({{\[Piecewise], ({{Cos[2 π 𝓉], 0≤𝓉<2}, {Cos[10 π 𝓉], 2≤𝓉<4}, {Cos[20 π 𝓉], 4≤𝓉<8}, {0, True}})}})],Directive[sty,18]],Top]]

In order to overcome the limitation of traditional Fourier transforms we break down the signal into smaller components.

Generally it is safe to say that a signal is stationary over a finite interval.

This adds an extra layer of complexity in choosing an appropriate window size.

ToPictureCell@Labeled[Block[{data,nst,sft,rect},
nst[t_]:=({{\[Piecewise], ({{Cos[2 π t], 0≤t<2}, {Cos[10 π t], 2≤t<4}, {Cos[20 π t], 4≤t<8}, {0, True}})}});
sft[t_,n_]:= Switch[n,1,({{\[Piecewise], ({{1, 0≤t<2}, {0, True}})}}),2,({{\[Piecewise], ({{1, 2≤t<4}, {0, True}})}}),3,({{\[Piecewise], ({{3, 4≤t<6}, {0, True}})}})];
rect[n_]:= {FaceForm[None],EdgeForm[{Thick,Darker[Green]}],Switch[n,1,Rectangle[{0,-1},{2,1}],2,Rectangle[{2,-1},{4,1}],3,Rectangle[{4,-3},{6,3}]]};
data[n_]:=Select[N@Table[sft[t,n]nst[t],{t,0, 6,(6/2047)}],#≠0&];
Grid[Join[{{Style["Data",FontFamily->"Verdana",Brown,Bold,16],Style["Fourier Transform",FontFamily->"Verdana",Brown,Bold,16]}}
,
Table[{
Plot[sft[t,n]nst[t],{t,0,6},ImageSize->400,AspectRatio->0.4,ImagePadding->{{30,Automatic},{Automatic,Automatic}},Epilog->rect[n]]
,
ListLinePlot[Take[fourier[data[n],(6/2047)],100],ImageSize->400,AspectRatio->0.4,ImagePadding->{{30,Automatic},{Automatic,Automatic}}]},{n,3}]]]],Style[TraditionalForm[({{\[Piecewise], ({{Cos[2 π t], 0≤t<2}, {Cos[10 π t], 2≤t<4}, {Cos[20 π t], 4≤t<8}, {0, True}})}})],Directive[sty,18]],Right]

Introduction: Wavelet Analysis

Consider a case where one encounters a significant short-lived event in a rather benign time series.

data=Table[Sin[4 π t]+2 Exp[-(10^5) ((((1/3)-t))^2)],{t,0,1,(1/511)}];
p1 = ListLinePlot[data]

Performing a Fourier analysis on the above signal, one realizes that the short-lived anomaly is completely missed. The Fourier analysis captures only the sine curve and is unable to resolve the short-lived event. Furthermore, even if it captures the event, it cannot tell us when it happened.

ListLinePlot[fourier[data,1/512], PlotRange -> All]

In real-world applications, such phenomena are often of critical importance (eg: earthquake, heart beats, etc.). Such time series can be efficiently analyzed using wavelet analysis.

cwt=ContinuousWaveletTransform[data,PaulWavelet[5],{7,12},SampleRate->512];
p2 = WaveletScalogram[cwt,ColorFunction->"AvocadoColors",Ticks -> Full,AxesLabel -> {"Time", "{Oct.,Voices}"}];
GraphicsRow[{p1,p2}]

Wavelets decompose a signal into approximations and details at different scales. This provides a very good understanding of the events and their location in the signal. Using a combination of scaling functions (ϕ) and wavelet functions (ψ) a large class of signals can be represented by

f(t)=(∑_(k ϵ ℤ))(Subscript[c, 0,k])(Subscript[ϕ, 0,k])(t)+(∑_(j=0)^∞)(∑_(k ϵ ℤ))(Subscript[d, j,k])(Subscript[ψ, j,k])(t),

where the coefficients are given by

(Subscript[c, j,k])(f)=⟨f,(Subscript[ϕ, j,k])⟩ = ∫f(t)(Subscript[ϕ, j,k])dt 
(Subscript[d, j,k])(f)=⟨f,(Subscript[ψ, j,k])⟩ = ∫f(t)(Subscript[ψ, j,k])dt

The signal is divided into different scales (j) of resolution, rather than frequencies. This form of analysis is called  multiresolution analysis and is fundamental to wavelet analysis.

The “time-frequency plane” in Fourier analysis is now taken by “time-scale plane” in wavelets.

Multiresolution analysis divides the frequencies into octave bands, from ω to 2ω (1, 2, 4, 8, 16, …) in contrast to Fourier analysis where frequencies split into uniform bands from ω to ω+Δω.

Frequency Detection in a Time Series

A single frequency shows up as a horizontal band at the equivalent scale and multiple frequencies show up as multiple bands at the equivalent scales.

Single Frequency

xval=N[Range[0,1,(1/1023)]];
generateData[fr_]:=Cos[2 π fr xval];
cwtfunc[fr_]:=ContinuousWaveletTransform[generateData[fr],GaborWavelet[6],{8,4},SampleRate->1023]
Manipulate[Column[{Grid[{{"Frequency (Hz):",fr}}],ListLinePlot[generateData[fr],AspectRatio->0.2,Axes->False],WaveletScalogram[cwtfunc[fr],Frame->True,PlotRangePadding->0.0001]}],{{fr,10,"Frequency"},10,200,5}]

Two Frequencies

generateData[fr1_,fr2_]:=N[Cos[2 π fr1 xval]+Cos[2 π fr2 xval]]
cwtfunc[fr1_,fr2_]:=ContinuousWaveletTransform[generateData[fr1,fr2],GaborWavelet[6],{8,4},SampleRate->1023]
Manipulate[Column[{Grid[{{"Frequency 1 (Hz):",fr1},{"Frequency 2 (Hz):",fr2}}],ListLinePlot[generateData[fr1,fr2],AspectRatio->0.2,Axes->False],WaveletScalogram[cwtfunc[fr1,fr2],Frame->True,PlotRangePadding->0.0001]}],{{fr1,10,"Frequency"},10,200,5},{{fr2,40,"Frequency"},10,200,5}]

Clear[cwtfunc,generateData,xval];

Higher frequencies are resolved at lower octaves and lower frequencies at higher octaves.

data=Table[({{\[Piecewise], ({{Sin[20 π t], 0≤t<(1/4)}, {Sin[50 π t], (1/4)≤t<(1/2)}, {Sin[100 π t], (1/2)≤t<(3/4)}, {Sin[200 π t], (3/4)≤t≤1}, {0, True}})}}),{t,0,1,(1/1023)}];

ListLinePlot[data,AspectRatio->0.2]

cwd=ContinuousWaveletTransform[data,DGaussianWavelet[5],{Automatic,12}]

WaveletScalogram[cwd,ColorFunction->"RustTones",Ticks->{Automatic,({#1,(Superscript[2, #1])}&)/@Range[9]}]

Analyzing a Time Series Using Wavelets

Wavelet analysis can be used for detecting discontinuities, isolating peaks, and non-stationary behavior.

Real wavelets can be used to locate discontinuities.

data=Table[HeavisideTheta[x],{x,-2,2,(4/1023)}];

ListLinePlot[data]

cwt=ContinuousWaveletTransform[data,DGaussianWavelet[4],{Automatic,12},Padding->"Fixed"]

WaveletScalogram[cwt]

Complex wavelets can be used to capture oscillatory behavior.

f[t_,ϵ_]:=(Sqrt[t (1-t)]) Sin[((2 π (1+ϵ))/(t+ϵ))]

data=Table[f[t,0.05],{t,0,1,(1/1023)}];

ListLinePlot[data]

cwt=ContinuousWaveletTransform[data,GaborWavelet[6],{Automatic,16}]

Amplitude of wavelet coefficients:

WaveletScalogram[cwt,ColorFunction->"BlueGreenYellow"]

Phase of wavelet coefficients:

WaveletScalogram[cwt,Automatic,Arg,ColorFunction->"BlueGreenYellow"]

ContinuousWaveletTransform can be used to detect certain features in a time series.

ClearAll[data,cwt]

data=Table[Sin[4 π t]+2 Exp[-(10^5) ((((1/3)-t))^2)],{t,0,1,(1/1023)}];

ListLinePlot[data]

cwt=ContinuousWaveletTransform[data,PaulWavelet[5],{7,12},SampleRate->1023]

WaveletScalogram[cwt,ColorFunction->"AvocadoColors",ImageSize->Large]

Filtering Frequencies from Time Series

Wavelet transforms can be used to filter frequencies.

d1 = Table[Sin[x]+Cos[2x],{x,-2π,2π,((4π)/2047)}];

d2 = Table[((ArcSin[Sin[20x]])/5),{x,-2π,2π,((4π)/2047)}];

ListLinePlot[d1+d2]

To filter out the two signals, first perform a wavelet transform.

dwd = DiscreteWaveletTransform[d1+d2,SymletWavelet[4],6];

Use WaveletListPlot to visualize frequency distribution.

WaveletListPlot[dwd,PlotLayout->"CommonXAxis",DataRange->{-2π,2π}]

To filter low frequencies, keep only the coarse coefficients.

rdwd1 = WaveletMapIndexed[(#*0.0)&,dwd,{___,1}];

GraphicsRow[{ListLinePlot[d1,PlotLabel->"data"],ListLinePlot[InverseWaveletTransform[rdwd1],PlotLabel->"filtered coarse data"]}, ImageSize -> 900]

To filter high frequencies, keep only the detail coefficients.

rdwd2 = WaveletMapIndexed[(#*0.0) &, dwd, {___, 0}];

GraphicsRow[{ListLinePlot[d2,PlotLabel->"data"],ListLinePlot[InverseWaveletTransform[rdwd2],PlotLabel->"filtered fine data"]}, ImageSize -> 900]

An example demonstrating frequency filtering.

f[x_]:=Sin[50 π x] Exp[-100 π (((x-0.5))^2)]+(Sin[50 π x]+2 Cos[140 π x]) Exp[-50 π (((x-0.2))^2)]+2 Sin[150 π x] Exp[-80 π (((x-0.8))^2)]

data=Table[f[x],{x,0.,1,(1/2047)}];

ListLinePlot[data]

cwt=ContinuousWaveletTransform[data,GaborWavelet[6],{10,16},Padding->0.0,SampleRate->2047,WaveletScale->Automatic]

WaveletScalogram[cwt,ColorFunction->"SolarColors"]

Filter the cosine with frequency 140π.

vec=ConstantArray[1,Length[data]];
vec[[1;;1000]]:=0.0

cwtThresh=WaveletMapIndexed[#1 vec&,cwt,{{4|5,_},{6,u_/;u<8}}];

WaveletScalogram[cwtThresh,ColorFunction->"SolarColors",ImageSize->Large]

rdata=InverseContinuousWaveletTransform[cwtThresh];

ListLinePlot[rdata]

ClearAll[data,rdata,cwtThresh,cwt,vec,f];

Directly analyze a sound using continuous wavelet transforms.

snd = ExampleData[{"Sound", "Apollo11ReturnSafely"}]

cwd=ContinuousWaveletTransform[snd,GaborWavelet[6],{Automatic,12}]

WaveletScalogram[cwd, {3|4 | 5 | 6 | 7, _},Ticks->{Automatic,{#,Superscript[2,#]}&/@Range[7]},AxesLabel->{Style["Time",sty],Style["Scale",sty]},ColorFunction->ColorData["SiennaTones"]]

Applying Wavelet Analysis to Finance

Extract the stock price trend for IBM since January 1, 2000.

price=FinancialData["IBM","Jan. 1, 2000"];

i1=DateListPlot[price]

dwd=DiscreteWaveletTransform[price[[All,2]],BiorthogonalSplineWavelet[3,3],4];

Find Financial Series Trend

The trend of the series is captured in the low-pass filter coefficients.

WaveletListPlot[dwd,PlotLayout->"CommonYAxis"]

Thresholding all detail coefficients and inverting the series gives the trend.

tr=InverseWaveletTransform[WaveletThreshold[dwd,{"Hard",100}]];

Row[{i1,Show[,ImageSize->70],DateListPlot[Transpose[{price[[All,1]],tr}]]}]

Detrend a Financial Series

Detail coefficients captured the detrended series.

WaveletListPlot[dwd]

Remove the trend by removing the coarse coefficients and inverting.

dtr=InverseWaveletTransform[WaveletMapIndexed[#1 0.0&,dwd,{___,0}]];

Row[{i1,Show[,ImageSize->70],DateListPlot[Transpose[{price[[All,1]],dtr}]]}]

Edge Detection on Images

Similar to detecting events and discontinuities in time series, wavelet analysis can be used to detect edges and discontinuities in images and two-dimensional data.

Example 1

dwd = DiscreteWaveletTransform[,SymletWavelet[2],4];

Visualize wavelet decomposition of image using WaveletImagePlot.

WaveletImagePlot[dwd,Automatic,ImageAdjust[ImageAdjust[#],{1,0.2,1.9}]&]

Set coarse coefficients to 0 and reconstruct using detail coefficients only.

imgEdge[img_,{0,0,0,0}]:= ImageApply[# 0.0&,img]
imgEdge[img_,___]:= ImageApply[# 2&,img]

Sharpen[InverseWaveletTransform[WaveletMapIndexed[imgEdge,dwd]]]

Example 2

dwd=LiftingWaveletTransform[,SymletWavelet[2],3]

WaveletImagePlot[dwd,Automatic,Sharpen[ImageAdjust[ImageApply[Abs,#1]]]&]

Set coarse coefficients to 0 and reconstruct using detail coefficients only.

rimg =Image[InverseWaveletTransform[dwd,Automatic,{___,1|2|3}],ImageSize->500];

Row[{Sharpen[rimg],"   ",Binarize[rimg]}]

Two-dimensional wavelet transform enables one to separate out horizontal and vertical discontinuities.

DynamicModule[{angle=0,p={1,0},lwt},Row[{LocatorPane[Dynamic[p,(angle=ArcTan@@#1;p={Cos[angle],Sin[angle]})&],Dynamic[Graphics[{{Thick,Red,Circle[]},{Blue,Line[{-p,p}]}},ImageSize->300]]],Dynamic[WaveletImagePlot[lwt=LiftingWaveletTransform[Image[Graphics[{{Red,Thick,Circle[]},{Blue,Line[{-p,p}]}},ImageSize->300],"Real"],Automatic,3],ImageSize->300,BaseStyle->Orange]]}]]

Images can be pre-processed using wavelets so that the image processing functions can captures features better.

img =;

dwd=LiftingWaveletTransform[img,SymletWavelet[4],3];

Plot transformed image coefficients.

WaveletImagePlot[dwd,Automatic,ImageAdjust[Sharpen[ImageAdjust[ImageApply[Abs,#]]],{0.2,0.8}]&,ImageSize->300]

Extract vertical and diagonal components from the wavelet-transformed image and reconstruct.

verimg=InverseWaveletTransform[dwd,Automatic,{___,1|3}];

Threshold image pixel values below 0.5.

preimg=Image[Threshold[ImageAdjust[verimg],0.5],ImageSize->300]

With vertical components intensified, apply ImageLines.

lines=ImageLines[preimg,0.12];

Show[img,Graphics[{Thick,Blue,Opacity[0.5],Line/@lines}],ImageSize->300]

Neat Examples

Shift-invariance is an important property used in applications such as:

Detecting changes in time series

Denoising

Pattern recognition

Clear[swtfunc];
swtfunc[s_,wavelet_]:= Module[{swt,data},
data =Table[SawtoothWave[s-t],{t,-2,2,4/31}];
swt= wavelet[data,Automatic,3];
Row[{Plot[SawtoothWave[s-t],{t,-2,2},PlotRange->{0,1},PlotLabel->Style["Input Data",16,FontFamily->"Verdana"],Axes->{True,False},Exclusions->None,ImageSize->450],Show[,ImageSize->50],WaveletListPlot[swt,Ticks->Full,Filling->Axis,ImageSize->450,PlotLabel->Style["Wavelet Transform",16,FontFamily->"Verdana"]]}]
]
Manipulate[swtfunc[s,wavelet],{s,0,5,0.1},{wavelet,{StationaryWaveletTransform,DiscreteWaveletTransform}},SaveDefinitions->True]

Apply wavelet transform at each decomposition level and instead of sub-sampling the data, up-sample the filter by padding with zeros. This can be useful in applications such as statistical analysis.

swd = StationaryWaveletTransform[,Automatic,3]

There is no sub-sampling of wavelet coefficients.

swd["Dimensions"]

Directly extract wavelet coefficients as an Image.

swd[{_,_},"Image"]

Customize WaveletImagePlot using an image function applied to a specified wavelet index.

Clear[imgFunc];
imgFunc[img_,{0,0,1}]:= EdgeDetect[ImageAdjust[ImageApply[Abs,img]],0.7,Method-> "ShenCastan"]
imgFunc[img_,{0,2}]:= Sharpen@Colorize[ImageAdjust[ImageApply[Abs,img]],ColorFunction->"SolarColors"]
imgFunc[img_,{___,0}]:= Composition[ImageAdjust,ImageApply[Abs,#1]&][img]
imgFunc[img_,wind_]:= Composition[Sharpen,ImageAdjust[#,{0,1}]&,ImageAdjust,ImageApply[Abs,#1]&][img]

WaveletImagePlot[swd,All,imgFunc,ImageSize->250]

Denoising Images

Denoising of images is analogous to removing higher frequency components from the image.

nimg = ImageEffect[,{"GaussianNoise",0.2}];

nimg

dwd = DiscreteWaveletTransform[nimg,BiorthogonalSplineWavelet[5,5]];

Perform "Soft" thresholding with threshold value "SURE" computed adaptively at each level.

wtdwd = WaveletThreshold[dwd,{"Soft","SURELevel"},{1|2|3}];

Invert thresholded coefficients.

{Image[InverseWaveletTransform[wtdwd],ImageSize->All],Image[nimg,ImageSize->All]}

Threshold specific components of an image.

dwd=DiscreteWaveletTransform[,BiorthogonalSplineWavelet[2,2],3]

Strongly threshold horizontal and vertical detail coefficients {___,1|2} only.

dwd2 = WaveletThreshold[dwd,{"Hard",1},{___,1|2}]

Compare coefficients.

{Labeled[WaveletImagePlot[dwd,ImageSize->Small],Text@"original"],Labeled[WaveletImagePlot[dwd2,ImageSize->Small],Text@"thresholded"]}

The inverse transform image has mainly diagonal features.

{InverseWaveletTransform[dwd],InverseWaveletTransform[dwd2]}

Half-toning is the process of representing an image with a series of dots of different sizes. Most notably used in newspapers.

img =ImageTake[ ,100]

Perform a stationary wavelet transform on the half-toned image.

swt =StationaryWaveletTransform[img,DaubechiesWavelet[4],2];

Apply GaussianFilter on the detail coefficients. The inverse wavelet transform of the filtered coefficients gives a inverse half-toned image.

InverseWaveletTransform[WaveletMapIndexed[GaussianFilter[#,4]&,swt,Except[{___,0}]]]

Wavelet Thresholding

Generate noisy data.

Remove[noisy,dwd,thr1,thr2];

noisy=Table[Sin[x]+RandomReal[0.2],{x,-2π,2π,0.01}];

ListLinePlot[noisy,BaseStyle->Automatic]

Perform a StationaryWaveletTransform on the noisy data.

dwd=StationaryWaveletTransform[noisy,SymletWavelet[4],5]

Use a specific threshold function.

thr1=WaveletThreshold[dwd]

Also specify how to choose the threshold value.

thr2=WaveletThreshold[dwd,{"Soft","SURE"}]

Perform InverseWaveletTransform and compare.

Row[{ListLinePlot[noisy,PlotLabel->Style["original",18],ImageSize->300,BaseStyle->Automatic],ListLinePlot[InverseWaveletTransform[thr1],PlotLabel->Style[{"Hard","Universal"},18],ImageSize->300,BaseStyle->Automatic],ListLinePlot[InverseWaveletTransform[thr2],PlotLabel->Style[{"Soft","SURE"},18],ImageSize->300,BaseStyle->Automatic]}]

Compare some named automatic thresholding methods.

ToPictureCell@DynamicModule[{data,σ=0.1,lwt,waveThresh,plot,wave=DaubechiesWavelet[2],dataplot,wfunc},
data[σ_]:= Table[N[10Cos[(x^2)]]+RandomReal[NormalDistribution[0,σ]],{x,0,4,0.005}];

lwt[wfunc_,σ_]:=wfunc[data[σ],wave,4];

waveThresh[σ_,tspec_,wfunc_]:=InverseWaveletTransform[WaveletThreshold[lwt[wfunc,σ],tspec]];

plot[σ_,tspec_,wfunc_]:= ListLinePlot[waveThresh[σ,tspec,wfunc],PlotLabel->Style[tspec,FontFamily->"Versana"],PlotRange->{-12,12},DataRange->{0,4},ImageSize->300,AspectRatio->0.5,BaseStyle->Automatic];

dataplot[σ_]:= ListLinePlot[data[σ],PlotLabel->Style["Data",FontFamily->"Versana"],PlotRange->{-12,12},DataRange->{0,4},ImageSize->300,AspectRatio->0.5,BaseStyle->Automatic];

Grid[{{Row[{Style["Noise (σ) : ",Bold,14],Manipulator[Dynamic[σ],{0.1,1},Appearance->"Open"]}],Dynamic[dataplot[σ]],Row[{Style["Transform : ",Bold,14],PopupMenu[Dynamic[wfunc],{LiftingWaveletTransform,StationaryWaveletTransform}]}]},
{Dynamic[plot[σ,"Universal",wfunc]],Dynamic[plot[σ,"VisuShrink",wfunc]],Dynamic[plot[σ,"SUREShrink",wfunc]]},{Dynamic[plot[σ,"SUREHybrid",wfunc]],Dynamic[plot[σ,"GCV",wfunc]],Dynamic[plot[σ,"FDR",wfunc]]}},Dividers->All,Alignment->Left]]

Compare named automatic adaptive thresholding methods.

ToPictureCell@DynamicModule[{data,σ=0.1,lwt,waveThresh,plot,wave=DaubechiesWavelet[2],dataplot,wfunc},
data[σ_]:= Table[2-5t+5 Exp[-500(((t-0.5))^2)]+RandomReal[NormalDistribution[0,σ]],{t,0,1,(1/1023)}];

lwt[wfunc_,σ_]:=wfunc[data[σ],wave,4];

waveThresh[σ_,tspec_,wfunc_]:=InverseWaveletTransform[WaveletThreshold[lwt[wfunc,σ],tspec]];

plot[σ_,tspec_,wfunc_]:= ListLinePlot[waveThresh[σ,tspec,wfunc],PlotLabel->Style[tspec,FontFamily->"Versana"],PlotRange->{-3,5},DataRange->{0,1},ImageSize->300,AspectRatio->0.5,BaseStyle->Automatic];

dataplot[σ_]:= ListLinePlot[data[σ],PlotLabel->Style["Data",FontFamily->"Versana"],PlotRange->{-3,5},DataRange->{0,1},ImageSize->300,AspectRatio->0.5,BaseStyle->Automatic];

Grid[{{Row[{Style["Noise (σ) : ",Bold,14],Manipulator[Dynamic[σ],{0.1,1},Appearance->"Open"]}],Dynamic[dataplot[σ]],Row[{Style["Transform : ",Bold,14],PopupMenu[Dynamic[wfunc],{LiftingWaveletTransform,StationaryWaveletTransform}]}]},
{Dynamic[plot[σ,"UniversalLevel",wfunc]],Dynamic[plot[σ,"SURELevel",wfunc]],Dynamic[plot[σ,"GCVLevel",wfunc]]}},Dividers->All,Alignment->Left]]

Image and Data Compression

Wavelets have excellent compression properties. Using wavelet packet transforms, wavelet best basis, and wavelet thresholding, one can perform image and signal compression.

Visualize default reconstruction tree coefficients for image data.

dwd=DiscreteWaveletPacketTransform[,Automatic];

WaveletImagePlot[dwd,ImageSize->Small]

Compute the reconstruction tree whose coefficients have the smallest total log energy.

best=WaveletBestBasis[dwd,"ShannonEntropy"]

Comparing the default bases needed for reconstruction vs. the one selected by WaveletBestBasis using Shannon entropy, we need far less bases to reconstruct the original signal.

WaveletImagePlot[best,ImageSize->Small]

Use WaveletBestBasis to choose a wavelet basis for threshold-based data compression.

data=Table[Sin[x+x^2],{x,0,10,0.01}];

WaveletBestBasis[DiscreteWaveletPacketTransform[data,DaubechiesWavelet[4],4]]

Use a hard threshold of 0.1.

dwd=WaveletThreshold[%,{"Hard",0.1}]

Count nonzero values in thresholded coefficients as a measure of compression.

Count[Last/@dwd[Automatic],Except[0.],∞]

Compare with compression by thresholding original data.

Count[Chop[data,0.2],Except[0],∞]

Plot original and reconstructed data.

{ListLinePlot[data],ListLinePlot[InverseWaveletTransform[dwd]]}

Compare threshold-based compression for various threshold values δ.

data=Table[Sin[x+x^2],{x,0,10,0.01}];

dwd=DiscreteWaveletPacketTransform[data,DaubechiesWavelet[],4];

Use the same threshold value δ for best basis selection and thresholding.

dwds=Table[WaveletThreshold[WaveletBestBasis[dwd,{"Threshold",δ}],{"Hard",δ}],{δ,{0.05,0.1,0.2,0.5}}]

Show reconstructed data along with the number of nonzero values as a measure of compression.

Table[ListLinePlot[InverseWaveletTransform[d],PlotLabel->Row[{"non-zero values: ",Count[Last/@d[Automatic],Except[0.],∞]}]],{d,dwds}]

Length of uncompressed original data.

Length[data]

Simple lossy image compression for various threshold values δ.

dwd=DiscreteWaveletPacketTransform[,CDFWavelet[],3];

dwds=Table[WaveletThreshold[WaveletBestBasis[dwd,{"Threshold",δ}],{"Hard",δ}],{δ,{0.02,0.05,0.2,0.5}}]

Show reconstructed data along with the number of nonzero values as a measure of compression.

Table[Labeled[InverseWaveletTransform[d],Row[{"non-zero values: ",Count[Last/@d[Automatic],Except[0.],∞]}]],{d,dwds}]

Total number of pixel values in all channels of original image.

ImageChannels[]*Apply[Times,ImageDimensions[]]

Image Fusion

Wavelet analysis can be used to combine two images in a variety of ways.

Take the texture of the grass from the second image and overlay it on the first image. This requires extracting the edge features of the grass from the second image.

{img,txt} ={ ,};

Perform a wavelet transform on both images.

dwd1 = StationaryWaveletTransform[img, CDFWavelet[], 3];
dwd2 = StationaryWaveletTransform[txt, CDFWavelet[], 3];

The textures of each image are in the detail coefficients.

{WaveletImagePlot[dwd1,ImageSize -> 100],WaveletImagePlot[dwd2,ImageSize -> 100]}

Combine detail coefficients of the two images by taking their mean.

rval1 = dwd1[{___, 1 | 2 | 3}, "Values"];
rval2 = dwd2[{___, 1 | 2 | 3}, "Values"];
wind = Cases[dwd1["WaveletIndex"], {___, 1 | 2 | 3}];
nimg=(Image[#1,Interleaving->False]&)/@((1/2) (rval1+2 rval2));
rnew = MapThread[Rule, {wind, nimg}];

Append the coarse coefficient of the first image.

rr = Append[rnew, First[dwd1[{0, 0, 0}, {"Image", "ImageFunction" -> Identity}]]];

Construct a new DiscreteWaveletData of the combined wavelet coefficients.

dwdnew = DiscreteWaveletData[rr, CDFWavelet[], StationaryWaveletTransform]

Reconstruct the combined image.

InverseWaveletTransform[dwdnew]

Use a mask image to vary between blurring and sharpening across an image.

WaveletMapIndexed[ImageMultiply[#,]&,StationaryWaveletTransform[,Automatic,3],{___,1|2|3}]

Compare with the original image.

{InverseWaveletTransform[%],}

Other Applications

Differentiate noisy data using a wavelet transform.

ClearAll[t]

f[x_]:=Log[2+Sin[3 π (Sqrt[x])]]

{t,fd}=Transpose[Table[{x,f[x]+RandomReal[NormalDistribution[0,0.01]]},{x,0.01,1,0.001}]];

ListLinePlot[fd]

To get more details on this application refer to the Application section in StationaryWaveletTransform.

WaveletD[t_,fd_]:=Module[{a,b,wd,j,du,dx,dudx},a=((Max[fd]-Min[fd])/(Max[t]-Min[t]));b=Min[fd]-a Min[t];wd=StationaryWaveletTransform[fd-(a t+b),HaarWavelet[],4];j=wd["Refinement"];du=((Sqrt[2])^j) First[wd[{0,0,0,1},"Values"]];dx=0.001;dudx=((4 du)/(dx ((((2^j)))^(3/2))))+a;Transpose[{t,dudx}]]

Compare wavelet-based numerical differentiation with exact differentiation.

Show[
ListLinePlot[WaveletD[t,fd][[21;;-1]]],
Plot[Evaluate[(∂_x)f[x]],{x,0.03,1},PlotStyle->{Thick,Red}],PlotRange->All]

Compare with standard Mathematica numerical differentiation.

ifunc = Interpolation[Transpose[{t, fd}]];

Plot[ifunc'[x], {x, 0.01, 1}]

Integer lifting maps integers to integer wavelet coefficients. Integer lifting is very important because it does not convert integers into real when doing the decomposition. This results in considerable savings in memory.

data = RandomInteger[10,{64}];

lwt = LiftingWaveletTransform[data,Method->"IntegerLifting"]

TableForm[lwt[Automatic]]

Norm[data-InverseWaveletTransform[lwt]]

Approximate data using few wavelet coefficients.

waveApprox[n_]:= Module[{data,dwd,rdata},
data=Table[Abs[Sin[2 π x]]+1.5 Abs[Cos[2 π x-π]],{x,0,1,(1/((2^6)-1))}];
dwd=LiftingWaveletTransform[data,HaarWavelet[]];
rdata=InverseWaveletTransform[WaveletThreshold[dwd,{"LargestCoefficients",n},Automatic]];
ListLinePlot[{data,rdata},PlotRange->All,PlotStyle->{Red,Green}]
]
Manipulate[waveApprox[n],{n,1,64,1},SaveDefinitions->True]

References

General References

Burrus, C. Sidney, Gopinath, Ramesh A., and Guo, Haitao. Introduction to Wavelets and Wavelet Transforms: A Primer. Prentice-Hall, 1998.

Daubechies, I. Ten Lectures on Wavelets. Society for Industrial and Applied Mathematics, 1992.

Jensen, A. and la Cour-Harbo, A. Ripples in Mathematics: The Discrete Wavelet Transform. Springer, 2001.

Kaiser, G. A Friendly Guide to Wavelets. Birkhäuser, 1994.

Mallat, S. A Wavelet Tour of Signal Processing. Academic Press, 1999.

Strang, G. and Nguyen, T. Wavelets and Filter Banks. Wellesley-Cambridge Press, 1996.

Van Fleet, Patrick J. Discrete Wavelet Transformations: An Elementary Approach with Applications. Wiley-Interscience, 2008.

Walden, Andrew T. and Percival, Donald B. Wavelet Methods for Time Series Analysis. Cambridge University Press, 2000.

Walnut, David F. An Introduction to Wavelet Analysis. Birkhäuser, 2002.

Application References

Addison, N. The Illustrated Wavelet Transform Handbook. The Institute of Physics, 2002.

Chemometrics 2004.

Gençay, R., Selçuk, F. and Whitcher, B. An Introduction to Wavelets and Other Filtering Methods in Finance and Economics. Elsevier, 2002.

Jansen, M. Noise Reduction by Wavelet Thresholding. Springer, 2001.

Ogden, R. Todd. Essential Wavelets for Statistical Applications and Data Analysis. Birkhäuser, 1996.

Salesin, David H. Stollnitz, Eric J., Derose, Tony D. Wavelets for Computer Graphics: Theory and Applications. Morgan Kaufmann Publishers, Inc., 1996.

Tang, Yuan Y. Wavelet Theory and Its Application to Pattern Recognition, 2nd edition. World Scientific, 2009.

Urban, K. Wavelet Methods for Elliptic Partial Differential Equations. Oxford University Press, 2009.

Vidakovic, B. Statistical Modeling by Wavelets. Wiley, 1999.

Initialization

fourier[data_,dt_]:=Module[{nn,fdata,freqlist,data1},
data1=data-Mean[data];
nn=Length[data1];
fdata=Abs[Fourier[data1,FourierParameters->{-1,1}]];
fdata=(((2 fdata[[1;;Round[(nn/2)]]]))^2);
freqlist=((Range[1,Round[(nn/2)]]-1)/((nn-1) dt));
Transpose[{2*freqlist,fdata}]];

imgsize=500;

sty =Directive[13,FontFamily->"Verdana"];

SetOptions[Plot,BaseStyle->Thick,PlotStyle->ColorData[1,"ColorList"],ImageSize->imgsize,PlotRange->All,FillingStyle->Automatic,TicksStyle->Directive[13,FontFamily->"Verdana"]];

SetOptions[Plot3D,ImageSize->imgsize,PlotRange->All,TicksStyle->Directive[13,FontFamily->"Verdana"],Mesh->None,PlotPoints->50];

SetOptions[ListPlot3D,ImageSize->imgsize,PlotRange->All,TicksStyle->Directive[13,FontFamily->"Verdana"],Mesh->None];

SetOptions[DateListPlot,ImageSize->imgsize,PlotRange->All,Joined->True,Frame->False,TicksStyle->Directive[13,FontFamily->"Verdana"]];

SetOptions[ListLinePlot,BaseStyle->Thick,PlotStyle->ColorData[19,"ColorList"],ImageSize->imgsize,PlotRange->All,FillingStyle->Automatic,TicksStyle->Directive[13,FontFamily->"Verdana"],FillingStyle->Lighter/@ColorData[19,"ColorList"]];

(SetOptions[#1,TicksStyle->Directive[13,FontFamily->"Verdana"],FrameTicksStyle->Directive[13,FontFamily->"Verdana",FrameTicks->Full],ImageSize->imgsize]&)/@{WaveletMatrixPlot,MatrixPlot,WaveletScalogram};

SetOptions[WaveletListPlot,PlotStyle->ColorData[19,"ColorList"],ImageSize->imgsize,PlotRange->All,TicksStyle->Directive[13,FontFamily->"Verdana"],FillingStyle->Lighter/@ColorData[19,"ColorList"]];

SetOptions[Manipulate,DefaultLabelStyle->Directive[13,Bold,FontFamily->"Verdana"]];

SetOptions[TabView,DefaultLabelStyle->Directive[12,FontFamily->"Verdana"]];

ToPictureCell[expr_]:=CellPrint[ExpressionCell[expr,"Picture",ShowStringCharacters->False]]

SetDirectory[NotebookDirectory[]];

$Line=0;