Wavelet Analysis: Concepts

Wolfram Research
wolfram-u@wolfram.com

MouseOver
MouseOver
MouseOver
MouseOver
MouseOver
MouseOver
MouseOver

Introduction: Fourier 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->550,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->550,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)],50],ImageSize->500,AspectRatio->0.4,ImagePadding->{{30,Automatic},{Automatic,Automatic}},AxesLabel -> {Style["2 frequency",Bold,18], ""}]},{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

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 ω+Δω.

Example: Continuous Wavelet Transform

ToPictureCell[Labeled[Block[{data,fourier,nst},
nst[t_]:=({{\[Piecewise], ({{Cos[32 π t], 0≤t<2}, {2 Cos[64 π t], 2≤t<4}, {3 Cos[128 π t], 4≤t<8}, {0, True}})}});
data=N[Table[nst[t],{t,0,6,(6/2047)}]];
Column[{
ListLinePlot[data,ImageSize->700,AspectRatio->0.15,PlotRangePadding->None,PlotStyle->Directive[ColorData["DarkRainbow",0.2]],PlotLabel->Style["Non Stationary Data",FontFamily->"Verdana",Brown,Bold,16],ImagePadding->{{30,Automatic},{Automatic,Automatic}},PlotRange->{{0,6},{-4,4}},AxesOrigin->{0,-4},DataRange->{0,6},BaseStyle->Automatic,AxesLabel->{Style["Time",Directive[13,FontFamily->"Verdana"]]}]
,
WaveletScalogram[ContinuousWaveletTransform[data,GaborWavelet[6],{6,12},SampleRate->(2047/6)],ImageSize->700,PlotRangePadding->None,PlotLabel->Style["Wavelet Decomposition",FontFamily->"Verdana",Brown,Bold,16],ImagePadding->{{30,Automatic},{Automatic,Automatic}},AxesLabel->{Style["Time",Directive[13,FontFamily->"Verdana"]]},AspectRatio->0.5]},Alignment->Right]],Style[TraditionalForm[({{\[Piecewise], ({{Cos[64 π 𝓉], 0≤𝓉<2}, {Cos[128 π 𝓉], 2≤𝓉<4}, {Cos[256 π 𝓉], 4≤𝓉<8}, {0, True}})}})],Directive[sty,18]],Right]]

Example: Discrete Wavelet Transform

ToPictureCell[Labeled[Block[{data,fourier,nst},
nst[t_]:=({{\[Piecewise], ({{Cos[64 π t], 0≤t<2}, {2 Cos[128 π t], 2≤t<4}, {3 Cos[256 π t], 4≤t<8}, {0, True}})}});
data=N[Table[nst[t],{t,0,6,(6/2047)}]];
Column[{
ListLinePlot[data,ImageSize->700,AspectRatio->0.15,PlotRangePadding->None,PlotStyle->Directive[ColorData["DarkRainbow",0.2]],PlotLabel->Style["Non Stationary Data",FontFamily->"Verdana",Brown,Bold,16],ImagePadding->{{30,Automatic},{Automatic,Automatic}},PlotRange->{{0,6},{-4,4}},AxesOrigin->{0,-4},DataRange->{0,6},BaseStyle->Automatic,AxesLabel->{Style["Time",Directive[13,FontFamily->"Verdana"]]}]
,
WaveletListPlot[StationaryWaveletTransform[data,ShannonWavelet[25],3],Automatic,DataRange->{0,6},PlotStyle->(ColorData["DarkRainbow"]/@Riffle[#1,#2]&)@@Partition[Range[0.05,0.85,(0.8/7)],4],ImageSize->700,PlotLayout->"CommonXAxis",PlotRangePadding->None,PlotLabel->Style["Wavelet Decomposition",FontFamily->"Verdana",Brown,Bold,16],ImagePadding->{{30,Automatic},{Automatic,Automatic}},AxesLabel->{Style["Time",Directive[13,FontFamily->"Verdana"]]},AspectRatio->0.5]},Alignment->Right]],Style[TraditionalForm[({{\[Piecewise], ({{Cos[64 π t], 0≤t<2}, {2 Cos[128 π t], 2≤t<4}, {3 Cos[256 π t], 4≤t<8}, {0, True}})}})],Directive[sty,18]],Right]]

Which transform to choose in analyzing a given signal is task-dependent. Following are some of the questions to ask:

What am I looking at in a signal?

Based on the purpose of analysis, what’s a suitable transform to use?

What are the advantages and disadvantages of a given analysis?

Introduction: Wavelet Analysis

Wavelets are short-lived wavelike oscillations whose amplitude starts at 0 increases and decreases back to 0.

Plot[WaveletPsi[MexicanHatWavelet[1],x],{x,-5,5},PlotRange->All,Axes->None]

Complex wavelet.

psi=WaveletPsi[GaborWavelet[6],x];
Plot[{Re[psi],Im[psi]},{x,-5,5},Axes->None,PlotStyle->{Darker[Red],{Dashed,Darker[Green]}}]

Wavelets are localized in space (time) and scale (frequency) domains. This localization feature makes the signal “sparse” in the wavelet domain, which can be used in a number of useful applications such as data compression, detecting features, and removing noise from signals.

Wavelets are based on a multi-resolution analysis (can analyze a signal at different scales).

Wavelets are characterized within a wavelet family by the number of vanishing moments. The higher the vanishing moment, the better the smooth signal approximation in the wavelet domain.

Manipulate[Plot[WaveletPsi[DGaussianWavelet[10],((x-k)/s)],{x,-10,10},PlotRange->All,PlotStyle->Directive[Darker[Red],Thick],ImageSize->Large],{{s,2,"Scaling"},1,5},{{k,0,"Translation"},-5,5}]

Introduction

The continuous wavelet transform of a function x(t) is given by

w(u,s)=(1/(Sqrt[s]))(∫_(-∞)^∞)x(t) (Conjugate[ψ])(((t-u)/s))dt

ContinuousWaveletData

The continuous wavelet transform of a function x(t) is given by

w(u,s)=(1/(Sqrt[s]))(∫_(-∞)^∞)x(t) (Conjugate[ψ])(((t-u)/s))dt

Perform a ContinuousWaveletTransform.

ListPlot[Sin[Range[40]/5]]

cwt=ContinuousWaveletTransform[Sin[Range[40]/5]]

The resulting ContinuousWaveletData is a symbolic representation of the transformed data. It represents arrays of coefficients at different scales {oct,voc}.

WaveletScalogram[cwt,Ticks->Full]

Extract properties including numerical scale corresponding to each pair, {oct,voc}.

cwt["Scales"]

Compute the inverse wavelet transform of arbitrary continuous wavelet transform coefficients.

Norm[InverseContinuousWaveletTransform[cwt]-Sin[Range[40]/5]]//Chop

cwd=ContinuousWaveletTransform[Sin[Range[8]/5],Automatic,{2,2}]

Use Normal to get all wavelet coefficients explicitly.

Normal[cwd]

Use All as an argument to get all coefficients.

cwd[All]

Extract specific coefficient arrays:

cwd[{_,1}]

Extract several wavelet coefficients corresponding to the list of wavelet index specifications.

cwd[{{1,1},{2,1}}]

Extract all coefficients whose wavelet indexes match a pattern.

cwd[{{_,1}}]

Plot coefficients as simple ListPlot icons.

cwd=ContinuousWaveletTransform[Table[Sin[(t^2)],{t,-3 π,3 π,((6 π)/1023)}],MexicanHatWavelet[],{8,1}]

cwd[All,"ListPlot"]

Get Dimensions of the input data.

cwd["DataDimensions"]

Use WaveletMapIndexed to operate on ContinuousWaveletData.

cwd=ContinuousWaveletTransform[{1,2,3,4}]

cwd=WaveletMapIndexed[c|->c*0.2,#]&[%]

InverseContinuousWaveletTransform[cwd]

Octaves, Voices, and Scales

The default number of octaves is Log2(n/2). The default number of voices is 4.

cwt=ContinuousWaveletTransform[Sin[Range[40]/5]];
{{cwt["Octaves"],cwt["Voices"]}, Round[Log2[40/2]]}

Each {oct,voc} is associated with a scale. Higher octaves have higher scales. The scale is computed as ((α2)^(oct-1))(2^(voc/nvoc)). The term α is the smallest wavelet scale.

cwt["Scales"]
cwt["WaveletIndex"]/. {oct_,voc_} -> cwt["WaveletScale"]*(2^(oct-1))*(2^(voc/4.))

The smallest wavelet scale depends on the type of wavelet used and the center frequency of that wavelet.

{1./cwt["Wavelet"]["FourierFactor"],cwt["WaveletScale"]}

Higher scales resolve lower frequencies. Lower scales resolve higher frequency. Scale is like using a microscope with different magnifications. Lower magnifications will show more with less detail and higher magnification will show more details locally.

data=Table[Sin[4 π t]+2 Exp[-(10^5) ((((1/3)-t))^2)],{t,0,1,(1/511)}];
p1 = ListLinePlot[data] ; 
cwt=ContinuousWaveletTransform[data,PaulWavelet[5],{7,12},SampleRate->511];
p2 = WaveletScalogram[cwt,ColorFunction->"AvocadoColors",Ticks -> Full,AxesLabel -> {"Time", "{Oct.,Voices}"}];
GraphicsRow[{p1,p2}]

Relating Scales and Frequency

There is a relationship between the wavelet scales and the frequency of a signal. Consider a signal with a single frequency.

xval=N[Range[0,1,(1/1023)]];
generateData[fr_]:=Cos[2 π fr xval];

The frequencies are proportional to the reciprocal of the wavelet scales and the center frequency of the wavelet.

myWavelet[order_]:= GaborWavelet[order];
cwtfunc[fr_]:=ContinuousWaveletTransform[generateData[fr],myWavelet[6],{8,4},SampleRate->1023]
cwt1 = cwtfunc[10];

Converting the scales into frequency.

freq = (#[[1]] -> 1023/(#[[2]]*GaborWavelet[6]["FourierFactor"]))& /@ cwt1["Scales"]

Manipulate[
cwt1 = ContinuousWaveletTransform[generateData[fr],myWavelet[order],{8,4},SampleRate->1023];
freq = (1023/(#*GaborWavelet[order]["FourierFactor"]))& /@ (Thread[{Range[8],1}]/.cwt1["Scales"]) ;
ticks = Transpose[{Range[Length[freq]],freq}];
Column[
{Grid[{{"Frequency (Hz):",fr}}],ListLinePlot[generateData[fr],AspectRatio->0.2,Axes->False],
WaveletScalogram[cwt1,Frame -> True, FrameTicks -> {{ticks,Automatic},Automatic},FrameLabel -> {"Time","Frequency"}]}],
{{order,6,"Order of the wavelet"},6,20,1},{{fr,10,"Frequency of signal"},10,100}]

Discrete Wavelets: Scaling Function

Scaling Function

Wavelet functions are basically band-pass filters and at each refinement level the frequency bandwidth is essentially halved. Theoretically, this corresponds to an infinite refinement levels decomposition in order to cover the entire spectrum. The scaling function captures the low-pass spectrum in order to fill in the gaps.

The scaling function ϕ(x) satisfies the following recursive equation:

ϕ(x)==2 (∑_n)(Subscript[a, n]) ϕ(2 x-n)

where (Subscript[a, n]) are the low-pass filter coefficients.

WaveletPhi[HaarWavelet[],x,WorkingPrecision->∞]

Clear[ϕ];
ϕ=WaveletPhi[DaubechiesWavelet[4]];
ϕ[x]

Plot[ϕ[x],{x,0,7}]

The low-pass filter coefficients can be found using WaveletFilterCoefficients.

lowpass = WaveletFilterCoefficients[DaubechiesWavelet[4],"PrimalLowpass"]

Using the translation and the scale, one can construct the scaling functions.

phi = WaveletPhi[DaubechiesWavelet[4]] ; 
phiX =2* Total[lowpass /. {translation_,coeff_} -> coeff*phi[2*x - translation]] ;

GraphicsRow[{Plot[phiX,{x,0,3}],Plot[phi[x],{x,0,3}]}]

Low-pass filter coefficients sum to unity.

Total[lowpass[[All,2]]]

Discrete Wavelets: Wavelet Function

Wavelet Function

The wavelet function ψ(x) is defined by

ψ(x)==2 (∑_n)(Subscript[b, n]) ϕ(2 x-n)

where (Subscript[b, n]) represents the high-pass filter coefficients.

WaveletPsi[HaarWavelet[],x,WorkingPrecision->∞]

Discrete wavelets are characterized by the number of coefficients and by the level of iteration.

Clear[ψ];
ψ=WaveletPsi[DaubechiesWavelet[4]];
ψ[x]

Plot[ψ[x],{x,-3,4},PlotRange->All]

The high-pass filter coefficients can be found using WaveletFilterCoefficients.

highpass = WaveletFilterCoefficients[DaubechiesWavelet[4],"PrimalHighpass"]

Using the translation, the coefficients and scaling functions, we can construct the wavelet functions.

phi = WaveletPhi[DaubechiesWavelet[4]] ; 
psiX =2* Total[highpass /. {translation_,coeff_} -> coeff*phi[2*x - translation]] ;

GraphicsRow[{Plot[psiX,{x,-1,2}],Plot[WaveletPsi[DaubechiesWavelet[4],x],{x,-1,2}]}]

High-pass filter coefficients sum to zero.

Total[highpass[[All,2]]]

Multivariate scaling and wavelet functions are outer products of univariate ones:

ToPictureCell@Grid[{{Plot3D[Evaluate[ϕ[x] ϕ[y]],{x,0,7},{y,0,7},ColorFunction->"DarkRainbow",ImageSize->350,PlotLabel->Style["ϕ(x) ϕ(y)",18]],
Plot3D[Evaluate[ϕ[x] ψ[y]],{x,0,7},{y,-3,4},ColorFunction->"DarkRainbow",ImageSize->350,PlotLabel->Style["ϕ(x) ψ(y)",18]]},{Plot3D[Evaluate[ψ[x] ϕ[y]],{x,-3,4},{y,0,7},ColorFunction->"DarkRainbow",ImageSize->350,PlotLabel->Style["ψ(x) ϕ(y)",18]],
Plot3D[Evaluate[ψ[x] ψ[y]],{x,-3,4},{y,-3,4},ColorFunction->"DarkRainbow",ImageSize->350,PlotLabel->Style["ψ(x) ψ(y)",18]]}}]

Wavelet Filter Bank

A wavelet filter bank is a set of low-pass and high-pass filters that separate the input signals into frequency bands.

The simplest example of a wavelet filter is the HaarWavelet. A lowpass filter performs a moving average over the time series.

f(t)=(1/2) x(t-1)+(1/2)x(t)

The low-pass filter coefficients are {(1/2),(1/2)}, which can be computed in Mathematica as follows:

WaveletFilterCoefficients[HaarWavelet[],"PrimalLowpass",WorkingPrecision->∞]

Frequency Response

To compute the frequency response to the filter we compute the components of the input vector to have pure frequency ω:

x(t)=(e^(i t ω))

For the moving average the response is:

f(t)=(1/2)x(t)+(1/2) x(t-1)
f(t)= (1/2)(e^(i t ω))+(1/2)(e^(i (t-1) ω))
f(t)= ((1/2)+(1/2)(e^(-i ω)))(e^(i t ω))

Low-pass frequency response function:

(Subscript[H, 0])(ω)=(1/2)+(1/2)(e^(-i ω))=cos((ω/2)) (e^(-i (ω/2)))

Plot[Abs[Cos[(ω/2)]],{ω,-π,π},Ticks->{{-π,π},{1}},TicksStyle->Directive[20,Bold],PlotStyle->Thick,Axes->True,AxesLabel->{Style[ω,20],Style[(Subscript[H, 0]) ω,20]}]

At ω=0, frequency response function (Subscript[H, 0])(ω) is unity, hence the name low-pass filter. 
The frequency response function for other filters can be computed in similar fashion.

The low-pass filter does an averaging operation, which in effect smooths out the spikes in the signal. The high-pass filter in contrast performs moving differences. It captures the spikes in the signal.

f(t)=(1/2)x(t)-(1/2) x(t-1)

The high-pass filter coefficients are {(1/2),-(1/2)}, which can be computed in Mathematica as follows:

WaveletFilterCoefficients[HaarWavelet[],"PrimalHighpass",WorkingPrecision->∞]

Frequency Response

To compute the frequency response to the high-pass filter we compute the components of the input vector to have pure frequency ω:

x(t)=(e^(i t ω))

For the moving average the response is:

f(t)=(1/2)x(t) - (1/2) x(t-1)
f(t)=(1/2)(e^(i t ω)) - (1/2)(e^(i (t-1) ω))
f(t)= ((1/2)-(1/2)(e^(-i ω)))(e^(i t ω))

High-pass frequency response function

(Subscript[H, 1])(ω)=(1/2)-(1/2)(e^(-i ω))=sin((ω/2))i (e^(-i (ω/2)))

Plot[Abs[Sin[(ω/2)]],{ω,-π,π},Ticks->{{-π,π},{1}},TicksStyle->Directive[20,Bold],PlotStyle->Thick,AxesLabel->{Style[ω,20],Style[(Subscript[H, 1]) ω,20]}]

At ω=0, frequency response function is zero. At the highest frequency (ω=π, ω=-π) , the response function (Subscript[H, 1])(ω)  is unity, hence the name high-pass filter.

Parametrization of Discrete Wavelets

Scaling functions ϕ(x) are generated from a single scaling function by simple scaling (j) and translation (k); these are the basic operations used in wavelet expansion:

(Subscript[ϕ, j,k])(x)=(2^(j/2)) ϕ ((2^j) x-k)

The factor (2^(j/2)) maintains a constant norm independent of the scale j.

Clear[ϕ];
ϕ[x_,j_,k_]:=(2^(j/2)) WaveletPhi[DaubechiesWavelet[4],(2^j) x-k]
ToPictureCell@
Row[{Plot[Evaluate[Table[ϕ[x,j,0],{j,0,4}]],{x,0,2},Filling->Axis,Axes->None,PlotRange->All,PlotLabel->Style["Dilations",sty,20,Brown]],Plot[Evaluate[Table[ϕ[x,2,k],{k,0,(2^2)-1}]],{x,0,2},Filling->Axis,Axes->None,PlotLabel->Style["Translations",sty,20,Brown]]}]

Similarly, a two-dimensional parametrization of a wavelet function ψ(x) is performed as

(Subscript[ψ, j,k])(x)=(2^(j/2)) ψ ((2^j) x-k)

Clear[ψ]
ψ[x_,j_,k_]:=(2^(j/2)) WaveletPsi[DaubechiesWavelet[4],(2^j) x-k]
ToPictureCell@Row[{Plot[Evaluate[Table[ψ[x,j,0],{j,0,4}]],{x,-0.5,2},Axes->None,Filling->Axis,PlotRange->All,PlotLabel->Style["Dilations",sty,20,Brown]],Plot[Evaluate[Table[ψ[x,2,k],{k,0,(2^2)-1}]],{x,-0.5,1.5},Axes->None,Filling->Axis,PlotRange->All,PlotLabel->Style["Translations",sty,20,Brown]]}]

DiscreteWaveletData

The forward transform is given by

(Subscript[c, j+1,n])=(Sqrt[2])(∑_m)(Subscript[a, m-2n+2]) (Subscript[c, j,m])

(Subscript[d, j+1,n])=(Sqrt[2])(∑_m)(Subscript[b, m-2n+2]) (Subscript[c, j,m])

where (Subscript[c, 0]) is the original signal. The coefficients a and b are the wavelet low-pass and high-pass filter coefficients. The inverse transform is given by

(Subscript[c, j,n])=(Sqrt[2]) ((∑_m)(Subscript[a, n-2 m+2]) (Subscript[c, j+1,n])+ (∑_m)(Subscript[b, n-2 m+2]) (Subscript[d, j+1,m]))

Simple Transform

dwt=DiscreteWaveletTransform[{0,0,1,0,0},WorkingPrecision->∞]

The resulting DiscreteWaveletData is a symbolic representation of the transformed data.

dwt["TreeView"]

The blue frames represent the basis used for reconstruction.

dwt["BasisIndex"]

Packet Transform

dwpt=DiscreteWaveletPacketTransform[{0,0,1,0,0},WorkingPrecision->∞]

The resulting DiscreteWaveletData is a symbolic representation of the transformed data.

dwpt["TreeView"]

The blue frames represent the basis used for reconstruction.

dwpt["BasisIndex"]

Coarse and Detail Coefficients in DWT

Use a Haar Wavelet to decompose a signal and perform four levels of decomposition.

data = Table[t*Sin[(t^2)],{t,-3π,3π,0.01}] ; 
dwd = DiscreteWaveletTransform[data,HaarWavelet[],4];
ListLinePlot[data]

The coarse coefficients resolve the lower frequencies (the center part of the signal) while the detail coefficients resolve the higher frequencies (the outer part of the signal).

GraphicsRow[{WaveletListPlot[dwd,Ticks -> Full,PlotLabel -> "Wavelet Coefficients"],WaveletListPlot[dwd,Method->"InverseTransform",Ticks -> Full, PlotLabel -> "Inverse transform of coefficients"]}]

At each refinement level, the length of the data is halved. So coefficients at level 1 have length (n/2), coefficients at level 2 have length (n/4), and so on. Therefore, the data is down sampled.

{Length[data],dwd["Dimensions"]}

In order to reconstruct the signal, one needs the highest-level coarse coefficient AND all the detail coefficients.

dwd["TreeView"]

The coarse coefficients c at level j are computed by using the coarse and detail coefficients at level j+1.

coarse1 = InverseWaveletTransform[dwd,Automatic,{0,0,0,0}] + InverseWaveletTransform[dwd,Automatic,{0,0,0,1}] ;

Norm[coarse1 - InverseWaveletTransform[dwd,Automatic,{0,0,0}]]

The sum of inverse transforms from individual coefficient arrays gives the original data.

reconData = Total[InverseWaveletTransform[dwd,Automatic,#]& /@ dwd["BasisIndex"]] ;

Norm[reconData - data]

Wavelet Best Basis

Wavelet best basis is a minimal/optimal representation of the data in the wavelet domain relative to a particular cost function. Best basis is primarily used in applications related to denoising and data compression.

Recall the difference between wavelet transform and wavelet packet transform. Also, recall that in order to get coefficients at level j, we need information from level j+1.

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

dwpt=DiscreteWaveletPacketTransform[data,BiorthogonalSplineWavelet[10,8],4,Padding->0];
dwt=DiscreteWaveletTransform[data,Automatic,4];
GraphicsRow[{dwpt["TreeView"], dwt["TreeView"]}, ImageSize -> 900]

The bases in DWT are fixed. We need these bases to reconstruct the signal.

dwt["BasisIndex"]

The bases in DWPT are not fixed. One could mix and match which bases to use in order to reconstruct the original signal.

GraphicsRow[{dwpt["BestBasisBlockView"], WaveletBestBasis[dwpt,{"Threshold",10^-5}]["BestBasisBlockView"]}, ImageSize -> 900]

WaveletBestBasis tells us which of the bases should be kept for optimal reconstruction of signal. Any coefficients below a user-defined threshold are ignored.

δ=(10^(-5));
best =WaveletBestBasis[dwpt,{"Threshold",δ}];
best["BestBasisBlockView"]

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

BarChart[{{Count[Flatten[data],u_/;Abs[u]>δ],Count[Flatten[dwpt[Automatic,"Values"]],u_/;Abs[u]>δ],Count[Flatten[best[Automatic,"Values"]],u_/;Abs[u]>δ]}},ChartLabels->Placed[{"Data","Basic Packet","Best Basis"},Below,Style[#,sty]&],LabelingFunction->(f|->Placed[f,Above,Style[#,sty]&]),ImageSize->Large]

Introduction

Thresholding generally is based on two core pieces of information:

Thresholding function to be applied to wavelet coefficients

Method used to compute the threshold value δ

:  Hard threshold

:  Soft threshold

:  Firm threshold

:  Hyperbola threshold

:  Piecewise Garrote threshold

:  Smooth Garrote threshold

Donoho and Johnstone’s Universal thresholding

Stein’s unbiased risk estimate (SURE) thresholding

Generalized cross validation (GCV) thresholding

Data adaptive thresholding methods

Hybrid thresholding methods based on SURE thresholding

Thresholding based on False Discovery Rate (FDR)

Thresholding based on retaining maximum energy coefficients

Support for user-defined thresholding methods

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 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]]

Wavelet Thresholds

Here we demonstrate the details of computing a "Universal" type thresholding. For other methods, please refer to the documentation on WaveletThreshold.

"Universal" thresholding value.

data = Table[Tanh[(t^2)],{t,-2π,2π,4π/1023}];
nd = data+RandomReal[NormalDistribution[0,0.03],Length[data]];
ListLinePlot[nd,PlotRange->All]

Perform a wavelet transform.

dwd = DiscreteWaveletPacketTransform[nd,SymletWavelet[4],4];

The mean absolute deviation is computed using detail coefficients at first refinement level {1}.

(Subscript[σ, MAD])=((MedianDeviation[{1}/.dwd[{1}]])/(Quantile[NormalDistribution[],(3/4)]))

Compute universal threshold value δ as (Subscript[σ, MAD])(Sqrt[2 log n])  and perform a "Hard" threshold.

δ = (Subscript[σ, MAD])(Sqrt[2Log[Length[data]]])

wd1 = WaveletThreshold[dwd,{"Hard",δ}];

Use WaveletThreshold to compute the "Universal" threshold value.

wd2 = WaveletThreshold[dwd,"Universal"];

Compare threshold values.

Row[{wd1["ThresholdTable"],"    ", wd2["ThresholdTable"]}]

Introduction

The forward transform is given by

(Subscript[c, j+1,n])=(∑_(m=1)^f)(Subscript[a, m]) (Subscript[c, j,n-(2^j) (1-m mod l)])

(Subscript[d, j+1,n])=(∑_(m=1)^f)(Subscript[b, m]) (Subscript[c, j,n-(2^j) (1-m mod l)])

where f is the filter length for the corresponding wavelet and l is the length of input data.

The inverse transform is given by

(Subscript[c, j,n])=(∑_(m=1)^f)(Subscript[a, m]) (Subscript[c, j+1,(2^j) ( m mod d-1)+n])+(Subscript[b, m]) (Subscript[d, j+1,(2^j) (m mod d-1)+n])

Stationary Wavelet Transform

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

Detecting changes in time series

Denoising

Pattern recognition

Clear[swtfunc];
swtfunc[s_]:= Module[{swt,data},
data =Table[SawtoothWave[s-t],{t,-2,2,4/31}];
swt= StationaryWaveletTransform[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],{s,0,5,0.1},SaveDefinitions->True]

Apply a 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"]

Compare with a discrete wavelet transform.

dwd = DiscreteWaveletTransform[, Automatic,3];
dwd["Dimensions"]

Comparison using WaveletImagePlot shows the difference in the decomposition.

{WaveletImagePlot[swd], WaveletImagePlot[dwd, ImageSize -> 450]}

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;