%Matlab file
clf
clear all
set(0,'defaultaxesfontsize',20);
%wc=7*pi/8; beta=0.15;
mend=3;
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));
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*(n3+mi*NL))./(pi*(n3+mi*NL)) + sin(pi*(n3-mi*NL))./(pi*(n3-mi*NL));
end
%
h1=h1.*cosw; h2=h2.*cosw; h3=h3.*cosw;
%
s0=getspeech('erf1s1t0.wav'); x0=s0(500:12500);
s1=getspeech('enf1s1t0.wav'); x1=s1(500:12500);
s2=getspeech('0wf1s1t0.wav'); x2=s2(500:12500);
s3=getspeech('0ef1s1t0.wav'); x3=s3(500:12500);
Fs=12500;
%
clear s0 s1 s2 s3
%
W4=dftmtx(4); W4I=W4';
X=[x0 ; x1 ; x2 ; x3];
S=W4I*X;
%
%S=4*ifft(X,4);
s0=S(1,:);
s1=S(2,:);
s2=S(3,:);
s3=S(4,:);
%
y0=conv(s0,h0); y1=conv(s1,h1); y2=conv(s2,h2); y3=conv(s3,h3);
l3=4*length(y0);
%y=zeros(1,l3);
y(1:4:l3)=y0; y(2:4:l3)=y1; y(3:4:l3)=y2; y(4:4:l3)=y3;
%
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)');
%
pause
%
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*(n3+mi*NL))./(pi*(n3+mi*NL)) + sin(pi*(n3-mi*NL))./(pi*(n3-mi*NL));
end
%
h1=h1.*cosw; h2=h2.*cosw; h3=h3.*cosw;
%
g0=conv([y(1:4:l3)],h0);
g1=conv([y(2:4:l3)],h1);
g2=conv([y(3:4:l3)],h2);
g3=conv([y(4:4:l3)],h3);
% %
% g0=conv([y(1:4:l3) 0],h0);
% g1=conv([0 y(2:4:l3)],h3);
% g2=conv([0 y(3:4:l3)],h2);
% g3=conv([0 y(4:4:l3)],h1);
%
G=[g0 ; g1 ; g2 ; g3];
Z=W4*G;
%Z=fft(G,3);
%
pause
soundsc(real(Z(1,:)),Fs)
pause
soundsc(real(Z(2,:)),Fs)
pause
soundsc(real(Z(3,:)),Fs)
pause
soundsc(real(Z(4,:)),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)));
hf3=abs(fftshift(fft(h3,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)));
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,'-g','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)
hold on
plot(omega,hf1,'-r','Linewidth',3)
plot(omega,hf2,'-k','Linewidth',3)
plot(omega,hf3,'-g','Linewidth',3)
axis([-pi pi 0 max(hf0)])
title('Magnitude of Each Filter');
xlabel('omega (radians/s)');
hold off
%SlopeH1=mean(diff(hfa1(24:end-24)))/domega
%SlopeH2=mean(diff(hfa2(24:end-24)))/domega

