窗函数法FIR带阻滤波器:从原理到工程落地
2026/9/19 13:18:21 网站建设 项目流程

简介:面向通信工程等专业数字信号处理课程设计的一份完整方案文档,针对基于窗函数法的FIR带阻滤波器设计需求,给出从指标分析、窗函数选择到MATLAB实现与频率响应验证的详细流程。包体为单个doc文档,大小312KB,包含课程设计任务书、摘要、MATLAB简介、窗函数设计法原理、线性相位分析、基本窗函数对比、方案设计程序与分析以及总结体会等内容,适合需要快速理解FIR带阻滤波器设计方法并撰写课程设计报告的本科生参考。已有306人学习下载,可帮助读者借助MATLAB完成滤波器设计、掌握矩形窗、三角窗、汉宁窗等不同窗函数的适用场景,并通过幅频响应曲线验证所设计滤波器是否满足通带最大衰减和阻带最小衰减等指标要求。整体内容结构清晰,既有理论推导又有程序实现,是数字信号处理实践入门和课设参考的实用资料。

1. 带阻 FIR 设计的起点:用窗函数法把指标写进系数里

在现场采集心电、振动或音频信号时,最常见的干扰不是宽带随机噪声,而是频率确定的窄带分量:50Hz 工频及其谐波、机械共振、载波泄漏。这类干扰幅度大、带宽窄,用低通或高通都切不干净,只有带阻滤波器能精确削掉一段。基于窗函数法的 FIR 带阻滤波器在软件实现层面成本低、相位线性、系数完全可预期,是工程上最不容易出错的方案之一。窗函数法的核心思路很反直觉:先构造一个理论上完美但无限长的理想滤波器,再用一条窗曲线把它截短成可计算的有限系数。截短带来的通带纹波和过渡带宽,全部由窗函数的选择决定,跟滤波器本身无关。这套逻辑一旦理清,设计带阻、带通、高通就只是改两个频率参数的事。本文按照“原理—设计—验证—落地”的顺序,把一条完整可跑的 FIR 带阻滤波器实现路径讲透,面向需要自己写信号处理代码,而不是只调库的工程师。

2. 窗函数法设计带阻 FIR 的原理:从理想冲激响应到加窗截断

2.1 理想带阻的冲激响应:全通减带通

窗函数法的第一步,是把频率域指标翻译成一条理想冲激响应序列。带阻滤波器在频域上等于“全通减去带通”:直流和低频、高频全部通过,只有中间一段被扣掉。全通滤波器的冲激响应就是单位冲激 δ[n],所以理想带阻的冲激响应可以写成:

h_BR[n] = δ[n] − h_BP[n]

其中 h_BP[n] 是理想带通的冲激响应。带通的推导从低通出发:截止频率为 fc 的理想低通,其冲激响应为 2fc·sinc(2fc·n),这里的 fc 是相对采样率的归一化频率,sinc(x) = sin(πx)/(πx)。两个低通相减就得到带通,于是带阻的完整表达式是:

h_BR[n] = δ[n] − 2fu·sinc(2fu·n) + 2fl·sinc(2fl·n)

n=0 时 sinc 取极限值 1,所以中心采样点的值是 1 − 2(fu − fl)。直接照这个公式生成系数即可,实现代码如下:

import numpy as np def bandstop_ideal_coeffs(N, fl, fu): # N: 系数个数;fl, fu: 归一化截止频率(相对采样率) m = np.arange(N) - (N - 1) / 2.0 h = np.zeros(N) zero = np.abs(m) < 1e-12 h[zero] = 1.0 - 2.0 * (fu - fl) h[~zero] = (2.0 * fl * np.sinc(2.0 * fl * m[~zero]) - 2.0 * fu * np.sinc(2.0 * fu * m[~zero])) return h

这段代码做了两件关键的事:一是把序列中心移动到 (N−1)/2 处,因为理想冲激响应是偶对称且非因果的,必须先平移才能截取;二是用 np.sinc 直接计算理想响应,避免了手动实现 sin(x)/x 在零点处的除零问题。注意 fl 和 fu 必须归一化到采样率,比如采样率 1000Hz、截止频率 42.5Hz,传入的就是 0.0425。平移量 (N−1)/2 决定最终滤波器的群延迟,这个值在实时系统中直接影响延迟预算。

2.2 加窗截断与阶数估算:过渡带宽度决定系数个数

理想冲激响应无限长,直接截断等于乘一个矩形窗,频域上会看到严重的 Gibbs 现象:阻带边缘出现约 21dB 的旁瓣,无论把系数取多长都压不下去。这就是要加窗的原因。窗函数的作用是在截断的同时,让序列两端平滑衰减到零,牺牲过渡带宽换取阻带衰减。窗函数法最核心的工程判断是:阻带衰减由窗的旁瓣水平决定,过渡带宽由窗的主瓣宽度决定,两者相互独立。想要 40dB 衰减就选汉宁窗,想要 53dB 就选海明窗,想要 74dB 就选布莱克曼窗,对应的过渡带会依次变宽。

