FIR滤波器设计

一、      滤波器设计概述

1.1 数字滤波器设计的基本要求

本章起将着重研究信号分析的逆问题，那就是给定了输入输出指标，设计系统。这种可以人为调整参数来满足给定要求的硬件或软件系统，称为滤波器。这种问题没有唯一的解，因为可能有多种类型的滤波器都能满足给定的要求。要人为地进行选择，就要考虑更多的因素,所以设计比之分析更加接近与工程实际。有两种类型的数字滤波系统。本章将讨论FIR数字滤波器的几种基本设计算法；下一章讨论IIR数字滤波器的设计问题。

在计算技术高度发展的今天，几乎所有的电子设计软件中都把滤波器设计作为其模块之一。那种用查表、手工计算的方法在工程实践中已经濒于淘汰。 

1.2数字滤波器的设计

数字滤波器的设计要经过三个步骤：

①确定指标：在设计一个滤波器前，必须有一些指标。这些指标要根据应用确定。

②模型逼近：一旦确定了技术指标，就可利用前面学习过的基本原理和关系式，提出一个滤波器模型来逼近给定的指标体系。这是滤波器设计所要研究的主要问题。

③实现：上面的两步的结果得到的滤波器，通常是以差分方程、系统函数或脉冲响应来描述的。根据这个描述，可以如第六章所讨论的那样，用硬件或计算机软件来实现它。

1.3数字滤波器的主要指标：

在很多实际应用，像语音或音频信号处理中，数字滤波器用来实现选频功能。因此，指标的形式通常为频域中的幅度和相位响应。幅度指标可以两种方式给出。第一种，叫做绝对指标，它提出了对幅度响应函数|H(jω)|的要求。这些指标一般可直接用于FIR滤波器。第二种方法叫做相对指标，它以分贝(dB)值的形式提出对幅频特性的要求。其值定义为：相对幅度指标：

由于定义中所包含的归一化，滤波器的相对幅频特性最高处的值为0dB，又由于定义中的负号，幅频特性小的地方，其dB值反而是正的。为了说明这些指标，以低通滤波器为例进行讨论。首先看下面两种指标图，它们分别用绝对指标和相对指标表示。

图1.1 FIR滤波器技术指标：（a）绝对指标（b）相对指标

绝对指标和相对指标间的变换关系：

可以用MATLAB语句fgp711准确地求出绝对指标和相对指标的对应数值。并画出如图中所示曲线。

图1.2滤波器绝对和相对指标之间的数据曲线

相位指标主要是线性相位条件。

理想的滤波器相位θ应该和频率ω成正比，也就是θ(ω)= -τω。因此满足此式的系统对信号中所有的频率分量都具有相同的时间延迟τ (无量纲)，因此无波形失真。

略低一点的要求是相位θ和频率ω成线性关系，即θ(ω)= θ0-τω。

IIR滤波器无法做到线性相位，通常用群迟延g(ω)来评价。其定义为相位θ对频率ω的导数，它愈接近常数愈好。

1.4 小结

上述指标是针对低通滤波器的，其他类型的选频滤波器，像高通和带通滤波器，它们的技术指标可作类似的规定，并且可以由低通滤波器演化而得。

本章讨论FIR数字滤波器的设计和逼近。FIR滤波器在设计和实现上具有如下的优越性：

      (a). 相位响应可为严格的线性；

        (b). 不存在稳定性问题；

        (c). 只包含实数算法，不需要递推运算。 

二、      线性相位滤波器的特性

2.1 线性相位滤波器的特性

线性相位条件要求滤波器分子系数bn满足对称性条件。设滤波器的系数长度为N，则这些系数应关于中心点τ=n0=(N-1)/2对称。偶对称时，h(n)=h(N-n-1)；而奇对称时，h(n)=-h(N-n-1)。再考虑到N可以为奇数或偶数，总共有四种类型的线性相位FIR滤波器。

