%Matlab file used for class demo conducted during 
%Session 25 illustrating the efficacy of a nonuniform 
%M=4 Channel PR Filter Bank derived from a
%Multiresolution dyadic tree-structured filter bank
%derived from a 2-channel QMF.  The M=4 nonuniform
%channel PR Bank is %applied to a word utterance sampled 
%at 12.5 KHz (but downconverted to 10 KHz)>
clf
clear all
set(0,'defaultaxesfontsize',20);
%Set up M-channel DFT filter bank:
N=25; beta=0.1;
%N=17; beta=0.35;
n=-N:N-1;  
n=n+0.5;
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));
h0=h/sqrt(h*h');
h1=(-1).^(0:length(n)-1).*h0;
h1=h1/sqrt(h1*h1');
%
%h=[1 1]/sqrt(2);
%h0=h; h1=[1 -1]/sqrt(2);
%
h02=zeros(1,2*length(h)); h12=h02;
h02(1,1:2:length(h02))=h0;
h12(1,1:2:length(h12))=h1;
h04=zeros(1,2*length(h02)); h14=h04;
h04(1,1:2:length(h04))=h02;
h14(1,1:2:length(h14))=h12;
%
h2z=conv(h0,h12);
ht3=conv(h0,h02);  h3z=conv(ht3,h14);  
h4z=conv(ht3,h04);
%
h00=h0; h11=h1;
%
%Rename the filters so that h0[n] is the lowpass
%filter and h3[n] is the highpass filter
%
h3=h1;
h2=h2z(1,1:(length(h2z)-1));
h1=h3z(1,1:(length(h3z)-4));
h0=h4z(1,1:(length(h4z)-4));
%
g0=h00;  g1=-h11;
g02=zeros(1,2*length(g0)); g12=g02;
g02(1,1:2:length(g02))=g0;
g12(1,1:2:length(g12))=g1;
g04=zeros(1,2*length(g02)); g14=g04;
g04(1,1:2:length(g04))=g02;
g14(1,1:2:length(g14))=g12;
%
g2z=conv(g0,g12);
gt3=conv(g0,g02);  g3z=conv(gt3,g14);  
g4z=conv(gt3,g04);
%
%
g3=-h11;
g2=g2z(1,1:(length(g2z)-1));
g1=g3z(1,1:(length(g3z)-4));
g0=g4z(1,1:(length(g4z)-4));
%
data=getspeech('00f1s2t0.wav');
Fs=12500;
clf
plot(data)
end1=3000;
end2=10000;
%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 x_0 [n] through x_M-1 [n]
%
w0=conv(x,h0);
x0=w0(1,1:8:length(w0));
%
w1=conv(x,h1);
x1=w1(1,1:8:length(w1));
%
w2=conv(x,h2);
x2=w2(1,1:4:length(w2));
%
w3=conv(x,h3);
x3=w3(1,1:2:length(w3));
%
%create y_0[n] and y_M-1 [n]
%
z0=zeros(1,8*length(x0));
z0(1,1:8:length(z0))=x0;
y0=conv(z0,g0);
%
z1=zeros(1,8*length(x1));
z1(1,1:8:length(z1))=x1;
y1=conv(z1,g1);
%
z2=zeros(1,4*length(x2));
z2(1,1:4:length(z2))=x2;
y2=conv(z2,g2);
%
z3=zeros(1,2*length(x3));
z3(1,1:2:length(z3))=x3;
y3=conv(z3,g3);
%
L1=round((length(y0)-length(y3))/2);
L2=round((length(y0)-length(y2))/2);
L1=L1+1;  L2=L2+1;
y=zeros(1,length(y0));
y(1,L1:L1+length(y3)-1)=y(1,L1:L1+length(y3)-1)+y3;
y(1,L2:L2+length(y2)-1)=y(1,L2:L2+length(y2)-1)+y2;
y(1,1:length(y1))=y(1,1:length(y1))+y1;
y(1,1:length(y0))=y(1,1:length(y0))+y0;
%
input('PR Bank 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',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 Nonuniform PR');
xlabel('omega (radians/s)');
%plot filter responses h[n]
h0f=abs(fftshift(fft(h0,512)));
h1f=abs(fftshift(fft(h1,512)));
h2f=abs(fftshift(fft(h2,512)));
h3f=abs(fftshift(fft(h3,512)));
input('Plot DTFT of h_0 [n] thru h_3 [n]')
clf
Hmax=max([max(h0f) max(h1f) max(h2f) max(h3f)]);
domega=2*pi/512;
omega=-pi:domega:pi-domega;
plot(omega,h0f,'b','Linewidth',4)
axis([-pi pi 0 Hmax])
hold on
plot(omega,h1f,'r','Linewidth',4)
plot(omega,h2f,'g','Linewidth',4)
plot(omega,h3f,'c','Linewidth',4)
legend('H_0 (w)','H_1 (w)','H_2 (w)','H_3 (w)');
title('Frequency Response of h0[n] and h1[n]');
xlabel('omega (radians/s)');
hold off