[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
%internal initialisation

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                     % 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;
ht = real(ifft(hf));             % phi(t) from inverse FFT
ht = ht(1:8*2^iterations); 
plot(ht(1:8*2^iterations));   % Plot phi(t)