简介在数字信号处理中低通滤波是去除高频噪声、保留低频分量的常用手段广泛应用于通信与音频降噪等场景。频率取样法为FIR滤波器设计提供了直观高效的路径。这份MATLAB源码完整演示了基于频率取样法设计并实现FIR低通数字滤波器的过程面向初学者、课程设计者及需要快速验证滤波算法的工程师。压缩包体积仅2KB包含1个m脚本文件结构紧凑便于阅读与修改参数。目前已有947人学习浏览说明这种小型可运行示例对动手理解滤波设计有实际帮助。脚本覆盖了指标设定、频率点采样、逆离散傅里叶变换求冲激响应以及系数调整等关键环节运行后可直观看到滤波器频率响应和输入信号低通滤波后的时域/频域变化通过修改变量即可调整截止频率快速适配不同需求有助于打通DTFT理论与工程实现之间的联系。1. 频率取样法设计 FIR 低通滤波器为什么直接在频域打点反而最容易翻车频率取样法可能是 FIR 低通滤波器设计里看起来最“直给”的方法把理想低通的频率响应按等间隔取 N 个点然后做一次 IDFT滤波器系数就出来了。但实际操作过的人都知道这里全是坑——取到的频率点没有控制好过渡带阻带衰减就只有 -20dB采样点虚部没清零ifft 出来的系数直接带复数根本没法用。这份 DTFT1.m 代码就是把这套流程完整跑了一遍从频率取样、IDFT 得到系数到用 filter 对输入信号做低通滤波并画图验证。对正在做数字信号处理课程设计、或者想手写 FIR 系数替代 fir1 的人这份脚本能直接当作入门模板。本文把每一步拆开讲清楚原理、代码、验证以及我踩过的五个坑。2. 频率取样法的原理N 个频域点怎么变成一组 FIR 系数2.1 从理想低通到可实现的 FIR为什么频率取样本身就能给出系数FIR 滤波器的长度为 M系数个数为 M其频率响应 H(e^{jω}) 在单位圆上是以 2π 为周期的连续函数。频率取样法的做法是把单位圆上的频率响应等间隔地取 N 个点每个点对应一个 H(k)然后通过 IDFT 得到 N 点冲激响应 h(n)。这一步的数学依据就是 IDFT 本身就是 DTFT 的一种内插重建DTFTH(e^{jω}) Σ h(n) e^{-jωn} IDFTh(n) (1/N) Σ H(k) e^{j 2πkn/N}当 N 足够大、且采样点设置合理时h(n) 就是目标 FIR 滤波器的系数。这就是为什么 DTFT1.m 里设计滤波器这一步并不需要调用 fir1 之类的高级函数——频率取样法的思路就是「先画理想频率响应再通过逆变换把频率响应翻译成时域系数」。这种设计方法特别适合教学场景因为每一步都能直观地对应到频域图上。实现时第一步是把 0~fs 映射到 0~2π把频率轴等分成 N 个点记录每个点落在通带、过渡带还是阻带。常见的设置是通带内取样值为 1阻带内取样值为 0过渡带根据宽度插入一两个中值点。这个中值点就是整个设计的“玄学”所在它决定了旁瓣到底能压到多少。2.2 线性相位约束取样点为什么不能全用实数值FIR 滤波器要做到线性相位冲激响应必须满足偶对称或奇对称条件而频率取样法设计时这个对称条件会映射成频域取样值 H(k) 的相位约束。简单说当 N 为偶数时k 0 和 k N/2 这两个点必须为实数其他点在构造时应按共轭对称的方式设置这样 IDFT 出来的 h(n) 才是实系数。很多人在这一步翻车直接把理想低通幅响应全是 0 或 1当作 H(k)然后调用 ifft得到的是复数系数。原因就是少了相位项 e^{-jπk(N-1)/N}。常见做法是给每个取样点乘上这个线性相位因子保证设计出的冲激响应是因果的、对称的。注意相位从 0 开始依次按固定步长增加频域取样点的虚部就自然被纳入计算。2.3 三种 FIR 设计方法选哪个频率取样法、窗函数法与等波纹法的对比做 FIR 低通滤波器设计常见的方法不止频率取样法一种。很多教材先讲窗函数法fir1再讲频率取样法最后讲最优化设计firpm。实际选型时频率取样法适合两种场景滤波器阶数不是特别高、希望直接控制某些指定频点的响应值以及用于教学演示——把设计过程可视化后学生能直观理解频率采样点与响应的关系。设计方法实现复杂度阻带衰减控制过渡带控制典型应用窗函数法fir1低可直接调用依赖窗类型可预测固定与窗函数有关常规滤波工程首选频率取样法中手动构造 H(k) 并 IDFT弱默认约 -20dB通过插入过渡带取样值手动调节教学、定制频率点的窄带滤波等波纹法firpm高需要迭代优化好可指定统一纹波最窄可逼近理论极限通信、高要求的精密滤波对比下来频率取样法的优点就是“直观”缺点也很明确对旁瓣的控制依赖过渡带取样点的插入而取样点的插入没有解析公式要靠试。DTFT1.m 选用这个方法就是为了把“从频域到系数”的流程完整展示而不是给出一个黑匣子式的最优解。3. DTFT1.m 逐段拆解从频率取样到滤波输出的完整代码路径3.1 设计参数采样频率、截止频率与点数 N 怎么选打开 DTFT1.m第一步通常是定义采样频率 fs、信号频率 f1、f2以及滤波器阶数或采样点数 N。以常用配置为例设 fs 1000通带截止频率 fc 100要求 120Hz 以上的信号被滤掉过渡带从 100Hz 到 140Hz。核心代码如下fs 1000; % 采样频率单位 Hz fc 100; % 通带截止频率 ft 140; % 过渡带右边界阻带起始 N 32; % 频率取样点数决定滤波器长度 k 0:N-1; % 频率索引 0~31 fk k * fs / N; % 每个取样点对应的实际频率 H zeros(1, N); % 初始化频域取样值 mid -1; % 过渡带取样值的占位 for i 1:N if fk(i) fc H(i) 1; % 通带取样值为 1 elseif fk(i) ft H(i) 0.5; % 过渡带取样值为 0.5可调 else H(i) 0; % 阻带取样值为 0 end end这段代码做的事是先把 0~fs 的频率轴切成 N 段然后逐点判断该频点在通带、过渡带还是阻带并赋目标值。N 的选择很关键N 越大取样点越密滤波器阶数越高过渡带越窄但计算量也增加。一般从 16 开始试不够再翻倍观察 freqz 的响应曲线决定是否继续加。mid 这个变量先占位后面如果要细调过渡带取样值直接改 0.5 为 0.5 和 0.25 串起来即可。3.2 加线性相位并执行 IDFTifft 之后为什么要 ifftshift频域取样值 H 构建完成后需要给每个点乘线性相位因子再调用 ifft 得到时域系数。这一步是频率取样法最容易出错的环节不乘相位因子直接 ifft得到的系数虽然数值上看起来是有限的但脉冲响应不是因果对称的滤波后的信号会多出莫名的相位失真。phase exp(-1j * pi * k * (N-1) / N); % 线性相位因子 H_lp H .* phase; % 带相位约束的取样值 h real(ifft(H_lp, N)); % IDFT 得到冲激响应取实部 % 常用的替代写法先用 ifftshift 再 ifft % h real(ifft(ifftshift(H)));逻辑说明ifft(H_lp) 会得到 N 个时域采样值其中绝大多数满足对称性但直接取出来的顺序是以 n0 为起点实际更常用的形式是 h 本身就应该对应于因果滤波器的系数。用 ifftshift 配合 ifft 的写法意义是先把频域数据按 0 频率居中再执行逆变换这样取出的 h 首元素对应脉冲响应 h(0)后续元素自然排列。两种写法结果一致只是中间语义不同初学时往往被这一步绕晕。如果直接用 real() 取实部说明设计者确认虚部是浮点误差而非真实的非零值。判断方式先看 sum(imag(h)).^2 是否小于 1e-10如果是虚部就是数值噪声可以放心丢弃如果虚部有明显大小说明相位因子构造错了要回头查 phase 的长度和方向。3.3 用 filter 对输入信号做低通滤波边界效应与延迟补偿系数 h 拿到之后滤波本身反而简单了。常见做法是构造一个低频加高频的测试信号然后用 filter 直接滤波。DTFT1.m 里对应的核心代码如下t (0:999) / fs; x sin(2*pi*50*t) sin(2*pi*200*t); % 50Hz 有用信号 200Hz 噪声 y filter(h, 1, x); % FIR 滤波h 为滤波器系数 % 也常用 conv 替代 % y_conv conv(x, h); % y_conv y_conv(1:length(x)); % 截断到和输入等长这里的 filter(h, 1, x) 表示分母为 1 的 FIR 结构h 直接作为分子系数执行的是差分方程 y(n) Σ h(m) * x(n-m)。因为 h 的长度为 N输出 y 和输入 x 等长但 FIR 滤波器有 N/2 个采样点的群延迟肉眼对比滤波前后波形时会看到 y 整体比 x 滞后 N/2 个点这是线性相位 FIR 的正常现象不是代码写错了。用 conv 实现时输出长度为 length(x) length(h) - 1需要手动截断到 length(x)。截断后两种方式结果一样但边界处的几个采样点会有差异这是卷积的边界效应引起的工程上通常丢弃开头和结尾各 N/2 个点再分析稳态波形。4. 滤波效果验证从幅频曲线到时域波形怎么确认这份代码真的达标4.1 freqz 幅频响应与相频响应一眼看出阻带衰减和线性相位设计完滤波器第一件事是用 freqz 看频率响应。freqz 在 MATLAB 中返回 0~π 弧度范围内的复频率响应配合 fs 参数换算成实际频率。它的本质是对滤波器系数做 FFT 后取前一半所以点数至少取 512 或 1024曲线才平滑。[Hf, f] freqz(h, 1, 1024, fs); figure; subplot(2,1,1); plot(f, 20*log10(abs(Hf))); % 幅频响应纵轴为 dB grid on; ylabel(幅度 (dB)); xlabel(频率 (Hz)); subplot(2,1,2); plot(f, unwrap(angle(Hf)) * 180 / pi); % 相频响应 grid on; ylabel(相位 (度)); xlabel(频率 (Hz));看幅频响应时重点确认三件事通带内是否接近 0dB阻带衰减是否达到预期过渡带是否落在 100~140Hz 之间。如果过渡带取样值设为 0.5阻带衰减通常在 -20dB 上下浮动若要求 -40dB 以下必须把过渡带取样点细化或增大 N。相位曲线在通带内应是一条直线这条直线说明滤波器保持了线性相位特性信号通过后各频率分量的相对时延一致波形不失真。4.2 时域对照滤波前后信号波形与频谱的对比只看频响还不够必须把实际信号喂进去。上一节构造的 x 包含 50Hz 有用分量和 200Hz 噪声分量滤波后 y 应该只保留 50Hz 的正弦。通过 FFT 对比滤波前后频谱能直观看到 200Hz 分量的衰减量。X abs(fft(x, 2048)); % 滤波前频谱 Y abs(fft(y, 2048)); % 滤波后频谱 f_axis (0:2047) / 2048 * fs; figure; plot(f_axis, 20*log10(X eps), b); hold on; plot(f_axis, 20*log10(Y eps), r); xlim([0 300]); grid on; legend(滤波前, 滤波后);这段代码把两个频谱画在同一个坐标系里。注意要加 eps 避免 log10(0) 产生 -Inf曲线出现断点。FFT 点数选 2048对应频率分辨率大约 0.49Hz足够看清 50Hz 和 200Hz 两个谱峰。如果滤波后在 200Hz 附近出现明显的残余谱峰说明阻带衰减不够如果 50Hz 幅值也比原始信号小说明通带内插损偏大需要检查通带取样值是否真的设成了 1。4.3 三个核心指标怎么量化纹波、衰减、过渡带宽判断一份频率取样法设计是否达标看的就是三个数通带纹波、阻带衰减、过渡带宽度。它们之间有强耦合关系在频率取样法里主要受 N 和过渡带取样值影响。指标定义频率取样法中的控制方式通带纹波通带内幅频响应的起伏单位 dB通带取样值必须严格为 1过渡带取样值越小通带边缘越容易下凹阻带衰减阻带内最大幅值与通带的比值单位 dB由过渡带取样点的个数和取值决定插一个 0.5 大约衰减 -20dB过渡带宽度通带边界到阻带边界的频率间隔约等于 2π/N × 过渡带取样点数N 越大过渡带越窄实际操作中通常先定过渡带宽度反推 N 的最小值然后调整过渡带取样值观察阻带衰减的变化。如果 N 已经很大但阻带衰减还是不够就要考虑在第 5 章提到的插入多个过渡带取样点的方法或者换窗函数法/等波纹法。5. 频率取样法五个常见坑现象、原因、解决办法5.1 系数带复数ifft 结果虚部不为零现象对 H 直接 ifft 后h 的虚部很大real(h) 后做滤波输出波形和预期完全对不上。原因频域取样值只设置了模值0 或 1没有附加线性相位因子IDFT 结果自然不是因果的实序列。N 为偶数时第 1 点k0和第 N/21 点没有成对的对称点必须手动保证它们是实数其他点需要满足共轭对称关系。解决按 3.2 节代码乘上 phase 因子或者使用 ifftshift 处理。另外检查 N 是否为偶数奇数 N 时对称法则不同相位因子的指数也要改成 exp(-1j2pik(N-1)/(2*N))两处对不上同样会产生复数系数。5.2 滤波输出整体延迟 N/2 点被误判为滤波失败现象把 x 和 y 直接叠加显示发现 y 波形和 x 错开了明显的一段看起来像相位错乱。原因线性相位 FIR 滤波器本身带有群延迟阶数为 N 时延迟恰好是 N/2 个采样点。这在频域上表现为相频曲线斜率为常数的直线是线性相位的自然属性不是代码错误。解决对齐时使用 y(N/21:end) 与 x(1:end-N/2) 对比或者做滤波前先把 x 补 N/2 个零。如果嫌延迟影响观察用 filtfilt 做零相位滤波但注意 filtfilt 要求信号长度至少是滤波器长度的 3 倍且会引入边界效应课程设计里一般不用。5.3 阻带衰减只有 -20dB过渡带取样值没插对现象幅频曲线显示阻带衰减约 -20dB想压到 -40dB 却怎么都压不下去增大 N 也没明显改善。原因频率取样法对阻带衰减的控制能力较弱默认只在理想响应边缘打一个点作过渡带值设 0.5 时 -20dB 左右基本是极限了。想要更深衰减需要增加过渡带的取样点数让频响从通带向阻带平滑过渡。解决最常用的是在过渡带先插入两个取样值 0.5925 和 0.1099这是等波纹设计经验值适用很广。N32 时把过渡带内的取样点从单点拆成两点阻带衰减可提升到 -40dB 附近。这个数值是工程经验值不是解析公式推导出来的换了截止频率也基本适用。5.4 filter 输出的前几个采样点异常边界效应现象滤波后波形开头几个点幅度明显异常甚至方向相反稳态后恢复正常。原因filter 初始条件默认为 0前 N-1 个输出点对应的输入还没填满整个滤波器窗口边界响应和稳态不一致。解决分析时直接丢弃前 N-1 个点或者观察稳态段。工程中这是标准做法如果必须处理短数据用 filter 的四参数版 filtic 设置初始条件但效果有限。5.5 增大 N 后阻带纹波反而更差频率取样法的吉布斯效应现象N 从 32 改成 64阻带衰减不但没好转边缘出现更密集的振荡纹波。原因频率取样法在理想频响的间断点处做等间隔采样采样点之间直接用插值重建而插值本身引入吉布斯振荡。N 越大振荡越是集中在间断点附近形成所谓的“过冲”阻带纹波峰值不一定随 N 单调改善。解决停掉一味增大 N 的惯性操作改插过渡带取样点。正确顺序是先按过渡带宽度估算 N 的最小值再用过渡带取样值去压阻带衰减最后微调通带边界观察 freqz 曲线的变化趋势。如果压纹波压不下说明该方法到了数学上的极限果断换 firpm 是明智做法不要在一个方法上硬耗。6. 进阶技巧过渡带插入多点取样把阻带衰减压到 -40dB频率取样法在基础流程跑通之后真正值得花时间的进阶操作是“多点过渡带取样”。原理很直接阻带衰减不够的根本原因是频响在通带边缘到阻带之间变化太陡插一个中值点不够平滑那就插两个。经验公式上过渡带插入 1 个点值 0.5时阻带衰减约 -20dB插入 2 个点值 0.5925、0.1099或近似 0.6、0.1时可到约 -40dB插入 3 个点可接近 -60dB。这套数值在很多教材中都能见到是经过大量仿真验证过的经验组合。% 在过渡带内额外定义取样值 % fc100, ft140, N64 时过渡带覆盖 k7~8 transition_vals [0.6 0.1]; % 两个过渡带取样值由内到外 H(7) transition_vals(1); H(8) transition_vals(2);这里的两个值不是随便拍的0.6 靠近通带0.1 靠近阻带它们让频响在过渡带内形成两段台阶近似地平滑了理想响应抑制了吉布斯振荡的峰值。改成 0.7、0.3 得到的效果是通带内插损增大、阻带衰减变化不大改成 0.5、0.1 则阻带衰减掉到 -30dB 左右所以最优组合还是有讲究的。从那以后我每次用频率取样法都会强制先跑一遍 freqz 确认阻带衰减再决定要不要插点绝不会把 N 往大了堆就完事。这个动作帮我省掉了很多“为什么滤波效果不行”的排查时间。这份 DTFT1.m 脚本本身就是很好的调试底稿下载后先从缺省参数跑一遍再逐个改过渡带取样值观察频率响应踩坑效率会低很多希望帮到你。本文还有配套的精品资源点击获取