简介面向通信工程等专业数字信号处理课程设计的一份完整方案文档针对基于窗函数法的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)n0 时 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 是过渡带宽Hzc 是窗函数对应的近似常数。换算得到 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要求阻带覆盖 4555Hz过渡带从 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, winhamming): 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/2Type 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, zizi) return ylfilter 的第二个参数是分母系数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_zeroTrue, fsfs, windowhamming) print(np.max(np.abs(h - h_ref)))pass_zeroTrue 表示直流在通带内对应带阻如果写 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, worN8192, fsfs) 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])))worN8192 是频响采样点数点数越多频率轴越细阻带边缘处的极值越容易被捕捉到。检查阻带区间时不要只看单个频点要用区间最大值因为窗函数法的阻带纹波不是单调的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.5Hz50Hz 正好落在整数 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, winhamming): 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。等波纹设计带阻时需要手动指定通带权重和阻带权重把衰减压力分配到指标更宽的频段上同样 4555Hz 阻带、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 结果逐拍对齐后再上板能省大量排查时间。本文还有配套的精品资源点击获取