FIRFIR数字滤波器的原理及数字滤波器的原理及设计设计令输入信号x(n) =(n)，代入上式，有： (7.2)于是得到： 又由(7.2)式可知，当nN-1时，h(n)=0, 即这个系统的冲激响应h(n) 是有限长度的。将ai =h(i)(i=0,1,N-1) 代入(7.1)式得到： (7.3)将(7.3)式的两边进行z变换后，可以得到FIR滤波器的系统函数： (7.4)又由(7.4)式有： n因此，FIR滤波器的系统函数H(z)的极点都位于z=0处，为N-1阶极点；而N-1个零点由冲激响应h(n)决定，一般来说，可以位于有限 z 平面的任何位置。n由于FIR数字滤波器的极点都集中在单位园内的原点z=0处，与系数h(n)无关，因此FIR滤波器总是稳定的，这是FIR数字系统的一大优点。nFIR数字滤波器的频率响应为： (7.5)n所谓线性相位滤波器，就是说此滤波器的相位特性，或者说其频率响应H(ej) 的幅角()，是频率的线性函数。7.2 7.2 线性相位线性相位FIRFIR滤波器滤波器 7.2.1 7.2.1 恒延时滤波恒延时滤波n数字滤波器的相延时为 (7.6)n数字滤波器的群延时为 (7.7)n所谓恒延时滤波就是要求p() 或g() 是不随变化的常量。要使p() 与g()都是不随变化的常量，() 的图象必定是一条过原点的直线，即有： ()= -, 为一常数 (7.8)7.2.2 7.2.2 线性相位线性相位FIRFIR滤波器满足的条件滤波器满足的条件7.2.2.1 7.2.2.1 要求恒相延时与恒群延时同时成立要求恒相延时与恒群延时同时成立() 0图图7.1 时的图象因为 故有： （7.9）由 (7.8) 式和 (7.9) 式有： 利用三角公式，由上式可以得到： （7.10）可以证明，当满足： (7.11)以及 0nN-1 (7.12)时，(7.10) 式成立。这就是说，如果 (7.11) 式和 (7.12) 式满足，便有：()=-，是的线性函数，而且有 ，即恒相延时与恒群延时同时成立。n(7.12)式说明冲激响应h(n)关于中心点偶对称，无论N为偶数还是奇数，对称中心都位于 。n若只要求群延时g() 为一常数，则相位特性是一条可以不经过原点的直线，即： （7.13）n并且有0=/2（这在下面会给予解释），即有 (7.14)7.2.2.2 7.2.2.2 只要求恒群延时成立只要求恒群延时成立0()图图7.3 时的图象时的图象由 (7.9) 式和 (7.14) 式可得： 利用三角公式，由上式可以得到： （7.15）可以证明，当满足: (7.16)以及 0nN-1 (7.17) 时，(7.15) 式成立。这就是说，如果 (7.16) 式和 (7.17) 式满足，便有 , 是的线性函数，而且有g()=，即恒群延时成立。n(7.17)式说明冲激响应h(n) 关于中心点奇对称，无论N为偶数还是奇数，对称中心都位于 。当N为奇数时有 。 n总的来说，当FIR滤波器的冲激响应h(n) 偶对称或者奇对称时，此滤波器的相位特性是线性的，而且群延时是恒定的，为 = 。 7.2.3 线性相位线性相位FIR滤波器的特性滤波器的特性n由冲激响应h(n)为偶对称或者奇对称的对称条件，可以导出线性相位FIR数字滤波器的一些特性。 7.2.3.1 7.2.3.1 网络结构网络结构n根据h(n)的对称性可以简化FIR滤波器的网络结构，详见下面8.3节。 7.2.3.2 7.2.3.2 频率响应频率响应nFIR滤波器的频率响应为:(7.18)n如果FIR滤波器是线性相位的，那末h(n)具有对称性，由此可以导出线性相位FIR数字滤波器频率响应的特有形式。 1. 1. 偶对称偶对称 ， N N为奇数为奇数n此时有h(n)=h(N-1-n)。对(7.18)式分段求和，得到：n令 ，则上式为：n其中： (7.20)n(7.19)式中求和号部分为实数，故H(ej) 的相位为 2. 2. 偶对称偶对称 ， N N为偶数为偶数n此时有h(n)=h(N-1-n)。对(7.18)式分段求和，得到：n于是得到： (7.21)其中： (7.22) n3 3奇对称奇对称 ， N N为奇数为奇数n此时有h(n)= -h(N-1-n)。对(7.18)式分段求和，得到：n于是得到： (7.23)n其中： (7.24)4 4奇对称奇对称 ， N N为偶数为偶数n此时有 h(n)= -h(N-1-n)。将(7.18)式分段求和，得到：n于是得到： (7.25)n其中： (7.26)n 上述四种情况有一个统一的形式，即： (7.27) n其中，H() 是的实函数，是三角函数的线性组合；因此H(ej)的相位由() 决定，而() 是的线性函数。当h(n)偶对称时， ；当h(n) 奇对称时， 。n现在可以解释为什么7.2.2.2节中的0只能够取 /2了。从上面讨论的第3、4种情况我们看到，只要h(n)是奇对 称 的 ， 所 推 导 出 的 频 率 响 应 的 表 达 式 (7.27)中 ， 必 然 有 ；另外，(7.27)式中的H()可能为负数，也就是与模值可能相差 -1=e-j，因为（/2）-=-/2，所以0 也可能为 /2。 就是说，0 只能取/2。n另外，幅度函数H()是三角函数的线性组合，在四种情况下各有不同的形式，但是，并不是每一种形式都能够用于低通、高通、带通、带阻等各种类型的滤波器。例如，在第4种情况下， ，由于是正弦函数的线性组合，故显然当=0时有H()=0，也就是说，=0不可以在相应的滤波器的通带，因此，这种形式不能够用于低通和带阻滤波器。 7.2.3.3 7.2.3.3 零点分布零点分布n如果FIR滤波器是线性相位的，则其N-1个零点在z平面上的分布是有一定的规律的。n对一线性相位FIR滤波器有： ， 0nN-1 因此有： n令m=N-1-n，则： n也即 (7.28)n因此，如果z=zi 是H(z)的零点，那末zi-1 也是H(z)的零点；此外，由于h(n)为实序列，故zi* 也是H(z)的零点，由此又得出(zi*)-1 也是零点。这四个零点构成了互为倒数、互为复共轭对的四点组。n几种特殊情况：若 ，则零点为单位圆上的复共轭对；若zi是不为0的实数，则零点为实轴上的倒数对；若zi= 1或zi= -1，则零点为单点。7.3 7.3 窗口法窗口法n从本节开始，讨论设计FIR数字滤波器的一些主要方法。注意，FIR数字滤波器不能够借助于模拟滤波器的设计方法来设计，而是直接逼近所要求的频率响应。窗口法是设计FIR滤波器的重要的基本方法。7.3.1 7.3.1 窗口法的基本思想窗口法的基本思想n用傅里叶反变换可以求得图7.6所示的理想低通滤波器的冲激响应： （7.29）n将此无限长的hd(n)截断就得到有限长的冲激响应h(n)： (7.30)n这里已经假设了h(n)的长度N为奇数。(7.30)式等价于： (7.31)其中 (7.32)7.3.2 7.3.2 理论分析理论分析n h(n) 是两个序列的乘积，故H(ej) 是这两个序列的傅里叶变换的卷积，即： (7.33)n其中Hd(ej) 是hd(n) 的傅里叶变换，也即理想的频响，而WR(ej) 是矩形窗wR(n) 的频谱。(7.34)WR(ej)是的偶函数。图图7.7矩形窗的频谱矩形窗的频谱n由 (7.33) 式有： （7.35）n式中积分等于由 -c到c区间曲线WRej(-)下的面积，如图7.8中阴影所示。当主瓣的中心变化时，此曲线左右移动，此面积也就发生变化。图 7.8 由 -c 到c 区间曲线WRej(-) 下的面积n当 =0时： n当逐渐增大，随着图中不同正负、不同大小的旁瓣移出和移入积分区间，使得H(ej) 的大小产生波动。n当 ： H(ej) 取最大值，约为1.0895 H(ej0)，此处称为上臂峰，或正肩峰。n当=c： 主瓣的中心移到了c处， 。n当 ：H(ej) 取最小值，约为 -0.0895 H(ej0)，此处称为下臂峰，或负肩峰。n继续增大到，H(ej) 随着积分区间内旁瓣的移动而在阻带内波动。n图7.9表示了H(ej)在由0到范围内变化的情况，图中已经假定H(ej0)=1。在 与 之间为过渡带。图图7.9加矩形窗后的频响与理想频响的比较加矩形窗后的频响与理想频响的比较n因此，加窗后所得到的FIR滤波器的频响H(ej) 出现了过渡带、肩峰以及通带和阻带内的波动。与这些特征有关的因素：1. 1. 过过渡渡带带：正、负肩峰之间的过渡带的宽度等于窗函数频谱的主瓣宽度。对于矩形窗频谱WR(ej)，此宽度为4/N。因此，过渡带宽度与所选窗函数有关；而对于一定的窗函数，增加窗口长度N可以使过渡带变陡。 2 2肩肩峰峰及及波波动动：是由窗函数频谱的旁瓣引起的。旁瓣越多，波动就越快；旁瓣相对值越大，波动就越厉害，肩峰也越强。因此，肩峰及波动与所选窗函数有关。长度N的增加能够使频响的波动加快，但是不能够改变肩峰和波动的相对大小。n因此，加窗法设计FIR滤波器，h(n) 之长度也即窗口长度N可以影响过渡带的宽度；而所选窗函数不仅可以影响过渡带的宽度，还能影响肩峰和波动的大小。选择窗函数应使其频谱：a）主瓣宽度尽量小，以使滤波器的过渡带尽量陡；b）旁瓣相对于主瓣越小越好，这样可使滤波器频响的肩峰和波动减小。n然而，一般情况是，若选择的窗口频谱旁瓣较小，其主瓣就会较宽，反之亦然。因此，常常要根据需要进行折衷的选择。 7.3.3 7.3.3 几种常用窗函数几种常用窗函数n这里介绍几种常用窗函数，它们的长度均设为N，并且都是因果窗，即定义在0nN-1区间，N可以是奇数或偶数，但w(n)都是偶对称的。由7.2.3.2节可以知道，w(n)的频谱可表示为： （7.36）n其中幅度函数W() 是的实函数，而相位函数： (7.37)n因此，对每种窗，只需考察其w(n) 和W() 的表示式。 7.3.3.1 7.3.3.1 矩形窗矩形窗n(7.32)式已给出以n=0为对称中心的非因果矩形窗wR(n)的定义，下面是实际使用的因果矩形窗： (7.38)n这个矩形窗的对称中心在(N-1)/2，它是wR(n)向右移位(N-1)/2的结果，即有： (7.39)n这两个窗函数的频谱之间的关系： (7.40) n如果序列以n=0为对称中心，则其延时= 0，同时相位函数()=0。因此，对称中心也即延时为0的非因果矩形窗的频谱是： (7.43)n而对称中心也即延时为 的因果矩形窗的频谱是： (7.44)7.3.3.2 7.3.3.2 升余弦窗升余弦窗-汉宁汉宁( (Hanning)Hanning)窗窗 (7.45) (N1) (7.46) 图 7.10 升余弦窗的频谱n从图7.10中可以看到，由于三部分频谱的相加，使总的频谱W()的旁瓣大大抵消，从而使能量有效地集中在主瓣内，但其代价是使主瓣与矩形窗相比加宽了一倍。 7.3.3.3 7.3.3.3 改进的升余弦窗改进的升余弦窗哈明哈明( (Hamming)Hamming)窗窗n对升余弦窗作一点调整，得到： （7.47） (7.48)n结果是，99.96% 的能量集中在主瓣内，而主瓣宽度仍与汉宁窗相同。n显然，汉宁窗和哈明窗可以统一表示为： (7.49)n对于汉宁窗， =0.5；对于哈明窗， =0.54。 7.3.3.4 7.3.3.4 二阶升余弦窗二阶升余弦窗布莱克曼布莱克曼( (Blackman)Blackman)窗窗n如果还要进一步抑制旁瓣，可以对升余弦窗再加一个二次谐波的余弦分量： n而 (7.51)（7.50）n以上几种窗函数的时域图像，除矩形窗之外，都是在对称中心 处取最大值，并且向两边逐渐减小，而在n=0和n=N-1处取最小值。 7.3.3.5 7.3.3.5 凯塞凯塞( (Kaiser)Kaiser)窗窗n这种窗通过调整一个参数来进行主瓣宽度与旁瓣衰减之间比重的选择。n凯塞窗利用零阶贝塞尔函数I0(x)构成： 0nN-1 (7.52)nw(n)从对称中心 向两边逐渐减小，越大，减小越快。n 参数值越大，其频谱的旁瓣越小，但主瓣宽度也相应增加，因而改变值就可以在主瓣宽度与旁瓣衰减之间进行选择。图 7.11 Kaiser 窗函数随参数而变化n以上所讨论的五种窗函数的主要性能归纳于表7.1中，可以看出，从上到下，旁瓣衰减越来越厉害，而主瓣却越来越宽。不过，哈明窗与汉宁窗的主瓣宽度相同，但是旁瓣衰减特性却好于汉宁窗。在实际应用中，经常采用哈明窗。表7.2则表示凯塞窗在不同值下的性能。7.3.4 7.3.4 设计方法小结设计方法小结n上面关于窗口法的基本思想以及理论分析虽然是关于非因果矩形窗的，但其基本原则和得出的那些结论对于采用其它窗、采用因果窗也完全适合。在实际应用中，滤波器的冲激响应h(n)都应该是因果序列；所加的窗当然也是因果窗，并且窗口长度N既可以是奇数，也可以是偶数。n线性相位FIR系统的冲激响应h(n)是偶对称或者奇对称的，其频率响应H(ej)是h(n)的傅里叶变换，并且可以表示为： H(ej)=H()ej() (7.53)n不管h(n)是有限长还是无限长，()=-只与h(n)的对称中心有关； 而h(n)是否是无限长的，只对应着H()是否是理想的矩形函数。n可以得到用窗口法设计线性相位因果FIR数字滤波器的步骤：1对给定的理想的频率响应（比如图7.6所示）进行傅里叶反变换，得到的是一个以n=0为对称中心的无限长序列，用hd(n)来表示。2如果所要求的滤波器冲激响应的长度为N(可以是奇数或者偶数)，则将hd(n)向右移位 ，于是得到一个以 为对称中心的无限长序列，用h(n)来表示，即有h(n)=hd(n- )。3用所选定的窗函数和所要求的长度N，对h(n)加窗，即有：h(n)=h(n)w(n)。由于h(n)和w(n)的对称中心都在 ，于是就得到了长度为N、对称中心在 的序列h(n)。 图7.13 非因果理想滤波器冲激响应的移位n有时，理想滤波器的频率响应是这样给出的： （7.54）n对这一Hd(ej)直接进行傅里叶反变换，就得到了以为对称中心的无限长序列h(n)。以低通滤波器为例，有 n这里 = = 。这样的做法相当于将上面的第1步和第2步合并为一步。例例7.17.1 用矩形窗设计一个线性相位因果FIR滤波器，它在0区间的理想频率响应为：n n例例7.17.1 用矩形窗设计一个线性相位因果FIR滤波器，它在0区间的理想频率响应为： n设所求的FIR滤波器的冲激响应h(n)长度为N，则 。 而 nh(n)还可以按照下面的方法求出：n将hd(n)向右移位 ，就得到h(n)： ，与上面得到的结果相同。这里要特别注意，虽然sin(n)=0，但是一定要保留hd(n)中的sin(n)，这样，移位后才能得到正确的h(n)。 n 例例7.27.2 一个带通FIR数字滤波器的衰减指标如图7.14所示，用窗口法设计这个线性相位因果滤波器。图7.14 一个带通滤波器的衰减指标解解：要求阻带最小衰减为45dB, 查表7.1, 可知应该选哈明窗。 n过渡带宽度：1 =p1 s1 =0.13, 2 =s2 p2 =0.15, 应选较窄的，故 =0.13。n由表7.1, 可知对于哈明窗，=8/N, 于是得到: N=8/()=8/0.13=61.54, 于是取 N=62。 n求截止频率：n先用傅里叶反变换求得以0为对称中心的无限长序列：n将hd(n)移位： -nn对h(n)加哈明窗，就得到要求的线性相位因果带通FIR数字滤波器的冲激响应： h(n)=h(n)w(n),其中： 7.6Maltab方法方法用用MatlabMatlab进行基于窗函数进行基于窗函数的的FIRFIR数字滤波器的设计数字滤波器的设计7.6.17.6.17.6.1.1 常用的窗函数常用的窗函数nMatlab中提供了很多常用的窗函数，其中的一些窗函数的调用形式为：矩形窗： w = boxcar ( N )汉宁窗： w = hanning ( N )哈明窗： w = hamming ( N )布莱克曼窗： w = blackman ( N )凯塞（Kaiser）窗： w = kaiser (N, beta) n其中输入参数N表示窗口的长度，返回的变量w是一个长度为N的列向量，表示窗函数在这N点的取值。Beta是控制Kaiser窗形状的参数。nMatlab还提供了函数 M,Wn,beta,ftype=kaiserord(f,a,dev)，以利用Kaiser窗来估计滤波器阶数M、参数Beta、截止频率Wn和需选用的滤波器类型ftype。如果ftype为空，表示滤波器为低通；如果ftype=high，滤波器为高通；ftype=stop，滤波器为带阻；ftype=DC-0，则表示多带滤波器第一个频带为阻带；ftype=DC-1，表示多带滤波器第一个频带为通带。输入参数f表示需设计的FIR滤波器的频带；a表示FIR滤波器在f定义的频带中的幅度值，一般对通带取值为1，阻带取值为0；dev表示FIR滤波器在f定义的频带内的波动值。7.6.1.2 窗函数法设计窗函数法设计FIR数字滤波器数字滤波器1. 1. b=fir1(N,Wn) b=fir1(N,Wn) nfir1函数用来设计 FIR 滤波器，其中N为滤波器的阶数，因此h(n)的长度为N+1；Wn是截止频率，其取值在01之间，这是以抽样频率为基准频率的标称值，故1对应抽样频率；b对应设计好的滤波器的系数h(n)。若Wn是一标量，则用来设计低通滤波器；如果Wn是12的向量，则可以用来设计带通滤波器；如果Wn是1L的向量，则可以用来设计L通带滤波器，注意这时的调用方式为： b=fir1(N,Wn,DC-1) 或 b=fir1(N,Wn,DC-0)n前者保证第一个带为通带，后者保证第一个带为阻带。nfir1 fir1 函数还有以下多种调用形式：函数还有以下多种调用形式： b = fir1(N,Wn,ftype) 当ftype 中的输入参数为 high 字串，即可用来设计高通滤波器；当ftype 中的输入参数为 stop 字串，即可用来设计带阻滤波器，此时Wn是12的向量。 b = fir1(N,Wn,window) b = fir1(N,Wn,ftype,window) 参量 window 表示设计 FIR 滤波器所采用的窗函数类型，以列向量形式表示。向量window 的长度必须为 N+1。若 window 缺省，则 fir1 默认使用哈明窗。2 2 b = fir2(N,f,m)b = fir2(N,f,m)n该函数采用窗函数法设计具有任意频率响应的FIR数字滤波器。其中 f是频率向量，其值在01（标称值）之间，1 对应抽样频率，其第一个点必须是 0，最后一个点必须是 1，而且频率点必须是递增的。m 是对应于频率点 f 处的期望的幅频响应， f和m的长度必须相等。 如同fir1函数，b = fir2(n,f,m,window)可以根据windows的值来选取不同的窗函数，缺省时自动选用哈明窗。nfir1和fir2函数可以用来设计低通、高通、带通、带阻和通用多带 FIR 滤波器。n例例7.37.3 利用矩形窗和哈明窗设计一个FIR低通滤波器，已知c =0.25，N=10。运行结果如图7.16所示。图7.16 分别加矩形窗和哈明窗后的幅频响应曲线例例7.4 7.4 利用窗函数法完成数字带通滤波器的设计，并画出所设计的滤波器的幅频响应图。滤波器的性能指标如下：n低端阻带边界频率：s1=0.2, 高端阻带边界频率：s2=0.8； 阻带最小衰减：As=60dB。n 低端通带边界频率：p1=0.35, 高端通带边界频率：p2=0.65； 通带最大衰减：Rp=1dB。 图7.17 例7.4中带通滤波器的幅频响应曲线解：解：因为 As=60dB，所以应该选用Blackman窗。程序运行结果如下：M = 75（滤波器阶数）Rp =0.0028（实际通带波动）As =75（最小阻带衰减）n从运行结果来看，75 阶 Blackman 窗的 FIR 数字滤波器的实际阻带衰减为 75dB，通带波动为 0.0028dB，显然满足上面所提的技术要求，其幅频响应曲线如上图所示。 例例7.5 7.5 利用Kaiser窗设计一个满足下列指标的数字高通滤波器，并画出幅频响应曲线图。n通带边界频率p=0.6, 阻带边界频率s=0.45， 阻带波纹=0.03。图7.18 例7.5的高通滤波器的幅频响应曲线7.6.2 用用Matlab进行等波纹进行等波纹FIR滤波器的设计滤波器的设计1、b = remez(N,f,a,w,ftype)b = remez(N,f,a,w,ftype)n它的几种调用形式为：（1 1） b = remez(N,f,a) b = remez(N,f,a) n用来设计一个N阶（h(n)的长度为 N+1）的FIR 数字滤波器。f是频率向量，其单位为，范围为 0 f 1，这些频率是顺序递增的；数组 a 对应f中为各指定频率上理想的幅频响应，f 和 a 的长度必须相等，且为偶数；每个频率带中所用的权函数等于 1，这说明在每个频率带中的容限（i）是相同的。数组 b 返回滤波器的系数（即冲激响应）。（2 2） b = remez(N,f,a,w) b = remez(N,f,a,w) n数组 w 是每个频率带的加权向量，其它参数与（1）中类似。 （3 3） b = remez(N,f,a,ftype)b = remez(N,f,a,ftype) n当 ftype 是 字 符 串 hilbert 或 differentiator 时，这个函数可以用来设计数字希尔伯特变换器或数字微分器。2、 N,fo,ao,w = remezord(f,a,dev) N,fo,ao,w = remezord(f,a,dev) n该函数用来确定当用切比雪夫最佳一致逼近设计FIR 滤波器时所需要的阶数、归一化边界频率、频带幅度以及加权系数。输入参量 f 是频率带边界频率，范围为 0 到 fs/2（抽样频率的一半）；a 是所期望的幅度值，f 的长度是 a 的长度的 2 倍减 2。参数 dev 是每个频率带允许的最大波动或误差。n需要注意的是，在某些情况下，remezord 可能高估或低估滤波器的阶数 N，因此，得到滤波器的系数 b 后，必须检查最小阻带衰减，并与给定的As 比较，如果滤波器不能满足技术指标要求，则应提高滤波器的阶数，例如 N+1、N+2 等。n若设计者事先不能确定要设计的滤波器的阶数，可以调用remezord 函数，根据得到的结果，再利用函数b = remez(n,fo,ao,w) 来设计等波纹逼近的 FIR 滤波器。n例例 7.6 7.6 设计满足下列指标的等波纹线性相位FIR低通滤波器： 图7.19 例7.6中等波纹线性相位FIR低通滤波器的幅频响应曲线程序生成的输出数据为：n从上面的结果可以看到，与预料的一样，滤波器系数满足对称约束条件 h(n)=h(N-n)，h(n)的长度是29（阶数N=28），它是一个线性相位 FIR 滤波器。例例7.7 7.7 设计满足例7.4中指标的等波纹线性相位FIR带通滤波器，即滤波器的性能指标如下：n低端阻带边界频率：s1=0.2, 高端阻带边界频率：s2=0.8； 阻带最小衰减：As=60dB。n低端通带边界频率：p1=0.35, 高端通带边界频率：p2=0.65； 通带最大衰减：Rp=1dB。 n程序生成的输出数据为：其幅频响应曲线如图7.20(a)所示。(a)N=23n从图7.20（a）的幅频响应曲线可以看出，该滤波器阻带衰减没有满足指标。（b）图是将阶数增加5之后的幅频响应曲线，满足了指标要求。 (b)N=28结束语结束语谢谢大家聆听！谢谢大家聆听！118