%Matlab file
clf
clear all
set(0,'defaultaxesfontsize',20);
wc=7*pi/8;
end1=75;
n=-end1:end1; 
beta=0.15;
NL=202;
n0=n+0.000001;
h0=sin(pi*n0)./(pi*n0); h0=zeros(1,1+2*end1); h0(1,end1+1)=1;
%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));
%
s0=getspeech('erf1s1t0.wav'); x0=s0(500:13500);
s1=getspeech('enf1s1t0.wav'); x1=s1(500:13500);
s2=getspeech('0wf1s1t0.wav'); x2=s2(500:13500);
Fs=12500;
%
hcmp=2*(sin((pi/2)*n0)./(pi*n0)).*exp(j*(pi/2)*n);
s0=conv(x0,hcmp); s1=conv(x1,hcmp); s2=conv(x2,hcmp);
%
% c1=exp(j*2*pi/3); c2=exp(j*4*pi/3); c3=exp(j*8*pi/3);
% s0=x0+x1+x2;
% s1=x0+c1*x1+c2*x2;
% s2=x0+c2*x1+c3*x2;
W3=dftmtx(3); W3I=W3';
%X=[x0 ; x1 ; x2];
%S=W3I*X;
X=[s0 ; s1 ; s2];
S=3*ifft(X,3);
s0=real(S(1,:)); s1=real(S(2,:)); s2=real(S(3,:));
% s0=S(1,:); s1=S(2,:); s2=S(3,:);
%
y0=conv(s0,h0); y1=conv(s1,h1); y2=conv(s2,h2);
l3=3*length(y0);
%y=zeros(1,l3);
y(1:3:l3)=y0; y(2:3:l3)=y1; y(3:3:l3)=y2;
%
domega=2*pi/16384;
omega=-pi:domega:pi-domega;
yf=abs(fftshift(fft(y(50:end-50),16384)));
plot(omega,yf,'Linewidth',3)
axis([-pi pi 0 max(yf)])
title('DTFT of Sum');
xlabel('omega (radians/s)');
%
beta=0.15;
n0=n+0.000001;
h0=sin(pi*n0)./(pi*n0); h0=zeros(1,1+2*end1); h0(1,end1+1)=1;
%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));
%
g0=conv([y(1:3:l3)],h0);
g1=conv([y(2:3:l3)],h1);
g2=conv([y(3:3:l3)],h2);
%
%c1=exp(-j*2*pi/3); c2=exp(-j*4*pi/3); c3=exp(-j*8*pi/3); 
%
%z0=g0+g1+g2;
%z1=g0+c1*g1+c2*g2;
%z2=g0+c2*g1+c3*g2;
G=[g0 ; g1 ; g2];
%Z=W3*G;
Z=fft(G,3);
%
pause
soundsc(real(Z(1,:)),Fs)
%
Zfilt2=conv(Z(2,:),hcmp);
pause
soundsc(real(Zfilt2),Fs)
%soundsc(real(Z(2,:)),Fs)
Zfilt3=conv(Z(3,:),hcmp);
pause
soundsc(real(Zfilt3),Fs)
%soundsc(real(Z(3,:)),Fs)
%
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)));
%
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,hf1,'-b','Linewidth',3)
axis([-pi pi 0 max(hf0)])
title('Magnitude of Each Filter');
xlabel('omega (radians/s)');
%SlopeH1=mean(diff(hfa1(24:end-24)))/domega
%SlopeH2=mean(diff(hfa2(24:end-24)))/domega
%
domega=2*pi/16384;
omega=-pi:domega:pi-domega;
Zf=abs(fftshift(fft(Z(2,:),16384)));
plot(omega,Zf,'Linewidth',3)
axis([-pi pi 0 max(Zf)])
title('DTFT');
xlabel('omega (radians/s)');
