1、数字信号处理上机实验 滤波器数字信号处理上机练习7.12（1）双线性变换法程序：%双线性变换法设计切比雪夫1型带通滤波器%w计算式：w=(f/ADfs)*2*pi 单位：radclc;clear;close all;ADfs=2000;%ADC的采样频率，单位：Hzws1=(100/ADfs)*2*pi;%下阻带截止频率，单位：radwp1=(200/ADfs)*2*pi;%通带下截止频率，单位：radwp2=(400/ADfs)*2*pi;%通带上截止频率，单位：radws2=(600/ADfs)*2*pi;%上阻带截止频率，单位：raddisplay(数字滤波器的边界频率 单位：rad);
2、ws1, wp1,wp2,ws2Rp=2;%通带最大衰减，单位：dBAs=20;%阻带最小衰减，单位：dBT=1;%模拟滤波器数字化时的采样间隔，单位：sWs1=(2/T)*tan(ws1/2);Wp1=(2/T)*tan(wp1/2);Wp2=(2/T)*tan(wp2/2);Ws2=(2/T)*tan(ws2/2);%转换为模拟滤波器指标Ws=Ws1,Ws2;Wp=Wp1,Wp2;display(转换后滤波器的边界频率 单位：rad/s);Ws1, Wp1,Wp2,Ws2 n,wn = cheb1ord(Wp,Ws,Rp,As,s);%计算切比雪夫1型滤波器阶数和截止频率display(滤
3、波器阶数);ndisplay(2dB 截止频率);wnb,a = cheby1(n,Rp,wn,s);% 计算模拟带通滤波器的传输函数Ha(s)display(模拟原型系统函数Ha(s)分子多项式系数：);bdisplay(模拟原型系统函数Ha(s)分母多项式系数);aHa,Wa = freqs(b,a,512);%计算并画出模拟原型滤波器的幅度频谱图subplot(211);plot(Wa/pi,20*log10(abs(Ha),LineWidth,2);axis(0,max(Wa/pi),min(20*log10(abs(Ha),5);xlabel(角频率rad/s();ylabel(对数
4、幅度dB);title(模拟原型带通滤波器的幅度谱（对数幅度）);grid;bz,az=bilinear(b,a,1/T);%用双线性变换法转换为数字滤波器display(数字滤波器的系统函数H(z)分子多项式系数：);bzdisplay(数字滤波器的系统函数H(z)分母多项式系数);azHz,Wz = freqz(bz,az,512,1/T); %计算并画出数字滤波器的幅度频谱图，Wz是以Hz为单位subplot(212);plot(Wz*2*T,20*log10(abs(Hz),LineWidth,2);axis(0,1,-100,5);xlabel(数字域频率();ylabel(对数幅度
5、dB);title(双线性变换法设计的数字滤波器的幅度谱（对数幅度）);grid;结果：数字滤波器的边界频率 单位：radws1 = 0.3142wp1 = 0.6283wp2 =1.2566ws2 = 1.8850转换后滤波器的边界频率 单位：rad/sWs1 = 0.3168Wp1 = 0.6498Wp2 =1.4531Ws2 = 2.7528滤波器阶数n = 22dB 截止频率wn = 0.6498 1.4531模拟原型系统函数Ha(s)分子多项式系数：b = 0 0 0.4218 0 0模拟原型系统函数Ha(s)分母多项式系数a = 1.0000 0.6457 2.4196 0.609
6、7 0.8916数字滤波器的系统函数H(z)分子多项式系数：bz = 0.0512 0.0000 -0.1024 0.0000 0.0512数字滤波器的系统函数H(z)分母多项式系数az = 1.0000 -2.0733 2.4881 -1.5944 0.6125（2）冲激响应不变法程序：%冲激响应不变法设计切比雪夫1型带通滤波器%w计算式：w=(f/ADfs)*2*pi 单位：radclc;clear;close all;ADfs=2000;%ADC的采样频率，单位：Hzws1=(100/ADfs)*2*pi;%下阻带截止频率，单位：radwp1=(200/ADfs)*2*pi;%通带下截止
7、频率，单位：radwp2=(400/ADfs)*2*pi;%通带上截止频率，单位：radws2=(600/ADfs)*2*pi;%上阻带截止频率，单位：raddisplay(数字滤波器的边界频率 单位：rad);ws1, wp1,wp2,ws2 Rp=2;%通带最大衰减，单位：dBRs=20;%阻带最小衰减，单位：dBT=1;%模拟滤波器数字化时的采样间隔，单位：sWs1=ws1/T;Wp1=wp1/T;Wp2=wp2/T;Ws2=ws2/T;%转换为模拟滤波器指标display(转换后滤波器的边界频率 单位：rad/s);Ws1, Wp1,Wp2,Ws2 Ws=Ws1,Ws2;Wp=Wp1,
8、Wp2;n,wn = cheb1ord(Wp,Ws,Rp,Rs,s);%计算切比雪夫1型滤波器阶数和截止频率display(滤波器阶数);ndisplay(2dB frequency);wnb,a = cheby1(n,Rp,wn,s);% 计算模拟带通滤波器的传输函数Ha(s)display(模拟原型系统函数Ha(s)分子多项式系数：);bdisplay(模拟原型系统函数Ha(s)分母多项式系数);aHa,Wa = freqs(b,a,512);%计算并画出模拟原型滤波器的幅度频谱图subplot(211);plot(Wa/pi,20*log10(abs(Ha),LineWidth,2);a
9、xis(0,1.5,-80,5);xlabel(角频率rad/s();ylabel(对数幅度dB);title(模拟原型带通滤波器的幅度谱（对数幅度）);grid;bz,az=impinvar(b,a,1/T);%用冲激响应不变法转换为数字滤波器display(数字滤波器的系统函数H(z)分子多项式系数：);bzdisplay(数字滤波器的系统函数H(z)分母多项式系数);azHz,Wz = freqz(bz,az,512,1/T); %计算并画出数字滤波器的幅度频谱图，Wz是以Hz为单位subplot(212);plot(Wz*2*T,20*log10(abs(Hz),LineWidth,2
10、);xlabel(数字域频率();ylabel(对数幅度dB);title(冲激响应不变法设计的数字滤波器的幅度谱（对数幅度）);axis(0,1,-80,5)grid;结果：数字滤波器的边界频率 单位：radws1 = 0.3142wp1 = 0.6283wp2 = 1.2566ws2 = 1.8850转换后滤波器的边界频率 单位：rad/sWs1 = 0.3142Wp1 = 0.6283Wp2 = 1.2566Ws2 = 1.8850滤波器阶数n = 32dB frequencywn = 0.6283 1.2566模拟原型系统函数Ha(s)分子多项式系数：b = 0 0 0 0.0811
11、0 0 0模拟原型系统函数Ha(s)分母多项式系数a = 1.0000 0.4636 2.7722 0.8132 2.1889 0.2890 0.4922数字滤波器的系统函数H(z)分子多项式系数：bz = 0.0000 0.0272 -0.0581 0.0109 0.0437 -0.0237 0数字滤波器的系统函数H(z)分母多项式系数az = 1.0000 -3.3030 6.0060 -6.7463 5.1356 -2.4093 0.62907.14程序：%双线性变换法设计切比雪夫1型带通滤波器%w计算式：w=(f/ADfs)*2*pi 单位：radclc;clear;close all
12、;ADfs=10000;%ADC的采样频率，单位：Hzwp1=(500/ADfs)*2*pi;%下通带截止频率，单位：radws1=(1000/ADfs)*2*pi;%阻带下截止频率，单位：radws2=(2000/ADfs)*2*pi;%阻带上截止频率，单位：radwp2=(3000/ADfs)*2*pi;%上通带截止频率，单位：raddisplay(数字滤波器的边界频率 单位：rad);wp1, ws1,ws2,wp2 Rp=3;%通带最大衰减，单位：dBAs=30;%阻带最小衰减，单位：dBT=1;%模拟滤波器数字化时的采样间隔，单位：sWs1=(2/T)*tan(ws1/2);Wp1=
13、(2/T)*tan(wp1/2);Wp2=(2/T)*tan(wp2/2);Ws2=(2/T)*tan(ws2/2);%转换为模拟滤波器指标display(转换后滤波器的边界频率 单位：rad/s);Wp1, Ws1,Ws2,Wp2 Ws=Ws1,Ws2;Wp=Wp1,Wp2;n,wn = cheb1ord(Wp,Ws,Rp,As,s);%计算切比雪夫1型滤波器阶数和截止频率display(滤波器阶数);ndisplay(3dB 截止频率);wnb,a = cheby1(n,Rp,wn,stop,s);% 计算模拟带阻滤波器的传输函数Ha(s)display(模拟原型系统函数Ha(s)分子多项
14、式系数：);bdisplay(模拟原型系统函数Ha(s)分母多项式系数);aHa,Wa = freqs(b,a,512);%计算并画出模拟原型滤波器的幅度频谱图subplot(211);plot(Wa/pi,20*log10(abs(Ha),LineWidth,2);axis(0,3,-120,5);xlabel(角频率rad/s();ylabel(对数幅度dB);title(模拟原型带阻滤波器的幅度谱（对数幅度）);grid;bz,az=bilinear(b,a,1/T);%用双线性变换法转换为数字滤波器display(数字滤波器的系统函H(z)分子多项式系数：);bzdisplay(数字滤
15、波器的系统函H(z)分母多项式系数);azHz,Wz = freqz(bz,az,512,1/T); %计算并画出数字滤波器的幅度频谱图，Wz是以Hz为单位subplot(212);plot(Wz*2*T,20*log10(abs(Hz),LineWidth,2);xlabel(数字域频率();ylabel(对数幅度dB);title(双线性变换法设计的数字滤波器的幅度谱（对数幅度）);axis(0,1,-120,5)grid;结果：数字滤波器的边界频率 单位：radwp1 = 0.3142ws1 = 0.6283ws2 = 1.2566wp2 = 1.8850转换后滤波器的边界频率 单位：r
16、ad/sWp1 = 0.3168Ws1 = 0.6498Ws2 = 1.4531Wp2 = 2.7528滤波器阶数n = 33dB 截止频率wn = 0.3430 2.7527模拟原型系统函数Ha(s)分子多项式系数：b = 1.0000 0 2.8326 0 2.6745 0 0.8418模拟原型系统函数Ha(s)分母多项式系数a = 1.0000 8.9270 16.6716 72.6944 15.7413 7.9584 0.8418数字滤波器的系统函H(z)分子多项式系数：bz = 0.0946 -0.3508 0.7174 -0.8802 0.7174 -0.3508 0.0946数字
17、滤波器的系统函H(z)分母多项式系数az =1.0000 -1.4601 0.3179 -0.3507 0.6885 0.2289 -0.38247.15程序：%双线性变换法设计切比雪夫2型带阻滤波器%w计算式：w=(f/ADfs)*2*pi 单位：radclc;clear;close all;ADfs=100000;%ADC的采样频率，单位：Hzwp1=(5000/ADfs)*2*pi;%下通带截止频率，单位：radws1=(10000/ADfs)*2*pi;%阻带下截止频率，单位：radws2=(20000/ADfs)*2*pi;%阻带上截止频率，单位：radwp2=(30000/ADfs
18、)*2*pi;%上通带截止频率，单位：raddisplay(数字滤波器的边界频率 单位：rad);wp1, ws1,ws2,wp2 Rp=3;%通带最大衰减，单位：dBAs=30;%阻带最小衰减，单位：dBT=1;%模拟滤波器数字化时的采样间隔，单位：sWs1=(2/T)*tan(ws1/2);Wp1=(2/T)*tan(wp1/2);Wp2=(2/T)*tan(wp2/2);Ws2=(2/T)*tan(ws2/2);%转换为模拟滤波器指标display(转换后模拟原型滤波器的边界频率 单位：rad/s);Wp1, Ws1,Ws2,Wp2 Ws=Ws1,Ws2;Wp=Wp1,Wp2;n,wn
19、= cheb2ord(Wp,Ws,Rp,As,s);%计算切比雪夫2型滤波器阶数和截止频率display(滤波器阶数);ndisplay(3dB 截止频率);wnb,a = cheby2(n,As,wn,stop,s);% 计算模拟带阻滤波器的传输函数Ha(s)display(模拟原型滤波器系统函数Ha(s)分子多项式系数：);bdisplay(模拟原型滤波器系统函数Ha(s)分母多项式系数);aHa,Wa = freqs(b,a,512);%计算并画出模拟原型滤波器的幅度频谱图subplot(211);plot(Wa/pi,20*log10(abs(Ha),LineWidth,2);axis
20、(0,1.5,min(20*log10(abs(Ha),5);xlabel(角频率rad/s();ylabel(对数幅度dB);title(模拟原型带阻滤波器的幅度谱（对数幅度）);grid;bz,az=bilinear(b,a,1/T);%用双线性变换法转换为数字滤波器display(数字滤波器系统函数H(z)分子多项式系数：);bzdisplay(数字滤波器系统函数H(z)分母多项式系数);azHz,Wz = freqz(bz,az,512,1/T); %计算并画出数字滤波器的幅度频谱图，Wz是以Hz为单位subplot(212);plot(Wz*2*T,20*log10(abs(Hz),
21、LineWidth,2);xlabel(数字域频率();ylabel(对数幅度dB);title(双线性变换法设计的数字滤波器的幅度谱（对数幅度）);grid;结果：数字滤波器的边界频率 单位：radwp1 = 0.3142ws1 = 0.6283ws2 = 1.2566wp2 = 1.8850转换后模拟原型滤波器的边界频率 单位：rad/sWp1 = 0.3168Ws1 = 0.6498Ws2 = 1.4531Wp2 = 2.7528滤波器阶数n = 33dB 截止频率wn = 0.5572 1.6946模拟原型滤波器系统函数Ha(s)分子多项式系数：b = 1.0000 -0.0000 3
22、.8028 -0.0000 3.5906 -0.0000 0.8418模拟原型滤波器系统函数Ha(s)分母多项式系数a = 1.0000 4.2458 12.8161 19.6444 12.1009 3.7851 0.8418数字滤波器系统函数H(z)分子多项式系数：bz = 0.2263 -0.7625 1.4500 -1.7406 1.4500 -0.7625 0.2263数字滤波器系统函数H(z)分母多项式系数az = 1.0000 -1.9477 1.5590 -1.0285 0.7650 -0.2894 0.02867.17（1）椭圆函数滤波器程序：%设计椭圆函数带通数字滤波器%w计
23、算式： w/pi=2f/ADfsclc;clear;close all;ADfs=25000;%ADC的采样频率，单位：Hzws1=(3500/ADfs)*2;%下阻带截止频率，单位：*pi radwp1=(5000/ADfs)*2;%通带下截止频率，单位：*pi radwp2=(7000/ADfs)*2;%通带上截止频率，单位：*pi radws2=(8500/ADfs)*2;%上阻带截止频率，单位：*pi radRp=0.5;%通带最大衰减，单位：dBAs=45;%阻带最小衰减，单位：dBws=ws1,ws2;wp=wp1,wp2;n,wc = ellipord(wp,ws,Rp,As);
24、%计算椭圆函数带通数字滤波器阶数和截止频率display(数字滤波器边界频率 单位：*pi rad);ws1 ,wp1 , wp2 ,ws2display(滤波器阶数);ndisplay(0.5dB 截止频率 单位：*pi rad);wcb,a = ellip(n,Rp,As,wc);% 计算椭圆函数带通数字滤波器的传输函数H（z）display(H(z)分子多项式系数：);bdisplay(H(z)分母多项式系数);ahn=impz(b,a); %计算单位冲激响应H,w=freqz(b,a);%画出系统的幅频特性曲线subplot(3,1,1);plot(12.5*w/pi,20*log10
25、(abs(H),LineWidth,2);%xlabel(f(khz);ylabel(dB);title(椭圆函数带通数字滤波器的幅度响应（db）);axis(0,12.5 ,-70,3);grid; %画出系统的相频特性曲线subplot(3,1,2);plot(w/pi,angle(H),LineWidth,2);%xlabel(w/);ylabel(rad);title(椭圆函数带通数字滤波器的相位响应);grid; %画出单位冲激响应subplot(3,1,3);stem(hn,.,Linewidth,2);title(系统单位冲激响应);xlabel(n);axis(0,120,-0
26、.2,0.2);结果：数字滤波器边界频率 单位：*pi radws1 = 0.2800wp1 = 0.4000wp2 = 0.5600ws2 = 0.6800滤波器阶数n = 40.5dB 截止频率 单位：*pi radwc = 0.4000 0.5600H(z)分子多项式系数：b = 0.0113 -0.0036 0.0110 -0.0045 0.0196 -0.0045 0.0110 -0.0036 0.0113H(z)分母多项式系数a = 1.0000 -0.4649 3.2535 -1.1443 4.1557 -0.9889 2.4381 -0.2974 0.5525（2）切比雪夫型滤
27、波器程序：%设计切比雪夫1型带通数字滤波器%w计算式： w/pi=2f/ADfsclc;clear;close all;ADfs=25000;%ADC的采样频率，单位：Hzws1=(3500/ADfs)*2;%下阻带截止频率，单位：*pi radwp1=(5000/ADfs)*2;%通带下截止频率，单位：*pi radwp2=(7000/ADfs)*2;%通带上截止频率，单位：*pi radws2=(8500/ADfs)*2;%上阻带截止频率，单位：*pi radRp=0.5;%通带最大衰减，单位：dBAs=45;%阻带最小衰减，单位：dBws=ws1,ws2;wp=wp1,wp2;n,wc
28、= cheb1ord(wp,ws,Rp,As);%计算切比雪夫1型带通数字滤波器阶数和截止频率display(数字滤波器边界频率 单位：*pi rad);ws1 ,wp1 , wp2 ,ws2display(滤波器阶数);ndisplay(0.5dB 截止频率 单位：*pi rad);wcb,a = cheby1(n,Rp,wc);% 计算切比雪夫1型带通数字滤波器的传输函数H（z）display(H(z)分子多项式系数：);bdisplay(H(z)分母多项式系数);ahn=impz(b,a); %计算单位冲激响应H,w=freqz(b,a);%画出系统的幅频特性曲线subplot(3,1,
29、1);plot(12.5*w/pi,20*log10(abs(H),LineWidth,2);%xlabel(f(khz);ylabel(dB);title(切比雪夫1型带通数字滤波器的幅度响应（db）);axis(0,12.5 ,-70,3);grid; %画出系统的相频特性曲线subplot(3,1,2);plot(w/pi,angle(H),LineWidth,2);%xlabel(w/);ylabel(rad);title(切比雪夫1型带通数字滤波器的相位响应);grid; %画出单位冲激响应subplot(3,1,3);stem(hn,.,Linewidth,2);title(系统单位冲激响应);xlabel(n);axis(0,120,-0.2