在讨论线性相位滤波器频率响应时需要引进幅特性正负号的概念。以往常设幅特性(Magnitude Response)为正数，因为幅特性的反号可以用相特性加减π来补偿。当相位特性要求线性，不得随便增减时，幅特性就必须分出正负,称为‘符幅特性’(Amplitude Response)。为了计算符幅特性，需要直接作解析计算，求出频率特性的解析式，并将它分解为幅度和相位特性，即，如果θ(ω)满足线性相位条件，这时的A(ω)就是符幅特性。对于长度为N的序列h(n)，先考虑两个对称的系数h(n)和h(N-1-n)，它们的傅立叶变换为（注意运用样本值相等的条件使两项合一）：

2.2  I型线性相位滤波器的特性

类型I，系数对称，N奇数。序列h的下标n从0起算，到N-1为止，中点位置是L=(N-1)/2。h(L)是一个孤项，其他可配成L对系数。将它们的符幅特性加起来，得到：

为了使四种类型滤波器的公式简明统一，引进两个参数。一个是令，它就是(7.1.5)式中的τ值，用以反映延迟的样本数或群延迟，它可以是分数。另一个是。即把τ值向下取整，因为要把L用作下标，它必须是整数。当N为奇数时，τ=L。

在推导中，还注意前一个求和号是从0到L-1,后一个求和号是从0到L,那是把孤项也当成成对项求和，然后再减去一个孤项h(L)。这样做的目的也是使四种类型滤波器的公式统一。

显然，为相角项，相角-ωτ,它与ω成严格的线性关系；A(ω)就是符幅项。

可见，其符幅特性由L 1个余弦项叠加组成。n=0时的符幅项为2h(0)cos(Lω)，在ω=[0~2π]频段间波动L个周期；余弦函数在ω=0，π和2π处都不等于零，因此类型I线性相位滤波器既可以用作低通滤波器(在ω=0处，幅特性不为零)；也可以用作高通滤波器(在ω=π处，幅特性不为零)；且可用作带通和带阻滤波器。

2.3 II型线性相位滤波器的特性

类型II，系数对称，N偶数。如果N为偶数，那末全部系数都可以配对，不会出现中心点的单项，一共有N/2组对称出现的系数，频率特性就成为

此时τ=(N-1)/2将不是整数, 也就是说，对称中心将在两个样本点的中间。其相角特性仍为θ(ω)= -ωτ,它仍与ω成严格的线性关系。

用整数来表示求和号的上限。得出其符幅特性的表达式

可见，它的符幅特性也由L 1个余弦分量构成。这样，不管N为奇数或偶数时，L都是整数，求和的项数都是L 1。由于τ=(N-1)/2是分数，所有余弦分量cos((τ-n)ω) 中都含一个0.5ω。在ω=π处，就会出现符幅特性必定为零的情况。所以类型II不能用作高通和带阻滤波器。

2.4 III型线性相位滤波器的特性

类型III，系数反对称，N奇数。当系数序列反对称时，在(7.2.1)式的推导过程中，h(N-1-n)项将反号而有h(N-1-n) = -h(n)，两个反对称项合成的结果是：

 所以只要把类型I 结果中的cos换成sin，得到：

当N为奇数时，τ=L=整数。相特性为θ(ω)= π/2-ω(N-1)/2=π/2 -ωτ,它虽然仍与ω成线性关系，但多了一个常数项π/2。注意反对称序列的对称中心n=L处，序列的值h(L)=0。所以可以放入求和号中，得到它的符幅特性为：

正弦函数在ω=0，π和2π处都等于零，因此类型III线性相位滤波器的符幅特性在ω=0，π和2π处都等于零。它不能用作低通，也不能用作高通滤波器。故不宜用于选频滤波器。 

2.5 IV型线性相位滤波器的特性

