基于布莱克曼窗的FIR低通滤波器设计与MATLAB实现

简介：面向数字信号处理初学者和 MATLAB 使用者，这份示例通过布莱克曼窗完成 FIR 低通滤波器设计，主程序直接调用 Blackman() 生成窗函数，并结合 ideal_lp() 理想低通函数与 freqz_m() 频率响应函数，完整演示了从理想低通特性构建、加窗处理到频响验证的过程，有助于理解窗函数法设计线性相位 FIR 滤波器的核心原理。包内文件总数为 4，含 3 个 .m 文件与 1 个 .txt 文本，.m 文件分别对应主程序及滤波、频率响应等辅助函数，.txt 为简要说明，压缩包整体仅 2KB，非常轻量，便于直接阅读和运行。代码逐行添加了注释，入口函数和附带子函数均有完整源码，设计结构清晰，既能用于课程作业和实验报告，也可在此基础上修改参数完成不同指标的低通滤波器设计。已有 3502 人学习下载，对刚接触 FIR 滤波器或 MATLAB 编程的读者来说，是一份能直接运行并快速上手的参考资料。

1. 用布莱克曼窗设计 FIR 低通滤波器，先算一笔过渡带的账

FIR 低通滤波器设计里，布莱克曼窗是最容易被低估的选项：它比汉明窗多付出约七成的过渡带，换来的却是接近 74 dB 的阻带衰减。很多人在 MATLAB 里第一次用 fir1 设计低通滤波器时不理解这笔账，把阶数设成和汉明窗一样，结果过渡带比预期宽很多，就开始怀疑代码写错了。避免这个坑，需要先搞清楚窗函数法为什么非截断不可、布莱克曼窗的频谱特征和阶数估算公式，再落到 MATLAB 的 fir1 实现、freqz 验证和混合信号测试，最后补齐可复用封装与系数导出。这套路径适合刚接触滤波器的 MATLAB 用户，也适合写过 fir1 但没深究过渡带与阶数关系的工程师。

2. 窗函数法设计 FIR 低通滤波器的原理与布莱克曼窗的定位

2.1 理想低通滤波器为什么只存在于公式里

理想低通滤波器的频率响应是矩形：|ω| ≤ ωc 时增益为 1，其余为 0。对频率响应做逆傅里叶变换，得到单位脉冲响应 h_d[n] = sin(ωc·n)/(π·n)，n 覆盖负无穷到正无穷。这个序列有两个工程上无法接受的属性：不因果——n < 0 时系数不为零，意味着输出先于输入出现；不绝对可和——直接求和发散，无法稳定实现。实用的 FIR 滤波器只能是它的某种近似，最直接的做法是截取中间一段点，再右移 M 个采样获得因果系数。

直接截断在数学上相当于给理想脉冲响应乘了一个矩形窗，频域里理想低通被矩形窗频谱做卷积，产生吉布斯效应：通带边缘出现约 9% 的过冲，靠近截止频率处振铃不断，阻带旁瓣只衰减约 13 dB，而且无论取多长的截断，旁瓣峰值都不再下降，只是变得更密。这就是窗函数法存在的根本原因——用一个频谱旁瓣很低的窗函数去乘理想脉冲响应，把泄漏能量压到主瓣附近，用可控的过渡带宽度换阻带衰减。理解这一点后，布莱克曼窗的三个系数 0.42、0.5、0.08 就不再是魔法数字，而是对旁瓣水平的显式设计。

2.2 布莱克曼窗的时域表达式与频域三要素

MATLAB 的 blackman(N) 生成 N 点对称窗：

w(n) = 0.42 - 0.5·cos(2πn/(N-1)) + 0.08·cos(4πn/(N-1))，n = 0, 1, …, N-1

系数由展开后让前三阶频谱分量在窗边界近似为零反推得到。和汉明窗（0.54 - 0.46·cos）相比，布莱克曼窗多了一项 4πn/(N-1) 的余弦，主瓣宽度相应从约 8π/N 拓宽到约 12π/N，旁瓣从约 -41 dB 压到约 -58 dB。对于 FIR 低通滤波器设计，这意味着同样的阶数下，布莱克曼窗的过渡带比汉明窗宽约七成，阻带衰减却好出 15 dB 以上。这个方向上的交换是窗函数法的核心，不存在既窄过渡带又低旁瓣的自由窗，除非引入 Kaiser 窗的 β 参数做连续调节。

频域上要同时看三个量：主瓣宽度决定滤波器的过渡带宽，旁瓣峰值决定阻带能达到的最小衰减，主瓣与旁瓣的能量分配决定通带纹波。布莱克曼窗的典型值是主瓣约 12π/(N-1) 宽、旁瓣约 -58 dB。后文验证章节会逐个实测，这里先记住一个结论：过渡带和阻带衰减主要由窗型决定，窗长只负责把过渡带压得有多窄。

