function [y,z]=AnalFilterBank(x,M,Fs)

% ANALFILTBANK performs analyzes input signal x to M bands
% each one corresponding to a frequency zone of Fs/(2*M) Hz
%   y=AnalFilterBank(x,M,Fs), x is the input signal, M is number of bands, 
%   M should be of length 2^n, where n is a positive integer (n>=1)
%   Fs is the sampling frequency and y is the output. 
%   The first row (samples K = length(x)/M) of output y correspond  
%   to band #1, the second row of y corresponds to band #2, etc.  
%
L=length(x);

if rem(L,M),
    L=ceil(L/M)*M;
end
w=round(L/M);
y = zeros(L,1);
y(1:length(x))=x;
x=y; 
y = reshape(y,M,w);

% Design the lowpass FIR filter
% In order for a cosine modulated filterbank to be designed the filter h 
% should be:
% (a) FIR
% (b) of order N = 2*m*M, where M is the number of filters in the
% filterbank and m is a positive integer (m>=1)
N = 2*2*M;
h = fir1(N,1/(2*M));
t=[0:N];
h1 = 2*h.*cos((2*1-1)*t(1:N+1)*pi/(2*M));
H=abs(fft(h1,1024));
mx=max(H);

z = zeros(1,L);

% Show input spectrum
%figure(1)
%Y1 = fft(x,1024);
%bar([0:511]*Fs/1024, abs(Y1(1:512))/1024,'r');
%axis([0 Fs/2 0 0.025])
%grid on

% zero pad at the end to account for delays due to filtering
x = [x; zeros(N/2,1)];

for k=1:M,
     hk = 2*h.*cos((2*k-1)*t*pi/(2*M))/mx;
     y1 = filter(hk,1,x);
     
     % Crop the delay
     y1 = y1(N/2+1:end);
     y(k,:) = y1(k:M:L);

    % Reconstruction
     y1 = zeros(L+N/2,1);
     y1(k:M:end-N/2) = M*y(k,:);
     y1 = filter(mx*hk,1,y1);
     y1 = y1(N/2+1:end);
     z = z + y1';
end
%figure(2)
%Y1 = fft(z,1024);
%bar([0:511]*Fs/1024, abs(Y1(1:512))/1024,'r');
%axis([0 Fs/2 0 0.025])
%grid on