Return to the homepage
Search the complete Wavelet Digest database
Help about the Wavelet Digest mailing list
About the Wavelet Digest
The Wavelet Digest
> Volume 4, Issue 4
Software: A fast matlab routine for calculating Daubechies filters

Previous :: Next

Author 
Message 
Kevin Amaratunga Guest

Posted: Mon Dec 02, 2002 6:07 pm Subject: Software: A fast matlab routine for calculating Daubechies filters




Software: A fast matlab routine for calculating Daubechies filters
Here is a matlab routine for computing Daubechies filter coefficients
which does the spectral factorization using cepstral analysis. This
method works well for long filter lengths and it tends to be faster
than methods which make use of root finding algorithms. With the
approximately 5.5 CPU seconds on a DECstation 5000/25. A quick check
shows that
sum(h(k)) = 2.00000000000044
sum(h(k)*g(k)) = 1.122455722880489e23
sum(h(k)*h(k2l)) = 1.871421275103674e13 (Max over l, l != 0)
Note that cepstral aliasing due to the FFT is not a significant problem.
References:
Oppenheim and Schafer  DiscreteTime Signal Processing.

% function h = daub(Nh)
%
% Generate filter coefficients for the Daubechies orthogonal wavelets.
%
% h = filter coefficients of Daubechies orthonormal compactly supported
% wavelets
% Nh = length of filter.
function h = daub(Nh)
K = Nh/2;
L = Nh/2;
N = 512; % Use a 512 point FFT by default.
k = 0:N1;
% Determine samples of the z transform of Mz (= Mz1 Mz2) on the unit circle.
% Mz2 = z.^L .* ((1 + z.^(1)) / 2).^(2*L);
z = exp(j*2*pi*k/N);
tmp1 = (1 + z.^(1)) / 2;
tmp2 = (1  z.^(1)) / 2;
Mz1 = zeros(1,N);
for l = 0:K1
Mz1 = Mz1 + binomial(L+l1,l) * (z).^l .* tmp2.^(2*l);
end
Mz1 = 4 * Mz1;
% Mz1 has no zeros on the unit circle, so use the complex cepstrum to find
% its minimum phase spectral factor.
Mz1hat = log(Mz1);
m1hat = ifft(Mz1hat); % Real cepstrum of ifft(Mz1). (= cmplx
% cepstrum since Mz1 real, +ve.)
m1hat(N/2+1:N) = zeros(1,N/2); % Retain just the causal part.
m1hat(1) = m1hat(1) / 2; % Value at zero is shared between
% the causal and anticausal part.
G = exp(fft(m1hat,N)); % Min phase spectral factor of Mz1.
% Mz2 has zeros on the unit circle, but its minimum phase spectral factor
% is just tmp1.^L.
Hz = G .* tmp1.^L;
h = real(ifft(Hz));
h = h(1:Nh)';
John R. Williams and Kevin Amaratunga
Intelligent Engineering Systems Group Department of Civil Engineering
MIT Room 1250 77 Massachusetts Avenue Cambridge, MA 02139.
(617) 253 7657 





All times are GMT + 1 Hour

Page 1 of 1 