2.3 常用窗函数对比与阶数估算公式

FIR 低通滤波器设计里最常见的五种窗对比如下，N 为窗长，c 为过渡带系数，实际过渡带宽 Δf ≈ c·fs/N：

注意最小阻带衰减并不等于旁瓣峰值。窗函数的旁瓣是它自己的频谱指标，而滤波器阻带里的最大纹波是窗谱与理想低通卷积后的结果，两者数值上有差异，工程上直接以阻带衰减列作为验收基准。用 c 值估算阶数的公式是 N ≥ c·fs/Δf，举一个实际例子：

fs = 1000; % 采样率，单位 Hz df = 80; % 要求的过渡带宽（比如 60 Hz 到 140 Hz） c = 5.5; % Blackman 窗的过渡带系数 N = ceil(c * fs / df); % 理论最小窗长，这里为 69 if mod(N, 2) == 0, N = N + 1; end % 强制为奇数，避免半采样群延迟

这里 N = ceil(5.5×1000/80) = 69，最后一行把偶数阶修正为奇数。为什么要奇数：窗长奇数时群延迟 (N-1)/2 是整数个采样周期，多路信号做时间对齐时直接按这个整数延迟补偿即可；偶数窗长也能设计，但半采样延迟在同步场景里要额外做插值，能避免就避免。把 c 换成 3.3（汉明窗），同样过渡带只需要 42 个系数，但阻带只能到 -53 dB，这就是后文所有调参决策的出发点。

3. matlab 中 fir1 与 blackman 组合实现 FIR 低通滤波器的完整代码

3.1 从采样率、截止频率和阶数到滤波器系数的三行代码

假设测试信号采样率 1000 Hz，需要保留 100 Hz 以下的低频，过渡带要求 80 Hz 左右。按上一章的估算公式，窗长取 69：

fs = 1000; % 采样率，所有频率参数都以此为基准 fc = 100; % 截止频率，fir1 内部以此为理想低通的截止参考点 N = 69; % 窗长 = 滤波器系数个数 b = fir1(N-1, fc/(fs/2), 'low', blackman(N)); % 得到 1×69 的行向量

逐项说明。fir1 的第一个参数是滤波器阶数，等于系数个数减一，这里写 N-1；第二个参数是归一化截止频率 fc/(fs/2)，100/500 = 0.2，范围必须落在开区间 (0,1)，1 对应奈奎斯特频率；第三个参数 'low' 指定低通；第四个参数 blackman(N) 是 69 点的布莱克曼窗向量，长度必须与系数个数一致。fir1 内部做的事情就是上一章的手工公式：生成以 (N-1)/2 为对称中心的理想低通脉冲响应，逐点乘上窗向量，输出 b。设计完成的 b 可以直接交给 filter(b, 1, x) 使用。

需要留意的是 fc 是理想低通的截止频率，实际幅频响应在 fc 处的增益在 -6 dB 附近、随 N 略有偏移；如果业务指标写的是「通带必须完整覆盖到 100 Hz」，应该把 fc 的取值调高一些，多留半条过渡带的余量。对应参数的含义和易错点整理如下：

提示：fir1 的第一个参数是阶数而不是窗长。曾经有同事传 N 进去、窗用 blackman(N+1)，滤波器比预期多一个系数自己没察觉，群延迟和预期差半个采样点，多路同步时怎么对都对不上。

3.2 截止频率、窗长对过渡带和阻带衰减的影响

把三档窗长的幅频响应画在同一张图上，能直观看到阶数的真实作用：

fs = 1000; fc = 100; Ns = [31 69 151]; % 短、中、长三档窗长 figure; hold on; for k = 1:3 b = fir1(Ns(k)-1, fc/(fs/2), 'low', blackman(Ns(k))); [H, f] = freqz(b, 1, 4096, fs); plot(f, 20*log10(abs(H)), 'LineWidth', 1.2); end legend('N=31', 'N=69', 'N=151'); grid on; xline(fc, '--', 'fc'); ylim([-100 5]); xlabel('频率 (Hz)'); ylabel('幅度 (dB)');

freqz(b, 1, 4096, fs) 用 4096 点 DFT 计算幅频响应，第四个参数 fs 让频率轴以 Hz 为单位输出。图上能读出两个规律。第一，过渡带随窗长增大同步收窄，-3 dB 到 -60 dB 的宽度分别约为 177 Hz、80 Hz、36 Hz，和 5.5×fs/N 的估算值吻合。第二，三条曲线的阻带底部都压在 -70 dB 附近，说明窗长几乎不影响阻带衰减水平，阻带指标是布莱克曼窗本身决定的。这解释了为什么很多人把 N 从 31 加到 151，阻带纹丝不动，只有过渡带在变——调参方向错了，加阶数永远等不到想要的结果。

