function fp = r_idpwt(w, l, h, level)
% FUNCTION:     r_idpwt
%
% SYNOPSIS:     FP = r_idpwt(W, L, H, LEVEL)
%
% DESCRIPTION:  This is the recursive portion of the inverse
%               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".
%               
%               Wavelet coefficient vector W has structure,
%               W  = [flevel, wlevel, ...w(p-1), wp]
%               where flevel and wlevel are the row vectors
%               that are the level LEVEL decomposition 
%               products of the DPWT applied to FP. The length
%               of W 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
%
% ORGANIZATION: University of California
%               Department of Electrical Engineering
%               and the Electronics Research Laboratory 
%               Cory Hall
%               Berkeley, CA 94720
%
% DATE: 12-11-92
% 

lenw = length(w); % lenw is 2^p;

if(lenw == 2^level)
  % Then we've hit bottom.
  fp = w;
else,
  % Haven't hit bottom yet.
  MATSHIFT = 1;	% Matlab indices start at 1, not 0. 
  
  % Initialization
  %
  fp = zeros(1,lenw);	% Place to put result.
  halflenw = lenw/2;	% Saves a couple of divisions. This is 2^(p-1)
  N = length(l)/2;
  
  %
  % DPWT filters l_tilde and h_tilde
  %
  [l_tilde,h_tilde] = lh_tilde(l, h, lenw);
  %
  
  % The recursion.
  %
  fpm1 = r_idpwt( ...
           w( (MATSHIFT + 0):(MATSHIFT + halflenw - 1) ), l, h, level ...
         );
  %
  
  % The IDPWT calculations
  %
  for i = 0:lenw - 1,
    firstcol = ceil((i+1)/2)-1;
    for k = mod([firstcol:-1:(firstcol - min([halflenw, N]) + 1)],halflenw),
      dexmod = mod(i-2*k,lenw);
      fp(MATSHIFT + i) = fp(MATSHIFT + i)                ...
                         + l_tilde(MATSHIFT + dexmod)    ...
                           * fpm1(MATSHIFT + k)          ...
                         + h_tilde(MATSHIFT + dexmod)    ...
                           * w(MATSHIFT + halflenw + k);
    end
  end
  %
end

% END of r_idpwt().