不明白的是双线性变换法得倒的滤波器没有冲激不变法的好。
%要求设计一chebyshev2低通数字滤波器，wp=30hz，ws=40hz，rp=0.5，rs=40，fs=100hz. wp=30;ws=40;rp=0.5;rs=40;fs=100; wp=30*2*pi;ws=40*2*pi;
[n,wn]=ellipord(wp,ws,rp,rs,'s'); [z,p,k]=ellipap(n,rp,rs); [num,den]=zp2tf(z,p,k);
[num1,den1]=impinvar(num,den); [num2,den2]=bilinear(num,den,100); [h,w]=freqz(num1,den1); [h1,w1]=freqz(num2,den2); subplot(1,2,1);
plot(w*fs/(2*pi),abs(h)); subplot(1,2,2);
plot(w1*fs/(2*pi),abs(h1)); figure(2);
subplot(1,2,1);
zplane(num1,den1); subplot(1,2,2);
zplane(num2,den2);
Butterworth带通
wp1=10;wp2=20;ws1=5;ws2=25;fs=100;rp=0.5;rs=50; wp1=10*2*pi;wp2=20*2*pi;ws1=5*2*pi;ws2=25*2*pi; wp=[wp1,wp2];ws=[ws1,ws2]; [n,wn]=buttord(wp,ws,rp,rs,'s');
bw=wn(2)-wn(1);wo=sqrt(wn(1)*wn(2)); [z,p,k]=buttap(n); [a,b,c,d]=zp2ss(z,p,k);
[a1,b1,c1,d1]=lp2bp(a,b,c,d,wo,bw); [num,den]=ss2tf(a1,b1,c1,d1); [h3,w3]=freqs(num,den); figure(4);
plot(w3/(2*pi),abs(h3)); title(‘模拟滤波器的幅频特性’); figure(5);
plot(w3/(2*pi),angle(h3));
title(‘模拟滤波器的相频特性’); figure(6);
zplane(num,den);
title(‘模拟滤波器的零极点图’); [num1,den1]=impinvar(num,den,fs); [h2,w2]=freqz(num1,den1);
[num2,den2]=bilinear(num,den,fs); [h1,w1]=freqz(num2,den2); figure(1);
subplot(1,2,1);
plot(w2*fs/(2*pi),abs(h2));
title(‘冲击不变法数字滤波器的幅频特性’); subplot(1,2,2);
plot(w1*fs/(2*pi),abs(h1));
title(‘双线性法数字滤波器的幅频特性’); figure(2);
subplot(1,2,1);
zplane(num1,den1);
title(‘冲击不变法数字滤波器的零极点图’); subplot(1,2,2);
zplane(num2,den2);
title(‘双线性法数字滤波器的零极点图’); figure(3);
subplot(1,2,1);
plot(w2*fs/(2*pi),angle(h2));
title(‘冲击不变法数字滤波器的相频特性’); subplot(1,2,2);
plot(w1*fs/(2*pi),angle(h1));
title(‘双线性法数字滤波器的相频特性’);
Chebshev1通带
wp1=10;wp2=20;ws1=5;ws2=25;fs=100;rp=0.5;rs=50; wp1=10*2*pi;wp2=20*2*pi;ws1=5*2*pi;ws2=25*2*pi; wp=[wp1,wp2];ws=[ws1,ws2]; [n,wn]=cheb1ord(wp,ws,rp,rs,'s');
bw=wn(2)-wn(1);wo=sqrt(wn(1)*wn(2)); [z,p,k]=cheb1ap(n,rp); [a,b,c,d]=zp2ss(z,p,k);
[a1,b1,c1,d1]=lp2bp(a,b,c,d,wo,bw); [num,den]=ss2tf(a1,b1,c1,d1);
[num1,den1]=impinvar(num,den,fs); [h2,w2]=freqz(num1,den1);
[num2,den2]=bilinear(num,den,fs); [h1,w1]=freqz(num2,den2); subplot(1,2,1);
plot(w2*fs/(2*pi),abs(h2)); subplot(1,2,2);
plot(w1*fs/(2*pi),abs(h1)); figure(1);
subplot(1,2,1);
zplane(num1,den1); subplot(1,2,2);
zplane(num2,den2);
Chebshev2通带
wp1=10;wp2=20;ws1=5;ws2=25;fs=100;rp=0.5;rs=50; wp1=10*2*pi;wp2=20*2*pi;ws1=5*2*pi;ws2=25*2*pi; wp=[wp1,wp2];ws=[ws1,ws2]; [n,wn]=cheb2ord(wp,ws,rp,rs,'s');
bw=wn(2)-wn(1);wo=sqrt(wn(1)*wn(2)); [z,p,k]=cheb2ap(n,rs); [a,b,c,d]=zp2ss(z,p,k);
[a1,b1,c1,d1]=lp2bp(a,b,c,d,wo,bw); [num,den]=ss2tf(a1,b1,c1,d1);
[num1,den1]=impinvar(num,den,fs); [h2,w2]=freqz(num1,den1);
[num2,den2]=bilinear(num,den,fs); [h1,w1]=freqz(num2,den2); subplot(1,2,1);
plot(w2*fs/(2*pi),abs(h2)); subplot(1,2,2);
plot(w1*fs/(2*pi),abs(h1)); figure(1);
subplot(1,2,1);
zplane(num1,den1); subplot(1,2,2);
zplane(num2,den2);
椭圆通带
wp1=10;wp2=20;ws1=5;ws2=25;fs=100;rp=0.5;rs=50; wp1=10*2*pi;wp2=20*2*pi;ws1=5*2*pi;ws2=25*2*pi; wp=[wp1,wp2];ws=[ws1,ws2]; [n,wn]=ellipord(wp,ws,rp,rs,'s');
bw=wn(2)-wn(1);wo=sqrt(wn(1)*wn(2)); [z,p,k]=ellipap(n,rp,rs); [a,b,c,d]=zp2ss(z,p,k);
[a1,b1,c1,d1]=lp2bp(a,b,c,d,wo,bw); [num,den]=ss2tf(a1,b1,c1,d1);
[num1,den1]=impinvar(num,den,fs); [h2,w2]=freqz(num1,den1);
[num2,den2]=bilinear(num,den,fs); [h1,w1]=freqz(num2,den2); subplot(1,2,1);
plot(w2*fs/(2*pi),abs(h2)); subplot(1,2,2);
plot(w1*fs/(2*pi),abs(h1)); figure(1);
subplot(1,2,1);
zplane(num1,den1); subplot(1,2,2);
zplane(num2,den2);