阶数 N 的估算公式由过渡带宽和窗函数决定:Δf ≈ c / N,其中 Δf 是过渡带宽(Hz),c 是窗函数对应的近似常数。换算得到 N ≈ c · fs / Δf。不同窗函数的 c 值差异很大,工程上常用的近似值如下表:

窗函数过渡带宽常数 c阻带衰减(典型)适用场景
矩形0.9~21 dB只做截断,不推荐用于带阻
汉宁3.1~44 dB通用,旁瓣衰减快
海明3.3~53 dB过渡带和衰减兼顾,最常用
布莱克曼5.5~74 dB需要深衰减,接受更高阶数
凯泽(A−7.95)/14.4连续可调阻带衰减有硬性指标时

表中 c 值是无量纲近似常数,乘上 fs 再除以目标过渡带宽就是所需阶数。

2.2.1 设计指标到设计参数的具体换算

以去除心电信号中的 50Hz 工频为例。假设采样率 fs = 1000Hz,要求阻带覆盖 45~55Hz,过渡带从 40Hz 到 45Hz、55Hz 到 60Hz,每侧过渡带宽 Δf = 5Hz。选用海明窗,阶数估算为 N ≈ 3.3 × 1000 / 5 = 660,取奇数 661。这里有个新手容易踩的坑:设计截止频率不要填阻带边缘 45Hz 和 55Hz,而要填过渡带的中点,即 42.5Hz 和 57.5Hz。因为加窗后滤波器频响的 −6dB 点大约落在理想截止频率处,过渡带会以这个点为中心向两侧展开。如果直接填阻带边缘,实际过渡带会整体向中心偏移,导致 45Hz 处的衰减不够,阻带指标超标。

提示:fir1、firwin 等设计工具的 cutoff 参数传入的也是过渡带中点,不是阻带边界。把阻带边缘当成截止频率传入,是窗函数法设计里最常见的指标偏移原因。

3. 用 Python 从零实现基于窗函数法的 FIR 带阻滤波器

3.1 完整设计函数:理想系数、加窗、归一化一次完成

上一章的理想系数生成只是半成品。把窗函数乘法、直流增益归一化整合到一个函数里,才是软件实现的标准形态。下面是完整可用的设计代码:

import numpy as np def window_coeffs(name, N): n = np.arange(N) M = N - 1 if name == 'hamming': return 0.54 - 0.46 * np.cos(2.0 * np.pi * n / M) if name == 'hanning': return 0.5 - 0.5 * np.cos(2.0 * np.pi * n / M) if name == 'blackman': return (0.42 - 0.5 * np.cos(2.0 * np.pi * n / M) + 0.08 * np.cos(4.0 * np.pi * n / M)) raise ValueError('unsupported window: %s' % name) def design_bandstop_fir(N, fc1, fc2, fs, win='hamming'): M = N - 1 m = np.arange(N) - M / 2.0 f1 = fc1 / fs f2 = fc2 / fs h = np.zeros(N) zero = np.abs(m) < 1e-12 h[zero] = 1.0 - 2.0 * (f2 - f1) h[~zero] = (2.0 * f1 * np.sinc(2.0 * f1 * m[~zero]) - 2.0 * f2 * np.sinc(2.0 * f2 * m[~zero])) h *= window_coeffs(win, N) h /= np.sum(h) return h

这个函数的参数设计遵循了工程惯例:N 是系数个数,不是滤波器阶数。严格说 FIR 的阶数是 N−1,但大多数函数库和文档里直接用 numtaps 代表系数个数,接口上保持一致可少踩一个坑。fc1 和 fc2 是过渡带中点频率,单位 Hz,在函数内部归一化。h /= np.sum(h) 这一行强制直流增益为 1,理论上带阻通带增益本身就是 1,但浮点累加和高阶数下数值误差会累积,显式归一化能保证通带增益不偏。窗函数单独抽成 window_coeffs,方便切换窗型做对比实验。

调用方式如下:

fs = 1000.0 # 采样率 1000 Hz N = 661 # 系数个数取奇数 fc1, fc2 = 42.5, 57.5 # 过渡带中点 h = design_bandstop_fir(N, fc1, fc2, fs, 'hamming') print("系数个数:", len(h)) print("关于中心对称:", np.allclose(h, h[::-1]))

系数个数取奇数不是习惯问题,而是线性相位 FIR 类型选择的问题。奇数 N 对应 Type I 滤波器,频率响应在奈奎斯特频率处没有约束;偶数 N 对应 Type II 滤波器,频响在 fs/2 处天然为零。如果带阻的上通带延伸到接近 fs/2,Type II 会把高频通带压出一个凹坑,指标直接作废。带阻滤波器统一用 Type I,即 N 取奇数。

