数字信号处理第6章.ppt

第6章 无限脉冲响应数字滤波器的设计,6.1 数字滤波器的基本概念 6.2 模拟滤波器的设计 6.3 用脉冲响应不变法设计IIR数字低通滤波器 6.4 用双线性变换法设计IIR数字低通滤波器 6.5 数字高通、 带通和带阻滤波器的设计 习题与上机题,6.1 数字滤波器的基本概念 所谓数字滤波器，是指输入、输出均为数字信号，通过数值运算处理改变输入信号所含频率成分的相对比例，或者滤除某些频率成分的数字器件或程序。因此，数字滤波的概念和模拟滤波相同，只是信号的形式和实现滤波方法不同。正因为数字滤波通过数值运算实现滤波，所以数字滤波器处理精度高、稳定、体积小、重量轻、灵活、不存在阻抗匹配问题，可以实现模拟滤波器无法实现的特殊滤波功能。如果要处理的是模拟信号，可通过A/DC和D/AC，在信号形式上进行匹配转换，同样可以使用数字滤波器对模拟信号进行滤波。 ,1数字滤波器的分类 按照不同的分类方法，数字滤波器有许多种类，但总起来可以分成两大类： 经典滤波器和现代滤波器。经典滤波器的特点是其输入信号中有用的频率成分和希望滤除的频率成分各占有不同的频带，通过一个合适的选频滤波器滤除干扰，得到纯净信号，达到滤波的目的。例如，输入信号x(t)中含有干扰，其时域波形和频谱图分别如图6.1.1(a)、(b)所示，由图可见，信号和干扰的频带互不重叠，可用图6.1.1(c)所示低通滤波器滤除干扰, 得到纯信号，如图6.1.1(d)所示。,图6.1.1 用经典滤波器从噪声中提取信号,但是，如果信号和干扰的频谱相互重叠，则经典滤波器不能有效地滤除干扰，最大限度地恢复信号，这时就需要现代滤波器，例如维纳滤波器、卡尔曼滤波器、自适应滤波器等最佳滤波器。现代滤波器是根据随机信号的一些统计特性，在某种最佳准则下，最大限度地抑制干扰，同时最大限度地恢复信号，从而达到最佳滤波的目的。本书仅介绍经典滤波器的设计分析与实现方法，而现代滤波器属于随机信号处理范畴，已超出本书学习范围。 ,经典数字滤波器从滤波特性上分类，可以分成低通、高通、带通和带阻等滤波器。它们的理想幅频特性如图6.1.2所示。这种理想滤波器是不可能实现的，因为它们的单位脉冲响应均是非因果且无限长的，我们只能按照某些准则设计滤波器，使之在误差容限内逼近理想滤波器，理想滤波器可作为逼近的标准。另外，需要注意的是, 数字滤波器的频率响应函数H(ej)都是以2为周期的，低通滤波器的通频带中心位于2的整数倍处，而高通滤波器的通频带中心位于的奇数倍处，这一点和模拟滤波器是有区别的。一般在数字频率的主值区, 描述数字滤波器的频率响应特性。,图6.1.2 理想低通、高通、带通和带阻滤波器幅度特性,根据滤波器对信号的处理作用又将其分为选频滤波器和其他滤波器。上述低通、高通、带通和带阻滤波器均属于选频滤波器，其他滤波器有微分器、希尔伯特变换器、频谱校正等滤波器。 滤波器可用于波形形成、调制解调器、从噪声中提取信号(见图6.1.1)、信号分离和信道均衡等。所以学习滤波器的设计与实现是必不可少的。运行本书程序集中的绘图程序fig611b.m可以清楚地观察用滤波器分离载波频率不同的两路双边带信号的原理。,2 数字滤波器的技术指标 常用的数字滤波器一般属于选频滤波器。假设数字滤波器的频率响应函数H(ej)用下式表示： H(ej)=|H(ej)|ej() 式中，|H(ej)|称为幅频特性函数; ()称为相频特性函数。幅频特性表示信号通过该滤波器后各频率成分振幅衰减情况，而相频特性反映各频率成分通过滤波器后在时间上的延时情况。因此，即使两个滤波器幅频特性相同，而相频特性不同，对相同的输入，滤波器输出的信号波形也是不一样的。,一般选频滤波器的技术要求由幅频特性给出，对几种典型滤波器(如巴特沃斯滤波器)，其相频特性是确定的，所以设计过程中，对相频特性一般不作要求。但如果对输出波形有要求，则需要考虑相频特性的技术指标，例如波形传输、图像信号处理等。本章主要研究针对幅频特性指标的选频滤波器设计。如果对输出波形有严格要求，则需要设计线性相位数字滤波器，这部分内容在第7章介绍。 ,对于图6.1.2所示的各种理想滤波器，我们必须设计一个因果可实现的滤波器去近似实现。另外，也要考虑复杂性与成本问题，因此实用中通带和阻带中都允许一定的误差容限，即通带不是完全水平的，阻带不是绝对衰减到零。此外，按照要求，在通带与阻带之间还应设置一定宽度的过渡带。,图6.1.3表示低通滤波器的幅频特性，p和s分别称为通带边界频率和阻带截止频率。通带频率范围为0|p，在通带中要求(11)|H(ej)|1，阻带频率范围为s|，在阻带中要求|H(ej)|2。从p到s称为过渡带，过渡带上的频响一般是单调下降的。通常，通带内和阻带内允许的衰减一般用分贝数表示，通带内允许的最大衰减用p表示，阻带内允许的最小衰减用s表示。对低通滤波器， p和s分别定义为：,片段常数特性： 对于选频型滤波器，一般对通带和阻带内的幅频响应曲线形状没有具体要求，只要求其波纹幅度小于某个常数，通常将这种要求称为“片段常数特性”。所谓片段，是指“通带”和“阻带”，常数是指“通带波纹幅度1”和“阻带波纹幅度2”，而通带最大衰减p和阻带最小衰减s是与1和2完全等价的两个常数。片段常数特性概念在选频型滤波器设计中很重要，尤其有助于理解IIR数字滤波器的双线性变换设计思想。,图6.1.3所示的单调下降幅频特性，p和s别可以表示为 如果将|H(ej0)|归一化为1，(6.1.3b)和(6.1.4b)式则表示为：,(6.1.3b),(6.1.4b),(6.1.5),(6.1.6),当幅度下降到 时，标记=c，此时 dB，称c为3 dB通带截止频率。p、c和s统称为边界频率，它们是滤波器设计中所涉及到的很重要的参数。对其他类型的滤波器，(6.1.3b)式和(6.1.4b)式中的H(ej0)应改成 ，0为滤波器通带中心频率。,图6.1.3 低通滤波器的幅频特性指标示意图,3 数字滤波器设计方法概述 IIR滤波器和FIR滤波器的设计方法完全不同。IIR滤波器设计方法有间接法和直接法，间接法是借助于模拟滤波器的设计方法进行的。其设计步骤是： 先设计过渡模拟滤波器得到系统函数Ha (s)，然后将Ha(s)按某种方法转换成数字滤波器的系统函数H(z)。这是因为模拟滤波器的设计方法已经很成熟，不仅有完整的设计公式，还有完善的图表和曲线供查阅； 另外，还有一些典型的优良滤波器类型可供我们使用。直接法直接在频域或者时域中设计数字滤波器，由于要解联立方程，设计时需要计算机辅助设计。FIR滤波器不能采用间接法，常用的设计方法有窗函数法、频率采样法和切比雪夫等波纹逼近法。,对于线性相位滤波器，经常采用FIR滤波器。可以证明，FIR滤波器的单位脉冲响应满足一定条件时，其相位特性在整个频带是严格线性的，这是模拟滤波器无法达到的。当然，也可以采用IIR滤波器，但必须使用全通网络对其非线性相位特性进行校正，这样增加了设计与实现的复杂性。 本章只介绍IIR滤波器的间接设计方法。为此，我们先介绍模拟低通滤波器的设计，这是因为低通滤波器的设计是设计其他滤波器的基础。,模拟高通、带通和带阻滤波器的设计过程是： 先将希望设计的各种滤波器的技术指标转换为低通滤波器技术指标，然后设计相应的低通滤波器，最后采用频率转换法将低通滤波器转换成所希望的各种滤波器。 应当说明，滤波器设计公式较多，计算繁杂。但是，在计算机普及的今天，各种设计方法都有现成的设计程序（或设计函数）供我们调用。所以，只要掌握了滤波器基本设计原理，在工程实际中采用计算机辅助设计滤波器是很容易的事。,6.2 模拟滤波器的设计 模拟滤波器的理论和设计方法已发展得相当成熟， 且有多种典型的模拟滤波器供我们选择，如巴特沃斯(Butterworth)滤波器、切比雪夫(Chebyshev)滤波器、椭圆(Ellipse)滤波器、贝塞尔(Bessel)滤波器等。这些滤波器都有严格的设计公式、现成的曲线和图表供设计人员使用，而且所设计的系统函数都满足电路实现条件。,这些典型的滤波器各有特点： 巴特沃斯滤波器具有单调下降的幅频特性； 切比雪夫滤波器的幅频特性在通带或者阻带有等波纹特性，可以提高选择性； 贝塞尔滤波器通带内有较好的线性相位特性； 椭圆滤波器的选择性相对前三种是最好的，但通带和阻带内均呈现等波纹幅频特性，相位特性的非线性也稍严重。设计时，根据具体要求选择滤波器的类型。,选频型模拟滤波器按幅频特性可分成低通、高通、带通和带阻滤波器，它们的理想幅频特性如图6.2.1所示。但设计滤波器时，总是先设计低通滤波器，再通过频率变换将低通滤波器转换成希望类型的滤波器。下面先介绍低通滤波器的技术指标和逼近方法，然后分别介绍巴特沃斯滤波器和切比雪夫滤波器的设计方法。椭圆滤波器的设计理论比较复杂，所以只介绍其MATLAB设计函数，并举例说明直接调用MATLAB函数设计椭圆滤波器的方法。其他滤波器的设计方法请参考文献9。,图6.2.1 各种理想模拟滤波器幅频特性,6.2.1 模拟低通滤波器的设计指标及逼近方法 本书中，分别用ha(t)、a(s)、Ha(j)表示模拟滤波器的单位冲激响应、系统函数、频率响应函数，三者的关系如下： 可以用ha(t)、Ha(t)、Ha(j)中任一个描述模拟滤波器，也可以用线性常系数微分方程描述模拟滤波器。但是设计模拟滤波器时，设计指标一般由幅频响应函数|Ha(j)|给出，而模拟滤波器设计就是根据设计指标，求系统函数Ha(s)。,工程实际中通常用所谓的损耗函数（也称为衰减函数）A()来描述滤波器的幅频响应特性， 对归一化幅频响应函数（本书后面都是针对该情况，特别说明的除外），A()定义如下（其单位是分贝，用dB表示）： （6.2.1）,应当注意，损耗函数A()和幅频特性函数|H(j)|只是滤波器幅频响应特性的两种描述方法。损耗函数的优点是对幅频响应|Ha(j)|的取值非线性压缩，放大了小的幅度，从而可以同时观察通带和阻带频响特性的变化情况。二者的特点如图6.2.2所示。图6.2.2(a)所示的幅频响应函数完全看不清阻带内取值较小（0.001以下）的波纹，而图6.2.2(b)所示的同一个滤波器的损耗函数则能很清楚地显示出阻带60 dB以下的波纹变化曲线。 另外，直接画出的损耗函数曲线图正好与幅频特性曲线形状相反，所以，习惯将A()曲线称为损耗函数(本书中也如此称谓)，如图6.2.2（b）所示。,图6.2.2 幅频响应与损耗函数曲线的比较,模拟低通滤波器的设计指标参数有 p、p、 s和s。其中p和s分别称为通带边界频率和阻带截止频率， p称为通带最大衰减（即通带0，p中允许A()的最大值）， s称为阻带最小衰减（即阻带s上允许A()的最小值）， p和 s的单位为dB。以上技术指标如图6.2.3所示, 图(a)以幅频特性描述, 图(b)以损耗函数描述。,图6.2.3 模拟低通滤波器的设计指标参数示意图,6.2.2 巴特沃斯低通滤波器的设计 1 巴特沃斯低通模拟滤波器设计原理 巴特沃斯低通滤波器的幅度平方函数|Ha(j)|2用下式表示： (6.2.7) 式中，N称为滤波器的阶数。当=0时，|Ha(j)|=1； =c时， ，c是3 dB截止频率。在=c附近，随加大，幅度迅速下降。,幅度特性与和N的关系如图6.2.4所示。幅度下降的速度与阶数N有关，N愈大，通带愈平坦，过渡带愈窄，过渡带与阻带幅度下降的速度愈快, 总的频响特性与理想低通滤波器的误差愈小。 以s替换j，将幅度平方函数|Ha(j)|2写成s的函数： (6.2.8),复变量s=+j，此式表明幅度平方函数有2N个极点，极点sk用下式表示： (6.2.9) 式中，k=0，1，2，2N-1。2N个极点等间隔分布在半径为c的圆上(该圆成为巴特沃斯圆)，间隔是/Nrad。例如N=3，极点间隔为/3 rad，如图6.2.5所示。,图6.2.4 巴特沃斯低通滤波器幅度特性与,图6.2.5 三阶巴特沃斯滤波器极点分布图,为形成因果稳定的滤波器，2N个极点中只取s平面左半平面的N个极点构成Ha(s)，而右半平面的N个极点构成Ha(s)。Ha(s)的表达式为 (6.2.10) 设N=3，极点有6个，它们分别为,取s平面左半平面的极点s0、s1、s2组成系统函数Ha(s)，即 由于不同的技术指标对应的边界频率和滤波器幅频特性不同，为使设计公式和图表统一，将频率归一化。巴特沃斯滤波器采用对3 dB截止频率c归一化，归一化后的系统函数为,将极点表示式(6.2.13)代入(6.2.12)式，得到Ga(p)的分母是p的N阶多项式，用下式表示： （6.2.15） 归一化原型系统函数Ga(p)的系数bk，k=0，1，N1，以及极点pk，可以由表6.2.1得到。另外，表中还给出了Ga(p)的因式分解形式中的各系数，这样只要求出阶数N，查表可得到Ga(p)及各极点, 而且可以选择级联型和直接型结构的系统函数表示形式，避免了因式分解运算工作。,表6.2.1 巴特沃斯归一化低通滤波器参数,将=s代入(6.2.7)式中，再将|Ha(j)|2代入(6.2.3)式中，得到： (6.2.17) 由(6.2.16)和(6.2.17)式得到：,请注意，如果采用(6.2.19)式确定c，则通带指标刚好满足要求，阻带指标有富余； 如果采用(6.2.20)式确定c，则阻带指标刚好满足要求，通带指标有富余。 总结以上，低通巴特沃斯滤波器的设计步骤如下： (1) 根据技术指标p、 p、s和s，用(6.2.18)式求出滤波器的阶数N。 (2) 按照(6.2.13)式，求出归一化极点pk，将pk代入(6.2.12)式，得到归一化低通原型系统函数Ga(p)。也可以根据阶数N直接查表6.2.1得到pk和Ga(p)。,(3) 将Ga(p)去归一化。将p=s/c代入Ga(p)，得到实际的滤波器系统函数 这里c为3 dB截止频率，如果技术指标没有给出c，可以按照(6.2.19)式或(6.2.20)式求出。 【例6.2.1】 已知通带截止频率fp=5 kHz，通带最大衰减 p=2 dB，阻带截止频率fs=12 kHz，阻带最小衰减 s=30 dB，按照以上技术指标设计巴特沃斯低通滤波器。,(2) 按照(6.2.13)式，其极点为 按照(6.2.12)式，归一化低通原型系统函数为 上式分母可以展开成五阶多项式，或者将共轭极点放在一起，形成因式分解式。这里不如直接查表6.2.1简单，由N=5直接查表得到： ,极点： 0.3090j0.9511， 0.8090j0.5878, 1.0000 归一化低通原型系统函数为 式中, b0=1.0000，b1=3.2361，b2=5.2361，b3=5.2361，b4=3.2361 分母因式分解形式为 以上公式中的数据均取小数点后四位。,(3) 为将Ga(p)去归一化，先求3 dB截止频率c。按照(6.2.19)式，得到： 将c代入(6.2.20)式，得到: ,此时算出的 比题目中给的s小，因此，过渡带小于指标要求。或者说，在s=212 krad/s时衰减大于30 dB，所以说阻带指标有富余量。 将p=s/c代入Ga(p)中, 得到：,2. 用MATLAB工具箱函数设计巴特沃斯滤波器 MATLAB信号处理工具箱函数buttap, buttord和butter是巴特沃斯滤波器设计函数。其5种调用格式如下。 1） Z，P, K=buttap(N) 该格式用于计算N阶巴特沃斯归一化（3 dB截止频率c=1）模拟低通原型滤波器系统函数的零、极点和增益因子。返回长度为N的列向量Z和P，分别给出N个零点和极点的位置，K表示滤波器增益。得到的系统函数为如下形式: （6.2.21）,式中，Z(k)和P(k)分别为向量Z和P的第k个元素。如果要从计算得到的零、极点得到系统函数的分子和分母多项式系数向量B和A,可以调用结构转换函数B, A=zp2tf(Z, P, K)。,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也是二元向量。N和wc作为butter函数的调用参数。,3） N, wc= buttord(wp, ws, Rp, As, s) 该格式用于计算巴特沃斯模拟滤波器的阶数N和3 dB截止频率wc。wp、ws和wc是实际模拟角频率（rad/s）。其他参数与格式2）相同。,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）,由于高通滤波器和低通滤波器都只有一个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阶高通滤波器级联。,【例6.2.2】 调用buttord和butter设计巴特沃斯低通模拟滤波器。要求与例6.2.1相同。 设计程序ep622.m如下: 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); subplot(2, 2, 1); plot(fk/1000, 20*log10(abs(Hk); grid on xlabel(频率(kHz); ylabel(幅度(dB) axis(0, 14, -40, 5) 运行结果: N=5，wc=3.7792e+004，B=7.7094e+022 A =1 1.2230e+005 7.4785e+009 2.8263e+014 6.6014e+018 7.7094e+022,将B和A代入（6.2.23）式写出系统函数为 与例6.2.1计算结果形式相同。滤波器的损耗函数曲线如图6.2.6所示。由图可以看出，阻带刚好满足指标要求，通带指标有富余。这就说明buttord函数使用（6.2.20）式计算3 dB 截止频率。,图6.2.6 程序ep622.m运行输出的损耗函数(例6.2.1的设计结果),6.2.3 切比雪夫滤波器的设计 1 切比雪夫滤波器的设计原理 巴特沃斯滤波器的频率特性曲线，无论在通带还是阻带都是频率的单调减函数。因此，当通带边界处满足指标要求时，通带内肯定会有较大富余量。因此，更有效的设计方法应该是将逼近精确度均匀地分布在整个通带内，或者均匀分布在整个阻带内，或者同时均匀分布在两者之内。这样，就可以使滤波器阶数大大降低。这可通过选择具有等波纹特性的逼近函数来达到。,切比雪夫滤波器的幅频特性就具有这种等波纹特性。它有两种形式： 振幅特性在通带内是等波纹的、在阻带内是单调下降的切比雪夫型滤波器； 振幅特性在通带内是单调下降、在阻带内是等波纹的切比雪夫型滤波器。采用何种形式的切比雪夫滤波器取决于实际用途。图6.2.7(a)和(b)分别画出不同阶数的切比雪夫型和型滤波器幅频特性。,图6.2.7 不同阶数的切比雪夫型和型滤波器幅频特性,我们这里仅介绍切比雪夫型滤波器的设计方法。其幅度平方函数用|Ha(j)|2表示： (6.2.24) 式中，为小于1的正数，表示通带内幅度波动的程度，愈大，波动幅度也愈大; p称为通带截止频率。令=/p，称为对p的归一化频率。CN(x)称为N阶切比雪夫多项式，定义为 ,当N=0时，C0(x)=1； 当N=1时，C1(x)=x； 当N=2时，C2(x)=2x21； 当N=3 时，C3(x)=4x33x。由此可归纳出高阶切比雪夫多项式的递推公式为 (6.2.25) 切比雪夫多项式的特性： (1) 切比雪夫多项式的过零点在|x|1的范围内； (2) 当|x|1时，CN(x)是双曲线函数，随x单调上升。,图6.2.8 四阶切比雪夫型和巴特沃斯低通滤波器的幅频特性比较,去归一化后的系统函数为 (6.2.40) 按照以上分析，归纳出切比雪夫型滤波器设计步骤： (1) 确定技术指标参数 p、p、 s和s。 p是=p时的衰减， s是 =s时的衰减，它们满足,(3) 求归一化系统函数Ga(p)。 为求Ga(p)，先按照(6.2.32)式求出归一化极点pk, k=1，2，N。 (6.2.46) 将极点pk代入(6.2.39)式，得到： ,(4) 将Ga(p)去归一化，得到实际的Ha(s)，即 (6.2.47) 【例6.2.2】 设计低通切比雪夫滤波器，要求通带截止频率fp=3 kHz，通带最大衰减 p=0.1 dB，阻带截止频率fs=12 kHz，阻带最小衰减 s=60 dB。 解 (1) 滤波器的技术要求：,(2) 由(6.2.44)和 (6.2.45)式求阶数N和：,(3) 将极点pk、N和代入(6.2.39)式求Ga(p)： 由(6.2.46)式求出N=5时的极点pi，代入上式，得到： (4) 将Ga(p)去归一化，得到：,2. 用MATLAB设计切比雪夫滤波器 MATLAB信号处理工具箱函数cheb1ap，cheb1ord和cheby1是切比雪夫型滤波器设计函数。其调用格式如下: 1） z, p, k= cheb1ap(N, Rp) 2） N, wpo= cheb1ord(wp, ws, Rp, As) 3） N, wpo= cheb1ord(wp, ws, Rp, As, s) 4） B, A= cheby1(N, Rp, wpo, ftype) 5） B, A= cheby1(N, Rp, wpo, ftype, s),切比雪夫型滤波器设计函数与前面的巴特沃思滤波器设计函数比较，只有两点不同。一是这里设计的是切比雪夫型滤波器； 二是格式2）和3）的返回参数与格式4）和5）的调用参数wpo是切比雪夫型滤波器的通带截止频率，而不是3 dB截止频率。其他参数含义与巴特沃思滤波器设计函数中的参数相同。系数向量B和A与数字和模拟滤波器系统函数的关系由（6.2.22）和（6.2.23）式给出。 MATLAB信号处理工具箱函数cheb2ap, cheb2ord和cheby2是切比雪夫型滤波器设计函数。其调用格式如下:,1） z, p, G= cheb2ap(N, Rs) 该格式用于计算N阶切比雪夫型归一化（阻带截止频率s=1）模拟低通滤波器系统函数的零、极点和增益因子。返回长度为N的列向量z和p，分别给出N个零点和极点的位置。G表示滤波器增益。Rs是阻带最小衰减（dB）。 2） N, wso= cheb2ord(wp, ws, Rp, As) 该格式用于计算切比雪夫型数字滤波器的阶数N和阻带截止频率wso。调用参数wp和ws分别为数字滤波器的通带边界频率和阻带边界频率的归一化值，要求0wp1，0ws1，1表示数字频率（对应模拟频率Fs/2）。,Rp和As分别为通带最大衰减和阻带最小衰减（dB）。当wswp时, 为高通滤波器； 当wp和ws为二元矢量时, 为带通或带阻滤波器，这时wso也是二元向量。N和wso作为cheby2的调用参数。 3） N, wso=cheb2ord(wp, ws, Rp, As, s) 该格式用于计算切比雪夫型模拟滤波器的阶数N和阻带截止频率wso。wp、 ws和wso是实际模拟角频率（rad/s）。其他参数与格式2）相同。,4） B, A=cheby2(N, Rs, wso, ftype) 该格式用于计算N阶切比雪夫型数字滤波器系统函数的分子和分母多项式系数向量B和A。调用参数N和wso分别为切比雪夫型数字滤波器的阶数和阻带截止频率的归一化值（关于归一化），一般调用函数cheb2ord计算N和wso。 5） B, A=cheby2(N, Rp, wso, ftype, s) 该格式用于计算N阶切比雪夫型模拟滤波器系统函数的分子和分母多项式系数向量B和A。调用参数N和wso分别为N阶切比雪夫型模拟滤波器的阶数和阻带截止频率（实际角频率）。 ftype的定义与巴特沃思滤波器设计函数中的ftype 相同。,【例6.2.3】 设计切比雪夫型和切比雪夫型模拟低通滤波器。要求与例6.2.2相同。 解 设计程序ep623.m如下: 例6.2.3设计程序: ep623.m 设计切比雪夫型模拟低通滤波器 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); 计算切比雪夫型模拟低通滤波器系统函数系数,subplot(2, 2, 1); fk=0:12000/512:12000; wk=2*pi*fk; Hk=freqs(B1, A1, wk); plot(fk/1000, 20*log10(abs(Hk); grid on xlabel(频率(kHz); ylabel(幅度(dB) axis(0, 12, -70, 5),运行结果: N=5 切比雪夫型模拟低通滤波器通带边界频率: wp1 =1.8850e+004 切比雪夫型模拟低通滤波器系统函数分子分母多项式系数: B=1.2187e+011 A = 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 椭圆滤波器的设计 椭圆（Elliptic）滤波器在通带和阻带内都具有等波纹幅频响应特性。由于其极点位置与经典场论中的椭圆函数有关，所以由此取名为椭圆滤波器。又因为在1931年考尔（Cauer）首先对这种滤波器进行了理论证明，所以其另一个通用名字为考尔（Cauer）滤波器。椭圆滤波器的典型幅频响应特性曲线如图6.2.10所示。,由图6.2.10(a)可见，椭圆滤波器通带和阻带波纹幅度固定时，阶数越高, 过渡带越窄； 由图6.2.10(b)可见，当椭圆滤波器阶数固定时，通带和阻带波纹幅度越小, 过渡带就越宽。所以椭圆滤波器的阶数N由通带边界频率p、阻带边界频率s、通带最大衰减 p和阻带最小衰减 s共同决定。后面对五种滤波器的比较将证实，椭圆滤波器可以获得对理想滤波器幅频响应的最好逼近，是一种性能价格比最高的滤波器，所以应用非常广泛。,图6.2.10 椭圆滤波器幅频响应特性曲线,椭圆滤波器逼近理论是复杂的纯数学问题，该问题的详细推导已超出本书的范围。只要给定滤波器指标，通过调用MATLAB信号处理工具箱提供的椭圆滤波器设计函数，就很容易得到椭圆滤波器系统函数和零极点位置。 MATLAB信号处理工具箱提供椭圆滤波器设计函数ellipap、 ellipord和ellip。其调用格式如下: 1） z, p, k= ellipap(N, Rp, As) 用于计算N阶归一化（通带边界频率wp=1）模拟低通椭圆滤波器的零点向量z、极点向量p和增益因子k。Rp和As分别为通带最大衰减和阻带最小衰减（dB）。返回长度为N的列向量z和p分别给出N个零点和N个极点的位置。 ,2） N, wpo= ellipord(wp, ws, Rp, As) 用于计算满足指标的椭圆数字滤波器的最低阶数N和通带边界频率wpo，指标要求由参数(wp, ws, Rp, As)给定。参数(wp, ws, Rp, As)的定义与巴特沃思滤波器设计函数buttord中的相应参数相同。 3） N, wpo= ellipord(wp, ws, Rp, As, s) 用于计算满足指标的椭圆模拟滤波器的最低阶数N和通带边界频率wpo。,4） B, A= ellip(N, Rp, wpo, ftype) 当wpo是表示滤波器通带边界频率的标量，而且缺省参数ftype时，该格式返回N阶低通椭圆数字滤波器系统函数的分子和分母多项式系数向量B和A，滤波器通带波纹为Rp dB； 当ftype=high时，返回N阶高通椭圆数字滤波器系统函数系数向量B和A。当wpo是表示带通滤波器通带边界频率的二元向量，而且缺省参数ftype时，该格式返回2N阶带通椭圆数字滤波器系统函数的分子和分母多项式系数向量B和A，滤波器通带波纹为Rp dB。当ftype=stop时，返回2N阶带阻椭圆数字滤波器系统函数系数向量B和A。二元向量参数wpo表示阻带上下边界频率。,5） B, A= ellip(N, Rp, wpo, ftype, s) 计算椭圆模拟滤波器系统函数系数向量B和A。当然，其中的边界频率均为实际模拟角频率值（rad/s）。 【例6.2.4】 设计椭圆模拟低通滤波器。要求与例6.2.2相同。设计程序ep624.m如下: % 椭圆滤波器设计程序 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);,%计算低通模拟滤波器系统函数系数 省去以下绘图部分 运行结果: 椭圆模拟低通滤波器阶数: N=4 模拟低通滤波器通带边界频率: wpo=1.8850e+004 椭圆模拟低通滤波器系统函数分子分母多项式系数: B=0.0010 8.3913e015 2.9126e+007 8.0051e004 1.0859e+017 A= 1 3.3792e+004 9.3066e+008 1.3646e+013 1.0984e+017,滤波器损耗函数如图6.2.11所示。虽然本例中椭圆滤波器阶数是4，但从图6.1.11可以看出，四阶椭圆模拟低通滤波器的过渡带宽度小于7 kHz，比指标要求（9 kHz）窄2 kHz。 而例6.2.3中需要五阶切比雪夫模拟低通滤波器，且其过渡带宽度大于7 kHz。对于本例的设计指标，如果用巴特沃斯模拟低通滤波器，计算所要求的阶数N=7。,图6.2.11 四阶椭圆模拟低通滤波器损耗,6.2.5 五种类型模拟滤波器的比较 前面讨论了四种类型的模拟低通滤波器（巴特沃思、切比雪夫型、切比雪夫型和椭圆滤波器）的设计方法，这四种滤波器是主要考虑逼近幅度响应指标的滤波器，第五种（贝塞尔滤波器）是主要考虑逼近线性相位特性的滤波器。为了正确地选择滤波器类型以满足给定的幅频响应指标，必须比较四种幅度逼近滤波器的特性。为此，下面比较相同阶数的归一化巴特沃思、切比雪夫型、切比雪夫型和椭圆滤波器的频率响应特性。,调用MATLAB滤波器设计函数，很容易验证: 当阶数相同时，对相同的通带最大衰减 p和阻带最小衰减s，巴特沃思滤波器具有单调下降的幅频特性，过渡带最宽。两种类型的切比雪夫滤波器的过渡带宽度相等，比巴特沃思滤波器的过渡带窄，但比椭圆滤波器的过渡带宽。切比雪夫型滤波器在通带具有等波纹幅频特性，过渡带和阻带是单调下降的幅频特性。切比雪夫型滤波器的通带幅频响应几乎与巴特沃思滤波器相同，阻带是等波纹幅频特性。椭圆滤波器的过渡带最窄，通带和阻带均是等波纹幅频特性。,相位逼近情况: 巴特沃思和切比雪夫滤波器在大约3/4的通带上非常接近线性相位特性，而椭圆滤波器仅在大约半个通带上非常接近线性相位特性。贝塞尔滤波器在整个通带逼近线性相位特性，而其幅频特性的过渡带比其他四种滤波器宽得多。 复杂性: 在满足相同的滤波器幅频响应指标条件下，巴特沃思滤波器阶数最高，椭圆滤波器的阶数最低，而且阶数差别较大。所以，就满足滤波器幅频响应指标而言，椭圆滤波器的性能价格比最高，应用较广泛。,由上述比较可见，五种滤波器各具特点。工程实际中选择哪种滤波器取决于对滤波器阶数（阶数影响处理速度和实现的复杂性）和相位特性的具体要求。例如，在满足幅频响应指标的条件下希望滤波器阶数最低时，就应当选择椭圆滤波器。,6.2.6 频率变换与模拟高通、带通、带阻滤波器的设计 高通、带通、带阻滤波器的幅频响应曲线及边界频率分别如图6.2.12(a)、 (b)和(c)所示。 低通、高通、带通和带阻滤波器的通带最大衰减和阻带最小衰减仍用 p和 s表示。图6.2.12中，ph表示高通滤波器的通带边界频率; pl和pu分别表示带通和带阻滤波器的通带下边界频率和通带上边界频率; sl和su分别表示带通和带阻滤波器的阻带下边界频率和阻带上边界频率。 ,图6.2.12 各种滤波器幅频特性曲线及边界频率示意图,从原理上讲，通过频率变换公式，可以将模拟低通滤波器系统函数Q(p)变换成希望设计的低通、高通、带通和带阻滤波器系统函数Hd(s)。在模拟滤波器设计手册中，各种经典滤波器的设计公式都是针对低通滤波器的，并提供从低通到其他各种滤波器的频率变换公式。所以，设计高通、带通和带阻滤波器的一般过程是: (1) 通过频率变换公式，先将希望设计的滤波器指标转换为相应的低通滤波器指标； (2) 设计相应的低通系统函数Q(p)； ,(3) 对Q(p)进行频率变换，得到希望设计的滤波器系统函数Hd(s)。 设计过程中涉及的频率变换公式和指标转换公式较复杂，其推导更为复杂。幸好一些学者已经开发出根据设计指标直接设计高通、带通和带阻滤波器的CAD程序函数，只要根据设计指标直接调用CAD程序，就可以得到高通、带通和带阻滤波器系统函数。前面所提到的MATLAB信号处理工具箱函数butter, cheby1, cheby2和ellip都具有这样的功能。 本节先简要介绍模拟滤波器的频率变换公式，再举例说明调用MATLAB信号处理工具箱函数直接设计高通、带通和带阻滤波器的方法。,后面的设计举例将说明, 如果低通滤波器Q(p)是关于某边界频率的 “归一化低通滤波器”, 则设计计算将大大简化。这里，“归一化低通滤波器”是指关于某个边界频率归一化的低通滤波器, 其系统函数就用Q(p)表示。归一化频率根据设计需要而定，对巴特沃斯滤波器关于3 dB截止频率归一化的系统函数称为巴特沃斯归一化低通原型（记为G(p)），而切比雪夫和椭圆滤波器的归一化低通原型一般是关于通带边界频率p归一化的低通系统函数（即Q(p)的通带边界频率为1）。,为了叙述方便，定义p=+j为Q(p)的归一化复变量，其通带边界频率记为p，称为归一化频率。用Hd(s)表示希望设计的模拟滤波器的系统函数，s=+j表示Hd(s)的复变量。例如，一阶巴特沃斯低通原型系统函数为 显然, 其3 dB截止频率p=1, 是关于3 dB截止频率归一化的。 模拟滤波器设计手册中给出了各种模拟滤波器归一化低通系统函数的参数（零、极点位置，分子、分母多项式系数等）。 ,在虚轴(频率轴)上该映射关系简化为如下频率变换公式: 式中，ph为希望设计的高通滤波器HHP(s)的通带边界频率。频率变换公式（6.2.51）意味着将低通滤波器的通带0, p映射为高通滤波器的通带, ph，而将低通滤波器的通带p, 0映射为高通滤波器的通带ph, 。,（6.2.51）,同样，将低通滤波器的阻带s, 映射为高通滤波器的阻带sh, 0，而将低通滤波器的阻带, s映射为高通滤波器的阻带0, sh。映射关系式（6. 2.50）确保低通滤波器Q(p)通带p, p上的幅度值出现在高通滤波器HHP(s)的通带ph|上。同样，低通滤波器Q(p)阻带s|上的幅度值出现在高通滤波器HHP(s)的阻带s, s上。,所以只要将（6.2.50）式代入（6.1.48）式，就可将通带边界频率为p的低通滤波器的系统函数Q(p)转换成通带边界频率为ph的高通滤波器系统函数: （6.2.52） 【例6.2.5】 设计巴特沃思模拟高通滤波器，要求通带边界频率为4 kHz，阻带边界频率为1 kHz，通带最大衰减为0.1 dB，阻带最小衰减为40 dB。,解 (1） 通过映射关系式（6.2.51），将希望设计的高通滤波器的指标转换成相应的低通滤波器Q(p)的指标。为了计算简单，一般选择Q(p)为归一化低通，即取Q(p)的通带边界频率p=1。则由（6. 2.51）式可求得归一化阻带边界频率为 转换得到低通滤波器的指标为： 通带边界频率p=1，阻带边界频率s=4，通带最大衰减 p=0.1 dB，阻带最小衰减 s=40 dB。 （2） 设计相应的归一化低通系统函数Q(p)。本例调用MATLAB函数buttord和butter来设计Q(p)。,p=1，,（3） 用（6.2.50）式将Q(p)转换成希望设计的高通滤波器的系统函数HHP(s)。本例调用MATLAB函数lp2hp实现低通到高通的变换。lp2hp函数的功能可用help命令查阅。BH, AH=lp2hp(B, A, wph)将系统函数分子和分母系数向量为B和A的低通滤波器变换成通带边界频率为whp的高通滤波器，返回结果BH和AH是高通滤波器系统函数分子和分母的系数向量。实现步骤（2）和（3）的程序ep616.m如下: ,例6.2.5 设计巴特沃思模拟高通滤波器程序: ep625.m wp=1; ws=4; Rp=0.1; As=40; 设置低通滤波器指标参数 N, wc=buttord(wp, ws, Rp, As, s); 计算低通滤波器Q(p)的阶数N和3 dB截止频率wc B, A=butter(N, wc, s); 计算低通滤波器系统函数Q(p)的分子分母多项式系数 wph=2*pi*4000; 模拟高通滤波器通带边界频率wph BH, AH=lp2hp(B, A, wph); 低通到高通转换,由系数向量B和A写出归一化低通系统函数为 由系数向量BH和AH写出希望设计的高通滤波器系统函数为 Q(p)和HHP(s)的损耗函数曲线如图6.2.13所示。,图6.2.13 例6.2.5 所得低通、高通滤波器损耗函数曲线,值得注意的是，实际上调用函数buttord和butter可以直接设计巴特沃思高通滤波器。设计程序ep616b.m如下: 例6.2.5 设计巴特沃思模拟高通滤波器程序: ep625b.m wp=2*pi*4000; ws=2*pi*1000; Rp=0.1; As=40; 设置高通滤波器指标参数 N, wc=buttord(wp, ws, Rp, As, s); 计算高通滤波器阶数N和3 dB截止频率 BH, AH=butter(N, wc, high, s); 计算高通滤波器系统函数HHP(s)分子分母多项式 系数,程序运行结果: N=5 BH=1 0 0 0 0 0 AH=1 5.1073e+004 1.3042e+009 2.0584e+013 2.0078e+017 9.7921e+020 分母多项式系数向量AH与程序ep616.m的运行结果相同。但是，分子多项式系数向量BH与程序ep616.m的运行结果有较大差别，这是由运算误差引起的。由于butter函数采用了归一化处理，所以计算误差小。本程序所得分子多项式系数向量BH与理论设计结果相同。所