FIR 数字滤波器的设计

上一章我们学习了 IIR 数字滤波器的设计方法--借助模拟低通原型(巴特沃斯、切比雪夫、椭圆),通过双线性变换把成熟的模拟设计"翻译"到数字域。IIR 滤波器可以用很低的阶次实现陡峭的过渡带,但它有一个挥之不去的痛点:相位非线性。在语音处理、图像处理、通信系统等对信号波形保真度要求高的场景中,非线性相位会造成波形失真。FIR 滤波器虽然阶次更高、计算量更大,却能实现严格的线性相位,并且天然绝对稳定。本节将系统介绍 FIR 滤波器的三种主流设计方法:窗函数法、频率采样法、等波纹(切比雪夫)最优逼近法。

FIR 和 IIR 的根本区别在于冲激响应是否有限长。这个看似简单的差异,直接决定了两者在结构、稳定性、相位特性和设计方法上的全面分野:

用一个最小例子来直观感受两者的区别。输入脉冲 $x(0)=1$,其余 $x(n)=0$,观察输出如何演化。

前 3 步输出完全一样,但第 4 步 FIR 归零了,IIR 还在继续衰减。原因是 IIR 差分方程里有 $y(n-1)$--反馈:输出反过来喂给自己。FIR 只看输入的历史 $x(n-k)$,不看自己的历史输出。

从这两个数值例子出发,可以推导出 FIR 与 IIR 在稳定性、相位和阶次上的全部差异:

Part 0 已经回答了"FIR 是什么"——冲激响应有限长、结构非递归、天然绝对稳定。接下来的问题是:拿到系数 $h(n)$ 之后,硬件或代码里具体怎么搭?

FIR 没有反馈路径,所以所有结构都围绕有限长冲激响应的不同组织方式展开:有时域的横截型、有转置型、有频域的频率抽样型、还有利用对称性省乘法的线性相位型。理解这些结构之后,我们再深入线性相位 FIR 的本质(第二组)与具体设计方法(第三组)。

FIR 的系统函数是 $z^{-1}$ 的 $N-1$ 次多项式,没有分母(除 $z=0$ 处),系统在有限 $z$ 平面处处收敛。

没有反馈意味着 FIR 天然稳定:只要输入有界,有限个输入样本的加权和一定有界。它还可以通过 $h(n)$ 的偶对称或奇对称设计出严格线性相位,因此在语音、通信、图像这类关心波形形状的场景中非常常用。代价是同样的幅频指标下,FIR 往往需要比 IIR 更高的阶数。

横截型直接用差分方程实现:$N-1$ 个延迟单元组成延时链,$N$ 个抽头分别乘 $h(0),h(1),\ldots,h(N-1)$,再相加。这就是 LTI 系统的卷积和公式。

横截型的优点是最直观:每个抽头对应一个冲激响应样本,改某个 $h(k)$ 就是改第 $k$ 条抽头的乘法系数。它常作为理解 FIR 的基准结构,也常用于系数需要频繁更新的自适应 FIR。

用转置定理把横截型变换为等价的转置横截型:将所有支路方向反转,输入输出互换,系统函数 $H(z)$ 不变。转置后,输入信号先分配到各个乘法器,再沿着延迟链逐级累加,延迟单元数仍为 $N-1$。

级联型把 $H(z)$ 分解为实系数一阶或二阶因子的乘积:

$N$ 为偶数时共 $N/2$ 个二阶节;$N$ 为奇数时有一个一阶节($\beta_{2k}=0$)。每节主要对应一对零点,适合从零点分布理解和调节 FIR;缺点是需要先因式分解,且节次顺序和系数缩放会影响有限字长误差。

频率抽样型利用 DFT 频域采样定理,把 $H(z)$ 用 $N$ 个频率采样值 $H(k)$ 表示:

这可以看成两部分级联:

• 梳状滤波器 $1-z^{-N}$:在单位圆上有 $N$ 个等间隔零点 $\omega_k=2\pi k/N$,频率响应为 $|2\sin(N\omega/2)|$。
• 谐振器组 $\sum H(k)/(1-W_N^{-k}z^{-1})$:每个一阶节在单位圆上有一个极点 $\omega_k=2\pi k/N$,正好与梳状滤波器的零点抵消,使得 $\omega=2\pi k/N$ 处的频率响应等于 $H(k)$。