类型IV，系数反对称，N偶数。与类型II滤波器相仿，若N为偶数，则共有N/2组数值相同反对称出现的系数，不出现单项，频率特性式就成为

此处仍有τ=(N-1)/2，因此τ将不是整数,带有小数部分0.5。故对称中心将在两个样本点的中间。它的相角特性为θ(ω)= π/2-ωτ,虽然仍与ω成线性关系，但多了常数项π/2，故不通过原点。

类型IV的滤波器符幅特性为：

因为L是整数，它的符幅特性也由L 1=N/2个正弦分量构成。因为sin(0)=0，在ω=0处，符幅特性必定为A(0)=0。这样的特性显然不能用作低通和带阻滤波器。实际上这种类型的滤波器也不适用于选频滤波。

把这四种类型的线性相位滤波器的相角特性和符幅特性用对比的方法画在一张图上如表2.1。  

表2.1四种类型的线性相位滤波器的相位与符幅特性表

2.6 符幅特性的计算

设系统的系数向量为b=[-1,2,4,2,-1]，a=1，求其幅频特性及符幅特性，并作比较。

先用手工计算，因为a=1,说明这是一个FIR滤波器，观察系数向量b，其长度N=5，阶次为M=N-1=4，传递函数为

                 H(z) = -1 2z -1 4z -2 2 z -3 – z –4

注意h(n)=bn，其频率响应为：

调用abs(H)函数只能求出常为正的幅特性。

要求出符幅特性，必须把对称系数配对计算。将n=0,4组成一对，n=1,3组成一对，n=2处是对称中心，h(2)独立计算，于是有：

由于对称配置，的高次和低次项搭配成为余弦，合成后的相位项都是，故其符幅特性是这三个幅度项的代数和：

把本题的幅特性和符幅特性画成曲线如下图

图2.2(a)是幅特性，图2.2(b)是符幅特性；可见符幅特性取绝对值，就得幅特性。

图2.2

图2.3(a)是相特性，图2.3(b)是线性相位，左图的第二个突跳与幅特性反号对应。

图2.3

2.7 计算符幅特性的MATLAB程序

编了一个通用的程序，它能自动判别滤波器的类型并能在(2.2a)~(2.5a)四个公式中，自动选择正确的符幅特性公式进行计算。把这个函数程序命名为amplres.m，其调用方法：

                 [A,w,type,tao]=amplres(h) 

其中： h=FIR滤波器的脉冲响应或分子系数向量

 A＝滤波器的符幅特性

 w＝取的频率向量，在0到pi之间取501个点

 type＝线性相位滤波器的类型

 tao＝符幅特性的群迟延

子程序amplres.m的核心语句为：

N=length(h);tao=(N-1)/2;

L=floor((n-1)/2);  

n=1:L 1;w=[0:500]*pi/500; % 滤波器频率向量

