第4章IIR数字滤波器的设计

1、第第4章章 无限单位脉冲响应数字无限单位脉冲响应数字滤波器设计滤波器设计 许多信息处理过程，如：过滤，检测、预测等都要用到滤波器。所谓数字滤波器是指输入、输出均为数字信号，通过一定运算关系改变输入信号的频率及或者相位特性的器件或者软件 。数字滤波器数字滤波器是数字信号处理中使用得最广泛的一种线性系统，是数字信号处理的重要基础。 4.1.1 数字滤波器的分类数字滤波器的分类 数字滤波器主要有以下三种分类方法：（1）根据滤波器的性能来分，可以分为经典数字滤波器和现代数字滤波器。经典数字滤波器，即一般滤波器，特点是有用信号和干扰信号各占不同的频带，通过一个合适的滤波器选择出有用的频率信号滤除干扰信号

2、，因此也称为选频滤波器。现代数字滤波器的特点是针对信号和干扰的频带相互重叠时，按照随机信号内部的一些统计分布规律，从干扰中提取有用信号。例如：维纳滤波器、卡尔曼滤波器、自适应滤波器等最佳滤波器。（2）根据数字滤波器的功能来分，可以分为：低通、高通、带通、带阻、多带数字滤波器。（3）根据实现的网络结构或者从单位脉冲响应来分，可以分成：IIR、FIR滤波器。01( )1MiiiNiiiazH zMNbz，一般0( )hMiiiH zzMiNiiiinybinxany01)()()(1. 数字滤波器的数学描述数字滤波器的数学描述1)差分方程NMzdzcAzbzaZHNiiMiiNiiiMiii一般1

3、11110)1 ()1 (1)(2)系统函数4.1.2 数字滤波器的技术指标 式中： 称为幅频特性， 称为相频特性。下面以低通IIR数字滤波器为例，介绍数字滤波器的技术指标。()() |()|jjjH eH ee 假设数字滤波器的传输函数用下式表示：|()|jH e( ) cjjeHeH)1lg(2011lg20| )(| )(|lg2011minmaxcjeH, 1| )(|11rjeH,| )(|2通带内波动阻带内最小衰减rjeHAt22maxlg201lg20| )(|1lg20问：其他类型滤波器的频率特性及技术指标？4.1.3 数字滤波器设计方法概述数字滤波器设计方法概述 i1）按照实

4、际需要确定滤波器的类型和性能指标要求；2）用一个因果稳定系统的 H(z) 或 h(n) 去逼近这个性能要求；3）用一个有限精度的运算去实现这个系统函数。包括选择运算结构：如级联型、并联型、卷积型、频率采样型以及快速卷积（FFT）型等。数字滤波器的设计思想都是采用逼近方法。设计的基本步骤是： IIR数字滤波器设计方法可以分为两类：直接设计法和间接设计法。 （1）直接设计法直接设计法是直接在频域或者时域进行逼近设计，有零极点累试法、频域逼近法、时域逼近法，通常必须借助计算机进行优化设计，因此也称为最优化设计方法。具体步骤分两步： 1) 确定一种最优准则，如：最小均方误差准则，即使设计出的实际频率响

5、应的幅度特性与所要求的理想频率响应的均方误差最小。min)()(21MijdjiieHeH2) 在此最佳准则下，通过不断地迭代运算求滤波器的系数，直到满足要求为止。此外还有其他多种误差最小准则。由于计算机技术的发展，最优化设计方法的使用也逐渐增多；但是，由于涉及到最优化设计算法，而且目前已经有商业设计程序可以使用，因此本章不介绍。（2）间接设计法间接设计法是借助于模拟滤波器的设计方法进行的，其设计步骤是：先设计模拟滤波器得到传输函数，然后将按某种方法转换成数字滤波器的系统函数。由于模拟滤波器的设计方法已经发展得很成熟，而且模拟滤波器有简单而严格的设计公式，设计起来方便、准确；因此可将这些理论推

6、广应用到数字域，作为设计数字滤波器的工具。本章后续内容将介绍这种设计方法。 模拟滤波器根据幅度特性可以分为低通、高通、带通和带阻滤波器，它们的理想特性如图所示。 4.2 4.2 模拟滤波器模拟滤波器的的设计设计模拟滤波器的设计就是根据一组设计规范设计模拟系统函数Ha(s)，使其逼近某个理想滤波器特性。通常都是先设计低通滤波器，再通过频率变换将低通变换为希望类型的滤波器；因此，下面先介绍模拟低通滤波器的技术指标和逼近方法，然后再介绍模拟滤波器的频率变换。4.2.1 模拟低通滤波器的技术指标及逼近方模拟低通滤波器的技术指标及逼近方法法模拟低通滤波器的设计指标用4个参数描述。 为通带截止频率 ； 为

