简介：这份资源面向学习数字信号处理、需要在MATLAB中实现FIR带阻滤波器的学生与工程人员，围绕长度N=45、阻带衰减AS=60dB的设计目标，给出凯塞-贝塞尔窗函数法的完整实现思路。压缩包内仅含1个doc文档，约60KB，以文字与源程序代码为主，便于直接阅读和复制调试。文档重点讲解窗函数参数beta对主瓣宽度、旁瓣大小与过渡带宽度的影响，并给出Beta=0.1102*(As-8.7)的计算依据；同时提供freqz.m子程序，用于求取相对振幅、绝对振幅、相位响应与延时群，配合ideal_lp函数生成理想低通响应，再与凯塞窗相乘得到实际冲激响应。读者可据此掌握从指标确定、窗函数选择到频率响应计算与图形化验证的完整流程，并通过stem与plot对比理想、窗函数及实际冲激响应和幅度响应，理解通带纹波、阻带衰减与实现复杂度之间的权衡。目前已有1376人学习下载，适合作为课程设计或滤波器实验的参考范例。
1. 从一段被 50 Hz 工频干扰毁掉的采集信号说起
做生物电、振动或音频采集的人大多遇到过这种场景：传感器本身没问题，放大器也干净，可一接上现场电源，频谱里就冒出一根又粗又稳的 50 Hz 尖峰，连它的二三次谐波都跟着起来。低通滤波器拦不住它，因为有用信号可能就在 40 Hz 到 60 Hz 之间；高通更没用，它只会把工频完整留下。这时候真正对口的工具是带阻滤波器，而 FIR 带阻因为能做到严格线性相位、系数稳定、不怕温漂，成了很多离线分析和嵌入式实现的首选。
MATLAB 里设计 FIR 带阻滤波器，核心就三件事：把阻带指标翻译成归一化频率，选一个合适的窗函数把理想冲激响应截断，再用freqz和filtfilt验证它到底有没有把那段频率压下去。凯塞—贝塞尔窗函数（Kaiser 窗）是这里最值得掌握的窗，因为它用一个可调的 β 参数就能在旁瓣衰减和过渡带宽度之间连续折中，比汉宁、汉明这些固定窗灵活得多。这篇就按“指标怎么定、窗怎么选、代码怎么写、结果怎么验”的顺序，把一套能直接复现的流程讲清楚。
2. FIR 带阻滤波器的指标换算与窗函数选型
2.1 把 Hz 指标翻译成归一化频率和阶数
FIR 设计的第一步不是写代码，而是把工程需求写成滤波器能读懂的四个数：通带边界、阻带边界、通带波纹、阻带衰减。假设采样率fs = 1000 Hz，要压掉 48~52 Hz 的工频，同时保留 0~40 Hz 和 60~500 Hz 的信号，那么：
- 阻带：48~52 Hz
- 下过渡带：40~48 Hz，上过渡带：52~60 Hz
- 阻带衰减：至少 60 dB
MATLAB 的滤波器设计函数统一使用归一化频率，即真实频率除以奈奎斯特频率fs/2。所以 48 Hz 对应48/500 = 0.096，52 Hz 对应0.104。这一步最容易出错的地方是有人直接除以fs，结果滤波器整体偏移一倍，阻带完全打偏。
过渡带宽度决定了阶数。对凯塞窗，经验公式是先用阻带衰减A反推 β，再由过渡带宽度Δω估算阶数N：
| 参数 | 含义 | 本例取值 |
|---|---|---|
fs | 采样率 | 1000 Hz |
f_pass1 | 下通带边界 | 40 Hz |
f_stop1 | 下阻带边界 | 48 Hz |
f_stop2 | 上阻带边界 | 52 Hz |
f_pass2 | 上通带边界 | 60 Hz |
A | 阻带衰减 | 60 dB |
2.2 凯塞窗的 β 与阶数估算公式
凯塞窗的定义里，β 控制窗的形状：β 越大，主瓣越宽、旁瓣越低，也就是阻带衰减越好但过渡带越宽。经典估算式（Kaiser 本人给出的经验式）是：
- 当
A > 50时，β = 0.1102 * (A - 8.7) - 当
21 <= A <= 50时，β = 0.5842*(A-21)^0.4 + 0.07886*(A-21) - 当
A < 21时，β = 0
阶数估算用N ≈ (A - 8) / (2.285 * Δω)，其中Δω是过渡带宽度（弧度）。本例过渡带 8 Hz，归一化后Δω = 2π * 8 / 1000 ≈ 0.0503，代入得N ≈ (60-8)/(2.285*0.0503) ≈ 452。这个数字说明：想用 8 Hz 过渡带换 60 dB 衰减，阶数会到四百多，FIR 的代价就在这里。
提示：如果阶数高到实时系统扛不住，优先放宽过渡带而不是降衰减。过渡带从 8 Hz 放到 20 Hz，阶数能掉到一百多。
2.3 为什么不用fir1直接一把梭
fir1确实能一行出系数，但它默认用汉明窗，β 不可调，遇到需要精确控制阻带衰减的场合就不够用。更可控的做法是显式构造凯塞窗再乘理想冲激响应，或者用fir1的kaiser参数形式。下面这段是显式版本，便于理解每一步在干什么：
fs = 1000; % 采样率 f = [0 40 48 52 60 fs/2] / (fs/2); % 归一化频带边界 a = [1 1 0 0 1 1]; % 对应期望幅度：通-阻-通 A = 60; % 阻带衰减 dB % 由衰减反推 Kaiser 窗 beta if A > 50 beta = 0.1102 * (A - 8.7); elseif A >= 21 beta = 0.5842*(A-21)^0.4 + 0.07886*(A-21); else beta = 0; end % 过渡带宽度（弧度）与阶数估算 dw = 2*pi*(48-40)/fs; N = ceil((A - 8) / (2.285 * dw)); if mod(N,2) == 0, N = N + 1; end % 保证奇数阶，类型 I 线性相位 b = fir1(N-1, f, a, kaiser(N, beta));逻辑说明：f和a是成对的频带-幅度描述，MATLAB 会在相邻点之间做理想过渡；fir1的第一个参数是阶数（比系数个数少 1），所以传N-1；kaiser(N, beta)生成长度 N 的窗，乘上去完成截断。参数上，N取奇数是为了得到类型 I 线性相位 FIR，群延迟是整数采样点，方便后续对齐。
3. 用 freqz 和 filtfilt 验证带阻效果
3.1 频响曲线要看哪几个点
系数出来不等于设计成功，必须看频响。freqz给出幅频和相频，重点核对三处：阻带最低衰减是否达到 60 dB、两个通带边缘有没有被削、过渡带是否落在 48~52 Hz 之外。
[H, w] = freqz(b, 1, 4096, fs); % 直接给采样率，横轴就是 Hz magdB = 20*log10(abs(H) + eps); % 检查阻带内最大增益 stopBand = (w >= 48) & (w <= 52); fprintf('阻带最大增益: %.2f dB\n', max(magdB(stopBand))); % 检查通带波纹 passBand = (w <= 40) | (w >= 60); fprintf('通带最大衰减: %.2f dB\n', -min(magdB(passBand))); plot(w, magdB); grid on; xlabel('频率 (Hz)'); ylabel('幅度 (dB)'); xline(48,'r--'); xline(52,'r--');逻辑说明：freqz(b,1,4096,fs)里第三个参数是频率点数，点越多曲线越平滑；eps防止对零取对数。参数上，stopBand和passBand的逻辑索引直接对应前面定的频带，如果这里算出来的阻带增益是正的，说明频带边界写反了。
3.2 用 filtfilt 做零相位滤波并对比
离线分析里我一般用filtfilt而不是filter，因为它前后各滤一遍，把相位抵消掉，群延迟为零，波形不会整体平移。代价是等效阶数翻倍，过渡带会略陡一点。
t = 0:1/fs:2; x = sin(2*pi*10*t) + 0.8*sin(2*pi*50*t) + 0.3*sin(2*pi*120*t); y_filt = filter(b, 1, x); % 单次滤波，有延迟 y_ff = filtfilt(b, 1, x); % 零相位 % 用 FFT 看 50 Hz 分量被压了多少 X = abs(fft(x)); Y = abs(fft(y_ff)); f_axis = (0:length(x)-1)*fs/length(x); idx50 = find(f_axis >= 49 & f_axis <= 51); fprintf('滤波前 50Hz 幅值: %.3f\n', max(X(idx50))); fprintf('滤波后 50Hz 幅值: %.3f\n', max(Y(idx50)));逻辑说明：filter的输出会滞后约N/2个采样点，做时域对齐时要补偿；filtfilt不需要补偿但要求信号长度大于 3 倍滤波器阶数，短信号会报错。参数上，filtfilt默认使用反射式边界延拓，信号两端有强瞬态时可以在末尾加'padlen'控制延拓长度。
3.3 常见翻车点与排查顺序
设计不达标时，按这个顺序查：先确认归一化频率除的是fs/2不是fs；再看N是不是奇数、fir1传的是不是N-1；然后核对f和a的长度是否一致、是否单调递增；最后检查kaiser的 β 有没有算错。多数“阻带压不下去”的问题，根源都在频带边界和阶数估算这两步。
4. 从凯塞窗到等波纹：进阶设计与工程落地技巧
4.1 用 firpm 换更短的长度
凯塞窗是“窗函数法”，阶数偏保守。如果对阶数敏感，可以换等波纹设计firpm（旧名remez），它在同样指标下通常能省 20%~30% 的阶数，代价是通带和阻带波纹等幅振荡，且设计可能不收敛。
% 等波纹带阻：f 与 a 含义同上，但需要指定各带权重 N_pm = 300; b_pm = firpm(N_pm, f, a, [1 10 1]); % 阻带权重给 10，压得更狠 [H2, w2] = freqz(b_pm, 1, 4096, fs); fprintf('firpm 阻带最大增益: %.2f dB\n', ... max(20*log10(abs(H2((w2>=48)&(w2<=52))) + eps)));逻辑说明：firpm第四个参数是权重向量，长度等于频带数的一半，给阻带更大权重会让它优先满足阻带衰减。参数上，N_pm必须偶数才能得到类型 II 线性相位，若要求奇数阶就传奇数并接受类型 I。和凯塞窗版本对比同一阻带增益下的阶数，就能判断哪种更划算。
4.2 定点化与嵌入式移植的注意点
FIR 系数最终要落到 MCU 或 FPGA 上时，浮点转定点是绕不开的一步。常见做法是先把系数归一化到最大绝对值 1，再乘2^Q取整，用 Q15 格式存储：
Q = 15; b_norm = b / max(abs(b)); b_fix = round(b_norm * (2^Q - 1)); b_fix(b_fix > 2^Q-1) = 2^Q-1; % 饱和处理 b_fix(b_fix < -2^Q) = -2^Q; fprintf('系数动态范围: %d ~ %d\n', min(b_fix), max(b_fix));逻辑说明：归一化保证不溢出，round后做饱和截断防止边界回绕。参数上，Q 越大精度越高但累加器位宽要求越高，Q15 配 32 位累加器是常见组合。移植后务必用同一段测试信号在 MATLAB 和硬件上各跑一遍，逐点比对输出，定点误差通常体现在阻带底部抬升几个 dB。
4.3 验证清单与参数速查
落地前过一遍这张表，能挡掉大部分返工：
| 检查项 | 期望 | 不达标时改什么 |
|---|---|---|
| 阻带最大增益 | ≤ -60 dB | 增大 N 或 β，或换 firpm |
| 通带最大衰减 | ≤ 0.5 dB | 减小 β，放宽过渡带 |
| 群延迟 | 约 N/2 采样点 | 用 filtfilt 消除 |
| 系数和 | 接近 0（带阻） | 检查频带边界 |
| 定点阻带抬升 | < 3 dB | 提高 Q 或加宽累加器 |
最后给一个实用技巧：设计完先把b存成.mat或文本，连同fs、N、beta一起记录，下次换采样率时按比例缩放频带边界即可复用整套流程，不必从头推公式。
本文还有配套的精品资源，点击获取
