%Matlab file used for class demo conducted during Session 21
%illustrating efficacy of Two-Channel PR Filter Bank
%based on Mother Wavelet having a 
%raised cosine spectrum.  The QMF is
%applied to a word utterance sampled at 12.5 KHz.
clf
clear all
set(0,'defaultaxesfontsize',20);
%Set up 2-channel QMF filter bank:
N=24; beta=0.1;
n=-N:(N-1);  
n=n+0.5;
%h=(sin(pi*n/2)./(pi*n/2)).*(cos(pi*beta*n/2)./(1-beta^2*n.^2));
h=2*beta*cos((1+beta)*pi*n/2)./(pi*(1-4*beta^2*n.^2));
h=h+sin((1-beta)*pi*n/2)./(pi*(n-4*beta^2*n.^3));
%h=h/2;
h0=h;
h1=(-1).^(0:(length(n)-1)).*h;
g0=h;
g1=-h1;
data=getspeech('erf1s1t0.wav');
Fs=12500;
clf
plot(data)
end1=input('end1?');
end2=input('end2?');
xr=data(end1:end2);
dl=end2-end1+1;
%change sampling rate to 10 KHz
x=resample(xr,4,5);
Fs=(4/5)*Fs;
input('Utterance played back at 10 KHz sampling rate'); 
soundsc(x,Fs)
%create x0[n] and x1[n]
w0=conv(x,h0); x0=w0(1:2:length(w0));
w1=conv(x,h1); x1=w1(1:2:length(w1));
%create y0[n] and y1[n]
z0=zeros(1,2*length(x0));
z0(1:2:length(z0))=x0;
z1=zeros(1,2*length(x1));
z1(1:2:length(z1))=x1;
y0=conv(z0,g0);
y1=conv(z1,g1);
y=y0+y1;
input('QMF Output played back at Fs=10 KHz'); 
soundsc(y,Fs)
%plot and compare DTFT's of x[n] and y[n]
domega=2*pi/8192;
omega=-pi:domega:pi-domega;
yf1=abs(fftshift(fft(x,8192)));
yf2=abs(fftshift(fft(y,8192)));
subplot(211)
plot(omega,yf1,'Linewidth',4)
axis([-pi pi 0 max(yf1)])
title('DTFT of original utterance');
%xlabel('omega (radians/s)');
subplot(212)
plot(omega,yf2,'Linewidth',4)
axis([-pi pi 0 max(yf2)])
title('DTFT of output of QMF');
xlabel('omega (radians/s)');
%plot filter response h[n]
hlf=abs(fftshift(fft(h0,512)));
hpf=abs(fftshift(fft(h1,512)));
input('Plot DTFT of h0[n] and h1[n]')
clf
domega=2*pi/512;
omega=-pi:domega:pi-domega;
plot(omega,hlf,'b','Linewidth',4)
ps=hlf.^2+hpf.^2;
axis([-pi pi 0 max(ps)])
hold on
plot(omega,hpf,'r','Linewidth',4)
plot(omega,ps,'k','Linewidth',4)
legend('H0(w)','H1(w)','|H0(w)|^2+|H1(w)|^2');
title('Frequency Response of h0[n] and h1[n]');
xlabel('omega (radians/s)');
hold off
pause
%plot and compare DTFT's of x0[n] and x1[n]
domega=2*pi/4096;
omega=-pi:domega:pi-domega;
xf1=abs(fftshift(fft(x0,4096)));
xf2=abs(fftshift(fft(x1,4096)));
subplot(211)
plot(omega,xf1,'Linewidth',3)
axis([-pi pi 0 max(xf1)])
title('DTFT of x0[n]');
%xlabel('omega (radians/s)');
subplot(212)
plot(omega,xf2,'Linewidth',3)
axis([-pi pi 0 max(xf2)])
title('DTFT of x1[n]');
xlabel('omega (radians/s)');