% function a = daub(N)
%
% Generate _lter coe_cients for the Daubechies orthogonal wavelets.
% a = _lter coe_cients of Daubechies' orthonormal compactly supported
% wavelets.
% Na = length of _lter.

% Hierarchical Wavelet-Based Numerical Models
% for Digital Data
% by
% Kevin S. Amaratunga
% B.Eng. (Hons. I), University of Southampton, UK (1991)
% S.M., Massachusetts Institute of Technology (1993)
% Submitted to the Department of Civil and Environmental
% Engineering
% in partial fulfillment of the requirements for the degree of
% Doctor of Philosophy in Computational Engineering
% at the
% MASSACHUSETTS INSTITUTE OF TECHNOLOGY
% June 1996
% c  Massachusetts Institute of Technology 1996. All rights reserved.

function a = hdaub(ord)
K = ord;
L = ord;
N = 512; % Use a 512 point FFT by default.
k = 0:N-1;
% Determine samples of the z transform of M1(z) = R(z) R(1/z) on the unit circle.
z = exp(j*2*pi*k/N);
tmp1 = (1 + z.^(-1)) / 2;
tmp2 = (-z + 2 - z.^(-1)) / 4; % sin^2(w=2).
M1z = zeros(1,N);
vec = ones(1,N);
for l = 0:K-1
    M1z = M1z + vec;
    vec = vec .* tmp2 * (L + l) / (l + 1);
end
M1z = 4 * M1z;
% M1(z) has no zeros on the unit circle, so use the complex cepstrum to _nd
% its minimum phase spectral factor.

M1zhat = log(M1z);
m1hat = ifft(M1zhat); % Real cepstrum of i_t(M1(z)).
% (= complex cepstrum since M1(z) 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.
Rz = exp(fft(m1hat)); % Min phase spectral factor of M1(z).
Az = Rz .* tmp1.^L;
a = real(ifft(Az));
a = 0.5 * a(1:(2*ord))';