news 2026/10/15 16:53:46

无FPU的MCU上实现纯整数二阶IIR滤波器:原理、定标与代码

作者头像

张小明

前端开发工程师

1.2k 24
文章封面图
无FPU的MCU上实现纯整数二阶IIR滤波器:原理、定标与代码

上手先讲个真实场景。我前阵子给一块用低成本MCU做的采集板调信号质量,输入是压力传感器的模拟量,经ADC采样后毛刺特别明显,高频抖动能把显示数字最后两位带得来回跳。首选的方案自然是经典二阶IIR低通,但芯片没有硬件浮点单元,软件浮点一次乘加就要好几条指令,传感器采样率一提到10kHz以上,CPU占用率直接飙到四成多,别说再做显示刷新和通信协议,连基本任务调度都开始抖。后来我把滤波器整个改成纯整数二阶滤波器,乘法全部用固定点整数运算,CPU占比一下降到5%以内,滤波效果和浮点版本几乎看不出差别。这篇文章就把我从原理、系数推导、定标到最终嵌入式代码的完整过程摊开讲,希望能给同样在资源受限平台上做信号处理的朋友一个可以直接抄作业的参考。

1. 为什么嵌入式里要死磕纯整数滤波

嵌入式工程师听到“浮点”两个字,第一反应往往是“贵”。这个贵不是指芯片价格,而是指运算代价。带硬件FPU的MCU,一次浮点乘加可能只要一两个周期,但这类芯片往往要多花不少预算。很多低成本场景用的是不带FPU的Cortex-M0或者老式8位机,浮点运算完全靠软件库模拟,一次乘法要几十个周期,跑起来相当心疼。

纯整数二阶级滤波器就能把这个问题绕开。它的所有系数预先算成整数,滤波过程里没有一位小数点,运算只有整数乘加和移位,几个周期就能完成一次输出。尤其适合传感器数据调理、电机电流环反馈、电源环路控制这类实时性要求高、又不方便外挂DSP的场景。

还有个容易忽略的优点:确定性。浮点运算在不同编译器、不同优化等级下,舍入行为可能有细微差异。整数运算的中间结果在所有平台上一模一样,这在多机联调、离线回放复现数据时特别有用。我调试某工业项目时,就是靠这份确定性在PC上复现了MCU里的每一个滤波输出点,排查问题省了不少时间。

说白了,纯整数滤波器不是“买不起FPU的将就方案”,而是一种稳定性、实时性、可移植性都更可控的工程选择。

2. 二阶滤波器原理与系数设计

2.1 从差分方程说起

二阶无限脉冲响应滤波器的离散域差分方程长这样:

y[n] = b0*x[n] + b1*x[n-1] + b2*x[n-2] - a1*y[n-1] - a2*y[n-2]

式子里x是输入序列,y是输出序列,四个历史量分别代表前两次的输入和输出。b0、b1、b2是前馈系数,a1、a2是反馈系数。滤波器的性格——是低通、高通还是带通,截止频率多少,有没有过冲——全部由这五个系数决定。

这里有一个新手容易绕晕的地方:符号。差分方程的反馈项写的是减号,这跟教科书传递函数H(z)的分母写法是对应的。分母是1 + a1*z^-1 + a2*z^-2,对应差分方程里就是- a1*y[n-1] - a2*y[n-2]。在你把系数存入结构体时,最好直接按差分方程的“减”来理解,否则程序里符号很容易搞混。我自己早年就吃过这个亏,后来统一约定:结构体里存的a1、a2已经带符号,运算时一律写成减。

2.2 双线性变换得到低通系数

设计滤波器系数的方法很多,双线性变换是嵌入式里最常用的一个。它的思路很简单:先写一个模拟低通滤波器的传递函数,再做一次从模拟域到数字域的映射。核心公式是:

s = K * (1 - z^-1) / (1 + z^-1)

其中:

K = tan(π * fc / fs)

fc是希望的截止频率,fs是采样率。注意这里用tan做预畸变。因为双线性变换在数字角频率和模拟角频率之间存在非线性压缩,如果不做预畸变,设计出来的实际截止频率会偏离目标值,高频偏差尤其明显。

