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 4


Software: A fast matlab routine for calculating Daubechies filters
 
images/spacer.gifimages/spacer.gif Reply into Digest
Previous :: Next  
Author Message
Kevin Amaratunga
Guest





PostPosted: Mon Dec 02, 2002 6:07 pm    
Subject: Software: A fast matlab routine for calculating Daubechies filters
Reply with quote

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.122455722880489e-23
sum(h(k)*h(k-2l)) = 1.871421275103674e-13 (Max over l, l != 0)


Note that cepstral aliasing due to the FFT is not a significant problem.

References:
Oppenheim and Schafer - Discrete-Time 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:N-1;

% 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:K-1
Mz1 = Mz1 + binomial(L+l-1,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 1-250 77 Massachusetts Avenue Cambridge, MA 02139.
(617) 253 7657
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.024604 seconds : 18 queries executed : GZIP compression disabled