当 FIR 滤波器的单位冲激响应满足偶对称 $h(n)=h(N-1-n)$ 或奇对称 $h(n)=-h(N-1-n)$ 时,相位响应具有严格的线性特性:

• 偶对称(情况 I、II):$\theta(\omega)=-\frac{N-1}{2}\omega$,纯线性相位。
• 奇对称(情况 III、IV):$\theta(\omega)=-\frac{N-1}{2}\omega+\frac{\pi}{2}$,线性相位加 $90°$ 相移。

利用对称性,可以把乘法次数减少约一半。根据 $N$ 的奇偶和对称性,线性相位 FIR 有四种情况:

$N$ 为奇数(情况 I)时,对称中心在整数样点 $n=(N-1)/2$,所有 $h(n)$ 都可配对。结构上,对称位置的两个延迟输出相加(偶对称取 $+1$)或相减(奇对称取 $-1$),再乘以系数 $h(n)$,乘法次数从 $N$ 减为 $(N+1)/2$。

$N$ 为偶数(情况 II)时,对称中心在 $n=(N-1)/2$ 处(非整数样点),没有单独的中心项。乘法次数从 $N$ 减为 $N/2$。

线性相位 FIR 在通信、图像处理中非常重要,因为线性相位意味着所有频率分量延迟相同,信号波形不会失真。下面的第二组会深入展开线性相位的本质与四类约束。

时域卷积 $y(n)=x(n)*h(n)$ 在频域等价于乘法 $Y(k)=X(k)\cdot H(k)$。当 FIR 滤波器的单位冲激响应 $h(n)$ 较长(比如 $N>64$)时,直接时域卷积需要 $N$ 次乘累加,而用 FFT 在频域做乘法的计算量反而更小。

当 $N$ 取 2 的幂时,FFT 的计算复杂度为 $O(N\log N)$,而直接卷积为 $O(N\cdot M)$。当 $M$ 较大时 FFT 方法显著更快。

实际应用中输入信号通常是无限长的流式数据(如音频),需要分段处理。两种经典方法:

重叠相加法(Overlap-Add)

将长序列 $x(n)$ 分成每段长 $L$ 的无重叠子段 $x_m(n)$。每段与 $h(n)$ 做 $N$ 点循环卷积($N\ge L+M-1$),结果长度为 $L+M-1$,相邻段的输出有 $M-1$ 个样本重叠,将重叠部分相加即可。

流程图:

graph LR
  A["x(n) 分段"] --> B["补零到 N 点"]
  B --> C["N 点 FFT"]
  C --> D["乘以 H(k)"]
  D --> E["N 点 IFFT"]
  E --> F["重叠相加"]
  F --> G["y(n)"]

重叠保留法(Overlap-Save)

每段输入取 $N$ 点,但相邻段有 $M-1$ 点重叠(每段前 $M-1$ 个点是上一段的后 $M-1$ 个点)。每段做 $N$ 点循环卷积后,丢弃前 $M-1$ 个受缠绕效应影响的样本,保留后 $L=N-M+1$ 个有效样本拼接输出。

至此,FIR 部分已经从时域横截型讲到频域快速卷积。可以看出,FIR 和 IIR 的结构选择标准不同:IIR 首先要管住极点和反馈误差,FIR 则更多是在乘法次数、频域调节、线性相位和分段卷积之间做取舍。接下来的第二组,我们将深入分析线性相位 FIR 的本质与约束——这是 FIR 区别于 IIR 的核心优势所在。

在深入具体设计方法之前,我们必须先回答一个根本问题:FIR 滤波器凭什么能做到线性相位?要做到这一点,它的冲激响应 $h(n)$ 必须满足什么对称条件?这些对称条件又把滤波器分成了哪几类,每类各自能做什么、不能做什么?

IIR 滤波器依靠递归结构(反馈)实现,系统函数 $H(z)$ 的极点可以靠近单位圆,从而用低阶次获得很窄的过渡带。但代价是:极点位置与相位特性深度耦合,无法独立控制。FIR 滤波器没有反馈,$H(z)$ 只有零点没有极点,这给它带来了两个天然优势:

