function  [w, lowpass, L1, L2]  =  analysis(f, J, alpha, tau)
%
% signal is first prefiltered, and then a J-level wavelet decomposition
% is performed. 
% wavelet coefficients are finally complexified.
%
% get projection filters
M        = length(f);
[P1, P2] = projectionFilters(M + 1, alpha, tau);

% make length M + 1 by repeating the last sample
% the odd length constraint comes up following certain symmetry properties
% of the Fourier transform (see references)
fpad  = [f  f(M)];
f1pad = real( ifft(P1.*fft(fpad)) );
f2pad = real( ifft(P2.*fft(fpad)) );

% extract the first M samples
f1    = f1pad(1 : M);
f2    = f2pad(1 : M);

% store last sample for synthesis
L1  = f1pad(M + 1);
L2  = f2pad(M + 1);

% (complex) wavelet coefficients
w = cell(J,1);

% analysis filters for the 2 channels 
shift    = 1/2;
M        = length(f);

[H1, G1] = filters1D(M, alpha, tau, 0);
[H2, G2] = filters1D(M, alpha, tau + shift, 0);

H1 = conj(H1);
G1 = conj(G1);
H2 = conj(H2);
G2 = conj(G2);

% perform analysis
for depth = 1 : J
    
   % channel 1
   f1 = fft(f1,M);
   Y1 = H1 .* f1;
   Z1 = G1 .* f1;
   Y1 = 1/2 * ( Y1(1:M/2) + Y1(M/2+(1:M/2)) );
   Z1 = 1/2 * ( Z1(1:M/2) + Z1(M/2+(1:M/2)) );
   f1 = real( ifft(Y1,M/2) );
   
   % channel 2 
   f2 = fft(f2,M);
   Y2 = H2 .* f2;
   Z2 = G2 .* f2;
   Y2 = 1/2 * ( Y2(1:M/2) + Y2(M/2+(1:M/2)) );
   Z2 = 1/2 * ( Z2(1:M/2)  + Z2(M/2+(1:M/2)) );
   f2 = real( ifft(Y2,M/2) );
   
   % complexify
   w{depth} = real(ifft(Z1,M/2)) + sqrt(-1) * real(ifft(Z2,M/2));
 
   % decimate filters 
   H1 = H1(1 : 2 : length(H1));
   G1 = G1(1 : 2 : length(G1));    
   H2 = H2(1 : 2 : length(H2));  
   G2 = G2(1 : 2 : length(G2));
  
   M = M/2;
end

% store lowpass (coarsest) signals
lowpass=[f1 ; f2];




