Daubechies
Daubechies
Daubechies
Algorithm in MATLAB
Brani Vidakovic
1
Table 1: Pollen parameterization for N = 2 (four-tap filters). [s = 2 2]
n hn for N = 2
0 (1 + cos sin )/s
1 (1 + cos + sin )/s
2 (1 cos + sin )/s
3 (1 cos sin )/s
Daubechies-Lagarias Algorithm
2
Table 2: Pollen parameterization for N = 3 (six-tap filters) [s = 2 2].
n hn for N = 3
0 (1 + cos 1 cos 2 cos 1 cos 2 + sin 1 cos 2 sin 1 sin 2
+ cos 1 sin 2 sin 1 sin 2 )/(2s)
1 (1 cos 1 + cos 2 cos 1 cos 2 + sin 1 + cos 2 sin 1 sin 2
cos 1 sin 2 sin 1 sin 2 )/(2s)
2 (1 + cos 1 cos 2 + cos 2 sin 1 cos 1 sin 2 + sin 1 sin 2 )/s
3 (1 + cos 1 cos 2 cos 2 sin 1 + cos 1 sin 2 + sin 1 sin 2 )/s
4 (1 cos 1 + cos 2 cos 1 cos 2 sin 1 cos 2 sin 1 + sin 2
+ cos 1 sin 2 sin 1 sin 2 )/(2s)
5 (1 + cos 1 cos 2 cos 1 cos 2 sin 1 + cos 2 sin 1 + sin 2
cos 1 sin 2 sin 1 sin 2 )/(2s)
Then
3
Pollen2 scaling function for phi1=0.6789 and phi2=0.12345 Pollen2 wavelet function for phi1=0.6789 and phi2=0.12345
2.5 1.5
1
2
0.5
1.5
0
1
0.5
0.5
1
0
1.5
0.5
2
1 2.5
0 1 2 3 4 5 6 2 1 0 1 2 3 4
(a) (b)
The convergence of ||Td1 Td2 Tdn Td1 Td2 Tdn+m || to zero, for
fixed m, is exponential and constructive, i.e., effective bounds, that decrease
exponentially to 0, can be established.
4
are given as:
1+ 3 3+ 3 1+ 3
0 0 0
4
4
4
T0 = 3 3 3+ 3 1+ 3 , and T1 = 1 3 3 3 3+ 3 .
4 4
4 2
4 4 4
1 3 3 3 1 3
0 4 4
0 0 4
The Daubechies and Lagarias algorithm gives only the values of the scal-
ing function. However, most of the evaluation needed involves the mother
wavelet. It turns out that another algorithm is unnecessary, due to Theorem
2.
If for some i the index i + 1 [2x] is negative or larger than N then the
corresponding component of u is zero.
Let the vector v be defined as
1 0 Y
v(x, n) = 1 Ti (4)
N
idyad({2x},n)
5
Then
Matlab Codes
2. This m-function gives low pass wavelet filter from the one parameter Pollen
family.
function h = MakePollen1(phi)
% One parameter Orthonormal Pollen Family
% y = MakePollen1(phi)
% Input: phi - angle in [0, 2 *pi)
% Output: h - low pass ON wavelet filter
% phi=pi/6 corresponds to Daubechies 4 tap ON filter
s = 2 * sqrt(2);
h(1) = (1 + cos(phi) - sin(phi))/s;
6
h(2) = (1 + cos(phi) + sin(phi))/s;
h(3) = (1 - cos(phi) + sin(phi))/s;
h(4) = (1 - cos(phi) - sin(phi))/s;
% B. Vidakovic, 12/8/2002
3. This m-function gives low pass wavelet filter from the two parameter Pollen
family. The spacial case 1 = 0.67890, and 2 = 0.12345 is plotted in Figure
1. The filter for this selection of parameters is: {-0.0153, 0.2469, 0.8404, 0.4675,
-0.1180, -0.0073}.
4. The functions Phijk.m and Psijk.m given below calculate the scaling
(wavelet) functions at the point with preassigned precision. Functions dec2bin.m,
t0.m, and t1.m are provided in each Phijk.m and Psijk.m for convenience.
7
function yy = Phijk(z, j, k, filter, n)
%--------------------------------------------------------------
% yy=Phijk(z, j, k, filter, n)
% Evaluation of the scaling function corresponding to an Orthogonal
% MRA by Daubechies-Lagarias Algorithm.
% inputs: z -- the argument
% j -- scale
% k -- shift
% filter -- ON finite wavelet filter, might be an
% output of WaveLabs: MakeONFilter
% n -- precision of approximation maesured by the number
% of Daubechies-Lagarias steps (default n=20)
%--------------------------------------------------------
% output: yy -- value of father wavelet (j,k) coresponding to
% filter at z.
%--------------------------------------------------------------
% Example of use:
% > xx = 0:0.01:6; yy=[];
% > for i=1:length(xx)
% > yy =[yy Phijk(x(i), 0, 1, MakeONFilter(Daubechies,4), 25)];
% > end
% > plot(x,yy)
%----------------------------------------------------------------
if (nargin == 4)
n=20;
end
daun=length(filter)/2;
N=length(filter)-1;
x=(2^j)*z-k;
if(x<=0|x>=N) yy=0;
else
int=floor(x);
dec=x-int;
dy=dec2bin(dec,n);
t0=t0(filter);
t1=t1(filter);
prod=eye(N);
for i=1:n
if dy(i)==1 prod=prod*t1;
else prod=prod*t0;
end
end
y=2^(j/2)*prod;
yyy = mean(y);
yy = yyy(int+1);
end
%--------------------functions needed----------------------------
%---------
function a = dec2bin(x,n)
a=[];
for i = 1:n
8
if(x <= 0.5) a=[a 0]; x=2*x;
else a=[a 1]; x=2*x-1;
end
end
%-----------
function t0 = t0(filter)
%
n = length(filter); nn = n - 1;
%
t0 = zeros(nn); for i = 1:nn
for j= 1:nn
if (2*i - j > 0 & 2*i - j <= n)
t0(i,j) = sqrt(2) * filter( 2*i - j );
end
end
end
%------------------
function t1 = t1(filter)
%
n = length(filter); nn = n - 1;
%
t1 = zeros(nn); for i = 1:nn
for j= 1:nn
if (2*i -j+1 > 0 & 2*i - j+1 <= n)
t1(i,j) = sqrt(2) * filter( 2*i - j+1 );
end
end
end
%---------------------- B. Vidakovic, 2002 --------------------
%--------------------------------------------------------------
9
% > plot(x,yy)
%----------------------------------------------------------------
if (nargin == 4) n=20;
end
N=length(filt)-1;
daun = (N+1)/2;
x=(2^j)*z-k;
if(x<=daun-N|x>=daun) yy=0;
%if(x<=0|x>=N) yy=0;
else
twox = 2 * x;
inti=floor(twox);
dec=twox-inti;
dy=dec2bin(dec,n);
t0=t0(filt);
t1=t1(filt);
prod=eye(N);
for i=1:n
if dy(i)==1
prod=prod*t1;
else prod=prod*t0;
end
end
uu=[];
for i=1:N
index = i + 1 - inti;
if ( index > 0 & index < N + 2 )
fi= (-1)^(-index)*filt(index);
else fi=0;
end
uu =[uu fi];
end
%----------------------------------------
v = 1/N * ones(1,N) * prod ;
yy=2^(j/2)* uu * v;
end
%--------------------functions needed----------------------------
%---------
function a = dec2bin(t,n)
a=[];
for i = 1:n
if(t <= 0.5) a=[a 0]; t=2*t;
else a=[a 1]; t=2*t-1;
end
end
%-----------
function t0 = t0(filt)
%
n = length(filt); nn = n - 1;
%
t0 = zeros(nn); for i = 1:nn
10
for j= 1:nn
if (2*i - j > 0 & 2*i - j <= n)
t0(i,j) = sqrt(2) * filt( 2*i - j );
end
end
end
%------------------
function t1 = t1(filt)
%
n = length(filt); nn = n - 1;
%
t1 = zeros(nn); for i = 1:nn
for j= 1:nn
if (2*i -j+1 > 0 & 2*i - j+1 <= n)
t1(i,j) = sqrt(2) * filt( 2*i - j+1 );
end
end
end
%--------------------- B. Vidakovic, 2002 --------------------
%-----------------------------------------------------------------
References
[1] Daubechies, I. and Lagarias, J. (1991). Two-scale difference equations I. Ex-
istence and global regularity of solutions, SIAM J. Math. Anal.,22 (5), 1388
1410.
[3] Wavelab802, Matlab toolbox for wavelet analysis. Donoho, D. et al., Stanford
University. http://www-stat.stanford.edu/~
wavelab/.
[4] Pollen, D. (1990). SU I (2, F [z, 1/z]) for F a subfield of C, J. Amer. Math.
Soc.,3, 611624.
11