function w = r_dpwt(f, l, h, level)
% FUNCTION:     r_dpwt			part of DPWT Toolbox
%
% SYNOPSIS:     W = r_dpwt(F, L, H, LEVEL)
%
% DESCRIPTION:  This is the recursive portion of the
%               discrete periodic wavelet transform using 
%               filters L and H, which are assumed to be  
%               a perfect reconstruction filter pair of even 
%               length. Decomposition is performed to level
%               LEVEL. For a description of the functionality 
%               of dpwt see the documentation for DPWT or 
%               enter "help dpwt".
%
%               The returned row vector W is W = [f0, w0, w1, ..., wp]
%               where f0 is length 1 and wp is a length 2^p row. These
%               are the wavelet coordinates of F.
%
%               The length of F must be an integer power of 2. 
%               Filter sequence L is not checked for satisfaction of 
%               perfect reconstruction criteria. It must be of even length.
%
% REFERENCE:    N. Getz, "A Fast Discrete Periodic Wavelet Transform",
%               Electronics Research Laboratory, U.C. Berkeley, 1992.
%

% AUTHOR:       Neil Getz
% DATE:         12-11-92
% 

lenf = length(f);

if (lenf == 2^level),
  % Then we've hit bottom.
  w = f;
else,
  % Haven't hit bottom yet.
  MATSHIFT = 1;  % matlab indices start at 1, not 0.
  
  % Initialization and local constants
  %
  lenfpm1 = lenf/2;	% Saves a couple of divisions.
  fpm1 = zeros(1,lenfpm1);
  wpm1 = fpm1;	% For size
  Nt2 = length(l);
  
  % Make DPWT filters l_tilde and h_tilde
  %
  [l_tilde,h_tilde] = lh_tilde(l, h, lenf);
  %
  
  % The DPWT calculations.
  %
  for i = 0:lenfpm1 - 1,
    two_i = 2*i;  % Saves a couple of multiplications.
    for n = mod([two_i:(two_i + min(lenf,Nt2) - 1)],lenf),
      dexmod = mod(n - two_i, lenf);
      fpm1(MATSHIFT + i) = fpm1(MATSHIFT + i) ...
                           + l_tilde(MATSHIFT + dexmod)*f(MATSHIFT + n);
      wpm1(MATSHIFT + i) = wpm1(MATSHIFT + i) ...
                           + h_tilde(MATSHIFT + dexmod)*f(MATSHIFT + n);
    end
  end
  %

  % The recursion
  %
  w = [r_dpwt(fpm1,l,h,level), wpm1];
  %

end

% END of r_dpwt().