% Example 6: Compute the samples of Daubechies scaling function and
% wavelet using the inverse DWT.

%p = 2                                     % Number of zeros at pi.
N = 2 * p - 1;                            % Support of the scaling function
numlevels = 5;                            % Number of iterations/levels.
M = 2^numlevels;
L = M * N;
f0 = daub(N+1) / 2;                       % Synthesis lowpass filter.
f1 = (-1).^[0:N]' .* flipud(f0);          % Synthesis highpass filter.

% For the scaling function, we need to compute the inverse DWT with a delta
% for the approximation coefficients.  (All detail coefficients are set
% to zero.)
y = upcoef('a',[1;0],f0,f1,numlevels);    % Inverse DWT.
phi = M * [0; y(1:L)];

% For the wavelet, we need to compute the inverse DWT with a delta for the
% detail coefficients.  (All approximation coefficients and all detail
% coefficients at finer scales are set to zero.)
y = upcoef('d',[1;0],f0,f1,numlevels);    % Inverse DWT.
w = M * [0; y(1:L)];

% Determine the time vector.
t = [0:L]' / M;

% Plot the results.
plot(t,phi,'-',t,w,'--')
legend('Scaling function','Wavelet')
title('Scaling function and wavelet by iteration of synthesis filter bank.')
xlabel('t')
pause

% Now compute the scaling function and wavelet by recursion.
% phivals (not part of the Matlab toolbox) does this.
[t1,phi1,w1] = phivals(daub(2*p),numlevels);

% Plot the results.
plot(t1,phi1,'-',t1,w1,'--')
legend('Scaling function','Wavelet')
title('Scaling function and wavelet by recursion.')
xlabel('t')
pause

% View the scaling functions side by side.
plot(t,phi,'-',t1,phi1,'--')
legend('Scaling function using iteration','Scaling function using recursion')
title('Comparison of the two methods (recursion is exact.)')
xlabel('t')