北理工数字信号处理实验3IIR数字滤波器设计.docx

《北理工数字信号处理实验3IIR数字滤波器设计.docx》由会员分享，可在线阅读，更多相关《北理工数字信号处理实验3IIR数字滤波器设计.docx（16页珍藏版）》请在得力文库 - 分享文档赚钱的网站上搜索。

1、实验3 IIR数字滤波器设计一、实验目的1、掌握利用脉冲响应不变法和双线性变换法设计IIR数字滤波器的原理及具体方法2、加深理解数字滤波器和模拟滤波器之间的技术指标转化3、掌握脉冲响应不变法和双线性变换法设计IIR数字滤波器的优缺点及适用范围二、实验设备与环境计算机、MATLAB软件环境三、实验理论基础IIR谑波器设计方法由己有的模拟淀波器得到数字造波器，我们将这些模拟漉波器称为 原型漉波常用的模拟原型漉波器有巴特沃斯(Butterworth)淀波器、切比雪夫(Chebyshev, I型和口型)漉波密和椭圆(Ellipse)谑波器等。2.1巴特沃斯滤波器N阶巴特沃斯低通谑波器的幅度平方响应为(

2、A1-3)其中，N为谑波器阶数，Q,悬截上频率(单位：rad9.对于式A1-3有几点需要 说明：对于所有的N,当Q = 0时，(2)对于所有的N,当时，这说明在Q处有3dB的衰 减。(3)风(/Q)是Q的单调递减国数。 Nf 8时，氏为趋于理想的低通滤波器(5) p/qC/Q)/在Q = 0处取得最大值，因为在这里各阶导数存在且为3我们日式A1-4来确定系统函数(s)儿玛(一5)=风(通吐”(A1-4)日式A13可以得到，分母多项式的极点为-150 00.51o o 5 o -1 号(cm-5 oOmega(n)10.5Omega(兀)oo 15 0 5 o o 3(cdho(3)切比雪夫2W

3、n=0.3*pi;As=15;N=6;b,a=cheby2(N,As,Wn,s,);bz,az=impinvar(b,a);H,w=freqz(bz,az);subplot(221);plot(w/pi,abs(H);grid on;xlabel(Omega(pi),);ylabel(|H(jOmega)r);subplot(222);plot(w/pi,20*logl0(abs(H)/max(abs(H);grid on;xlabel(,Omega(pi),);ylabel(|H(jOmega)|,dB,);subplot(223);plot(w/pi,angle(H)/pi);grid o

4、n;xlabel(Omega(pi),);ylabel(,Phase of H(jOmega)(pi),);t=0:0.01:30;h=impulse(b,a,t);subplot(224);plot(t,h);grid on;xlabel(,t(s),);ylabel(Impulse Response1);00.51Omega(7r)00.51Omega(7r) suods cc s-ndlu-COPWCD工(4)椭圆Wn=0.2*pi;Rp=l;As=15;N=6;b,a=ellip(N,Rp,As,Wn/s);bz,az=impinvar(b,a);| H,w|=freqz(bz,az)

5、;subplot(221);plot(w/pi,abs(H);grid on;xlabel(Omega(pi),);ylabel(,|H(jOmega)r);subplot(222);plot(w/pi,20*log 10(abs(H)/max(abs(H);grid on;xlabel(,Omega(pi),);ylabel(|H(jOmega)|,dB,);subplot(223);plot(w/pi,angle(H)/pi);grid on;xlabel(Omega(pi),);ylabel(,Phase of H(jOmega)(pi),);t=0:0.01:30;h=impulse(

6、b,a,t);subplot(224);plot(t,h);grid on;xlabel(t(s),);ylabel(Impulse Response1);00.51o o O 12 3 _ _ _ spwcdh-5* o-40 suods CC s-ndlu_5 m /( O。0.5Omega（兀）1OO 15 0 5 o O s(crO SPLid-10五、实验心得与体会通过本次实验，掌握了利用脉冲响应不变法和双线性变换法设计IIR数字滤波器的原理及具体方 法；加深理解了数字滤波器和模拟滤波器之间的技术指标转化；掌握了脉冲响应不变法和双线性变换 法设计IIR数字滤波器的优缺点及适用范围，M

7、ATLAB操作更加熟练。白式Al-5可知/4(-s)有2N个极点，等间隔分布在半径为Q.的圆上，其角度间隔为产/N,且极点不会落在虚轴上。我们选择左半平面的极点，得到一个稳定、因果的 谑波耨的(Al-6)QvHs) = fJt-1在MATLAB中，设计阶数为N,截止频率为Q,的巴特沃斯模拟原型滤波器可以逋过 两种方法实现：方法1；在一般设计中,先把Q。选择为1厂ad/s,使频率得到归一化，归一化后的巴 特沃斯谑波器系统函数的极点分布及分母多项式系数等都有现成的表格可查n然后将归一化 巴特沃斯流波器系统函数中的变量S用S/Q,替换以后，即可得到任意Q,的非归一化的巴 特沃斯港波器口由于巴特沃斯谑

8、波器系统函数不存在零点，则将归一化巴特沃斯滤波器系统 函数的极点倍乘以Q,增益倍乘以Q.即可实现上述由归一化巴特沃斯造波器到非归一化 巴特沃斯漉波器的变换。MATLAB提供了函数乙.幻=buttap(N)来设计N阶归一化巴特沃斯模拟原型滤波器， 它返回数阻z和p(零点和积点)以及增益k.方法2：在MATLAB中还可以直接利用butter的数设计巴特沃斯模拟原型漉波系butter 函数段可用于设计模拟巴特沃斯港波it也可,用于设计数字巴特沃斯滤波器，用于设计模拟巴 特沃斯谑波维时，其茨用格式为b,a = butterfNn/s) N表示滤波器的阶数，Wn表示截止频率，父表示设计模拟谑波巴特沃斯模

9、拟原型原波器设计方程：模拟低通滤波器是白参数%, Q.和A描述的，设计巴特沃斯滤波器其实就悬 根据这些参数求得港波盘阶数N和截止频率Q小为了满足设计要求，我们希望Q = Q0时,-101og10pfflQQ)|2=,即-101os10 j iu-101os10 j iu=Rp(Al-7)并且Q = Q门时，-101吗0风=a，即-101og10=A, *(Al-8)1+求解上述两个方程，得到N =(Al-9)一般情况下，上式中的A?不为整数，但N应该为整数，所以我们取N为比实际值大的最接近实际值的整数，例如N=43 则取N=5,我们用b表示选取比x大的最小整数，则白于A，比实际需要的要大，所以

10、在或处将满足或超过指标，为了在处满足 指标，则Q, =_Ko% 一1为了在Q口处满足指标Q,=Q,=喻0%-12.2切比雪夫低通浊波器切比雪夫淀波就有两种，切比雪夫I型谑波器在通带中具有等波动响应，而切比雪夫n 型滤波塞在阻带中具有等波动响应.相比于单道特性的漉波雪，选择等波动特性的漉波器, 可以得到较低阶数的避波器，例如巴特沃斯滤波塞在逋带和阻带上的响应均是单道的，因此 对于相同的指标，切比雪夫滤波益的阶数比巴特沃斯滤波熊的阶数低。切比雷夫I型滤波器的幅度平方晌应为其由A，为谑波器阶数，Q,为截止频率，为逋带波动系数，它与R,的关系如式A1-2 所示，J，（x）是N阶切比雪夫多项式Gv（X）

11、=Gv（X）=cos(Acos x)chNchAx)确定系统函数Ha（s）的过程与前面巴特沃斯谑波器里介绍的方法类似，先求出式5）式-5）的极点，然后把左半平面的极点分配给通过求下式的根得到玛日/-0的极点如果Pi =+ J。9左= OJ,N 1是上式的根（左半平面）,则/小、 n (2左+ 1)开q =(aQJcos + -2 2Nc /工八、乃(2左+ 1)7Q. =(6Q)sin + -*2 2N其中4 =；版)1 I f;其中 口 = 一+ 1 + F般+酝 丫铲这些根落在长轴为6Q,短轴为aQ,的椭圆上，最后可以得至1系统函数为Hs) =Hs) =n（sf）其中K是归一化因子，可以由

12、下式决定N为奇数N为偶数N为奇数N为偶数(A1-19)1区(川)=1 5+在MATLAB中，设计阶数为N,截止频圭为Q。，通带波动为&的切比雪夫I型模 拟原型滤波器可以通过两种方法实现：方法;MATLAB提供了函数zpk = cheblap（NRp）来设计一个阶数为N ,逋带波动 为及的归一化（数止频率为1）切比雪夫I型模拟原型漉波若，它返回数阻z和p （妾点和 积点）以及增益匕切比雪夫I型流波器系统函数不存在零点，将归一化切比雪夫I型漉波 器系统函数的极点倍乘以Q,增益倍乘以s = 0时非归一化漉波器系统函数分母多项式和归一化滤波器系统函数多项式的比值,即可实现白归一化切比雪夫I型谑波器到非

13、归一化 切比雪夫I型滤波器的变换n方迭2_在位71，.包中还可以直接利用chebyl函数设计切比雪夫I型模拟原型谑波器, chebyl函数既可用于设计模拟切比雪夫I型谑波境也可用于设计数字切比雪夫I型谑波器, 用于设计模拟切比雪夫I型波器时，其遁用格式为b:a = cheby 1 （N?Rp; Wil s） N表示谑波器的阶数，Rp表示通带波动，Wn表示截止频 率，s表示设计模拟滤波器。切比雪夫I型模拟原型漉波器设计方程: 给定参数Q,，R。, Q.和4,设计切比雪夫I型滤波舒=首先根据式A1-2可以得到 = J10% - 14 = 10%对于切比雪夫I型滤波泰有根据式A143,由过零起点Q:

14、r处幅度平方函数的关系可以求得N 0arccoarccoarccos 也/%)(Al-22)切比切夫n型漉波器有一个单个的通带和等波动的阻带,这就意味着这种谑波器在s 平面上既有零点又有极点，因此通带中的群延迟特性比切比雪夫I型滤波恭要好。切比雪夫 n型建波器的幅度平方里应为11(A1-23)在MATLAB中，设计阶数为N,截止频率为。一阻帝衰减为4的切比雪夫口型模拟 tr,1型漉波琴可以通过两种方法实现:1型漉波琴可以通过两种方法实现:方法L MATLAB提供了函数z：p,k = cheb2即(N,As)来设计一个阶数为N,阻带衰减为过的归一化(截止频率为1)切比雪夫II型模拟原型淀波器，它

15、返叵数阻z和p (零点和 积点)以及甯益匕 切比雪夫H里谑波器系统函数存在零点，可以将归一化切比雪夫II型滤 波器系统函数的极点和零点倍乘以Q. 增益倍乘以s = 0时非归一化的有理函数和归一化的有理函数的比值,却可实现白归一化切比雪夫II型港波器到非归一化切比雪夫II型港波器 的变换.方法2：在MATLAB中还可以直接利用cheby2函数设计切比雪夫II型模拟原型谑波器. Cheby2函数SE可用于设计模拟切比雪夫D型滤波器也可用于设计数字切比雪夫n型滤波 劈，用于设计模拟切比雪夫口型波器时，其建用格式为b：a = cheby2(N,AsJVn；s) N表示涯波器的阶数，As表示过帝衰减，W

16、n表示截止频 率，表示设计模拟滤波器：切比雪夫n型模拟模型原波器设计方程除了=。二：之外，与切比雪夫I型类似，这是因为波动在阻带中。2.3椭圆漉波器椭圆滤波琴的逋带和在帝均具有等波动里应，因此相对于巴特沃斯谑波器和切比雪夫途 波罂，对于给定相司的指标，它可使阶数最小，换言之对于给定的阶数N,它能使过渡带 最陡。椭圆遮波器的幅度平方响应为|七=Al-24）1 + 识其中N为谑波器阶数，Q,为截止频率，为通带波动系数，它与的关系如式A1-2 所示，L1（x）是N阶雅可比（Jacobian）椭圆函数n在MATLAB中，设计阶数为N,截止频率为Q,逋笨波动为R”阻带衰减为4的 椭圆模拟原型滤波釜可以通

17、过两种方法实现：方法：MATLAB提供了函数乙p,k = ellipap（N,Rp,As）来设计一个阶数为N ,通带波动为火一阻带衰减为乱的归一化（截止频率为1）椭圆模拟原型滤波器，它返回数里z和p（零点和积点）以及年基k将归一化椭圆漉波器系统函数的极点和零点倍乘以增益 倍乘以s = Q时非归一化的有理函数和归一化的有理函数的比值,即可实现自归一化椭圆港波景到非归一化插圆滤波器的变换,方法2_左MATLAB中还可以直接利用ellip函数设计椭圆模拟原型涯波器。ellip函数正可用于设计模拟椭圆港波琴也可.用于设计数字椭圆谑波泰，目于设计模拟椭圆波泰时，其 追月格式为b:a = ellip（N:

18、Rp;As:Wn/s,） N表示淀波器的阶数，Rp表示通波动，As表示阻带衰 减，Wn表示欲止频率，管表示设计模拟遮波器。(A1-25)(A1-26)对于椭圆谑波若的分析垓困荒，在这里我们不作具体分析，但是桶圆模拟原型滤波器设 让方程在很多文献上都己经给出了K（左）K（正左；）K6）K（J1-二）其中K（x）为完全I型椭圆积分，MATLAB提供了函数ellipke用于计算该积分,四、实验内容设采样频率为=10狂Zz,设计数字低通淀波器，满足如下指标通带截止频率：f=lkHz,通带波动：RP = dB阻带截止频率：九=L5kHz,阻笔衰减：4=15dB 要求分别设计巴特沃斯、切比雪夫I型、切比雪

19、夫II型和椭圆模拟原型速波器，并分别结合 脉冲响应不变法和双线性变换法进行设计.结合实验结果，分别讨论采用上述方法设计的数 字漉波器是否都能满足给定指标要求，分析脉冲响应不变法和双线性变换法设计IIR数字谑 波器的优景点及适用范围n解：(1)巴特沃斯Wn=0.2*pi;Rp=l;As=15;N=6;b,a=butter(N, Wn/s*);bz,az=impinvar(b,a);H,w=freqz(bz,az);subplot(221);plot(w/pi,abs(H);grid on;xlabel(,Omega(pi);ylabel(|H(jOmega)r);subplot(222);plo

20、t(w/pi,20*log 10(abs(H)/max(abs(H);grid on;xlabel(,Omega(pi),);ylabel(|H(jOmega)|,dB,);subplot(223);plot(w/pi,angle(H)/pi);grid on;xlabel(,Omega(pi),);ylabel(,Phase of H(jOmega)(pi);t=0:0.01:30;h=impulse(b,a,t);subplot(224);plot(t,h);grid on;xlabel(,t(s),);ylabel(,lmpulse Response);0.5Omega(兀)150一 G

21、Q一(cm一5 o-100 00.51Omega(兀)oo 15 0 5 o o 3(cdho-10(2)切比雪夫1Wn=0.2*pi;Rp=l;N=6;b,a=cheby 1 (N,Rp,Wn,s)；bz,az=impinvar(b,a);H,w=freqz(bz,az);subplot(221);plot(w/pi,abs(H);grid on;xlabel(Omega(pi),);ylabel(|H(jOmega)r);subplot(222);plot(w/pi,20*logl0(abs(H)/max(abs(H);grid on;xlabel(,Omega(pi),);ylabel(|H(jOmega)|,dB,);subplot(223);plot(w/pi,angle(H)/pi);grid on;xlabel(Omega(pi),);ylabel(,Phase of H(jOmega)(pi),);t=0:0.01:30;h=impulse(b,a,t);subplot(224);plot(t,h);grid on;xlabel(,t(s),);ylabel(Impulse Response1);