• 绝对稳定:没有极点,不存在稳定性问题。
• 线性相位可实现:只要 $h(n)$ 满足对称条件,群延迟就是常数。

线性相位的物理意义是:信号通过滤波器后,所有频率分量被延迟的时间完全相同,波形只是整体平移,不会发生色散失真。这在图像滤波、数据通信、雷达信号处理中至关重要。

这句话的数学表述是:如果滤波器的相位响应为 $\varphi(\omega)=-\omega\tau$(关于 $\omega$ 的线性函数,斜率为常数 $\tau$),那么频率为 $\omega$ 的分量被延迟 $\tau$,频率为 $2\omega$ 的分量也被延迟 $\tau$--所有频率的延迟完全相同。结果就是波形形状不变,只是整体晚了 $\tau$ 个采样点。反过来,如果 $\varphi(\omega)$ 不是线性的,不同频率延迟不一样,波形散了,这就是相位失真。

线性相位的数值验证。上图:因果 5 点滑动平均(输出写在 $i+4$),系数不对称,橙线相对蓝线明显右移(相位非线性)。下图:对称 5 点滑动平均(输出写在 $i+2$ 即窗口中心),系数对称,橙线和蓝线完全重合,只去掉噪声没有位移(线性相位=恒定延时)。

图中对比了同一个 5 点滑动平均的两种写法:

• 上图(因果实现):输出写在 $i+4$,系数序列 $[\frac{1}{5},\frac{1}{5},\frac{1}{5},\frac{1}{5},\frac{1}{5}]$ 关于中心不对称(因果约束),相位非线性,滤波结果相对原始信号有明显位移。
• 下图(对称实现):输出写在 $i+2$(窗口中心),同样的系数但输出位置对齐了窗口中心,等效于 $h(n)=h(N{-}1{-}n)$,相位线性,滤波结果与原始信号完全重合,只有平滑效果没有位移。

这正是 $h(n)=h(N{-}1{-}n)$ 的物理后果:对称系数保证所有频率分量被延迟相同的时间 $\tau=\frac{N-1}{2}$,波形保真。

FIR 滤波器的系统函数为有限项多项式:

其频率响应 $H(e^{j\omega})$ 可以写成幅度与相位的乘积。为了得到线性相位,需要相位 $\varphi(\omega)$ 与 $\omega$ 成严格的线性关系。数学上可以证明,这要求冲激响应 $h(n)$ 关于中心点 $\alpha=\frac{N-1}{2}$ 满足对称性:

为什么对称能保证线性相位?直观理解:对称的 $h(n)$ 相当于把信号按中心"镜像折叠"后再加权求和,这种操作对每个频率分量的延迟完全一致,因此相位斜率是常数。

这里有一个常见的误解需要澄清。对称条件 $h(n)=h(N{-}1{-}n)$ 要求滤波器的"中心"在第 $\frac{N-1}{2}$ 个采样点,也就是说要等到输入信号的"未来"样本到达后才能计算当前输出。有人把这叫作"非因果"--但这个词容易误导。更准确的说法是非即时:

• 频域视角:线性相位意味着 $\varphi(\omega)=-\omega\cdot\tau$,即相位随频率线性变化,斜率恒为 $\tau=\frac{N-1}{2}$。这不是"相位出问题了",而是所有频率分量统一偏移了一个固定斜率。
• 时域视角:这个固定斜率 $\tau$ 就是恒定延时--输出波形整体晚 $\frac{N-1}{2}$ 个采样点出现,形状完全不变。

所以"非因果"不是真的违反因果律,而是为了对称而付出的延迟代价。因果系统完全可以实现线性相位,只要接受输出不是"即时"的,而是固定延迟 $\frac{N-1}{2}$ 个采样。在实时处理中,这个延迟是可预测、可补偿的;在离线处理中,直接把输出对齐回来就行。前面实验图里对称滑动平均的橙线和蓝线完全重合,就是因为输出写在窗口中心 $i+2$,自动补偿了这个延迟。

