function [phi,psi,alpha]=wave_cascade(g,h)
%
% A.R. Butz, March 2012
%
% Iterates the cascade algoritm, plotting the resulting scaling  function
% at each iteration corresponding to low pass filter coefficients g,
% asking whether or not you want to continue.
%
% h is the high pass filter coefficients
%
% When you are done phi is the scaling function and psi is the wavelet.
% alpha is the scale of the independent variable, its length being the
% max of the lengths of phi and psi. These can be plotted.
%
fac=2/sum(g);
c=fac*g;
rake=sign(cos(pi*[1:length(h)]));
[m,n]=size(h);
if n < m
    rake=rake';
end
fad=2/abs(sum(h.*rake));
d=fad*h;
phi=c;
psi=d;
iter=1;
ans=1;
while (ans ~= 0)
    phi=upsample(phi,2);
    phi=conv(c,phi);
    psi=upsample(psi,2);
    psi=conv(d,phi);
    deltaalpha=1/2^(iter+1);
    L=max(length(phi),length(psi));
    alpha=deltaalpha*[1:L];
    plot(alpha(1:length(phi)),phi);
    fprintf('iteration %d\n',iter)
    iter=iter+1;
    ans=input('Continue ( Yes(1) or No(0) )? ');
end
phi=[phi,zeros(1,length(alpha)-length(phi))];
psi=[psi,zeros(1,length(alpha)-length(psi))];
return

function xup=upsample(x,L)
[m,n]=size(x);
if min(m,n) ~= 1
    fprintf('\nInput vector to upsample should be row or column matrix.\n')
    return
end
len=max(m,n);
xup=zeros(1,L*len);
if n == 1
    xup=xup';
end
xup(1:L:L*len)=x;
return