function y=cascade(x,l,h,options)
% Computes the scaling function or wavelet by the cascade algorithm.
%
% x       evaluation points
% l       recursion depth
% h       scaling filter
%
% options 'phi' scaling function 
%         'psi' wavelet
%
% Course in Wavelets, TU Chemnitz, 2012

y=zeros(size(x));
if(strcmp(options,'psi'))
  if(l==0)
    y((-0.5<=x)&(x<0))=1;
    y((0<=x)&(x<0.5))=-1;
  else
    g=h(end:-1:1);
    g=g.*((-1).^(0:length(h)-1));
    for n=0:length(h)-1
      y=y+g(n+1)*cascade(2*x-n,l-1,h,'phi');
    end;
    y=sqrt(2)*y;
  end;
else
  if(l==0)
    y((-0.5<=x)&(x<0.5))=1;
  else
    for n=0:length(h)-1
      y=y+h(n+1)*cascade(2*x-n,l-1,h,'phi');
    end;
    y=sqrt(2)*y;
  end;
end;