对称条件有"偶/奇"两种,滤波器长度 $N$ 有"奇/偶"两种,组合起来正好四类。每一类的频率响应形式和适用滤波器类型各不相同,这是 FIR 设计中最核心的约束表,做题时必须首先查表确定类型。

第一类与第二类线性相位 FIR 的对比总结。共同点是都满足偶对称 $h(n)=h(N-1-n)$，因此相位都是严格线性的：$\phi(\omega)=-\tau\omega$。区别在于：第一类 $N$ 为奇数，对称轴正好落在中心采样点上，因此幅度函数里有一个独立的中心项 $h(\tau)$；第二类 $N$ 为偶数，对称轴落在两个采样点之间，因此没有独立中心项，只剩成对余弦项，这也导致它必有 $H(\pi)=0$。

最近几张课件图还可以和这里形成一条完整的理解链：

• 理想低通图给出窗函数法的目标：矩形幅频响应 $|H(e^{j\omega})|$ 在时域对应无限长的 sinc 序列 $h_d(n)=\frac{\sin[\omega_c(n-\tau)]}{\pi(n-\tau)}$，这解释了为什么后面必须“加窗”。
• 5 点滑动平均例题说明“偶对称 + 奇数长度”就是第一类线性相位 FIR，也是最简单的线性相位低通。
• $h(n)=\{2,2,3,0,-3,-2,-2\}$ 例题说明“奇对称 + 奇数长度”属于第三类，并且天然有 $H(0)=0$，所以能滤除直流。
• 汉明窗例题说明只要构造出关于中心 $\tau=\frac{N-1}{2}$ 偶对称的窗序列，就能得到线性相位；当 $N=31$ 时，它属于第一类。

线性相位 FIR 的零点不是随意分布的,它们具有严格的共轭倒数对称结构。若 $z_i$ 是 $H(z)$ 的零点,则 $z_i^*$、$1/z_i$、$1/z_i^*$ 也必然是零点--四个零点组成一组"四重对称"。

这种结构来源于线性相位条件的代数要求:$H(z)=\pm z^{-(N-1)}H(z^{-1})$。它强制 $H(z)$ 的零点关于单位圆镜像对称,同时关于实轴共轭对称。

• 当零点恰好落在单位圆上($|z_i|=1$)时,$z_i$ 与 $1/z_i^*$ 重合,退化为二重零点。
• 当零点恰好落在实轴上($z_i=\pm 1$)时,四个点全部重合为单点,且 $\omega=0$ 或 $\omega=\pi$ 处必有过零。

明确了线性相位的分类约束后,核心问题转化为:如何由给定的频率指标 $H_d(e^{j\omega})$,求出满足线性相位条件的有限长 $h(n)$? 三种方法对应三条不同的路径:时域截断(窗函数法)、频域采样(频率采样法)、误差最小化(等波纹最优法)。

窗函数法：理想响应的时域截断

最直观的设计思路是:先写出理想滤波器的频率响应 $H_d(e^{j\omega})$,通过 IDTFT 得到无限长的理想冲激响应 $h_d(n)$,然后用一个有限长的"窗" $w(n)$ 把它截断为 $N$ 点,得到可实现的 $h(n)=h_d(n)\cdot w(n)$。

理想低通滤波器的频域与时域对应关系。上图是“砖墙式”理想低通幅频响应：$|H_d(e^{j\omega})|=1$（$|\omega|\le \omega_c$），其余频率完全抑制；下图是对应的单位冲激响应 $h_d(n)$，呈无限长的 sinc（Sa）序列：$h_d(n)=\frac{\sin[\omega_c(n-\tau)]}{\pi(n-\tau)}$，在 $n=\tau$ 处取极限值 $h_d(\tau)=\omega_c/\pi$。这正是窗函数法的起点：先得到理想目标序列，再用窗把它截成有限长 FIR。

