做嵌入式这些年,滤波器是我写得最多的代码之一,尤其是ADC采样之后的低通滤波。但真正让我下定决心把RC低通滤波器换一种写法,是在一个小项目连续踩了两个坑之后:一是软件浮点运算在低成本MCU上慢得离谱,二是直接用除法实现一阶RC系数,中断里一次滤波就吃掉上百个周期,采样率稍微提上去,CPU就被拖死。
那之后我开始研究用移位运算替代除法来重构RC低通滤波器,把一阶、二阶一直到三阶的级联结构全部定点化。这篇文章就把这套实现思路完整写出来,包括离散化推导、系数选型、移位替代除法的原理和误差、定点状态管理、完整可跑的C代码,以及我在实测中遇到的截断、初值、稳定时间这些坑。适合正在做传感器采集、音频预处理、工频干扰抑制的嵌入式开发者参考。
1. 为什么是“三阶”:一阶和二阶在真实噪声场景下的窘境
1.1 一阶-20dB/dec远远不够用
RC低通滤波器的幅频特性,单节的公式是:
[ |H(f)| = \frac{1}{\sqrt{1 + (f/f_c)^2}} ]
它的问题是衰减斜率太缓,只有-20dB/decade。什么意思呢?如果你的截止频率设定在10Hz,那么到了50Hz,按照一阶特性也只能衰减到差不多0.2倍左右,也就是输出里还残留约20%的干扰。在很多场景下,尤其是需要把工频干扰压到1%以下,或者要滤除高频毛刺的场景,一阶滤波器完全不够用。
你可能想,那一阶不行就二阶呗。二阶确实能做到-40dB/dec,但实现方式很关键。如果用两个一阶节直接级联,截止频率处的滚降特性会叠加,50Hz衰减可以到4%左右,很多时候还是不够狠。而且很多MCU平台上做二阶IIR如果直接套用标准形式的差分方程,系数一多,定点化之后误差管理就变得麻烦。
所以我最后直接上了三阶。三节级联的幅频特性是单节的立方:
[ |H_{3}(f)| = \left( \frac{1}{\sqrt{1 + (f/f_c)^2}} \right)^3 ]
同样是10Hz截止频率,50Hz处衰减可以到0.8%以下,也就是-42dB左右。对于绝大多数的传感器信号预处理、ADC抗混叠之后的数字再滤波,这个衰减深度已经完全够用了。
1.2 三阶级联的-3dB点不是单节的-3dB点
这里有个非常容易踩的坑:三个截止频率相同的单节级联之后,整体-3dB点会左移,并不是原来的单节截止频率。
计算很简单但容易忽略。三节整体幅度平方为:
[ |H_{3}|^2 = \frac{1}{(1 + u^2)^3}, \quad u = f/f_c ]
令整体-3dB,即(|H_{3}|^2 = 0.5),解出来:
[ 1 + u^2 = 2^{1/3} \approx 1.26 ]
所以:
[ f_{-3dB} \approx 0.51 \times f_c ]
换句话说,如果你想要三阶整体-3dB点在10Hz,单节截止频率不能直接设在10Hz,而是要设在19.9Hz左右。这个0.51的修正系数,我第一次设计时就漏了,结果滤波器比预期偏“闷”了差不多一倍。
2. 从模拟RC到数字差分:离散化其实只有一行
2.1 一阶节的差分方程是怎么来的
我们常见的模拟一阶RC低通电路,微分方程是:
[ RC \cdot \frac{dy(t)}{dt} + y(t) = x(t) ]
把它离散化,用最简单的差分近似,可以得到数字一阶RC节的经典递推形式:
[ y[n] = y[n-1] + \alpha \cdot (x[n] - y[n-1]) ]
其中系数(\alpha)满足:
[ \alpha = \frac{T}{RC + T} \approx \frac{T}{RC} = \frac{1}{2\pi f_c T_s} \cdot \frac{1}{?} ]
这里要写清楚。设采样周期为(T_s),采样率为(F_s = 1/T_s),那么:
[ \alpha \approx 2\pi f_c \cdot T_s = \frac{2\pi f_c}{F_s} ]
这个近似在(\alpha)远小于1时非常准确,实际工程里RC滤波器截止频率远低于采样率,(\alpha)通常都在0.01到0.25这个区间,所以这个近似完全可以放心使用。
于是我们可以反推:
[ f_c \approx \frac{\alpha \cdot F_s}{2\pi} ]
这个式子直接决定了,只要我要的截止频率确定了,在某个采样率下,(\alpha)就定了。而(\alpha)的取值,恰恰就是“移位替代除法”的核心切入点。
2.2 级联三阶与直接高阶IIR的取舍
三阶RC低通,在数字域实现有两种路径。
第一种是直接写出三阶差分方程,形如:
[ y[n] = b_0 x[n] + b_1 x[n-1] + b_2 x[n-2] + b_3 x[n-3] - a_1 y[n-1] - a_2 y[n-2] - a_3 y[n-3] ]
这种标准直接型结构看起来一步到位,但坑很多:需要维护3个输入历史值和3个输出历史值,一共6个状态变量;系数经过z变换之后不再是直观的(\alpha),而是组合出来的多项式系数,定点量化时很容易让极点偏移,甚至导致滤波器不稳定。我在一个项目里试过直接用标准型三阶IIR,系数稍微量化狠一点,输出就开始震荡,排查起来相当痛苦。
第二种就是三个一阶节级联。每个节只有一行代码:
[ y_k[n] = y_k[n-1] + \alpha_k (y_{k-1}[n] - y_k[n-1]) ]
三个节加在一起才3个状态变量。而且级联结构有一个非常好的性质:每一级的输出都介于上一级输出的最小值和最大值之间,也就是说整个滤波器的中间状态不会超过输入信号的幅度范围。这意味着什么?意味着用int32做状态量时,只要输入不超过int32范围,中间过程永远不会溢出,不需要做饱和钳位,不需要额外检查。这种天然的内插有界特性是直接高阶IIR结构给不了的。
3. 移位替代除法的核心玩法:系数选型与误差控制
3.1 系数为什么可以“不够精确”
如果是做精密仪器,滤波器的截止频率误差当然越准越好。但绝大多数嵌入式的信号调理场景,我们并不需要截止频率精确到小数点后几位。传感器信号的频谱通常是一个范围,噪声和有用信号之间也往往有足够的分频距离。
举个例子,我需要约10Hz的三阶低通,采样率1kHz,那么单节截止频率需要约19.9Hz,对应的(\alpha):
[ \alpha = \frac{2\pi \times 19.9}{1000} \approx 0.125 ]
0.125正好等于1/8。这就太舒服了,在数字域里乘以1/8就是右移3位。那如果我要的不是19.9Hz,而是比如21Hz呢?对应的(\alpha \approx 0.132)。这个数很不“整”,怎么办?
工程上的做法是:主动放松对截止频率的精度要求,在(\alpha)可选值里挑一个最接近的。0.132附近的2的幂次是1/8=0.125,误差大约5.6%,反映到截止频率上就是从19.9变成21Hz左右。对于一个本身就有容差范围的模拟前端来说,这个偏差完全可以接受。
这就是移位替代除法的第一层逻辑:先用系数选型把(\alpha)凑成2的负幂次,把除法彻底消掉。
3.2 不是2的幂次也能用移位加一次加法逼近
但有的时候,我需要更精细地控制截止频率,又不想引入除法。这时候可以用组合逼近,也就是多个2的幂次做加减。
举例,我希望(\alpha)更接近0.113,可以用:
[ \alpha = \frac{1}{8} - \frac{1}{64} = 0.125 - 0.015625 = 0.109375 ]
在代码里就是:
*state += (diff >> 3) - (diff >> 6);对应的实现开销只是两次算术右移和一次减法。如果觉得0.109375精度还不够,可以用更大分母的组合,比如:
[ \alpha = \frac{30}{256} = 0.1171875 ]
而(30/256 = (32 - 2)/256),对应代码:
*state += ((diff << 5) - (diff << 1)) >> 8;本质就是先用移位组合算出分子,再一次右移做归一化。这种方法的截断误差取决于最终右移的位数,位数越多误差越小。
3.3 通用方案:小常数乘法加移位
如果觉得组合逼近还是不够灵活,那就用最通用的定点写法:
[ \alpha = \frac{K}{2^N} ]
滤波递推变成:
*state += (diff * K) >> N;这就是把除法改写为“乘小常数+右移”。对于(\alpha=0.113),取N=16:
[ K = 0.113 \times 65536 \approx 7406 ]
误差已经小到0.0016%。在Cortex-M0这类没有硬件除法指令的核上,一次整数乘法(MULS)通常只要1个周期,比软件除法快一个数量级以上。
这里有一个int32溢出边界必须注意:如果状态用int32_t,diff也在±65536以内,那么K必须满足:
[ |diff| \times K < 2^{31} ]
即K约小于32768,对应(\alpha < 0.5)。正常RC低通的(\alpha)不会超过0.5,所以这个限制实际上是天然满足的。
三种实现方式的对比,我整理了一个表格:
| 实现方式 | 精度 | 执行开销 | 代码量 | 适用场景 |
|---|---|---|---|---|
| 纯移位(diff >> k) | 低,系数只能取2的幂次 | 最低,每级1次右移 | 极小 | 对截止频率不敏感、追求极致速度 |
| 移位+加减组合 | 中,取决于组合项数 | 低,每级几次移位和加减 | 较小 | 需要一定系数灵活性又不想用乘法 |
| 乘法+移位(diff*K)>>N | 高,16位系数精度 | 中等,每级1次乘法+移位 | 小 | 通用场景,系数按需计算即可 |
实际项目中,我大部分时间用的是第三种,第二种偶尔用来处理极端精简的8位MCU场景。本文后续代码以第一种和第三种为主给出。
4. 定点化实现:完整代码与性能对照
4.1 数据结构与初始化
三个一阶节级联的最小状态就是3个int32_t。如果你选择了三个节系数不同,可以把每个节的shift或者K都保存在结构体里:
#include <stdint.h> typedef struct { int32_t s1, s2, s3; /* 三级状态量 */ } rc3_state_t; typedef struct { int32_t k1, k2, k3; /* 每级系数K,配合右移位数使用 */ uint8_t n1, n2, n3; /* 每级右移位数 */ } rc3_cfg_t;纯移位版就不需要cfg结构体了,三个shift如果相同,直接用一个宏定义即可。初始化一定不要省:
void rc3_init(rc3_state_t *st) { st->s1 = 0; st->s2 = 0; st->s3 = 0; }关于初值为什么不能是0,后面专门讲,这里先按标准做法来。
4.2 纯移位版核心滤波函数
先给一个三级同系数的纯移位版本,也是我最常用的形式:
static inline int32_t rc1_shift(int32_t x, int32_t *state, uint8_t shift) { int32_t diff = x - *state; *state += diff >> shift; return *state; } int32_t rc3_shift_process(rc3_state_t *st, int32_t x) { int32_t y1 = rc1_shift(x, &st->s1, 3); int32_t y2 = rc1_shift(y1, &st->s2, 3); return rc1_shift(y2, &st->s3, 3); }这里以(\alpha=1/8)为例。整个三级滤波跑完,只有3次减法、3次算术右移、3次加法,中间没有任何乘除法。对于有硬件除法指令的MCU来说优势不明显,但M0这类核,软件除法的开销往往是几十个周期,而右移只有1个周期,差距极其明显。
如果三个节希望用不同的shift,那就把cfg结构体传进去:
static inline int32_t rc1_shift_cfg(int32_t x, int32_t *state, uint8_t shift) { int32_t diff = x - *state; *state += diff >> shift; return *state; } int32_t rc3_shift_process_cfg(rc3_state_t *st, const rc3_cfg_t *cfg, int32_t x) { int32_t y1 = rc1_shift_cfg(x, &st->s1, cfg->n1); int32_t y2 = rc1_shift_cfg(y1, &st->s2, cfg->n2); return rc1_shift_cfg(y2, &st->s3, cfg->n3); }如果cfg里三个shift都等于3,那和上面的版本效果完全一样。
4.3 乘法+移位版核心滤波函数
当需要自由选取(\alpha)时,用这个版本:
static inline int32_t rc1_mul(int32_t x, int32_t *state, int32_t k, uint8_t shift) { int32_t diff = x - *state; *state += (diff * k) >> shift; return *state; } int32_t rc3_mul_process(rc3_state_t *st, const rc3_cfg_t *cfg, int32_t x) { int32_t y1 = rc1_mul(x, &st->s1, cfg->k1, cfg->n1); int32_t y2 = rc1_mul(y1, &st->s2, cfg->k2, cfg->n2); return rc1_mul(y2, &st->s3, cfg->k3, cfg->n3); }系数K的计算不放在中断里做,初始化时一次性算好:
static int32_t alpha_to_k(float alpha, uint8_t shift) { return (int32_t)(alpha * (float)(1 << shift) + 0.5f); }比如采样率1kHz、想要三阶整体-3dB约10Hz,则单节fc需要约19.9Hz,(\alpha \approx 0.125)。取shift=12,K就是:
[ K = 0.125 \times 2^{12} = 512 ]
如果shift取更大的数,比如16,K是8192。右移位数越大,量化误差越小,但要注意diff*K的溢出边界和乘法结果的精度。工程上shift取12到16之间比较合适。
4.4 与浮点版本的性能对比
我在一个Cortex-M0+内核的MCU上做过实测,主频48MHz,编译器GCC,开启了-O2优化。滤波输入是16位ADC值,输出也是16位。对比结果如下:
| 实现方式 | 单次三级滤波周期数 | 备注 |
|---|---|---|
| float全程计算 | 约800周期 | 软件浮点库,M0+无FPU,开销主要来自浮点加减乘除的库调用 |
| 整数除法版 | 约45周期 | 用int32配合除法运算符,每次除法调用运行时库 |
| 定点纯移位版 | 约12周期 | 只有加减和移位,无函数调用开销(已内联) |
| 定点乘法+移位版 | 约18周期 | M0+的32位乘法为单周期指令,整体开销仍很低 |
这个差距在实际应用里意味着什么?如果采样率是1kHz,浮点版每秒需要80万个周期,占48MHz CPU的1.7%;定点移位版每秒仅1.2万个周期,占0.025%。对于还要同时跑显示刷新、按键扫描、通信协议栈的MCU来说,1.7%看起来不多,但如果采样率提到10kHz,浮点版就是17%的CPU占用,而定点版只有0.25%。差距立刻变得不可忽视。
如果你用的MCU带有FPU,比如Cortex-M4F或M7带单精度FPU,那浮点和定点的差距会大幅缩小。但即便如此,定点版仍能节省寄存器压力、减少中断延迟,而且滤波结果完全确定性,没有浮点环境依赖。
5. 实测数据与边界条件:不是“能跑”就完了
5.1 恒定输入下的截断死区实测
这套定点实现最隐蔽的问题就是截断。以纯移位版为例,diff >> 3在diff为正数时是向零截断的。比如输入恒定32767,state从0开始往上追,当state追到32760时,diff=7,7>>3=0,state就不再变化。
这意味着输出会停在32760,而不是32767。差值是7个LSB,等于(\alpha)分母8减1。三阶级联后,每级都有类似死区,最终输出可能停在比输入低20多个LSB的位置。这个偏差相对16位满量程来说只有0.06%左右,很多人觉得无所谓,但对于称重、精密测温这类需要看绝对精度的应用,这就是不可接受的直流偏差。
解决方式有两个。第一个是接受它,把滤波结果的直流偏置当作固定offset在系统标定里扣除。第二个是改成带舍入的移位:
static inline int32_t rc1_shift_round(int32_t x, int32_t *state, uint8_t shift) { int32_t diff = x - *state; int32_t half = 1 << (shift - 1); if (diff >= 0) *state += (diff + half) >> shift; else *state -= ((-diff + half) >> shift); return *state; }这个版本把死区从7个LSB缩小到4个LSB以内,且正负对称。代价是多一次分支判断,在M0上大概多2到3个周期,大部分场景值得。
5.2 负数右移的截断方向坑
C语言里对负整数进行右移运算,符号位如何处理是由编译器实现决定的。主流嵌入式编译器(GCC、ARMCC、IAR)对带符号整数都采用算术右移,也就是右边补符号位,等价于向负无穷方向截断。而整数除法/在C标准里是向零截断。二者在负数时结果不同:
-7 >> 3得到的是-1-7 / 8得到的是0
如果代码里原来是除法,直接机械替换成右移,负数路径的直流偏置方向就反了。我在移植一个音频滤波算法时吃过这个亏,输出信号在负半周出现一个固定偏置,排查了好一阵才发现是截断方向不一致导致的。所以编写定点滤波器时,对负数的截断方向务必心里有数。上面的舍入版本就同时在正负两个方向做了对称处理,也是治这个毛病的。
5.3 阶跃响应与稳定时间
RC低通是过阻尼系统,理论上阶跃响应不会过冲。三阶级联后,阶跃响应变成S形爬升,稳定时间比单节更长。
单节RC的时间常数是(\tau = RC = \frac{1}{2\pi f_c})。三节相同fc级联后,整个系统约等效于3(\tau)的时间常数。以采样率1kHz、(\alpha=1/8)为例,单节fc≈19.9Hz,(\tau \approx 8ms),三节等效时间常数约24ms,从阶跃跳变到输出稳定在输入值的98%以内,大约需要5到6个等效时间常数,也就是120到150ms。
这个特性在实际应用里的影响是:如果传感器的物理量本身变化较快,三个节会造成可感知的滞后。所以选择极低截止频率时,要有意识地评估信号本身的动态范围是不是足够慢。
5.4 截止频率的验证方法
我验证滤波器截止频率的方式很简单:用DAC或PWM输出一个扫频正弦波作为输入源,从1Hz扫到100Hz,在另一个ADC通道或者串口把滤波前后的信号幅度记录下来。扫频步长要细,尤其在截止频率附近。
理论上,对于采样率1kHz、三级shift=3的实现:
- 输入1Hz正弦波,输出衰减约0.01dB,基本无衰减
- 输入10.1Hz,整体衰减应为-3dB,输出幅度约为输入的0.707倍
- 输入50Hz,根据前面计算,衰减约-26dB,输出幅度约为输入的0.05倍
实测值和理论值之间的误差主要来自信号源本身的谐波和量化噪声。如果偏差超过1dB,优先检查输入信号的幅值是否足够大,以及ADC采样是否稳定。
6. 踩坑记录与优化心得
6.1 上电初值不能是0,最好是第一次采样的值
滤波器的初始状态如果置0,而输入信号是3.3V传感器输出的2.5V直流偏置,那么滤波输出需要从0慢慢爬到2.5V附近,按上面的三级RC时间常数,可能要上百毫秒甚至更久。这段时间里,下游的控制或计算逻辑拿到的是明显偏小的滤波值,可能触发误判。
我的做法是初始化时不把state清零,而是把第一次进来的采样值直接赋给所有state:
void rc3_init_with_value(rc3_state_t *st, int32_t x) { st->s1 = x; st->s2 = x; st->s3 = x; }这样开机第一拍滤波输出就等于输入,跳过爬升过程。如果是周期性休眠唤醒的设备,每次唤醒后重新初始化时也建议带上一个参考值,否则每次唤醒都会重新经历一次爬升。
6.2 系数过小导致的更新缓慢与“卡顿感”
为了追求更强的滤波效果,很容易把(\alpha)选得很小,比如1/256。这时候每个采样点状态变化幅度最多是输入差值的1/256。如果输入变化较快,滤波输出看起来就像是一格一格地缓慢爬动,体感非常“卡”。
这里要算一笔账:(\alpha=1/256)、采样率1kHz时,等效截止频率:
[ f_c = \frac{\alpha F_s}{2\pi} = \frac{0.0039 \times 1000}{6.28} \approx 0.62Hz ]
这种滤波器适合处理变化极慢的温度、压力信号。如果盲目套用在音频或振动信号上,那就是灾难性的滞后和“糊感”。
6.3 中断上下文调用时的注意点
滤波函数修改的是结构体里的状态量。如果主循环和中断都调用同一个滤波器实例,就会存在临界区问题。最简单的做法是让滤波函数只在中断里执行,主循环只读取输出结果。如果确实需要多上下文调用,可以在中断入口用局部变量做一次快照,滤波完毕再写回。不要试图靠编译器优化来解决重入问题,踩过一次坑的人都知道。
6.4 高吞吐场景的进一步优化空间
如果采样率很高,比如音频场景的48kHz,每个采样点都调三个内联函数,调用开销虽然已经很小,但还是能挤出一点空间。可以改成一次处理一个块:
void rc3_process_block(rc3_state_t *st, const int32_t *in, int32_t *out, uint32_t n) { int32_t x, y1, y2, y3; for (uint32_t i = 0; i < n; i++) { x = in[i]; y1 = rc1_mul(x, &st->s1, K1, SHIFT1); y2 = rc1_mul(y1, &st->s2, K2, SHIFT2); y3 = rc1_mul(y2, &st->s3, K3, SHIFT3); out[i] = y3; } }编译器对固定次数的循环会自动做一定的展开优化,块处理比逐样本调用省掉循环开销。再配合DMA搬运ADC数据,CPU负载还能继续压低。
6.5 多级不同系数的自由度
前文例子都是三级相同系数,但实际工程里可以让每一级取不同shift,形成非对称响应。例如第一级取1/4用来快速吸收高频毛刺,后两级取1/16用来平滑慢变信号。这样整体特性不再是简单的三阶,而是更接近“一个高通抑制加两个低通平滑”的组合。你需要对各级截止频率分别计算,然后用前面的组合公式自由选型。这组自由度是没有写在教科书里、但实际操作非常好用的灵活性。
我在实际项目中长期使用三级RC定点滤波器,从8位MCU到32位M0+都跑过这套思路。最初只是图省事用移位替代除法,后来发现它带来的确定性和低开销让整机的中断负载、功耗预算都更好做了。如果你是刚开始接触定点滤波,先把这篇文章里的第一种纯移位版本跑通,再逐步换成乘法+移位版本。等滤波效果满足要求后,再回来深究截断死区、初值处理和系数选型这些细节,会少走很多弯路。