function lo = dblei(W,L,h0,h1,h2)
%  Double-density DWT (Inverse Transform)
%
%  y = dblei(W,L,h0,h1,h2);
% 
%  See also dblef (forward transform)

N = length(h0);
h0 = h0(N:-1:1);
h1 = h1(N:-1:1);
h2 = h2(N:-1:1);

lo = W{L+1};
for j = L:-1:1
   lo = pupfir(lo,h0,1-N);
   lo = lo + pupfir(W{j}{1},h1,1-N);
   lo = lo + pupfir(W{j}{2},h2,1-N);
end


% --------------------------------------------------
%      subroutine
% --------------------------------------------------


function y = pupfir(x,h,C)
% y = pupfir(x,h,C)
%
% Periodic up-sampling and FIR filtering.
%
% x: input
% h: fir filter
% C: index of first value of h (default C = 0)
% (The support of h is C:C+Nh-1)
% (This parameter is useful for alignment)
%
% It is assumed that the support of x is 0:Nx-1.
%
% % Example
% x = [1 2 5 3 1 2 4 5];
% h = [1 2 3 2];
% C = 1;
% y = pupfir(x,h,C);

if nargin < 3
   C = 0;
end

   Nx = length(x);
   Nh = length(h);
   Ny = 2*Nx;

   % upsample and filter
   x = upfirdn(x,h,2,1);

   % periodize at boundary (add tail to front)
   x(1:Nh-2) = x(1:Nh-2) + x(Ny+1:Ny+Nh-2);

   % alignment
   C = mod(C,Ny);
   y = zeros(1,Ny);
   y(1:C) = x(Ny-C+1:Ny);
   y(C+1:Ny) = x(1:Ny-C);