这里有一个根本问题：为什么要加窗？理想低通滤波器的冲激响应是无限长的 sinc 函数 $h_d(n)=\frac{\sin[\omega_c(n-\alpha)]}{\pi(n-\alpha)}$，向两侧无限延伸。一个无限长的 $h(n)$ 在物理上无法实现——你不能存储无限个历史样本，也不能做无限次乘法。FIR 滤波器的定义就是冲激响应有限长：$h(n)$ 只在 $n=0,1,\ldots,N-1$ 处非零。所以必须截断。窗函数就是截断的工具：$w(n)$ 在 $0\le n\le N-1$ 内非零、外面为零，乘以 $h_d(n)$ 后就把无限长序列"剪"成了 $N$ 点。

反过来，正因为 FIR 允许我们自由选择保留哪 $N$ 个点、怎么过渡到零，所以才有了"加窗"这个设计空间。不同窗函数的区别就在于截断的方式不同：矩形窗是硬切，汉宁/汉明窗是软着陆。截断方式不同，频域的"涂抹"效果就不同——这就是窗函数法设计的全部内容。

从频域看,时域乘积对应频域卷积:

因此,$H(e^{j\omega})$ 的形状由窗函数频谱 $W(e^{j\omega})$ "涂抹"理想响应得到。窗函数的频谱特性直接决定了设计结果的三大指标:

• 过渡带宽度:由窗频谱主瓣宽度决定,近似为 $\Delta\omega\approx c\cdot\frac{2\pi}{N}$。
• 通带波纹:由窗频谱主瓣与理想响应的卷积造成。
• 阻带衰减:由窗频谱旁瓣的峰值相对幅度决定。

五种常用窗函数

窗的选择是窗函数法的核心。以下五种窗从“最粗糙”到“最平滑”递进，主瓣越窄则过渡带越陡，旁瓣越低则阻带衰减越大——二者不可兼得。

图 7.11 五种常用窗函数的时域波形。所有窗关于中心 $n=\frac{N-1}{2}$ 对称，峰值 $w(n)=1$。矩形窗（最宽、截断最硬）→ 三角窗 → 汉宁窗 → 海明窗 → 布莱克曼窗（最窄、边缘衰减最柔和）。底部公式：$h(n)=h_d(n)\cdot w(n)$，加窗后的频率响应 $H(e^{j\omega})=\sum_{n=0}^{N-1}h(n)e^{-j\omega n}$。

从图中可以看出窗函数设计的核心权衡：矩形窗的“剪刀”最锋利（主瓣最窄，过渡带最陡），但截断太硬导致频谱旁瓣很高；越往右的窗边缘衰减越柔和，旁瓣越低，但主瓣越宽（过渡带越宽）。五种窗的关键参数对比见下表。

工程上最常用的是汉明窗:它在主瓣宽度与旁瓣衰减之间取得了最佳平衡,阻带衰减约 53 dB,足以满足大多数应用。

上图展示了 $N=32$ 时矩形窗(红色)与汉明窗(蓝色)的时域波形。汉明窗在两端逐渐衰减到接近零,这种"软着陆"显著抑制了频谱旁瓣,代价是主瓣变宽一倍。

窗函数法设计步骤

窗函数法的完整设计流程可以归纳为四步,考试时按此顺序书写,不容易遗漏。

其中 $\omega_c$ 通常取通带边缘 $\omega_p$ 与阻带边缘 $\omega_{st}$ 的算术中点:$\omega_c=\frac{\omega_p+\omega_{st}}{2}$。这是因为加窗后通带和阻带各向中间"侵蚀"了一段。

频率采样法：从频域直接构造

窗函数法走的是"时域截断"路线。另一条思路是:既然 $h(n)$ 是 $N$ 点有限长序列,它的 DTFT 在单位圆上采样 $N$ 点正好对应 DFT。那么,何不直接在频域"画"出想要的 $H_d(k)$,再反变换回时域?这就是频率采样法。

频率采样法的优点是直接:你把想要的频响"画"在 $N$ 个采样点上,IDFT 后就自动得到 $h(n)$。缺点也很明显:采样点之间的频率响应由内插函数"连接",如果采样点处跳变太剧烈(比如从通带 1 直接跳到阻带 0),内插会在中间产生较大的起伏,导致阻带衰减不够。

线性相位约束下的采样点相位

为保证线性相位,$H(k)$ 的相位不能随意取,必须满足与四类对称对应的约束。以最常见的第一类(偶对称、$N$ 奇数)为例,幅度采样 $H_k$ 必须为偶对称 $H_k=H_{N-k}$,相位采样则为

