Discrete Wavelets with Applications

071 - Solving for the Coefficients of Daubechies Filters

Helmut  Knaust - Department of Mathematical Sciences, University of Texas at El Paso, hknaust@utep.edu
10/26/2010; last edit 4/5/2022

SetDirectory[NotebookDirectory[]];
<<DW.m

D-4: The Orthogonality Conditions for the Low Pass Filter

Our low pass filter will have the form [h0, h1, h2, h3] with corresponding matrix

LP[8]=SparseArray[{{i_,j_}/;j==2i-1->h0,{i_,j_}/;j==2i->h1,{i_,j_}/;j==2i+1->h2,{i_,j_}/;j==2i+2->h3,{4,1}->h2,{4,2}->h3},{4,8}]//Normal

Note the wrapping of the last row to accommodate the length 4 filter.

To obtain orthogonality we need (1) that the scalar product of the first row with itself is 1, and that (2) the scalar product of two consecutive rows is 0.

O1=h0^2+h1^2+h2^2+h3^2==1

O2=h0 h2+h1 h3==0

Here is a quick orthogonality check :

LP[8].Transpose[LP[8]]

LP[8].Transpose[LP[8]]/.{h0^2+h1^2+h2^2+h3^2->1,h0 h2+h1 h3->0}

D-4: The High Pass Filter

We will use an earlier orthogonality "trick" to choose the coefficients for the high pass filter :

HP[8]=SparseArray[{{i_,j_}/;j==2i-1->-h3,{i_,j_}/;j==2i->h2,{i_,j_}/;j==2i+1->-h1,{i_,j_}/;j==2i+2->h0,{4,1}->-h1,{4,2}->h0},{4,8}]//Normal

HP[8].Transpose[HP[8]]

Corresponding rows of the low and high pass filters thus are automatically orthogonal. Overlapping rows of the high and low pass filters are also orthogonal.

LP[8].Transpose[HP[8]]
HP[8].Transpose[LP[8]]

D-4: Filter Conditions for the Lowpass Filter

Here is the Fourier expansion of the LP filter  lp[t]::

Clear[t];
lp[t_]={h0,h1,h2,h3}.Table[Exp[k I t],{k,0,3}]

We need to impose the filter condition 
		lp[π] = 0  :

lp[π]
F1=lp[π]==0//FullSimplify

Recall the plot of the Haar filter :

FourierFilterPlot[{1/Sqrt[2],1/Sqrt[2]}]

To obtain a better filter we will impose the the new condition that the derivative of lp[t] vanishes at π :
   	lp'[π] = 0

lp'[π]==0

Dividing by i yields :

