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