[LO_D,HI_D,LO_R,HI_R] = WFILTERS('db2')
h=LO_R
g=HI_R
%h(1)=h(1)+0.01
%h(2)=h(2)-0.01
iterations=8

% [hf,ht] = pf(h,kk) computes and plots hf, the Fourier transform
%  of the scaling function  phi(t)  using the freq domain
%  infinite product formulation with kk iterations from the scaling
%  function coefficients h.  Also calculates and plots ht = phi(t)
%  using the inverse FFT               csb 5/19/93

L = 2^(iterations+3); P = L;     % Sets number of sample points
omega=0:2*3.14159/(L-1):2*3.14159;
hp = fft(h,L); hf = hp;          % Initializes iteration
plot(omega,abs(hf));
pause
for k = 1:iterations-1                    % Iterations
  hp = [hp(1:2:L), hp(1:2:L)];   % Sample
  hf = hf.*hp/sqrt(2);           % Product
  plot(omega,abs(hf(1:L)));pause;    % Plot Phi(omega) each iteration
     P=P/2;                      % Scales axis for plot
end;

%init paralell iterations for HighPart
gp = fft(g,L);
for k = 1:iterations                    % Iterations
  %paralell iterations for HighPart
  gp = [gp(1:2:L), gp(1:2:L)];   % Sample   
end
hf = hf.*gp/sqrt(2);           % Product

ht = real(ifft(hf));             % psi(t) from inverse FFT
ht = ht(1:8*2^iterations); 

plot(ht(1:8*2^iterations));   % Plot psi(t)