3.2 用 lfilter 做实时滤波:初始条件决定瞬态长度

系数设计好之后,滤波本身是一个卷积运算。逐点卷积写法直观但效率低,工程上使用直接 I 型差分方程实现,即 scipy.signal.lfilter。很多人直接调用 lfilter 却忽略初始条件,导致输出开头多出一段建立过程,这在离线处理里看不出问题,在实时系统里会表现为前几十个采样点幅度异常。正确的做法是用 lfilter_zi 初始化状态:

from scipy.signal import lfilter, lfilter_zi def apply_fir(x, h): zi = lfilter_zi(h, 1.0) * x[0] # 按首采样点缩放初始状态 y, _ = lfilter(h, 1.0, x, zi=zi) return y

lfilter 的第二个参数是分母系数,FIR 滤波器分母为 1.0,所以直接传标量。zi 是滤波器内部延迟线的初始状态,乘以 x[0] 表示假设滤波前信号已经稳定在第一个采样值,这是“稳态启动”的近似,能显著缩短瞬态。如果环境不允许预填充,也可以接受输出前 (N−1)/2 个点作废,群延迟是 330 个采样点,这点事先算进延迟预算即可。

3.3 与 scipy.signal.firwin 交叉验证

自己写的设计函数需要对照验证,最直接的方法是跟 scipy.signal.firwin 的结果做差。firwin 同样是窗函数法,内部用频率采样构造理想滤波器,理论上与我们的解析公式等价:

from scipy.signal import firwin h_ref = firwin(N, [fc1, fc2], pass_zero=True, fs=fs, window='hamming') print(np.max(np.abs(h - h_ref)))

pass_zero=True 表示直流在通带内,对应带阻;如果写 False 则变成带通。两者的系数差通常在 1e-12 量级,只有浮点舍入差异。若差异很大,优先检查 fc1、fc2 是否归一化,以及窗函数实现里是否把 N 和 N−1 弄混。窗函数的分母必须是 N−1,工程上写成 N 会导致端点不为零,阻带衰减直接掉十几 dB。

4. 频响验证与滤波效果评估:从设计到上信号的完整流程

4.1 用 freqz 检查阻带深度与通带纹波

系数设计出来不等于指标达标,一定要先看频率响应再上信号。scipy.signal.freqz 返回离散时间傅里叶变换的采样点,配合 fs 参数可以直接用 Hz 读坐标:

from scipy.signal import freqz w, H = freqz(h, worN=8192, fs=fs) mag_db = 20.0 * np.log10(np.maximum(np.abs(H), 1e-12)) sb = (w >= 45.0) & (w <= 55.0) pb = ((w >= 0.0) & (w <= 40.0)) | ((w >= 60.0) & (w <= 500.0)) print("阻带最大增益: %.2f dB" % np.max(mag_db[sb])) print("通带最大纹波: %.3f dB" % np.max(np.abs(mag_db[pb])))

worN=8192 是频响采样点数,点数越多频率轴越细,阻带边缘处的极值越容易被捕捉到。检查阻带区间时不要只看单个频点,要用区间最大值,因为窗函数法的阻带纹波不是单调的,45Hz 处达标不代表 47Hz 处也达标。在本例参数下,海明窗设计的滤波器阻带衰减通常能压到 50dB 以上,通带纹波在 0.1dB 以内。

4.2 构造测试信号,定量测量 50Hz 的衰减量

频响曲线只能说明线性系统特性,最终要确认的是真实信号经过滤波后的效果。构造一个三段叠加的测试信号:5Hz 有效信号、50Hz 工频干扰、120Hz 高频噪声,幅度比例模仿真实场景:

import numpy as np t = np.arange(0, 2, 1 / fs) x = (1.0 * np.sin(2 * np.pi * 5 * t) + 0.5 * np.sin(2 * np.pi * 50 * t) + 0.1 * np.sin(2 * np.pi * 120 * t)) y = apply_fir(x, h) def amp_at_freq(sig, freq): X = np.fft.rfft(sig) / len(sig) f = np.fft.rfftfreq(len(sig), 1 / fs) idx = np.argmin(np.abs(f - freq)) return 2.0 * np.abs(X[idx]) for f in (5, 50, 120): print("%3d Hz 幅度: 滤波前 %.3f -> 滤波后 %.3f" % (f, amp_at_freq(x, f), amp_at_freq(y, f)))

FFT 幅度谱里每个 bin 的频率间隔是 fs / len = 0.5Hz,50Hz 正好落在整数 bin 上,可以直接读幅度。如果干扰频率不是整数 Hz,要用 argmin 找到最近的 bin,但会引入频谱泄漏误差,更严谨的做法是加窗或用 Goertzel 算法精确测量单频幅度。滤波后 50Hz 分量应衰减 50dB 量级,5Hz 和 120Hz 分量的幅度基本不变,滤波前后相位差正好对应群延迟 330 个采样点。

