%Uniform Filter Bank
%Written By: Ebrahim Soroush
% Dept. of Electrical Engineering
% Amirkabir of University Technology
%the bandwidths,as defined,are 4pi/N the filter pass-bands overlap
%by 50% A superposition of the bandpass frequency responses for N
%Also shown is the frequency-response sum, which we will show to
%be exactly constant and equal to N This gives our filter bank the
%perfect reconstruction property. We can simply add the outputs of
%the filters in the filter bank to recreate our input signal exactly.
%This is the source of the name Filter-Bank Summation (FBS)
% I used direct DTFT formula for implemntig it, but first of all
% for single filter of MA with N=15 [-7,7] plotted with FFT func
% for clarifying our figure
clc,clear all,close all
N=15;N2=(N-1)/2;k=0:N-1;Wk=2*k*pi/N;%define N and N/2 and wk shift key
n=0:N-1;x=ones(1,N);lfft=512;
stem(n,x);title('time series of MA');% plot x(n)
xlabel('time index');ylabel('amplitude');
WT=[1:256]/256;WT=[-fliplr(WT),0,WT(2:end)]*3;%define WT for fft plot
X1=real(fft(x,lfft));X=[(X1(lfft/2:lfft)),X1(1:lfft/2-1)];%plot X [-pi,pi]
figure,plot(WT,X);%fft compute between[0,2pi] but we change to [-pi,pi]
k=-200:200;m=k/200*pi;%frequincy between [-pi,pi]
X=x*exp(-1j*(n')*m);% DTFT of x
fig=figure;plot(m,real(X));hold on % plot filter bank for m=0
for i=1:N-1 % plot filter bank for m= 1 to N-1
X=x*exp(-1j*(n')*(m+2*i*pi/N));
plot(m,real(X),'--');
end
plot(m,15*ones(1,length(m)));
axis([-pi pi -3 17])
title(['Modulated-Running-Sum Filter Bank for N= ',num2str(N)]);
xlabel('Normalized Frequency WT');ylabel('magnitude');
hold off
print(fig,'-depsc','uniform_filter15') % for more resolution print
%% plot frequency response of filter directly from it's formula
% see table 2.3 from DSP oppenhim second edition
w=-pi:.01:pi;f=sin(w*(N/2))./sin(w/2).*exp(-1j*w*(N-1)/2);
figure,plot(w,real(f));axis([-pi pi -1 5])