第6章无限脉冲响应数字滤波器的设计ppt课件.ppt

第6章 无限脉冲响应数字滤波器的设计,6.1 数字滤波器的基本概念6.2 模拟滤波器的设计6.3 用脉冲响应不变法设计IIR数字低通滤波器6.4 用双线性变换法设计IIR数字 低通滤波器,6.1 数字滤波器的基本概念,本章讲述滤波器的基本概念，讲述无限脉冲响应数字滤波器（低通、高通、带通和带阻）的设计方法转换法。,例1,0.滤波器的基本概念,设,分贝(dB)的定义:,当,-3dB带宽:,例2,系统函数为,频率函数为,1.数字滤波器的分类,(1)经典滤波器与现代滤波器,(2)低通,高通,带通和带阻,(3)无限脉冲响应(IIR)滤波器和有限脉冲响应(FIR)滤波器。,数字滤波器的低频频带位于2的整数倍处，高频频带位于的奇数倍附近，,IIR和FIR 系统函数分别为：,(6.1.1),(6.1.2),(NM)称为N阶IIR滤波器函数,称为N-1阶FIR滤波器函数,2.数字滤波器的技术要求,(1)理想滤波器,图6.1.1 理想低通、高通、带通、带阻滤波器幅度特性,带通,带阻,假设数字滤波器的频率函数H(e j)用下式表示：,()为相频特性，表示通过滤波器后各频率成分在时间上的延迟情况。,(2)物理可实现滤波器,为幅频特性，表示通过滤波器后各频率成分衰减情况。,p 通带截止频率 c 3dB通带截止频率 s 阻带截止频率,图6.1.2 低通滤波器的技术要求,通带(0,p)内要求,阻带(s,)内要求,p 到s为过渡带。,例如：某滤波器的频率特性为,p和s分别定义为：,通带内允许的最大衰减用p表示，阻带内允许的最小衰减用s表示.,如将|H(ej0)|归一化为1，(6.1.3)和(6.1.4)式则表示成：,(6.1.5),(6.1.6),与理想低通比较:,设计步骤：先设计模拟滤波器得到传输函数Ha(s)，然后将Ha(s)按某种方法转换成数字滤波器的系统函数H(z)。,3.数字滤波器设计方法概述,IIR滤波器设计方法,经常用的一类设计方法是借助于模拟滤波器的设计方法进行的。,另一类是直接在频域或时域中进行设计。,FIR滤波器的设计方法,常用的有窗函数法和频率采样法。,6.2 模拟滤波器的设计,6.2.1模拟低通滤波器的设计指标及逼近方法,(1)技术指标,模拟低通滤波器的设计指标有p,p,s和s。,图6.2.2 低通滤波器的幅频特性,p和s分别称为通带截止频率和阻带截止频率.,损耗函数：,用对数表示幅频特性,图中c称为3dB截止频率。,p是通带(=0p)中的最大衰减系数，s是阻带s的最小衰减系数，p和s一般用dB数表示。,损耗函数曲线：,对于单调下降的幅度特性，可表示成：,如果=0处幅度已归一化到1，即|Ha(j0)|=1,p和s表示为：,由给定的指标p,p,s和s 求出,然后得到Ha(s).,(2)逼近方法,称为幅度平方函数.,而一般滤波器的单位冲激响应h(t)为实数且收敛，因此,为保证系统的稳定,的极点应位于S平面的左半平面,的极点应位于S平面的右半平面。,具体做法：,因为,将 中位于左半平面的极点组成。,巴特沃斯低通滤波器的幅度平方函数|Ha(j)|2用下式表示：,(6.2.6),6.2.2 巴特沃斯低通滤波器的设计方法,N称为滤波器的阶数.,(1)巴特沃斯低通滤波器的传输函数Ha(s),1.巴特沃斯低通滤波器的设计原理,图6.2.3 巴特沃斯幅度特性和N的关系,将幅度平方函数,(6.2.7),此式表明幅度平方函数有2N个极点.,写成s的函数：,得,极点sk用下式表示：,(6.2.8),为形成稳定的滤波器，2N个极点中只取s平面左半平面的N个极点构成Ha(s)，而右半平面的N个极点构成Ha(-s)。Ha(s)的表示式为,比如N=3,通过下式可以计算出6个极点.,同理，得：,图6.2.4 三阶巴特沃斯滤波器极点分布,当N=3时，6个极点中位于左半平面的三个分别为:,取s平面左半平面的极点s0,s1,s2组成Ha(s)：,(2)归一化传输函数Ha(p),将 对3dB截止频率c归一化；,(6.2.10),归一化后的Ha(s)表示为,=/c，称为归一化频率。,上式中，,令,p称为归一化拉氏复变量.,式中，pk=sk/c为归一化极点，sk为位于左半平面的极点;用下式表示：,(6.2.11),归一化巴特沃斯的传输函数为,（6.2.12）,将极点表示式(6.2.12)代入(6.2.11)式，得到的Ha(p)的分母是p的N阶多项式，用下式表示：,表6.2.1 巴特沃斯归一化低通滤波器参数,续表6.2.1,续表6.2.1,(3)N的确定,由技术指标p,p,s和s确定.,整理得：,由(6.2.14)和(6.2.15)式得到：,令,则N由下式表示：,(6.2.16),取大于等于N的最小整数。,关于3dB截止频率c，如果技术指标中没有给出，可以按照(6.2.14)式或(6.2.15)式求出，,(4)总结低通巴特沃斯滤波器的设计步骤如下：,根据技术指标p,p,s和s，用 求出滤波器的阶数N。,按照，求出归一化极点pk，将pk代入，得到归一化传输函数Ha(p)。,将Ha(p)去归一化。将p=s/c代入Ha(p)，得到实际的滤波器传输函数Ha(s)。,例6.2.1 已知通带截止频率fp=5kHz，通带最大衰减p=2dB，阻带截止频率fs=12kHz，阻带最小衰减s=30dB，按照以上技术指标设计巴特沃斯低通滤波器。,解(1)确定阶数N。,技术指标:fp=5kHz p=2dB fs=12kHz s=30dB 代入,取 N=5,得到：,(2)确定归一化传输函数.,由N=5，直接查表6.2.1，得到：极点：-0.3090j0.9511,-0.8090j0.5878，-1.0000,将共轭极点放在一起,形成因式分解形式:,上式分母也可以展开成为五阶多项式,式 中 b0=1.0000,b1=3.2361,b2=5.2361,b3=5.2361,b4=3.2361,先求3dB截止频率c。按照(6.2.17)式，得到：,(3)将Ha(p)去归一化,将p=s/c代入Ha(p)中得到：,例6.2.1图1 所设计滤波器的幅频特性,close all,N=256;t=linspace(0,1,N);dt=t(2)-t(1);xt=cos(2*pi*4*t)+cos(2*pi*10*t)+cos(2*pi*20*t);f=(0:(N/2-1)/(dt*N);Xt=fft(xt,N);subplot(2,1,1),plot(t(1:128),xt(1:128),subplot(2,1,2),plot(f(1:64),abs(Xt(1:64),xlabel(f(kHz)Qc=5.2775;b0=1;b1=3.2361;b2=5.2361;b3=b2;b4=b1;Q=f%Q=linspace(0,25,N);Ha=Qc5./(j*Q).5+b4*Qc*(j*Q).4+b3*Qc2*(j*Q).3+b2*Qc3*(j*Q).2+b1*Qc4*(j*Q)+b0*Qc5);L=length(Ha)Has=20*log10(abs(Ha);figure,plot(Q(1:64),Has(1:64),Q,-30,r*,12,Has,*,5,Has,*),axis(0 30-70 2),xlabel(f(kHz),ylabel(20lg(abs(H_a(jOmega)(dB);Yt=Xt(1:L).*Ha;yt=ifft(Yt);figure,subplot(2,1,1),plot(t(1:128),abs(yt(1:128),subplot(2,1,2),plot(f(1:64),abs(Yt(1:64);,xlabel(f(kHz),附例题的绘图程序,例6.2.1图2 假定的输入信号的时域波形和频谱,例6.2.1图3 滤波器的输出信号的时域波形和频谱,(4)简单讨论,将c代入(6.2.18)式,得到：,第一,此值小于题目给定的值,说明设计出的滤波器在fs=12kHz处的衰减,超过30dB,过渡带小于要求的.,fp=5kHz p=2dB fs=12kHz s=30dB,要求设计的滤波器,设计的滤波器一,阻带性能比设计要求好.,第二,如果用 计算c,则,用 计算p,则,此值大于题目给定的值,这样设计出的滤波器在fp=5kHz处的衰减,小于2dB.,设计的滤波器二,fp=5kHz p=2dB fs=12kHz s=30dB,要求设计的滤波器,通带性能比设计要求好.,2.用MATLAB工具箱函数设计巴特沃斯滤波器,MATLAB信号处理工具箱函数buttap,buttord和butter是巴特沃斯滤波器设计函数。,调用格式如下。,Z,P,K=buttap(N)%计算归一化（以3dB截止频率归一化）系统函数的零极点和增益。,N,wc=buttord(wp,ws,Rp,As,s),B,A=butter(N,wc,s),B,A=zp2tf(Z,P,K),式中，Z(k)和P(k)分别为向量Z和P的第k个元素。,1）Z,P,K=buttap(N)该格式用于计算N阶巴特沃斯归一化（3 dB截止频率c=1）模拟低通原型滤波器系统函数的零、极点和增益因子。得到的系统函数为如下形式:（6.2.21）,Z,P,K=buttap(5)Z=P=-0.3090+0.9511i-0.3090-0.9511i-0.8090+0.5878i-0.8090-0.5878i-1.0000 K=1,B,A=zp2tf(Z,P,K)B=0 0 0 0 0 1A=1.0000 3.2361 5.2361 5.2361 3.2361 1.0000,2）N,wc=buttord(wp,ws,Rp,As)该格式用于计算巴特沃斯数字滤波器的阶数N和3 dB截止频率wc。调用参数wp和ws分别为数字滤波器的通带边界频率和阻带边界频率的归一化值，要求0wp1，0ws1,1表示数字频率（对应模拟频率Fs/2，Fs表示采样频率）。Rp和As分别为通带最大衰减和阻带最小衰减（dB）。当wswp时，为高通滤波器；当wp和ws为二元矢量时，为带通或带阻滤波器，这时wc也是二元向量。,3）N,wc=buttord(wp,ws,Rp,As,s)该格式用于计算巴特沃斯模拟滤波器的阶数N和3 dB截止频率wc。wp、ws和wc是实际模拟角频率（rad/s）。其他参数与格式2）相同。,wp=2*pi*5000;ws=2*pi*12000;Rp=2;As=30;N,wc=buttord(wp,ws,Rp,As,s)N=5wc=3.7792e+004,4）B,A=butter(N,wc,ftype)计算N阶巴特沃斯数字滤波器系统函数分子和分母多项式的系数向量B和A。调用参数N和wc分别为巴特沃斯数字滤波器的阶数和3 dB截止频率的归一化值（关于归一化），一般按格式2）调用函数buttord计算N和wc。由系数向量B和A可以写出数字滤波器系统函数：,（6.2.22）,式中，B(k)和A(k)分别为向量B和A的第k个元素。,5）B,A=butter(N,wc,ftype,s)计算巴特沃斯模拟滤波器系统函数的分子和分母多项式的系数向量B和A。调用参数N和wc分别为巴特沃斯模拟滤波器的阶数和3 dB截止频率（实际角频率）。由系数向量B和A写出模拟滤波器的系统函数为（6.2.23）,（6.2.21）、（6.2.22）和（6.2.23）式也适用于后面要介绍的切比雪夫和椭圆滤波器的MATLAB设计函数。,B,A=butter(N,wc,s)B=1.0e+022*0 0 0 0 0 7.7094A=1.0e+022*0.0000 0.0000 0.0000 0.0000 0.0007 7.7094,Z,P,K=buttap(5);B,A=zp2tf(Z,P,K)B=0 0 0 0 0 1A=1.0000 3.2361 5.2361 5.2361 3.2361 1.0000,由于高通滤波器和低通滤波器都只有一个3 dB截止频率wc，因此仅由调用参数wc不能区别要设计的是高通还是低通滤波器。当然仅由二维向量wc也不能区分带通和带阻。所以用参数ftype来区分。ftype=high时，设计3 dB截止频率为wc的高通滤波器。缺省ftype时默认设计低通滤波器。ftype=stop时，设计通带3 dB截止频率为wc的带阻滤波器，此时wc为二元向量wcl,wcu，wcl和wcu分别为带阻滤波器的通带3 dB下截止频率和上截止频率。缺省ftype时设计带通滤波器，通带为频率区间wclwcu。应当注意，设计的带通和带阻滤波器系统函数是2N阶的。这是因为带通滤波器相当于N阶低通滤波器与N阶高通滤波器级联。,ftype选项：,【例6.2.2】调用buttord和butter设计巴特沃斯低通模拟滤波器。要求与例6.2.1相同。,设计程序如下:wp=2*pi*5000;ws=2*pi*12000;Rp=2;As=30;%设置滤波器参数N,wc=buttord(wp,ws,Rp,As,s)%计算滤波器阶数N和3 dB截止频率 B,A=butter(N,wc,s)%计算滤波器系统函数分子分母多项式系数 k=0:511;fk=0:14000/512:14000;wk=2*pi*fk;Hk=freqs(B,A,wk);plot(fk/1000,20*log10(abs(Hk);grid onxlabel(频率(kHz);ylabel(幅度(dB),axis(0,14,-40,5)fc=wc/2/pi,运行结果:N=5，wc=3.7792e+004fc=6.0148e+003B=1.0e+022*0 0 0 0 0 7.7094A=1.0e+022*0.0000 0.0000 0.0000 0.0000 0.0007 7.7094,低通滤波器的幅频特性,将B和A代入（6.2.23）式写出系统函数为与例6.2.1计算结果形式相同。滤波器的损耗函数曲线如图6.2.6所示。由图可以看出，阻带刚好满足指标要求，通带指标有富余。这就说明buttord函数使用（6.2.20）式计算3 dB 截止频率。,6.2.3 切比雪夫滤波器的设计,1 切比雪夫滤波器的设计原理,切比雪夫滤波器有两种形式：,一、振幅特性在通带内是等波纹的、在阻带内是单调下降的切比雪夫型滤波器；二、振幅特性在通带内是单调下降、在阻带内是等波纹的切比雪夫型滤波器。采用何种形式的切比雪夫滤波器取决于实际用途。,图6.2.7(a)和(b)分别画出不同阶数的切比雪夫型和型滤波器幅频特性。,图6.2.7 不同阶数的切比雪夫型和型滤波器幅频特性,型,型,图中的归一化角频率定义为：,型：对通带截止频率归一化,型：对阻带截止频率归一化,我们这里仅介绍切比雪夫型滤波器的设计方法。,式中，为小于1的正数，表示通带内幅度波动的程度，愈大，波动幅度也愈大;p称为通带截止频率。令=/p，称为对p的归一化频率。,其幅度平方函数用|Ha(j)|2表示：,CN(x)称为N阶切比雪夫多项式，定义为,MATLAB信号处理工具箱函数cheb1ap，cheb1ord和cheby1是切比雪夫型滤波器设计函数。其调用格式如下:1）Z,P,G=cheb1ap(N,Rp)%计算归一化（以通 带截止频率归一化）系统函数的零极点和增益。2）N,wpo=cheb1ord(wp,ws,Rp,As,s)3）B,A=cheby1(N,Rp,wpo,ftype,s),2.用MATLAB设计切比雪夫滤波器,wpo是切比雪夫型滤波器的通带截止频率。,MATLAB信号处理工具箱函数cheb2ap,cheb2ord和cheby2是切比雪夫型滤波器设计函数。其调用格式如下:,1）Z,P,G=cheb2ap(N,As)%计算归一化（以阻带 截止频率归一化）系统函数的零极点和增益。2）N,wso=cheb2ord(wp,ws,Rp,As,s)3）B,A=cheby2(N,Rp,wso,ftype,s),wso是切比雪夫型滤波器的阻带截止频率。,设计切比雪夫型模拟低通滤波器wp=2*pi*3000;ws=2*pi*12000;Rp=0.1;As=60;设置指标参数N1,wp1=cheb1ord(wp,ws,Rp,As,s);计算切比雪夫型模拟低通滤波器阶数和通带边界频率B1,A1=cheby1(N1,Rp,wp1,s);计算切比雪夫型模拟低通滤波器系统函数系数fk=0:14000/512:14000;wk=2*pi*fk;Hk=freqs(B1,A1,wk);%在指定的频率点上计算频率响应plot(fk/1000,20*log10(abs(Hk);grid onxlabel(频率(kHz);ylabel(幅度(dB)axis(0,14,-70,5),【例6.2.3】设计切比雪夫型和切比雪夫型模拟低通滤波器。要求与例6.2.2相同。,运行结果:N=5切比雪夫型模拟低通滤波器通带边界频率:wp1=1.8850e+004切比雪夫型模拟低通滤波器系统函数分子分母多项式系数:B=1.2187e+011A=1 3.2873e+004 9.8445e+008 1.6053e+013 1.8123e+017 9.7448e+020,滤波器损耗函数如图6.2.9所示。,图6.2.9 五阶切比雪夫型模拟低通滤波器损耗函数（例6.2.3的设计结果）,6.2.4 椭圆滤波器的设计,椭圆滤波器的典型幅频响应特性曲线如图6.2.10所示。,椭圆滤波器的阶数N由通带边界频率p、阻带边界频率s、通带最大衰减 p和阻带最小衰减 s共同决定。,图6.2.10椭圆滤波器幅频响应特性曲线,图中的归一化角频率定义为：,对通带截止频率归一化。,MATLAB信号处理工具箱提供椭圆滤波器设计函数ellipap、ellipord和ellip。其调用格式如下:,1）z,p,k=ellipap(N,Rp,As)%计算归一化（以通带截止频率归一化）系统函数的零极点和增益。,2）N,wpo=ellipord(wp,ws,Rp,As,s),3）B,A=ellip(N,Rp,wpo,ftype,s),%椭圆滤波器设计程序wp=2*pi*3000;ws=2*pi*12000;Rp=0.1;As=60;%设置指标参数N,wpo=ellipord(wp,ws,Rp,As,s);%计算椭圆低通模拟滤波器阶数和通带边界频率B,A=ellip(N,Rp,As,wpo,s);%计算低通模拟滤波器系统函数系数fk=0:14000/512:14000;wk=2*pi*fk;Hk=freqs(B,A,wk);plot(fk/1000,20*log10(abs(Hk);grid onxlabel(频率(kHz);ylabel(幅度(dB)axis(0,14,-70,5),【例6.2.4】设计椭圆模拟低通滤波器。要求与例6.2.2相同。,运行结果:椭圆模拟低通滤波器阶数:N=4模拟低通滤波器通带边界频率:wpo=1.8850e+004椭圆模拟低通滤波器系统函数分子分母多项式系数:B=0.0010 8.3913e015 2.9126e+007 8.0051e004 1.0859e+017A=1 3.3792e+004 9.3066e+008 1.3646e+013 1.0984e+017,滤波器损耗函数如图6.2.11所示。,图6.2.11 四阶椭圆模拟低通滤波器损耗,6.2.5 四种类型模拟滤波器的比较,巴特沃思、切比雪夫型、切比雪夫型和椭圆滤波器,wp=2*pi*5000;ws=2*pi*12000;Rp=2;As=30;,滤波器阶数N1=5 N2=3 N3=3 N4=3,滤波器参数为N=5,Rp=2;As=30,高通、带通、带阻滤波器的幅频响应曲线及边界频率分别如下图所示。,6.2.6 模拟高通、带通、带阻滤波器的设计,用MATLAB实现：,N,wc=buttord(wp,ws,Rp,As,s),B,A=butter(N,wc,ftype,s),高通：ftype=high,带通：wp,ws均为二元向量，ftype缺省。,带阻：wp,ws均为二元向量，ftype=stop。,wp=2*pi*12000;ws=2*pi*5000;Rp=2;As=30;%设置滤波器参数N1,wc=buttord(wp,ws,Rp,As,s);%计算滤波器阶数N和3 dB截止频率 B,A=butter(N1,wc,high,s);%计算滤波器系统函数分子分母多项式系数,wp=2*pi*5000,10000;ws=2*pi*3000,12000;Rp=2;As=30;%设置滤波器参数N1,wc=buttord(wp,ws,Rp,As,s);%计算滤波器阶数N和3 dB截止频率B,A=butter(N1,wc,s);%计算滤波器系统函数分子分母多项式系数,ws=2*pi*5000,10000;wp=2*pi*3000,12000;Rp=2;As=30;%设置滤波器参数N1,wc=buttord(wp,ws,Rp,As,s);%计算滤波器阶数N和3 dB截止频率B,A=butter(N1,wc,stop,s);%计算滤波器系统函数分子分母多项式系数,6.3 用脉冲响应不变法设计IIR数字低通滤波器,主要讲述:将模拟滤波器Ha(s)转换为数字滤波器H(z),(1)因果稳定的模拟滤波器转换成数字滤波器，仍是因果稳定的。转换关系应实现 s平面的左半平面映射到z平面的单位圆内。,(2)数字滤波器的频率响应模仿模拟滤波器的频响.转换关系应实现 s平面的虚轴映射到z平面的单位圆，相应的频率之间成线性关系。,为了保证转换后的H(z)稳定且满足技术要求，对转换关系提出两点要求：,设模拟滤波器的传输函数为Ha(s),相应的单位冲激响应是ha(t)，由拉普拉斯变换，,设模拟滤波器Ha(s)只有单阶极点，且分母多项式的阶次高于分子多项式的阶次；,一.脉冲响应不变法,(6.3.1),式中si为Ha(s)的单阶极点。,将Ha(s)用部分分式表示：,(6.3.1),式中si为Ha(s)的单阶极点。,将Ha(s)用部分分式表示：,将Ha(s)进行逆拉氏变换得到ha(t)：,式中u(t)是单位阶跃函数。,(6.3.2),对ha(t)进行等间隔采样，采样间隔为T，得到：,(6.3.3),对上式进行Z变换，得到数字滤波器的系统函数H(z)：,(6.3.1),比较,二.标准映射关系,所以,即,由此得到标准映射关系:,(6.3.6),设,按照z=esT式，得到:,因此得到:,(6.3.10),即，当=0时,r=1;,S平面上的虚轴映射到z平面的单位圆.,当0时,r1;,当0时,r1;,S平面上的左半平面映射到z平面的单位圆内.,S平面上的右半平面映射到z平面的单位圆外.,另外，z=esT是一个周期函数，可写成,所以,可见数字频率与模拟频率是线性关系但不是一一对应的.,=0,r=1 0,r1,图6.3.1 z=esT,s平面与z平面之间的映射关系,s平面,z平面,采样信号 的频谱与原模拟信号的关系：,离散信号x(n)的频谱与原模拟信号的关系：,比较：,有,在 中，令=0（r=1），也可得到这个结论。,1、脉冲响应不变法（变换方法）,采样,若,则,2、标准映射关系（是否满足第一点要求）,3、,（是否满足第二点要求）,0,0,模拟滤波器,希望的数字滤波器,转换后,0,实际的数字滤波器,三.频谱混叠现象,设模拟滤波器的频谱为,(6.3.7),将(6.3.7)式代入上式,可得:,假设,在这种情况下,用脉冲响应不变法设计的数字滤波器可以重现原模拟滤波器的频响.,则,图6.3.2 脉冲响应不变法的频率混叠现象,脉冲响应不变法不适合设计高通滤波器和带阻。,四.滤波器的传输函数,一般Ha(s)的极点si是一个复数，且以共轭成对的形式出现.,在 中将一对复数共轭极点放在一起，形成一个二阶基本节。,相应的数字滤波器二阶基本节(只有实数乘法)的形式为,(6.3.12),如果模拟滤波器的二阶基本节的形式为,(6.3.11),如果模拟滤波器二阶基本节的形式为,(6.3.14),相应的数字滤波器二阶基本节的形式为,(6.3.13),脉冲响应不变法的优点与缺点:,优点：1、频率变换线性。2、时域特性好。,缺点：频谱混叠，不适合设计高通和带阻。,用脉冲响应不变法将Ha(s)转换成数字滤波器的系统函数H(z)。,例6.3.1 已知模拟滤波器的传输函数Ha(s)为,极点为,解 首先将Ha(s)写成部分分式：,式中T是采样间隔,T的选取应满足抽样定理.,方法一:按照，并经过整理，得到,那么H(z)的极点为,设T=1s时用H1(z)表示，T=0.1s时用H2(z)表示，则,方法二:转换时，也可以直接按照(6.3.13),(6.3.14)式进行转换。,首先将Ha(s)写成(6.3.13)式的形式，如极点 s1,2=1j1,则,再按照(6.3.14)式，H(z)为,图6.3.3 例6.3.1的幅度特性,图(a)表示模拟滤波器的幅度特性,图(b)表示数字滤波器的幅度特性,模拟频率与数字频率 满足=T,如T=0.1s(a)图中C点=15rad/s=4.77rad/s,(b)图中c点=4.77 0.1=0.477 rad,令,(a)图的纵坐标为,令,(b)图的纵坐标为,P194 6.,解:,收敛域包括虚轴。,另解：直接对 采样，,数值上,因此有,则,数字系统为低通滤波器.,五.脉冲响应不变法的MATLAB实现,Bz,Az=impinvar(B,A,Fs)%用脉冲响应不变法将模拟滤波器转换为数字滤波器,B,A分别是模拟滤波器系统函数Ha(s)分子分母多项式的系数。,Bz,Az分别是数字滤波器系统函数H(z)分子分母多项式的系数。,Fs为采样频率。,impluse invariable,%例题6.3.1的转换程序b=0.5012;a=1,0.6449,0.7079;h,w=freqs(b,a);%计算模拟滤波器的频响subplot(2,1,1),plot(w(1:180),20*log10(abs(h(1:180)/max(h),grid onxlabel(归一化角频率),ylabel(幅度),title(模拟滤波器)bz,az=impinvar(b,a,1)h,w=freqz(bz,az);%计算数字滤波器的频响subplot(2,1,2),plot(w/pi,20*log10(abs(h/max(h),grid onxlabel(以pi为单位),ylabel(幅度),title(数字滤波器),hold onbz,az=impinvar(b,a,10)h,w=freqz(bz,az);%计算数字滤波器的频响subplot(2,1,2),plot(w/pi,20*log10(abs(h/max(h),grid onxlabel(以pi为单位),ylabel(幅度),title(数字滤波器),例题6.3.1的转换结果,例：用脉冲响应不变法设计数字滤波器。要求通带和阻带具有单调下降特性，技术指标如下：,解：,要求数字滤波器的通带和阻带具有单调下降特性，可采用巴特沃思模拟滤波器作为原型；,设计步骤：1、利用关系 将数字频率转换为模拟频率；2、先设计巴特沃思模拟滤波器；3、用脉冲响应不变法将模拟滤波器转换为数字滤波器。,将数字频率转换为模拟频率（假设T=1s）；,第二步，第三步用MATLAB 实现，程序如下：,结果如下：,wp=0.2*pi;ws=0.35*pi;Rp=2;As=30;%设置滤波器参数N,wc=buttord(wp,ws,Rp,As,s);%计算滤波器阶数N和3 dB截止频率B,A=butter(N,wc,s);%计算滤波器系统函数分子分母多项式系数 NHk,wk=freqs(B,A);subplot(2,1,1);plot(wk(1:150),20*log10(abs(Hk(1:150);grid onxlabel(频率(Hz);ylabel(幅度(dB)bz,az=impinvar(B,A,1);h,w=freqz(bz,az);%计算数字滤波器的频响subplot(2,1,2),plot(w/pi,20*log10(abs(h/max(h),grid onxlabel(以pi为单位),ylabel(幅度),title(数字滤波器),如果频率指标是用数字频率给的，采样间隔可任取。,6.4 用双线性变换法设计IIR数字低通滤波器,为了克服频谱混叠的缺点,采用非线性频率压缩方法.,即:将整个频率轴上的频率范围压缩到 之间,再用 转换到z平面上.,一.非线性频率压缩,图6.4.1 双线性变换法的映射关系,当1从-/T经过0变化到/T时，则由-经过0变化到+，实现了s平面上整个虚轴完全压缩到s1平面上虚轴的/T之间的转换.。,用正切变换实现频率压缩：,(6.4.1),式中T是采样间隔.,设,经非线性频率压缩后,用 表示.,这样,由s=j和s1=j1得到:,(6.4.2),(6.4.3)(6.4.4)称为双线性变换.,将 代入,得:,由,将s平面映射为s1平面,图6.4.1 双线性变换法的映射关系,二.模拟频率和数字频率之间的关系。,令s=j,z=e j，并代入 中，有,图6.4.2 双线性变换法的频率变换关系,(6.4.5),所以,模拟频率与数字频率是一一对应的，但非线性。,图6.4.3 双线性变换法幅度和相位特性的非线性映射,三.幅度特性和相位特性的失真,设,四.数字滤波器的传输函数:,则,表6.4.1 系数关系表,小结：用双线性变换法将模拟滤波器Ha(s)转换为数字滤波器H(z):,1、转换方法：,2、频率映射关系：,3、优点：没有频谱混叠,4、缺点：频率关系非线性。,例6.4.1试分别用脉冲响应不变法和双线性不变法将图6.4.4所示的RC低通滤波器转换成数字滤波器。,解:首先按照图6.4.4写出该滤波器的传输函数Ha(s)为,图6.4.4 RC低通滤波器,利用脉冲响应不变法转换，数字滤波器的系统函数H1(z)为,利用双线性变换法转换，数字滤波器的系统函数H2(z)为,图6.4.5 例6.4.1图H1(z)和H2(z)的网络结构(a)H1(z);(b)H2(z),脉冲响应不变法,双线性变换法,图6.4.6例6.4.1 图数字滤波器H1(z)和H2(z)的幅频特性,总结:利用模拟滤波器设计IIR数字低通滤波器的步骤。,(1)确定数字低通滤波器的技术指标：通带截止频率p、通带衰减p、阻带截止频率s、阻带衰减s。,(2)将数字低通滤波器的技术指标转换成模拟低通滤波器的技术指标。,如果采用双线性变换法，边界频率的转换关系为,T可任选,=T,如果采用脉冲响应不变法，边界频率的转换关系为,具体转换：,(3)按照模拟低通滤波器的技术指标设计模拟低通滤波器Ha(s)。,(4)将模拟滤波器Ha(s)，从s平面转换到z平面，得到数字低通滤波器系统函数H(z)。,采样间隔T的选择:,（1）脉冲响应不变法要求,选择T满足;s是阻带截止频率。,（2）双线性变换法T可任选,一般选T=1.,若先给定数字滤波器的技术指标,脉冲响应不变法和双线性变换法的T可任选,一般选T=1.,例6.4.2 设计低通数字滤波器，要求在通带内频率低于0.2rad时，容许幅度误差在1dB以内；在频率0.3到之间的阻带衰减大于15dB。指定模拟滤波器采用巴特沃斯低通滤波器。试分别用脉冲响应不变法和双线性变换法设计滤波器。,解:(1)用脉冲响应不变法设计数字低通滤波器。,p=0.2rad/s,p=1dB;s=0.3rad/s,s=15dB。,模拟低通的技术指标为,选T=1s,由=T 可得,数字低通的技术指标为 p=0.2rad,p=1dB;s=0.3rad,s=15dB。,设计巴特沃斯低通滤波器。先计算阶数N及3dB截止频率c。,取N=6,此值满足通带技术要求，同时给阻带衰减留一定余量，这对防止频率混叠有一定好处。,为求3dB截止频率c,将p和p代入下式，得到,为去归一化，将p=s/c代入Ha(p)中，得到实际的传输函数Ha(s),根据阶数N=6，查表6.2.1，得到归一化传输函数为,用脉冲响应不变法将Ha(s)转换成H(z),图6.4.7 例6.4.2图用脉冲响应不变法设计的数字低通滤波器的幅度特性,wp=0.200*pi;ws=0.30*pi;Rp=1;As=15;%设置滤波器参数N,wc=buttord(wp,ws,Rp,As,s)%计算阶数N和3 dB截止频率（用阻带参数）B,A=butter(N,wc,s)%计算系统函数分子分母多项式系数 bz,az=impinvar(B,A)%脉冲响应不变法转换为数字系统,N=6wc=0.7087 B=0 0 0 0 0 0 0.1266A=1.0000 2.7380 3.7484 3.2533 1.8824 0.6905 0.1266bz=0.0000 0.0007 0.0105 0.0167 0.0042 0.0001 0az=1.0000-3.3443 5.0183-4.2190 2.0725-0.5600 0.0647,模拟低通的技术指标为,(2)用双线性变换法设计数字低通滤波器。,数字低通技术指标仍为 p=0.2rad,p=1dB;s=0.3rad,s=15dB。,设计巴特沃斯低通滤波器。,这样阻带技术指标满足要求，通带指标已经超过。,阶数N计算如下：,N=6,用 计算c,得到c=0.7662rad/s。,根据N=6，查表6.2.1得到的归一化传输函数Ha(p)与脉冲响应不变法得到的相同。为去归一化，将p=s/c代入Ha(p)，得实际的Ha(s)，,用双线性变换法将Ha(s)转换成数字滤波器H(z)：,图6.4.8 例6.4.2图用双线性变换法设计的数字低通滤波器的幅度特性,wp=0.65*pi/3.14;ws=1.019*pi/3.14;Rp=1;As=15;%设置滤波器参数N,wc=buttord(wp,ws,Rp,As,s)%计算阶数N和3 dB截止频率B,A=butter(N,wc,s)%计算系统函数分子分母多项式系数 S,G=tf2sos(B,A)%将直接型转换为级联型bz,az=bilinear(B,A,1)%用双线性变换Sz,Gz=tf2sos(bz,az)%将直接型转换为级联型,N=6 wc=0.7666B=0 0 0 0 0 0 0.2029A=1.0000 2.9618 4.3862 4.1181 2.5775 1.0228 0.2029S=0 0 1.0000 1.0000 1.4809 0.5876 0 0 1.0000 1.0000 0.3968 0.5876 0 0 1.0000 1.0000 1.0841 0.5876G=0.2029bz=0.0007 0.0044 0.0111 0.0148 0.0111 0.0044 0.0007az=1.0000-3.1824 4.6194-3.7763 1.8118-0.4794 0.0544Sz=1.0000 2.0260 1.0262 1.0000-0.9040 0.2154 1.0000 1.9997 1.0000 1.0000-1.0102 0.3581 1.0000 1.9743 0.9745 1.0000-1.2682 0.7050Gz=7.3937e-004,Bz,Az=bilinear(B,A,Fs)%用双线性变换法将分子分母多项式系数分别为B和A的模拟滤波器系统函数Ha(s)转换为数字滤波器H(z)，Bz和Az分别为H(z)的分子分母多项式系数，Fs为抽样频率。,五.双线性变换法的MATLAB实现,N,wc=buttord(wp,ws,Rp,As),B,A=butter(N,wc,ftype),直接设计数字滤波器的命令：,N,wc=ellipord(wp,ws,Rp,As),B,A=ellip(N,wc),默认用双线性变换法将模拟滤波器转换为数字滤波器。,%用巴特沃斯模拟滤波器作为原型设计数字滤波器,%用椭圆模拟滤波器作为原型设计数字滤波器,用MATLAB实现：,N,wc=buttord(wp,ws,Rp,As,s),B,A=butter(N,wc,ftype,s)