% function h = daub(Nh) % % Generate filter coefficients for the Daubechies orthogonal wavelets. % % Kevin Amaratunga % 9 December, 1994. % % 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 = (-z + 2 - z.^(-1)) / 4; % sin^2(w/2) Mz1 = zeros(1,N); vec = ones(1,N); for l = 0:K-1 % Mz1 = Mz1 + binomial(L+l-1,l) * tmp2.^l; Mz1 = Mz1 + vec; vec = vec .* tmp2 * (L + l) / (l + 1); 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)';