以最简单的二阶巴特沃斯低通为例。模拟归一化传递函数的分母是1 + 2ζs + s^2,对二阶巴特沃斯来说阻尼比ζ=0.7071,对应的Q=0.7071。代入双线性变换,整理后可以得到一组比较干净的系数公式:

a0 = 1 + 2ζK + K^2 a1 = 2*(K^2 - 1) a2 = 1 - 2ζK + K^2 b0 = K^2 b1 = 2*K^2 b2 = K^2

再把a0除到分母外面,让分母首项变成1,就得到差分方程里需要的那五个最终系数:

b0' = b0 / a0 b1' = b1 / a0 b2' = b2 / a0 a1' = a1 / a0 a2' = a2 / a0

顺手算一组。设fc=50Hz,fs=1000Hz,ζ=0.7071。K=tan(π*50/1000)=0.1584,K^2=0.0251,2ζK=0.224。于是:

a0=1.2491,a1=-1.9498,a2=0.8011,b0=b2=0.0251,b1=0.0502。

归一化得到:

b0'=0.0201,b1'=0.0402,b2'=0.0201,a1'=-1.5611,a2'=0.6414。

这组数后面做整数化还会用到。

如果你想设计其他响应,比如巴特沃斯之外的切比雪夫、贝塞尔,只需改变模拟原型里的ζ值,或者引入纹波系数,公式结构不变,所以这套计算流程可以覆盖大多数场景。

2.3 为什么二阶就够用

一阶低通的滚降只有20dB/十倍频,衰减能力偏弱。传感器信号一旦混入宽频噪声,一阶滤波后高频部分依然残留不少。二阶巴特沃斯的滚降是40dB/十倍频,比一阶强得多,实现代价却只多了一个状态变量和几次乘加。

三阶、四阶衰减厉害,但在嵌入式定点环境下,成本上升很快。系数敏感度变高,量化误差更容易把极点推到单位圆外。而且很多传感器噪声频谱的主要能量集中在几百赫兹以内,二阶低通配合适中的截止频率,已经能把信噪比改善几个数量级。这就像盖房子,一阶是个棚子,二阶是个正经平房,三阶四阶是楼房——但平房在大多数情况下已经完全够住。

3. 整数化的核心:定标、动态范围与溢出控制

3.1 把小数系数变成整数

上一节算出来的系数全是小数,比如b0'=0.0201。MCU处理这种数,要么用浮点,要么用一个技巧——定标。定标的思路很朴素:既然你小数位数太多,我干脆把整个差分方程乘以2的N次方,把所有系数变成整数,算完再右移N位把结果缩回来。

差分方程变成:

y[n] = (B0*x[n] + B1*x[n-1] + B2*x[n-2] - A1*y[n-1] - A2*y[n-2]) >> N

其中:

