%Matlab file used for class demo conducted during clas
%illustrating upsampling by a factor of 3 using the
%efficient polyphase implementation of upsampling.
%A word utterance is sampled at 12.5 KHz, upsampled by
%a factor of 3 and played back at 37.5 KHz
clf
clear all
set(0,'defaultaxesfontsize',20);
end1=100;
n=-end1:end1;
NL=202;
wm=2*pi/(NL+1);
cosw=0.5*ones(size(n))+0.5*cos(wm*n);
%cosw=ones(size(n));
mend=3;
n0=n+0.000001;
h0=sin(pi*n0)./(pi*n0); h0=zeros(1,1+2*end1); h0(1,end1+1)=1;
n1=n+(1/4);
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
n2=n+(2/4);
h2=sin(pi*n2)./(pi*n2);
for mi=1:mend,
h2=h2+sin(pi*(n2+mi*NL))./(pi*(n2+mi*NL)) + sin(pi*(n2-mi*NL))./(pi*(n2-mi*NL));
end
n3=n+(3/4);
h3=sin(pi*n3)./(pi*n3);
for mi=1:mend,
h3=h3+sin(pi*(n2+mi*NL))./(pi*(n2+mi*NL)) + sin(pi*(n2-mi*NL))./(pi*(n2-mi*NL));
end
%
h1=h1.*cosw; h2=h2.*cosw; h3=h3.*cosw;
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
x=getspeech('0wf1s1t0.wav');
x=x(500:13500);
Fs=12500;
input('Utterance played back at orignal sampling rate, Fs = 12.5 KHz'); 
soundsc(x,Fs)
y1=conv(x,h0); y2=conv(x,h1); y3=conv(x,h2); y4=conv(x,h3);
l4=4*length(y1);
y=zeros(1,l4);
y(1:4:l4)=y1; y(2:4:l4)=y2; y(3:4:l4)=y3; y(4:4:l4)=y4;
input('Utterance interpolated by a factor played at thrice sampling rate, 4*Fs = 50 KHz');
soundsc(y,4*Fs)
%
N1=16384;
N2=65536;
% domega=2*pi/16384;
% omega=-pi:domega:pi-domega;
omega1=linspace(-pi,pi,N1);
omega2=linspace(-pi,pi,N2);
yf1=abs(fftshift(fft(x,N1)));
yf2=abs(fftshift(fft(y(75:end-75),N2)));
subplot(211)
plot(omega1,yf1,'Linewidth',3)
axis([-pi pi 0 max(yf1)])
title('DTFT of original utterance');
%xlabel('omega (radians/s)');
subplot(212)
plot(omega2,yf2,'Linewidth',3)
axis([-pi pi 0 max(yf2)])
title('DTFT of utterance interpolated by a factor of 3');
xlabel('omega (radians/s)');
pause
domega=2*pi/1024;
omega=-pi:domega:pi-domega;
hf0=abs(fftshift(fft(h0,1024)));
hf1=abs(fftshift(fft(h1,1024)));
hf2=abs(fftshift(fft(h2,1024)));
subplot(211)
plot(omega,hf1,'Linewidth',3)
axis([-pi pi 0 max(hf1)])
title('DTFT of Filter 1');
%xlabel('omega (radians/s)');
subplot(212)
plot(omega,hf2,'Linewidth',3)
axis([-pi pi 0 max(hf2)])
title('DTFT of Filter 2');
xlabel('omega (radians/s)');
pause
domega=2*pi/1024;
omega=-pi:domega:pi-domega;
hfa0=angle(exp(j*end1*omega).*fftshift(fft(h0,1024)));
hfa1=angle(exp(j*end1*omega).*fftshift(fft(h1,1024)));
hfa2=angle(exp(j*end1*omega).*fftshift(fft(h2,1024)));
hfa3=angle(exp(j*end1*omega).*fftshift(fft(h3,1024)));
subplot(211)
plot(omega,hfa0,'Linewidth',3)
axis([-pi pi -pi pi])
hold on 
plot(omega,hfa1,'-r','Linewidth',3)
plot(omega,hfa2,'-k','Linewidth',3)
plot(omega,hfa3,'-m','Linewidth',3)
legend('Phase of H0(w): slope=0','Phase of H1(w): slope=1/4','Phase of H2(w): slope=2/4','Phase of H3(w): slope=3/4')
hold off
title('Phase of 4 Polyphase Filters');
%xlabel('omega (radians/s)');
subplot(212)
plot(omega,hf0,'-b','Linewidth',3)
axis([-pi pi 0 max(hf0)])
title('Magnitude of Filter 0');
xlabel('omega (radians/s)');
SlopeH1=mean(diff(hfa1(24:end-24)))/domega
SlopeH2=mean(diff(hfa2(24:end-24)))/domega
SlopeH3=mean(diff(hfa3(24:end-24)))/domega

