实验7 用MATLAB设计IIR数字滤波器
一、实验目的：
1、加深对IIR数字滤波器的基本设计方法的理解。 2、掌握用模拟滤波器原型设计IIR数字滤波器的方法。
3、了解MATLAB有关IIR数字滤波器设计的子函数的调用方法。 二、实验原理：
1、双线性变换法的基本知识
双线性变换法是将整个s平面映射到z平面，其映射关系为
1+sT/221-z-1z=s= 或 -11-sT/2T1+z双线性变换法克服了脉冲响应不变法从s平面到z平面的多值映射的缺点，消除了频谱
混叠现象。但其在变换过程中产生了非线性畸变，在设计IIR数字滤波器的过程中需要进行一定的修正。
用双线性变换法设计IIR数字滤波器的步骤如下： ① 输入给定的数字滤波器的设计指标；
② 根据公式Ω=(2/T)tan(ω/2)进行预修正，将数字滤波器设计指标转换为模拟滤波器设计指标；
③ 确定模拟滤波器的最小阶数和截止频率； ④ 计算模拟低通原型滤波器的系统传递函数；
⑤ 利用模拟域频率变换法求解实际模拟滤波器的系统传递函数； ⑥ 用双线性变换法将模拟滤波器转换为数字滤波器。 2、用双线性变换法设计IIR数字低通滤波器
例8-1 设计一个巴特沃斯数字低通滤波器，要求：ωp=0.25П，Rp=1dB；ωs=0.4П，As=15dB，滤波器采样频率Fs=100Hz。
程序清单如下：
wp=0.25*pi; %滤波器的通带截止频率 ws=0.4*pi; %滤波器的阻带截止频率 Rp=1;As=15; %滤波器的通阻带衰减指标 ripple=10^(-Rp/20); %滤波器的通带衰减对应的幅度值 Attn=10^(-As/20); %滤波器的阻带衰减对应的幅度值 %转换为模拟滤波器的技术指标 Fs=100;T=1/Fs;
Omgp=(2/T)*tan(wp/2);%原型通带频率的预修正 Omgs=(2/T)*tan(ws/2);%原型阻带频率的预修正 %模拟原型滤波器计算
[n,Omgc]=buttord(Omgp,Omgs,Rp,As,'s') %计算阶数n和截止频率 [z0,p0,k0]=buttap(n); %设计归一化的巴特沃思模拟滤波器原型 ba1=k0*real(poly(z0)); %求原型滤波器的系数b aa1=real(poly(p0)); %求原型滤波器的系数a [ba,aa]=lp2lp(ba1,aa1,Omgc); %变换为模拟低通滤波器
%也可将以上4行替换为[bb,aa]=butter(n,Omgc,'s');直接求模拟滤波器系数 %用双线性变换法计算数字滤波器系数 [bd,ad]=bilinear(ba,aa,Fs)
[sos,g]=tf2sos(bd,ad) %转换成级联型 %求数字系统的频率特性 [H,w]=freqz(bd,ad);
dbH=20*log10((abs(H)+eps)/max(abs(H))); subplot(2,2,1);plot(w/pi,abs(H));
ylabel('|H|');title('幅度响应');axis([0,1,0,1.1]); set(gca,'XTickMode','manual','XTick',[0,0.25,0.4,1]); set(gca,'YTickMode','manual','YTick',[0,Attn,ripple,1]);grid subplot(2,2,2);plot(w/pi,angle(H)/pi); ylabel('\\phi');title('相位响应');axis([0,1,-1,1]); set(gca,'XTickMode','manual','XTick',[0,0.25,0.4,1]); set(gca,'YTickMode','manual','YTick',[-1,0,1]);grid subplot(2,2,3);plot(w/pi,dbH);title('幅度响应(dB)'); ylabel('dB');xlabel('频率(\\pi)');axis([0,1,-40,5]); set(gca,'XTickMode','manual','XTick',[0,0.25,0.4,1]); set(gca,'YTickMode','manual','YTick',[-50,-15,-1,0]);grid subplot(2,2,4);zplane(bd,ad);
axis([-1.1,1.1,-1.1,1.1]);title('零极点图'); 程序运行结果如下： n = 5 Omgc = 103.2016
bd = 0.0072 0.0362 0.0725 0.0725 0.0362 0.0072 ad = 1.0000 -1.9434 1.9680 -1.0702 0.3166 -0.0392 sos =
1.0000 0.9956 0 1.0000 -0.3193 0 1.0000 2.0072 1.0072 1.0000 -0.6984 0.2053 1.0000 1.9972 0.9973 1.0000 -0.9257 0.5976 g = 0.0072
频率特性如图8-1所示：
幅度响应10.89131相位响应|H|00.1778000.250.4幅度响应(dB)0-111-100.250.4零极点图1Imaginary Part00.250.4频率(?)1?dB-150.50-0.5-1-10Real Part1
图8-1
由频率特性曲线可知，该设计结果再通阻带截止频率处能满足Rp≤1dB、As≥15dB的设计指标要求，系统的极点全部在单位圆内，是一个稳定系统。由n=5可知，该滤波器是一个5阶系统，原型Ha(s)在s=-∞处有5个零点，映射到z=-1处。该滤波器的传递函数为
0.0072+0.0362z-1+0.0725z-2+0.0725z-3+0.0362z-4+0.0072z-5H(z)=(直接型)-1-2-3-4-51-1.9434z+1.9680z-1.0702z+0.3166z-0.0392z -150.0072(1+z)H(z)=（级联型）-1-1?2-1-2(1-0.3193z)(1-0.6984z?0.2053z)(1-0.9257z+0.5976z)3、用双线性变换法设计IIR数字高通滤波器
例8-2 设计一个椭圆数字高通滤波器，要求：通带fp=250Hz，Rp=1dB；阻带fs=150Hz，
As=20dB，滤波器采样频率Fs=1000Hz。
程序清单如下：
fs=150;fp=250;Fs=1000;T=1/Fs;
wp=fp/Fs*2*pi; %滤波器的通带截止频率 ws=fs/Fs*2*pi; %滤波器的阻带截止频率 Rp=1;As=20; %滤波器的通阻带衰减指标 ripple=10^(-Rp/20); %滤波器的通带衰减对应的幅度值 Attn=10^(-As/20); %滤波器的阻带衰减对应的幅度值 %转换为模拟滤波器的技术指标
Omgp=(2/T)*tan(wp/2);%原型通带频率的预修正 Omgs=(2/T)*tan(ws/2);%原型阻带频率的预修正 %模拟原型滤波器计算
[n,Omgc]=ellipord(Omgp,Omgs,Rp,As,'s') %计算阶数n和截止频率 [z0,p0,k0]=ellipap(n,Rp,As); %设计归一化的椭圆滤波器原型 ba1=k0*real(poly(z0)); %求原型滤波器的系数b aa1=real(poly(p0)); %求原型滤波器的系数a [ba,aa]=lp2hp(ba1,aa1,Omgc); %变换为模拟高通滤波器 %用双线性变换法计算数字滤波器系数 [bd,ad]=bilinear(ba,aa,Fs) %求数字系统的频率特性 [H,w]=freqz(bd,ad);
dbH=20*log10((abs(H)+eps)/max(abs(H))); subplot(2,2,1);plot(w/2/pi*Fs,abs(H),'k'); ylabel('|H|');title('幅度响应');axis([0,Fs/2,0,1.1]); set(gca,'XTickMode','manual','XTick',[0,fs,fp,Fs/2]); set(gca,'YTickMode','manual','YTick',[0,Attn,ripple,1]);grid subplot(2,2,2);plot(w/2/pi*Fs,angle(H)/pi*180,'k'); ylabel('\\phi');title('相位响应');axis([0,Fs/2,-180,180]); set(gca,'XTickMode','manual','XTick',[0,fs,fp,Fs/2]); set(gca,'YTickMode','manual','YTick',[-180,0,180]);grid subplot(2,2,3);plot(w/2/pi*Fs,dbH);title('幅度响应(dB)'); ylabel('dB');xlabel('频率(\\pi)');axis([0,Fs/2,-40,5]);
相关推荐：
loading