B0 = (int)(b0' * 2^N) B1 = (int)(b1' * 2^N) B2 = (int)(b2' * 2^N) A1 = (int)(a1' * 2^N) A2 = (int)(a2' * 2^N)

取整可以选择就近取整或截断,我习惯用(int)floor(x*scale + 0.5)做就近取整,误差更小。

上面那组50Hz系数,选N=12:

B0=82,B1=165,B2=82,A1=-6394,A2=2627。

差分方程就变成:

y[n] = (82*x[n] + 165*x[n-1] + 82*x[n-2] + 6394*y[n-1] - 2627*y[n-2]) >> 12

计算过程里没有任何小数,全是int或int64乘法。

3.2 定标位数怎么选:中间量的动态范围

选N不是越大越好,也不是越小越好,关键看中间累加结果会不会溢出。

还是上面那组。假设输入x是16位ADC的满幅输出,范围±32768,滤波器的输出y也控制在±32768量级。累加器里每一项的绝对值上限是:

系数绝对值乘32768。最大的是A1对应项,约6394乘32768,约2.1e8。五项累加,最坏情况可能到8e8到1e9。这个量级已经逼近int32的极限,2^31约等于2.147e9。

所以用32位int做中间累加不是不行,但余量不宽。只要输入信号动态范围再大一点,或者系数因为截止频率靠近奈奎斯特频率而变大,累加器就可能在瞬时极端输入下翻车。

我的习惯是直接用int64_t做累加器。嵌入式编译器对64位乘法会生成优化代码,虽然比32位慢一点,但在低通滤波这种每样本只做五次乘加的场合,开销完全可以接受。求稳,别省这一个类型。

如果实在想用32位累加,可以把定标位数降到N=10或N=8。但系数量化误差会变大,滤波器的实际幅频响应会和设计值有偏差。取舍关系如表所示:

定标位数N系数量化误差中间累加器安全度适用场景
8较大,截止频率可能偏1%-3%很安全,int32随便用对频率精度不敏感的粗滤波
12较小,适合多数传感器int32需小心,推荐int6416位信号常见选择
15很小,接近浮点精度必须int64高精度信号处理

3.3 一个很实用的小验证:直流增益

定标换算做完之后,怎么快速判断系数没算错?有个非常实用的检查方法。一个稳定的低通滤波器在直流处的增益应当是1,也就是说,输入一个恒定不变的值,输出最终要收敛到输入值。

用整数系数验证起来很简单,把五个整数系数全部代入公式:

dc_gain = (B0 + B1 + B2 - A1 - A2) / 2^N

代回上面的例子,82+165+82+6394-2627=4096,正好等于2^12。这个结果如果等于1(乘以定标后等于2^N),说明系数换算的直流增益正确。如果不是,说明至少某个系数取整错了一位数,或者符号写反了。

这个办法我每次写滤波器都用,两分钟就能排除一大半低级错误。

4. 嵌入式代码实现

4.1 结构体与核心滤波函数

话不多说,直接上代码。这里用纯C实现,任何嵌入式平台都能用。

#include <stdint.h> typedef struct { /* 定标后的整数系数,a1/a2 已带符号 */ int32_t b0; int32_t b1; int32_t b2; int32_t a1; int32_t a2; /* 定标位数,对应差分方程的右移量 */ uint32_t shift; /* 历史状态 */ int32_t x1, x2; int32_t y1, y2; /* 输出限幅用,可选 */ int32_t min_out; int32_t max_out; } biquad_i_t;

核心滤波函数长这样:

int32_t biquad_i_process(biquad_i_t *f, int32_t x) { int64_t acc; acc = (int64_t)f->b0 * x + (int64_t)f->b1 * f->x1 + (int64_t)f->b2 * f->x2 - (int64_t)f->a1 * f->y1 - (int64_t)f->a2 * f->y2; /* 可选:加舍入偏置,再右移 */ #if 1 if (acc >= 0) { acc += (int64_t)1 << (f->shift - 1); } else { acc -= (int64_t)1 << (f->shift - 1); } #endif /* 状态更新 */ f->x2 = f->x1; f->x1 = x; f->y2 = f->y1; f->y1 = (int32_t)(acc >> f->shift); /* 可选:限幅 */ if (f->y1 < f->min_out) { f->y1 = f->min_out; } if (f->y1 > f->max_out) { f->y1 = f->max_out; } return f->y1; }

这里有两个细节值得展开说。

第一,为什么中间累加器用int64_t。看最坏情况,乘法项里A1可能到六千多,乘上y1的32768,单这一项就超过两个亿,再累加五项容易突破int32的2.1亿上限。int64_t范围大得多,这是最简单可靠的办法。

第二,右移前为什么要加舍入偏置。acc >> shift是向下取整,对正数来说相当于直接砍掉小数部分。这会造成一个系统性的负偏差,尤其在信号接近零的位置会产生可感知的波形台阶。加上偏置后变成就近舍入,误差更小。注意我对acc>=0和acc<0分开处理,避免负数的算术右移偏置方向搞反。

4.2 初始化与硬编码系数

初始化有两条路线。

路线A:初始化函数在运行时接收浮点系数,内部完成定标换算。好处是同一套代码可以在不同截止频率间切换,坏处是库函数里还得带浮点运算和数学库,代码体积增大。

路线B:把定标后的整数系数直接硬编码写进结构体。这是我最推荐嵌入式场景使用的。系数在PC或主机端用公式算好,写死到C文件里,目标MCU上完全不需要浮点计算,程序连数学库都不必链接。切换滤波器特性时,改几个常量重编译即可。

以fs=1kHz、fc=50Hz、Q=0.707、定标N=12为例:

const biquad_i_t biquad_lpf_50 = { .b0 = 82, .b1 = 165, .b2 = 82, .a1 = -6394, .a2 = 2627, .shift = 12, .x1 = 0, .x2 = 0, .y1 = 0, .y2 = 0, .min_out = -32767, .max_out = 32767, };

上电后将结构体清零,或者在每次输入值发生变化前手动reset:

void biquad_i_reset(biquad_i_t *f) { f->x1 = 0; f->x2 = 0; f->y1 = 0; f->y2 = 0; }

注意reset后前几个输出存在过渡过程。如果系统对启动瞬态敏感,可以把历史状态初始化成当前输入值,比如f->y1 = f->y2 = f->x1 = f->x2 = 首采样值,这样输出从第一个点就基本贴住输入,没有爬升过程。

4.3 代码再往细里抠一抠

类型转换的坑值得多说一句。(int64_t)f->b0 * x,为什么要显式把b0转成int64_t再乘?因为b0是int32,x也是int32,C语言会把两者都提升为int32做乘法,结果溢出后才赋值给int64,到时候已经晚了。写成(int64_t)f->b0 * x强制让乘法在64位下进行,是必须的。

状态更新顺序也有讲究。必须先更新x2和y2,再更新x1和y1,顺序搞反会串位导致输出变成另一个滤波器。我自己用代码生成器批量产生滤波器时,专门给这个写过静态断言。

5. 实战:给压力传感器信号去噪

5.1 参数选型思路

场景是一块采集板,压力传感器输出0.5V到4.5V,经12位ADC采样,采样率fs=1000Hz,信号里主要噪声来自工频耦合和高频开关毛刺,目标是把噪声压下去,同时保留200Hz以内有用的脉动波形。

截止频率选择要遵循一个原则:截止频率要比有用信号最高频率高,但又要比噪声主频低。这里有用信号脉动大概在几十赫兹,工频干扰是50Hz附近但不算太强,因此我把fc定在80Hz比较合适。当然这也意味着50Hz工频还是会漏过来一些,想要更狠的抑制,需要加陷波器,那是另一个话题。

Q值这里选0.7071,即巴特沃斯响应,通带最平坦,没有明显过冲。如果后面跟的环节怕过冲,可以选Q=0.5,也就是贝塞尔响应,群延迟更平但下降更缓。对压力信号,巴特沃斯是最稳的起点。

5.2 系数换算全过程

还是上面的公式,先把fc=80Hz、fs=1000Hz代入:

K=tan(π*80/1000)=tan(0.2513)=0.2572。

K^2=0.0662,2ζK=1.4142*0.2572=0.3637。

a0=1+0.3637+0.0662=1.4299 a1=2*(0.0662-1)=-1.8676 a2=1-0.3637+0.0662=0.7025 b0=b2=0.0662 b1=2*0.0662=0.1324

归一化:

b0'=0.0463 b1'=0.0926 b2'=0.0463 a1'=-1.3061 a2'=0.4913

定标N=12,取整:

B0=190 B1=379 B2=190 A1=-5350 A2=2012

检查直流增益:190+379+190+5350-2012=4097,相比2^12=4096偏差1。这个偏差来自系数取整的累积误差,完全可以忽略,实际直流增益为0.9998。

写进C代码:

const biquad_i_t biquad_lpf_80 = { .b0 = 190, .b1 = 379, .b2 = 190, .a1 = -5350, .a2 = 2012, .shift = 12, };

这样整个滤波器就完成了。整个过程不需要在MCU上做任何浮点运算,所有设计计算都在PC上完成。

5.3 实测效果与波形表现

实际跑起来,先给一个1kHz采样率的伪方波,叠加约满幅5%的白噪声。滤波后明显看到高频毛刺被压平,上升沿还能保持锐利,没有明显的高频损失。再给一个固定直流输入,输出大约在几十个采样点内稳定下来,稳定后波动在±1个LSB以内(Q12定标下),这就是16位整数量化的理论地板。

对比浮点仿真结果,两者在主波形上几乎不可区分,最大瞬时误差在输出满幅的0.1%以内。对绝大多数传感器应用来说,这个精度足够了。如果你的项目要求更高精度,把定标位数提到N=15,同时所有中间累加保持int64,量化误差还能再降一个数量级。

6. 常见问题与调优经验

6.1 问题速查表

现象原因排查与解决
输出缓慢漂移或不稳定系数符号写反,a1/a2存储与差分方程不一致用直流增益验证公式检查,确认B0+B1+B2-A1-A2近似等于2^N
输出前几十点严重偏小启动时历史状态为零,滤波器需要建立时间用首个采样值初始化全部状态,或丢弃前N个输出样本
突然出现大幅跳变中间累加器溢出,通常是int32累加且输入接近满幅换成int64_t累加器,或降低定标位数N
固定输入时输出周期性小幅振荡定点滤波器的极限环现象提高定标位数,使用舍入而非截断,必要时加1LSB抖动
滤波后仍有很多高频残留截止频率选太高,或阶数不够降低fc,检查是否真的需要更高阶滤波器
输出与浮点仿真差别明显定标位数太低,系数量化误差大提高N到15,状态更新保持int64

6.2 溢出的深水区:什么时候必须上64位

一个容易被忽略的溢出场景是y1、y2不再局限于16位范围。比如你在做闭环控制,输出作为误差量叠加到控制器,y的幅值可能超过ADC量程,这时候A1*y1单一项就能轻松到几十亿。另一种场景是滤波器串联。两个二阶滤波器的中间量如果都放在16位,很可能截断,正确做法是每个滤波器的输出仍然保留int32,状态变量全部int32,累加器一律int64。

更隐蔽的陷阱:输入信号本身如果来自32位累加器如定时器计数,幅值范围可能是±2^31。这时候无论定标N取多大,累加器都必须使用int64,否则输入那一步就先爆了。我一般这样判断——写代码之前先估算:

max_acc = (|B0|+|B1|+|B2|+|A1|+|A2|) * max_amplitude

这个值如果明显小于2^31,int32可行;接近或超过,直接int64。五项的绝对值之和我算过,很多常见配置下大概在2e9附近,刚好逼近int32的边界,所以保守派直接int64。

6.3 量化误差、极限环与稳定性

量化误差的来源有两个:一是系数的取整误差,二是输出每次右移时丢掉的小数部分。系数取整误差可以靠提高N来减小,但会挤占累加器的动态范围。输出量化误差则与舍入策略强相关,我测试过截断与就近舍入,在中低频段差别不大,但在输入信号幅度较小时,就近舍入的波形明显干净得多。

极限环是个很经典的现象:输入固定不变时,理想滤波器输出应该收敛到固定值,但定点实现可能出现几个LSB的周期性振荡。根因是状态变量的小数部分丢失,使得系统在平衡点附近形成了一个不收敛的小环。缓解办法不多,一是提高定标位数,二是在状态更新时用舍入而不是截断,三是在系统层面接受这个误差,因为几个LSB的振荡在大多数控制回路里都在死区范围内。

稳定性的检查我习惯给滤波器喂单位阶跃信号,看输出是否收敛。判断条件很简单:|A2| < 2^N且|A1| < 2^N + A2时要小心,整数环境下靠近边界时极容易越过。如果你的系数经过取整后触发了这个判断,稳妥起见可以把极点往圆心方向挪一点,即把A1、A2的绝对值各减小万分之几,换来的是稳定余量。

6.4 一些个人调试心得

第一,永远先测直流增益。每次改了系数,先用恒定输入跑一遍,输出如果收敛不到输入值,系数换算或者定标一定有问题,别急着接传感器找原因。

第二,做一个扫频测试。嵌入式上不方便用标准库的sin,可以预先算一张正弦波表,频率从1Hz到400Hz步进扫,记录每个频率点输出的幅值,画出来和理想幅频曲线对比。这个测试能同时发现系数错误、极限环、截止频率偏移三个问题,性价比很高。

第三,把滤波器的中间变量用调试器全部拉出来,加入一个只用于测试的峰值记录。跑一轮后检查峰值有没有逼近int64边界,这比事后看输出波形更早发现隐患。

第四,代码里尽可能把shift做成编译期常量,不要留成运行期变量。编译器遇到常量移位会生成更紧凑的指令,我实测大概能省4到6个周期每样本,在高速采样里这个收益不小。

我自己现在做MCU上的信号调理,第一反应就是列一列能不能用纯整数搞定,而不是急着挂浮点库。这个思路从低通滤波器出发,扩展到高通、带通、带阻也只是改系数的问题,整套框架可以一直复用。如果你在项目中试出了更极端的定标方案,或者踩到了我没提到的坑,欢迎交流,这种问题往往聊一聊就能少走一大段弯路。

版权声明: 本文来自互联网用户投稿,该文观点仅代表作者本人,不代表本站立场。本站仅提供信息存储空间服务,不拥有所有权,不承担相关法律责任。如若内容造成侵权/违法违规/事实不符,请联系邮箱:809451989@qq.com进行投诉反馈,一经查实,立即删除!
网站建设 2026/10/15 16:52:34

MATLAB图形模式识别实战:圆、正方形、矩形自动分类

简介&#xff1a;这份资源聚焦图形模式识别中的基础几何形状检测&#xff0c;面向正在学习数字图像处理、计算机视觉或MATLAB编程的本科生、研究生及工程技术人员&#xff0c;帮助解决圆形、正方形与矩形自动识别与分类的算法实现问题。压缩包共218个文件&#xff0c;约1.07MB&…

作者头像 李华
网站建设 2026/10/15 16:52:30

HART温度变送器直连以太网的工业落地实践

1. 这不是普通温度采集——高温管道变送器直连以太网的底层逻辑在某建材集团下属的水泥熟料生产线车间里&#xff0c;我第一次见到那台安装在回转窑尾气管道上的HART温度变送器时&#xff0c;它表面温度计显示外壳已达85℃&#xff0c;而内部测点正实时监控着420℃以上的高温烟…

作者头像 李华
网站建设 2026/10/15 16:52:27

AT24C02读写实战:从I2C时序到蓝桥杯掉电保存全解析

接到一个做蓝桥杯单片机赛项准备的朋友问起AT24C02&#xff0c;他在这一块卡了好几天&#xff0c;读写数据死活不对&#xff0c;有时候读回来全是0xFF&#xff0c;有时候写完掉电再上电又恢复成初始值。这玩意儿确实坑不少&#xff0c;尤其是竞赛里要求在几分钟内调通、在几秒内…

作者头像 李华
网站建设 2026/10/15 16:51:50

RISC-V E300时钟系统详解:PLL配置与低功耗设计

调过的板子里&#xff0c;十个有八个“看起来是软件问题”的怪毛病&#xff0c;最后都能在时钟树上找到根源。UART出来全是乱码&#xff0c;SPI时序对不上&#xff0c;PWM频率差一大截——第一反应怀疑代码&#xff0c;第二反应怀疑信号完整性&#xff0c;折腾一圈回头才发现&a…

作者头像 李华
网站建设 2026/10/15 16:51:47

通过LiveNVR实现摄像头集中服务器上录像存储备份与回放、以及远程回放现场摄像头sd卡录像或录像机现有录像

通过LiveNVR实现摄像头集中服务器上录像存储备份与回放、以及远程回放现场摄像头sd卡录像或录像机现有录像1、视频流媒体2、录像类型3、录像回看4、查看录像4.1、时间轴视图4.2、列表视图5、相关接口5.1、云端录像接口5.2、设备录像接口6、录像计划7、相关问题7.1、录像存储位置…

作者头像 李华
网站建设 2026/10/15 16:51:43

Qt无触摸屏按键交互:从GPIO到QKeyEvent的完整方案

简介&#xff1a;本资源面向嵌入式开发初学者与树莓派爱好者&#xff0c;提供一套基于QT框架的无触摸屏系统按键交互方案&#xff0c;重点解决无触控设备下如何通过物理按键或模拟按键控制硬件的问题。包内共7个文件&#xff0c;包含2个cpp源文件、1个h头文件、1个ui界面文件、…

作者头像 李华