%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); wc=7*pi/8; end1=100; n=-end1:end1; NL=202; beta=0.15; n0=n+0.000001; h0=sin(pi*n0)./(pi*n0); %h0=(sin(wc*n0)./(wc*n0)).*(cos(beta*wc*n0)./(1-(4/pi^2)*beta^2*(wc*n0).^2)); n1=n+(1/3); h1=sin(pi*n1)./(pi*n1); h1=h1+sin(pi*(n1+NL))./(pi*(n1+NL)) + sin(pi*(n1-NL))./(pi*(n1-NL)); h1=h1+sin(pi*(n1+2*NL))./(pi*(n1+2*NL)) + sin(pi*(n1-2*NL))./(pi*(n1-2*NL)); %h1=(sin(wc*n1)./(wc*n1)).*(cos(beta*wc*n1)./(1-(4/pi^2)*beta^2*(wc*n1).^2)); n2=n+(2/3); h2=sin(pi*n2)./(pi*n2); h2=h2+sin(pi*(n2+NL))./(pi*(n2+NL)) + sin(pi*(n2-NL))./(pi*(n2-NL)); h2=h2+sin(pi*(n2+2*NL))./(pi*(n2+2*NL)) + sin(pi*(n2-2*NL))./(pi*(n2-2*NL)); %h2=(sin(wc*n2)./(wc*n2)).*(cos(beta*wc*n2)./(1-(4/pi^2)*beta^2*(wc*n2).^2)); 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); l3=3*length(y1); y=zeros(1,l3); y(1:3:l3)=y1; y(2:3:l3)=y2; y(3:3:l3)=y3; input('Utterance interpolated by a factor played at thrice sampling rate, 3*Fs = 37.5 KHz'); soundsc(y,3*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))); 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) legend('Phase of H0(w): slope=0','Phase of H1(w): slope=1/3','Phase of H2(w): slope=2/3') hold off title('Phase of 3 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