很多做音频处理的人第一次碰均衡器习惯直接拖一个现成的FIR设计工具把频响曲线拉出来一看挺漂亮可一放进实时音频链路就傻眼延迟大、CPU吃紧、调参还不跟手。我早先在做一个模拟简易调音台项目时就把均衡器改成了“串联二阶IIR”结构用MATLAB把原型写出来后后面不管是批量处理音频素材还是移植到实时环境都轻松了很多。这篇文章想把那套串联二阶IIR实现音频均衡器的MATLAB代码完完整整拆开讲为什么选这个结构双二阶系数怎么生成逐样本滤波怎么写频响怎么验证又有哪些容易踩的坑。代码可以直接复制跑通适合想自己写音频算法原型的人也适合准备把MATLAB原型搬到C/C的嵌入式音频开发者。1. 均衡器为什么非要用串联二阶IIR1.1 均衡器本质上只是在“按频段改幅度”先说均衡器到底在干什么。你在面板上看到的每一个旋钮背后都是一个滤波器它决定了“哪个频率附近能改多少dB”以及“这个改动影响多宽”。这三个数就是中心频率fc、增益gain和Q值。Q值越大影响的频带越窄Q值越小曲线越平缓。比如想突出人声就选1k~3kHz附近的一个峰想消掉房间里的嗡嗡声就在50~80Hz做一个窄带陷波。把这些频段的改动叠加起来最终就是一条总的幅度响应曲线。这里的关键词是“叠加”。一种做法是把每个滤波器的输出直接相加另一种做法是把它们一个个串起来前一个的输出进后一个。本文用的是后者串联二阶IIR。总传递函数就是每一段双二阶的乘积频响上表现为每段频响的逐点相乘画到dB坐标上就是逐段相加。这种方式最大的优点是调试直觉很强我改某一段的参数主要影响这一段对应的频点附近虽然级联的幅频严格说是乘积不会完全互不影响但通常比并联相加时的相位叠加要可控得多。可能你会问为什么不用FIRFIR当然能做任意幅频响应而且能做到线性相位但代价是抽头数量大。一个48kHz采样率下、中心频率1kHz、Q值1的峰值滤波器用FIR实现等效精度往往需要几百甚至上千个抽头而一个双二阶IIR只要5个系数每个样本每个频段只需要几次乘加。做实时音频、嵌入式效果器省CPU是硬需求。IIR的相位不是线性的但对均衡器来说人耳对均衡的相位变化敏感度远低于幅度变化所以工程上普遍接受IIR方案。1.2 为什么是“二阶”不是一阶也不是高阶一阶滤波器的频率响应是固定的6dB/oct斜坡做不了真正的峰和陷波要想做出峰或者陷波至少需要一对共轭极点和一对共轭零点也就是二阶。二阶IIR又叫双二阶Biquad传递函数是这样一个标准形式H(z) (b0 b1z^-1 b2z^-2) / (1 a1z^-1 a2z^-2)分子决定零点位置用来形成陷波或补偿分母决定极点位置用来形成谐振峰或共振特性。一组共轭极点加一组共轭零点正好覆盖均衡器最常见的两种需求峰值提升和窄带衰减。那为什么不直接写一个高阶IIR呢高阶直接型滤波器虽然也能拉出复杂的频响但系数量级差异极大量化误差很容易把极点推离单位圆稍微动一下、量化一下就发散。拆成多个二阶IIR再串联每一段只负责一对极点稳定性检查、系数调试都是独立的这也是很多专业音频软件内部的实际做法不追求单段有多陡而是用合理的二阶段组合出整个目标曲线。1.3 串联和并联的区别为什么这里选串联并联结构同样能实现均衡器思路是把每个滤波器的输出直接相加整体是求和。但并联在实现时需要注意各段输出同步而且相邻频段的相位响应会产生复杂干涉调试曲线时经常出现“这里明明设了3dB叠加之后却变成-1dB”的情况。串联则是频响相乘曲线形状更可控也更容易从模拟电路里找到对应——很多经典模拟均衡器本质上就是滤波级的级联结构。所以我这个例子选串联一段双二阶的输出进入下一段双二阶一路串下去。每一段可以单独调节整条链的稳定性只需要逐段检查极点不需要处理一个高阶多项式。后面的代码也全部围绕“按顺序过一遍每一段”来写。2. 可以直接跑的MATLAB代码参数表、系数函数、滤波循环、频响2.1 先把均衡器参数表定义出来写均衡器原型的第一步不是写滤波器而是把你要的均衡曲线定义清楚。我的习惯是用三个长度相同的数组中心频率 fc、Q 值、增益值 gaindB每个位置对应一个频段。比如一个五段均衡器可以写成fs 44100; fc [100, 350, 1000, 3000, 8000]; Q [1.1, 1.1, 1.0, 1.0, 1.0]; gdB [5, -3, 2, 4, -2];这段配置的意思是100Hz附近提升5dB350Hz衰减3dB1kHz提升2dB3kHz提升4dB8kHz衰减2dB。用参数表而不是在代码里写死每段滤波器好处是后期很容易改成配置文件或者从界面上的旋钮直接映射进来。你甚至可以把这个数组导出成CSV由另一个脚本自动生成不同预设的均衡器。2.2 系数生成双二阶均衡器系数公式有了参数表接下来要把每个频段的增益、频率、Q值换算成5个滤波系数。这里我用的是一套音频领域广泛流传的双二阶均衡系数公式网上常叫audio EQ cookbook专门给实时音频用的。下面这个函数生成单段峰值式系数也就是Peaking EQ可以提升也可以衰减function sos peaking_coefficients(fc, Q, gaindB, fs) % 返回单段二阶IIR的传递函数系数按a0归一化 % sos [b0, b1, b2, a1, a2]对应的差分方程 % y[n] b0*x[n] b1*x[n-1] b2*x[n-2] - a1*y[n-1] - a2*y[n-2] A 10^(gaindB/40); w0 2*pi*fc/fs; alpha sin(w0)/(2*Q); b0 1 alpha*A; b1 -2*cos(w0); b2 1 - alpha*A; a0 1 alpha/A; a1 -2*cos(w0); a2 1 - alpha/A; sos [b0/a0, b1/a0, b2/a0, a1/a0, a2/a0]; end里面 A 和 alpha 是公式的中间变量后面我会单独解释。这里你只要知道传入中心频率、Q、增益dB和采样率返回一行系数。有个细节值得留意当 gaindB 0 时A 10^(0/40) 1b0/a0、a1、a2 会和分母对应系数完全相同这段滤波器就变成纯直通不会改变声音。这是一个很有用的特性方便你在调试时把某一段临时旁路。2.3 逐样本滤波转置直接II型结构系数算出来之后把它们串起来逐样本处理。我不直接用MATLAB的 filter 函数而是写一个显式的逐样本循环因为这才是实时处理的真实形态来一个样本过第一段再过第二段一路串到最后。结构用转置直接II型Transposed Direct Form II每个双二阶只需要两个状态变量 z1、z2function y biquad_cascade_process(x, sos) % 串联二阶IIR均衡器逐样本处理整段信号 % sos: 每一行是 [b0,b1,b2,a1,a2]行数等于滤波器段数 N length(x); M size(sos,1); y zeros(size(x)); z zeros(M,2); % 每段两个状态变量 for n 1:N u x(n); for k 1:M b0 sos(k,1); b1 sos(k,2); b2 sos(k,3); a1 sos(k,4); a2 sos(k,5); v b0*u z(k,1); z(k,1) b1*u - a1*v z(k,2); z(k,2) b2*u - a2*v; u v; end y(n) u; end end这段代码把第k段的输出直接当成第k1段的输入从第1段一直走到第M段实现了真正的串联。z(k,1)和z(k,2)分别是第k段在转置直接II型结构里的两个状态。要注意更新顺序先算当前输出v再更新z1最后更新z2。如果把z2先更新了再回头算z1会用到错误的状态值代码里这个顺序是有讲究的后面讲移植时还会再提到。2.4 频响展示与一段整体可运行脚本有了系数和处理函数就可以画频响曲线验证了。整体频响是各段的乘积我习惯在对数频率点上逐点计算。下面是一个完整示例定义五段均衡、画整体频响、生成简单测试音频并处理% eq_demo.m fs 44100; fc [100, 350, 1000, 3000, 8000]; Q [1.1, 1.1, 1.0, 1.0, 1.0]; gdB [5, -3, 2, 4, -2]; sos zeros(length(fc), 5); for k 1:length(fc) sos(k,:) peaking_coefficients(fc(k), Q(k), gdB(k), fs); end % 生成一段2秒测试音440Hz 5kHz混合 t (0:fs*2-1) / fs; x sin(2*pi*440*t) 0.3*sin(2*pi*5000*t); y biquad_cascade_process(x, sos); % 绘图整体频响 f logspace(log10(20), log10(fs/2), 2000).; z exp(-1j*2*pi*f/fs); H ones(size(f)); for k 1:size(sos,1) b [sos(k,1), sos(k,2), sos(k,3)]; a [1, sos(k,4), sos(k,5)]; H H .* polyval(b, z) ./ polyval(a, z); end figure; semilogx(f, 20*log10(abs(H)eps)); grid on; xlabel(Frequency (Hz)); ylabel(Magnitude (dB)); title(Cascaded Biquad EQ Response); xlim([20 fs/2]);绘图代码里用 polyval 对复数 z 做多项式求值MATLAB 会逐点返回结果再用点乘逐段累积画出来的就是整条串联均衡器的总响应。如果想单独看某一段把H初始化换成那一行的b/a再画一次就行。我建议你跑通这一段脚本后先把gdB全部改成0这条曲线应该是一条0dB直线再把其中某一段改成6dB看看峰是不是出现在你要求的中心频率上。这一步能验证你手里的系数函数工作正常。3. 系数公式背后为什么双二阶公式这么写3.1 A和alpha的来路以及那个容易搞错的dB参数那个 A 10^(gaindB/40) 经常让人迷惑增益是dB为什么不直接 10^(gaindB/20)其实这是经典双二阶均衡公式里的一个约定A是公式里的中间增益因子峰值滤波器在中心频率附近的增益和A的平方相关所以公式里用dBgain/40而不是dBgain/20。如果你自己改成 10^(gaindB/20)同样的增益值会得到完全失真的系数曲线幅度直接不对。我第一次自己推导时就栽在这里后来干脆记住这个约定cookbook公式原样写先照抄再理解。alpha sin(w0)/(2Q) 则可以粗略理解为阻尼项经离散化后的产物。整组公式可以看成先在连续时间域设计一个带谐振峰的滤波器再用双线性变换离散化使数字频率 w0 对应模拟原型的目标频率。sin和cos的组合就是离散化后自然出现的量。由于这套公式在设计时已经处理了从模拟到数字的映射你直接用 w0 2pi*fc/fs 是安全的不需要额外再手动预畸变。3.2 Peaking、Low Shelf、High Shelf怎么选专业均衡器的两端通常不用Peaking型。频率太低时峰值曲线形状很难做得自然频率接近Nyquist时Peaking又容易在高频边缘翘起一个不想要的尾巴。更常见的做法是最下面一段用Low Shelf最上面一段用High Shelf中间频段用Peaking。Peaking对中心频率附近进行提升或衰减Shelving则整体抬升或压低某一段频响形似缓坡。上面的示例代码只写了Peaking是为了让代码更短、更容易抓住核心。正式做均衡器时建议在最底和最顶各补一段Shelf。判断标准也很简单当中心频率低于100Hz或高于10kHz时如果用了Peaking出现“过冲”、频段边界噪声变大先怀疑是不是应该换成Shelf型而不是硬调Q值。3.3 采样率、频率范围边界和被忽略的物理带宽自己写系数最容易被忽略的是采样率。同一个fc、同样的Q和增益在44.1kHz和8kHz采样率下对应的模拟频带宽度完全不同听起来也就不一样。公式里 w0 是数字角频率已经和fs绑定所以设定参数时一定要清楚当前工程用的是什么采样率。另外fc必须严格大于0且小于fs/2。我在写参数表时经常设一个接近20kHz的峰值段如果采样率只有44.1kHz曲线会在高频剧烈抖动一旦fc超过Nyquistsin和cos的组合就失去物理意义了。建议在参数加载时统一检查一下fc不要超过0.45*fs给奈奎斯特频率留一点安全余量尤其是后续要插值或做变采样处理时在线末端很容易出现高频杂讯。4. 级联顺序、数值稳定性和状态量初始化那些坑4.1 级联顺序到底要不要管线性时不变系统串联的顺序不影响总频响这是教科书结论。但在浮点实现里顺序会影响中间信号的动态范围和误差积累。我的经验是正增益段尽量往后放负增益段尽量往前放。比如某段要9dB某段要-6dB把-6dB放前面信号先被压下去再经过9dB中间峰值不大不容易削波反过来先9dB再-6dB中间会出现一个比最终输出大得多的信号在定点或低headroom环境里就危险了。如果所有段都是double精度输入信号本身余量也足够顺序的影响可以忽略这时候我更建议按频率从低到高排不是为了更好听而是为了查参数、看图方便。我自己的标准做法是低频、中频、高频升序排列同时把增益大的段尽量往链路靠后的位置放。你可以在代码里做一个简单的重排根据gdB的符号和数值调整sos的行顺序直观感受一下中间信号峰值的变化。4.2 验算极点别让系数悄悄发散串联二阶IIR稳定性的条件是每一段分母的极点都在单位圆内。MATLAB里可以写个循环检查对每段的分母系数求根然后判断模长。for k 1:size(sos, 1) r roots([1, sos(k,4), sos(k,5)]); if any(abs(r) 1) error(第 %d 段有不稳定极点, k); end end实际使用中我遇到发散的情况并不多除非Q设得特别大比如大于50或者采样率特别高。当极点离单位圆非常近时double精度下还能勉强工作换成float精度就可能因为四舍五入把极点推到圆外。所以如果你要面向嵌入式定点或float平台建议在界面层把Q限制在0.3~10这个范围系数计算保留double只有滤波运算转成float。另一个边界情况是fc非常低比如20Hz且Q很小alpha可能接近0这时b0和b2几乎相等差分方程里的高频分量很小状态变量几乎不更新听觉上就是这段完全没反应。这不是代码错了而是Q和fc的组合物理效果太弱可以考虑降低fc或增大增益也可以用Low Shelf来替代低频Peaking。4.3 状态变量清零、爆音和流式处理biquad_cascade_process函数每次调用都会把z初始化为零适合一次性处理完整音频。但要做实时音频插件或连续处理流这就出问题了下一帧开始必须保留上一帧末尾的z状态否则每帧开头都会出现明显的“咔哒”声。正确做法是把z存成类成员或持久变量每次处理只更新它不在函数内部清零。参数在运行中变化也会引入爆音。一种简单有效的方法是把增益、频率这些参数本身做平滑而不是让系数突变。例如用户把100Hz增益从0dB拉到6dB你可以让gaindB在几毫秒内按线性或对数曲线插值不断重新计算系数。由于参数变化是渐进的系数变化也是渐进的爆音会明显减少。想彻底消除爆音还需要做输出crossfade但多数场景下参数平滑已经够用。这个小技巧对后面做自动化控制特别有帮助直接改系数会有很多“啪嗒”声平滑参数基本没有。5. 怎么验证这个均衡器没有白写5.1 理论响应和实测响应交叉验证写完代码第一件事是对比设计目标和实际频响。先看第二节里的绘图代码这叫理论曲线然后还要做一次“实际测量”对滤波器的输入给一个单位冲激把输出做FFT得到的就是真实频响或者直接给一段正弦扫频看输出包络。我在调试某个模拟均衡器原型时就遇到过理论曲线和实测曲线在3kHz附近差1.5dB的情况。查到最后是参数表里有一段Q误写成10另一段中心频率设成了超过Nyquist的35kHz两者叠加出来一个异常峰。频响对比能把这些低级错误一次性暴露出来所以不要只信仿真图一定要有步骤地和实测结果对一遍。5.2 用对数扫频信号做完整的实测验证推荐用对数扫频信号过一遍算法再做短时傅里叶变换看包络。生成对数扫频并处理的MATLAB代码可以这样写T 5; n round(T*fs); tt (0:n-1)/fs; f_start 20; f_stop 20000; inst_f f_start * (f_stop/f_start).^(tt/T); % 瞬时频率对数扫频 phase 2*pi*cumsum(inst_f)/fs; x_chirp sin(phase); y_chirp biquad_cascade_process(x_chirp, sos); % 对y_chirp做短时傅里叶变换观察包络即可 spectrogram(y_chirp, hann(1024), 512, 1024, fs, yaxis);扫频的好处是能直观听出来从低频一路扫到高频哪个频段被抬、哪个频段被压和频响图完全对应。我习惯在每次调完参数后先用耳朵扫一遍再回头看频响图调整细节比纯盯图快很多。如果扫频结果里某个频点出现断续或异常明显的“呜呜”声基本可以判断是极点太靠近单位圆或者参数表里有重复峰叠加。5.3 常见问题排查表我把常见现象、可能原因和解决方向整理成一个小表方便你对着查现象可能原因检查方向某频段完全没效果中心频率超出范围或增益为0检查参数表fc应小于0.45fs频段交界处出现凹陷或凸起相邻峰靠太近、Q值偏大适当减小Q或重新规划频段间隔输出出现NaN或爆音极点发散状态未正确传递运行极点检查检查帧间状态高频发虚、噪声大高频段用了大Q、或Shelf没做好改用High Shelf限制Q上限整体增益比预期低级联顺序导致中间削波正增益段后置、负增益段前置预留headroom遇到曲线异常时第一件事是把某一段单独拿出来画频响比如用sos(2,:)单独算而不是盯着一堆叠加后的曲线猜。单段正常了再检查两段组合、三段组合很快就能定位到出问题的那一节。6. 从MATLAB原型到实时系统的移植思路6.1 MATLAB里提高吞吐量的写法本文的biquad_cascade_process用双循环是为了还原实时逐样本流。如果只是离线批处理几百万个样本这个循环在MATLAB里偏慢可以改用内置filter函数按段串联y x; for k 1:size(sos,1) y filter([sos(k,1), sos(k,2), sos(k,3)], [1, sos(k,4), sos(k,5)], y); end这段代码效果和双循环在数值精度内一致速度会快很多。但真要移植到实时环境反而是前面那种逐样本逐段更新的写法更贴近硬件模型每个采样点依次经过所有双二阶也就是一个样本进、一个样本出。所以你理解算法时用循环批量验证时用filter两者思路完全一样。6.2 C/C移植时的关键点把上述函数搬成C并不难关键是状态管理。建议为每段定义一个结构体typedef struct { double b0, b1, b2, a1, a2; double z1, z2; } Biquad; double process_biquad(double x, Biquad *f) { double y f-b0 * x f-z1; f-z1 f-b1 * x - f-a1 * y f-z2; f-z2 f-b2 * x - f-a2 * y; return y; }处理一帧音频时对每个样本依次调用M个process_biquad。帧与帧之间要保持状态千万别在函数内部把z清零只在初始化时清零。系数计算可以在初始化阶段或参数变化时重新算不要每个样本都算系数那样开销太大。嵌入式平台上我一般用float保存滤波状态用double只做系数计算。编译开最高优化等级后M段级联的开销就是每样本M次结构体更新对主流MCU来说几段、十几段都能实时跑。如果还不够可以把多个Biquad结构体放进连续数组用数组索引代替结构体指针缓存命中率会好一些。6.3 这套代码还能怎么扩展从这个基础结构出发能扩展出不少东西把中心频率按ISO 1/3 octave排布就能做成31段图形均衡器在参数表里加一列滤波器类型就能同时支持Peaking、Low Shelf、High Shelf把每一段的增益参数做几毫秒的线性插值就能支持自动化控制和旋钮渐变。更进阶的可以在每个频段输出端加一个软削波器防止多个正增益段叠加后超出0dBFS。如果你要做真正的参数均衡器另一个重要扩展是“参数联动修改”当用户同时拖动频率和Q值时新系数不能瞬间跳变。我习惯在参数变更时重新计算系数再做glide插值或者先对增益值做低通滤波再更新系数。这个扩展在实时效果器里非常实用可以从最简单的一阶低通平滑开始先让声音不爆音再慢慢追求更干净的过渡。最后说点个人体会。我在实际做这个原型时最常犯的错不是数学公式而是忘了做“前端参数检查”fc设置超出奈奎斯特频率、Q值填了个极端值、级联顺序没规划。所以现在我的每个均衡器代码都会在初始化时做一次参数合法性检查并且把正增益段放在链路后部这几乎成了肌肉记忆。另一个小技巧是调试时先用扫频信号听一遍再瞄频响图耳朵往往比曲线更快暴露问题。如果你也打算做自己的第一版均衡器建议先用这篇文章的代码跑通五段再慢慢加Shelf和参数平滑后面会顺很多。