《数字信号处理教程——MATLAB释义与实现》第八章

1、1第八章第八章 IIRIIR滤波器设计滤波器设计 2IIRIIR滤波器设计滤波器设计本章主要内容8.1 关于模拟滤波器的预备知识 8.2 模拟滤波器原型的设计 8.3 由模拟滤波器变换为数字滤波器 8.4 模拟滤波器的频带变换 8.5 数字滤波器的频带变换 8.6 IIR数字滤波器的直接设计函数 8.7 IIR和FIR数字滤波器的比较 8.8 MATLAB中的滤波器设计工具 38.1 8.1 模拟滤波器的预备知识模拟滤波器的预备知识 模拟滤波器以是平方幅度响应作为指标。通常把通带最大幅特性归一化为一，故以相对指标的形式给出。对于低通滤波器，在通带和阻带内分别为：22221()1,110(),a

2、pasHjHjA 4关于模拟滤波器的预备知识关于模拟滤波器的预备知识图中为通带波动系数，p和s是通带和阻带边缘频率。A为阻带衰减参数。参数和A化成dB值后可与第七章中的数字滤波器指标Rp和As建立关系。其关系式为： 11011log1010/210RppR20/210101log10AssAAA5关于模拟滤波器的预备知识关于模拟滤波器的预备知识模拟滤波器技术指标以平方幅度响应的形式给出，所以它不包含相位信息。为了从平方幅度响应计算s域系统函数，数学上就是要由复平面频率轴上的函数，推算出全复平面上的函数，那是没有唯一解的。从可实现性考虑，要作两点限制：一是规定函数的结构形式为s的实有理分式；二是