第二类(偶对称、$N$ 偶数)还需要额外满足 $H_{N/2}=0$,因为 $H(\pi)=0$。

频率采样法通常分为 I 型和 II 型,区别仅在于第一个采样点的位置:

• I 型:第一个采样点在 $\omega=0$ 处($k=0$),即标准的 $N$ 点 DFT 采样。
• II 型:第一个采样点在 $\omega=\frac{\pi}{N}$ 处,采样点整体偏移半个间隔。

两种类型在通带平坦度和设计灵活性上略有差异,但核心思想完全相同。

过渡带采样的技巧

频率采样法最经典的问题出现在低通设计:通带最后一个采样点为 1,阻带第一个采样点为 0,中间没有任何过渡。内插函数在跳变边缘会产生类似吉布斯现象的振荡,阻带衰减通常只有 20 dB 左右,远不能满足工程需求。

解决方法是:在通带与阻带之间故意插入一个或几个过渡带采样点,取介于 0 和 1 之间的值。这相当于在频域"修坡",让跳变变缓。代价是过渡带变宽,但阻带衰减可以大幅提高。

上图示意了三种情况:理想矩形频响(蓝色实线)、无过渡带采样时的频响起伏(红色虚线,吉布斯振荡明显)、插入过渡带采样后的平滑结果(绿色)。过渡带采样"牺牲"了过渡带的陡峭度,换取了阻带内更干净的高衰减。

等波纹逼近法：最优设计

窗函数法和频率采样法都有一个共同缺陷:它们无法精确控制通带波纹 $\delta_1$ 和阻带衰减 $\delta_2$,设计结果是"碰运气"的--阶数 $N$ 选大了浪费计算,选小了指标不达标。有没有一种方法,能在给定阶数 $N$ 的约束下,让最大逼近误差最小化?

答案是等波纹最优逼近法,也称切比雪夫逼近法或 Parks-McClellan 算法。它的数学目标非常清晰:

为什么叫"等波纹"?因为最优解有一个惊人的性质:误差函数 $E(\omega)=W(\omega)[H_d(\omega)-P(\omega)]$ 在逼近区间 $A$ 上的极值点数目达到最大,且所有极值点的绝对值相等,形成均匀分布的波纹。直观上,这就像把一块弹性膜绷紧在理想响应上,膜上的波峰波谷高度完全一致--没有任何一处"特别突出"。

交错定理与 Remez 交换算法

等波纹最优解的存在性和唯一性由交错定理保证。

基于交错定理,Parks 和 McClellan 提出了高效的迭代算法--Remez 交换算法。它的核心思想是:

1. 在逼近区间 $A$ 内猜测一组 $r+2$ 个极值点频率 $\{\omega_i\}$。
2. 解线性方程组,求出这组极值点上的等波纹误差值 $\delta$ 和多项式系数。
3. 计算整个区间上的误差函数 $E(\omega)$,找出真正的极值点。
4. 用新的极值点替换旧的猜测,重复迭代直到收敛。

Remez 算法收敛非常快,通常只需 5-10 次迭代即可达到机器精度。MATLAB 的 firpm(原 remez)和 SciPy 的 scipy.signal.remez 都是这一算法的实现。

现在我们已经掌握了 FIR 滤波器的三种设计方法。窗函数法概念最直观、手算最方便;频率采样法适合需要精确控制频域某些点响应的场景;等波纹法在给定阶数下性能最优。下面从多个维度进行系统对比。

无论用哪种方法,FIR 滤波器设计都可以归纳为以下决策流程:

flowchart TD
    A[给定滤波器指标] --> B{需要线性相位?}
    B -->|是| C[FIR]
    B -->|否| D[IIR]
    C --> E{阶次受限且指标苛刻?}
    E -->|是| F[等波纹最优法
Parks-McClellan]
    E -->|否| G{需要手算/快速原型?}
    G -->|是| H[窗函数法
汉明/汉宁/布莱克曼]
    G -->|否| I[频率采样法
配合过渡带优化]
    F --> J[量化系数
实现为 DSP/FPGA]
    H --> J
    I --> J