公司动态
嵌入式IIR滤波器从设计到STM32落地全指南
做嵌入式音频和传感器信号处理的朋友一定绕不开IIR滤波器。它能在单片机这种资源有限的平台上只用几个乘法器就实现非常陡峭的频率响应音频均衡、心率信号去噪、电源纹波滤除、ADC采样预处理全都有它的身影。我最早接触IIR是在一个低成本的心率血氧模块上当时需要在不换MCU的前提下把工频干扰和基线漂移压下去FIR阶数太高跑不动最终就是靠级联双二阶IIR解决了问题。也是从那时起我意识到IIR滤波器不只是一个“公式套用”系数怎么设计、结构怎么选择、定点怎么处理、为什么一个符号看错就会自激每一项都是坑。这篇文章我会完整梳理IIR滤波器的核心概念、工作原理、设计工具和嵌入式落地全流程重点讲清楚FIR和IIR怎么选、直接I型/SOS矩阵是怎么回事以及STM32上从系数生成到实际跑通的步骤。既适合刚接触数字滤波的开发者建立整体认知也适合正在做工程应用的同学直接抄作业。文章里所有代码和步骤都是我在实际项目中验证过的思路你可以照着跑一遍再根据自己信号的频段去调整参数。1. 为什么我选择IIR滤波器1.1 IIR到底解决了什么问题IIR的全称是Infinite Impulse Response无限冲激响应滤波器。这个“无限”并不是说它永远不停地在工作而是说它在时域上受到一个单位冲激激励之后输出会无限地延续下去不会真正归零。之所以会这样是因为IIR的输出不仅依赖当前和过去的输入还依赖过去已经计算出来的输出也就是存在反馈回路。反馈是这个结构最核心的差异它让IIR可以用很低的阶数实现很高的频率选择性。举个直观的例子你要在2kHz采样率下去掉50Hz工频要求过渡带很窄。用FIR做可能需要128阶甚至256阶每来一个样本就要做几百次乘加运算而IIR用4阶、也就是两个双二阶节就能达到接近的滚降特性。对于MCU来说这个运算量的差距可能就是“跑得动”和“跑不动”的区别。我经常用一个生活化的类比来理解IIRFIR像是一堆人排队传话每个人只根据队伍前面几个人说的话来复述传完一遍就结束IIR则像一个装了弹簧的机械系统你推它一下它会按自身振动特性来回晃很久这些晃动的频率就是滤波器的极点位置决定的。反馈让历史输出“残留”下来所以我们才能用更少的系数去构造更复杂的频率响应。1.2 FIR和IIR滤波器到底有什么区别很多新手在选择滤波器类型时最纠结的就是FIR还是IIR这个选择没有绝对的对错完全看应用场景。我整理了一个对比表格方便你直接对照对比维度FIR滤波器IIR滤波器冲激响应有限长度总会归零无限长度理论不归零输出计算只用输入的历史值同时用输入和输出的历史值传递函数H(z)B(z)只有零点H(z)B(z)/A(z)有零点和极点实现同滚降的阶数高通常数十到数百阶低通常2到8阶相位特性可严格线性相位非线性相位群延迟不均稳定性天生稳定需要约束极点位置计算开销高低数值敏感性较低对系数量化敏感典型应用通信、音视频要求线性相位工业控制、传感器滤波、音频EQ看这张表就能明白IIR的核心优势是“性价比”——用更少的计算量换取更陡峭的滚降。代价则是相位非线性、稳定性设计和数值实现都要更小心。我的实际选择经验是如果信号是传感器采样、电机控制反馈、音频EQ这类对实时性和MCU资源敏感、又不关心波形绝对相位的场景优先用IIR。而如果你的系统需要多个滤波器级联后波形不能有任何相位失真比如某些高精度测量、心电波形分析、通信基带处理那FIR的线性相位特性就是不可替代的优势。换句话说在“资源受限”和“相位失真可接受”这两个前提下IIR是更聪明的方案。1.3 相位非线性的影响有哪些IIR最隐蔽的一个问题就是相位非线性和群延迟波动。信号经过IIR之后不同频率成分的延迟时间不一样反映在时域波形上就是“波形的形状变了”。对单一频点的正弦波你只会看到幅值变化和整体时移对包含多个频率成分的复杂波形比如心电信号里的QRS波群、音频里的瞬态鼓点IIR会让波形在时间轴上“拉偏”极端情况下甚至出现明显的前后振铃。我在做音频EQ的时候就遇到过这个问题。一个4段参数均衡器全部用IIR调出来频响曲线很漂亮但人声一经过就感觉“糊”了一层原因是不同频段的群延迟差异把瞬态轮廓抹掉了。后来我采用了折中方案低频段用IIR中高频段改用FIR或者把Q值限制住在音色和资源之间找平衡。所以不要只看幅频响应图。如果你的应用对波形保真度有高要求在设计完成之后一定要输入一个方波或脉冲信号来观察输出有没有明显的前后振铃、有过冲。这类问题在幅频曲线上完全看不出来。2. 从差分方程到系统函数弄清IIR的工作原理2.1 从一条最简单的公式看IIR的差分方程在做什么IIR滤波器在数字域的数学描述是差分方程一阶IIR低通滤波器的公式很简单y[n] b0 * x[n] b1 * x[n-1] - a1 * y[n-1]一条公式两个输入项、一个输出反馈项这就是最原始的IIR。实际工程中一般会用到二阶甚至更高阶但本质结构不变。通用的N阶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]注意反馈项前面的符号在绝大多数教科书和DSP库中差分方程里反馈是减号但实现时系数存的往往是正的a1、a2程序里做的是减法。很多人移植代码时符号搞反滤波器直接变成震荡器这个坑后面排查章节我会再讲。前向路径上的b系数决定“零点”影响哪些频率被压掉反馈路径上的a系数决定“极点”影响哪些频率被放大或谐振。两者共同作用构成完整的频率响应。对一阶IIR来说b0b10.5、a1-0.5时就是一个最简单的移动平均低通反馈越接近1滤波越平滑但响应也越慢。2.2 Z变换和零极点为什么反馈会发散理解IIR稳定性的关键在Z变换。把差分方程两边做Z变换得到系统函数H(z) (b0 b1*z^-1 ... bM*z^-M) / (1 a1*z^-1 ... aN*z^-N)令分子为0得到零点令分母为0得到极点。零点让某些频率的响应跌到谷底极点让某些频率的响应抬起来形成山峰。极点位置决定了系统是否稳定如果所有极点都落在Z平面的单位圆内部也就是极点模值小于1系统就是稳定的一旦任何一个极点的模值大于等于1输出就会发散或等幅振荡。这里可以继续用弹簧类比极点相当于系统本身的固有振动输入信号只是不断在“喂”这个系统。固有振动衰减极点在圆内那输出最终稳定固有振动发散极点在圆外就算输入停止系统自己也会越振越大这就是所谓的“滤波器自激”。在理想数学环境下设计工具给出来的系数必然稳定但当你把系数转成有限位数的定点数、或者用浮点但中间误差累积时极点可能被硬生生推到单位圆外。我在STM32上做定点滤波时就亲身踩过一次在MATLAB里算好的4阶巴特沃斯系数一切正常用Q15格式量化后直接跑输出发散了。逐项比对量化前后的极点才发现一个原本距离单位圆只有0.002的极点量化后漂移到了1.0001。这个案例让我明白高阶IIR在嵌入式里不能“直接用bh”一定要转成二阶节的级联结构并且对系数做合理的增益分配。2.3 直接I型、直接II型和SOS矩阵说到IIR的实现结构排在第一位的概念就是直接型结构。直接I型结构就是严格按照差分方程来排运算先算前向路径零点部分再算反馈路径极点部分中间用两个延迟链分别存输入和输出历史值。这种结构好理解但缺点是动态范围不好控制——如果滤波器的增益峰值很高中间节点数值会很大容易在定点实现中溢出。直接II型结构通过合并前向和反馈延迟链来减少内存占用理论上最少只需要N个状态变量就能实现N阶IIR。但直接II型中间节点是零极点合并后的总响应动态范围问题比直接I型更严重。在浮点平台上这些问题不明显但在定点处理器上一旦中间节点饱和输出就会产生严重的非线性失真。高阶IIR滤波器最稳妥的实现方式是SOS矩阵也就是Second Order Sections二阶节级联。思路很简单把一个4阶、6阶甚至8阶的系统拆成多个二阶环节联起来每一节只有一对共轭极点和一对共轭零点滤波器系数范围更可控极点的数值敏感性也大幅降低。SOS矩阵通常是一个二维数组每一行代表一个二阶节[b0, b1, b2, a0, a1, a2]其中a0通常是1实际运算时所有系数除以a0归一化。更规范的设计工具会额外输出一个增益参数g可以把这个增益乘到第一节的b0上也可以乘到DC增益上目的是防止中间级输出过大。CMSIS-DSP库里的Biquad结构就是典型的SOS级联实现后面的实操章节我会专门讲。3. 设计工具与系数获取从Matlab/Python到STM323.1 用MATLAB快速生成IIR系数的完整流程如果你手上有MATLABFilter Designer是最快的系数设计工具。在命令行输入fdatool老版本或filterDesigner新版本打开图形界面选择“Bandpass”或“Lowpass”指定采样率Fs、通带边界频率Fpass、阻带边界频率Fstop、通带纹波和阻带衰减。设计算法一般选Butterworth巴特沃斯或者Elliptic椭圆前者过渡带平坦后者滚降陡但对相位影响更大。设计完成后在菜单里选择“Convert to Second-Order Sections”软件会把高阶传递函数拆成SOS矩阵。导出时选“Export to C Header File”可以生成一个包含SOS系数和增益的数组。这里有一个很容易忽略的细节MATLAB导出的SOS矩阵里a0不一定是1而CMSIS-DSP等库默认a01所以移植时要把每一行都除以a0。我自己后来的习惯是不再用GUI而是写脚本批量化处理因为设计参数经常要反复微调GUI点来点去效率太低。把设计流程用脚本固定下来新项目改几个参数就能出系数还能自动画频响曲线建议你也这么干。3.2 实战演算用Python设计一个4阶巴特沃斯低通并导出系数没有MATLAB licence时Python的scipy.signal库就是最顺手的设计工具。下面我用一个具体案例来演示采样率选2kHz想滤除200Hz以上的高频噪音设计一个4阶巴特沃斯低通滤波器。import numpy as np from scipy.signal import butter, sosfreqz, sosfilt fs 2000.0 fc 200.0 # 截止频率 order 4 # 归一化频率截止频率 / (采样率/2) Wn fc / (fs / 2) sos butter(order, Wn, btypelow, outputsos) print(sos)这里最关键的坑就是Wn的计算公式。很多人直接传200/2000得到的结果完全不对因为数字滤波器设计用的归一化频率是相对于奈奎斯特频率的也就是采样率的一半。200Hz截止、2kHz采样率Wn应该是0.2而不是0.1。设计结果是一个3行6列的SOS矩阵每条记录对应一个二阶节。第一行大概是这样的数据[[ 0.02008335 0.0401667 0.02008335 1. -1.56101808 0.64135178] [ 1. 2. 1. 1. -1.06091543 0.41095086]]看到这个结果不要慌它和我们直接用butter(order, Wn)返回的b/a形式不一样但滤波效果是一样的。SOS格式更适合数值稳定的级联实现。设计完成后可以用sosfreqz检查频率响应确认-3dB点在200Hz附近import matplotlib.pyplot as plt w, h sosfreqz(sos, worN8000, fsfs) plt.semilogx(w, 20*np.log10(np.maximum(abs(h), 1e-5))) plt.axvline(fc, colorred, linestyle--) plt.grid(True) plt.show()再用一段仿真信号验证滤波效果t np.arange(0, 1, 1/fs) x np.sin(2*np.pi*50*t) 0.5*np.sin(2*np.pi*500*t) y sosfilt(sos, x)如果频谱上50Hz成分保留、500Hz成分被明显衰减滤波器设计就通过了。3.3 定点量化系数转成Q格式之前必须做的事STM32F4以上或Cortex-M7带FPU的MCU可以直接用float跑IIR但如果你用的是F103这种不带FPU的M3核浮点运算会让CPU负载爆表必须走定点路线。定点化系数不是简单地乘一个缩放因子那么简单有几个大坑必须提前规避。第一先做SOS归一化确保每一节的a01。第二分配增益。把设计工具输出的总增益g合理地分配到各个二阶节上通常做法是让每一节在通带内的增益不超过1避免中间级输出过大。第三把浮点系数乘以2^(N-1)四舍五入取整存成Q15或Q31格式。以Q15为例系数的取值范围是[-1, 1)但实际设计出的a1、a2可能接近甚至超过1。对于超过1的系数你需要先把这一节的整体增益缩小算完再补回来或者用Q31扩大动态范围。我早期做定点时就因为a21.32量化失败输出直接飞了。后来学乖了设计完系数后统一检查模值超过0.99的系数一律做增益拆解。另外定点实现时状态变量也要对齐。输入ADC是12位整型的话先左移到Q15格式再做滤波输出再右移回12位。这个过程看似简单但位宽处理不对信号会变得面目全非。关于定点实现的完整代码我放在下一章结合STM32一起说。4. 在STM32上实现IIR滤波器的完整过程4.1 硬件环境和工程准备我用的开发环境是STM32CubeIDE芯片选的是STM32F407主频168MHz带FPU所以本章示例用浮点实现。如果你的芯片没有FPU可以借助CMSIS-DSP的Q15接口改造思路一样。具体外设配置如下ADC用定时器触发采样率2kHzDMA把采样结果搬运到内存数组滤波处理在主循环里批量执行处理完的数据通过DAC输出或串口发送到上位机绘图。定时器触发采样加DMA搬运的好处是采样间隔稳定不占用CPU中断时间这是做实时滤波的基础。工程里需要添加CMSIS-DSP库。在STM32CubeIDE里可以通过“Software Packs”安装或者直接把CMSIS/DSP目录加入源码路径。如果你用的是旧版本库还要注意在C/C设置里加上ARM_MATH_CM4宏定义和浮点库-mfpufpv4-sp-d16 -mfloat-abihard否则编译器不启用FPU指令。4.2 直接I型C语言实现与逐行讲解在不依赖库的情况下一个4阶IIR直接I型实现的核心就是维护8个状态变量4个输入历史、4个输出历史。下面的代码展示了一个双二阶节的直接I型浮点实现如果要多阶就把多个节串起来。typedef struct { float b0, b1, b2; float a1, a2; float x1, x2; float y1, y2; } biquad_df1_t; float biquad_df1_process(biquad_df1_t *f, float x) { // 直接I型差分方程 // y[n] b0*x[n] b1*x[n-1] b2*x[n-2] - a1*y[n-1] - a2*y[n-2] float y f-b0 * x f-b1 * f-x1 f-b2 * f-x2 - f-a1 * f-y1 - f-a2 * f-y2; // 状态更新注意顺序x2 -- x1 -- xy2 -- y1 -- y f-x2 f-x1; f-x1 x; f-y2 f-y1; f-y1 y; return y; }这里最容易出错的就是状态更新顺序。必须先计算y再用到旧状态绝对不能先刷新状态再算y否则整个滤波器就变成另一个完全不同的系统。我见过很多人把这段拼错出来的滤波效果乱七八糟。直接I型的好处是直观、好调试每节的中间结果都能单独打印。缺点是在定点实现中动态范围不理想所以生产项目我更推荐CMSIS-DSP的级联Biquad结构。4.3 使用CMSIS-DSP库五步配置级联双二阶CMSIS-DSP库把IIR滤波器封装得很完整你不需要自己推导差分方程填对参数就能跑。核心分五步第一定义结构体和缓冲区#include arm_math.h #define NUM_STAGES 2 // 4阶 2个二阶节 arm_biquad_casd_df1_inst_f32 S; float32_t coeffs[5 * NUM_STAGES]; // 每节5个系数b0,b1,b2,a1,a2 float32_t state[4 * NUM_STAGES]; // 每节4个状态x1,x2,y1,y2 float32_t input[BLOCK_SIZE]; float32_t output[BLOCK_SIZE];第二把Python或MATLAB导出的SOS系数填进coeffs数组。假设第一节是[b0, b1, b2, a1, a2]第二节同样。这里注意CMSIS-DSP不存a0默认a01。如果设计工具导出的a0不是1先归一化。第三初始化结构体arm_biquad_cascade_df1_init_f32(S, NUM_STAGES, coeffs, state);第四处理一个数据块arm_biquad_cascade_df1_f32(S, input, output, BLOCK_SIZE);第五重置功能切换参数或重新初始化时很有用arm_biquad_cascade_df1_reset_f32(S);这个接口是按块处理的BLOCK_SIZE建议取32或64。块处理的好处是充分利用CPU缓存、减少函数调用开销实测比逐样本调用快很多。4.4 实时数据流从定时采样到滤波再到输出实际工程里滤波不是孤立跑一个函数而是嵌在数据采集和输出链路里。我常用的框架是这样的volatile bool dma_done false; float adc_buf[BLOCK_SIZE]; float dac_buf[BLOCK_SIZE]; void HAL_ADC_ConvCpltCallback(ADC_HandleTypeDef *hadc) { dma_done true; } int main(void) { // ...初始化ADC、DAC、定时器、DMA... HAL_TIM_Base_Start_IT(htim6); HAL_ADC_Start_DMA(hadc1, (uint32_t*)adc_buf, BLOCK_SIZE); while (1) { if (dma_done) { dma_done false; // ADC是12位无符号先转成有符号浮点 for (int i 0; i BLOCK_SIZE; i) { adc_buf[i] adc_buf[i] - 2048.0f; } // 滤波 arm_biquad_cascade_df1_f32(S, adc_buf, dac_buf, BLOCK_SIZE); // 输出到DAC或串口传送数据 HAL_DAC_Start_DMA(hdac, DAC_CHANNEL_1, (uint32_t*)dac_buf, BLOCK_SIZE); } } }这里有几个实践细节ADC原始数据是0~4095直接滤波会有直流偏移所以先减2048。滤波器输出如果继续走DAC需要判断DAC输出格式必要时还要叠加回直流电平。如果是串口回传改成打包float数据发送即可。性能方面用STM32F407跑4阶IIR滤波每个样本的耗时大约只有几个微秒2kHz采样率下CPU占用几乎可以忽略不计。如果跑8阶IIR同时做三路也完全在预算内。这就是IIR在嵌入式中最大的价值。4.5 调参与验证波形对不上设计值怎么办滤波跑通之后第一时间要做的是验证频率响应是否和设计值一致。我通常用信号发生器输出扫频信号给ADC输入然后在某个频点用示波器看输出幅值或者用串口把滤波前后的数据发到PC端用Python画频谱图对比。别偷懒这一步能避免很多“差不多但不对”的隐患。如果实测的滚降点明显偏离设计值最常见的原因是采样率不对定时器配置的采样率和你设计时用的Fs不一致或者ADC采样触发频率设错了。还有一种情况是你设计时用了单精度浮点但CMSIS-DSP库默认也是单精度两者应该一致不需要担心。如果用的是Q15定点那要重点检查系数量化误差导致的极点漂移。5. 常见问题与排查技巧实录5.1 滤波后信号发散了先从这三个地方排查这是IIR项目里最让人头疼的现象输入一加上输出就直接拉满、饱和或爆音撤销输入输出还在持续增长典型的不稳定表现。必须按顺序排查以下三点。第一查系数符号。差分方程里反馈项是减号但设计工具返回的a1、a2本身就是带符号的。如果你手写实现时把y b0*x - a1*y1写成了y b0*x a1*y1极点符号反了滤波器立刻变成振荡器。我之前就干过一次花了两小时才意识到是手滑多加了个负号。第二查a0归一化。MATLAB导出的SOS矩阵a0不一定是1如果直接把b、a填进CMSIS-DSP的coeffs数组而忘了除以a0相当于每个节都多乘了一个未知增益输出大概率发散。把脚本里加上coeffs[i] / a0就没这个问题。第三查定点量化。把浮点系数直接转成Q15如果某些极点离单位圆很近量化误差极有可能把极点推出去。解决方法是改用SOS级联结构或者用Q31提高精度再不行就重新设计滤波器给极点留出更多余量。5.2 输出有直流偏移或突然跳变状态数组未初始化是这类现象的第一嫌疑。CMSIS-DSP库的init函数会把state数组清零但在部分版本中如果你只定义数组不调初始化就处理数据状态里是随机值输出前几百个样本会有一段奇怪的“过渡”或跳变。解决方式很简单上电后先调用一次reset接口。第二个常见原因是有无符号数混用。ADC输出是12位无符号0~4095你把它直接当作有符号float送进滤波器等于给滤波器叠加了约2048的直流偏置输出自然有直流分量。前面代码里先减2048就是这个道理。第三如果你在系统运行中切换了系数新的状态变量是从旧值继续还是清零会产生不同的瞬态响应。CMSIS-DSP提供的reset接口会在切换参数时把状态清零但如果你不清零输出会有短暂冲击。设计允许的话在参数切换时加一点淡入淡出体验会好很多。5.3 定点实现出现严重噪声或饱和定点IIR的噪声主要来自两部分系数量化误差和乘法累加舍入误差。前者在设计阶段就能预判后者则要靠动态范围控制。Q15格式下每个乘法结果都是两个17位有符号数相乘需要32位保存累加时如果直接用32位很容易溢出。标准做法是用Q15乘法配合Q31累加器每算完一节后缩放到Q15再进入下一节。CMSIS-DSP的Q15接口内部就是这么实现的所以能用库就别自己写。我早期在STM32F103上自己写Q15版IIR没有及时做中间缩放结果滤波信号在幅值接近最大时出现明显的“削顶”失真。后来改成每节输出后右移15位中间节点用32位累加问题就好了。噪声和饱和是所有定点DSP应用的永恒主题没有捷径只能逐级检查中间变量的范围。5.4 输出频响和设计图对不上最后一种典型问题滤波器确实在“工作”但实测的截止频率、增益和设计图对不上。如果不是采样率和硬件问题大部分原因出在归一化频率的计算上。设计工具里的Wn必须是0到1之间的归一化频率算法是截止频率除以Nyquist频率。很多脚本书写正确但一旦采样率或截止频率其中一个用错了单位比如把kHz当Hz整个滤波器的设计目标就全偏了。排查时把设计脚本里所有参数打印出来和实际板卡的采样时钟对照一遍。还有一种容易被忽略的情况CMSIS-DSP处理的是分块数据如果你的BLOCK_SIZE数和实际ADC DMA配置不一致部分数据块会包含旧采样值频谱上会出现时域伪影。检查DMA回调里处理的数组长度是否和传给滤波函数的blockSize一致这也是我踩过的坑。5.5 快速排查速查表为了方便日常调试我把上面所有问题整理成一个速查表现象可能原因排查方法输出发散/自激反馈符号写反逐节打印系数对比设计工具输出发散/自激a0未归一化检查SOS每行a0是否为1输出发散/自激极点量化移出单位圆用SOS级联Q31格式直流偏置大ADC数据未去掉偏移输入先减2048开始时有跳变state数组未复位调用reset接口频响偏移归一化频率算错检查Wn fc/(fs/2)定点时噪声大累加溢出用Q15乘法Q31累加器输出有失真啸叫中间节点饱和逐节分配增益限制级间峰值最后再分享一个我个人的工作习惯任何工程里我接手IIR滤波相关改动一定先把各节系数和状态变量打印出来和设计工具里的数值逐项比对然后分步验证——先用纯理想信号灌进滤波函数看输出是否和MATLAB一致再接入真实ADC数据。这个过程虽然多花十分钟但能省掉后面好几个小时的排查时间。IIR滤波器入门不难真正难的是在不稳定的数值环境和紧张的算力预算里做出一个稳定可靠的结果希望这篇文章能帮你少走一些我走过的弯路。