3.3 手工构造窗向量验证 fir1 内部细节

fir1 是封装好的，想确认它的内部逻辑，可以用 MATLAB 的 sinc 函数手工复现：

n = 0:N-1; mid = (N-1) / 2; hd = 2*fc/fs * sinc(2*fc/fs * (n - mid)); % 理想低通脉冲响应 w = 0.42 - 0.5*cos(2*pi*n/(N-1)) + 0.08*cos(4*pi*n/(N-1)); % 布莱克曼窗 b_manual = hd .* w; % 窗函数法核心一步 b_fir1 = fir1(N-1, fc/(fs/2), 'low', blackman(N)); max(abs(b_manual - b_fir1(:)')) % 输出约 1e-10，验证两者等价

sinc(x) = sin(πx)/(πx)，x = 0 时 MATLAB 直接返回 1，中心点不需要单独处理除零。对比结果通常落在 1e-10 量级，差异来自 fir1 内部双精度运算的舍入顺序。手工实现的意义在于确认两件事：窗函数法的本质就是「理想低通 × 有限窗」；b 没有额外的增益归一化，DC 增益约等于 1 是因为窗中心采样点权值恰好为 1，而其他点偏离 1 的部分正是通带纹波的来源。如果项目要求 DC 增益严格等于 1，可以执行 b = b / sum(b)。

4. 用 freqz 和混合信号验证 FIR 低通滤波器的真实性能

4.1 幅频、相频和群延迟的读取方法

滤波器系数拿到手，第一步是看频率响应的三个维度：

[H, f] = freqz(b, 1, 4096, fs); % 幅频与相频 gd = grpdelay(b, 1, 4096, fs); % 群延迟，单位：采样周期 figure; subplot(3,1,1); plot(f, 20*log10(abs(H))); ylim([-120 5]); grid on; xline(fc, '--', 'fc'); ylabel('幅度 (dB)'); subplot(3,1,2); plot(f, unwrap(angle(H))*180/pi); grid on; ylabel('相位 (度)'); subplot(3,1,3); plot(f, gd); grid on; ylabel('群延迟 (sample)'); xlabel('频率 (Hz)');

幅频图看三个位置：通带 0 到约 60 Hz 的纹波，布莱克曼窗通常压在 ±0.1 dB 以内；fc 处的增益，大约在 -6 dB 附近；180 Hz 以上的阻带峰值，应该低于 -70 dB。相频图为一条过原点的直线，对应恒定群延迟；群延迟图是水平线，数值等于 (N-1)/2 = 34 个采样点，在 fs = 1000 Hz 时是 34 ms。恒定群延迟意味着通带内的波形只延迟不畸变，这是线性相位 FIR 滤波器相对 IIR 滤波器的根本优势，做数据采集同步时可以直接按 34 ms 做时间补偿。

4.2 构造混合信号做端到端滤波验证

频率响应是理论值，真实数据上的表现还要用混合信号实测。构造一个同时包含通带与阻带成分的测试信号：

t = (0:4095)' / fs; x = 1.0*sin(2*pi*30*t) ... % 30 Hz，通带内 + 0.8*sin(2*pi*50*t) ... % 50 Hz，通带内 + 0.6*sin(2*pi*200*t) ... % 200 Hz，阻带内 + 0.7*sin(2*pi*350*t); % 350 Hz，阻带内 y = filter(b, 1, x); f_test = [30 50 200 350]; for k = 1:numel(f_test) wk = 2*pi*f_test(k)/fs; % Hz 转 rad/sample Hk = freqz(b, 1, [wk wk]); % 传两个相同频率点 Hk = Hk(1); fprintf('%3d Hz: %7.2f dB, 输出幅度约 %.4f\n', ... f_test(k), 20*log10(abs(Hk)), abs(Hk)); end

freqz(b, a, w) 的第三种调用形式要求 w 是至少两个元素的频率向量，标量会被解析成点数 n 而不是角频率，所以这里传 [wk wk] 再取第一个值，顺便绕开这个文档里不容易注意到的歧义。30 Hz 和 50 Hz 处增益接近 0 dB，200 Hz 和 350 Hz 处应该在 -70 dB 以下，理论上 0.6 × 10^(-70/20) ≈ 0.00019，时域上完全不可见。滤波前后波形对比时，把 y 初始的 N 个采样点去掉再画——filter 从零状态启动，前 N-1 个输出被起始瞬态污染，这是数字滤波器不可避免的边界效应。

提示：如果想验证过渡带，把 80 Hz 也加进测试信号。它落在 60~140 Hz 的过渡区之间，输出幅度会介于通带和阻带之间，而且随窗长变化，这是观察窗函数法过渡带宽最直接的办法。

4.3 三个可量化的验收指标与调参方向

三个指标互相牵制：固定窗型时，增大 N 只压窄过渡带，阻带衰减几乎不动；想同时改善两项，必须换窗型。Kaiser 窗是这里最常用的替代品，designfilt 里指定 'Window', 'kaiser' 配合 β 参数，β 越大阻带越好、过渡带越宽，本质是把布莱克曼窗的固定权衡变成连续可调。辅助验证还可以看阶跃响应：

stepz(b, 1, 256, fs); % 阶跃响应：上升沿与过渡带宽度成反比，无过冲说明阻带抑制充分

stepz(b, 1, 256, fs) 画出前 256 个采样点的阶跃响应，布莱克曼窗的阶跃响应没有过冲但上升沿比汉明窗慢，这是宽主瓣的直接表现。如果产品对时延敏感，要在设计阶段就换算清楚：34 ms 群延迟能否被系统预算接受，不能只看幅频图。

5. 把布莱克曼窗 FIR 低通滤波器封装成可复用函数并导出系数

5.1 用 designfilt 重写同一设计，把参数校验交给 MATLAB

fir1 是低层接口，参数写错时可能静默产出一组看着正常、实际偏离目标的系数。新版 Signal Processing Toolbox 推荐用 designfilt 以名称-值对表达设计意图：

d = designfilt('lowpassfir', ... 'FilterOrder', N-1, ... 'CutoffFrequency', fc, ... 'Window', 'blackman', ... 'SampleRate', fs); b = d.Coefficients;

designfilt 会直接校验 CutoffFrequency 必须小于 SampleRate/2、FilterOrder 必须为非负整数，约束违反时立即报错，不会像 fir1 那样把错误埋进系数里。返回的 d 是 DFILT 对象，filter(d, x)、freqz(d)、stepz(d) 都直接接受对象作为第一个参数，不必每次手动取出系数。如果项目里有多个不同截止频率的滤波器，把这段封装成一个函数，输入 fs、fc、df，输出 d，所有设计逻辑收拢到一个文件里，后续维护只改一处。

5.2 流式数据用 dsp.FIRFilter，别再用 filter 逐块调用

filter(b, 1, x) 每次调用都从零状态开始，适合离线处理整段数据。数据分块到达的实时系统，比如串口或 DAQ 采集，逐块调用 filter 会在每块开头产生一次瞬态。正确做法是用 dsp.FIRFilter 维护内部延迟线：

hfir = dsp.FIRFilter('Numerator', b); y1 = hfir(x(1:512)); % 第一块数据 y2 = hfir(x(513:1024)); % 内部状态自动延续，块间无缝衔接

dsp.FIRFilter 在对象内部保存滤波器状态，第二块输出的前几个采样点不再有起始瞬态。离线处理超长数组时，如果内存允许，优先一次性 filter 整段；只有数据本来就是分块到达、或者单块内存受限时，才需要切到 System object 方案。两种做法的系数 b 完全一致，区别只在状态管理方式。

5.3 导出定点系数：Q15 量化前先做归一化和频谱复测

把 b 带到没有浮点单元的嵌入式平台前，系数要量化为定点数。推荐顺序：先归一化 DC 增益，再乘 2^Q 四舍五入，最后用 freqz 复测量化后的频谱：

b = b / sum(b); % 归一化，保证 H(0) = 1 Q = 15; bq = round(b * 2^Q); % 量化到 Q15，峰值约 0.2×32768，int16 放得下 [H_float, f] = freqz(b, 1, 4096, fs); [H_q] = freqz(bq/2^Q, 1, 4096, fs); % 反缩放后对比 figure; plot(f, 20*log10(abs(H_float)), f, 20*log10(abs(H_q))); legend('浮点系数', 'Q15 系数'); grid on;

量化误差是逐点引入的，69 个系数的误差累加后，阻带通常从 -74 dB 抬升到 -65 dB 量级，通带纹波几乎不变。对比图上阻带如果被抬到 -60 dB 以上，说明该窗长的系数对量化太敏感，优先换 Q31 或直接导出双精度浮点系数；只有确定目标平台的乘累加器是 16 位时才保留 Q15。导出 C 数组用单个 fprintf 就能完成：

fprintf('%d,\n', bq); % 输出可直接粘贴进 C 的 int16_t coefs[] 数组

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