4.3 群延迟与通带边缘振荡的处理

线性相位 FIR 的群延迟恒定为 (N−1)/2,本例是 330 个采样点,在 1000Hz 采样率下就是 0.33 秒。这个数字对离线分析毫无影响,但对实时控制系统可能不可接受。如果延迟超预算,不要急着降阶数,先确认是不是采样率取得过高——把采样率从 1000Hz 降到 250Hz,同样的 5Hz 过渡带,阶数从 661 降到 165,延迟同步缩到 0.33 秒以内。

注意:filtfilt 零相位滤波会让延迟翻倍,因为正向反向各过一次。离线分析用它没问题,实时系统里必须用 lfilter 配合精确的延迟补偿。

通带边缘的振荡是窗函数法固有的,表现为 40Hz 和 60Hz 附近的微小幅值起伏,属于设计指标的一部分,不是数值 bug。振荡幅度由窗的旁瓣水平决定,换成布莱克曼窗可以压得更深,但过渡带会从 5Hz 拓宽到约 8Hz,阻带边缘就要相应外扩。振荡如果出现在敏感频段,说明指标给得太紧,优先调整通带边缘频率而不是修改窗函数。

5. 工程落地的三个进阶操作:多带阻、等波纹与定点导出

5.1 单带阻扩展到多带阻

很多时候干扰不止一个频点,比如 50Hz 工频和 100Hz 二次谐波要同时滤除。窗函数法做多带阻几乎零成本:全通减去多个带通之和,系数做一个循环累加即可。复用前面 design_bandstop_fir 的理念,把每条阻带的带通响应算出来累加,再用同一个窗加窗:

def design_multi_bandstop_fir(N, bands, fs, win='hamming'): M = N - 1 m = np.arange(N) - M / 2.0 zero = np.abs(m) < 1e-12 h = np.zeros(N) h[zero] = 1.0 for fl, fh in bands: f1, f2 = fl / fs, fh / fs bp = np.zeros(N) bp[zero] = 2.0 * (f2 - f1) bp[~zero] = (2.0 * f2 * np.sinc(2.0 * f2 * m[~zero]) - 2.0 * f1 * np.sinc(2.0 * f1 * m[~zero])) h -= bp h *= window_coeffs(win, N) h /= np.sum(h) return h

多带阻的阶数以最窄那条阻带的过渡带为准,其他阻带即使设得宽一点也不会增加总阶数。但要注意,多个阻带累加后通带纹波会叠加,海明窗下两条阻带的通带纹波比单条多出约 1 倍,设计指标要留余量。

5.2 等波纹设计会在什么时候替代窗函数法

窗函数法的频率响应误差是均匀分布在通带和阻带的,无法针对某个频段局部优化。当阻带要求超过 60dB、或过渡带窄到窗函数法需要上千阶时,就该换 Parks-McClellan 等波纹算法,在 scipy 里对应 scipy.signal.remez。等波纹设计带阻时需要手动指定通带权重和阻带权重,把衰减压力分配到指标更宽的频段上,同样 45~55Hz 阻带、40/60Hz 通带,等波纹方案往往能省下 20%~30% 的阶数。代价是设计复杂度明显上升,权重的反复调整是常态。

5.3 从浮点系数到 FIR Compiler:定点化和多相落地

软件仿真通过后,滤波器最终可能落到 FPGA 平台。浮点系数要先量化为定点,常见做法是转成 Q15 格式,即 16 位有符号数,范围 −1 到 1−2⁻¹⁵:

coeff_q15 = np.round(h * 32768.0).astype(np.int16) print("最大绝对值:", np.max(np.abs(coeff_q15))) print("系数和:", coeff_q15.sum())

定点化后的系数和通常不再是 2 的整数次幂,直流增益会偏移万分之几,在 FIR Compiler IP 核里要勾选重新归一化,或在软件端预先补偿。把 coeff_q15 写入 Xilinx 的 .coe 文件格式后,可在 Vivado 的 FIR Compiler 里加载,滤波器类型选 Bandstop,系数宽度匹配 Q15。针对高采样率场景,FIR Compiler 支持自动多相分解,把一次长卷积拆成多路并行短卷积,降低时钟频率,系数本身不需要手工拆相位,IP 核会按抽取率或通道数自动重排。验证时先用仿真激励对比 IP 输出与 Python 的 lfilter 结果,逐拍对齐后再上板,能省大量排查时间。

本文还有配套的精品资源,点击获取

需要专业的网站建设服务?

联系我们获取免费的网站建设咨询和方案报价,让我们帮助您实现业务目标

立即咨询