function W = dblef(lo,L,h0,h1,h2)
%  Double-Density DWT (Forward transform)
%
%  W = dblef(x,L,h0,h1,h2)
%  x : data
%  L : number of scales
%  h0, h1, h2 : filters
%  W{scale}{subband} : transform coefficients
% 
%  See also dblei (inverse transform)
%
%  % Example:
%  K422                       % design filters
%  L = 3;                     % number of scales
%  N = 2^7;                   % data length
%  x = rand(1,N);             % data
%  W = dblef(x,L,h0,h1,h2);   % forward transform
%  y = dblei(W,L,h0,h1,h2);   % inverse transform
%  max(abs(x-y))              % check perferct reconstruction


% STAGE 1:L
for j = 1:L
   W{j}{1} = pfirdn(lo,h1);
   W{j}{2} = pfirdn(lo,h2);
   lo = pfirdn(lo,h0);
end
W{L+1} = lo;


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


function y = pfirdn(x,h,C)
% y = pfirdn(x,h,C)
%
% Periodic FIR filtering and down-sampling.
%
% 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.
% It is assumed that x is of even length.
%
% % Example
% x = [1 2 5 3 1 2 4 5];
% h = [1 2 3 2];
% C = 1;
% y = pfirdn(x,h,C);

if nargin < 3
   C = 0;
end

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

   % shift
   C = mod(C,Nx);
   g = zeros(1,Nx);
   g(1:C) = x(Nx-C+1:Nx);
   g(C+1:Nx) = x(1:Nx-C);

   % filter and down-sample
   y = upfirdn(g,h,1,2);

   % periodize at boundary (add tail to front)
   Nt = length(y)-Ny; % (length of tail)
   y(1:Nt) = y(1:Nt) + y(Ny+1:Ny+Nt);
   y = y(1:Ny);