FIR滤波器和IIR滤波器是数字信号处理里躲不开的两个大类。最近整理自己多年攒下来的滤波器设计笔记发现项目里最常用的还是窗函数法。这个方法算是数字滤波器设计的入门第一课就算后面等纹波设计、最小二乘设计讲得再多不少工程项目的首版原型仍然会跳回窗函数法——快、直观、可控。这篇文章围绕窗函数法设计FIR数字滤波器把原理、公式、代码、容易踩的坑一起过一遍。内容偏入门但坑有很多是老手也会犯的所以后面那些实操记录部分值得慢慢看。1. 窗函数法的底层逻辑用最小成本换取可用的滤波器1.1 设计一个FIR滤波器首先要回答为什么不用IIR在动手算系数之前我想先讲清楚一个最基础也最容易被忽略的问题到底什么场景该选FIR什么场景该选IIR。这个选择如果做错了后面所有设计流程都白搭。IIR滤波器的优点是阶数低、计算量小同样的衰减指标IIR可能只需要三五阶FIR可能要用几十甚至上百阶。但IIR的相位响应是严重非线性的信号经过滤波器以后各个频率分量之间的相对时间关系会发生扭曲。很多刚接触数字信号处理的人不太在乎这一点因为用示波器观察幅度响应时滤波效果好像也挺好。但一旦处理的是心电信号、语音信号或者通信基带信号相位非线性带来的波形畸变就会变得非常明显。FIR滤波器虽然阶数高、计算量大但它有一个无可替代的优点只要系数对称就能实现严格的线性相位。线性相位意味着滤波器只对信号做整体延时不改变每个频率分量的相对关系。这一条在工程上的价值极高。举个例子我曾经在一个振动信号采集项目里同时试过IIR和FIR做同一组数据的带通滤波IIR滤波后的波形在时域上看起来完整但和原始信号对照时上升沿的位置发生了微小的偏移而且不同频段的偏移量不一样。换用线性相位FIR以后整个波形的轮廓被完整保留只是整体平移了一个固定时长做后续特征提取就变得非常靠谱。1.2 窗函数法到底解决了一个什么困难问题FIR设计有很多方法窗函数法不是唯一的但它是理解所有其他方法的基础。它的基本思路可以概括成一句话先构造一个理想滤波器的冲激响应然后在时域用窗函数截断它。理想滤波器在频域上是一个“矩形”——通带内增益为1阻带内增益为0过渡带宽度为0。这种滤波器物理上无法实现因为它对应的冲激响应在时域是无限长的而且非因果。窗函数法做的事情就是把这个无限长的响应切出一段有限的来再想办法让“切”的动作尽量温和使频率响应在工程可接受的范围内逼近理想指标。这个方法最大的好处是直观而且不需要迭代优化。等纹波设计比如常见的Parks-McClellan算法需要反复迭代调整系数最小二乘设计要解矩阵方程这些方法在掌握之前很容易让人一头雾水。窗函数法则不一样选择窗函数、估算阶数、乘系数三步走完滤波器的雏形就出来了。对很多实际工程来说这个“雏形”的性能已经足以直接投入使用。我在自己的项目里总结过一个经验窗函数法适合做需求明确、指标在常规范围内的滤波器设计。如果你只需要一个截止频率几百赫兹到几千赫兹的低通、带通或高通阻带衰减要求不高过七八十分贝窗函数法几乎是最省时间的选择。反过来如果你要求极窄的过渡带同时还要极高的阻带衰减那窗函数法可能会因为阶数过大而显得笨重这时才需要转而考虑多级实现或者直接用IIR。2. 三个关键步骤理想冲激响应、截断与加窗2.1 理想低通滤波器的冲激响应是怎么来的窗函数法的一切都从一个数学结果出发理想低通滤波器的冲激响应是一个sinc函数。假设采样频率为 (f_s)理想低通滤波器的截止频率为 (f_c)单位Hz那么这个滤波器在数字域的频率响应从 (-\pi) 到 (\pi) 的范围里在 (|f| f_c) 内是1其余是0。对它做傅里叶反变换就得到时域的序列[ h_d[n] \frac{2f_c}{f_s} \cdot \text{sinc}\left(\frac{2f_c}{f_s}(n - M)\right) ]其中 (\text{sinc}(x) \frac{\sin(\pi x)}{\pi x})(M \frac{N-1}{2})(N) 是我们最终确定的滤波器长度。这个公式里最关键的一点是那个 ((n - M)) 的偏移量。它表示理想冲激响应的对称中心放在了 (M) 处这样截取出来的序列在 (0 \sim N-1) 范围内左右对称滤波器天然满足线性相位条件。初学者最容易犯的错是直接把 (n 0, 1, 2, ...) 代入公式而不作平移结果得到一个看起来七扭八歪的序列频率响应完全不对。sinc函数本身向两端无限延伸衰减速度只有 (1/x)这意味着它的尾巴拖得很长。我们只能截取其中一段这个截取动作直接影响最终滤波器的质量。2.2 矩形窗截断为何会引发吉布斯现象如果我们直接用矩形窗去截等同于在理想冲激响应上乘一个长度为N的全1序列。这个操作在频域里相当于把理想矩形的频谱和一个sinc函数频谱做卷积。结果就是过渡带变宽通带边缘出现振荡阻带衰减始终被限制在大概 -21dB 左右。这个现象叫吉布斯现象。它的本质是傅里叶级数在间断点附近的收敛问题——矩形频谱这个信号本身有突变用有限项傅里叶级数逼近时突变处一定会出现过冲。很多教科书上都画过那张经典图矩形窗截断后的频率响应在截止频率附近上下振荡而且无论N取多大过冲的峰值幅度不会消失只会越来越靠近突变点。也就是说单纯靠增加长度是没办法把阻带衰减压下去的。实际设计滤波器时-21dB的阻带衰减对绝大多数应用都不够。哪怕做一个简单的音频低通这种滤波器也会在阻带留下明显的泄漏成分。2.3 窗函数扼住了吉布斯现象的咽喉吉布斯现象的根源是“突然切断”。既然突然切断会带来麻烦那就换一种切断方式让序列在两端平缓地过渡到0。这就是窗函数法的核心思想。窗函数 (w[n]) 是一个长度为 (N)、关于中心对称、两端向0靠近的序列。我们用 (w[n]) 乘以理想冲激响应 (h_d[n])得到最终的滤波器系数[ h[n] h_d[n] \cdot w[n] ]由于窗函数的两端是逐渐过渡的截断的“切口”被柔化了频域卷积中的旁瓣被压低。代价是主瓣变宽反映在频率响应上就是过渡带变宽。这里有一个工程上必须记住的权衡关系窗函数的主瓣宽度决定过渡带宽旁瓣高度决定阻带衰减。想压低旁瓣主瓣就会变宽想让过渡带窄旁瓣就压不下来。不同的窗函数本质就是在这一对矛盾里取不同的平衡点。3. 常用窗函数横向对比主瓣、旁瓣与阶数估算3.1 各类窗函数的工程特征速览窗函数家族里最常见的几个成员在实际使用中各有各的脾气。矩形窗是理想但不可用的基准线。它主瓣最窄在同等长度下过渡带最窄但旁瓣只压到 -13dB阻带衰减只有约 -21dB。它的意义更多是作为理论参考而不是实际设计工具。汉宁窗Hanning的表达式是 (w[n] 0.5 - 0.5\cos(\frac{2\pi n}{N-1}))旁瓣峰约 -31dB衰减到约 -44dB。它用在功率谱估计里比较多做滤波器也可以但表现中规中矩。海明窗Hamming是工程里最常见的通用窗表达式 (w[n] 0.54 - 0.46\cos(\frac{2\pi n}{N-1}))旁瓣峰约 -41dB阻带衰减约 -53dB。这个衰减水平覆盖了相当多的应用场景。和汉宁窗相比海明窗在靠近主瓣的几片旁瓣压得更低代价是远处旁瓣衰减不再继续下降而是维持在一个水平附近。布莱克曼窗Blackman的表达式里多了二次谐波项 (0.08\cos(\frac{4\pi n}{N-1}))旁瓣更低阻带衰减约 -74dB。但它的主瓣比海明窗宽不少需要在滤波器长度上做出补偿。它适合对阻带衰减要求高于60dB的场景。凯泽窗Kaiser是特殊的一类它带一个可调参数 (\beta)通过调整 (\beta) 可以在主瓣宽度和旁瓣高度之间灵活选择。(\beta) 越大窗越“尖”旁瓣越低但过渡带也越宽。这使得凯泽窗能覆盖从接近矩形窗到极低旁瓣的很大范围是工程上最灵活的选择。3.2 滤波器长度的经验估算公式有了窗函数的特性下一步就是估算滤波器长度。这是窗函数法里最实用的一个技能我希望每个看文章的人都能默写下来。设过渡带宽度为 (\Delta f)单位Hz即阻带起始频率减去通带边缘频率采样率为 (f_s)那么滤波器长度 (N) 约等于[ N \approx \frac{C \cdot f_s}{\Delta f} ]不同窗函数对应的 (C) 值和可达阻带衰减如下表窗函数C值典型阻带衰减矩形窗0.9约 -21 dB汉宁窗3.1约 -44 dB海明窗3.3约 -53 dB布莱克曼窗5.5约 -74 dB凯泽窗可调可按公式精确估算N计算出来后一般取奇数方便得到第一类线性相位的对称滤波器。例如算出52.8我会取53而不是52因为奇数长度的系数关于中间点对称群延迟是整数个采样周期逻辑上最好理解实现上也最清晰。凯泽窗有专门的估算公式。给定阻带衰减 (A_s)单位dB先计算 (\beta)[ \beta \begin{cases} 0.1102(A_s - 8.7), A_s 50 \ 0.5842(A_s - 21)^{0.4} 0.07886(A_s - 21), 21 \le A_s \le 50 \ 0, A_s 21 \end{cases} ]再计算 (N)[ N \approx \frac{A_s - 7.95}{2.285 \cdot \Delta \omega} 1 ]其中 (\Delta \omega \frac{2\pi \Delta f}{f_s})是归一化后的过渡带宽度。这几条公式不管手算还是写进脚本里都很好用。我在很多项目里都会先用它们快速估算阶数再结合频响验证微调。4. 完整实操手算并实现一个1500Hz截止的低通滤波器4.1 指标转换先算理想截止频率我把具体的计算过程完整走一遍这样拿到现成代码就能直接替换参数。假设采样率 (f_s 8000) Hz要求通带边缘频率 (f_p 1500) Hz阻带起始频率 (f_s 2000) Hz阻带衰减 (A_s \ge 50) dB。过渡带宽度[ \Delta f f_s - f_p 2000 - 1500 500 \text{ Hz} ]这里要特别注意一个容易被忽略的细节窗函数法公式里的“截止频率”(f_c)不是通带边缘而是理想低通滤波器的突变点。这个突变点通常会放在过渡带的中间也就是[ f_c \frac{f_p f_{st}}{2} \frac{1500 2000}{2} 1750 \text{ Hz} ]很多教材为了简化会直接用通带边缘当截止频率这样做出来的滤波器实际通带边缘会向阻带一侧偏移。严格的设计应该把理想截止频率放在过渡带中心最后得到的通带和阻带范围才会准确落在我们想要的位置。4.2 阶数估算窗函数法里决定性的一步阻带衰减要求50dB海明窗的典型衰减约 -53dB刚好够用而且海明窗过渡带系数C为3.3比布莱克曼窗更节省阶数。[ N \approx \frac{3.3 \times 8000}{500} 52.8 ]取奇数得到 (N 53)。对称中心[ M \frac{N-1}{2} \frac{52}{2} 26 ]海明窗的实际过渡带宽度约为[ \Delta f \approx \frac{3.3 \times 8000}{53} 498 \text{ Hz} ]和要求的500Hz几乎一致。这说明阶数选择是合适的。4.3 冲激响应与加窗一步步得到最终系数理想低通滤波器的冲激响应为[ h_d[n] \frac{2 \times 1750}{8000} \cdot \text{sinc}\left(\frac{2 \times 1750}{8000}(n - 26)\right) ]化简一下[ h_d[n] 0.4375 \cdot \text{sinc}\left(0.4375(n - 26)\right) ]其中 (n 0, 1, 2, ..., 52)。海明窗为[ w[n] 0.54 - 0.46 \cdot \cos\left(\frac{2\pi n}{52}\right) ]最终滤波器系数[ h[n] h_d[n] \cdot w[n] ]由于 (h_d[n]) 和 (w[n]) 都关于 (n26) 对称最终系数也严格对称。这个对称性是线性相位的基础也是窗函数法最可靠的一点——只要公式写对线性相位天然成立。4.4 用Python验证设计频响曲线怎么读我用一段精简的Python代码展示完整的实现过程。脚本里尽量不加多余的依赖核心是NumPy和SciPy。import numpy as np import matplotlib.pyplot as plt from scipy.signal import freqz fs 8000 # 采样率 fp 1500 # 通带边缘频率 fst 2000 # 阻带起始频率 delta_f fst - fp # 海明窗阶数估算 N int(np.ceil(3.3 * fs / delta_f)) if N % 2 0: N 1 M (N - 1) // 2 fc (fp fst) / 2 # 理想截止频率放在过渡带中心 n np.arange(N) hd 2 * fc / fs * np.sinc(2 * fc / fs * (n - M)) w 0.54 - 0.46 * np.cos(2 * np.pi * n / (N - 1)) h hd * w # 归一化直流增益为1便于后续级联使用 h h / np.sum(h) # 频率响应 w_axis, H freqz(h, 1, worN2048, fsfs) H_db 20 * np.log10(np.abs(H) 1e-12) plt.figure(figsize(10, 5)) plt.plot(w_axis, H_db) plt.axvline(1500, colorred, ls--, lw1) plt.axvline(2000, colorred, ls--, lw1) plt.ylim(-80, 5) plt.grid(True) plt.ylabel(幅度 (dB)) plt.xlabel(频率 (Hz)) plt.show()这段代码输出的频响图里应该看到三条重要信息第一通带边缘大约在1500Hz附近增益只波动零点几个dB这表明通带纹波很小。第二2000Hz之后曲线迅速跌到 -50dB 以下最终停在约 -53dB说明阻带衰减符合海明窗的理论值。第三整条曲线没有出现在阻带内上升到很高再跌下去的情况说明阶数估算合适没有严重的频谱泄漏。4.5 MATLAB版本两条命令就能微调实际项目里很多团队还是在用MATLAB做信号处理原型。MATLAB的滤波器设计工具箱提供更快捷的接口但理解背后的手动过程仍然重要。fs 8000; fp 1500; fst 2000; df fst - fp; N ceil(3.3 * fs / df); if mod(N, 2) 0 N N 1; end fc (fp fst) / 2; M (N - 1) / 2; n 0 : N - 1; hd 2 * fc / fs * sinc(2 * fc / fs * (n - M)); w 0.54 - 0.46 * cos(2 * pi * n / (N - 1)); h hd .* w; % 归一化直流增益 h h / sum(h); freqz(h, 1, 2048, fs);MATLAB的sinc函数与Python中numpy.sinc的定义一致都是 (\sin(\pi x)/(\pi x))所以公式可以直接平移。实际开发中如果想快速做一个滤波器原型我经常还会用fir1这条命令h fir1(N-1, fc/(fs/2), hamming(N));fir1内部就是窗函数法它要求截止频率用Nyquist频率归一化。注意它内部的截止频率也是放在过渡带中心的所以传入fc/(fs/2)是合理的。有些人不理解这一点直接传入通带边缘1500Hz做出来的滤波器经过验证就会发现实际通带边缘偏了。5. 从低通扩展到高通、带通设计5.1 类型I与类型II线性相位的区别FIR滤波器按系数的对称性可以分成几类其中最常用的是第一类线性相位系数偶对称长度取奇数。这种滤波器在零频和Nyquist频率都没有额外约束低通、高通、带通都能用。如果长度取偶数那就是第二类线性相位频率响应在Nyquist频率处必然有零点。这个约束意味着它天然不适合做高通滤波器因为高通在Nyquist频率处应该有最大增益结果被强制压成0。因此设计高通、带阻这类在(f_s/2)附近仍有通带的滤波器时尽量把N取成奇数。这也是我前面一直强调“取奇数”的原因。5.2 高通、带通、带阻的系数推导从低通系数出发可以通过频谱搬移得到其他类型。高通滤波器的理想冲激响应可以看成“全通”减去“低通”[ h_{d_hp}[n] \delta[n - M] - h_{d_lp}[n] ]其中(\delta[n-M])是在对称中心位置的单位脉冲。直观理解所有频率都通过再把低频部分抠掉剩下的就是高频。带通滤波器可以看成两个低通的差。设下截止频率为(f_{c1})上截止频率为(f_{c2})[ h_{d_bp}[n] \frac{2f_{c2}}{f_s} \cdot \text{sinc}\left(\frac{2f_{c2}}{f_s}(n - M)\right) - \frac{2f_{c1}}{f_s} \cdot \text{sinc}\left(\frac{2f_{c1}}{f_s}(n - M)\right) ]带阻滤波器则是在全通基础上扣除带通[ h_{d_bs}[n] \delta[n - M] - h_{d_bp}[n] ]得到理想冲激响应后后续加窗过程与低通完全一致直接用选好的窗函数相乘即可。在这个部分我额外提醒一点带通滤波器有两个过渡带阶数估算时要以更窄的那个过渡带为准。比如一个滤波器下边过渡带400Hz、上边过渡带600Hz算阶数时用400Hz来算否则上边下边会有一个指标不达标。6. 实际开发中常见问题与排查经验6.1 过渡带不够窄先算N再验证不要盲目加长我见过很多人设计完滤波器频率响应一看过渡带比预期宽出一截第一反应就是把N翻倍甚至翻三倍重算。这当然有效但不是最优解。正确流程是先确认阶数估算公式里用的(\Delta f)对不对。如果理想截止频率没放在过渡带中心实际表现就会整体平移被误判成“过渡带太宽”。还有一种情况是过度追求阻带衰减选用了布莱克曼窗C值从海明窗的3.3变成5.5那相同过渡带下N自然要大得多这不是“滤波器变差”而是性能指标的必然换取。如果在阶数正确的前提下仍然觉得过渡带偏宽我非常不建议继续无限加长N。FIR阶数一上去计算量和资源占用成倍增长。这时更务实的方案是换用多级滤波器或者重新审视指标里有没有一条真正需要那么窄的过渡带。6.2 频率响应通带边缘和相位曲线出现异常我在调试中遇到过好几次通带边缘位置偏了、相位曲线不直线的情况。排查下来原因高度集中在两个地方。第一理想冲激响应公式里的对称中心算错。比如N53时M(53-1)/226如果模板代码里写成了MN//2那对称中心偏移半个点相位自然是非线性的。第二归一化截止频率传参错误。用现成的firwin、fir1时如果不清楚它们内部对截止频率的约定把通带边缘和理想截止频率混着传最后得到的通带边缘就会和设计目标不一致。我在我的代码里习惯把所有频率先换算成Hz不在代码里写“归一化到Nyquist的浮点数”这样至少能减少一类错误。6.3 嵌入式定点实现的精度坑滤波器系数算出来是浮点数直接放进嵌入式设备经常出问题。我用过一个32位定点芯片做实时滤波一开始把系数简单地转成定点结果在阻带衰减要求较高的频段实测幅频特性比浮点仿真差了将近10dB。后来排查发现问题出在系数动态范围上。窗函数法得到的系数有大有小直接定点量化会损失小系数的精度而恰恰是那些小系数决定了阻带的细微结构。解决办法有两个方向。一个是在量化前对系数做整体缩放把系数范围放大到接近定点满量程再量化。另一个是把滤波器拆成多相结构或者用高一点的中间位宽做累加避免乘累加过程中的截断误差。具体用哪种要看芯片算力和存储空间。但有一点是确定的浮点仿真结果不能直接当作定点实现的最终指标必须在硬件平台上重新测一遍频响。6.4 遇到苛刻指标时如何妥协窗函数法有一个绕不开的缺点指标一旦特别严苛N会变得很大。比如说要求过渡带只有几十赫兹、阻带衰减超过80dB用布莱克曼窗或凯泽窗算出来的N可能上千。直接把上千阶的FIR跑在实时系统里延迟和计算量都可能受不了。这时候我会老实评估两种情况其一能不能接受IIR滤波器的相位非线性把设计切换到IIR其二能否把任务拆成两级前级用一个宽过渡带的滤波器粗滤后级用一个窄过渡带的滤波器精滤两个滤波器的总阶数通常比单个滤波器小很多。我在一个信号处理Demo里遇到过类似情况单个1600阶FIR滤波器的延迟让整套系统的反馈环不稳定拆成三级以后每级只有200300阶总延迟反而降到了可接受范围。窗函数法不是终点是工具箱里的一个底板实际工程要灵活搭配其他方法。6.5 直流增益和归一化问题最后一个非常常见的坑是忘记检查直流增益。要知道理想低通的冲激响应在(nM)处取值为(2f_c/f_s)加窗后所有系数的和理论上接近1但由于截断效应实际可能有1%左右的偏差。这个偏差在单级滤波时无伤大雅一旦多个滤波器级联每个都偏一点点级联总增益就可能偏差好几个百分点。我在实际项目里的处理习惯是生成系数后先求sum(h)然后统一除以它把直流增益严格归一化为1。这样做的代价是通带增益可能不再是精确的1但换来的是系统总增益可控后续再做增益补偿时更容易。低通滤波器归一化直流增益带通滤波器也可以归一化中心频率处的增益方法类似。7. 最后聊一点个人体会窗函数法用了这么多年给我的最大感受就是即使有了各种现成设计工具手动把理想冲激响应、窗函数、阶数估算这些步骤走一遍仍然很有必要。因为它让我始终清楚一个滤波器系数背后对应着什么样的频域效果指标不达标时能快速判断是哪个环节出了问题。另外我也越来越推崇“先手动算一遍再用工具快速验证”的工作方式。手动算虽然慢但每一步的逻辑都在掌控之中工具验证虽然快但内部做了很多自动处理很容易掩盖错误。把这两者结合起来既能保证设计质量又能提高迭代效率。这一篇讲的是窗函数法设计FIR滤波器也是整个FIR和IIR设计系列的第一篇。后面我会继续写更多关于滤波器设计的实战内容包括等纹波设计、IIR设计、以及不同实现结构对硬件资源的影响。如果这篇文章对你有帮助或者你有踩过其他我没提到的坑欢迎私下交流互相补全经验库。