%Matlab file used for class demo conducted during Session 22
%illustrating efficacy of M-Channel PR Filter Bank
%based on a 2-channel QMF.  The 4 channel PR Bank is
%applied to a word utterance sampled at 12.5 KHz.
clf
clear all
set(0,'defaultaxesfontsize',20);
%Set up Efficient M-channel DFT filter bank:
M=8 %no. of channels
%
xdata=getspeech('enf1s1t0.wav');
Fs=12500;
%plot(xdata)
end1=500;
end2=13503+(M-4);
%end1=input('end1?');
%end2=input('end2?');
x=xdata(end1:end2);
dl=end2-end1+1;
%
input('Utterance played back at 12.5 KHz sampling rate'); 
soundsc(x,Fs)
%
mend=3;
end1=100;
n=-end1:end1; 
NL=201;
wm=2*pi/(NL+1);
cosw=0.5*ones(size(n))+0.5*cos(wm*n);
%cosw=ones(size(n));
n0=n+0.000001;
h0=sin(pi*n0)./(pi*n0); h0=zeros(1,1+2*end1); h0(1,end1+1)=1;
Hp(1,:)=h0;
for mk=1:M-1,
n1=n+(mk/M);
h1=sin(pi*n1)./(pi*n1);
for mi=1:mend,
h1=h1+sin(pi*(n1+mi*NL))./(pi*(n1+mi*NL)) + sin(pi*(n1-mi*NL))./(pi*(n1-mi*NL));
end
Hp(mk+1,:)=h1.*cosw;
end 
%
for m=1:M;
Tx(m,:)=x(1,m:M:end);
end
% 
Fp(1,:)=conv(Hp(1,:),Hp(1,:));
for m=2:M,
Fp(m,:)=conv(Hp(m,:),Hp(M-m+2,:));
end
%
Gp(1,:)=Hp(1,:);
for m=2:M,
Gp(m,:)=Hp(M-m+2,:);
end
%
%create y_0[n] and y_M-1 [n]
YpL=max(size(Tx))+max(size(Fp))-1;
Yr(1,:)=zeros(1,M+M*YpL);
for m=1:M;
Zp(m,:)=conv(Tx(m,:),Hp(m,:));
end
Wp=fft(Zp,M);
Vp=ifft(Wp,M);
for m=1:M;
Yp(m,:)=conv(Vp(m,:),Gp(m,:));
%Yp(m,:)=conv(Tx(m,:),Fp(m,:));
end
%
Yr(1,1:M:end)=[Yp(1,:) 0];
for m=2:M,
Yr(1,m:M:end)=[0 Yp(m,:)];
end
%
input('QMF Output played back at Fs=12.5 KHz'); 
soundsc(real(Yr),Fs)
%plot and compare DTFT's of x[n] and y[n]
Nfft=2^16;
omega=linspace(-pi,pi,Nfft);
yf1=abs(fftshift(fft(x,Nfft)));
%yf2=M*abs(fftshift(fft(y,Nfft)));
yf2=abs(fftshift(fft(Yr,Nfft)));
subplot(211)
plot(omega,yf1,'Linewidth',3)
axis([-pi pi 0 max(yf1)])
title('DTFT of original utterance');
%xlabel('omega (radians/s)');
subplot(212)
plot(omega,yf2,'Linewidth',3)
axis([-pi pi 0 max(yf2)])
title('DTFT of output of QMF');
xlabel('omega (radians/s)');
Nfft=2048;
%
h0f=2*abs(fftshift(fft(Fp(1,:),Nfft)));
h1f=2*abs(fftshift(fft(Fp(2,:),Nfft)));
h2f=2*abs(fftshift(fft(Fp(3,:),Nfft)));
h3f=2*abs(fftshift(fft(Fp(4,:),Nfft)));
%
input('Plot DTFT of h0[n] thru h3[n]')
clf
domega=2*pi/Nfft;
omega=-pi:domega:pi-domega;
subplot(2,2,1)
plot(omega,h0f,'b','Linewidth',4)
axis([-pi pi 0 max(h0f)])
legend('H0(w)'); xlabel('omega (radians/s)');
%hold on
subplot(2,2,2)
plot(omega,h1f,'r','Linewidth',4)
axis([-pi pi 0 max(h0f)])
legend('H1(w)'); xlabel('omega (radians/s)');
subplot(2,2,3)
plot(omega,h2f,'g','Linewidth',4)
axis([-pi pi 0 max(h0f)])
legend('H2(w)'); xlabel('omega (radians/s)');
subplot(2,2,4)
plot(omega,h3f,'c','Linewidth',4)
axis([-pi pi 0 max(h0f)])
legend('H3(w)'); xlabel('omega (radians/s)');
title('Frequency Response of h0[n] and h1[n]');
xlabel('omega (radians/s)');
%hold off
pause
clf
domega=2*pi/Nfft;
omega=-pi:domega:pi-domega;
hfa1=angle(exp(j*end1*omega).*fftshift(fft(Hp(2,:),Nfft)));
hfa2=angle(exp(j*end1*omega).*fftshift(fft(Hp(3,:),Nfft)));
hfa3=angle(exp(j*end1*omega).*fftshift(fft(Hp(4,:),Nfft)));
plot(omega,hfa1,'Linewidth',3)
axis([-pi pi -pi pi])
hold on 
plot(omega,hfa2,'-r','Linewidth',3)
plot(omega,hfa3,'-k','Linewidth',3)
legend('Phase of H1(w): slope=1/M','Phase of H2(w): slope=2/M','Phase of H3(w): slope=3/M')
hold off
title('Phase of 3 Polyphase Filters');


