FIR与IIR低通滤波器在语音降噪中的Matlab实现与对比分析

做语音信号处理方向的课程设计或毕业设计时FIR和IIR数字滤波器组合起来做一套Matlab程序加设计报告是出现频率非常高的选题。这个项目表面上是“写两个滤波器滤一下语音”真做起来却发现指标怎么定、阶数怎么算、相位失真怎么解释、报告里的图和表怎么组织每一个环节都能把人卡住。我最近刚陪某高校的同学把一个类似项目完整过了一遍把带噪语音分别用FIR和IIR低通滤波做了处理对比频谱、波形和主观听感最后整理出一套可以照着复现的完整思路。这篇文章就把从WAV读取、预处理、滤波器设计到报告撰写的全过程展开讲讲顺带把那些课本上不会写、但实验里一定会遇到的坑都交代清楚。1. 项目的最初定位FIR和IIR为什么必须放一起做1.1 两种滤波器的本质差异数字滤波器按单位冲激响应的长度分成有限冲激响应FIR和无限冲激响应IIR两大类。这个分类不是纯理论游戏它直接决定了你在设计时拥有哪些自由度、要付出什么代价。FIR滤波器的输出只依赖当前和历史输入样本没有反馈回路因此冲激响应是有限长。它的突出优点是可以用对称系数实现严格的线性相位也就是说信号通过滤波器后所有频率分量的时间延迟完全相同波形不会发生相位失真。缺点是达到同样的频率选择性能时需要的阶数往往很高计算量大、延迟长。IIR滤波器的输出除了依赖输入还依赖之前的输出存在反馈项冲激响应理论上是无限长的。它最大的优势是用很低的阶数就能实现非常陡峭的过渡带计算效率高。代价是相位响应是非线性的不同频率分量通过滤波器时被延迟的程度不一样会导致波形失真。此外因为有反馈设计不当还可能不稳定。我用一个生活类比FIR像所有人必须排队依次走过安检门节奏稳定、秩序好但队伍长、耗时长IIR像有人从侧门抄近路速度快但路径绕来绕去到达时间不一致队形容易乱。1.2 语音场景下的职责划分语音信号本身有几个特点能量主要集中在300Hz到3400Hz这个窄带范围人耳对频率强度变化敏感但对相位变化的感知相对迟钝语音处理任务通常对实时性有要求比如通话系统。正因为这些特点语音信号处理场景中FIR和IIR都大量使用但侧重点不同。FIR的线性相位对波形保真要求高的场景有价值比如语音波形分析、基音检测前的预处理IIR的高效率适合实时降噪、响度控制这类对延迟敏感、人耳主要关心幅度谱的应用。课程设计里把两者放在同一个项目里其实是考察你能否针对同一组指标用两种不同的设计路线得到结果并正确解释差异。我的建议很直接项目里用相同或相近的频带指标分别设计一个FIR低通和一个IIR低通同时作用在同一段语音上然后从幅频响应、相频响应、计算量、滤波后的波形细节四个维度做对比。这样报告里的分析有依据答辩时也有讲的底气。2. 语音信号的预处理链路拿到WAV之后先做什么2.1 读文件与基本信息检查很多人拿到一个wav文件就急着设计滤波器这是最容易翻车的地方。我在项目中养成的习惯是先读文件把采样率、声道数、位深这三个基础信息打印出来确认。[y, fs] audioread(speech_noisy.wav); info audioinfo(speech_noisy.wav); disp(fs); disp(info.NumChannels); disp(info.BitsPerSample);采样率是最关键的参数后面所有频率设计都要除以奈奎斯特频率做归一化。常见语音信号的采样率有几种电话语音8kHz宽带语音16kHzCD音质44.1kHz。如果你拿到的文件采样率是44.1kHz却把滤波器参数按8kHz场景设计效果直接错得离谱。如果文件是多声道一般要转成单声道再处理if size(y, 2) 1 y mean(y, 2); end这里有个细节不要直接取某一声道因为左右声道如果存在相位差单声道信息会有损失取平均是最保守的做法。当然如果项目本身就是做立体声处理则另当别论。2.2 去直流、预加重、幅值归一化预处理环节我固定做三步每一步都有实际意义。第一步是去直流y y - mean(y);语音采集过程中麦克风偏置或模数转换器可能引入直流分量。直流分量不去掉的话FFT频谱的0Hz附近会出现一条很高的谱线掩盖低频细节滤波后也可能出现一段缓慢的基线漂移。第二步是预加重y filter([1, -0.95], 1, y);语音信号的高频部分在声带激励和口唇辐射过程中衰减很快预加重通过一阶高通提升高频段能量让频谱在高频区域的细节更清晰。系数取0.95到0.97比较常见数值越大高频提升越明显。这一步不是必须的如果你的任务只是验证滤波器效果可以不做但做了之后频谱图会更好看分析也更方便。第三步是幅值归一化y y / max(abs(y));这样做有两个好处一是避免信号幅值接近浮点上限导致削波二是后面计算信噪比时不同幅值尺度下的结果更具可比性。2.3 分帧加窗和滤波的边界这里有一个很多人混淆的问题语音识别里的分帧加窗和滤波处理到底什么关系实际上滤波器处理是对整段信号流做的不需要先把信号分成一帧一帧。分帧加窗是为了做短时傅里叶变换STFT观察信号频谱在时间轴上的变化。我见过不少同学把这两件事混在一起先分帧再逐帧滤波结果帧边缘出现大量拼接噪声还以为是滤波器设计错了。正确做法是整段信号直接过滤波器需要画时频分布图时再用帧长20到30ms、帧移10到15ms的参数对滤波前后的信号分别做短时谱分析。比如16kHz采样率下帧长取512点帧移取256点加汉明窗。这样得到的对比图既能体现滤波前后的频带变化又能在时间轴上看到噪声压制的动态效果。3. FIR滤波器设计窗函数、阶数估算与Matlab实现3.1 设计指标怎么定FIR低通滤波器的设计要先有明确的频带指标。以我参与的这个项目为例场景是某段采样率为16kHz的带噪语音噪声主要分布在4kHz以上的高频区我们希望保留语音主要频段同时压低高频噪声。指标定为通带截止频率 fp 3400 Hz阻带截止频率 fst 4000 Hz通带最大衰减 Rp 1 dB阻带最小衰减 As 40 dB过渡带宽度是600Hz归一化到采样率就是Δf_norm (fst - fp) / fs 600 / 16000 0.0375注意这里归一化除以的是采样率本身不是奈奎斯特频率。用Matlab的fir1时截止频率则是按奈奎斯特频率归一化的Wn (fpfst)/2 除以 fs/2这一点后面会用到。指标定得好不好直接影响后续阶数。过渡带越窄滤波器阶数越大实时成本越高。如果项目没有特殊要求过渡带不要拍脑袋定参考实际噪声分布和应用需求来定。3.2 窗函数的参数计算FIR设计最直观的方法是窗函数法。思路很清楚理想低通滤波器的冲激响应是无限长的sinc序列直接截断必然引起吉布斯现象所以用一个有限长窗函数去平滑截断控制通带和阻带的纹波。常用的汉明窗阻带衰减大约53dB主瓣宽度约8π/N过渡带宽度和阶数N的关系可以估算为N ≈ 3.3 / Δf_norm代入数据N ≈ 3.3 / 0.0375 88所以用88阶汉明窗理论上可以达到约53dB的阻带衰减超过40dB的指标要求。如果想让阻带衰减更可控可以用凯塞窗。凯塞窗的阶数估算公式是N ≈ (As - 7.95) / (2.285 × Δf_norm)代入As 40dBN ≈ (40 - 7.95) / (2.285 × 0.0375) ≈ 374这个结果比汉明窗大很多原因是凯塞窗可以独立控制主瓣宽度和旁瓣高度公式给出的阶数是保守估计。实际项目中如果追求灵活性更推荐直接用Matlab的kaiserord函数[n, Wn, beta, ftype] kaiserord([3400 4000], [1 0], [0.891 0.01], 16000);其中通带波纹1dB换算成幅值容差是10^(-1/20)0.891阻带衰减40dB换算成幅值容差是10^(-40/20)0.01。kaiserord会自动返回满足指标的阶数和窗参数。不过我个人在实际项目中用得最多的是等波纹优化法也就是Parks-McClellan算法它的优势是在同样的阶数下通带和阻带的纹波分布更均匀通常可以用更低的阶数达到同样的最坏情况指标。Matlab里就是firpm函数。3.3 用fir1和firpm实现并验证汉明窗版本fs 16000; Wp 3400 / (fs/2); % 0.425 Ws 4000 / (fs/2); % 0.5 N_hamm 88; if mod(N_hamm, 2) 1 N_hamm N_hamm 1; end b_fir fir1(N_hamm, (Wp Ws) / 2, low, hamming(N_hamm 1));这里阶数N和窗函数长度差1是Matlab对滤波器阶数的定义造成的初学者很容易有误解fir1(N, ...)生成的是N阶滤波器也就是长度为N1的系数向量窗函数长度也必须是N1。如果这里写错了设计出来的滤波器过渡带位置会有偏差。验证频响[H, f_axis] freqz(b_fir, 1, 4096, fs); plot(f_axis, 20*log10(abs(H))); xline(3400, --, fp); xline(4000, --, fst);freqz返回的频率轴单位如果指定了fs单位就是Hz不指定时是归一化rad/sample。画图之前确认一下单位避免把横轴标错。滤波操作y_fir filter(b_fir, 1, y);这里用filter而不是conv因为filter处理的是与输入等长的输出且符合滤波器实时处理的语义。如果你用conv做卷积要自己处理长度关系和延迟补偿容易出错。3.4 高阶数带来的延迟问题FIR线性相位的代价是群延迟恒定为N/2个采样点。汉明窗88阶延迟44个采样点在16kHz下约2.75ms这个数值不算大。但如果你听了kaiserord的建议用374阶延迟就是187点约11.7ms。对于语音通话应用ITU标准里单向端到端延迟建议在150ms以内所以11.7ms完全可接受。但如果项目后续要接到反馈控制系统里延迟就成了必须纳入考虑的因素因为相位延迟等效于把信号整体后移闭环系统对相位裕度很敏感。在报告里讨论这一点非常加分FIR的设计自由度和相位特性是优点但高阶数延迟不是免费午餐。对比IIR时这个数值差异就是最有说服力的论据。4. IIR滤波器设计从模拟原型到数字实现的完整过程4.1 原型滤波器选型IIR滤波器通常从模拟原型变换而来常见原型有巴特沃斯、切比雪夫I型、切比雪夫II型和椭圆滤波器。它们之间的取舍可以用一个表格说清楚原型通带特性阻带特性过渡带宽度相位特性巴特沃斯单调无波纹单调无波纹较宽四种中最平缓切比雪夫I型等波纹单调较窄有波纹切比雪夫II型单调等波纹较窄有波纹椭圆等波纹等波纹最窄非线性最强课程设计我推荐巴特沃斯理由很实际代码简单数值稳定性好通带无波纹意味着语音频谱不会被额外调制相位特性相对平缓主观听感最不容易出怪声。椭圆滤波器虽然阶数最低但对系数量化误差最敏感实验里容易出现设计指标达标但实际滤波效果发飘的情况。4.2 双线性变换与频率预畸变IIR设计里最难理解的环节是双线性变换。简单说双线性变换把模拟s域映射到数字z域映射关系为s (2/Ts) × (z - 1) / (z 1)这个变换有一个重要性质模拟频率轴会被压缩到数字Nyquist频率范围内频率映射关系是ω 2 × arctan(Ω × Ts / 2)这意味着模拟频率和数字频率不是线性对应的。如果你直接把模拟滤波器的截止频率拿来做数字归一化设计出来的数字滤波器截止频率会偏。解决办法是先对模拟截止频率做预畸变Ω 2 × fs × tan(ωc / 2)其中ωc是期望的数字截止角频率rad/sample。好在这套流程Matlab已经封装好了butter函数内部自动完成双线性变换不需要手动预畸变。但答辩时老师大概率会问原理所以这个公式要能默写出来而且要能解释清楚为什么要预畸变而不是直接把频率等比缩放。4.3 butterbuttord的代码实现用buttord先计算最小阶数fs 16000; Wp 3400 / (fs/2); % 0.425 Ws 4000 / (fs/2); % 0.5 Rp 1; As 40; [n, Wn] buttord(Wp, Ws, Rp, As); [b, a] butter(n, Wn, low);buttord返回的Wn是实际设计时使用的3dB截止频率归一化值它不一定等于3400对应的Wp而是一个位于通带边缘附近的值。理解这一点对正确分析很重要滤波器的-3dB点并不是通带边缘的3400Hz而是由满足指标的最小阶数反推出来的。验证幅频和相频[H_iir, f_axis] freqz(b, a, 4096, fs); subplot(2, 1, 1); plot(f_axis, 20*log10(abs(H_iir))); subplot(2, 1, 2); plot(f_axis, unwrap(angle(H_iir)) * 180 / pi);相位图最好用unwrap展开不然角度在±π之间跳变看不出真实相位趋势。IIR的相位响应是一条明显弯曲的曲线不同频率的延迟不一致这就是非线性相位。4.4 线性相位补偿filtfilt的使用与注意针对IIR非线性相位带来的波形失真Matlab里常用filtfilt做零相位滤波。它的做法是把信号正向过一遍滤波器再把结果反转过一遍正向滤波的相位延迟被反向滤波抵消掉最终输出的相位延迟近似为零但幅度响应会变成原滤波器幅频响应的平方等效阻带衰减翻倍。y_iir_linear_phase filtfilt(b, a, y);这个函数好用但不能无脑用。几个实际注意点第一filtfilt输出结果的幅频特性不是butter设计出的那个响应而是它的平方。所以报告里的滤波器幅频响应图和滤波后的信号频谱图不是直接对应的。第二filtfilt存在起始瞬态问题。信号开头的一小段在没有足够数据时会被滤波器强行处理导致边界处出现较大误差。信号很短时这个误差会扩散到整个信号建议至少保证信号长度是滤波器阶数的十倍以上。第三filtfilt是非实时的因为它需要整段信号都已知才能做双向滤波。如果项目讨论的是实时通话场景filtfilt只能作为离线分析工具不能用于实时处理流程。我的习惯是实时方案用filter离线对比分析用filtfilt两者结果都放在报告里分别解释适用场景。这样既展示了技术覆盖面又体现对工程边界的理解。5. 滤波效果的三维验证频谱、波形与听感5.1 频域验证信号频谱对比滤波器设计出来到底有没有用不能只看幅频响应图。截取滤波前后的语音信号分别做FFT观察频谱变化NFFT 2^nextpow2(length(y)); Y_before fft(y, NFFT); Y_after fft(y_fir, NFFT); f_axis (0:NFFT/2-1) * fs / NFFT; P_before abs(Y_before(1:NFFT/2)).^2 / NFFT; P_after abs(Y_after(1:NFFT/2)).^2 / NFFT; plot(f_axis, 10*log10(P_before eps)); hold on; plot(f_axis, 10*log10(P_after eps)); legend(滤波前, 滤波后);在4kHz以上频率区滤波后的频谱应该有明显下降。这里一个实用技巧是把滤波器幅频响应曲线缩放后叠加在信号频谱图上。因为滤波器响应和信号频谱量纲不同所以要对滤波器响应做适当增益调整但这样可以直观看出信号频谱变化是否与滤波器频率选择一致。5.2 时域波形与信噪比计算频域图够了还要有时域证据。对比原始带噪语音和滤波后语音的波形最明显的变化是高频毛刺减少波形变得更平滑。如果要量化效果可以计算滤波前后的信噪比SNR。由于这类课程设计通常是人工叠加噪声有干净的原始语音参照可以这样算snr_before 10*log10(sum(speech_clean.^2) / sum((speech_noisy - speech_clean).^2)); snr_after 10*log10(sum(speech_clean.^2) / sum((speech_filtered - speech_clean).^2));注意这里分母中的speech_filtered要先与speech_clean对齐。FIR滤波和filtfilt都会引入延迟直接相减会得到一个包含时间错位的大误差SNR计算结果会非常难看。必须先做对齐对FIR把speech_clean也平移N/2个采样点后再比较或者用xcorr做互相关对齐。这是一个非常常见的错误我见过太多同学在报告里给出离谱的负SNR原因就是对齐没做。5.3 实际听感测试的主观评价客观指标再漂亮最后也要过耳朵这一关。我的项目里做了三个维度的主观评价噪声是否明显下降特别是高频噪声带语音清晰度是否保持有没有出现闷的感觉是否有不自然的金属声或空洞感。实测中FIR滤波后的语音听起来更自然因为线性相位保证了波形结构不被破坏IIR滤波后降噪更彻底但某些辅音开头会有轻微的钝化感这与相位失真和通带边缘的具体频率响应有关。主观听感部分不要只写声音变干净了一句话最好给出评级标准比如1到5分让几位同学盲听打分然后列表汇总。这样报告的分析部分就更有说服力。6. 程序与报告的规范化让答辩老师挑不出毛病6.1 程序模块划分与命名规范Matlab课程设计最忌脚本里几百行堆在一起。规范的做法是把功能拆成独立文件main.m主流程按步骤调用各函数preprocess.m读取、去直流、预加重、归一化design_fir.mFIR设计返回滤波器系数design_iir.mIIR设计返回滤波系数evaluate_fir_iir.m滤波、SNR计算、频谱分析plot_results.m汇总绘图命名全部用英文小写不用拼音缩写。每个函数顶部加注释块说明输入输出、参数含义、设计指标。关键代码行加简短中文注释但注释不要写满每一行那样反而影响阅读。6.2 报告的关键章节与图表引用报告结构我推荐按工程项目的逻辑组织而不是按教科书目录摘要一段话说清楚做了什么、得到什么结论。引言项目背景语音降噪的需求为什么需要滤波器。FIR和IIR数字滤波器原理各自数学模型、设计原理、Matlab实现方法。设计方案包括设计指标、整体流程图用文字描述、预处理方案、滤波器选型。实验结果与分析频谱对比图、波形对比图、SNR对比表、主观听感评分。结论与心得两种滤波器的适用场景项目经验教训。参考文献。图表编号要和正文引用对应比如图3-2展示了88阶汉明窗FIR滤波器的幅频响应。正文中每一张图都必须有引用否则这张图就是多余的。表格要有表头对比类表格建议放在结果分析章节不要散乱插入。6.3 结果分析部分的写作技巧结果分析是报告中拉开差距的核心部分。建议用表格对比FIR和IIR的实际测试结果对比项FIR汉明窗/等波纹IIR巴特沃斯阶数约88到100约12相位特性线性非线性实测阻带衰减约53dB满足40dB群延迟44个采样点频率相关不恒定计算量较高较低主观听感自然略闷分析时不要只说FIR好或IIR好而是说明取舍。比如IIR用较低的阶数实现同样的阻带衰减计算效率高但它的非线性相位导致波形失真在基音检测等对波形敏感的场景下可能影响后续算法。FIR虽然阶数高但相位线性波形保真好实现结构始终稳定工程上更可控。如果能够把理论差异和实测现象对应起来答辩时即使问得再细你也有话可讲。7. 实操中真正折磨人的几个坑7.1 边界效应滤波起点的瞬态用filter做IIR滤波时信号开头会有一段明显的瞬态响应样本数大约与滤波器阶数相当。这是因为滤波器起始状态全零信号突然输入时滤波器需要几个周期的响应才能进入稳定状态。FIR也有这个问题只是表现不那么明显。常规解决方法是前补零再裁剪pad_len 1000; y_padded [zeros(pad_len, 1); y]; y_filtered_padded filter(b, a, y_padded); y_filtered y_filtered_padded(pad_len 1 : end);或者直接用filtfilt消除相位效应但前面提过它自身也有边界瞬态问题。对于课程设计我更推荐前补零再裁剪的方法逻辑简单便于讲解。7.2 相位失真对语音波形的影响这个坑经常被忽视。语音波形对时序要求严格的地方不止是人耳还有后续算法。我在项目里对比过同段语音经过IIR和FIR滤波后再提取基音周期结果IIR滤波后的基音检测在某些音节上出现了跳变原因是相位失真造成了波形形状的改变。原因是IIR滤波器对高频分量和低频分量的延迟不同导致语音中的辅音起始位置发生偏移。对于只做人耳主观评价的任务这种偏移不易察觉但对于自动语音识别、语速估计这类任务相位失真会直接影响特征提取精度。报告里谈到IIR的局限性时这一点比相位非线性四个字更有说服力。7.3 参数整定的调试思路如果设计出来的滤波器不满足指标不要盲目增大阶数。先明确是哪一项指标不达标通带纹波太大说明通带边缘控制不够可以加大阶数或改用等波纹设计过渡带太宽说明主瓣宽度不够需要加大阶数阻带衰减不足说明旁瓣还没压下去需要换更好的窗函数或提高阶数。调试时用freqz同时观察幅频和相频用零极点图检查稳定性。IIR滤波器阶数过高时极点会非常靠近单位圆系数量化后可能越过单位圆导致不稳定所以IIR的参数并不是越高越好。7.4 数据类型与数组维度问题最后说几个程序层面很低级但实际经常遇到的错误audioread读入的数字默认是double但如果原文件是16位PCM实际幅值范围是[-1,1]做归一化或SNR计算前要确认量纲双声道信号和滤波器系数做filter时filter会按列独立处理如果你的系数向量是行向量而信号是列向量结果可能完全不对freqz不指定fs时频率轴是rad/sample指定fs后是Hz代码里一旦混用画出来的图会让人完全看不懂如果用filter处理一个列向量和一个行向量Matlab有时会报维度错误最好统一用列向量。我做这个项目时还遇到过一个经典问题某段信号是44.1kHz采样率我却按16kHz设计了滤波器通带截止频率完全错位滤波后语音变得非常闷。检查半天才发现是采样率变量没改只是复制代码时忘了更新。这件事提醒我在项目初期就把采样率作为全局参数写死并注明来源所有后续计算都从同一个变量引用可以省掉很多低级错误。实际上做完这个项目我最大的体会是数字滤波器设计的答案永远不止一个同一个指标FIR和IIR给的是两种风格完全不同的解法。把两种解法的结果量化出来理解差异背后的原因并且能在报告里用数据和图表支撑你的结论这项能力比背住任何一个设计公式都更有价值。如果你也在做类似的项目建议不要只交一个能跑的程序把调试过程中那些为什么记录不对、为什么效果发闷、为什么频谱对不上的问题也写进报告里这些真实的排查过程反而是最能体现工程素养的地方。