Notebook on wavelet theory

belonging to the master course on Fourier Theory and Signal Analysis by T.H. Koornwinder (T.H.Koornwinder@uva.nl)

University of Amsterdam, February-June 2006

Last modified 7 June 2014

Reference: I. Daubechies, Ten lectures on wavelets, Regional Conference Series in Applied Math. 61, SIAM, 1992

Load packages

no spelling warnings:

Off[General::spell, General::spell1]

Pyramid (subband) algorithm conected to Daubechies wavelets with N=7

Compute the wavelet coefficients:

pp[n_,y_]:=Sum[Binomial[n-1+j,j]y^j,{j,0,n-1}]
qq[n_]:=Expand[z^(n-1)*pp[n,1/2-(z+z^(-1))/4],z]
roots[n_]:=Module[{sol},sol=NSolve[qq[n]==0,z];z/.sol]
diskroots[n_]:=Select[roots[n],Abs[#]<1&]
bb[n_]:=Sqrt[Abs[Coefficient[qq[n],z,2n-2]]]*Sqrt[Abs[Times@@diskroots[n] ]^(-1)]*2^(-n+1/2)*(1+z)^n z^(-2n+1)Product[(z-diskroots[n][[j]]),{j,1,n-1}]
hh[n_]:=Table[Chop[Coefficient[Expand[bb[n],z],z,-k]],{k,0,2n-1}]
hhh=hh[7]
ggg=Table[(-1)^k  hhh[[14-k]],{k,0,13}]

Normalization checks of the wavelet coefficients:

Sum[hhh[[i]],{i,1,14}]/Sqrt[2]

Chop[Sum[ggg[[i]],{i,1,14}]/Sqrt[2]]

Orthogonality checks of the wavelet coefficients:

Sum[hhh[[i]]^2,{i,1,14}]

Sum[ggg[[i]]^2,{i,1,14}]

Chop[Sum[hhh[[i]] ggg[[i]],{i,1,14}]]

For[j=1,j≤6,j++,
Print[j,", ",
Chop[Sum[hhh[[i]] hhh[[i+2j]],{i,1,14-2j}]],", ",
Chop[Sum[ggg[[i]] ggg[[i+2j]],{i,1,14-2j}]],", ",
Chop[Sum[hhh[[i]] ggg[[i+2j]],{i,1,14-2j}]]
]
]

f is the function to be processed:

f=Exp[-Pi x^2] (Sum[Random[Real,{-1,1}]/(k+1) Cos[k Pi x],{k,0,20}]+Sum[Random[Real,{-1,1}]/(k+1) Sin[k Pi x],{k,1,20}]);
Plot[f,{x,-2,2}, PlotStyle->Red,PlotRange->All,AxesOrigin->{-2,0}]

Initialize the c and d arrays:

c=Table[0,{j,-7,-1},{k,-2*2^7,3*2^7}];
d=Table[0,{j,-7,-1},{k,-3*2^7,3*2^7}];

c at resolution 2^(-7) is obtained by sampling the function f:

For[k=-2*2^7,k≤2*2^7,k++,c[[1,k+2*2^7+1]]=2^(-7/2)f/.x->k*2^(-7)];

Obtain c and d at resolutions 2^(-6), 2^(-5),  2^(-4), 2^(-3). Below, the red graphs are c at resolutions 2^(-7), 2^(-6), 2^(-5),  2^(-4), 2^(-3), the green graphs are d at resolutions 2^(-6), 2^(-5),  2^(-4), 2^(-3).

ListPlot[Table[c[[1,k+257]]*2^(7/2),{k,-256,256}],PlotStyle->Red,Joined->True,PlotRange->All,Ticks->{Table[{128i +257,i},{i,-2,2}],Automatic},PlotLabel->{"c",-7}]
For[j=-6,j≤-3,j++,
For[k=-2^(1-j),k≤2^(1-j),k++,
c[[j+8,k+2^(1-j)+1]]=
Sum[c[[j+7,2k+2^(2-j)+1+i]]  hhh[[i+1]],{i,0,13}];
d[[j+8,k+2^(1-j)+1]]=Sum[c[[j+7,2k+2^(2-j)+1+i]]  ggg[[i+1]],{i,0,13}]
];
pc=ListPlot[Table[c[[j+8,k+2^(1-j)+1]]2^(-j/2),{k,-2^(1-j),2^(1-j)}],PlotStyle->Red,Joined->True,PlotRange->All,Ticks->{Table[{2^(-j)i +2^(-j+1)+1,i},{i,-2,2}],Automatic},PlotLabel->{"c",j},DisplayFunction->Identity];
pd=ListPlot[Table[d[[j+8,k+2^(1-j)+1]]2^(-j/2),{k,-2^(1-j),2^(1-j)}],PlotStyle->Green,Joined->True,PlotRange->All,Ticks->{Table[{2^(-j)i +2^(-j+1)+1,i},{i,-2,2}],Automatic},PlotLabel->{"d",j},DisplayFunction->Identity];
Show[{pc,pd}]//Print
]

Initialize the array cr:

cr=Table[0,{j,-7,-3},{k,-2*2^7,3*2^7}];

cr will reconstruct c from c at resolution 2^(-3) and all the d. The graphs compare c (red) with cr (blue).

For[k=-16,k≤16,k++,cr[[5,k+17]]=c[[5,k+17]]
];

pc=ListPlot[Table[c[[5,k+17]]2^(3/2),{k,-16,16}],PlotStyle->Red,Joined->True,PlotRange->All,Ticks->{Table[{8i +17,i},{i,-2,2}],Automatic},PlotLabel->{"c",-3},DisplayFunction->Identity];
pcr=ListPlot[Table[cr[[5,k+17]]2^(3/2),{k,-16,16}],PlotStyle->Blue,Joined->True,PlotRange->All,Ticks->{Table[{8i +17,i},{i,-2,2}],Automatic},PlotLabel->{"cr",-3},DisplayFunction->Identity];
Show[{pc,pcr}]//Print;
For[j=-4,j>=-7,j--,
For[n=-2^(-j),n≤2^(-j),n++,cr[[j+8,2n+2^(-j+1)+1]]=Sum[(cr[[j+9,n+2^(-j)+1-i]]hhh[[2i+1]]+d[[j+9,n+2^(-j)+1-i]]ggg[[2i+1]]),{i,0,Min[n+2^(-j),6]}];
cr[[j+8,2n+2^(-j+1)+2]]=Sum[(cr[[j+9,n+2^(-j)+1-i]]hhh[[2i+2]]+d[[j+9,n+2^(-j)+1-i]]ggg[[2i+2]]),{i,0,Min[n+2^(-j),6]}]
];
pc=ListPlot[Table[c[[j+8,k+2^(-j+1)+1]]*2^(-j/2),{k,-2^(-j+1),2^(-j+1)}],PlotStyle->Red,Joined->True,PlotRange->All,Ticks->{Table[{2^(-j)i +2^(-j+1)+1,i},{i,-2,2}],Automatic},PlotLabel->{"c",j},DisplayFunction->Identity];
pcr=ListPlot[Table[cr[[j+8,k+2^(-j+1)+1]]*2^(-j/2),{k,-2^(-j+1),2^(-j+1)}],PlotStyle->Blue,Joined->True,PlotRange->All,Ticks->{Table[{2^(-j)i +2^(-j+1)+1,i},{i,-2,2}],Automatic},PlotLabel->{"cr",j},DisplayFunction->Identity];Show[{pc,pcr}]//Print
]

Cascade algorithm for computing Daubechies wavelets

Reference: §6.5 in I. Daubechies, Ten lectures on wavelets, Regional Conference Series in Applied Math. 61, SIAM, 1992

nn=8;

Initialize the arrays cr and dr:

cr=Table[0,{j,1,nn},{k,0,2^(nn+4)}];
dr=Table[0,{j,1,nn},{k,0,2^(nn+4)}];

cr will generate the wavelet φ at resolution 2^(-8) from an initial value and using the wavelet coefficients hhh[[i]] (red graphs).
dr will generate the wavelet ψ at resolution 2^(-8) from an initial value and using the wavelet coefficients ggg[[i]] in the first step and next the wavelet coefficients hhh[[i]] (blue graphs).

For[i=0,i≤13,i++,cr[[nn,i+1]]=hhh[[i+1]];dr[[nn,i+1]]=ggg[[i+1]]
];
pc=ListPlot[Table[cr[[nn,k+1]]*2^(1/2),{k,0,26}],PlotStyle->Red,Joined->True,Ticks->{Table[{2i +1,i},{i,0,13}],Automatic},PlotRange->All,PlotLabel->-1,DisplayFunction->Identity];
pd=ListPlot[Table[dr[[nn,k+1]]*2^(1/2),{k,0,26}],PlotStyle->Blue,Joined->True,Ticks->{Table[{2i +1,i},{i,0,13}],Automatic},PlotRange->All,PlotLabel->-1,DisplayFunction->Identity];
Show[{pc,pd}]//Print
For[j=-2,j>=-nn,j--,
For[n=0,n≤2^(-j+3)-1,n++,cr[[j+nn+1,2n+1]]=Sum[cr[[j+nn+2,n-i+1]]hhh[[2i+1]],{i,0,Min[n,6]}];
cr[[j+nn+1,2n+2]]=Sum[cr[[j+nn+2,n-i+1]]hhh[[2i+2]],{i,0,Min[n,6]}];
dr[[j+nn+1,2n+1]]=Sum[dr[[j+nn+2,n-i+1]]hhh[[2i+1]],{i,0,Min[n,6]}];
dr[[j+nn+1,2n+2]]=Sum[dr[[j+nn+2,n-i+1]]hhh[[2i+2]],{i,0,Min[n,6]}]
];
pc=ListPlot[Table[cr[[j+nn+1,k+1]]*2^(-j/2),{k,0,13*2^(-j)}],PlotStyle->Red,Joined->True,Ticks->{Table[{2^(-j)i +1,i},{i,0,13}],Automatic},PlotRange->All,PlotLabel->j,DisplayFunction->Identity];
pd=ListPlot[Table[dr[[j+nn+1,k+1]]*2^(-j/2),{k,0,13*2^(-j)}],PlotStyle->Blue,Joined->True,Ticks->{Table[{2^(-j)i +1,i},{i,0,13}],Automatic},PlotRange->All,PlotLabel->j,DisplayFunction->Identity];
Show[{pc,pd}]//Print
]

The above graphs stabilize as resolution becomes bigger. Compare with the graph of (Subscript[, 7])φ and (Subscript[, 7])ψ on p.197 of Daubechies.

We can also compute ϕ(x)  and ψ(x) by their  expressions in terms of the ϕ(2x-i) and the wavelet coefficients hhh[[i]] resp. ggg[[i]]:

ccr=Table[0,{k,0,2^(nn+4)}];
ddr=Table[0,{k,0,2^(nn+4)}];
For[n=0,n≤2^(nn+4)-1,n++,
ccr[[n+1]]=Sum[cr[[2,n-i*2^(nn-1)+1]]hhh[[i+1]],{i,0,Min[n*2^(-nn+1),13]}];ddr[[n+1]]=Sum[cr[[2,n-i*2^(nn-1)+1]]ggg[[i+1]],{i,0,Min[n*2^(-nn+1),13]}]
];
pc=ListPlot[Table[ccr[[k+1]]*2^(nn/2),{k,0,13*2^nn}],PlotStyle->Red,Joined->True,Ticks->{Table[{2^nn i +1,i},{i,0,13}],Automatic},PlotRange->All,DisplayFunction->Identity];
pd=ListPlot[Table[ddr[[k+1]]*2^(nn/2),{k,0,13*2^nn}],PlotStyle->Blue,Joined->True,Ticks->{Table[{2^nn i +1,i},{i,0,13}],Automatic},PlotRange->All,DisplayFunction->Identity];
Show[{pc,pd}]//Print

Numerical results for the successive approcimations of ϕ:

For[i=1,i≤12,i++,Print[i,", ",cr[[3,i*2^(nn-2)+1]]*2^((nn-2)/2),", ",cr[[2,i*2^(nn-1)+1]]*2^((nn-1)/2),", ",cr[[1,i*2^nn+1]]*2^(nn/2)]]