if all(abs(h(n)-h(N-n 1))<1e-10)   % 若系数对称

    A=2*h(n)*cos(((N 1)/2-n)'*w) - mod(N,2)*h(L 1); 

% 对称条件下计算A（两种类型）

elseif all(abs(h(n) h(N-n 1))<1e-10)&(h(L 1)*mod(N,2)==0)   

% 系数为反对称条件下计算A （两种类型）

A=2*h(n)*sin(((N 1)/2-n)'*w); end       

三、      设计FIR滤波器的窗函数法

3.1 理想低通滤波器的实际滤波器

为了构造特性接近理想低通滤波器的实际滤波器，应当使它具备如下特性：

(1).因果序列，在n<0处无样本，起点在n=0处；

(2).有限长度，长度为N(阶次为N-1),系数下标从0到N-1;

(3).系数对称，以保证线性相位。对称中心应在n0=(N-1)/2处；也即把理想脉冲响应右移n0;

用长N对称中心在n0=(N-1)/2的矩形窗RN(n)与hd(n)相乘，截取出的实际系数序列h(n)， ，如下图3.1：

图3.1

图3.2为理想滤波器的脉冲响应hd(n)和其频谱；

图3.2

图3.3为矩形窗函数R_N(n)和其频谱；

图3.3

图3.4时域为hd(n)和R_N(n)的乘积，其频谱为上两频谱的卷积。出现波动。

图3.4

求理想脉冲函数加窗后的频谱特性(图7.3.1(f))可以利用频域卷积定理。时域相乘h(n)= hd(n)RN(n),对应于频域卷积，做卷积必须考虑函数的正负号，所以其频谱需要用符幅特性来表示：

其中Hd(jω)为hd(n)的频谱

WRN(jω)为RN(n)的频谱，用符幅特性表示。

分析其幅度，可得截断后滤波器的频率特性

这一积分的几何意义见图3.5。

实际的H(ω)和理想的Hd(ω)有两点差别：

（1）通带内增加了波动，最大的峰值在ω=ωc-2π/N处。阻带内产生了余振，最大的负峰在ω=ωc 2π/N处；

（2）在理想特性不连续点ω=ωc附近形成过渡带。过渡带的宽度近似等于ARN(ω)主瓣宽度。

图3.5矩形窗与理想低通滤波器符幅特性的频域卷积

(a).理想滤波器频谱；

(b).矩形窗函数符幅特性；此时相当于ω=0，它与(a)卷积的结果为矩形窗截断区的正面积减负面积，设为1。

(c ).ω=ωc时卷积输出应为阴影面积的积分，此处应为0.5。

(d). ω <ωc时少一块负面积，卷积输出大于1；

(e) ω <ωc时多一块负面积，卷积输出小于0。

(f).卷积输出的波动波形。

3.2 矩形窗截取的实际特性

仔细地分析一下过渡带和波动的特性。我们关心的不是ARN(ω)本身的波动，而是它与理想幅特性卷积后的H(ω)在ω=ωc附近的波动。ARN(ω)从-∞逐步右移，进入-ωc~ωc积分区间，形成H(ω)的过程实际上就是对ARN(ω)进行积分，积分的结果称为累积符幅函数。

因为不能用解析方法计算，这里就用MATLAB进行数值积分。所用的函数为y=cumtrapz(x)，它是按梯形积分方法编程的。

程序hc731如下： 

N=?;dx=0.01;               % 设置初值       

x=[-pi/2:dx:pi/2] 1e-10;    % 自变量数组

AR=sin(N*x)./sin(x);  % 矩形窗符幅特性

H=cumtrapz(AR)*dx/pi;  % 数字积分

subplot(2,1,1),plot(x,AR), % 画符幅特性曲线

subplot(2,1,2),plot(x,H) % 画累积符幅特性曲线

delta1=max(H)-1，% 通带最大波动（绝对值）

delta2=abs(min(H)),  % 阻带最大波动（绝对值）

% 求过渡带宽度

xtrans=dx*(find(H==max(H))-find(H==min(H)))

N=?;dw=0.01;               % 设置初值       

w=[-pi/2:dw:pi/2] 1e-10;           % 自变量数组

AR=sin(N*w)./sin(w);  % 矩形窗符幅特性

H=cumtrapz(AR)*dw/pi;  % 数字积分

subplot(2,1,1),plot(w,AR), % 画符幅特性曲线

subplot(2,1,2),plot(w,H) % 画累积符幅特性曲线

delta1=max(H)-1,  % 通带最大波动（绝对值）

delta2=abs(min(H)),  % 阻带最大波动（绝对值）

% 求过渡带宽度

wtrans=dw*(find(H==max(H))-find(H==min(H)))

3.3 结果

        比较N=10及100的结果，可以看出N只影响过渡带宽度(分别为0.628和0.063)，不影响吉布斯效应波动的幅度(As分别为22及25dB)。

图3.6 N=10（左）和N=100（右）时矩形窗的符幅特性和累积符幅特性