开场:别被"老技术"这三个字骗了
IIR滤波器,无限脉冲响应滤波器,大学数字信号处理课上最劝退的那一章,也是我工作这些年里用得最多的滤波器。你可能觉得FIR才是万金油,什么场合都能套一个窗函数上去,但我告诉你,在很多资源受限的嵌入式场景里,IIR才是真正的保命方案。
为什么?一句话:同样的滤波效果,IIR用的阶数更低,算得快,省内存。一个3阶的巴特沃斯低通,性能大致能抵得上十几阶甚至几十阶的FIR。对跑在STM32这种Cortex-M内核上的程序来说,这差别不是一点点,是实实在在的算力开销和RAM占用。
这篇文章我不会从Z变换的严格定义开始讲,那套东西教材里多得是。我按自己实际用下来的思路来:先搞明白IIR到底是个什么东西,它和FIR的核心差异在哪,然后讲清楚SOS矩阵和直接I型这两种最常用的实现方式,再落到STM32上系数怎么算、代码怎么写、怎么避开那些害死人的坑。读完你能直接上手,至少不会再用错结构、算错系数。
1. IIR滤波器的本质:输出不仅取决于输入,还取决于过去的输出
1.1 从差分方程看IIR的"反馈"本质
FIR滤波器的差分方程长这样:
y[n] = b0·x[n] + b1·x[n-1] + ... + bN·x[n-N]
注意,输出只和当前及过去的输入有关,没有输出反馈。所以一个单位脉冲进去,经过N个采样点之后输出就归零了,脉冲响应是有限长的,这就是"Finite Impulse Response"名字的由来。
IIR滤波器的差分方程多了一项:
y[n] = b0·x[n] + b1·x[n-1] + ... + bM·x[n-M] - a1·y[n-1] - a2·y[n-2] - ... - aN·y[n-N]
后面这一串带a系数的项,是把过去的输出又加权加回来,这就是反馈。因为输出会不断反馈到输入端,所以一个单位脉冲进去,理论上输出永远不会完全归零(虽然实际上因为有限精度会逐渐衰减到零),脉冲响应是无限长的,这就是"Infinite Impulse Response"。
这个反馈是IIR的灵魂,也是它所有优缺点的根源。反馈让同样的滤波效果只需要很少的系数——阶数低、运算量小、内存占用小。但反馈也带来了稳定性问题:如果反馈系数不合适,输出可能发散,直接飘到天上去。FIR是绝对稳定的(只要系数是有限值),IIR则必须检查极点位置。
1.2 频域视角:IIR的优势从哪来
IIR滤波器在频域上可以做到极陡的过渡带。比如你有一个50Hz的工频干扰,旁边就是你要的100Hz信号,你用FIR想把这个干扰压下去,可能需要60阶甚至80阶,但用IIR,2阶到4阶就差不多能做到。
原因是IIR的系统函数可以写成:
H(z) = B(z) / A(z)
分子多项式B(z)决定零点,分母多项式A(z)决定极点。FIR只有分子,相当于只有零点;IIR既有多项式又有分母,可以实现更复杂的频率形状。极点可以看成是让滤波器在某些频率上"谐振",从而实现高增益的窄带特性,或者急剧变化的相位特征。
用个生活化的类比:FIR是纯粹靠"记住更多历史数据"来平滑结果的算法,像一个记性好但反应慢的人;IIR是"记得住过去的结果并据此调整当前判断"的算法,像一个经验丰富但偶尔会"一根筋"(不稳定)的老手。
1.3 一个看得见的例子:同参数下直接对比
我用一个具体的例子来说明。假设采样率Fs = 1000Hz,截止频率Fc = 50Hz的低通滤波器,MATLAB/Octave里用同样的设计规格来对比。
如果用切比雪夫I型IIR,3阶就能做到通带波纹0.5dB、阻带衰减40dB。但如果用FIR,要达到差不多的过渡带宽度和阻带衰减,用窗函数法设计,估算需要的阶数大约是:
N ≈ (Attenuation_dB - 8) / (2.285 × Δω)
其中Δω是归一化过渡带宽度。算下来至少需要20多阶。在STM32F103这种72MHz的MCU上,每秒钟采样1000次,每个采样点要算20次乘加运算,IIR只需要算6次乘加。高下立判。
2. 从传递函数到实际代码:IIR的三种常见实现结构
这一节是实操的基石。很多新手直接拿高阶IIR系数往代码里一塞,发现输出全是NaN,根本不知道自己错在哪。问题往往出在实现结构上。
2.1 直接I型(Direct Form I):最直观但也最浪费
直接I型是最容易理解的实现方式,就是把差分方程照抄成代码:
y[n] = b0·x[n] + b1·x[n-1] + ... + bM·x[n-M] - a1·y[n-1] - ... - aN·y[n-N]
在代码里你需要维护两个数组:一个是输入的历史缓冲x_buffer,一个是输出的历史缓冲y_buffer。处理每个新样本的伪代码如下:
// Direct Form I 伪代码 y = b0*x + b1*x_buf[0] + b2*x_buf[1] - a1*y_buf[0] - a2*y_buf[1]; // 更新缓冲 x_buf[1] = x_buf[0]; x_buf[0] = x; y_buf[1] = y_buf[0]; y_buf[0] = y;这样有两个历史缓冲,各占M和N个位置。看起来逻辑简单,但要注意两个问题:
第一,存储浪费。虽然IIR阶数低,但每个滤波器都要维护两个缓冲。如果是多通道的音频处理或者多路传感器信号,内存消耗会成倍增加。
第二,数值灵敏度。直接I型在系数取值范围很大时,中间结果可能非常大,在定点DSP或单片机上很容易溢出。这是它最大的问题。
但直接I型也有好处:它对系数量化误差的敏感度相对较低(相比直接II型来说),而且移植简单、容易查错。STM32这类带FPU的MCU上,用浮点运算实现直接I型,性能基本不是问题,出不了大乱子。
2.2 直接II型(Direct Form II):省内存但更矫情
直接II型(也叫Canonical Form)是一种更节省内存的结构。它首先计算一个中间变量w[n]:
w[n] = x[n] - a1·w[n-1] - ... - aN·w[n-N]
然后再计算输出:
y[n] = b0·w[n] + b1·w[n-1] + ... + bM·w[n-M]
这样做的好处是你只需要维护一组状态变量w_buffer,而不是两组,内存开销减少了一半。在很多老的DSP课程里它被大大推崇,但实际工程中用得反而少,原因在于它对系数误差更敏感。特别是在把高阶滤波器系数直接量化成16位定点数时,直接II型的极点位置偏移可能比直接I型更大,稍不留神滤波器就不符合设计规格了。
2.3 级联型(SOS):高阶滤波器的正确打开方式
如果你搜索过"iir滤波器 sos 矩阵",你看到的SOS(Second-Order Sections)就是级联型的标准格式。它的核心思想是:不把高阶IIR滤波器作为一个整体去实现,而是把它拆成多个二阶滤波器,串行级联起来。
为什么这么做?因为高阶多项式在数值计算上非常脆弱。假设你设计了一个8阶的巴特沃斯滤波器,分母多项式A(z)有8个系数,这些系数的动态范围可能非常大,小到10的负几次方,大到几百。在浮点运算中可能还好,但在定点处理器上,量化误差会迅速放大,极点的实际位置和理论位置偏差很大,滤波器可能变成"振荡器"——输出持续抖动甚至发散。
如果把8阶滤波器拆成4个二阶节,每个二阶节的系数范围小得多、数值稳定性好得多,依次计算,每一级的输出是下一级的输入,整体效果不变但数值表现优秀得多。
SOS矩阵的格式通常是这样的,每一行表示一个二阶节:
[b0, b1, b2, 1, a1, a2]举个例子,一个4阶巴特沃斯低通滤波器(Fs=1000Hz,Fc=50Hz)在MATLAB/Octave中设计后,用sos函数输出的就是一个Nx6的矩阵:
sos = 0.0201, 0.0402, 0.0201, 1.0000, -1.5606, 0.6414 1.0000, 2.0000, 1.0000, 1.0000, -1.3643, 0.5098每一行代表一个二阶节。第一行增益较低(b系数小),第二行增益较高(b系数接近1)。两级串联在一起,整体的传递函数就是这两个二阶节的乘积。
级联SOS是实际工程中最推荐的选择。TI的DSP库、CMSIS-DSP里的arm_biquad_cascade_df1_f32函数,以及各种CMSIS滤波器库,都是基于SOS级联结构实现的。
3. 用SOS矩阵串联还是直接算:代码实现的关键区别
现在问题来了:已知SOS矩阵,怎么在代码里实现?
3.1 逐级处理的实现方法
最直接的方法是按顺序处理每一级。注意每一级的输出就是下一级的输入,处理完第一级得到中间信号,再传给第二级。伪代码如下:
// 假设有numSections个二阶节,sos是numSections x 6的系数矩阵 // x_in是当前采样值,state是numSections x 4的状态缓冲 float filter_process(float x_in) { float y = x_in; for (int i = 0; i < numSections; i++) { y = biquad_process(y, &sos[i], &state[i]); } return y; }其中biquad_process输入一个值,输出这个二阶节的输出。它的内部实现(直接I型二阶节):
float biquad_process(float x) { float y = b0*x + b1*state[0] + b2*state[1] - a1*state[2] - a2*state[3]; state[1] = state[0]; state[0] = x; state[3] = state[2]; state[2] = y; return y; }这就是CMSIS-DSP中arm_biquad_cascade_df1_f32函数的思路。这种逐级处理的方法也有额外的调试好处:你可以在任意一级输出处加打印或者断点,确认是哪一级出了问题,排错方便很多。
3.2 直接实现一个"塞满系数"的高阶滤波器为什么危险
有人会想:既然我有8阶滤波器的全部系数,为什么不直接把差分方程里所有a和b系数写进代码?
理论上可以,但实际工程中我强烈不建议。原因还是之前提到的数值稳定性问题。用MATLAB/Octave设计一个8阶切比雪夫II型滤波器,它的极点和零点可能非常接近单位圆,任何微小的量化误差都可能把极点推出单位圆外,导致滤波器不稳定。就算不推到单位圆外,极点位置偏移也会让实际频响曲线和设计规格差异很大,可能原来要求-40dB处实际上只有-32dB。
还有一个隐患:如果输入的x信号太大,中间各级的输出可能很大,在级联结构里你可以在每一级之间重新归一化(很多库会自动做),但直接实现的高阶结构没有这个灵活度。如果中间变量超出浮点数范围,直接就NaN了。
3.3 从MATLAB/Octave导出SOS矩阵的操作流程
这段操作值得仔细看,因为它是搜索"iir滤波器 sos 矩阵"的人最想找到的答案。
在MATLAB/Octave中的标准流程是:
- 设计原型滤波器,比如一个5阶的巴特沃斯低通,采样率Fs=1000Hz,截止频率50Hz:
[z, p, k] = butter(5, 50/(1000/2), 'low');注意50/(1000/2)是归一化截止频率,数字信号处理中频率必须用奈奎斯特频率(Fs/2)归一化,这是最容易搞错的地方。
- 把零极点增益模型转换成SOS矩阵:
[sos, g] = zp2sos(z, p, k);sos就是二级节矩阵,g是全局增益。这里有一个重要技巧:zp2sos可以带参数'up'或'down'指定级联顺序,通常用'up'把最靠近单位圆的极点放在最后,数值稳定性更好:
[sos, g] = zp2sos(z, p, k, 'up');- 把全局增益分配进第一级和第二级。许多库的biquad结构自带增益系数(每个节都有b0/b1/b2),所以你可以把增益
g乘到第一级的b系数上,也可以平均分配到所有级。为了缩小每级的动态范围,通常建议g乘到第一级:
sos(1, 1:3) = sos(1, 1:3) * g;这样后面每一节都不需要额外乘增益。但如果g特别大或特别小,建议在级间观察信号幅度,必要时手动调整增益分配。
- 如果要在STM32上用定点实现,还需要把系数转换成Q格式。比如用16位Q15表示系数(取值范围-1到1之间),但要注意IIR系数的负数范围可能超过-1,这需要额外处理。通常建议直接用浮点(STM32F4以上带FPU)或者用32位定点,而不是用16位。后面我会详细讲这个问题。
4. FIR和IIR到底怎么选:从几个实际工程场景看
4.1 一张表站在需求角度做对比
先给一个实用对比表,这是我做选型时必看的:
| 对比维度 | IIR | FIR |
|---|---|---|
| 相同滤波效果的阶数 | 低(通常2~8阶) | 高(通常几十阶) |
| 计算量 | 小 | 大 |
| 内存占用 | 小 | 大 |
| 线性相位 | 不支持(相位非线性) | 天然支持(关于中心对称) |
| 稳定性 | 需要检查极点 | 始终稳定 |
| 数值敏感性 | 高,需小心实现 | 低 |
| 适合场景 | 资源受限、实时性要求高的嵌入式 | 音频处理、需要无相位失真的数据采集 |
表格里最关键的指标是线性相位。
4.2 相位敏感场景:FIR胜出
如果你的应用是对信号做高精度测量,比如振动分析、心电信号处理、音频滤波,波形的时域形状很重要,那么IIR的非线性相位可能是个问题。IIR滤波器在某些频率点会有较大的群延迟波动,信号经过滤波器后不同频率成分的延迟不同,时域波形会发生畸变。
举个例子,你用IIR滤波器处理心电信号,它的P波、QRS波群、T波含有不同的频率成分,IIR滤波器会导致这些波的相对时间关系发生偏移,医生一看波形就觉得不对。你以为滤波后信号变干净了,实际上它已经被"扭曲"了。这种情况下,FIR的线性相位特性就是刚需。
但如果你只是从传感器数据里滤掉一些噪声,不关心波形的具体相位,比如做温控的PID控制里滤掉高频抖动、电池管理系统里读取电流电压的平均值,那IIR完全没有问题。
4.3 资源有限场景:IIR胜出
这又要说回STM32了。假设你用STM32F103(主频72MHz,没有FPU,只有单精度浮点仿真库)做一个电机电流环的采样滤波。电流环的采样频率可能要到10kHz甚至20kHz。如果你用FIR滤波器,阶数50,每次采样要做50次乘加运算,在无FPU的芯片上这已经是很重的开销了。但如果用IIR 2阶滤波器,每次采样只需要大约6次乘加运算,少了一个数量级。更关键的是状态缓冲只需要4个浮点数,而FIR需要50个。
所以在硬实时系统中,IIR往往是唯一的合理选择。你可能觉得"反正计算量也不大",但你要考虑整个系统的预算——中断里除了滤波还有PID计算、通信协议处理、状态机判断。每多1微秒的开销,在高速控制回路里都是要命的。
5. STM32上实现IIR:从系数获取到实测
这是搜索"stm32 iir滤波器直接i型 系数"的人最关心的部分,我直接按完整流程来。
5.1 在STM32上用浮点直接I型实现2阶IIR
先提供一个可以直接用的函数,这是最经典的2阶直接I型实现,系数从MATLAB/Octave导出后手动填入:
typedef struct { float b0, b1, b2; float a1, a2; float x1, x2; // 输入历史 float y1, y2; // 输出历史 } iir_biquad_t; float iir_biquad_process(iir_biquad_t *f, float x) { float y = f->b0 * x + f->b1 * f->x1 + f->b2 * f->x2 - f->a1 * f->y1 - f->a2 * f->y2; // 更新历史 f->x2 = f->x1; f->x1 = x; f->y2 = f->y1; f->y1 = y; return y; }注意这里的符号约定:在MATLAB中,滤波器系数形式一般是y[n] = b0*x[n] + ... - a1*y[n-1] - ...,所以代码里面减法要对应好。很多人就是在这里搞错了符号导致滤波器完全不对。
5.2 从MATLAB/Octave获取系数并检查
以STC/STM32项目中常用到的一个50Hz工频陷波器为例,采样率Fs=500Hz,希望滤除50Hz干扰。可以用iirnotch函数:
Fs = 500; fo = 50; bw = 5; % 3dB带宽 [b, a] = iirnotch(fo/(Fs/2), bw/(Fs/2));得到类似这样的系数:
b = [0.97551, -1.10060, 0.97551] a = [1.00000, -1.10060, 0.95102]注意这里a的第一个元素是1,代码里不需要用它,但需要确认。把b0、b1、b2、a1、a2分别填到结构体的对应字段里就行。
在把系数烧到单片机之前,强烈建议先在电脑上做一次快速仿真验证。可以在MATLAB/Octave里用freqz(b, a, 1024, Fs)看频响曲线,确认陷波点位置和带宽正确;再用stepz或impz看时域响应是否收敛。这一步花5分钟,能省去后面在板子上调试几小时。
5.3 定点还是浮点:STM32上怎么选
STM32F0、F1这类不带FPU的芯片,用浮点数做实时的IIR是有代价的。虽然编译器有软件浮点库,但乘法运算会扩展到几十条汇编指令,速度慢。这种情况下有两种方案:
方案一:用Q15或Q31格式实现在STM32上不太推荐自己搞。STM32官方库里的arm_biquad_cascade_df1_f32是单精度浮点版本,需要用带FPU的Cortex-M4/M7/M33。CMSIS-DSP里也有arm_biquad_cascade_df1_q15和arm_biquad_cascade_df1_q31定点版本,可以直接用,省心很多。
方案二:如果你必须用STM32F1系列,说实话我建议先算算实际负载再决定。一个2阶IIR每秒钟跑1000次,软件浮点也就需要几百微秒,对于慢速采样场景其实无所谓。只有当采样率很高或多路滤波时才需要考虑软件浮点性能瓶颈。
方案三:用定点手动实现系数,非常不建议初学者做,因为会遇到饱和、舍入误差、极限环振荡等问题。
我个人的偏好是:如果成本允许,直接上STM32G4或STM32F4,带FPU,用单精度浮点实现IIR,开发速度和调试体验远超定点方案,价格差距也就几块钱。
5.4 一个简易测试方法:板子上怎么确认滤波器没写错
代码写完了,板子跑起来了,怎么确认滤波器正常工作?
我的做法是:在滤波函数入口处注入一个已知信号,观察输出是否和MATLAB对同一信号的仿真结果一致。
具体操作:
- 在单片机里生成一个固定频率的正弦波数组(比如放在Flash里的const数组,或者直接在代码里用简单的sin函数生成)。
- 给定一个包含50Hz低频和300Hz高频混合的信号,经过滤波器后观察输出波形是否只留下低频部分。
- 把UF问题:在输出端加一个调试变量,通过串口或者JTAG调试器把数据点抓出来画波形。
如果输出波形形状和MATLAB仿真基本一致(幅值允许有微小误差),滤波器就通了。这一步别跳,曾经我身边不止一个人跳过这一步直接接到实际传感器上,结果滤波器系数里符号弄错,折腾了一天最后发现低频全被滤掉了。
6. 实际踩坑记录:IIR滤波器工程化的五个教训
最后这部分是我自己踩坑踩出来的,每一件都付出过时间成本。
6.1 初始瞬态问题:滤波器的"启动冲击"
IIR滤波器有反馈,所以它有一个"建立时间"。假设初始化时状态变量全部清零,突然输入一个大的阶跃信号,输出会有一个很大的过冲,可能持续几毫秒甚至更长。
在很多控制系统中,这个过冲会导致执行机构猛地动一下,这是不可接受的。
我的处理办法是:系统启动时先不要马上让滤波器输出控制量,而是先让传感器信号稳定一段时间,或者在初始化时把滤波器的输入历史状态赋值为当前采样值的合理估计,让滤波器"从工作点开始"。比如传感器刚上电时输出接近0,就把x1、x2、y1、y2都初始化成0,如果预期传感器稳定在某个偏置电压,你可以先采N个样本求平均然后设置初始状态。
6.2 系数量化误差:16位定点的灾难
我第一次在STM32F103上做IIR时,天真地把浮点系数直接转成Q15格式,结果滤波器特性完全乱套。原因很简单:IIR系数的动态范围经常超过Q15能表示的范围,比如系数可能在0.998到1.002之间,Q15量化后直接变成一个小整数,极点的精度损失太严重。
后来我改用浮点加上CMSIS-DSP的库,问题就消失了。如果确实只能用定点,至少要选Q31,并且考虑用级联SOS结构。16位定点做IIR,尤其是高阶IIR,基本等于自找麻烦。
6.3 采样率不稳定的影响
IIR系数是基于固定采样率设计的。如果一个系统用定时器中断触发采样,但是定时器优先级设置不当,导致采样周期抖动,滤波器实际表现会和设计差很多。
最典型的例子是:你用HAL_Delay或者简单的while循环来做ADC采样触发,时钟不准确,采样率漂移,你以为是5kHz采样并计算截止频率,实际可能是4.7kHz,滤波器特性偏移还不算大,但如果触发被其他中断打断,采样间隔忽长忽短,IIR的反馈会累积误差,输出出现抖动。
解决办法是用硬件定时器触发ADC采样,或者用DMA,保证采样间隔精确。这是嵌入式工程师做任何数字滤波都必须刻进DNA的规范。
6.4 状态变量精度:float还是double
Cortex-M4的FPU有单精度浮点,速度极快但精度只有大约7位有效数字。对于2阶IIR来说,这精度足够。但对于8阶以上的高阶IIR,单浮点可能不够用,状态变量累积误差会导致输出噪声增加。
我的建议:8阶以下放心用float,更高阶建议用double或者保证SOS级联顺序合理。还有,如果你需要极窄带滤波器(比如Q值很高的带通或陷波器),float可能不够,因为极点离单位圆太近,对精度极其敏感。我的经验阈值是:如果滤波器的极点模值超过0.999,用double。
6.5 滤波器稳定性:必须在真实运行范围内验证
很多滤波器的极点位置是在设计阶段算好的,但实际运行过程中,信号频率和幅度的变化可能导致滤波器中间变量过大,出现饱和等非线性现象,使得系统不稳定。特别是定点实现中,饱和溢出后状态变量变成极大值,反馈恢复不过来,滤波器就卡在"抱死"状态。
所以稳定性验证不能只看理论,要加到最恶劣的输入信号进行测试,比如信号满幅值、接近截止频率的阶跃信号。如果出现异常,首先检查各节输出是否触及数据范围上界,必要时加入饱和保护逻辑。
结尾:一些关于IIR的个人心得
IIR滤波器不是银弹,但它绝对是嵌入式工程师和算法工程师工具箱里不可缺少的一把锤子。我用它做过电机电流环的陷波、心率传感器的噪声抑制、电源纹波数据的平滑、音频回声消除的前置滤波,每个场景都踩过不同的坑,但核心方法论是一样的:先算清需求指标,选对阶数和类型,用MATLAB/Octave设计并用SOS级联实现,最后在板子上用已知信号验证。
如果你想从IIR开始入手,我建议先从2阶巴特沃斯低通开始,把差分方程、直接I型实现、频响验证这套链路走通,再尝试高阶的切比雪夫、椭圆滤波器。在这个过程中要特别留意那个符号约定的问题——MATLAB里a系数带着负号出现在差分方程里,很多初学者在这里翻车不是少数。
最后分享一个小技巧:如果你的系统里信号频率范围比较宽,又担心IIR的相位失真影响波形,可以在信号链路里先做一次正向滤波,再做一次反向滤波(所谓零相位滤波),但这只适合离线数据处理,实时系统里就别想了。实时系统的话,要么接受IIR的相位特性,要么踏踏实实上FIR。搞清楚自己的需求边界,比纠结"哪个技术更高级"重要得多。