FIR数字滤波器设计和实现.ppt

第五章FIR 数字滤波器设计和实现,概述：IIR 和 FIR 比较,IIR与FIR性能比较IIR数字滤波器幅频特性较好，但相频特性较差 FIR数字滤波器可以严格线性相位，又可任意幅度特性因果稳定系统可用 FFT 计算(计算两个有限长序列的线性卷积)但阶次比 IIR 滤波器要高得多,概述：IIR 和 FIR 比较,IIR 与 FIR 设计方法比较IIR DF无限冲激响应，H(Z)是 z-1 的有理分式，借助于模拟滤波器设计方法，阶数低（同样性能要求）。其优异的幅频特性是以非线性相位为代价的。缺点：只能设计特定类型的滤波器，不能逼近任意的频响。,FIR DF有限冲激响应，系统函数 H(Z)是 z-1 的多项式，采用直接逼近要求的频率响应。设计灵活性强。缺点：设计方法复杂；延迟大；阶数高。(运算量比较大，因而在实现上需要比较多的运算单元和存储单元)FIR DF 的技术要求通带频率p，阻带频率s 及最大衰减p，最小衰减s很重要的一条是保证 H(z)具有线性相位。,概述：FIR DF 设计方法,FIR 数字滤波器设计 FIR 滤波器的任务给定要求的频率特性，按一定的最佳逼近准则，选定 h(n)及阶数 N。三种设计方法 窗函数加权法 频率采样法 FIR DF 的 CAD-切比雪夫等波纹逼近法,概述：FIR DF 零极点,FIR滤波器的I/O 关系：,FIR 滤波器的系统传递函数：,在 Z 平面上有 N-1 个零点；在原点处有一个（N-1）阶极点，永远稳定。,FIR 系统定义：一个数字滤波器 DF 的输出 y(n)，如果仅取决于有限个过去的输入和现在的输入x(n),x(n-1),.,x(n-N+1)，则称之为 FIR DF。FIR 滤波器的单位冲激响应：,FIR DF 的频率响应为：,FIR 滤波器的最重要特点是能实现线性相位。具有线性相移特性的 FIR 滤波器是 FIR 滤波器中应用最广泛的一种。,Hr()：振幅响应，它是一个取值可正可负的实函数。,()=arg H(ejw)为数字滤波器的相位响应。,概述：FIR DF 频率响应,信号通过线性滤波器时，其幅度和相位可能会发生改变，滤波器幅频特性|H()|和相频特性()可能会随频率的变化而改变。如：输入正弦信号 Acos(n0)则：输出为|H(0)|Acos(n0)，其中相移(0)输出频率和输入频率相同，但幅度和相位都发生了变化输出信号比输入信号滞后的样点数 n(位移)可由下式求得:设：n00,滤波器在数字频率0 处的相位延迟（位移）,由于相位延迟 n 的不同，最终产生了相位失真。,确保不产生相位失真的办法：使不同频率的信号通过滤波器时有相同的延迟 n。,概述：相位失真,对不同的频率有恒定的相移，不同的相位延迟 n，会产生相位失真.,如：方波 y(t)可以用无数奇次谐波的正弦波的叠加来得到：,若每个正弦波相移/2 弧度：,确保所有频率具有相同相位延迟的简单方法：随着频率的变化而改变相位，使滤波器具有线性相位特性，即使所有频率的相位延迟保持恒定，这种方法可通过使系统的相位函数()为频率的线性函数来实现。,概述：相位失真,可见相移之后正弦波之和已不再是方波。,线性相移FIR DF 约束条件和频率响应,三个内容 约束条件恒延时滤波h(n)偶对称：恒相延时和恒群延时同时成立 h(n)奇对称：仅恒群延时成立 频率响应Type I：h(n)偶对称、N 为奇数Type II：h(n)偶对称、N 为偶数Type III：h(n)奇对称、N 为奇数Type IV：h(n)奇对称、N 为偶数 FIR DF 零极点分布,相延时：,群延时：,线性相移FIR DF 约束条件：恒延时滤波,恒延时滤波 滤波器的延时分为相延时和群延时两种,令,恒延时滤波器：p()或g()是不随变化的常量，这时滤波器具有线性相位特性。,（负号是因为系统必有时延）,由于 FIR 滤波器的传递函数为:,故:,线性相移FIR DF 约束条件：恒延时,恒相延时和恒群延时同时成立要使p、g 都不随 变化,()必须是一条过原点直线,于是:,线性相移FIR DF 约束条件：恒延时,可以证明，当,线性相移FIR DF 约束条件：恒延时,上式成立，此时,恒相延时和恒群延时同时成立时，线性相位滤波器的必要条件是：不管 N 为偶数，还是 N 为奇数，系统冲激响应 h(n)都关于中心点(N-1)/2 偶对称。当 N 为奇数时对称中心轴位于整数样点上;当 N 为偶数时对称中心轴位于非整数样点上。,于是有：,线性相移FIR DF 约束条件：恒群延时,只要求恒群延时成立 若只要求群延时g()为一常数，则相移特性为不过原点的直线。,故,可以证明，当,上式成立，此时,故,线性相移FIR DF 约束条件：恒群延时,FIR滤波器单独满足恒定群延时的必要条件为：冲激响应 h(n)对中心点(N-1)/2 成奇对称。此时，无论 N 为奇数或偶数，滤波器的相频特性均为线性，并包含有/2 的固定相移：因此，信号通过此类滤波器时不仅产生(N-1)/2 个取样点的延迟，还将产生 90o 的相移，通常这类滤波器又被称为 90o 移相器，并具有很好的应用价值。,当 N 为奇数时,，故,线性相移FIR DF 约束条件：恒群延时,奇对称：()对所有的频率成分都有一个 90相移。,因此，有四种类型的 FIR DF：,线性相移 FIR DF 约束条件,线性相位约束条件对于任意给定的值 N，当 FIR 滤波器的 h(n)相对其中心点(N-1)/2 是对称时，不管是偶对称还是奇对称，此时滤波器的相移特性是线性的，且群延时都是=(N-1)/2。偶对称：()为过原点的，斜率为-的一条直线,线性相移 FIR DF 频率响应：Type I,h(n)偶对称，N 为奇数（恒相时延、恒群时延此时，由于 h(n)序列的长度为奇数，因此滤波器的频率响应函数可进行以下拆分（前后对称部分、中心点）：,对上式的第二和式作变量替换（n=N-1-m)后得到:,由对称条件,则 H(ej)表示为：,线性相移 FIR DF 频率响应：Type I,令,则上式为,由此可以看出其线性相位特性。由于 cos(n)对于=0、2都是偶对称，所以幅度响应 Hr()对=0、2也是偶对称。,线性相移 FIR DF 频率响应：Type I,其中,振幅响应：,相频响应：,N=9,Hr(w),h(n)偶对称，N 为偶数（恒相时延、恒群时延由于h(n)序列的长度为偶数，因此滤波器的频率响应函数可拆分成如下两部分（前后对称部分，中心点处无值）：,线性相移 FIR DF 频率响应：Type II,对上式的第二和式作变量替换（n=N-1-m)后得到:,由对称条件,则 H(ej)表示为：,线性相移 FIR DF 频率响应：Type II,令,，则上式为：,其中,（注意 n 从1 开始，即 b(0)=0，或没有定义）,线性相移 FIR DF 频率响应：Type II,与所设计的 b(n)或 h(n)无关，恒为 0。这种类型（即 h(n)偶对称，N为偶数）不能用于高通或带阻滤波器。（2）由于 cos(n-1/2)对于=是奇对称，所以，Hr(w)对=也是奇对称；以=0、2为偶对称。,振幅响应：,相频响应：,Hr(w),h(n)奇对称，N 为奇数（恒群时延h(n)长度为奇数，拆分成前后两部分：,线性相移 FIR DF 频率响应：Type III,对上式的第二和式作变量替换，并利用对称条件 h(n)=-h(N-1-n)，得:,Hr(w),线性相移 FIR DF 频率响应：Type III,，则上式为：,其中,令,振幅响应：,相频响应：,n 从1开始,与 c(n)或 h(n)的值无关，因此，这种类型的滤波器不适用于低通、带阻或高通滤波器设计，而且，这说明 jHr(w)是纯虚数，对于逼近理想数字希尔伯特变换和微分器，它是很有用的。理想的希尔伯特变换是一个全通滤波器，它对输入信号产生 90 度的相移，它频繁用于通信系统中的调制。微分器广泛用于模拟和数字系统中对信号求导。(2)由于 sin(n)对于=0、2 都是奇对称，所以，Hr(w)以=0、2为奇对称。,注意：,(1)在=0 和 处，有：,线性相移 FIR DF 频率响应：Type III,线性相移 FIR DF 频率响应：Type IV,h(n)奇对称，N 为偶数（恒群时延,其中,线性相移 FIR DF 频率响应：Type IV,Hr(w),与 d(n)或 h(n)的取值无关，因此传输函数 H(z)在 z=1 处为零点。显然，这种类型不能用于实现低通滤波器。又有，所以这类滤波器适用于设计希尔伯特变换和微分器。,注意：,(1)在=0 处，有：,(2)由于 sin(n-1/2)在=处偶对称，在0、2 是奇对称，所以，Hr(w)以=偶对称，0、2为奇对称。,一般形式：,偶对称：,奇对称：,（两个恒时延条件）,（一个恒时延条件）,(Hr()为 的实函数),线性相移 FIR DF 频率响应：小结,四类线性相位FIR滤波器,第一类FIR 系统是 的线性组合，在 时，易取得最大值，因此这一类滤波器易体现低通特性，且是偶函数。通过频率移位，又可体现高通、带通、带阻特性。所以，经典的低通、高通、带通和带阻滤波器的 都是偶对称的。第三、四类FIR 系统是 的线性组合，在 时，的值为零，且是奇函数。这一类滤波器都是作为特殊形式的滤波器，如Hilbert变换器、差分器等。请使用时注意。N最好取为奇数，以便以中心点为对称。,一般的 FIR DF 的零、极点：,在z=0处，有一个（N-1）阶的极点，故滤波器稳定；,其零点要求 f(z)=0，根据代数理论，它为 N-1阶多项式，应有 N-1 个根，所以有 N-1 个零点。如果 h(n)为实数值，其根肯定是共轭对称的。,线性相移 FIR DF 零极点分布,令：m=N-1-n,于是：,线性相移 FIR DF 零极点分布,线性相移 FIR DF 的零极点：,如果 zi 是 H(z)的零点，即 H(zi)=0 则 H(z-1)=0，即 zi-1 亦为 H(z)的零点。,上面提到 zi 肯定是共轭的，故 zi*亦必为其零点于是零点有：,线性相移 FIR DF 零极点分布,总结:,(1)一般情况，,，有四个零点：,(2)r=1，单位圆上的零点：,(共轭对),(3)位于实轴上的实数：b,1/b(实轴上的倒数对)。,(4)zi=1：单零点,例：设FIR滤波器的系统函数为：求出该滤波器的单位取样响应h(n)，判断是否具有线性相位，求出其幅度特性和相位特性。,解：对FIR数字滤波器，其系统函数为：,所以，该FIR滤波器具有第一类线性相位特性，设其频率响应函数为,思路：,理想数字滤波器,设计的 FIR 数字滤波器,要求：,线性相位,尽可能降低逼近误差,FIR DF 窗口法（傅里叶级数法）,hd(n)无限长，且非因果,h(n)有限长，且因果,设所要求的 DF 的频率响应是 Hd(ej)，需要注意：它可能是低通、高通、带通和带阻 FIR DF，没有特指某种类型的数字滤波器。不管是何种 FIR DF，它的频率响应是频域中的周期函数，周期为 2，所以它可以展开为傅氏级数形式：,窗口法：基本原理,式中 hd(n)是傅里叶系数，也是单位取样响应序列。由傅里叶级数理论可得：,因此，所要求的 DF 的系统函数便可求得：显然，Hd(z)是非因果的，且 hd(n)的持续时间为-+，物理上不可实现。我们可以采用逼近 Hd(ej)的方法 首先把 hd(n)先截短为有限项，把 hd(n)截为2M+1项，得：,窗口法：基本原理,然后把截短后的 hd(n)右移，使之变成因果性的序列。令 H(z)等于 H1(z)乘以 z-M 得：令 h(n)=hd(n-M),n=0,1,2,.,2M，则,频率响应 z=ej,窗口法：基本原理,显然H(z)是物理可实现的其冲激响应 h(n)的持续时间也是有限的选择 hd(n)=hd(N-1-n)，保证H(z)具有线性相位。,对 hd(n)的截短必然产生误差，即以|H(ej)|近似|Hd(ej)|。,定义逼近误差为均方误差：,而 Hd(ej)可以展开为：,式中:,窗口法：性能分析,|H(ej)|对|Hd(ej)|的逼近,因为|H(ej)|是对 hd(n)截短而产生的，假定：,即当|n|M 时，An=0,Bn=0。,所以把上述两式代入逼近误差中，利用三角函数的正交性可得：,由于上式中每一项都是正的，所以，只有当,最小。,窗口法：性能分析,说明：,当用|H(ej)|Hd(ej)|时，要使 2=min，|H(ejw)|的截短后的单位取样响应 h(n)的系数必须等于所要求的幅频响应|Hd(ejw)|展成傅里叶级数的系数 hd(n)。,有限项傅氏级数是在最小均方意义上对原信号的最佳逼近,其逼近误差为：,截短的长度 M 越大，逼近误差2 愈小（因为 hd(n)值愈小）。,窗口法：性能分析,将 hd(n)截短：,相当于将 hd(n)与一窗函数 wR(n)相乘，即,窗口法：Gibbs 效应,其中,在一定意义上来看，窗函数决定了我们能够“看到”多少个原来的冲激响应，“窗”这个用词的含义也就在此。,窗函数的频谱,窗口法：Gibbs 效应,此矩形窗谱为一钟形偶函数，在+2/N 之间为其主瓣，主瓣宽度=4/N，在主瓣两侧有无数幅度逐渐减小的旁瓣,见图所示。,截短，根据时域相乘映射为频域卷积，得：,窗口法：Gibbs 效应,为便于分析，我们假定|Hd(ejw)|是理想低通滤波器 LPF。,式中积分等于由 c 到 c 区间内 WNej(w-)下的面积，随着变化，窗函数的主瓣和不同正负、不同大小的旁瓣移入和移出积分区间，使得此面积发生变化，也即|H(ejw)|的大小产生波动。,卷积,窗口法：Gibbs 效应,现在分析几个特殊频率点的滤波器性能=0 时：由于一般情况下都满足 c 2/N，因此，H(0)的值近似等于窗谱函数 WR(ejw)与轴围出的整个面积。,窗口法：Gibbs 效应,=c 时：此时窗谱主瓣一半在积分区间内一半在区间外，因此，窗谱曲线围出的面积，近似为=0 时所围面积的一半，即。,=c-2/N 时，正肩峰 此时窗谱主瓣全部处于积分区间内，而其中一个最大负瓣刚好移出积分区间，这时得到最大值，形成正肩峰。之后，随着值的不断增大，H(ejw)的值迅速减小，此时进入滤波器过渡带。=c+2/N 时，负肩峰 此时窗谱主瓣刚好全部移出积分区间，而其中一个最大负瓣仍全部处于区间内，因此得到最小值，形成负肩峰。之后，随着值的继续增大，H(ejw)的值振荡并不断减小，形成滤波器阻带波动。,理想滤波器的不连续点演化为过渡带,通带与阻带内出现起伏,过渡带：正负肩峰之间的频带。其宽度等于窗口频谱的主瓣宽度。对于矩形窗 WR(ejw),此宽度为 4/N。,肩峰及波动：这是由窗函数的旁瓣引起的。旁瓣越多，波动越快、越多。相对值越大，波动越厉害，肩峰越强。肩峰和波动与所选窗函数的形状有关，要改善阻带的衰减特性只能通过改变窗函数的形状。,加窗处理对理想矩形频率响应的影响,窗口法：Gibbs 效应,Gibbs 现象,在对 hd(n)截短时，由于窗函数的频谱具有旁瓣，这些旁瓣在与 Hd(ejw)卷积时产生了通带内与阻带内的波动，称为吉布斯现象。长度 N 的改变只能改变 坐标的比例及窗函数 WR(ejw)的绝对大小，但不能改变肩峰和波动的相对大小（因为不能改变窗函数主瓣和旁瓣的相对比例，波动是由旁瓣引起的），即增加 N，只能使通、阻带内振荡加快，过渡带减小，但相对振荡幅度却不减小。,结论：过渡带宽度与窗的宽度 N 有关，随之增减而变化。阻带最小衰减（与旁瓣的相对幅度有关）只由窗函数 决定，与 N 无关。,窗口法：Gibbs 效应,Gibbs现象；,窗口法：Gibbs 效应,设计FIR DF时，窗函数不仅可以影响过渡带宽度，还能影响肩峰和波动的大小(阻带的衰减)，因此，选择窗函数应使其频谱：主瓣宽度尽量小，以使过渡带尽量陡。旁瓣相对于主瓣越小越好，这样可使肩峰和波动减小，即能量尽可能集中于主瓣内。对于窗函数，这两个要求是相互矛盾的，要根据需要进行折衷的选择，,窗口法：常用窗函数,为了定量地比较各种窗函数的性能，给出三个频域指标：3db 带宽 B，单位为（最大可能的频率分辨力）最大旁瓣峰值 A(dB)，A 越小，由旁瓣引起的谱失真越小旁瓣谱峰渐进衰减速度 D（dB/oct）一个好的窗口，应该有最小的 B、A 及最大的 D。,以下介绍的窗函数均为偶对称函数，都具有线性相位特性。设窗的宽度为N，窗函数的对称中心点在(N-1)/2处。因此，均为因果函数。矩形窗最简单的窗函数，从阻带衰减的角度看，其性能最差。它的频率响应函数为：,窗口法：基本窗函数_矩形窗,为了对过渡带和阻带衰减进行精确分析，对窗振幅响应进行连续积分（或累积振幅响应），即矩形窗函数 w(n)以及它的振幅响应、累积振幅响应如下图所示。,窗口法：基本窗函数_矩形窗,性能指标3dB 带宽 B=0.89最大旁瓣峰值 A=-13dB旁瓣谱峰渐进衰减速度 D=-6dB/oct,在 Matlab 中，实现矩形窗的函数为 w=boxcar(n)。,振幅响应在=1 处具有第一个零点：因而主瓣的宽度为 2，所以过渡带宽也近似为 2。大约在=3/N 处，出现第一个旁瓣（即主旁瓣），其幅度为：将它与主瓣振幅 N 比较，则最大旁瓣峰值A(dB)为 A=-13db。累积振幅响应第一个旁瓣为 21dB，这个 21dB 的阻带衰减与窗长度 N 无关。根据最小阻带衰减，可以精确地计算出过渡带宽为：它大约是近似带宽的一半。,窗口法：基本窗函数_矩形窗,三角窗（或 巴特利特 Bartlett 窗）由于矩形窗从 0 到 1（或 1 到 0）有一个突变的过渡带，这造成了吉布斯现象。Bartlett 提出了一种逐渐过渡的三角窗形式，它是两个矩形窗的卷积。B=1.28,A=-27dB,D=-12dB/oct，近似过渡带宽 8/N，精确过渡带宽 6.1/N，最小阻带衰减 25dB。与矩形窗来比较，阻带衰减性能有所改善，但代价是过渡带的加宽。,窗口法：基本窗函数_三角窗,在 Matlab 中，函数 bartlett(n)和 triang(n)用来计算相似的三角窗，但它们有两个重要的区别：bartlett 函数返回的序列两端总是 0，因此，对于奇数 n，语句 bartlett(n+2)的中间部分等于 triang(n)；对于偶数 n，bartlett 仍然是两个矩形序列的卷积，但 n 为偶数时的三角窗没有标准定义。,5.3.3 窗口法：基本窗函数_三角窗,余弦窗B=1.2,A=-23dB,D=-12dB/oct。近似过渡带宽 8/N，精确过渡带宽6.5/N，最小阻带衰减 34dB。,窗口法：基本窗函数_余弦窗,或,其中,频率响应,窗口法：基本窗函数_余弦窗,升余弦窗函数汉宁窗、汉明窗、布莱克曼窗都是升余弦窗的特例。它们都是频率为 0 2/(N-1)和 4/(N-1)的余弦序列的组合。升余弦窗的频率特性比矩形窗有很大改善。其中 A、B、C 为常数。当 A=0.5，B=0.5，C=0 时，为汉宁(Hanning)窗。Matlab 中，w=hanning(n)当 A=0.54，B=0.46，C=0 时，为汉明(Hamming)窗。Matlab 中，w=hamming(n)当 A=0.42，B=0.5，C=0.08 时，为布莱克曼窗。Matlab 中，w=blackman(n),窗口法：升余弦窗函数,窗口法：升余弦窗函数,Hanning 窗（升余弦窗）,窗口法：升余弦窗函数_汉宁窗,B=1.44,A=-32db,D=-18db/oct，近似过渡带宽 8/N，精确过渡带宽 6.2/N，最小阻带衰减 44dB。与矩形窗来比，最小阻带衰减性能明显提高，但过渡带也明显增大。,窗口法：升余弦窗函数_汉宁窗,Hamming 窗（改进的升余弦窗）,窗口法：升余弦窗函数_汉明窗,B=1.3,A=-43dB,D=-6dB/oct，近似过渡带宽 8/N，精确过渡带宽6.6/N，最小阻带衰减 53dB。通过这一系数调整，使能量的 99.963%都集中在了窗谱的主瓣内。,窗口法：升余弦窗函数_汉明窗,Blackman 窗（二阶升余弦窗）,窗口法：升余弦窗函数_布莱克曼窗,B=1.68,A=-58db,D=-18db/oct，近似过渡带宽 12/N，精确过渡带宽11/N，最小阻带衰减 74dB。通过增加余弦的二次谐波分量，能够进一步抑制旁瓣，但主瓣宽度却比矩形窗谱的主瓣宽度大三倍。,窗口法：升余弦窗函数_布莱克曼窗,比较以上窗函数，可以看到，矩形窗函数具有最窄的主瓣B，但也有最大的旁瓣峰值 A 和最慢的衰减速度 D。汉宁窗主瓣稍宽，但有着较小的旁瓣和较大的衰减速度，因而被认为是较好的窗口。,窗口法：凯瑟窗函数,凯瑟（Kaiser）窗上面讨论的几种窗函数以牺牲主瓣宽度，换取旁瓣抑制；Kaiser 窗全面反映了这种主瓣和旁瓣衰减之间的互换关系；定义了一组可调的由零阶贝塞尔 Bessel 函数构成的窗函数；通过调整参数可以在主瓣宽度和旁瓣衰减之间自由选择它们的比重。从而实现以同一种窗类型来满足不同窗性能需求的目的。Kaiser 窗函数由 J.F.Kaiser 提出，由下式给出：其中 I0 是修正过的零阶贝塞尔 Bessel 函数是用来调整窗形状的参数，依赖于参数 N。,对于相同的 N，Kaiser 窗可以提供不同的过渡带宽，这是其他窗函数做不到的。通过调整参数，就可以方便地完成对过渡带宽度和阻带衰减的调整。参数越高，其频谱的旁瓣越小，但主瓣宽度也相应增加。,零阶贝塞尔函数,Kaiser 窗函数依参数 而变化,Matlab 中，函数 w=kaiser(n,beta)实现 Kaiser 窗,窗口法：凯瑟窗函数,弗里德里希威廉贝塞尔（Friedrich Wilhelm Bessel）1784年1846 年德国天文学家及数学家,下面图固定，当窗的长度变化时，相应的旁瓣的高度保持不变。,窗口法：凯瑟窗函数,=5.658，则过渡带宽等于 7.8pi/N，最小阻带衰减为 60dB，如下图所示：,窗口法：凯瑟窗函数,凯瑟窗的计算由于 Bessel 函数的复杂性，这种窗的设计公式很难推导，为此，Kaiser 提出了经验公式。给定 p、s、Rp 和 As，参数定义如下：对于过渡带宽=s-p(rad/s)，滤波器阶数为需要强调的是：阶数为 N 的滤波器大致能满足要求，但最后的结果还必须要演算以便证明这一点。在 Matlab 中，函数 w=kaiser(n,beta)实现 Kaiser 窗。,5.3.3 窗口法：凯瑟窗函数,窗口法：常用窗函数的性能指标,FIR DF 窗口法设计步骤 性能要求 Hd(e j)把 Hd(e j)展成傅里叶级数，得到 hd(n)；把 hd(n)自然截短到所需的长度 N=2M+1；将截短后 hd(n)右移 M 个取样间隔，得 h(n)；将 h(n)乘以合适的窗函数，即得所需的滤波器的冲激响应，这时窗函数以 n=M对称（当然窗函数也可直接加在 hd(n)上，这时窗函数以原点为对称）；利用 h(n)，既可用硬件构成滤波器的系统函数 H(z)，也可直接用计算机软件实现滤波。,窗口法：设计步骤,数字低通滤波器的设计例：一个理想低通数字滤波器的频率响应如图所示，为：假定 c=0.25，分别取 N=11、21、31的线性相位 FIR，观察加窗后对滤波器幅频特性的影响。,窗口法：设计步骤_ 数字低通,解:由于 Hd(e j)是一个实周期函数，把它展成为傅氏级数：,窗口法：设计步骤_ 数字低通,式中,将 hd(n)截短为 N=2M+1，并将截短后的 hd(n)移位，得：,窗口法：设计步骤_ 数字低通,然后乘以窗函数 wR(n)，得到最后得 h(n)。,对于 c=0.25，由上式得：,（当 c=，就得到一个全通滤波器）,当 N=11 时，M=5，求得,h(0)=h(10)=-0.045，h(1)=h(9)=0，h(2)=h(8)=0.075，h(3)=h(7)=0.1592，h(4)=h(7)=0.2251，h(5)=0.25。,当 N=11 时，乘以汉明窗：,窗口法：设计步骤_ 数字低通,当 N 取不同值时，H(ej)都不同程度上近似于Hd(ej)。N 过小时，通带过窄,且阻带内波纹较大,当 N 增加时，通频带接近于0.25pi,阻带内波纹减小,但在通带内出现了波纹，随着 N 的继续增加，这些波纹并不能消失。由图中还可以看到，使用汉明窗后，通带内的振荡基本消失，阻带内的纹波也大大减小，从这一点上来说，滤波器的性能得到了改善，但是，这是以过渡带的加宽为代价的。,窗口法：设计步骤_ 数字低通,适用范围对于能用解析式表达，且傅里叶级数的系数容易求解的滤波器：此时，窗函数法是设计 FIR DF 较为方便的一个方法如果 hd(n)不易求，则使用该方法较为困难窗函数用窗口法设计 FIR DF，一个重要问题是选用何种窗函数 w(n)进行截短，以及截短的长度 N。窗函数的选择：阻带衰减指标满足设计给定的阻带衰减和其它滤波器性能要求；能量尽量集中于主瓣内；个人的经验及喜好有关。,窗口法：设计步骤_ 数字低通,窗函数长度 N 的选择：过渡带宽指标采用试验方法，即逐渐增大 M，直至 H(ej)在通带和阻带内部达到指标要求。若对|Hd(ej)|的过渡带提出了具体要求，因为 FIR DF 的过渡带等于窗函数的主瓣宽度，那么通过查表计算 N：例如 5.26，5.86。注意：（1）根据所要设计线性相位 FIR DF类型来决定最终 N 取奇数还是偶数。（2）一般选择 N 为奇数。,5.3.4 窗口法：设计步骤_ 数字低通,窗口法：设计步骤_ 数字低通,截止频率 c 的确定：截止频率 c 对应于明确的 0.5 增益点，而不再标志某个增益点。对于非理想滤波器，其截止频率 c 不采用通带边缘频率 p 或阻带边缘频率 s，而使用过渡带的中点（即通带边缘和阻带边缘之间的中点。因此，窗函数法不能精确确定其通带和阻带的边缘频率）：,例：根据下列指标设计一个线性相位FIR低通滤波器 通带边缘频率 fp=2kHz 阻带边缘频率 fstop=3kHz 阻带衰减 40dB 取样频率 fs=10kHz,解：(1)求对应的理想数字频率：过渡带宽=3kHz 2kHz=1kHz。转换为数字频率过渡带：截止频率：数字截止频率,(2)设理想线性相位滤波器为：由此可得脉冲响应：,(4)由过渡带宽确定窗口长度：,(3)由阻带衰减确定窗函数：因为阻带衰减 40dB，通过查表知道，可以 选择 Hanning 窗：,(N-1)/2=15,则此滤波器的脉冲响应为：,例:用矩形窗、汉宁窗和布莱克曼窗设计FIR低通滤波器，设N=11,c=0.2rad。,解:用理想低通作为逼近滤波器，有,用汉宁窗设计：,用布莱克曼窗设计：,低通幅度特性,例：请选择合适的窗函数及N来设计一个线性相位低通滤 波器 要求其最小阻带衰减为-45dB，过渡带宽为。求出。（设）,例:一段乐曲中夹杂着高频噪声，严重影响收听质量。下图给出了噪声污染后的信号及其频谱。系统的取样频率为 16kHz。设计一个滤波器来提高声音质量。,从图中频谱可以看出，噪声从 2kHz 的频率点开始占据支配地位，因为这个频率很低，消除噪声的同时也会损失一部分音乐信息。可以用通带边缘频率为 2kHz 的低通滤波器对此段音乐进行滤波。由于没有特殊的阻带衰减要求，任何具有合理特性的窗函数即可。在此选择 Hamming 窗。本例中，对过渡带宽也没有特殊要求。为了得到合理的陡峭滚降，选择 N=101，则过渡带宽为,此低通滤波器的参数总结如下：取样频率：16kHz通带边缘频率：2kHz过渡带宽：545Hz截止频率：数字截止频率：窗函数：Hamming 窗，即,窗口长度：N=101。则此滤波器的脉冲响应为：得到的脉冲响应、幅频特性如下图所示。,伴有噪声的音乐经过上述滤波器滤除噪声后的信号及其频谱如下图所示。滤波后的信号几乎没有噪声，但歌曲听起来有点压抑，这是因为歌曲中高频分量也随噪声一起被滤除的缘故。可以想象，如果噪声存在于所有的频率分量上，则不可能在不严重降低信号质量的情况下滤除噪声。,以上设计的是数字低通滤波器，若希望设计数字高通、带通和带阻滤波器，只需要改变付氏级数系数中积分的上、下限即可。数字高通滤波器,窗口法：设计步骤_数字高通,令其时域右移 M 位后的幅频特性为：,（频域为 e-jM，表示时域右移 M 位）,窗口法：设计步骤_数字高通,则,求得,窗口法：设计步骤_数字高通,从这个结果可以看出：一个高通滤波器相当于用一个全通滤波器（即c=）减去一个低通滤波器。传输函数：脉冲响应：,数字带通滤波器,窗口法：设计步骤_数字带通,令其幅频特性为：,（频域为 e-jM，表示时域右移 M 位）,则,窗口法：设计步骤_数字带通,求得,从这个结果可以看出：一个带通滤波器相当于两个低通滤波器相减，其中一个截止频率为 h，另一个为 l。,传输函数：,脉冲函数：,窗口法：设计步骤_数字带通,或者一个带通滤波器相当于一个低通滤波器和一个高通滤波器相乘，即先经过一个 LP DF，再经过一个 HP DF。,传输函数：,脉冲函数：,（频域相乘，时域卷积）,例：用布拉克曼窗设计一个线性相位的理想带通滤波器 求出 的表达式。,解:可求得此滤波器的时域函数为,采用布拉克曼窗设计时(N=51)：,解：,因h(n)为偶对称且N为奇数，其线性相位结构如图表示。,例：为取样频率为 22kHz 的系统设计一个 FIR 带通滤波器，中心频率为 4kHz，通带边缘在 3.5kHz 和 4.5kHz，过渡带宽为 500Hz，阻带衰减 50dB。,数字带阻滤波器,窗口法：设计步骤_数字带阻,令其幅频特性为：,（频域为 e-jM，表示时域右移 M 位）,则,窗口法：设计步骤_数字带阻,求得,从这个结果可以看出：一个带阻滤波器相当于一个低通滤波器加上一个高通滤波器，低通滤波器的截止频率为 l，高通在 h。,传输函数：,脉冲函数：,例：利用窗函数法设计满足下列指标的最低阶数的线性相位 FIR数字带阻滤波器。该数字带阻滤波器的性能指标如下：低端阻带边缘：；低端通带边缘：；高端阻带边缘：；高端通带边缘：；阻带衰减不小于52dB。,窗口法：设计步骤,窗口法：设计步骤,设理想线性相位带阻滤波器的频率响应为：,其中,=(N-1)/2=33(或=26),窗口法：设计步骤,则,汉明窗,FIR DF Matlab 设计函数b=fir1(n,wn,options)，单带 FIR 滤波器b=fir2(n,f,m,options)，多带 FIR 滤波器两者可设计低通、高通、带通、带阻和通用多带 FIR 滤波器Fir1 具有以下多种形式：b=fir1(n,wn)b=fir1(n,wn,ftype)b=fir1(n,wn,window)b=fir1(n,wn,ftype,window)b=fir1(.,noscale)参数向量 b 是 n 阶 FIR 滤波器的系数截止频率 wn 是从 0 到 1 之间的数，1 对应着奈氏频率。对于高通滤波器，ftype 为 high；带阻滤波器 ftype 为 stop,窗口法：Matlab 实现,对于带通或带阻滤波器，wn 为包含通频带边带频率的一个二元素向量 wn1,wn2。参量 window 表示所采用的窗函数类型。window 的长度必须为 n+1（n 为滤波器的阶数），若 window 却省，则 fir1 使用汉明窗。注意因为奇数阶的 II 型滤波器（h(n)为偶对称，长度N为偶数）在高频段的频率响应为零，所以 fir1 函数在高通和带阻情况下不设计 II型滤波器，因此，如果 n 为奇数时，fir1 将阶次加 1 并返回 I 型滤波器。函数 b=fir2(n,f,m)也可设计加窗 FIR 滤波器，但它针对任意形状的分段（piece-wise）线性频率响应。向量 f 由从 0 到 1 的频率点组成，其中 1 表示奈氏频率，第一个点必须是 0，最后一个点必须是 1，频率点必须是递增的。m 是对应于频率点 f 处的期望的频率幅值响应。f 和 m 的长度必须相等。,窗口法：Matlab 实现,Matlab 频率响应函数在 Matlab 中提供了一个 freqz 函数，利用这个函数开发一个新的函数 freqz_m，它给出了绝对的和相对的 dB 值幅度响应、相位响应以及群延时响应。,窗口法：Matlab 实现,function db,mag,pha,grd,w=freqz_m(b,a);%Modified version of freqz subroutine%-%db,mag,pha,grd,w=freqz_m(b,a);%db=Relative magnitude in dB computed over 0 to pi radians%mag=absolute magnitude computed over 0 to pi radians%pha=Phase response in radians over 0 to pi radians%grd=Group delay over 0 to pi radians%w=501 frequency samples between 0 to pi radians%b=numerator polynomial of H(z)(for FIR:b=h)%a=denominator polynomial of H(z)(for FIR:a=1)%H,w=freqz(b,a,1000,whole);H=(H(1:1:501);w=(w(1:1:501);mag=abs(H);db=20*log10(mag+eps)/max(mag);pha=angle(H);grd=grpdelay(b,a,w);,例5.5 根据下列技术指标，设计一个数字 FIR 低通滤波器 wp=0.2,Rp=0.25dB ws=0.3，As=50dB 采用汉明窗，确定脉冲响应，并给出所设计的滤波器的频率响应图。解:在设计中，没有使用通带波动值 Rp=0.25dB，但必须检查设计的实际波动，验证它是否确实在给定容限内。,窗口法：Matlab 实现,%Lowpass filter design-Hamming windowwp=0.2*pi;ws=0.3*pi;tr_width=ws-wpN=ceil(6.6*pi/tr_width)+1 n=0:1:N;wc=(ws+wp)/2h=fir1(N,wc/pi);db,mag,pha,grd,w=freqz_m(h,1);delta_w=2*pi/1000;Rp=-(min(db(1:1:wp/delta_w+1)%Passband RippleAs=-round(max(db(ws/delta_w+1:1:501)%Min Stopband attenuation,运行结果如下：tr_width=0.3142（过渡带宽）N=67（滤波器的阶数，长度为 68，Type II 偶对称偶数）wc=0.7854（理想 LPF截止频率）Rp=0.0364（实际通带波动）As=53（最小阻带衰减）,窗口法：Matlab 实现,从结果看，67 阶 Hamming 窗的 FIR 数字滤波器的实际阻带衰减为 53dB，通带波动为 0.0364dB，显然满足上面所提的技术要求，其时域和频域响应曲线如下图所示。,例5.6 利用例 5.5 给出的设计技术指标，选择 Kaiser 窗，设计所需的低通滤波器。解：,窗口法：Matlab 实现,%Lowpass filter design-Kaiser windowwp=0.2*pi;ws=0.3*pi;As=50;tr_width=ws-wp;N=ceil(As-7.95)/(14.36*tr_width/(2*pi)+1)+1n=0:1:N;beta=0.1102*(As-8.7)wc=(ws+wp)/2;h=fir1(N,wc/pi,Kaiser(N+1,beta);db,mag,pha,grd,w=freqz_m(h,1);delta_w=2