% 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)';