F2=1/I (lp'[π])==0//FullSimplify

D-4: Solving this System of 4 Equations

Here is our complete systemof equations :

{O1,O2,F1,F2}//TableForm

We introduce a new variable c by rewriting the second equation as
	h3/h0 = -h2/h1 = c

With this identification, the first equation becomes

Factor[O1/.{h3->c h0,h2->-c h1}]/(1+c^2)

So solving for (h0^2)+(h1^2) yields 

	(E1)         (h0^2) + (h1^2)=(1/(1+(c^2)))

Substituting c in the third equation yields

F1/.{h3->c h0,h2->-c h1}

Solving for h1 yields 
	(E2)	 h1=((1-c)/(1+c))h0

Solve[F1/.{h3->c h0,h2->-c h1},h1]

For a given c, the last two equations (E1) and (E2) represent a circle and line :

cc=0.8;
Plot[{Sqrt[1/(1+cc^2)-h0^2],-Sqrt[1/(1+cc^2)-h0^2],(1+cc)(1-cc)h0},{h0,-1,1},AxesLabel->{h0,h1},AspectRatio->Automatic]

We still have not looked at the consequences of equation (F2) :

F2

Solving for h1 yields

	(E3)      h1=-((3c)/(1+2c))h0

Solve[F2/.{h2->-c h1,h3-> c h0},h1]

Let' s solve (E2) and (E3) for c :

sol=Solve[-3c/(1+2c)==(1-c)/(1+c)]

Discarding the first solution, and using the second solution in (E2) we obtain :

E4= h1==((1-c)/(1+c))h0/.sol[[2]]

Rationalizing the denominator :

sol2=Solve[E4,h1]//FullSimplify

Substituting in (E1) yields :

(h0^2) + (h1^2)==(1/(1+(c^2)))/.sol2[[1]]

Solving for h0 yields :

Solve[(h0^2) + (h1^2)==(1/(1+(c^2)))/.sol2[[1]]/.sol[[2]],h0]

By computing the numerical value of the second solution, we can see that this coincides with the answer given in the book.

N[(1/(4 (Sqrt[2-(Sqrt[3])])))]

N[1/4/Sqrt[2] (1+Sqrt[3])]

Finally h1, h2, and h3 are then obtained by backsolving.

D-4: Solving this System of 4 Equations Following the Lazy Man's Rule

Here is our complete system of equations :

{O1,O2,F1,F2}//TableForm

Apply the Lazy Man' s Rule : Let' s solve the (linear) filter conditions for h2 and h3 first:

sol1=Solve[{F1,F2},{h2,h3}]

Next, we substitute h2 and h3 in the orthogonality conditions ...

sol2={O1,O2}/.sol1[[1]];TableForm[sol2]

... and simplify :

sol3=FullSimplify[sol2];
TableForm[sol3]

These two equations represent an ellipse, and two lines, respectively:

ContourPlot[{14 (h0^2)+6 (h1^2)==16 h0 h1+1,3 (h0^2)==(h1^2)},{h0,-2,2},{h1,-2,2},AspectRatio->Automatic]

Here are the equations of the two lines  - we select the one with the positive slope:

sol4=Solve[sol3[[2]],h1]
sol4[[2]]

Next we substitute h1 in the ellipse equation to obtain h0 ...

sol5=sol3[[1]]/.sol4[[2]]

... and finally solve for h0, once again selecting the positive solution.

sol6=Solve[sol5]
sol7=sol6[[2]]

The rest is backsolving :

sol8=sol4[[2]]/.sol7

sol9=sol1[[1]]/.sol8/.sol7

Here is the final filter, both in symbolic and numerical form :

h={h0,h1,h2,h3}/.sol7/.sol8/.sol9;
FullSimplify[h]
N[h]

FourierFilterPlot[h]

Daubechies-4 Filter

Clear[h];
sols=Solve[{
h[0]^2+h[1]^2+h[2]^2+h[3]^2==1,
h[0] h[2] + h[1] h[3]==0,
h[0]-h[1]+h[2]-h[3]==0,
h[1]-2h[2]+3h[3]==0
}]//FullSimplify

N[%]//TableForm

sols[[2]]

Daub4=Array[h,4,0]/.sols[[2]]

N[Daub4]

d4=FourierFilterPlot[Daub4]

Daubechies-6 Filter

eqns[n_]:=Module[{L=n-1,k=2},Clear[h];Join[
	{Sum[h[i]^2,{i,0,L}]==1},
	Table[Sum[h[i] h[i+2k],{i,0,L-2k}]==0,{k,1,(L-1)/2}],
	{Sum[(-1)^i h[i],{i,0,L}]==0},
	Table[Sum[(-1)^(i+1)i^k h[i],{i,1,L}]==0,{k,1,(L-1)/2}]	
	]]

eqns[6]//TableForm

Clear[sols];
sols=NSolve[eqns[6]];
sols=Table[Sort[sols[[i]]],{i,1,Length[sols]}];
sols//TableForm

Length[sols]

sol=N[sols[[1]]]

Daub6=sol[[All,2]]

d6=FourierFilterPlot[Daub6]

Daubechies-10 Filter

eqns[10]//TableForm

sols=Solve[eqns[10]];

Length[sols]

N[sols]//TableForm

Picking the real solutions :

sols=Select[sols,(#[[1,2]]==Conjugate[#[[1,2]]]&)];
sols=Select[N[sols],(Abs[#[[1,2]]]<10&)];
sols=Table[Sort[sols[[i]]],{i,1,Length[sols]}];
sols//TableForm

Length[sols]

Daub10=sols[[1,All,2]]//N

ListPlot[Transpose[{Table[i,{i,0,9}],Daub10}],PlotStyle->AbsolutePointSize[10]]

d10=FourierFilterPlot[Daub10]