7、幅度下降到0.707，也就是下降了3dB时的频率，因此通常被称为3dB截止频率； 为阻带截止频率。 pcr22maxmaxmin2minmin|()|()|20lg10lg10lg|()|()|()|jjjjjH eH eH eH eH e2maxmax20lg |()|10lg |()|jjAtH eH e 0)()(dtethjHtjaa0sincos)()(dttjtthjHaa)()(jHjHaa对于一般滤波器的单位脉冲响应 为实数的因果系统有： ( )ah t定义振幅平方函数222()()()()()()()( )()aaaaaaasjAHjHjHjAHjHjHs Hs 又 S=j，

8、2= - S2 A(2)=A(-S2)|S=j问题：如何由A(-S2)| Ha(S)？对于给定的A(-S2)，先在S复平面上标出A(-S2)的极点和零点，由于A(-S2)的极点和零点总是“成对出现”，且对称于S平面的实轴和虚轴；故选用选用A(-S2)在在S左半平面的极点作为左半平面的极点作为Ha(s)的极点，零的极点，零点可选用任一半点可选用任一半。 可见只要给定幅度平方函数模型，模拟滤波器就可以设计出来。 4.2.2巴特沃兹（巴特沃兹（Butterworth）滤波器逼近）滤波器逼近NcajjjHA22211)()(特点：具有通带内最大平坦的振幅特性，且随f ，幅频特性 单调 。巴特沃兹滤波器

9、 (Butterworth 滤波器) 幅度平方函数：巴特沃兹滤波器幅度平方函数图1中，N增加，通带和阻带的近似性越好，过渡带越陡。通带内，分母/c1， ( /c)2N 1, 增加， A(2)快速减小。=c, 幅度衰减相当于3db衰减点。21)0()(2AAc21)(2A可以证明：理想幅度响应及其逼近响应的差值，与=0处的导数等于0的阶数有关，应尽量使导数等于0的阶次高些。N增加，通带内频率响应变得更为平坦，阻带内的衰减更大；过渡带内频响更趋于斜率为过渡带内频响更趋于斜率为-6NdB/倍频程的渐倍频程的渐近线近线。 振幅平方函数的极点：NcaajSSHSH2)(11)()(可见，Butter w

10、orth滤波器 的振幅平方函数有2N个个极点极点，它们均匀对称地分布在|S|=c的圆周上。 例：为N=3阶BF振幅平方函数的极点分布，如图。1121()222( 1)()kjNNPkccSje 3阶A(-S2)的极点分布 考虑到系统的稳定性，可知DF的系统函数是由S平面左半部分的极点（平面左半部分的极点（SP0，SP1，SP2）组成的，它们分别为：2233012,jjpcpcpcSeSSe30120121( )()()()()()()caPPPPPPccccccHsSSSssssSsSsSccsjcpj101( )()aNpkkHppS121()22kjNpkSe 称为归一化的极点，它只与N有

11、关；因此，当N确定时，归一化的极点可以事先算出并建立一个参数表。 ( )aHp2101211( ).aNNNHpbb pb pbpp可以写成以下分母多项式形式 当N确定时，通过查表就可以得到滤波器的归一化原型，例如：当N=3时，查表得： 321( )221aHppppcsp 321( )( /)2( /)2( /) 1acccHssss设计的巴特沃兹模拟低通滤波器传递函数为 ：问题：如何确定滤波器的问题：如何确定滤波器的N、 c ？ 根据滤波器的技术指标来确定的。NcpjjjeHp22)(11lg10|)(|lg102maxmax| )(|lg10| )(|lg20jjeHeHAt2min2m

12、in2maxminmax| )(|lg10| )(| )(|lg10| )(| )(|lg20jjjjjeHeHeHeHeHNcrjj2)(11lg10查表得H(p)，则：Ncp21 . 0)(110lglgspspkNNrpAt21 . 01 . 0)(110110NcrAt21 .0)(1100.120.1101101spAtkcspapHsH/| )()(pprrfspfNpc211.0)110(NAtrc211.0)110(【例4.2.1】 已知模拟低通滤波器通带截止频率fp=0.5Hz，通带最大衰减=2dB，阻带截止频率fr=1.2Hz，阻带最小衰减=30dB，试采用巴特沃斯模型设计

13、该模拟低通滤波器。 解：利用上述的方法设计，设计步骤为： (1) 确定模拟低通滤波器的阶数N。 0.10.11010.0242101spAtk0.4167pprrfspflg4.25lgspspkN（2）查表得到5阶归一化低通原型滤波器：b0=1.0000，b1=3.2361，b2=5.2361，b3=5.2361，b4=3.2361取N=5。 （3）求3dB截止频率 23451( )1 3.23615.23615.23613.2361aHpppppp10.12(101)20.52755/Ncpards 10.12(101)20.6015/AtNcrards （4）去归一化，得到所设计的低通滤

14、波器：554233245( )3.23615.23615.23613.2361cacccccHssssss MATLAB提供三个函数用于ButterWorth模拟滤波器的设计。 （1）阶数及截止频率的求取函数buttordbuttord函数的语句格式为：n,Wn = buttord( Wp, Wr, Rp, Rr, s) 其中：Wp为通带截止频率、Wr为阻带截止频率，单位是rad/s；Rp为通带最大衰减、Rr为阻带最小衰减，单位是dB；选项s表示设计模拟滤波器，省略此参数为设计数字滤波器。 当设计带通、带阻滤波器时，Wp，Wr为2元数组。函数的返回值Wn为截止频率，n为满足技术指标的滤波器的最

15、小阶数。（2）设计函数butter有两种语句格式，当设计低通或者带通（带通的 Wn为2元数组）滤波器时，语句格式为： b,a = butter(n, Wn, s) 当设计其他类型滤波器时，采用ftype说明滤波器的类型，语句格式为： b,a = butter(n, Wn, ftype, s) 式中：ftype可以为high,low , stop，分别表示设计的滤波器为高通、低通、带阻滤波器。函数的返回参数b,a为模拟滤波器的传递函数模型：11(0)(1)( )( )(0)(1)( )nnannbsbsb nHsasasa n （3）归一化模拟低通原型设计函数 z,p,k=buttap(n)，n

16、为阶数，z为零点，p为极点，k为放大系数。【例4.2.2】例4.2.1的MATLAB设计代码为： 方法（1）：利用模拟低通原型设计，程序如下： fp=0.5*2*pi;ap=2;fr=1.2*2*pi;as=30; n,Wn = buttord(fp,fr,ap,as,s) z,p,k = buttap(n);%低通原型设计 b0=k*real(poly(z); %poly函数是把多项式根转换为多项式 系数，real是复数取实部运算 a0=real(poly(p); b,a=lp2lp(b0,a0,Wn) ； h,f=freqs(b,a); mag = abs(h);mag = 20*log1

17、0(mag); phase = angle(h); subplot(2,1,1);plot(f/(2*pi),mag); title(N=5 Butterworth Lowpass Filter); axis(0 1.5 -35 0);xlabel(f(Hz);ylabel(幅度(dB); grid; subplot(2,1,2);plot(f/(2*pi),phase); axis(0 1.5 -4 4);xlabel(f(Hz);ylabel(相位(ard);方法（2）：直接设计模拟滤波器，修改程序如下： fp=0.5*2*pi;ap=2;fr=1.2*2*pi;as=30; n,Wn =

18、 buttord(fp,fr,ap,as,s) b,a = butter(n,Wn,s) h,f=freqs(b,a); 程序运行结果： （1）滤波器的阶数N = 5， 与上述计算结果相同。 （2）截止频率Wn = 3.7792rad/sec= 0.6015Hz，与计算结果相同。 （3）程序运行可得： b = 770.9440 a = 1.0000 12.2299 74.7850 282.6305 660.1398 770.9440 则所设计的滤波器的传递函数为：5432770.9440( )s12.2299s74.7850s282.6305s66 0.1398s770.9440aHs （4）

19、滤波器的特性分析结果可以看出在通带截止频率0.5Hz处的衰减小于2dB，阻带截止频率1.2Hz处阻带最小衰减大于30dB，满足指标要求。相位特性曲线在通带内相位特性接近线性。切比雪夫滤波器具有波动性。它有两种类型：（1）切比雪夫型，在通带中具有等波动响应；（2）切比雪夫型，在阻带中具有等波动响应。切比雪夫滤波器在通带范围内是等幅起伏的，所以同样的通带衰减，其阶数较巴特沃兹同样的通带衰减，其阶数较巴特沃兹滤波器滤波器要小要小。可根据需要对通带内允许的衰减量（波动范围）提出要求，如：要求波动范围小于1dB等。4.2.3切比雪夫滤波器切比雪夫滤波器切比雪夫型滤波器的振幅平方函数为 )(11)()(2

20、222cNaVjHA1)coshcosh(1)coscos()(11xxNxxNxVN)(, 11)(,1xVxxxVxNN时VN（x）为N阶切比雪夫多项式，定义为式中：表示|Ha(j)|波动范围的参数； 当 =0时，)2(cos11)0arccos(cos11)(22220NNjHa22011)(jHa1)(02jHa1)2(cos2N0)2(cos2N2111c ，随/c ， 0 (迅速趋于零)N为偶数，N为奇数，通带内变化范围切比雪夫滤波器的振幅平方特性 切比雪夫型的平方幅度响应为：切比雪夫型滤波器通带幅频特性是单调的，而阻带是等波动的。这种滤波器在s平面上既有极点，又有零点。切比雪夫模

21、拟滤波器的设计计算比较复杂，因此只介绍用Matlab设计。MATLAB提供多个函数用于切比雪夫模拟滤波器的设计，下面以型为例进行介绍；对于切比雪夫型滤波器设计，只要把切比雪夫型设计函数中的“1”改为“2”，即为切比雪夫型滤波器设计函数。21221()()1aNcJHjT （1）阶数求取函数cheb1ord，语句格式为：n,Wn = cheb1ord( Wp, Wr, Rp, Rr, s) 其中：Wp、Wr为角频率，单位是rad/s。当所设计的滤波器为带通、带阻滤波器时，Wp，Wr为2元数组；返回值Wn为截止频率，n为满足指标滤波器的最小阶数。（2）设计函数cheby1，语句格式有两种：设计低通

22、或者带通（带通的 Wn为2元数组）滤波器时格式：b,a = cheby1(n,Rp,Wn,s) 设计的是其他类型的滤波器时格式为：b,a = cheby1(n, Rp, Wn, ftype, s)式中：ftype可以为high,low , stop，分别表示设计的滤波器为高通、低通、带阻滤波器。（3）归一化模拟低通原型设计函数 z,p,k=cheb1ap(n,Rp)，n为阶数，Rp为通带最大衰减，z为零点，p为极点，k为放大系数。【例4.2.3】例4.2.1的切比雪夫型设计。方法（1）：先设计低通原型，再转换成所要设计的模拟低通滤波器。程序代码： fp=2*pi*0.5;Ap=2;fr=2*p

23、i*1.2;Ar=30; N,Wn = cheb1ord(fp,fr,Ap,Ar,s) z,p,k = cheb1ap(N,Ap);%低通原型 A,B,C,D = zp2ss(z,p,k);%转换为状态空间表达式 At,Bt,Ct,Dt=lp2lp(A,B,C,D,Wn);%低通原型转换 b,a=ss2tf(At,Bt,Ct,Dt); %转换为传输函数表达式 h,f=freqs(b,a); mag = abs(h);mag = 20*log10(mag); phase = angle(h); subplot(2,1,1); plot(f/(2*pi),mag); title(Cheby1 Lo

24、wpass Filter); axis(0 1.5 -35 0); xlabel(f(Hz);ylabel(幅度(dB); grid; subplot(2,1,2); plot(f/(2*pi),phase); axis(0 1.5 -4 4); xlabel(f(Hz);ylabel(相位(ard);程序运行结果：（1）滤波器的阶数N = 3（小于巴特沃斯的阶数）。（2）截止频率fn = 3.1416 ard/s=0.5 Hz 。（3）传递函数：b = 0 0 0 10.1356 a = 1 2.3179 10.0886 10.1356 得到所设计的滤波器的传递函数为：3210.1356(

25、)s2.3179s10.0886s10.1356aHs （4）幅频特性曲线：通带内等波纹波动并且在截止频率0.5Hz处的衰减小于2dB，阻带截止频率1.2Hz处阻带最小衰减大于30dB，满足指标要求；相位特性曲线：在通带内相位特性接近线性。 方法（2）：直接采用Matlab函数设计，代码设计部分改为： fp=0.5*2*pi;ap=2;fr=1.2*2*pi;as=30 n,Wn = cheb1ord(fp,fr,ap,as,s) b,a = cheby1(n,fp,Wn,s) h,f=freqs(b,a); 请读者在上述程序，加入绘制幅频特性、相频特性曲线的代码，调试、运行程序，并与前面结果

26、比较，看是否相同。 椭圆滤波器也称考尔滤波器（Cauer Filter），其振幅平方函数为：),(11)()(2222LRjHANa式中：RN（，L）为N阶雅可比椭圆函数 L是表示波纹性质的参量4.2.4椭圆滤波器（Elliptic Filter） 可见，在归一化通带内（-11）， 在（0，1）间振荡，而超过L后， 在 间振荡。 。 这一特点使滤波器同时在通使滤波器同时在通带和阻带具有任意衰减量带和阻带具有任意衰减量。 ),(25LRLL,2L),(25LR),(25LR的特性曲线 椭圆滤波器的振幅平方函数 图中和A的定义同切比雪夫滤波器特点：幅值响应在通带和阻带内都是等波纹的幅值响应在通带和

27、阻带内都是等波纹的，对于给定的阶数和波纹要求，椭圆滤波器能获得较其它滤波器更窄的过渡带宽更窄的过渡带宽；就这点而言，椭圆滤波器是最优的。但是其设计计算最复杂，因此通常采用MATLAB设计，设计设计函数函数有：有： n, Wn = ellipord( Wp, Wr, Rp, Rr, s) z, p, k = ellipap( n, Rp, Rr) b, a = ellip( n, Rp, Rr, Wn, s) b, a = ellip( n, Rp, Rr, Wn, ftype ,s) 实际应用中：在给定相同指标下，三种模型设计都能满足指标要求，但椭圆滤波器阶次最低，切比雪夫次之，巴特沃兹最高，

28、参数的灵敏度则恰恰相反。设计时需根据要求，合理选用。4.2.5 模拟高通、带通、带阻滤波器设计模拟高通、带通、带阻滤波器设计模拟高通、带通、带阻滤波器设计是通过设计模拟低通滤波器再通过频率变换实现的。假设： 模拟低通滤波器的传输函数用 表示， 归一化频率用 表示， 令： ，称为归一化拉斯复变量， 为归一化低通滤波器； 所要设计类型滤波器的传输函数用 表示， 归一化频率用 表示， 令： ，称为归一化拉斯复变量， 称为归一化传输函数。 ( )G spj( )G p( )H sqj( )H q4.2.5.1 低通到高通的频率变换 低通到高通的频率变换公式： 11H(j )()G j4.2.5.2.

29、低通到带通的频率变换低通到带通的频率变换带通滤波器归一化边界频率用以下式子计算： 11/rrB22/rrB/llB /uuB 22200/luB ， 带通滤波器频率转换为归一化低通滤波器频率的关系式： 220220B 或 归一化低通传输函数到带通传输函数的转换公式 220qpq/qs B2220()ululsspssB 2()H(s)( )ululspsG p 4.2.5.3 低通到带阻的频率变换低通到带阻的频率变换220220B 或 2220()ululssBpss 2()H(s)( )ululspsG p 4.2.5.4 由模拟低通原型设计模拟高通、带通、带阻由模拟低通原型设计模拟高通、带

30、通、带阻滤波器的滤波器的Matlab函数函数Matlab提供的函数如下：1.低通原型转换为低通：NUMT,DENT = LP2LP(NUM,DEN,Wo)，Wo为截止频率(rad/sec)；2.低通原型转换为高通：NUMT,DENT = LP2HP(NUM,DEN,Wo)，Wo为截止频率(rad/sec)；3.低通原型转换为带通：NUMT,DENT = LP2BP(NUM,DEN,Wo,Bw)，Wo为中心频率(rad/sec)，Bw为通带带宽(rad/sec)；4.低通原型转换为带阻：NUMT,DENT = LP2BS(NUM,DEN,Wo,Bw) ，Wo为中心频率(rad/sec)，Bw为阻

31、带带宽(rad/sec)；4.2.5.5 模拟高通、带通、带阻滤波器直接设计的模拟高通、带通、带阻滤波器直接设计的Matlab函数函数Matlab提供了直接设计模拟滤波器函数：（1）b,a = butter(N, Wn, ftype, s)ftype选项确定滤波器类型；可选high、low、stop缺省时，设计的是低通、带通滤波器。当设计带通、带阻滤波器时，Wn为2元数组，结果为2*N阶。（2）切比雪夫型滤波器：b,a = cheby1(N, Rp,Wn, ftype, s)，Rp为通带最大波动（dB）。（3）切比雪夫型滤波器：b,a = cheby2(N, Rr,Wn, ftype, s)，

32、Rr为阻带最小波动（dB）。（4）椭圆滤波器：b,a = ellip(N, Rp，Rr，Wn, ftype, s)，Rp为通带最大波动（dB），Rr为阻带最小波动（dB）。例4.2.4 3阶巴特沃斯型的归一化低通滤波器的传递函数为试设计下列模拟滤波器：（1）通带为10Hz的低通滤波器；（2）通带下边频为10Hz的高通滤波器；（3）中心频率为10Hz，带宽为2Hz的带通滤波器；（4）中心频率为10Hz，带宽为2Hz的带阻滤波器；又设采样周期为0.01秒，求相应的数字滤波器。解：Matlab程序如下：b=1;a=1 2 2 1;Wn=10*2*pi; B=2*2*pi; 3211( )221H s

33、sss b1,a1=lp2lp(b,a,Wn)；b2,a2=lp2hp(b,a,Wn) b3,a3=lp2bp(b,a,Wn,B)；b4,a4=lp2bs(b,a,Wn,B) 运行结果： b1 = 2.4805e+005 a1 = 1 125.66 7895.7 2.4805e+005 b2 = 1 1.2095e-014 -4.2767e-013 -2.7202e-020 a2 = 1 125.66 7895.7 2.4805e+005 b3 = 1984.4 -5.8081e-012 -7.7026e-009 -2.9869e-013 a3 = 1 25.133 12159 2.0042e

34、+005 4.8003e+007 3.9171e+008 6.1529e+010 b4 = 1 3.5545e-013 11844 2.8096e-009 4.6756e+007 5.5449e-006 6.1529e+010 a4 = 1 25.133 12159 2.0042e+005 4.8003e+007 3.9171e+008 6.1529e+01053252.4805 101( ) s 125.66s7895.7s 2.4805 10H s33252( ) s 125.66s7895.7s2.4805 10sHs 36545 37 28101984.43( ) s 25.133s1

35、2159s2.0042 104.8003 103.9171 106.1529 10sH ssss647 2106545 37 2810118444.676 106.1529 104( ) s25.133s12159s2.0042 104.8003 103.9171 106.1529 10sssH ssss4.3 利用模拟滤波器设计利用模拟滤波器设计IIR数字滤波器数字滤波器 利用模拟滤波器设计IIR数字滤波器， Ha(s) H（z） S平面 Z平面 这种映射变换应遵循两个基本原则：1）H（z）的频响要能模仿Ha(s)的频响，即S平面的虚轴应映射到Z平面的单位圆上。2）Ha(S) 的因果稳定性映

36、射成 H（z）后保持不变，即S平面的左半平面 ReS0 应映射到Z平面的单位圆以内|Z|1。脉冲响应不变法是从滤波器的脉冲响应出发，使:数字滤波器的单位脉冲响应序列数字滤波器的单位脉冲响应序列h(n)h(n)正好等于模拟正好等于模拟滤波器的冲激响应滤波器的冲激响应h ha a(t)(t)的采样值的采样值，即: h(n)=ha(nT)T为采样周期 假设：Ha(s) 及 H（z）分别表示 ha(t) 的拉氏变换及 h(n) 的 Z 变换，即: Ha(s)=Lha(t) H(z)=Zh(n)计算 H(z) ： 4.3.1脉冲响应不变法脉冲响应不变法4.3.1.1变换关系变换关系假设模拟滤波器的传递函

37、数若只有单阶极点，且分母的只有单阶极点，且分母的阶数高于分子阶数阶数高于分子阶数 NM，则可表达为部分分式形式；NiiissAsHa1)(NitsiatutueAthi1)(),()(NiNinTsinTsinueAnueAnThanhii11)()()()()(其拉氏反变换为： 对ha(t)采样得到数字滤波器的单位脉冲响应序列 再对h(n)取Z变换，得到数字滤波器的传递函数：01011)()(nNinnTsNiinnTsizeAzeAzHiikTskTszezeii111)(1, 0)(1kkTszeiNiTsizeAzHi111)(NiiissAsHa1)(比较看到： （1）S平面上的极点

38、 S=Si 变换到Z平面上是极点Tsiez 1)(TsRTsieiee（2）Ha(s)与H(Z)中部分分式所对应的系数不变（3）这种Ha(s)到H(Z)的对应变换关系，只有将Ha(s)表达表达为部分分式形式才成立为部分分式形式才成立。4.3.1.2稳定性分析稳定性分析 如果模拟滤波器是稳定的，则所有极点 Si 都在S左半平面，即 Resi0 , 那么变换后H(z)的极点也都在单位圆以内，即： 因此，数字滤波器保持稳定。4.3.1.3频率映射关系频率映射关系根据理想采样序列拉氏变换与模拟信号拉氏变换的关系。maamTjsHTsH21)( nnsTnnsTnststnenhenThdtenTtth

39、dtenTtthsH)()()()()()()(aaaasTez s平面与z平面的映射关系 以上表明，采用脉冲响应不变法将模拟滤波器变换为数字滤波器时，它所完成的 S 平面到 Z 平面的变换，正是以前所讨论的拉氏变换到Z变换的标准变换关系，即：（1）首先对）首先对Ha(s)作周期延拓，作周期延拓，（2）然后再经过映射关系映射到）然后再经过映射关系映射到 Z 平面上。平面上。sTzejT2jsrezMj,)2(令TMerT2,:则TMT)2(S平面上每一条横带，都将重叠地映射到Z平面的整个平面上: 每一横带的左半部分映射到每一横带的左半部分映射到Z平面单位圆以内，平面单位圆以内， 每一横带的右半

40、部分映射到每一横带的右半部分映射到Z平面单位圆以外，平面单位圆以外，TjSTeez)(轴上每一段轴上每一段都对应于绕单位圆一周。都对应于绕单位圆一周。j0TT3T3T)Im( zj)Re( z0S 平面Z 平面: 4.3.1.4频率响应分析频率响应分析数字滤波器的频响并不是简单的重现模拟滤波器的频响，而是模拟滤波器频响的周期延拓：TmwjHTeHamj21Ts2 TmjsHTsHzHamaezsT21)(如果模拟滤波器的频响带限于折叠频率折叠频率S/2 以内，20)(sajH 这时数字滤波器的频响才能不失真地重现模拟滤波器的频响（存在于折叠频率S/2以内）)(1)(TjHTeHaj但任何一个实

41、际的模拟滤波器，其频响都不可能是真正带限的，因此不可避免地存在频谱的交叠，即混淆。)(jeH)(TjHa00脉冲响应不变法中的频响混淆【例4.3.1】将一个具有如下传递函数的模拟滤波器变换为数字滤波器。 3111) 3)(1(2)(sssssHTTezezzH3111111)(243131)(1)(zeeezeezTTTTT模拟滤波器的频率响应为: 4)3(2) 3)(1(2)()(2jjjsHjHajs数字滤波器的频率响应为: 2433)(1)()()(jTjTTjTTezjeeeeeeeezHeHj 4.3.1.5脉冲响应不变法的特点脉冲响应不变法的特点（1）频率变换是线性的，即：与是线性

42、关系。因此如果模拟滤波器的频响带限于折叠频率以内的话，通过变换后数字滤波器的频响可不失真地反映原响应与频率的关系。（2）如果模拟滤波器是稳定的，映射后得到的数字滤波器也是稳定的。（3）缺点：有频谱周期延拓效应，会产生频率混叠。因此只能用于带限的频响特性，如衰减特性很好的低通或带通滤波器。至于高通和带阻滤波器，不适用。因此，脉冲响应不变法只有在一定要满足频率线性关系或严格要求保持瞬态响应时才采用。 脉冲响应不变法的主要缺点是频谱交叠产生的混淆，这是从S平面到Z平面的标准变换z=esT的多值对应关系导致的，为了克服这一缺点，设想变换分为两步变换分为两步： 第一步：将整个S平面压缩到S1平面的一条横

43、带里； 第二步：通过标准变换关系将此横带变换到整个Z平面上去。 由此建立S平面与Z平面一一对应的单值关系，消除多值性，也就消除了混淆现象。4.3.2 双线性变换法双线性变换法4.3.2.1变换关系变换关系s平面s1平面z平面双线性变换法的映射关系 为了将S平面的j轴压缩到S1平面j1轴上的/T到/T 一段上，可通过以下的正切变换正切变换实现：0)2(tan1Tc取C=2/T，经过这样的频率变换， 当由时, 1由-/T经过变化到/T ，即S平面的整个j轴被压缩到S1平面的。将这一关系解析扩展至整个S平面，则得到S平面到S1平面的映射关系：TsTseecTscs1111)2(tanh1zzTsTs

44、ez1sTsTz)/()/(表明：（1）从s平面到z平面之间的映射关系都为一一对应的线性变换关系，因此该变换称为双线性变换。（2）当s左半平面映射在单位圆内，右半平面映射到单位圆外；因此，稳定的模拟滤波器通过双线性变换后，所得到的数字滤波器也是稳定的。 S平面的虚轴(整个j)对应于Z平面 单位圆的一周，S平面的=0处对应于Z平面的=0处，对应的数字滤波器的频率响应终止于折迭频率处，所以双线性变换不不存在混迭效应存在混迭效应。)2/tan(2T考虑 z = ej , jjTjTeeTsjj)2(tan22cos)2/sin(21124.3.2.2特点特点 （1）s平面与z平面是单值的一一对应关系

45、，即整个j轴单值的对应于单位圆一周，不存在频率混叠。（2） 模拟域频率（）与数字域频率（）之间的变换成非线性关系：2tan2T这将导致如下问题： a. 数字滤波器的幅频响应相对于模拟滤波器的幅频响应有畸变。例如：一个模拟微分器，它的幅度与频率是直线关系，但通过双线性变换后，就不可能得到数字微分器bkjHeHbkjHj2tan)()()(2tanb. 线性相位模拟滤波器经双线性变换后，得到的数字滤波器为非线性相位非线性相位。 c.要求模拟滤波器的幅频响应必须是分段恒定分段恒定的，故双线性变换只能用于设计低通、高通、带通、带阻等选频滤波器。d. 设计时需要预畸变即将模拟滤波器的临界频率事先加以畸变

46、事先加以畸变，然后通过双线性变换后正好映射到所需要的频率上。 2tan2iiT双线性变换时频率的预畸变（3）双线性变换比脉冲响应法的设计计算更直接和简单 从模拟滤波器的传递函数可直接通过代数置换得到数字滤波器的系统函数。 变换关系式为11112112 1( )( )1aazsTzzH zHsHTz2tan22()()tan2jaaTH eHjHjT 【例4.3.2】设有一模拟滤波器21( )1aHsss试用双线性变换法将它转变为数字系统函数H(z) ，设采样周期T = 2。解： T = 2时变换公式变为： 1111zsz1121111111()( ) |11111azszHzHszzzz121

47、222122(1)12330 .3 3 3 30 .6 6 6 70 .3 3 3 310 .3 3 3 3zzzzzzzz4.3.3 用用Matlab实现模拟滤波器变换到实现模拟滤波器变换到IIR数字滤波器数字滤波器（1）脉冲响应不变法的Matlab变换函数语句格式为： BZ,AZ = impinvar(B,A,Fs)（2）双线性变换法的Matlab变换函数语句格式为： BZ,AZ = bilinear(B,A,Fs)其中： B为模拟滤波器系统函数的分子多项式系数矢量， A为模拟滤波器系统函数的分母多项式系数矢量， BZ为数字滤波器系统函数的分子多项式系数矢量， AZ为数字滤波器系统函数的分

48、母多项式系数矢量， Fs为采样频率（单位为Hz）。【例4.3.3】用Matlab实现例4.3.1具有如下传递函数的模拟滤波器变换为数字滤波器，设采样频率1Hz。 解：Matlab代码为： b=2;a=1 4 3;Fs=1; bz,az=impinvar(b,a,Fs); 程序运行结果： bz = 0 0.3181 az = 1.0000 -0.4177 0.0183 程序运行结果与例4.3.1的计算结果一致。222( )(1)(3)43H sssss 【例4.3.4】用Matlab实现例4.3.2。 解：Matlab代码： b=1;a=1 1 1;Fs=0.5; bz,az=bilinear(

49、b,a,Fs); 程序运行结果： bz = 0.3333 0.6667 0.3333 az = 1.0000 0.0000 0.3333 与例4.3.2的计算结果一致。4.4 从模拟低通原型到各种数字滤波器从模拟低通原型到各种数字滤波器的频率变换的频率变换 原型变换 映射变换 原型变换模拟原型模拟低通、高通带通、带阻数字低通、高通带通、带阻也可把前两步合并成一步，直接从模拟低通归一化原型通过一定的频率变换关系，完成各类数字滤波器的设计。下面举例讨论设计各种数字滤波器的变换方法 通过模拟低通原型设计数字低通滤波器的步骤： 1）确定数字滤波器的性能指标临界频率k、通带阻带衰减。2）由变换关系将k映

50、射到模拟域，得出模拟低通滤波器的临界频率值k。3）根据k设计模拟滤波器的Ha(s)4） 把Ha(s) 变换成 H(z)（数字滤波器传递函数）变换方法有两种：脉冲响应不变法和双线性变换法，下面举例说明。 4.4.1 低通变换低通变换【例4.4.1】设采样周期 ，已知数字低通滤波器通带截止频率 ，通带最大衰减2dB，阻带截止频率 ，阻带最小衰减20dB，试采用巴特沃斯模型设计该数字低通滤波器。250(4)sTs fkhz解法一：脉冲响应不变法1）确定数字低通滤波器的技术指标。0.5pfkHz1.2rfkHz ，通带最大衰减 =2dB； ，阻带最小衰减 =20dB。2）由于脉冲响应不变法的频率变换关

51、系是，因此有，所以相应的模拟低通滤波器的指标为： ，通带最大衰减 =2dB； ，阻带最小衰减 =20dB。20.25ppf T20.6rrf TAtAt2/1ppf T Tk 22.4rrf Tk 因此，取N=3。查表得到模拟低通原型滤波器：3）设计巴特沃兹模拟低通滤波器0.10.11010.0769101spAtk0.4167pprrfspflg2.93lgspspkN231( )122Ha ssss10.12(101)20.5468/Ncpkards 10.12(101)20.5579/AtNcrkard s 去归一化得模拟低通滤波器： 231( )1 2( /) 2( /)( /)acc

52、cH ssss 4）脉冲响应不变法变换可见，H（z）与采样周期T有关，T越小，H（z）的相对增益越大，这是不希望的。改进：H（z）乘以因子T，使H（z）只与 有关。/6/6/3/3( )(13)/ 2(13)/ 2jjcceceHa sscscjscj11( )1iNisTiAH zez/6/61(13)/2(13)/211/(/3 )(/3 )( )111cccjjcccjjTT eT eH zezezez111210.87650.87650.6349( )1 0.41621 0.93610.4162zH zTzzz/csff 解法二：双线性变换法 1）首先确定滤波器数字域指标，结果与脉冲响

53、应不变法的步骤1）相同。 2）确定模拟低通滤波器的技术指标。根据频率的非线性关系，进行频率预畸变。11120.87650.87650.6349( )1 0.41621 0.93610.4162zH zzzz22tan0.41422ppTT 22tan1.37642rrTT 3) 设计模拟低通滤波器因此，取N=3。3dB截止频率：0.10.11010.0769101spAtk0.3009pprrfspflg2.1360lgspspkN10.122(101)0.6399/AtNcrkard sT 231( )12( /)2( /)( /)acccHssss4） 进行双线性变换 11112 11/1

54、.562711cczzsTzz 231111111( )111121.562721.56271.5627111H zzzzzzz31123112.8257 10.20736.43911.0575zzzz 1、从模拟低通原型到数字低通滤波器的程序： fs=4000;Ts=1/fs;fp=500;ap=2;fr=1200;ar=20; %以下用脉冲响应不变法设计 wp=2*pi*fp/fs;wr=2*pi*fr/fs;% wp、wr为数字频率 Wp=wp/Ts;Wr=wr/Ts; % Wp、Wr为模拟频率 N1,Wn = buttord(Wp,Wr,ap,ar,s) z,p,k=buttap(N1

55、); %模拟低通原型设计 b0=k*real(poly(z); a0=real(poly(p); B,A=lp2lp(b0,a0,Wn); %模拟低通原型变换到模拟低通滤波器 num1,den1=impinvar(B,A,fs); %模拟滤波器经脉冲响应不变法变换到数字滤波器 h1,w=freqz(num1,den1);解法三：用MATLAB完成设计%以下用双线性变换法设计Wp=2/Ts*tan(wp/2);Wr=2/Ts*tan(wr/2); %频率预畸变N2,Wn = buttord(Wp,Wr,ap,ar,s)；z,p,k=buttap(N2); b0=k*real(poly(z); a

56、0=real(poly(p); B,A=lp2lp(b0,a0,Wn) num2,den2=bilinear(B,A,fs); %双线性变换 h2,w=freqz(num2,den2); f=w/pi*fs/2; figure; plot(f,abs(h1),-.,f,abs(h2),-); text(600,0.66,leftarrow脉冲响应不变法); text(880,0.46,leftarrow双线性变换法); grid; xlabel(频率/Hz )； ylabel(幅度) 2、从模拟滤波器到数字滤波器的设计程序： fs=4000;Ts=1/fs; fp=500;ap=2; fr=1

57、200;ar=20; %以下用脉冲响应不变法设计 wp=2*pi*fp/fs;wr=2*pi*fr/fs;%wp、wr为数字频率 Wp=wp/Ts;Wr=wr/Ts; %Wp、Wr为模拟频率 N1,fn = buttord(Wp,Wr,ap,ar,s) B,A=butter(N1,fn,s)； num1,den1=impinvar(B,A,fs); h1,w=freqz(num1,den1);%以下用双线性变换法设计Wp=2/Ts*tan(wp/2);Wr=2/Ts*tan(wr/2); N2,fn = buttord(Wp,Wr,ap,ar,s)B,A=butter(N2,fn,s)num2

58、,den2=bilinear(B,A,fs); h2,w=freqz(num2,den2); f=w/pi*fs/2; figure; plot(f,abs(h1),-.,f,abs(h2),-); text(600,0.66,leftarrow脉冲响应不变法); text(880,0.46,leftarrow双线性变换法); grid; xlabel(频率/Hz )； ylabel(幅度) 设计结果满足指标要求。 设计高通数字滤波器时，只能采用双线性变换法。在模拟高通滤波器的设计中，低通至高通的变换是s变量的倒置，这一关系应用于双线性变换，则从模拟低通到数字高通的变换，只要将双线性变换式中的

59、s代之以1/s，就可得到模拟低通到数字高通滤波器的变换关系式即：11112zzTs4.4.2 高通变换高通变换 由于倒数关系不改变模拟滤波器的稳定性，因此，也不会影响双线变换后的稳定条件，而且 轴仍映射在单位圆上，只是方向颠倒了。即jjjTeeTseZjjj2cot2112,时2cot2T如图22Tcot 01z1z01.01.00高通变换频率关系这一曲线的形状与双线性变换时的频率非线性关系曲线相对应，只是将 坐标倒置，因而通过这一变换后可直接将模拟低通变为数字高通，如图。模拟低通原型到数字高通变换 应当明确： 所谓高通DF，并不是高到2由0由02cot2kkT由于数字频域存在 折叠频 率，对

60、于实数响应的数字滤波器因此有效的数字域仅是有效的数字域仅是 高通变换的计算步骤和低通变换一样。但在确定模拟原型预畸的临界频率时，应采用是镜象部分高通也仅指这一段的高端高通也仅指这一段的高端例例【例4.4.2】 设计一数字高通滤波器，它的通带为0.61kHz(fs=2kHz)，通带允许波动0.5dB，阻带衰减在400Hz的频带内至少为20dB。解法一：利用切比雪夫滤波器模型的双线性变换设计1）确定数字高通滤波器的技术指标2）预畸变模拟边界频率确定模拟低通原型指标260020.62000ppf T240020.42000rrf T0.5dB20AtdB0.6cot()0.7265222pTT 0.