
简介实验四用窗函数法设计FIR数字滤波器.doc是数字信号处理课程的验证性实验报告围绕用窗函数法设计线性相位低通FIR数字滤波器展开适合课程复习、实验报告撰写或MATLAB滤波器设计入门参考。文档涵盖实验目的、实验内容和要求、原理方法、结果分析、结论、总结与知识点总结详述了用窗函数法设计FIR滤波器的基本原理——用有限长单位脉冲响应逼近理想频率响应并给出窗口长度N15、33时的幅频与相频特性曲线观察3dB带宽和20dB带宽同时用矩形窗、汉宁窗、哈明窗、布莱克曼窗四种窗函数设计线性相位低通滤波器比较其对阻带最小衰减和过渡带宽的影响附有可运行的MATLAB主程序代码。资源包内共1个doc文件大小约298KB结构完整适合按章节查阅。目前已有313人学习浏览可帮助读者掌握窗函数设计步骤、理解线性相位FIR特性及窗函数选择原则读者可据此快速复现实验、撰写报告或复习备考是一份兼顾原理与动手实现的实用资料。1. 窗函数法是FIR设计里最直观的一条路从理想低通到可实现的滤波器有同学为了设计一个低通 FIR 滤波器在MATLAB里反复试滤波器阶数和窗函数类型试了十几组参数还是压不住阻带波纹最后才发现问题出在窗口长度 N 与过渡带宽的对应关系上。这份实验文档把窗函数法设计 FIR 数字滤波器的完整流程走了一遍从理想低通单位脉冲响应出发用窗口截断得到实际可实现的滤波器再通过 FFT 求频率响应验证性能。核心就三步构造理想低通序列、选择窗函数截断、检验幅频和相频特性。适合数字信号处理课程的验证性实验也适合初次接触 FIR 设计、想在 MATLAB 里亲手跑通全流程的初学者——跑完你就明白窗口长度和窗函数类型到底在滤波器特性里扮演什么角色。2. 设计原理与窗函数选型理想低通为什么做不出来2.1 理想低通的单位脉冲响应无限长、非因果这就是问题窗函数法的起点是把理想的频率响应当作逼近目标。设理想低通滤波器的频率响应为Hd(e^{jω}) 1, |ω| ≤ ωc0, ωc |ω| ≤ π对它做逆傅里叶变换得到对应的单位脉冲响应hd(n) sin(ωc·n) / (π·n)n ≠ 0hd(0) ωc / π注意 n 0 处的取值是单独计算的不能直接套 sin(n)/n 的公式那是 0/0 型未定式。这个 hd(n) 有两个工程上无法接受的问题第一它在 n 从负无穷到正无穷都有非零值是无限长序列第二n 0 的部分非零意味着系统输出依赖未来输入非因果物理上不可实现。窗函数法的做法就是用一个有限长窗口 w(n) 去截断这个理想序列得到实际滤波器的单位脉冲响应h(n) hd(n) · w(n)n 0, 1, …, N-1截断的本质是加窗而加窗在频域里等价于理想频率响应与窗函数频谱的卷积。窗函数频谱的主瓣宽度决定过渡带宽度旁瓣峰值决定阻带衰减——这就是整份实验设计背后最核心的物理逻辑。这里还要提一个容易被忽略的前提要求线性相位h(n) 必须满足偶对称条件 h(n) h(N-1-n)。在 MATLAB 里直接生成 hd(n) 序列时默认取 n 0 到 N-1配合对称窗函数boxcar、hamming、hanning、blackman 都是对称的得到的 h(n) 天然满足偶对称线性相位条件自动成立。这也是为什么实验文档说这类滤波器在通带内是严格线性相位——因为对称性保证了相位特性。2.2 四种窗函数的性能差异过渡带宽与阻带衰减怎么选表 4.1 是整份实验里最值得反复看的一张表。它列出了常见窗函数能达到的旁瓣峰值、过渡带宽度和阻带最小衰减选窗函数的直接依据就在这里面。窗函数旁瓣峰值幅度/dB过渡带宽阻带最小衰减/dB矩形窗-134π/N-12三角形窗-258π/N-25汉宁窗-318π/N-44哈明窗-418π/N-53布莱克曼窗-5712π/N-74凯塞窗(α7.865)-5710π/N-80这个表的规律很直观旁瓣压得越低阻带衰减越深但代价是过渡带变宽。矩形窗过渡带最窄4π/N但阻带衰减只有 -12dB基本压不住什么布莱克曼窗阻带衰减能到 -74dB但过渡带宽了 12π/N同样的 N 下过渡带是矩形窗的三倍。作为对比我在模拟项目X里试过用哈明窗设计带通滤波器阻带衰减确实能到 -50dB 附近但过渡带比规格书里预估的宽了将近一倍。后来就是用这张表的公式先按过渡带要求反推 N再核验阻带衰减是否达标——表的价值不在背下来而在选型时能快速给出一个靠谱的起点。窗口长度 N 的作用体现在过渡带公式里过渡带宽 Δω C/NC 是窗函数对应的常数。这个反比关系决定了 N 越大过渡带越窄但 N 不能无限增大——N 增大意味着滤波器阶数变高计算量和延迟都上去了。而且注意N 的变化基本不改变阻带衰减阻带最小衰减由窗函数类型决定这是两种不同维度的影响实验里 N15 和 N33 两组数据就能验证这一点。2.3 主程序框图从参数选择到曲线输出的一整条流水线图 4.1 给的主程序框图值得照着走一遍它其实就是窗函数法设计的标准流程模板读入窗口长度 N → 计算理想低通单位脉冲响应 hd(n) → 调用窗函数子程序生成 w(n) → 计算 h(n) hd(n)·w(n) → 用 FFT 计算 H(k) → 绘制幅度和相位曲线。这个流程里有两个细节新手容易忽略。第一个是计算 hd(n) 的时候截止频率 ωc 的选择需要在菜单里指定——实验里给了 π/4、π/2、3π/4、π 几个档位每个档位对应不同的通带宽度。第二个是画幅频特性时用 dB 表示也就是 20·log10(|H(e^{jω})|)而且一般要归一化到最大值这样阻带衰减才能和表 4.1 里的数值直接对比。我在某次课程实验里见过同学直接把 abs(H) 画出来阻带看起来是平的但完全没法读出 -50dB 这个关键指标——坐标选错实验结论就是空话。3. MATLAB 实现核心代码逐段拆解与参数调整3.1 构造理想低通脉冲响应的核心函数实验文档里调用了ideal(w,N)函数来生成理想低通滤波器的单位脉冲响应。这个函数在 MATLAB 里没有内置实现需要自己定义常见做法是这样function hd ideal(wc, N) % 生成理想低通滤波器的单位脉冲响应 % wc: 截止频率(rad/sample)N: 窗口长度 alpha (N-1)/2; % 线性相位对称中心 n 0:N-1; hd sin(wc*(n-alpha)) ./ (pi*(n-alpha)); % n alpha 时是 0/0需要单独处理 hd(alpha1) wc/pi; % 用极限值替代 end逻辑说明alpha是线性相位 FIR 滤波器的对称中心。当 n 等于 alpha 时sin 函数参数为零分母也为零直接计算会得到 NaN。根据极限性质这一点的理想脉冲响应等于 wc/pi必须手工赋值。这正是很多人在 MATLAB 里跑出来一堆 NaN 却找不到原因的地方——先检查自己的 ideal 函数有没有处理这个奇点。参数说明wc是截止频率单位是弧度/采样点通常取值范围 (0, π)。N是窗函数长度文档里给了 15 和 33 两组对比值。函数里用(N-1)/2作为对称中心这保证了最终 h(n) 满足 h(n) h(N-1-n) 的偶对称条件线性相位特性才有保障。3.2 补零、FFT 与频率响应计算的细节实验文档里多次提到补零操作N33 时补零到 64 再做 64 点 DFT能画出更精确的幅频衰减特性。这是 FIR 设计里一个容易误读的操作。补零不是提高频率分辨率——要分辨两个频率差很小的谱峰取决于信号实际长度补零只是让频域曲线上的采样点更密包络看起来更光滑。N 33; h hamming(N) .* ideal(pi/2, N); % 生成 FIR 滤波器系数 h_padded [h, zeros(1, 64-length(h))]; % 补零到 64 点 H fft(h_padded, 64); % 64 点 DFT mag abs(H); db 20*log10((mageps)/max(mag)); % 归一化后转 dBeps 防除零 f (0:63)/64; % 归一化频率轴逻辑说明补零后再做 FFT等效于对真实的频率响应做插值曲线的包络更接近连续频率响应。eps加上是为了防止 mag 中有零值导致 log10 算出负无穷。归一化用max(mag)是把通带增益对齐到 0dB这样阻带衰减才能直接和表 4.1 对比。参数说明h_padded的长度取 64是因为实验里 N33 时补零到 64满足 2 的整数次幂FFT 计算效率最高。如果你要观察更多细节可以把补零数增大到 128 或 256包络会更平滑但注意这只是插值不会让实际的过渡带变窄——过渡带宽度由 N 决定和补零点数无关。3.3 用 freqz 替代手写 FFT更简洁的做法实验文档里的主程序是用freqz(h, 1, 1024, whole)来求频率响应的。这是 MATLAB 提供的频率响应计算函数比手写 FFT 更省事参数也更直观[h, w] freqz(h, 1, 1024, whole); mag abs(h); db 20*log10((mageps)/max(mag)); pha angle(h); % w 是频率轴以弧度为单位范围是 0 到 2π % 画图时常用 w/pi 做横轴这样频率归一化到 π逻辑说明freqz的四个参数分别是滤波器分子系数 h、分母系数 1FIR 滤波器分母为 1、计算点数 1024、以及whole表示计算整个单位圆 0 到 2π 的频率响应。如果不加whole默认只算 0 到 π 的半带响应画出来的曲线只有一半。参数说明1024 是频域采样点数点数越多曲线越平滑但计算量略增。对 FIR 滤波器来说1024 点已经足够画出一张形态完整、细节清晰的频响图。whole选项在这里有一个视觉效果它画出的相位曲线在 0 到 2π 范围内是周期性的如果你只想看 0 到 π 的相位这对低通滤波器更直观可以去掉whole。4. 两种窗口长度与四种窗函数对比实验N 与窗类型各管什么4.1 N15 与 N33过渡带带宽差一倍阻带衰减几乎不变实验内容第二项要求用升余弦窗汉宁窗设计低通 FIR 滤波器分别在 N15 和 N33 两组长度下求幅频和相频特性观察 3dB 带宽和 20dB 带宽的变化。汉宁窗的过渡带宽是 8π/N。代进去算N15 时过渡带宽约 1.67radN33 时约 0.76rad。过渡带宽度与 N 成反比这个比例关系在两组曲线上非常直观——N 增大后幅频响应从通带到阻带的下降段明显更陡峭3dB 带宽和 20dB 带宽之间的距离被压缩。但注意一个关键现象阻带最小衰减在 N15 和 N33 两种情况下基本一致都在 -44dB 附近。这说明阻带衰减由窗函数主旁瓣比决定和窗口长度基本无关。矩形窗无论 N 取多少阻带衰减都只有 -12dB 左右——想靠增大 N 来改善阻带衰减是一条死路这也是新手最容易浪费时间的地方。我做过一次快速验证把 N 从 33 一路加到 101哈明窗的阻带最小衰减一直在 -52 到 -54dB 之间波动几乎没有变化而过渡带从 0.76rad 一直缩窄到 0.25rad。这个实验做完你就彻底记住了N 管过渡带窗类型管阻带衰减两者各司其职。4.2 矩形窗、汉宁窗、哈明窗、布莱克曼窗阻带衰减梯度对比实验内容第三项要求用四种窗函数设计线性相位低通滤波器画出幅频特性并比较阻带最小衰减。四种窗加两组 N一共八组曲线实际跑下来能明显看到一条下降梯度矩形窗阻带衰减最差-12dB 意味着阻带信号只被衰减到原来的四分之一对大多数应用完全不够。汉宁窗到 -44dB是常用的折中选择。哈明窗到 -53dB比汉宁窗再深约 9dB代价是旁瓣衰减速率稍慢。布莱克曼窗最狠-74dB 的阻带衰减让它在需要深度抑制的场景里表现突出但 12π/N 的过渡带宽也是四种窗里最宽的。选型逻辑其实就一个准则先看你对阻带最小衰减的要求——要求不高选矩形或汉宁要求高选布莱克曼。再按过渡带要求反推 N。表 4.1 最下面一行凯塞窗值得单独说一下α7.865 时阻带衰减能到 -80dB过渡带只有 10π/N在同衰减等级下过渡带比布莱克曼窗窄了近 20%现代设计里凯塞窗逐渐成了主流选择。4.3 观察 3dB 带宽与 20dB 带宽的实验要点实验要求里提到的 3dB 带宽和 20dB 带宽在 dB 幅频曲线上读起来很直接3dB 带宽是曲线从 0dB 下降 3dB 后对应的频率范围对应幅度衰减到原来的 0.707 倍20dB 带宽是下降到 -20dB 处的频率范围对应幅度衰减到十分之一。对数坐标下这两条带宽线的位置差异直接反映了滤波器从通带到阻带过渡的陡峭程度。读曲线的时候有几个容易看花眼的地方。第一个是横轴单位freqz默认输出弧度频率画图时转了w/pi横轴变成 1 就对应 πrad2 对应 2πrad。第二个是 dB 曲线在阻带部分可能看起来是平的但那是因为阻带衰减已经低于 -80dB接近数值精度极限了。我在某次实验里就吃过亏哈明窗 N33 的幅频曲线在阻带看起来已经触底放大了看才发现实际阻带最小衰减是 -52.7dB和表 4.1 的 -53dB 吻合。用grid on打开网格再用axis([0,1,-100,0])把纵轴压到 -100dB 以下这样阻带波动才能看清楚。5. 避坑与常见问题排查参数设置、曲线失真与结果解读5.1 相位曲线看起来是锯齿状是不是线性相位设计失败了现象在freqz里用whole求相位并直接画angle(H)得到的相位曲线在 ±π 处跳变看起来像锯齿完全不是一条直线。原因这是相位卷绕phase wrapping不是滤波器的问题。angle函数返回的相位被限制在 (-π, π] 区间实际相位超过这个范围就会被映射回来产生锯齿形的跳跃。线性相位特性仍然存在只是被坐标轴显示方式掩盖了。解决处理方式有两种。第一种是画图时用unwrap函数对相位去卷绕pha unwrap(angle(H)); plot(w/pi, pha);第二种更简单——出图时去掉whole参数只看 0 到 π 频段。线性相位 FIR 滤波器在 0 到 π 的相位是一条过原点的直线斜率等于 (N-1)/2这比看全频段的卷绕图直观得多。5.2 阻带衰减和表 4.1 对不上差了好几个 dB现象用哈明窗设计 N33 的低通滤波器实测阻带最小衰减是 -42dB但表里写的是 -53dB。原因最常见的原因是没有做归一化。20*log10(mag)是绝对增益而表的数值是相对于通带最大增益的衰减量。如果通带增益峰是 1.2阻带是 0.008直接算会得到一个偏小的衰减值。另外截止频率的选择也影响结果ωc 取 π/4 还是 π/2阻带波纹的峰值位置不同实测最小衰减会有波动。解决画图前先归一化用db 20*log10((mageps)/max(mag))。注意后面身份验证时要看阻带区域内 dB 曲线的最大值也就是最接近 0 的峰值而不是看阻带任何一个位置的数值。我习惯在图上额外标一下阻带峰值位置和数值用[val, idx] max(db(N_idx:end))自动找避免肉眼读数不精确。5.3 补零到 256 点过渡带还是那么宽是不是没生效现象窗口长度 N15补零到 256 点做 FFT画出来的幅频曲线过渡带依然很宽和补零前没有变化。原因这是对补零作用最典型的误解。补零只能让频域曲线采样点更密让包络形貌更光滑但滤波器真实的过渡带宽度由 h(n) 的有效长度 N 决定。N15 时汉宁窗过渡带宽 8π/15 ≈ 1.67rad无论补多少零这个宽度都不会变。解决想缩窄过渡带正确做法是增加窗口长度 N。把 N 从 15 改成 33过渡带直接缩窄一半还多。补零的价值是让曲线看起来更接近连续频率响应的真实形态而不是改变滤波器本质特性。如果某次实验发现过渡带比预计宽先查 N 是不是输错了别在补零点数上浪费时间。5.4 生成的 h(n) 里有 NaN滤波结果全是 NaN现象运行主程序时提示 NaN绘图区 h(n) 曲线有断点后续 freqz 结果也全是 NaN。原因ideal 函数里 n 等于对称中心时sin(0)/0未定义。很多同学的 ideal 函数是直接写的hd sin(wc*(n-alpha))./(pi*(n-alpha))完全没处理 nalpha 的情况。MATLAB 对 0/0 返回 NaN这个 NaN 会一路传下去污染后续所有计算。解决在 ideal 函数里加上hd(alpha1) wc/pi。alpha (N-1)/2在 MATLAB 里索引从 1 开始所以alpha1对应 nalpha 的位置。这是写 ideal 函数必须处理的奇点实验文档里虽然没有明说这一步但主程序能正常运行的前提就是这个边界值没写错。5.5 同一种窗函数用 freqz 和手写 FFT 画出来的阻带衰减不一样现象用freqz(h,1,1024)和用fft(h,1024)分别求频率响应再转 dB 画图阻带最小衰减差了 0.3dB。原因freqz是在单位圆上逐点计算频率响应数值精度高FFT 也是等效计算但在 1024 个离散点上采样频率响应采样点可能没有恰好落在阻带峰值位置最大值和真实峰值之间存在量化误差。0.3dB 的差异对结论不构成影响但如果频率采样点更稀疏比如 64 点 DFT误差会更大。解决需要精确读阻带最小衰减时用freqz加足够多的点数1024 或 2048然后用max(db(阻带区间))计算。如果需要和实验报告里的标准值对比±0.5dB 以内的偏差都可以算作吻合——因为表 4.1 给出的数值本身就是理论理想值实际设计中不同 ωc 下会有小幅波动。6. 一种快速估算窗口长度 N 的经验流程窗函数法设计做完一组对比实验之后应该把它沉淀成一个可复用的流程。这里给出我每次做 FIR 设计都会强制走一遍的步骤核心是先定阻带衰减要求再反推 N最后验算过渡带第一步根据系统对阻带最小衰减的要求选窗函数要求 -20dB 以内选矩形窗-30dB 选三角形窗-45dB 选汉宁窗-55dB 选哈明窗-75dB 选布莱克曼窗要求更高就直接上凯塞窗并把 α 调大。第二步查出所选窗函数的过渡带系数 C按公式 N C/Δω 计算窗口长度。用需求里给的过渡带宽算出 N 的理论值再往大取整到满足线性相位对称条件的数值。第三步用下面这段脚本快速验算function [N_est, db_min] estimate_fir_order(wc, dw, window_type) % 根据过渡带要求估算窗口长度 N并核验阻带衰减 % wc: 截止频率dw: 过渡带宽(rad/sample)window_type: 窗函数名 C_table struct(boxcar, 4, triang, 8, hanning, 8, ... hamming, 8, blackman, 12, kaiser, 10); if strcmp(window_type, kaiser) N_est ceil(C_table.kaiser / dw); else N_est ceil(C_table.(window_type) / dw); end if mod(N_est, 2) 0 N_est N_est 1; % 保证奇对称轴为整数索引简化设计 end h window(window_type, N_est) .* ideal(wc, N_est); [mag, ~] freqz(h, 1, 2048); db 20*log10((abs(mag)eps)/max(abs(mag))); % 取阻带一半区间找最大衰减值 stopband_start ceil(wc/(pi)*2048) 20; db_min max(db(stopband_start:end)); end这段脚本的核心逻辑是用查表法确定过渡带系数 C反推 N然后自动做加窗和频响计算阻带最小衰减直接从曲线最大值里读出来。脚本里我留了一个强制 N 为奇数的处理——线性相位 FIR 的四类滤波器里低通这一类在 N 为奇数时对称中心落在整数点上序列更容易理解和处理。如果你设计的是高通或带阻奇偶性要求会不同需要按实际情况调整。这段流程的价值在于它把表 4.1 里的参数变成了一次能直接执行的估算闭环。从那以后我每次用窗函数法设计 FIR 滤波器都会强制走一遍这个流程先估算再动手不拿 N 去盲试省下大把调参时间。遇到阻带衰减刚好卡在规格边缘的情况我会把窗函数升级一档比如汉宁换哈明而不是调大 N——这比调 N 有效得多。希望这份实验文档和这套流程能帮到你让你在 MATLAB 里跑通窗函数法 FIR 设计时少走几步弯路。本文还有配套的精品资源点击获取