数字信号处理89人文科技.ppt
第8章 数字滤波器的结构 在得到H(z)之后,还应该将其实现。H(z)实际上是一个数学表达式,它的不同的数学表示形式对应着不同的算法结构。这就是说,同一个系统函数H(z),可以用不同的算法结构来实现。不同的算法结构,使得滤波器具有不同的性能。与算法结构有关的滤波器性能主要指以下几个方面。,1相乘、相加、延迟等运算的运算量,这关系到运算的复杂程度和运算速度。2所用的加法器、乘法器、延迟器的数目,这关系到系统的复杂度和系统成本,而延迟器的数目直接影响系统所需存储器的多少。,3系统频率特性对于乘法器系数变化的灵敏度,这是数字滤波器的非常重要的性能。在设计H(z)时,并没有考虑系数的精度问题,即是在无限精度的条件下进行设计的。但是,在实现H(z)时,这些系数都要化为二进制数以有限的字长即有限的精度来进行运算。于是,实现时H(z)产生的系数误差使得所设计的零极点位置发生变化,系统的频率特性也就发生变化。,在相同的二进制字长的条件下,或者说在相同的系数误差范围的情况下,不同的滤波器结构所引起的频率特性的变化程度不同,频率特性变化大的结构对于乘法器系数变化的灵敏度就高,或者说这种结构的灵敏度特性差,反之则是灵敏度特性好。,对于IIR滤波器,结构的灵敏度特性不仅影响系统的频率特性,还影响系统的稳定性,这是因为滤波器实现时系数的误差若使得极点的位置变化较大,就有可能从单位园内移出,这将使系统不稳定。,综上所述,算法结构的选择对于一个系统的实现非常重要。一个确定的数字系统有其确定的差分方程、确定的单位抽样响应h(n)以及确定的系统函数H(z),而同一个系统函数H(z),却可以画出不同的算法结构。用什么样的结构来实现H(z),也是数字滤波器的一个很重要的问题,而信号流图是数字网络结构的一种非常方便、直观、有效的表示。,8.1 数字网络的信号流图 8.1.1 信号流图及其有关概念信号流图是由连接节点的有向线段构成的网络,它是表示信号流通的几何图形。信号流图清楚地表示了系统的算法结构,通过信号流图,可以对系统进行有效的分析,还可以方便地求出系统函数。下面结合图8.1来介绍信号流图的有关概念。,图8.1 一个信号流图,1节点:信号流图中每一节点都对应一个变量,或者说代表一个信号。节点又叫做节点变量。2支路与支路传输:支路是连接两个节点的有向线段。支路旁边标注的系数(或称加权)叫做支路传输,它们起着相乘的作用。一种支路传输就是乘法器系数;另一种实际上表示延迟,如图中支路X0X1的传输z-1,有X1=z-1 X0,即x1(n)=x0(n-1)。,3.源(节)点:对于一个节点,流入该节点的信号叫输入,流出该节点的信号叫输出。若一个节点只有输出支路与之相连接,则称之为源节点,或输入节点。4.汇点:若一个节点只有输入支路与之相连接,则称之为汇点,或输出节点。,5.混合节点:若一个节点既有输入支路与之相连接,又有输出支路与之相连接,则称之为混合节点。6.开路径:也叫通路,即是从某一节点出发,沿支路方向连续经过一些支路而中止到另一节点上的路径。注意,每一节点只通过一次。,7闭路径:也叫自环或者环路,即是从某一节点出发,沿着支路方向,连续经过一些支路又中止在出发节点的路径。注意,途中各节点只通过一次。8节点变量的值:设连接节点Xi 与Xj的支路XjXi的支路传输为tji,则节点变量Xi的值为:(8.1)即节点变量的值等于流入该节点的全部信号的叠加。要特别注意,计算节点变量值时不要考虑输出支路。,图8.1中,节点变量X1、X2、X3、X4之值分别为:X1=z-1(8.2)X2=aX0+bX1+eX3(8.3)X3=cX1+dX2+fX3(8.4)X4=X2(8.5)而源节点X0 无值可言。,8.1.2 解代数方程组求节点变量之值将流图中各节点变量的值都表示出来,就组成了一代数方程组。如果流图中除了源点之外有M个节点,那末可以得到含有M个未知数由M个方程组成的线性方程组。但是为了简化方程组的求解过程,应当尽量减少需要求解的节点个数。,一般情况下,这样的线性方程组有且只有一组解,另外,往往也不需要解出每个节点变量,因此,采用克莱姆法则来求解比较方便。,例8.1 用代数方程组求解法求图8.1的系统函数 H=X4/X0。解:由于X4=X2,因此只需要解出节点变量X2。由上面已有的节点之值表示式(8.2)、(8.3)和(8.4),可以得到方程组的矩阵表示形式:,根据克莱姆法则,可求得:X2=2/,这里为系数矩阵行列式:,而,于是可求得系统函数:,看得出,从X2 到X4这条传输为1的支路是从X2延伸出来的,实际上,可以直接从X2 获取输出,并且直接求出节点X2与源点X0之间的系统函数。以此类推,可以根据信号流图求出从源点到其它任何一个节点的传输函数也即系统函数。,8.1.3 化简信号流图求系统函数如果一个信号流图最终被化简为图8.2所示的最简形式,那末支路传输H显然就是系统函数,因为Y=HX0,H=Y/X0。信号流图化简的依据是节点变量的表示式和代数方程的恒等关系,主要有下面三种情况。,图8.2 信号流图的最简形式,Y,X0,H,8.1.3.1 合并支路串联的支路通过支路传输相乘来合并,并联的支路通过支路传输相加来合并。例如图8.3中:H=a(bdeg+cfh)ij。,图 8.3 支路的合并,8.1.3.2 消除节点 图8.4 节点的消除 图8.4中,由于 X2=bX1=b(aX0)=abX0 以及 X3=cX1=c(aX0)=acX0,因此消去了节点X1。,8.1.3.3 消除自环,图8.5 自环的消除,由图8.5左边的流图有:X2=bX1=b(aX0+cX2),于是X2=abX0+bcX2,即为中间的流图;再将含X2的项移到等式左边,得到:,即为右边的流图,此时消除了自环,流图化为最简。,例8.2 如图8.6所示,利用自环消除的规则首先消除了左边流图的自环和节点X2,然后经支路合并后得到了最简形式。此例的流图还可以首先消除节点X3来进行简化,如图8.7所示。,图 8.6 例8.2流图的化简(一),图 8.7 例8.2流图的化简(二),例8.3 信号流图如图8.8(a)所示,中间的三个节点X2、X3、X4以及三个自环按图中所示步骤逐一消除。应注意,在流图(c)中,节点X3后面有一条输入支路的权值要受X3上自环消除的影响,事实上,在流图(b)中,有:X3=abX1+bcX3+eX4由此式得:这即为流图(c)中的情形。后面的化简也类似处理。,图 8.8 例8.3流图的化简,最后得:显然,D=X5/X1 就是这个系统的系统函数。,8.1.4 Mason公式 在给出Mason公式之前,首先解释以下名词。1通路传输:通路边界间(即开始节点与终止节点之间)各支路传输之乘积。2环路传输:绕环路一周各支路传输之乘积。,3不接触:两条通路间或两个环路间或一条通路与一个环路之间若无公共节点则称它们互不接触。Mason公式给出了一个信号流图中,从源点到其余任一节点的传输函数(即系统函数)的表达式:(8.6),其中,为流图的行列式:=1-(所有环路传输之和)+(每两个互不接触的环路传输乘积之和)-(每三个互不接触的环路传输乘积之和)+而 gi 是从源点到这一节点的第 i 条通路的通路传输,i 则是此通路流图的余子式:,i=1-(与此通路不接触的各环路传输之和)+(与此通路不接触的每两个互不接触的环路传输乘积之和)-(与此通路不接触的每三个互不接触的环路传输乘积之和)+例8.4 由一电路网络所得到的信号流图如图8.9所示,试用Mason公式求其系统函数:H=V4/Vg。,图 8.9 一个电路网络的信号流图,解:这实际上是一个模拟网络的信号流图,但是,从流图的观点,并不需要考虑节点和支路传输的物理意义,也就是说,模拟网络的信号流图可以与数字网络的信号流图同样地处理。图中,每条支路用其编号来代表,如、等等;另外,对于同一条支路,有:ZjYj=1,RjGj=1。,为求得Mason公式中的分母,应找出流图中所有环路:a.2,13;b.10,11;c.3,12;d.8,9;e.6,7,8,13;f.8,3,4,11。于是=1-(-Z1Y2-G3R4+-Z1G3+Z1G3-G3R4)+(G3R4Z1Y2-G3R4)为求得Mason公式中的分子,需要找出由Vg到V4的所有通路。各条通路的具体情况如下:,1).,;g1=G3R4;1=1-(-Z1Y2+Y1Z1)=1+Z1Y2-.2).,;g2=Y2Z1Y1R4=Y2R4;2=1.3).,;g3=-Y2Z2G3R4=-G3R4;3=1-Y1Z1=1-.,4).,;g4=G3Z1Y1R4=G3R4;4=1.5).,;g5=Y2Z1(-G3)R4=-Z1Y2G3R4;5=1.6).,;g6=Y2Z2(-G3)Z1Y1R4=-G3R4;6=1.,故可得到:采用Mason公式可以快捷地从一个信号流图得到系统函数,也可以快捷地由已知的系统函数画出信号流图。,8.1.5 信号流图的转置将信号流图转置是对其形式的一种变换,这种变换包括三项操作:1.将输入变量(源点)和输出变量(汇点)交换位置;2.将各条支路反向;3.保持每条支路的支路传输不变。,当信号流图中只有一个源点和一个汇点时,转置后的流图与原流图有相同的系统函数。例8.5 将图8.10(a)的信号流图进行转置,并根据Mason公式说明转置前后的流图具有相同的系统函数。,图 8.10 信号流图的转置,解:在图8.10中,对流图(a)进行转置的三项操作就得到了流图(b),再按习惯将输入画在左边,输出画在右边,就得到流图(c),图(b)和图(c)完全等价。,根据Mason公式,很容易得到这三个流图的系统函数。在这三个流图中,都只有一个环路,而且环路传输都是acz-1;又,三个流图中从源点到汇点都只有一条通路,通路传输为c,并且这条通路与环路接触。于是,三个信号流图具有相同的系统函数,即为:,8.2 IIR数字滤波器的结构 8.2.1 直接型IIR数字滤波器是递归型的线性时不变因果系统,其差分方程可以表示为:(8.7),对上式两边进行z变换,可得:(8.8)图8.11所示的信号流图是对差分方程(8.7)式或者其z变换式(8.8)式的直接实现,因此这种形式的IIR滤波器结构叫做直接型。,图 8.11 直接 I 型,由(8.8)式可以得到IIR滤波器的系统函数:(8.9)显然,分子的系数ai 确定了H(z)的零点,而分母的系数bi 确定了H(z)的极点。,图8.11所示的直接型结构可以看成是两个独立网络的级联,也即第一个网络的输出y1(n)正是第二个网络的输入。第一个网络是:,故可得到第一个网络的系统函数:(8.10),H1(z)与H(z)的分子多项式对应,也就是说,第一个网络实现了滤波器的零点。第二个网络是:故其系统函数是:(8.11),H2(z)的分母正是H(z)的分母,因此第二个网络实现了滤波器的极点。整个滤波器是这两个网络的级联:(8.12)现在,将图8.11中的两个网络的级联次序交换,如图8.12所示。,图8.12 直接型,。,这时,前面的网络即反馈网络为:对这个式子进行z变换后可以得到其系统函数:(8.13),图8.12后面的网络为:经z变换后得到系统函数:(8.14),图8.12的系统函数为:(8.15)(8.15)式与(8.9)式最后结果完全相同。,图8.12也是IIR滤波器的直接型结构,叫直接型。上面的讨论同时也说明了,级联网络总的输入输出关系即总的系统函数与各子网络级联的次序无关;级联网络总的系统函数等于各级联子网络系统函数的乘积。8.2.2 正准型将图8.12中两列传输比为z-1的支路合并为一列(假设M=N),就得到正准型结构,图8.13 IIR滤波器的正准型结构,观察图8.11、图8.12以及图8.13(a)、(b),容易看出,这四个流图中环路的情况以及由输入到输出的通路情况都完全相同(设M=N),根据Mason公式,很快就可以知道,这四个流图的系统函数都如(8.9)式所示;并且,四个流图中的前向通路与系统函数H(z)的分子很好地对应,而流图中的反馈回路与H(z)的分母很好地对应。因此,在这样的情况下,我们很容易根据信号流图写出其系统函数,也很容易根据系统函数画出相应结构的信号流图。,例8.6 已知系统函数:,用正准型结构实现。解:将H(z)化成IIR滤波器系统函数的标准形式:流图如图8.14所示。这里应该特别注意,H(z)的分母的第一项应该化为1,否则容易出错。,图8.14 例8.6的信号流图,IIR滤波器的直接型结构和正准型结构都是对系统差分方程的直接实现,它们的共同缺点是灵敏度特性差,系数误差会引起零极点位置的较大的变化,而正准型比直接型优越的地方是节省一半的存储单元,因此,直接型结构一般不采用,而正准型结构当阶次高(N大于3或者4)时也不采用。,8.2.3 级联型 将(8.9)式的分子、分母多项式进行因式分解:(A=a0),因系数ai、bi都为实数,故零点ci与极点di或是实根,或是共轭复根,而每一对共轭因式又可以合并为一个实系数的二次三项式,于是有:(8.16),单实根因式可以看作二阶因式的特例,设NM,则上式又可写成:(8.17),这里L是由N/2到N范围内的一个整数。子网络Hi(z)的结构如图8.15所示,而整个系统就是由L个这样的二阶子网络级联而成。当然,可能有的子网络中有些系数等于0,此时图8.15所示的结构会得到简化。,图 8.15 一个子网络的结构,由(8.17)式可以清楚地看到,1i和2i确定第i对零点,1i和2i确定第i对极点,也就是说,级联型结构是将零点和极点分散到各个子网络中来实现的,这就使得这种结构的灵敏度特性优于直接型和正准型结构。并且,调整任何一个子网络的零点和极点都不影响其它的零点和极点,即零极点具有独立性。,应该特别指出,零点对与极点对之间共有L!种不同的搭配方式,而对于每一种搭配所得到的L个子网络,又有L!种不同的排列次序。这些不同的方案总的系统函数都相同,但系统特性或者说零极点位置对于系数变化的灵敏度不同。因此,可以选择最优的零极点搭配和级联次序,以降低系统频率特性对于系数变化的灵敏度,提高系统的稳定性。,8.2.4 并联型还可以将(8.9)式的H(z)表示成如下的部分分式展开式:(8.19),将式中每一对共轭的一次因式合并为实系数的二次三项式,并且假设M=N,就得到:(8.20),图8.16 IIR滤波器的并联型结构,由于并联支路的极点也是整个网络的极点,而并联支路的零点却不一定是整个网络的零点,因此并联型的网络结构,可以独立地调整极点位置,但不能独立地控制零点。当然,与级联型一样,并联型结构的灵敏度特性优于直接型和正准型。另外,并联型结构的运算误差比级联型结构小。,8.3 FIR 数字滤波器的结构 FIR数字滤波器的差分方程为:(8.21)并且系统函数为:(8.22),8.3.1 横截型由于横截型结构实际上是对差分方程(8.21)式的直接实现,所以横截型又叫做直接型。另外,由于差分方程:所以横截型还可以叫做卷积型。,图 8.17 横截型结构(一),图 8.18 横截型结构(二),运用Mason公式很容易将系统函数H(z)与图8.17或者图8.18对应起来。显然,这两个流图中都没有环路,因此H(z)的分母为1;而分子是从输入到输出的N条通路传输之和。如果FIR滤波器是线性相位的,那末其冲激响应h(n)满足对称条件,这时,图8.17和图8.18都可以简化。,1.偶对称的情形,有:h(n)=h(N-1)-n 0nN-1 当N为偶数时,有:(8.23)信号流图如图8.19(a)所示,此结构与FIR滤波器一般的直接型结构相比,乘法器减少了一半。,图8.19 线性相位FIR滤波器的直接型结构,当N为奇数时,有:(8.24)信号流图如图8.19(b)所示,乘法器数目也减少了近一半。,2.奇对称的情形,有:h(n)=-h(N-1)-n 0nN-1当N为偶数时,可得:(8.25)只需在图8.19(a)中,在由z-(N-1-n)来的信号旁写一个减号即可得到这种情况下的网络结构。,当N为奇数时,可得:(8.26)此时的网络结构,除了在图8.19(b)中,在由z-(N-1-n)来的信号旁写一减号以外,还应去掉乘法器系数为h(N-1)/2)的支路。,8.3.2 级联型一般情况下,冲激响应h(n)为实数,所以将(8.22)式中的多项式分解因式后,H(z)可以写成若干实系数的二次三项式的乘积,即:(8.27)这里K是(N-1)/2到N-1范围内的某一整数,于是可以得到FIR滤波器的级联型结构。,图8.20 FIR滤波器的级联型结构,这种形式的结构中,每一个子网络控制一对零点,即零点可以独立调整,而且零点位置变化的灵敏度优于横截型结构。但是这种结构所需的乘法运算次数比横截型多。,8.3.3 频率抽样型 非递归型的FIR滤波器也可以用递归算法来实现,这就是频率抽样型结构。事实上,如果按照用频率抽样法设计的系统函数(7.56)式来实现FIR滤波器,那末得到的就是频率抽样型结构。,显然,这样的结构中出现了反馈回路,当然同时也出现了极点。这又如何解释FIR滤波器的极点都集中在z=0呢?(8.28)显然,这个结构主要由两个子网络级联而成,它们的流图分别如图8.21和图8.22所示。,图 8.21 Ha(z)的信号流图,图 8.22 Hb(z)的信号流图,可以看到,Ha(z)在单位圆上均匀分布的N个零点正是Hb(z)的N个极点,因此,当这两个网络前后级联时,它们正好相互抵消;另外,级联结构还使得Ha(z)在z=0处的极点抵消了Hb(z)在z=0的一阶零点。因此,最终的结果是保留了FIR滤波器原有的零极点,即在z=0的N-1阶极点和有限z平面上的N-1个零点。,但是实际上,Ha(z)在单位园上的N个零点是靠延时来实现的,因此能够准确实现;而Hb(z)在单位圆上的极点却是靠复数乘法来实现的,故不能准确实现。既然单位圆上的极点不能被零点完全抵消,滤波器就会出现不稳定现象,因此,应当对上面所述的网络结构进行修正。,令0r1并且r1,用rz-1 来代替(8.28)式的H(z)中的z-1,即有:(8.31)于是,的零点移到了半径为r的圆上:,而并联网络的极点也移到了半径为r的圆上:这样,即使极点不能完全被零点抵消,由于是在单位圆内,也不会导致整个系统不稳定。下面,还要将并联网络部分所含的运算都化为实数运算,为此要利用以下周期性:,还可以得到:设0kN-1,那么也有0N-kN-1,此时有:H*(k)=H(N-k);又令H(k)=|H(k)|ej(k)。现在将并联网络中的第k及第(N-k)子网络合并为一个二阶网络,即有:,(8.33)这里(8.34),若N为偶数,则(8.33)和(8.34)式中k=1,2,N/2-1,故N-k=N-1,N-2,N/2+1。当k=0:(8.35)当k=N/2:(8.36),若N为奇数,则(8.33)和(8.34)式中k=1,2,(N-1)/2,故N-k=N-1,N-2,(N+1)/2,而H0(z)如(8.35)式所示。因此,由(8.31)式,有:当N为偶数,(8.37)当N为奇数,(8.38),当N为偶数时,FIR滤波器的整个频率抽样型的网络结构就如图8.27所示。在这样的结构中,He(z)的N个零点也即并联网络的N个极点都移到了半径为r(r1,r1)的圆上,而且所有的运算都是实数运算。关于最后一点,只需要补充说明H(0)和H(N/2)是实数,实际上,由于h(n)为实序列,故很容易证明 H(0)和H(N/2)也都是实数。,图8.23 He(z)的网络结构,图8.24 Hk(z)的网络结构,图8.25 H0(z)的网络结构 图8.26 HN/2(z)的网络结构,图 8.27 频率抽样型的网络结构(N为偶数),8.4 FIR数字滤波器与IIR数字滤波器的比较 第三部分主要讨论了数字滤波器的原理、设计方法和实现的结构。我们之所以要分别讨论IIR和FIR这两大类滤波器以及它们的各种设计方法,是因为没有一种滤波器也没有一种设计方法在所有的情况下都是最佳的。作为这部分的最后一小节,有必要对这两大类滤波器从各方面进行比较。,关于设计方法,IIR滤波器可以借助于模拟滤波器来进行设计,因此设计过程相对来说比较简单;而FIR滤波器不能借助于模拟滤波器来进行设计,虽然窗口法较简单,但有时还需要作一些迭代运算,实际上,设计FIR滤波器的大多数方法都需要进行大量的迭代运算;但是,FIR滤波器的设计方法具有灵活性,能够逼近更加任意的频率响应;此外,IIR滤波器的设计一般忽略了滤波器的相位响应,而FIR滤波器却可以有精确的线性相位特性。,由于FIR滤波器的极点都集中在z=0处,因此不会出现不稳定的情况,这点比IIR滤波器优越。但是,对于一定的频率响应指标,若用IIR滤波器来实现,一般来说所需的阶数明显低于用FIR滤波器实现所需的阶数。,关于滤波器的结构,无论是IIR还是FIR,直接实现时滤波器的频率特性对于系数变化的灵敏度都比较差,但是直接实现的结构所需的运算次数一般都比级联型和并联型等结构所需的要少。不过,IIR滤波器的极点灵敏度不但影响其频率特性,还影响系统的稳定性,,因此,一般情况下IIR滤波器当其阶次在3、4以上就应该采用级联型或者并联型结构来实现。而FIR滤波器在阶次N为几十甚至上百时也常常用横截型结构直接实现,这是因为对于非递归结构的FIR滤波器,只存在零点灵敏度问题,不会影响稳定性;并且本来FIR滤波器的阶次就比较高,如果再分解为级联型来实现,所需的运算次数会更多。,8.5 用matlab实现数字滤波器的结构 8.5.1 IIR 数字滤波器的结构实现 8.5.1.1 直接型可以直接调用函数filter(b,a,x)来实现直接型结构,其中b和a 分别表示系统函数分子和分母多项式的系数。,8.5.1.2 级联型1sos=zp2sos(z,p,k)zp2sos函数可以用来从已给定的系统函数H(z)直接求出二阶因式,它产生以零、极点形式确定的等效系统函数H(z)的各个二阶因式的系数矩阵sos,sos是一个L6 的矩阵。,该矩阵的第j行包含了第j个二阶因式的分子和分母多项式的系数bij 和aij,i=0、1、2,L表示级联的二阶子系统的数目。因此整个系统函数为:,2sos,g=tf2sos(b,a)该函数用来实现将(8.39)分解成一系列二阶子系统Hj(z)级联的形式,其中Hj(z)的表达式为,3.z,p,k=tf2zp(b,a)函数 z,p,k=tf2zp(b,a)完成将系统函数转换为零极点的形式。注意函数 tf2zp(b,a)要求分母 a 的长度要大于等于分子b的长度。即可以将形如(8.39)式的系统函数变为如下的零极点形式:,例8.7 已知离散系统函数为,将其转换成零极点形式解:b=2 3;a=1 0.4 1;b,a=eqtflength(b,a);%让b和a长度相等z,p,k=tf2zp(b,a)%得到零极点增益形式,z=0-1.5000k=2从运行结果可以写出零极点增益形式的系统函数为:,4.b,a=sos2tf(sos)函数 sos2tf(sos)的功能和tf2sos相反,它用来由二阶子系统构成直接形式的系统函数。,8.5.1.3 并联型利用我们在2.9节介绍的函数residuez就可以实现IIR滤波器的并联型结构,这个函数的调用请见2.9节。这里要补充说明的是,如果系统函数有共轭复数极点,则利用b1,a1=residuez(R1,p1,0)语句可以获得其所对应的实系数二阶分式的分子、分母多项式的系数。其中R1为共轭复数留数所构成的向量,p1为共轭复数极点所构成的向量,b1、a1分别为有理分式分子和分母多项式的系数向量。,例8.9 一个滤波器的差分方程如下:求它的并联型结构。解:由差分方程可以得到该滤波器的系统函数,程序运行的结果为:r=0.5000-0.5669i 0.5000+0.5669i 2.0000 p=-0.2500+0.6614I-0.2500-0.6614i 0.3333 k=b1=1.0000 1.0000 0a1=1.0000 0.5000 0.5000,由b1、a1、r(3)和p(3)可以写出并联形式的系统函数为:,8.5.2 FIR 数字滤波器的结构实现 8.5.2.1 直接型FIR直接型结构由行向量b来描述,b即系数。该结构也由函数filter来实现,只需要将矢量a的位置设为1,即调用方式为filter(b,1,x)。线性相位FIR数字滤波器的结构在本质上还是直接型,只是减少了乘法的运算量,因此在Matlab实现上可以与直接型结构相同。,8.5.2.2 级联型FIR滤波器的级联型结构可以继续使用前面在IIR 级联型中描述的 Matlab 函数z,p,k=tf2zp(b,a)和 sos=zp2sos(z,p,k),只需要把分母矢量 a 置为 1。,但这里应该注意,因为函数 tf2zp(b,a)要求分母 a 的长度要大于等于分子b的长度,且第一个系数不能为 0,因此,如果分子 b 的长度为 N 的话,则分母 a 需要在 1 后面补上 N-1 个零,即 a=1,0,0,0。类似地,用函数 sos2tf 可以从级联型得到直接型的表达式。FIR滤波器系统函数级联的各个二阶节的系数矩阵sos为:,例8.10 已知FIR滤波器的冲激响应为试用其级联型结构实现。解:由题目可以写出该滤波器的系统函数为:,程序运行结果为:p=0 0 0 0 k=1,sos=1.0000 0.1000 0.4000 1.0000 0 0 1.0000 0.2000 0.3000 1.0000 0 0所以级联型结构为:,8.5.2.3 频率抽样型结构 在8.3.3节中最后所得到的FIR滤波器的频率抽样型结构为:当N为偶数:当N为奇数:,在Matlab中,用函数C,B,A=fir_freq_structure(h)来实现FIR滤波器的频率抽样型结构。但是,这个函数没有考虑零极点位置的调整,也就是说,并没有将级联网络的前面的子网络的零点及后面的并联网络的极点从单位园上移到单位园内。因此,在这里,%FIR频率抽样结构实现functionC,B,A=fir_freq_structure(h)%h是FIR滤波器的单位抽样响应向量%C是频率抽样结构所包含的并联部分增益的行向量%B是按照行排列的分子系数向量%A是按照行排列的分母系数向量,Len=length(h);%得到单位抽样响应向量的长度H=fft(h,Len)%求出单位抽样响应的DFTMag_H=abs(H);%得到H(k)的幅度Phase_H=angle(H);%得到H(k)的辐角%判断Len的奇偶if(mod(Len,2)=0),L=Len/2-1;%如果Len为偶数,则L=Len/2-1%如果Len为偶数,则包含()项,所以A1、C1赋值如下A1=1,-1,0;1,1,0;C1=real(H(1),real(H(L+2);,else L=(Len-1)/2;%如果Len为奇数,则L=(Len-1)/2%如果Len为奇数,则仅包含 项,所以A1、C1赋值如下 A1=1,-1,0;%奇数,不包含H(N/2)项 C1=real(H(1);%奇数,不包含H(N/2)项,endk=1:L;B=zeros(L,2);%初始化B,B为长度为2的行向量A=ones(L,3);%初始化A,A为长度为3的行向量,A(1:L,2)=-2*cos(2*pi*k/Len)%A=A;A1;B(1:L,1)=cos(Phase_H(2:L+1);%B1=B(1:L,2)=-cos(Phase_H(2:L+1)-(2*pi*k/Len);%B2=-C=2*Mag_H(2:L+1),C1;%,例8.11 已知FIR滤波器的的冲激响应为 设抽样点数 N=5,求出其频率抽样结构。解:h=1,-1,0,0,1;C,B,A=fir_freq_structure(h);,程序运行结果为:C=4.2979 3.0867 1.0000 B=0.4653-0.9856 0.6479 0.0765 A=1.0000-0.6180 1.0000 1.0000 1.6180 1.0000 1.0000-1.0000 0于是有:,第9章 数字信号处理中的有限字长效应9.1 概述 在这部分之前,我们讨论的实际上是抽样数据信号的处理,因为前面的内容,都只是涉及信号在时间上是离散的这一特征,并没有涉及数值上离散的特征。而对于真正的数字信号的处理,只需要在前面所讨论的离散时间信号处理的原理和方法的基础上,加入字长效应的影响。,9.1.1 数字系统与有限字长效应 一个线性、时不变、因果系统的差分方程为:,到目前为止,在讨论问题时,对于系统的输入序列x(n)、输出序列y(n),对于方程中的系数ai、bi等等,都认为它们的数值有无限精度。但是当具体实现一个离散系统时,无论用软件方式还是硬件方式,实际处理的数都是有限字长的二进制数。当认为一个离散系统的数据是无限精度的,这样的系统就叫做抽样数据系统;若离散系统的数据是有限字长的,则此系统就是数字系统。,对于一个数字系统,由于本应为无限精度的数据变为有限字长来进行处理,因此肯定会对系统的特性产生一定影响,这就是有限字长效应问题。这个问题本来是数字信号处理中的一个重要问题,但是,随着计算机和微处理器技术的飞速发展,运算速度和运算精度都在不断提高,使得有限字长效应的重要性已逐渐降低。不过,在数字信号处理的一些实际应用中,这个问题还是存在的,因此,有必要了解它的影响以及降低影响的一些方法。,9.1.2 关于数的表示 进行数字信号处理时,数的表示有定点制和浮点制两种。浮点制运算比定点制运算的动态范围大,处理精度高,但实现较复杂而且运算速度较慢,因而常用于计算机上的软件实现,进行非实时处理。,在实时处理中定点制运算得到广泛应用,因为它运算速度较快而且硬件实现较经济,但是由于定点运算的动态范围和处理精度受限制较大,因而有限字长效应问题比较突出。本章主要讨论定点制算法的有限字长效应。,9.1.3 量化误差 数的定点表示:设寄存器长L+1位,则除了一位符号位外,可表示的最小数为q=2-L,这个值称为量化间距。若要处理的数有M+1位(含符号位),且ML,则这个数要存储于寄存器中就必须被量化。有两种量化方法:截尾和舍入。截尾就是将寄存器容纳不下的低位数截断;舍入是在数据的L+1位上加1,然后截断为L位。,当数x被量化时,就引入误差e,有:(9.1)其中Qx为x的量化值,即经截尾或者舍入后的值。,由图9.1可知,定点制截尾处理的量化误差et的范围为:补码:-q0时,-qet 0 当 x0时,0etq由图9.2可知定点制舍入处理的量化误差er的范围为:-q/2erq/2,图 9.1 定点制截尾处理的量化特性,图 9.2 定点制舍入处理的量化特性,9.2 A/D变换的字长效应所谓A/D变换即由模拟到数字的变换,一般可分为两步,即抽样与量化编码。抽样数据信号x(n)=xa(nTs)的每个抽样值的精度是无限的,经过量化编码之后,成为有限精度的数字信号。这里不考虑如何进行编码的问题,只讨论其量化效应问题。,9.2.1 量化效应的统计分析A/D变换的结果一般都用定点制补码来表示。量化方法无论采取截尾还是舍入,其误差都可以表示为:e=Qx-x。因此,量化后的抽样值可以表示为:x(n)=Qx(n)=x(n)+e(n)(9.2),(请将图中 改为x),图 9.3 A/D 变换的模型,为了对此模型进行统计分析,要对量化误差序列e(n)作如下假设:1)e(n)是一个平稳随机序列;2)e(n)与信号x(n)不相关;3)e(n)本身样值间不相关,即为白噪声过程;4)e(n)具有等概率密度分布(在一定的量化间距上)。,也就是说,我们将e(n)作为量化噪声,它是白噪声,量化后的信号可以等效为无限精度信号与一噪声相叠加。量化噪声的均值和方差:舍入时:均值 补码截尾时:均值,图 9.4 量化噪声的概率密度函数,而信号功率与噪声功率之比即信噪比为:用对数表示:(9.4)因此,寄存器长度每增加一位(L加上1),信噪比约提高6db。,以上所假设的统计模型对于实际遇到的大多数信号都适合,但若输入信号xa(t)为规则信号,如直流或周期性方波等等,上述假设就不成立。,9.2.2 线性时不变系统对量化噪声的响应当已量化的信号通过一LTI系统H(z)时,由于实际的输入信号如(9.2)式所示,故输出信号为:y(n)=y(n)+f(n)(9.5)其中y(n)是此线性系统对无限精度信号x(n)的响应,f(n)是系统对量化噪声e(n)的响应,故f(n)为输出噪声。输出噪声的功率为:(9.6),这个积分可以用留数定理来计算,其中积分围线C是在H(z)与H(z-1)的公共收敛域内的一条围绕原点的闭合曲线。如果H(z)是稳定系统,则可选单位园为围线C,将z=ej 代入(9.6)式,可以得到:(9.7),9.3 乘积误差的影响 在数字网络中,典型的乘法运算可以表示为:y(n)=a x(n)这里x(n)为数据值,a为乘法器系数。一般来说,每次相乘后要对乘积作舍入或截尾处理。本节讨论在定点实现的相乘运算中乘积的舍入量化误差问题。,图 9.5 相乘运算的统计模型,相乘后的实际结果为:y(n)=Qy(n)=y(n)+e(n)=ax(n)+e(n)(9.8)由9.1节,舍入误差范围为:(9.9)e(n)的统计特性可利用9.2节的假设,故其均值为0,方差为:(9.10),9.3.1 IIR滤波器中乘积误差的影响 首先来分析一阶IIR滤波器:(9.11)其中含有乘积项ay(n-1)。可以将与系数a相乘后乘积的舍入误差所产生的影响等效为存在噪声源e(n),如图9.6所示。,(请将图中 改为y),图 9.6 一阶IIR 滤波器的统计模型,图 9.7 量化噪声通过一阶 IIR 系统,滤波器实际的输出为:y(n)=y(n)+f(n)(9.12)根据9.2.2节,可以求出输出噪声f(n)的方差(功率)为:(9.13),这里e2 如式(9.10)所示。H(z)为系统的传递函数,有:(9.13)式中的(9.14)可用留数定理计算,被积函数为:,选单位园为积分围线C,故被积函数在C内只有一个极点,即z=a,于是有:代入(9.13)式可得:,9.3.1.1 乘积项的有限字长效应与IIR滤波器的结构的关系 1.直接型(正准型)结构 例9.1 已知一个LTI系统如下,用直接型结构实现,求乘积的舍人误差所产生的输出噪声。(9.15),解:信号流图如图9.8(a)所示,图中实线表示H(z)的直接型(正准型)网络结构,e0(n)、e1(n)、e2(n)分别表示与系数0.04、1.7、-0.72相乘后的舍入噪声,它们都经过相同的传输网络H0(z)=1/A(z),A(z)为(9.15)式中的分母多项式,这个等效的噪声网络H0(z)如图9.8(b)所示。因此输出噪声f(n)的功率为:(9.16),被积函数为:(9.17)围线C选单位圆,B(z)在C内有两个极点:z=0.9,z=0.8,故(9.18)因此 f2=(q2/4)90=22.5q2(9.19),2 级联型结构例9.2 H(z)仍如(9.15)式所示,但是用级联型结构实现,求输出噪声功率。解:将H(z)表示为:(9.20),图 9.9 级联型结构的相乘误差,图9.9中e0(n)、e1(n)、e2(n)分别表示与系数0.04、0.9、0.8相乘后的舍入噪声。由图可知,e0(n)和e1(n)通过网络,而e2(n)只通过网络。故输出噪声功率为:,(9.21)可以算得:I1=90 I2=2.78 因此(9.23),2.并联型结构例9.3 将(9.15)式的H(z)用并联型结构实现,求输出噪声功率。解:将H(z)分解为部分分式之和:(9.24)图9.10中e0(n)、e1(n)、e2(n)、e3(n)分别为与系数0.36、0.9、-0.32、0.8相乘后的舍入噪声。,图 9.10 并联型结构的相乘误差,由图可知,e0(n)和e1(n)只通过网络,e2(n)和e3(n)只通过网络,