polyphase.m

10channel的polyphase滤波器

function filter_ten(flag)
% filter_ten(flag) flag=0 for flat sidelobes, flag=1 for falling sidelobes
hh1=remez(169,[0 40 60 500]/500,[1 1 0 0],[1 100]);
frq=[0 40 60 99 100 149 150 199 200 249 250 299 300 349 350 399 400 449 450 500]/500;
gn= [1 1 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0];
pn= [ 1 100 140 180 220 260 300 340 380 420];
hh2=remez(169,frq,gn,pn);
hh=hh1;
flag=0;
if flag==1
hh=hh2;
end
figure(1)
subplot(2,1,1)
plot(hh)
grid
title('Impulse Response: Prototype Filter')
xlabel('Normalized time nT/T')
ylabel('Amplitude')
subplot(2,1,2)
plot((-0.5:1/1024:.5-1/1024)*1000,fftshift(20*log10(0.000001+abs(fft(hh,1024)))))
grid
axis([-500 500 -90 10])
title('Frequency Response: Prototype Filter')
xlabel('Frequency (kHz)')
ylabel('Log-Magnitude (dB)')
pause
x1=1*cos(2*pi*(0:2299)*10/1000);
x2=2+cos(2*pi*(0:2299)*15/1000);
x2=x2.*cos(2*pi*(0:2299)*100/1000);
x3=(3*cos(2*pi*(0:2299)*22/1000)+5*sin(2*pi*(0:2299)*6/1000));
x3=x3.*sin(2*pi*(0:2299)*300/1000);
xx=x1+x2+x3;
xx=[xx zeros(1,200)];
figure(2)
subplot(2,1,1)
plot(xx(1:200));
grid
title('Real Input Time Series')
xlabel('Normalized time nT/T')
ylabel('Amplitude')
subplot(2,1,2)
ww=kaiser(1024,8)';
ww=ww/sum(ww);
plot((-0.5:1/1024:.5-1/1024)*1000,fftshift(20*log10(abs(fft(xx(1:1024).*ww,1024)))))
grid
axis([-500 500 -90 10])
title('Spectrum of Real Input Series')
xlabel('Frequency (kHz)')
ylabel('Log-Magnitude (dB)')
pause
hold
plot((-0.5:1/1024:.5-1/1024)*1000,fftshift(20*log10(0.00001+abs(fft(hh,1024)))),'r')
pause
gg1=hh.*exp(j*2*pi*(-84.5:84.5)*100/1000);
plot((-0.5:1/1024:.5-1/1024)*1000,fftshift(20*log10(abs(fft(gg1,1024)))),'r--')
gg2=hh.*exp(j*2*pi*(-84.5:84.5)*200/1000);
plot((-0.5:1/1024:.5-1/1024)*1000,fftshift(20*log10(abs(fft(gg2,1024)))),'r--')
gg3=hh.*exp(j*2*pi*(-84.5:84.5)*300/1000);
plot((-0.5:1/1024:.5-1/1024)*1000,fftshift(20*log10(abs(fft(gg3,1024)))),'r--')
gg4=hh.*exp(j*2*pi*(-84.5:84.5)*400/1000);
plot((-0.5:1/1024:.5-1/1024)*1000,fftshift(20*log10(abs(fft(gg4,1024)))),'r--')
hold
pause
hh2=reshape(hh,10,17);
reg=zeros(10,17);
n2=1;
for nn=1:10:2500
reg(:,2:17)=reg(:,1:16);
reg(:,1)=flipud(xx(nn:nn+9)');
for mm=1:10
vv(mm)=reg(mm,:)*hh2(mm,:)';
end
yy(:,n2)=fft(vv)';
n2=n2+1;
end
figure(3)
subplot(5,2,1)
plot(real(yy(1,:)))
grid
title('Time Series, Channels 1-5')
rr=axis;
rr(2)=250;
axis(rr);
subplot(5,2,2)
ww=kaiser(200,8)';
ww=ww/sum(ww);
plot((-0.5:1/256:.5-1/256)*100,fftshift(20*log10(abs(fft(yy(1,20:219).*ww,256)))))
axis([-50 50 -80 10])
grid
title('Spectra, Channels 1-5')
subplot(5,2,3)
plot(real(yy(2,:)))
grid
rr=axis;
rr(2)=250;
axis(rr);
xlabel('Frequency (kHz)')

ylabel('Log-Magnitude (dB)')
pause
hold
plot((-0.5:1/1024:.5-1/1024)*1000,fftshift(20*log10(0.00001+abs(fft(hh,1024)))),'r')
pause
gg1=hh.*exp(j*2*pi*(-84.5:84.5)*100/1000);
plot((-0.5:1/1024:.5-1/1024)*1000,fftshift(20*log10(abs(fft(gg1,1024)))),'r--')
gg2=hh.*exp(j*2*pi*(-84.5:84.5)*200/1000);
plot((-0.5:1/1024:.5-1/1024)*1000,fftshift(20*log10(abs(fft(gg2,1024)))),'r--')
gg3=hh.*exp(j*2*pi*(-84.5:84.5)*300/1000);
plot((-0.5:1/1024:.5-1/1024)*1000,fftshift(20*log10(abs(fft(gg3,1024)))),'r--')
gg4=hh.*exp(j*2*pi*(-84.5:84.5)*400/1000);
plot((-0.5:1/1024:.5-1/1024)*1000,fftshift(20*log10(abs(fft(gg4,1024)))),'r--')
hold
pause
hh2=reshape(hh,10,17);
reg=zeros(10,17);
n2=1;
for nn=1:10:2500
reg(:,2:17)=reg(:,1:16);
reg(:,1)=flipud(xx(nn:nn+9)');
for mm=1:10
vv(mm)=reg(mm,:)*hh2(mm,:)';
end
yy(:,n2)=fft(vv)';
n2=n2+1;
end
figure(3)
subplot(5,2,1)
plot(real(yy(1,:)))
grid
title('Time Series, Channels 1-5')
rr=axis;
rr(2)=250;
axis(rr);
subplot(5,2,2)
ww=kaiser(200,8)';
ww=ww/sum(ww);
plot((-0.5:1/256:.5-1/256)*100,fftshift(20*log10(abs(fft(yy(1,20:219).*ww,256)))))
axis([-50 50 -80 10])
grid
title('Spectra, Channels 1-5')
subplot(5,2,3)
plot(real(yy(2,:)))
grid
rr=axis;
rr(2)=250;
axis(rr);
subplot(5,2,4)
plot((-0.5:1/256:.5-1/256)*100,fftshift(20*log10(abs(fft(yy(2,20:219).*ww,256)))))
axis([-50 50 -80 10])
grid
subplot(5,2,5)
plot(real(yy(3,:)))
grid
rr=axis;
rr(2)=250;
axis(rr);
subplot(5,2,6)
plot((-0.5:1/256:.5-1/256)*100,fftshift(20*log10(abs(fft(yy(3,20:219).*ww,256)))))
axis([-50 50 -80 10])
grid
subplot(5,2,7)
plot(real(yy(4,:)))
grid
rr=axis;
rr(2)=250;
axis(rr);
subplot(5,2,8)
plot((-0.5:1/256:.5-1/256)*100,fftshift(20*log10(abs(fft(yy(4,20:219).*ww,256)))))
axis([-50 50 -80 10])
grid
subplot(5,2,9)
plot(real(yy(5,:)))
grid
rr=axis;
rr(2)=250;
axis(rr);
subplot(5,2,10)
plot((-0.5:1/256:.5-1/256)*100,fftshift(20*log10(abs(fft(yy(5,20:219).*ww,256)))))
axis([-50 50 -80 10])
grid
pause
figure(4)
subplot(2,2,1)
plot(0,0)
hold
for mm=1:10
plot(0:1/64:1-1/64,(10*abs(fft(hh2(mm,:),64))))
end
hold
grid
title('Spectral Magnitude Response of Ten Polyphase Filters')
xlabel('Normalized Frequency')
ylabel('Amplitude')
axis([0 0.5 0.0 1.2])
subplot(2,2,3)
plot(0,0)

hold
for mm=1:10
plot(0:1:16,0.5*mm+10*hh2(mm,:))
end
plot([7.55 8.45],[5.628 1.128],'r')
hold
grid
title('Impulse Response of Ten Polyphase Filters')
xlabel('Normalized Time')
axis([-1 17 -0.2 6])
subplot(2,2,2)
plot(0,0)
hold
for mm=1:10
plot(0:1/64:1-1/64,unwrap((angle((fft(hh2(mm,:),64)))))/pi)
end
hold
grid
title('Spectral Phase Response of Ten Polyphase Filters')
axis([0 0.5 -9 0])
xlabel('Normalized Frequency (f/fs)')
ylabel('Phase Shift (\theta/2\pi)')
subplot(2,2,4)
plot(0,0)
hold
for mm=1:10
vv=unwrap(angle(fftshift(fft(hh2(mm,:),64))))/(2*pi);
vv2=64*filter([1 -1],1, vv);
plot(0:1/64:1-1/64,fftshift(vv2))
end
hold
grid
title('Spectral Group Delay Response of Ten Polyphase Filters')
xlabel('Normalized Frequency (f/fs)')
ylabel('Group Delay (d\theta/d\omega) (Samples)')
%axis([0 0.5 -9 0])
axis([0 0.5 -9 -7])
% pause
% figure(5)
% plot(0,0)
% hold
% for mm=1:10
% plot(0:1/64:1-1/64,unwrap((angle((fft(hh2(mm,:),64)))))/pi)
% end
% hold
% grid
% title('Spectral Phase Response of Ten Polyphase Filters')
% axis([0 0.5 -9 0])
% xlabel('Normalized Frequency (f/fs)')
% ylabel('Phase Shift (\theta/2\pi)')
% figure(6)
% plot(0,0)
% hold
% for mm=1:10
% vv=unwrap(angle(fftshift(fft(hh2(mm,:),64))))/(2*pi);
% vv2=64*filter([1 -1],1, vv);
% plot(0:1/64:1-1/64,fftshift(vv2))
% end
% hold
% grid
% title('Spectral Group Delay Response of Ten Polyphase Filters')
% axis([0 0.5 -9 -7])
% xlabel('Normalized Frequency (f/fs)')
% ylabel('Group Delay (d\theta/d\omega) (Samples)')
%
% gg=get(gca);
% set(gca,'gridlinestyle','-')