% This function computes the compactly supported 
% scaling and wavelet functions based on the 
% impulse response of a lowpass wavelet fileter h^(z)
% using the cascade algorithm.
% Inputs:  h -- coefficents of the lowpass wavelet filter
%          g -- coefficients of the highpass wavelet filter
%       epsi -- tolerance (for convergence)
% Outputs: phi -- scaling function on [0, N-1]
%          psi -- wavelet function on [-N/2+1, N/2]
% Example: [h,g,rh,rg] = daub(4); [phi,psi] = cascade(rh,rg,1e-3);
% Written by W.-S. Lu, University of Victoria
% Last modified: Feb. 18, 1999.

function [phi,psi] = cascade(h,g,epsi)

% Initialize iteration using the Haar scaling function

N = length(h);
phw = [ones(1,1024) zeros(1,(N-2)*1024)];
e = 1;
% Compute scaling function using the cascade iteration algorithm

while e > epsi,
   phw20 = [phw(1:2:(N-1)*1024-1) zeros(1,(N-1)*512)];
   phw2 = h(1)*phw20;
   for k = 1:N-1,
      w =  h(k+1)*[zeros(1,k*512) phw20(1:(N-1)*1024-k*512)];
      phw2 = phw2 + w;
   end
   phw2 = sqrt(2)*phw2;
   e = norm(phw2-phw)/N;
   phw = phw2;
end
phi = phw;

% Compute wavelet function using the wavelet equation

psw20 = [phw(1:2:(N-1)*1024-1) zeros(1,(N-1)*512)];
psw2 = g(1)*psw20;
for k = 1:N-1,
  w = g(k+1)*[zeros(1,k*512) phw20(1:(N-1)*1024-k*512)];
  psw2 = psw2 + w;
end
psi = sqrt(2)*psw2;

% Plot the scaling and wavelet functions computed

t = 0:(N-1)/((N-1)*1024-1):N-1;
t1 = -N/2+1:(N-1)/((N-1)*1024-1):N/2;
figure(1)
subplot(121)
plot(t,phi)
axis([0 N-1 1.1*min(phi) 1.1*max(phi)])
axis('square')
grid
title('scaling function')
xlabel('time')
subplot(122)
plot(t1,psi)
axis([-N/2+1 N/2 1.1*min(psi) 1.1*max(psi)])
axis('square')
grid
title('wavelet function')
xlabel('time')