3、此函数(也即其零极点)要关于虚轴对称，满足 2/)()()(jsaaajHsHsH6关于模拟滤波器的预备知识关于模拟滤波器的预备知识于是模拟滤波器设计问题可以表述如下。先找到一个满足幅特性要求的2的有理真分式，然后取出它的左半平面零点zj和极点pk，用这些零极点可以组成一个稳定的模拟滤波器。 0),()()()(pkzjrealkjapszssHNNNNNNcccdddCDJJ2212022120222)()()()(7关于模拟滤波器的预备知识关于模拟滤波器的预备知识在这样的意义下，我们说模拟滤波器的传递函数Ha(s)可由平方幅度函数指标直接设计出来。 图 8.1.2 典型的 Ha(s)Ha(

4、-s)零极点模式。 让左半 平面的零极点构成 Ha(s)，就得到最小相位系统 8IIRIIR数字滤波器的设计方法数字滤波器的设计方法 IIR滤波器设计的基本技术就是利用复变函数映射把模拟滤波器转换成数字滤波器。这种技术的优点是：模拟滤波器设计方法已很成熟。而模拟到数字滤波器的复值映射也是成熟的。设计过程中的关键是进行两种不同的变换：(1)频带变换（低通到高通，带通，带阻等的变换）；(2)模拟/数字变换；根据这两种变换的先后次序，引出两种设计方法：9IIRIIR数字滤波器的两种设计方法数字滤波器的两种设计方法1. 模拟低通原型频带变换模拟/数字变换;2. 模拟低通原型模拟/数字变换频带变换;10

5、8.2 8.2 模拟滤波器原型设计模拟滤波器原型设计 巴特沃思低通滤波器巴特沃思低通滤波器这种滤波器的特征是其通带和阻带都有平坦的幅度响应。N阶滤波器的平方幅度表达式为：其中，N为滤波器的阶数，c是边缘频率（单位：弧度/秒）。NcajHJ2211)()(11模拟滤波器原型的设计模拟滤波器原型的设计取N=1,2,10，得到的平方幅度响应见图8.2.1。 图8.2.1 巴特沃思低通滤波器的平方幅频特性 12模拟滤波器原型的设计模拟滤波器原型的设计为确定系统函数Ha(s)，对J()进行零极点分解，特征方程 的根就是系统的极点，且此滤波器没有零点。 这2N个极点为 它们关于j对称，并绝不会落在虚轴上。

6、222( )()NcaaNNcjHs Hssj220NNcsj1/2(21)/2() ( 1),0,1,21Njk NNkccpjekN 13模拟滤波器原型的设计模拟滤波器原型的设计巴特沃思滤波器平方幅特性的2N个极点分布图。其左半平面的N个极点构成了它的系统函数。 图 8.2.2 巴特沃思低通滤波器平方幅特性的极点分布 14模拟滤波器原型的设计模拟滤波器原型的设计例 8.2.1 给定 试确定模拟滤波器的系统函数。解：根据给定的条件：N3，c0.5。Ha(s)Ha(-s)的极点如上图左所示。取出左半平面的三个极点，可写出系统函数,6411)(62jHa320.5( )(0.250.433)(0

7、.5)(0.250.433)0.125(0.5)(0.50.25)aHssjssjsss15模拟滤波器原型的设计模拟滤波器原型的设计MATLAB提供了一个叫做buttap的函数，来设计N阶归一化（即c=1）巴特沃思模拟原型滤波器，其调用格式为z0,p0,k0=buttap(N)。 去归一化的零极点及增益为z=z0c、p=p0c 及k=k0cN 。要求出系统函数分子分母的系数b和a，可用b=poly(z); p=poly(p)这使得高阶巴特沃斯滤波器的设计计算非常方便，不必用查表的老办法。16模拟滤波器原型的设计模拟滤波器原型的设计计算本例的程序hc821应为：N=input(N= ); % 输

8、入滤波器参数OmegaC=input(OmegaC= );z0,p0,k0 = buttap(N) % 调用buttap函数, p = p0*OmegaC; % 非归一化极点z = z0*OmegaC; % 非归一化零点k = k0*OmegaCN; % 非归一化增益b = real(poly(z);b = k*b % 分子系数向量a = real(poly(p) % 分母系数向量由此程序可求b,a，并得到系统函数的任何形式。17模拟滤波器原型的设计模拟滤波器原型的设计设计过程中，需要由给定的滤波器指标参数 p,Rp, s和As，求得滤波器的阶数N和边缘频率c，即要把滤波器公式变为它的设计方程

9、。令滤波器在通带和阻带边缘频率的幅特性满足指标Rp和As，建立两个联立方程。解出N和c：)/(log2) 110/() 110(log1010/10/10spAsRpN22/10/10101101pscpcsNNRpAs，或18模拟原型的设计例模拟原型的设计例8.2.28.2.2例8.2.2 设通带边缘频率fp=5kHz，最大衰减Rp为1dB，阻带边缘频率fs=12kHz，最大衰减As为30dB，要求设计巴特沃思低通滤波器。解：用手工运算的方法可按下列步骤：（1）。用式(8.2.7)和式(8.2.8)确定阶数N及边缘频率c；（2）。用(8.2.3)式，根据N值求出极点pk分布，取出其左半平面全

10、部极点；（3）。用(8.2.4)式求出滤波器的系数向量。19模拟原型设计模拟原型设计程序程序hc822OmegaP=2*pi*5000; OmegaS=2*pi*12000;Rp=1;As=30;% 确定滤波器的阶次和边缘频率N=ceil( log10(10 . (0.1*abs(Rp)-1)./(10 .(0.1*abs(As) - 1) ) / (2*log10(OmegaP/OmegaS) )OmegaC= OmegaP / ( (10(.1*abs(Rp) - 1)(1/(2*N) ); % 上两条语句可用下面一条语句代替% N, OmegaC =buttord(OmegaP, Ome

11、gaS, Rp, As, s) z0,p0,k0 = buttap(N) %得到左半平面零极点% 求归一化系数b0,a0 b0=k0*real(poly(z0); a0=real(poly(p0)20模拟原型设计模拟原型设计程序程序hc822程序运行的结果为:c=37792弧度?秒。分子分母上出现的cN达到1022的数量级，以致向量a中其它较小的系数无法显示。因此最好用归一化的滤波器系数， b0=1,a0=1,3.2361,5.2361,3.2361,1而把变量s写成（s0=s/c），这时自变量仍为s，所以整个滤波器系统函数仍是非归一化形式。 1 3.2361s025.2361s0 35.23

12、61s0 43.2361s05s0 1)(saHs037792s其中21模拟滤波器原型的设计模拟滤波器原型的设计为了检验此滤波器在p和s处的幅特性，可以用连续系统频率特性计算函数freqs，它是freqz函数在连续域的对等物。典型调用格式为：H =freqs(b,a,w)用归一化参数时：H =freqs(b0,a0,w/OmegaC)其中w为要取的频率向量，不允许取标量，故至少要同时算两个频点。 若对频率范围不清楚，可以用另一种调用格式：H,w =freqs(b,a)系统会自动选200个适当的频点进行计算 。22模拟滤波器原型的设计模拟滤波器原型的设计用以下语句可以检验设计结果：wx0=500

13、0, 12000*2*pi/37792; % 两个频点Hx=freqs(b0,a0,wx0) % 计算该频率上的幅特性dbHx=-20*log10(abs(Hx)/max(abs(H) % 分贝运行程序得到 dbHx = 1.0000 32.1555可见滤波器在OmegaP处的Rp=1分贝，在OmegaS处的As=32分贝，满足了设计要求。键入无左端变量的语句freqs(b,a)，MATLAB会自动给出系统的波德图 。23模拟滤波器原型的设计模拟滤波器原型的设计 图 8.2.3 例 8.2.2 的滤波器幅频和相频特性 24切比雪夫切比雪夫I I模拟原型的设计模拟原型的设计切比雪夫滤波器具有波动

14、性。它有两种类型：切比雪夫型在通带中具有等波动响应，而切比雪夫型在阻带中具有等波动响应。切比雪夫型滤波器的平方幅度响应为：其中N为滤波器的阶数，为通带波动系数，它与Rp有关,TN(x)是N阶切比雪夫多项式 。cNaTjHJ22211)()(xxxxNxTN1),(cosh(cosh10),(coscos()(1125切比雪夫切比雪夫I I模拟原型的设计模拟原型的设计切比雪夫滤波器的等波动响应是由多项式TN(x)造成的。它的主要性质是:（a）0 x1时，TN(x) 在 -1至1之间振荡；（b）0 x时，TN(x)单调增至。为了确定一个因果的和稳定的系统，必须找出J(2)的极点，它们应该是的根。然

15、后把左半平面的极点分配给Ha(s)。0122cNjsT26切比雪夫切比雪夫I I模拟原型的设计模拟原型的设计切比雪夫型滤波器的这些极点落在长轴为bc，短轴为ac的椭圆上。如右图所示。它也没有零点，因此滤波器的归一化系统函数为： 图 8.2.5 三阶切贝雪夫滤波器的极点分布 NiipNapssH11)(21)(27切比雪夫切比雪夫I I模拟原型的设计模拟原型的设计MATLAB提供了切比雪夫型的设计函数：z0,p0,k0=cheblap(N,Rp)用来设计一个阶数为N，通带波动为Rp的归一化切比雪夫型原型滤波器，它在数组z0中返回零点，数组p0中返回极点，并且返回增益k0。如果要得到非归一化的滤波

16、器参数，可把p0和z0乘以c。k0-则乘以cN。但最好还是求出b0=k0poly(z0),a0=poly(p0),及s0=s/OmegaC)来表示系统函数。要调用cheb1ap函数，必须先知道滤波器的阶数N，通带波动Rp和边缘频率c。 28切比雪夫切比雪夫I I模拟原型的设计模拟原型的设计在实际设计时，需要导出切比雪夫型原型滤波器的设计公式，以便根据给定的滤波器技术指标p,Rp,s和As，求得其阶数N，波动系数和边缘频率c，见式 (8.2.16) (8.2.19)。MATLAB已把这些设计公式也归纳成一个函数cheb1ord，其调用格式为：N, OmegaC = cheb1ord(OmegaP

17、, OmegaS, Rp, As, s)它与buttord调用格式相似，用以求得切比雪夫型原型的阶数N和截止频率c，s表示求的是模拟原型（不是数字滤波器）的阶次。29切比雪夫切比雪夫I I原型设计例原型设计例8.2.38.2.3按例8.2.2的指标，fp=5kHz，Rp=1dB，fs=12kHz，As=30dB，设计切贝雪夫I型低通滤波器。解：首先计算需要的参数,按公式（8.2.16）到（8.2.19）依次算得N,c,a,b等，将非常繁琐。我们直接调用MATLAB语句如下：OmegaP=2*pi*5000;OmegaS=2*pi*12000; Rp=1; As=30; N, OmegaC =

18、cheb1ord(OmegaP, OmegaS, Rp, As, s)z0,p0,k0 = cheb1ap(N,Rp) b0 =k0* real(poly(z0); a0 = real(poly(p0) 30切比雪夫切比雪夫I I原型设计例原型设计例8.2.38.2.3运行程序(注意b0,a0语句中的real的作用)得到：N=4, OmegaC =31416, b0 = 0.2457a0 = 1.0000 0.9528 1.4539 0.7426 0.2756因此此切比雪夫型模拟滤波器的传递函数为：画出其幅频特性和极点分布如图所示。可以计算出其通带、阻带的波动各为0.9998和41dB。131

19、4167426. 0314164539. 1314169528. 0314162457. 0)(234sssssHa31切比雪夫切比雪夫I I原型设计例原型设计例8.2.38.2.3 图 8.2.6 切比雪夫型模拟低通滤波器原型的对数幅频特性和极点分布 可以看出此滤波器通带内的波动(1dB)和呈椭圆分布的极点(右图）。32切比雪夫切比雪夫IIII模拟原型的设计模拟原型的设计切比雪夫型的平方幅度响应为：切比雪夫型滤波器通带幅特性是单调的，而阻带是等波动的.这种滤波器在s平面上既有极点，又有零点。它的公式推导比I型更加复杂，故从略。可以直接利用MATLAB函数z0,p0,k0=cheb2ap(N,

20、As)来设计阶数为N，阻带衰减为As的切比雪夫原型滤波器 。同样，要先由cheb2ord函数根据给定的指标p,Rp,s和As，求得滤波器的阶数N和边缘频率c。 122211)()(cNaTjHJ33切比雪夫切比雪夫IIII原型设计例原型设计例8.2.48.2.4按例8.2.2的指标，fp=5kHz，Rp=1dB，fs=12kHz，As=30dB，设计切贝雪夫II型低通滤波器。解：程序hc824几乎和hc823相同，只改变了几个变元： 输入设计指标语句相同N, OmegaC = cheb2ord(OmegaP, OmegaS, Rp, As, s) % 调用切贝雪夫II设计函数% 调用切贝雪夫I

21、I原型函数,得左半平面零极点z0,p0,k0 = cheb2ap(N, As) 34切比雪夫切比雪夫IIII原型设计例原型设计例8.2.48.2.4根据算出的系数，滤波器的传递函数为写成级联型为： 其中 s0=(s/57146)为归一化拉普拉斯算子根据所得结果可以绘出滤波器的幅频特性曲线和零极点分布图如下。2530. 008149. 005653. 107648. 10100)(23424sssssssHa6003. 003673. 1 . 004214. 003976. 001716. 106848. 00)(2222sssssssHa35切比雪夫切比雪夫IIII原型设计例原型设计例8.2.

22、48.2.4从图上可以看出切比雪夫型滤波器在阻带的等幅波动的特性，其峰值为30分贝。并看出它的四个零和四个极点的分布。 图 8.2.8 例 8.2.4 中切比雪夫 II 型模拟低通滤波器原型的对数幅频特性和零极点分布 36椭圆型模拟原型的设计椭圆型模拟原型的设计椭圆滤波器的通带和阻带均展现出等波动响应，它们的幅度响应特性与FIR等波动滤波器类似。因此对于给定的指标，它们可使阶数N最小。椭圆滤波器的平方幅度响应为：其中N是阶数，是通带波动（它与Rp有关），UN( )是N阶雅可比椭圆函数，这种滤波器数学形式复杂，设计时必须现成的程序或表格。先用名为ellipord的设计函数求N和c 再用名为ell

23、ipap的函数计算滤波器参数， 。 cNaUjHJ22211)()(37椭圆滤波器原型设计例椭圆滤波器原型设计例8.2.58.2.5按例8.2.2的指标，fp=5kHz,Rp=1dB,fs=12kHz，As=30dB，设计椭圆型低通滤波器。解：解题程序名为hc825。它和前几个程序的语句相仿，只有下面两句不同。N, OmegaC = ellipord(OmegaP, OmegaS, Rp, As, s)z0,p0,k0 = ellipap(N,Rp,As)算得N=3,它有三个极点,两个零点。系统函数为：画出的幅频特性和零极点分布如下图5687. 002460. 109701. 00816. 3

24、0)(232sssssHa38椭圆滤波器原型设计例椭圆滤波器原型设计例8.2.58.2.5由图可见，椭圆滤波器在通带和阻带内的幅特性都是等波动的。它只要三阶就达到了指标。 它也兼有零点和极点，它的零点也在虚轴上。 图8.2.10 椭圆型模拟低通滤波器原型的对数幅频特性和零极点分布 39模拟滤波器原型的设计模拟滤波器原型的设计滤波器原型的相位响应滤波器原型的相位响应椭圆滤波器提供了最优的平方幅度响应，但通带上的相位响应非线性较大，在某些情况下这是个重要的问题。巴特沃思滤波器的幅度响应是单调的，往往需要更高的阶数。然而它的相位响应比较线性。切比雪夫滤波器的相位特征介于两者之间。因此，在实际应用中，

25、巴特沃思滤波器、切比雪夫滤波器和椭圆滤波器都各有其适用的场合。选择哪种类型取决于实际问题的要求。408.3 8.3 滤波器由模拟变换为数字滤波器由模拟变换为数字 为了利用模拟滤波器的设计结果来求相应的数字滤波器，可以用把s平面映射到z平面的方法,它使Ha(s)变换成等价的数字系统函数H(z)，这种映射(变换)关系，必须满足三条基本要求：第一， H(z)的频率响应要能模仿Ha(s)的频率响应，即s平面虚轴应映射到z平面的单位圆上。第二，因果稳定的Ha(s)应能映射成因果稳定的H(z)。也就是s域的左半平面必须映射到z平面单位圆的内部。第三，变换前后的滤波器在时域和频域的主要特征应尽可能接近。 4

26、1脉冲响应不变法脉冲响应不变法脉冲响应不变法是使数字滤波器的单位脉冲响应序列h(n)模仿模拟滤波器的冲激响应ha(t)。2.1节中讨论过这两者之间量纲的不同，要使两者建立等价关系，应当满足:连续系统脉冲响应由s平面上的各个极点pak的分量组成，而离散系统脉冲响应由z平面上的各个极点pdk的分量组成。(下标a和d分别表示模拟和数字域)。因为s与z之间的映射关系为z=esT,所以有 。( )()(8.3.1)ah nT h nTTpdkakep42脉冲响应不变法脉冲响应不变法复平面s与z之间的映射关系 如右图。 图8.3.1 脉冲响应不变法的映射关系 1,ln ,sTzeszTs平面上每一条宽度为

27、2/T的横条都将映射到整个z平面上，每一横条的左半边映射到z平面单位圆以内，右半边映射到z平面单位圆以外，而s平面虚轴j映射到z平面单位圆上。 T243脉冲响应不变法脉冲响应不变法szsz映射图映射图右图给出了s平面的等实部线和等虚部线到z平面的映射图。图中有意把等实部线的长度取为从 -33，比 -小一些，其相应的z平面上的圆就缺了一个角。 图 8.3.2 脉冲响应不变法映射关系的内部细节 44脉冲响应不变法脉冲响应不变法szsz映射图映射图s平面上的等实部线映射到z平面上的同心圆。s平面上长为2/T的虚轴则映射到z平面的单位圆上。虚轴的长度向正负无限延长时，z平面上的映射将在单位圆上不断转圈

28、，反映出多点对一点的非单值映射。 左半s平面全部映射到单位圆的内部。说明因果性和稳定性在映射中保持不变。s平面上的等虚部线，也即是等线，它们在z平面上的映射是等相角的径向直线。这说明模拟域中的频率，是以线性关系=T映射到数字域的。45脉冲响应不变法脉冲响应不变法脉冲响应不变法是一种时域响应等价的方法。它以变换前后的模拟和数字滤波器脉冲响应等价性为基础，得到相应的数字滤波器。实现这种等价映射的步骤如下：（1）。根据已知的模拟滤波器传递函数，用拉普拉斯反变换求出它的冲激响应；（2）。对此冲激响应进行采样，再乘以T，得到等价的脉冲响应序列，（3）。对此脉冲序列求z变换，得到其等价数字系统的传递函数。

29、 46脉冲响应不变法脉冲响应不变法用脉冲响应不变法设计数字滤波器的步骤为： 按给定的数字指标p, s，确定采样周期T并选择模拟频率p /T, s /T 。根据指标，设计模拟滤波器Ha(s)。 展开 ，作反变换求出ha(t)对ha(t)进行采样，并乘以T，得到等价的脉冲响应序列he(n)=Tha(nT)对he(n)求z变换,1.得到 NkakkapsRsH1)(11( )1akNkep TkT RHzez47脉冲响应不变法脉冲响应不变法实际计算时从第3步求出Rk和pdk后，可免去中间正反z变换，直接写出最后一步的结果。为了研究它们在频率特性方面是否等价，对(8.3.7)两端做傅立叶变换，得到可见

30、，等价数字滤波器的频率响应是模拟滤波器频率响应的周期平移之叠加。只有当模拟滤波器的频率响应是限带的，且频带限于折叠频率以内时，才能使等价数字滤波器的频率响应重现滤波器频率响应， ksaTjejejkjHeHeH)()(48脉冲响应不变法算例脉冲响应不变法算例8.3.18.3.1 利用脉冲响应不变法，把转换成等价的数字滤波器H(z)，设T0.1秒。解： 首先，用部分分式展开Ha(s)：极点p1= -3 和p2= -2，按（8.3.4）和T0.1，得：手工计算是麻烦易错的，需要开发一个MATLAB函数来实现脉冲响应不变映射。 651)(2ssssHa2132651)(2ssssssHa211121

31、36065. 05595. 1108966. 01 . 0112)(zzzzeTzeTzHTTe49脉冲响应不变法子程序脉冲响应不变法子程序脉冲响应不变法子程序的核心语句如下function bd,ad=impinvar0(ba,aa,Fs)% 将模拟滤波器系数向量变为模拟极点和留数Ra,pa, Ca=residue(ba, aa);T=1/Fs; pd=exp(pa*T);%模拟极点变为数字极点% 用原留数Ra和数字极点pd求得数字滤波器系数bd,ad=residuez(T*Ra, pd, Ca)bd=real(bd), ad=real(ad)% 去除微小虚数 50脉冲响应不变法子程序脉冲响

32、应不变法子程序第一行语句把模拟域的系数行向量ba和aa转变为留数Ra、极点pa和直接项Ca。第三行语句的输入变元为三个，residuez认为它们是数字滤波器的留数Rd、极点pd和直接项Cd矩阵，并把它们转变为系数向量bd和ad。 将例8.3.1 的数据b=1,1;a=1,5,6;T=0.1; 分别作为子程序impinvar0的输入变元，可得与上面手工计算的结果一致的bd和ad。 可以用下面的程序检验变换前后的滤波器的脉冲响应和幅频特性， 51例例8.3.18.3.1频域和时域特性比较频域和时域特性比较subplot(2,1,1), t=0:0.1:3;% 计算模拟系统的脉冲响应，乘以T，画出曲

33、线ha=impulse(ba,aa,t);plot(t,ha*T),hold on% 数字系统的脉冲响应，同样时间轴画出点列h=impz(bd,ad,31);stem(0.1*0:30,h) subplot(2,1,2), w=0:0.1:10*2*pi;% 计算模拟系统的幅频特性，画曲线Ha=freqs(ba,aa,w); plot(w,abs(Ha),hold on% 计算数字系统的幅频特性，画曲线H=freqz(bd,ad,w/Fs); plot(w/Fs,abs(H)52例例8.3.18.3.1频域和时域特性比较频域和时域特性比较 把模拟系统的脉冲响应采样，乘以T后与数字系统的脉冲响应

34、是重合的，而数字滤波器的幅频特性却高于变换前的模拟滤波器。 图8.3.3 脉冲不变法变换前后的系统的时域和频域特性的比较 53例例8.3.18.3.1频域和时域特性比较频域和时域特性比较 本题中的采样频率Fs为10Hz，乃奎斯特频率为31.4弧度/秒。根据采样定理，高于这个频率的模拟频谱将折叠回来叠加在主频谱上，所以在该频率上，数字频率的幅特性约为模拟频率幅特性的两倍。本例中模拟滤波器在Fs/2处的幅特性有0.03左右，所以混叠后的幅度就近0.06。这是脉冲响应不变法的致命缺点，除非原来的模拟滤波器具有很好的低通性能，在乃奎斯特频率处幅特性极小，才能使用。54脉冲响应不变法算例脉冲响应不变法算

35、例8.3.28.3.2用巴特沃思原型设计低通数字滤波器，满足：p=0.2,Rp = 1dB；s=0.3,As = 15dB，设滤波器的采样频率Fs为1000Hz。解：模拟滤波器设计的过程： % 第一段程序：设定数字滤波器指标； % 第二段程序：转换为模拟原型指标； % 第三段程序：模拟巴特沃思原型滤波器计算:由这些程序得出模拟滤波器系数b,a。在上节已介绍过，此处不再重复。55脉冲响应不变法算例脉冲响应不变法算例8.3.28.3.2只要加一条语句就得出数字滤波器的系数bd,ad = impinvar0(ba,aa,1/T); % 也可用工具箱函数bd,ad = impinvar(ba,aa,1

36、/T); C,B,A=tf2par(bd,ad)% 变换为并联二阶环节由算出的数据bd,ad，知设计的是一个六阶的巴特沃思滤波器，其系统函数H(z)的并联形式为：在下图中给出它的频率特性图。 1111212121.87 0.632.16 1.150.29 0.45( )1 0.990.251 1.060.371 1.290.69zzzH zzzzzzz56脉冲响应不变法算例脉冲响应不变法算例8.3.28.3.2 图 8.3.4 用脉冲响应不变法设计的巴特沃斯数字滤波器频率特性 578.3.2 8.3.2 双线性变换法双线性变换法 双线性变换法是使数字滤波器的频率响应模仿模拟滤波器的频率响应的一

37、种变换方法，它克服了多值映射引起频率混叠的缺点。 双线性变换法把整个s平面变换到整个z平面上去，而且也使s的左半平面映射到z平面的单位圆内。使s平面与z平面建立一一对应的关系，消除了多值变换性，也就消除了频谱混叠现象。 它同时也保持了数字滤波器与原来的模拟滤波器有同样的稳定性和因果性。58双线性变换法双线性变换法这种映射的映射函数如下：其中T是参数。通常取采样周期，也可以取其它值。这种变换的另一个名字叫做双线性分式变换，这是因为分式展开后，得到：若固定其中一个变量，则另一个是线性的。或者说，这种变换对于s和z是双线性的。它的复平面映射如图8.3.5所示。 2/12/1,11211sTsTzzz

38、Ts0122zsTszT59双线性变换法双线性变换法s平面的等实部线在z平面上的映射是一组圆；当它们在s平面上向j延伸时，其在z平面上的映射都向 (-1,0)点汇聚。其等虚部线在z平面上的映射是一组向外辐射的弧线。 图8.3.5 双线性变换法映射关系的内部细节 60双线性变换映射的特征双线性变换映射的特征1.s平面虚轴映射到 z平面单位圆上；2.整个左半平面映射成单位圆的内部，因此这是一个稳定的变换。3.虚轴上的(-)点一一对应地映射成单位圆(-/T/T)点，对应关系为非线性的：11222010111222TTTjjzzTTTjj及,2tan2T12tan,2T及逆函数61双线性变换法双线性变

39、换法4. 模拟滤波器和经过采样后的数字滤波器之间的关系，可以写出： H的下标e表示等价，以后将把下标e省略。这个变换关系非常简明，只要把Ha(s)中的s换成 即可。下面用同一个例子来研究它与脉冲响应不变法的不同。11112)()(zzTsaesHzH11112zzTs62双线性变换法算双线性变换法算例例8.3.3 利用双线性变换，把转换成数字滤波器，设T1。解：根据（8.3.11），得到整理后得到651)(2ssssHa11121111112111( )211125 2611aazzzHsHzzzzz1211212 . 0105. 01 . 015. 042023)(zzzzzzzH63双线性

40、变换法的子程序双线性变换法的子程序要找到双线性变换的系数变化规律，才能编出其计算程序。假如系统为N阶，它的分子分母都按s正幂降序排列，最后以常数项结束。将 代入通分，得到(8.3.19)式。它的分子分母都是由多项式的高次乘积构成，可以用conv函数的多次递归调用来实现。根据这样的思路，可以编出双线性变换的子程序。101211121)(NNNMMMMabasaasaasaabasbasbasbasH112(1)1sFzNzszDz64双线性变换法子程序双线性变换法子程序bilinear0编程要点：步骤1：把系数ba,aa变得同长，使分子分母程序公用；步骤2：逐项计算(8.3.18)分子分母中各个

41、N次多项式乘积（内层的两个for循环）；分子分母都有N+1项要计算（外层的for循环),每项算完后就叠加在总结果中。步骤3：将系数归一化。65双线性变换法子程序双线性变换法子程序bilinear0调用方法为：function bd,ad = bilinear0(ba,aa,Fs)其中:Fs为采样频率ba,aa为模拟滤波器分子分母（s正幂）系数多项式；bd,ad为等价数字滤波器分子分母（z负幂）系数多项式；66双线性变换法算例双线性变换法算例8.3.48.3.4 用bilinear函数对例8.3.3中的系统函数Ha(s)进行变换，但把T改为0.1。解：MATLAB程序hc834ba=1,1; a

42、a=1,5,6; T=0.1; Fs=1/T; bd,ad = bilinear(ba,aa,Fs); 也可用bd,ad = bilinear0(ba,aa,Fs); 根据算出的bd,ad,写出滤波器传递函数为：21216047. 05573. 110375. 00040. 00415. 0)(zzzzzH67双线性变换法算例双线性变换法算例8.3.48.3.4把原模拟滤波器与变换后数字滤波器在时域和频域进行比较。注意在进行脉冲响应的比较时，同样要把模拟脉冲响应ha(t)的采样值乘以采样周期T。编成下列程序，其运行结果曲线为：subplot(2,1,1),t=0:0.1:3;ha=impuls

43、e(ba,aa,t);plot(t,ha*T), hold onh=impz(bd,ad,31);stem(0.1*0:30,h,.) subplot(2,1,2), w=0:0.1:10*2*pi;Ha=freqs(ba,aa,w); plot(w,abs(Ha),hold onH=freqz(bd,ad,w/Fs); plot(w/Fs,abs(H)68双线性变换时域和频域比较双线性变换时域和频域比较 图 8.3.6 双线性变换前后的系统的时域和频域特性的比较 69双线性变换法双线性变换法从曲线可以看出，两者的脉冲响应在起始点重合得不好，而在以后的各点是很相符的。两者的幅特性在低频区域形状

44、和峰值大小很相似，只是它们的峰值频率发生了变化。把数字频率除以T，即放大十倍，它的低频段和模拟幅特性几乎重合，而到高频区就差别很大了。在乃奎斯特频率处的数字幅特性为零，这说明频率混叠的现象得到了消除。把图8.3.6与图8.3.4作一比较，容易看出，双线性变换得到的数字滤波器基本保持了原模拟滤波器的幅特性，比脉冲响应不变法优越得多。70双线性变换法双线性变换法需要注意的是双线性变换存在一个频率的非线性偏移的问题。数字频率的最大值为，而模拟频率可以向延伸，两者之间又要保持一一对应的映射。两者之间就必然是非线性的关系。其结果是，一个线性相位的模拟滤波器经双线性变换后就得到非线性相位的数字滤波器，不再

45、保持原有的线性相位了。另外，为了使变换后的数字滤波器具有规定的频带，其模拟原型的频带要根据非线性映射关系预先进行修正。 2tan,(8.3.17)2T 71双线性变换法频率预修正双线性变换法频率预修正例如要求数字带通滤波器的四个边缘频率为1,2,3,4，所对应的模拟滤波器的四个边缘频率应分别为： 利用这组做过修正的模拟频率来设计模拟带通滤波器原型。以后对此模拟原型滤波器采用双线性变换，变成数字滤波器后。它的边缘频率应该正是原先所要求的那组1,2,3,4。 2tan,(1,2,3,4)2iiiT 72双线性变换法频率预修正双线性变换法频率预修正 利用这组做过修正的模拟频率来设计模拟带通滤波器作为

46、原型。然后对此模拟原型滤波器采用双线性变换，变成数字滤波器。这时它的边缘频率应该正是我们原先所设定的。这一预修正过程可从图8.3.8看出。73双线性变换法设计步骤双线性变换法设计步骤用双线性变换法设计IIR数字低通的步骤如下。设技术指标p,s, Rp, 和As已经给定：1. 选择T值：它对设计结果没有影响，但会影响计算的精度。一般，使cT数量级在1左右。2. 对边缘频率p和s预修正；计算p和 s：3. 设计满足指标p, s, Rp, As的模拟原型Ha(s);4. 最后，令 ， 得到H(z)。2tan2,2tan2ssppTT11112)(zzTHzHa74双线性变换法算双线性变换法算例例8.

47、3.5 设计例8.3.2的数字巴特沃思滤波器。设滤波器的采样频率为1000Hz.解：MATLAB程序hc835要点(1) 列出数字滤波器指标:(2) 转换为模拟原型指标OmegaP = (2/T)*tan(wp/2); % 通带频率预修正OmegaS = (2/T)*tan(ws/2); % 阻带频率预修正(3) 模拟巴特沃思原型滤波器计算: （见上节）75双线性变换法算双线性变换法算例例8.3.5(4) 用双线性变换转换为数字滤波器系数；bd,ad = bilinear(ba,aa,Fs); %双线性变换:sos,G = tf2sos(bd,ad) % 二阶级联结构 得到一个六阶的滤波器，根

48、据其系统函数画出频率响应曲线如图所示。注意：bilinear函数也可用归一化参数，此时Fs用归一化频率Fs/OmegaC代入即可。bd1,ad1 = bilinear(ba0,aa0,Fs/OmegaC);76双线性变换法算双线性变换法算例例8.3.5特性特性 图 8.3.7 用双线性变换法设计的巴特沃斯数字滤波器频率特性 77双线性变换法算双线性变换法算例例8.3.6用切比雪夫I型为原型设计数字滤波器，指标同例8.3.5。解：MATLAB程序hc836要点% 滤波器指标变换和予修正语句从略% 模拟切贝雪夫I原型滤波器系数ba,aa计算;% 用双线性变换求出数字滤波器系数bd,ad;bd,ad

49、 = bilinear(ba0,aa0,Fs/OmegaC); % 双线性变换:sos,G = tf2sos(bd,ad) % 变为二阶环节级联结构 得出一个四阶滤波器，系统函数及频率特性见后：)8482. 04996. 11)(6493. 05548. 11 ()1 (0018. 0)(212141zzzzzzH78双线性变换法算双线性变换法算例例8.3.6特性特性 图 8.3.8 用双线性变换法设计的切贝雪夫 I 型数字滤波器频率特性 79双线性变换法算双线性变换法算例例8.3.7用切比雪夫II原型，设计数字滤波器，指标同例8.3.5。解：MATLAB程序hc837% 滤波器指标变换和予修

50、正语句从略% 模拟切贝雪夫II原型滤波器系数ba,aa计算;% 用双线性变换求出数字滤波器系数bd,ad;bd,ad = bilinear(ba0,aa0,Fs/OmegaC); % 双线性变换:sos,G = tf2sos(bd,ad) % 变为二阶环节级联结构 得出一个四阶滤波器，系统函数及频率特性见后：)7476. 03087. 11)(1884. 06039. 01 ()2989. 11)(2133. 01 (1653. 0)(21212121zzzzzzzzzH80双线性变换法算双线性变换法算例例8.3.7特性特性 图 8.3.9 用双线性变换法设计的切贝雪夫 II 型数字滤波器频率

51、特性 81双线性变换法算双线性变换法算例例8.3.8 用椭圆型作原型，设计数字滤波器，指标同例835。解：MATLAB程序hc838% 滤波器指标变换和予修正:同例835，程序语句从略% 模拟椭圆原型滤波器计算，得出系数ba0,aa0；% 双线性变换求数字滤波器系数bd,ad；bd,ad = bilinear(ba0,aa0,Fs/OmegaC); sos,G = tf2sos(bd,ad) % 变为二阶环节级联结构得出一个三阶滤波器，系统函数及频率特性见后：)8612. 04928. 11)(6183. 01 ()4211. 11)(1 (4 .121)(211211zzzzzzzH82双线

52、性变换法算双线性变换法算例例8.3.8特性特性 图 8.3.10 用双线性变换法设计的椭圆型数字滤波器频率特性 838.48.4模拟滤波器的频带变换模拟滤波器的频带变换 本节讨论在模拟域把低通滤波器变换为高通、带通或带阻滤波器的方法，也就是图8.1.3中方法1中的ss的变换。讨论的都是模拟滤波器。 高通、带通和带阻滤波器的传递函数可以通过频带变换，由低通滤波器的传递函数求得。因此不论设计哪一种滤波器，都可以先将该滤波器的技术指标转换为低通滤波器的技术指标，先设计低通滤波器。再通过频带变换，将低通的传递函数转换成所需类型的滤波器传递函数。 84模拟低通到低通的频带变换模拟低通到低通的频带变换低通

53、到低通变换：如果已经有了一个设计好的边缘频率为0c=1的归一化低通滤波器，现在要得到另一个边缘频率为c，而其它指标相同的滤波器，那就只要把原来的滤波器进行一个频带变换，使因此，如果归一化低通传递函数为G(s0)，则去归一化低通滤波器的传递函数H(s)用下式转换： csssGsH/0)0()(0ssc85模拟低通到高通的频带变换模拟低通到高通的频带变换低通滤波器和高通滤波器的幅频特性是反的。要由前者变换成后者，就要把通带从低频区换到高频区，把阻带由高频区换到低频区，如右图所示。不难想象，这可以通过对自变量作倒数变换来实现。86模拟低通到高通的频带变换模拟低通到高通的频带变换令 （8.4.3）其中

54、k由给定的高通滤波器通带边缘频率c决定。将归一化滤波器的通带边缘频率0c=1代入（8.4.3），要求左端为c，由此得知k=c。s=c/s0。因而，低通到高通的变换为 (8.4.4)设k=1，按(8.4.4)式画出的0曲线见图8.4.2。 0/ ksscsGsH0)0()(87模拟低通到高通的频带变换模拟低通到高通的频带变换横坐标为归一化低通滤波器频率0，纵坐标是变换后的高通滤波器频率。从-101的低通通带变换为-1-和1的两段高通频段。低通和高通的边缘频率c此处取为1。 图 8.4.2 低通到高通的频率变换关系 88模拟低通到高通的频带变换模拟低通到高通的频带变换由此，模拟高通滤波器的设计步骤

55、如下： （1）确定高通滤波器的技术指标：通带边缘频率c，通带最大波动Rp，阻带最大衰减As。 （2）按照上述Rp，As及通带边缘频率0c=1 ，设计归一化低通滤波器G(s0)。 （3）按照(8.4.4)式，令G(s0)中s0=c/s,转换成高通滤波器的H(s)。89模拟低通到带通的频带变换模拟低通到带通的频带变换 归一化低通与带通滤波器的幅度特性如右图所示。图中c 1和c2分别是带通滤波器的通带下边频和上边频边缘频率。带通滤波器一般用通带中心频率z和通带带宽B两个参数来表征。其关系式为： 12ccB212,zcc 90模拟低通到带通的频带变换模拟低通到带通的频带变换要寻找一种变换，使得低通滤波

56、器频率轴上0=0的点变为带通滤波器的两个中心频率z，而0=1处分别变为带通滤波器两个通带的上下边频c1和c2，并使其带宽为B。带通滤波器频率特性是正负对称的，故这个变换必须是一对二的映射，它应该是的二次函数，在开方时取正频率和负频率都能满足。研究说明，下列映射能满足要求：(8.4.7) Bz22091模拟低通到带通的频带变换模拟低通到带通的频带变换在B=1时， 按这个函数画出的0曲线如右图所示。这函数由上下两段曲线组成，每一段曲线映射出正负半轴上的一个通带。其中心频率为z,带宽为B。 图8.4.4 带通和低通滤波器的频率映射关系 92模拟低通到带通的频带变换模拟低通到带通的频带变换因此，模拟带

57、通滤波器的设计步骤如下：（1）确定模拟带通滤波器的技术指标：包括 通带中心频率z，通带宽度B，通带最大波动Rp，阻带最大衰减As。（2）按照通带宽度0c=1 ，通带最大波动Rp，阻带最大衰减As。设计归一化低通滤波器G(s0)。（4）由下式直接将G(s0)转换成带通滤波器H(s)BszsssGsH220)0()(93模拟低通到带阻的频带变换模拟低通到带阻的频带变换 归一化低通与带阻滤波器的幅特性如右图。图中c1和c2分别是带阻滤波器的阻带下边频和上边频。带阻滤波器用阻带中心频率z和阻带带宽B两个参数来表征。其关系式为：12ccB212,zcc 94模拟低通到带阻的频带变换模拟低通到带阻的频带变

58、换要寻找一种变换，使得低通滤波器频率轴上0=0的点变为带阻滤波器的两个中心频率z，而0=1处分别变为带阻滤波器两个阻带的上下边频c1和c2，并使其带宽为B。带阻滤波器频率特性也是正负对称的，故这个变换必须是一对二的映射，它也应该是的二次函数，在开方时取正频率和负频率都能满足。研究说明，下列映射能满足要求：(8.4.12) 220zB95模拟低通到带阻的频带变换模拟低通到带阻的频带变换在B=1时，按这个函数画出的0曲线如右图所示。这函数由三段曲线组成，上半平面的曲线映射出正半轴上的一个阻带；下半平面曲线映射出负频率阻带。 图8.4.6 带阻和归一化低通的频率映射关系 96模拟低通到带阻的频带变换

59、模拟低通到带阻的频带变换因此，模拟带阻滤波器的设计步骤如下：（1）确定模拟带阻滤波器的技术指标：包括阻带中心频率z，阻带宽度B，通带最大波动Rp，阻带最大衰减As。（2）按照通带宽度0c=1 ，通带波动Rp，阻带衰减As。设计归一化低通滤波器G(s0)。（3）由下式直接将G(s0)转换成带通滤波器H(s)220)0()(zsBsssGsH97模拟滤波器的频带变换模拟滤波器的频带变换上述的这四种频带变换，都要将变量s0代换为s的多项式分式，手工运算非常麻烦。所以要编成子程序，其编程的原理和双线性变换相仿。MATLAB提供了四个函数lp2lp, lp2hp, lp2bp, lp2bs来完成这些频带

60、变换。以lp2bs（低通到带阻）为例，调用格式为：numT,denT = lp2bs(num,den,OmegaZ,B) 其中：num,den为低通原型的分子分母系数 OmegaZ,B 为带阻的中心频率z和阻带宽度BnumT,denT为带阻滤波器的分子分母系数98模拟滤波器频带变换算例模拟滤波器频带变换算例8.4.18.4.1设已根据指标要求，设计好了归一化低通滤波器为三阶巴特沃斯型，其传递函数为试设计下列模拟滤波器： (1）通带为10Hz的低通滤波器；（2）通带下边频为10Hz的高通滤波器；（3）中心频率为10Hz，带宽为2Hz的带通滤波器；（4）中心频率为10Hz，带宽为2Hz的带阻滤波器