The Wavelet Digest Homepage
Return to the homepage
Search the complete Wavelet Digest database
Help about the Wavelet Digest mailing list
About the Wavelet Digest
The Digest The Community
 Latest Issue  Back Issues  Events  Gallery
The Wavelet Digest
   -> Volume 4, Issue 3


Software: Accurate Daubechies filter coefficients in Mathematica
 
images/spacer.gifimages/spacer.gif Reply into Digest
Previous :: Next  
Author Message
paul@earwax.pd.uwa.edu.au (Paul Abbott)
Guest





PostPosted: Mon Dec 02, 2002 5:45 pm    
Subject: Software: Accurate Daubechies filter coefficients in Mathematica
Reply with quote

Software: Accurate Daubechies filter coefficients in Mathematica

Accurate computation of the Daubechies filter coefficients using Riesz
factorisation (as described in Chapter 6 of Daubechies' Ten Lectures on
Wavelets) is straightforward to implement in Mathematica:

P[n_, y_] := Sum[Binomial[n-1+k,k] y^k, {k,0,n-1}]

Filter[n_, prec_:20] :=
Module[{poly, r, small, z, l},
poly = P[n/2,1/2 - 1/(4z) - z/4];
r = z /. N[Solve[poly == 0, z], prec];
small = Select[r, Abs[#] < 1 &];
poly = (z+1)^(n/2) Times @@ (z - small);
l = CoefficientList[poly,z] // Reverse;
l / Sqrt[l.l]
] /; EvenQ[n]

The precision, prec, is the working precision not the final precision of
the resulting filter coefficients. For example, here are the D30
coefficients computed using 40 decimals working precision:

(Daubechies[30] = Filter[30, 40]) // Timing

{15.5333 Second,
{0.00453853736157889888145939491, 0.04674339489276627189170969335,
0.2060238639869957315398915009, 0.4926317717081396236067757074,
0.645813140357424358176420912, 0.3390025354547315276912641144,
-0.1932041396091454287063990534, -0.288882596566965646248412501,
0.0652829528487728169228310792, 0.1901467140071229823484893117,
-0.039666176555790944483843668, -0.111120936037231693365671032,
0.033877143923507686208548178, 0.054780550584507612689137903,
-0.025767007328439962585945258, -0.0208100501696930816778848342,
0.0150839180278359023632927446, 0.0051010003604075431697088602,
-0.0064877345603157449951816831, -0.00024175649076162428116672253,
0.00194332398038221154176491233, -0.000373482354137616992009809421,
-0.000359565244362468812164962008, 0.0001558964899205997479471658241,

0.00002579269915531893680925862418,-0.00002813329626604781364755324777,
-6 -6
3.36298718173757980312484521 10 , 1.811270407940577083768510912
10 ,
-7 -8
-6.3168823258816644212015973 10 , 6.13335991330575202905629946 10 }}

It is easy to check how well the conditions on the filter coefficients are
satisfied. Checking the Normalisation:

(Plus @@ Daubechies[30]) - N[Sqrt[2], 50]

-33
0. 10

Daubechies[30].Daubechies[30] - 1

-33
0. 10

Checking the Approximation Conditions:

Daubechies[30] . Table[(-1)^i,{i,30}]

-27
0. 10

Table[Daubechies[30] . Table[i^l (-1)^i,{i,n}], {l,0,30/2-1}]

-27 -39 -37 -36 -35 -34
{0. 10 , 0. 10 , 0. 10 , 0. 10 , 0. 10 , 0. 10 ,
-32 -30 -21 -21 -21 -21
0. 10 , 0. 10 , 0. 10 , 0. 10 , 0. 10 , 0. 10 ,
-22 -22 -20
0. 10 , 0. 10 , 0. 10 }

It is clear that the highest moment condition leads to the greates
cancellation error.

Checking the Orthogonalisation Conditions:

Table[Take[Daubechies[30],k].Take[Daubechies[30],-k],{k,2,30,2}]

-47 -40 -40 -40 -40 -40
{0. 10 , 0. 10 , 0. 10 , 0. 10 , 0. 10 , 0. 10 ,
-39 -39 -38 -38 -38 -38
0. 10 , 0. 10 , 0. 10 , 0. 10 , 0. 10 , 0. 10 ,
-37 -38
0. 10 , 0. 10 , 1.}

Paul Abbott
Department of Physics Phone: +61-9-380-2734
The University of Western Australia Fax: +61-9-380-1014
Nedlands WA 6907 email: paul@earwax.pd.uwa.edu.au
AUSTRALIA
www: http://www.pd.uwa.edu.au/Paul/abbott.html

----
Note from the editor: Using this software, Chris Heil was able to
find new and confirm known typos in the tables in Ingrid Daubechies'
"Ten Lectures". I have summarised these below.

p. 195, N = 4, n = 2, 11th digit should be 2 i.s.o. 3
p. 196, N = 8, 5th coefficient, first digit should be 9 i.s.o. blank
p. 277, N~ = 3, N = 7, coeff of z^{-3} is 363 i.s.o. 336
All times are GMT + 1 Hour
Page 1 of 1

 
Jump to: 
 


disclaimer - webmaster@wavelet.org
Powered by phpBB

This page was created in 0.027893 seconds : 18 queries executed : GZIP compression disabled