1. 为什么嵌入式环境要跟浮点较劲纯整数方案的真实动机先抛个大多数人都踩过的场景小封装 MCU 上跑传感器采集ADC 出来的原始数据带噪声很自然地想到写一个二阶低通滤波器。顺手用 float 把差分方程抄了一遍编译下载后发现事情不对——这块芯片根本没有 FPU每个 float 乘法背后都是一长串软浮点库函数一个 1kHz 的采样中断里跑一次滤波要吃掉几百个周期还要看编译器脸色。更麻烦的是浮点库调用时机不可控中断里偶尔蹦出个大延迟采样节拍一乱滤波器性能肉眼可见地变差。“纯整数二阶滤波器”就是干这个的把滤波器系数从 (-1, 1) 区间的小数放大成整数用定点乘加代替浮点运算整个滤波过程只有整数乘法、整数加法和移位。没有软浮点库调用没有隐式函数跳转执行时间固定到几十个周期以内这在很多 M0/M0 单核小芯片上属于“降维打击”级别的优化。这套东西解决的核心问题有三个没有 FPU 的 MCU 上浮点滤波器运行时间不可控、开销大滤波处理必须和采样中断严格同步允许的抖动窗口很窄内存和代码空间受限能省则省每个栈字节都金贵。适合看这篇的读者基本是这几类做 ADC 信号调理的嵌入式工程师、写电机或电源控制环的固件开发者、用单片机做简单音频处理的爱好者以及刚接触定点 DSP、想把滤波器真正落到板子上的学生。文章不会绕理论绕太久重点在怎么算系数、怎么定标、怎么写代码、上线前怎么验证。纯整数方案也不是没有代价。它的精度受字长限制动态范围不如 float系数需要经过量化搞不好会让滤波器的极点移动、甚至越过单位圆导致发散累加过程还可能溢出。这些坑后面每一处都会展开讲。2. 二阶 IIR 滤波器的结构为什么嵌入式里默认选直接 I 型2.1 差分方程与三种直接结构二阶 IIR 滤波器更常见的叫法是 biquad也就是双二阶滤波器。它的输入输出关系可以用一条差分方程描述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 控制极点。滤波器的频率响应、稳定性、群延迟全部由这五个系数决定。同一个差分方程在代码里可以翻译成好几种结构最常见的是三种直接 I 型DF1把输入历史和输出历史都存下来按顺序乘加。状态变量是 x[n-1]、x[n-2]、y[n-1]、y[n-2]物理意义非常直观——“前两个输入、前两个输出”。直接 II 型DF2先算递归部分再算前馈部分状态变量只有两个省内存但内部节点动态范围大定点实现时容易溢出。转置直接 II 型DF2T状态更新公式和 DF2 反过来稳定性更好是很多专业 DSP 库的默认选择但状态变量的物理含义不直观调试时容易懵。嵌入式定点开发我推荐直接 I 型。原因是它的每个状态变量都对应着真实的信号历史溢出分析和排查问题都能直接套用信号的幅值范围而 DF2、DF2T 的内部状态可能比输入输出大得多字长设计要额外小心。省那 4 字节内存往往不值得那个调试成本。2.2 系数从哪来双线性变换与 RBJ 计算表二阶滤波器的系数不是拍脑袋写的一般路径是从模拟原型出发经过双线性变换映射到数字域。双线性变换的核心思想是把模拟频率轴“卷”到数字频率轴保证了模拟滤波器的稳定性在变换后依然成立。这里有个关键概念叫频率预畸变。模拟频率和数字频率的关系不是线性的如果直接拿模拟截止频率去套数字滤波器实际截止点会偏移。预畸变的公式是ωa 2 * fs * tan(ωd / 2)其中 ωd 2π * fc / fsfc 是想要的数字截止频率fs 是采样率。所有正经的系数计算工具包括音频行业广泛使用的 RBJ EQ 手册公式内部都隐含了这一步。我平时算系数很少手推传递函数直接用 RBJ cookbook 的低通套路ω0 2π * fc / fs α sin(ω0) / (2 * Q) 低通 b0 (1 - cos(ω0)) / 2 b1 1 - cos(ω0) b2 (1 - cos(ω0)) / 2 a0 1 α a1 -2 * cos(ω0) a2 1 - α顺手把高通、带通、陷波也列在下面做嵌入式信号调理常用。类型b0b1b2a1a2低通(1-cosω0)/21-cosω0(1-cosω0)/2-2cosω01-α高通(1cosω0)/2-(1cosω0)(1cosω0)/2-2cosω01-α带通α0-α-2cosω01-α陷波1-2cosω01-2cosω01-α注意上表还没归一化。计算完要统一除以 a0让分母常数项变成 1才能得到最终可用的 b0、b1、b2、a1、a2。这一步在浮点世界里只是除一下到了整数世界直接影响定标结果。2.3 一组可以直接抄的系数示例理论公式看着干巴巴直接给两个我几乎每版项目都会用到的例子先跑通代码再说别的。例一采样率 1kHz截止 50HzQ0.707这就是二阶 Butterworth 低通。代入公式算完并归一化浮点系数是b0 0.02008 b1 0.04017 b2 0.02008 a1 -1.56105 a2 0.64135例二采样率 500Hz截止 10HzQ0.707。截得很低系数值更“贴近”极点边界正好用来后面演示量化误差b0 0.00362 b1 0.00724 b2 0.00362 a1 -1.82270 a2 0.83718记住这两组数的量级规律系数基本都是小于 1 的小数但 a1 的绝对值很容易超过 1甚至逼近 2。这个事实直接决定了我们后面选 Q 格式时不能盲目用 Q15。3. 纯整数化实现Q 格式、累加器与 C 代码全解3.1 Q 格式是什么把小数“放大”成整数再算浮点滤波器和整数滤波器的本质差异就是小数怎么表示。Q 格式的思路是用一个整数加一个固定缩放因子来表示小数。Q14 的意思是每个系数都乘上 2^14 16384存成整数实际值等于整数除以 16384。选 Q 格式的核心矛盾在于动态范围和精度Q 值越大步进越小精度越高但能表示的数值范围越小Q15 能表示的最大值接近 1.0放不下 a1可能到 1.56 甚至 1.9Q14 最大值接近 2.0覆盖绝大多数稳定 biquad 的分母系数。所以嵌入式定点 biquad 的甜点位是 Q14系数精度约 6.1e-5对应量化噪声在 -84dB 量级对 16 位传感器数据绰绰有余同时能容纳 a1 接近 2 的窄带滤波器。以例一低通滤波器为例把五个系数乘 16384 再四舍五入系数浮点值Q14 整数b00.02008329b10.04017658b20.02008329-a11.5610525576-a2-0.64135-10508这里有个不少人第一次写会蒙圈的符号约定。差分方程里分母项是减 a1y[n-1] 减 a2y[n-2]在代码里为了避免减法、也为了利用乘加的流水线我习惯把 c1 -a1、c2 -a2 直接存进结构体。这样递推式就变成了标准的五项乘加acc b0*x b1*x1 b2*x2 c1*y1 c2*y2例一中存入结构体的 c1 是 25576c2 是 -10508。这个“取负号再存”的做法是很多定点 DSP 库通行的约定理解以后查问题会顺畅很多。3.2 系数定标、舍入与量化后的稳定性检查浮点系数转整数最土也最稳妥的做法是四舍五入round(coeff * 16384)。不要用截断截断会引入系统性偏置直流增益很容易飘掉 1%-2%。量化以后必须做一次稳定性复核。二阶 IIR 的极点都在单位圆内才稳定量化后的判定条件可以总结成两条直接在整数域上验证1. |c2| S S 16384 2. |c1| S - c2拿例一验证|c2| 10508 16384满足|c1| 25576S - c2 16384 - (-10508) 2689225576 26892满足。滤波器稳定。例二就有意思了量化后分子求和是 59 119 59 237分母等效值是 S - c1 - c2 16384 - 29870 13717 231直流增益变成了 237/231 ≈ 1.026比理论值偏了大约 2.6%。这就是系数量化带来的直流增益误差。大多数传感器场景下这点偏差无所谓但如果后面接的是高精度测量链路就要注意预先补偿或者在系数微调时手动找一组让直流增益更接近 1 的整数组合。3.3 一套可直接移植的 C 实现完整实现放出来直接建工程就能跑。输入输出都是 int16_t系数用 Q14累加器用 int32_t#include stdint.h typedef struct { int16_t b0, b1, b2; /* Q14 分子系数 */ int16_t c1, c2; /* Q14 分母系数注意已取负号 */ int16_t x1, x2; /* 历史输入 */ int16_t y1, y2; /* 历史输出 */ } biquad_q14_t; int16_t biquad_df1_process(biquad_q14_t *f, int16_t x) { int32_t acc; acc (int32_t)f-b0 * x; acc (int32_t)f-b1 * f-x1; acc (int32_t)f-b2 * f-x2; acc (int32_t)f-c1 * f-y1; acc (int32_t)f-c2 * f-y2; /* Q14 还原加半 LSB 做四舍五入再右移 14 位 */ acc (1 13); acc 14; /* 状态移位与输出限幅 */ f-x2 f-x1; f-x1 x; f-y2 f-y1; if (acc 32767) acc 32767; else if (acc -32768) acc -32768; f-y1 (int16_t)acc; return f-y1; }结构体的 x1、x2、y1、y2 就是滤波器状态。处理完当前样点后把状态往后挪一拍等价于数学表达式里的 z^-1 延迟。整段代码没有除法、没有浮点、没有条件分支依赖数据真正做到了定时确定。初始化的时候要把整个结构体清零尤其是 y1、y2 不能是随机值否则上电头几百个样点会吐出吓人的毛刺。初始化代码一句话biquad_q14_t lpf {329, 658, 329, 25576, -10508, 0, 0, 0, 0};3.4 累加器溢出分析什么时候必须上 64 位int32_t 累加器到底够不够是个必须较真的问题。每个乘积项最大是 int16 输入32767乘以 Q14 系数最多 32767单乘积最大约 1.07e9。五个乘积全取最大值相加约 5.4e9已经超过 int32 上限 2.147e9。所以理论上最坏情况确实可能溢出。但工程上要分场景。一个最常用的边界分析办法是算“三角形不等式上界”最大累加和 ≤ 输入满幅 × (|b0| |b1| |b2| |c1| |c2|)例一的系数绝对值之和 329 658 329 25576 10508 37400。如果输入是 12 位 ADC 出来的 ±2048最大累加和约 2048 × 37400 76.6e6离 int32 上限差两个数量级随便用。如果输入满幅是 16 位 ±32767最大约 1.23e9还在 int32 范围内但余量只有 40% 左右此时阶跃、谐振峰等瞬态响应可能把状态推到更恶劣的位置不再安全。我的实践规则是这样传感器场景信号来自 12 位或 16 位 ADC 但预处理后限制在 ±16 位以下int32 累加器有余量放心用输出一路级联多级滤波、或者输入是 16 位满幅音频信号、或者滤波器 Q 值很高导致 a1 接近 2直接换 int64 累加器稳妥任何拿不准的边界用后面第 4 节的仿真办法先做最坏情况扫描。int64 版本改动很小把 acc 类型换成 int64_t乘法的强制转换改成(int64_t)移位常数加LL后缀即可int64_t acc; acc (int64_t)f-b0 * x; acc (int64_t)f-b1 * f-x1; acc (int64_t)f-b2 * f-x2; acc (int64_t)f-c1 * f-y1; acc (int64_t)f-c2 * f-y2; acc (1LL 13); acc 14;3.5 性能账本M0 和 M3/M4 上的实测量级精度之外性能是整数方案的卖点。我按工程经验给一组量级参考不同工具链会有差异但数量级不会变。Cortex-M3/M4 有单周期的 32x32 乘 64 位乘法指令编译器会把整段乘加翻译成几条 SMLAL 类指令。一个 Q14 biquad 处理一个样点大概在 20-40 个周期主频 168MHz 的 M4 上换算不到 0.25 微秒48kHz 双声道音频也才占几个点 CPU余量非常大。Cortex-M0 没有 32x32 乘 64 位指令int64 累加会退化成库函数调用一个样点可能到 80-150 周期。但主频 48MHz 上依然不到 3 微秒跑 1kHz 采样绰绰有余。如果 M0 上用 int32 累加器周期数能砍到 40 以内代价是需要按 3.4 节的规则做满幅分析。作为对照同一块 M0 上软浮点跑一个 biquad我一个量产项目里实测量级在 5-15 微秒之间取决于编译器浮点库的优化程度。纯整数方案在无 FPU 芯片上的优势不是“好一点点”而是差了一个数量级。4. 系数生成、验证与上线流程4.1 用 Python 直接生成 Q14 整数系数手算系数容易算错我把生成 Q14 整数系数的脚本固化下来了每次设计新滤波器直接改参数。脚本实现 RBJ 低通公式不依赖 scipy任何环境都能跑import math def rbj_lpf(fc, fs, q0.7071, bits14): w0 2.0 * math.pi * fc / fs cw math.cos(w0) sw math.sin(w0) alpha sw / (2.0 * q) b0 (1.0 - cw) / 2.0 b1 1.0 - cw b2 (1.0 - cw) / 2.0 a0 1.0 alpha a1 -2.0 * cw a2 1.0 - alpha # 归一化 s b0 / a0, b1 / a0, b2 / a0, a1 / a0, a2 / a0 b0n, b1n, b2n, a1n, a2n s scale 1 bits b0i int(round(b0n * scale)) b1i int(round(b1n * scale)) b2i int(round(b2n * scale)) c1i int(round(-a1n * scale)) # 注意负号 c2i int(round(-a2n * scale)) return b0i, b1i, b2i, c1i, c2i, scale # 例一1kHz 采样50Hz 截止 print(rbj_lpf(50.0, 1000.0)) # 例二500Hz 采样10Hz 截止 print(rbj_lpf(10.0, 500.0))输出就是可以直接填进 C 结构体的初始化值。高通、带通用同样的套路把 2.2 节的公式对应替换即可。如果想设计四阶、六阶巴特沃斯手写公式很烦建议直接上 scipyfrom scipy.signal import butter fs 1000.0 fc 100.0 sos butter(4, fc, btypelow, fsfs, outputsos) print(sos)sos返回的每一行就是一个二阶节的 [b0, b1, b2, a0, a1, a2]再用上面的脚本逐行转成 Q14 整数级联使用。4.2 定点仿真两类误差都要看系数量化以后必须把“浮点设计”和“Q14 定点实现”摆在一起对比不能只看波特图理论值。两类误差需要分别验证幅度响应误差。量化后的滤波器在截止频率附近的增益可能偏零点几分贝窄带滤波器尤其敏感。我一般直接在仿真里生成白噪声分别过浮点滤波器和整数滤波器对比输出频谱。整数滤波器不要直接在 C 里跑先在 Python 里按 Q14 规则模拟一遍改参数成本为零。瞬态与极限环误差。给两个滤波器喂同一个阶跃信号对比稳定后的残差。如果整数输出在零输入时还残留小幅振荡那就是第 5 节要讲的极限环问题要在仿真阶段就发现不要拖到硬件上查。量化对极点位置的扰动有个简单经验Q14 的步进是 6.1e-5二阶节的极点半径偏移通常在 1e-4 量级对 Q 值不高的滤波器完全无感但 Q 值超过 10、带宽极窄的滤波器极点本来就压在单位圆边缘量化后可能越界。这种情况要么降 Q 值要么用 Q13 甚至 Q12 换取更大的系数表示范围。4.3 上线前必做的三步检查代码合入固件之前我的固定流程就三步第一步直流增益检查。算一下 (b0i b1i b2i) / (S - c1i - c2i)确认在可接受偏差内。这一步能拦住绝大多数系数符号搞反的错误。第二步阶跃响应与最坏情况扫描。输入按预期最大物理量换算成整数再叠加一个正弦扫频信号跑 10 万个样点观察 acc 是否有逼近 int32 边界的时刻。顺手把滤波器初始化状态设成满幅测启动瞬态会不会冲过头。第三步稳定性复核。用 3.2 节的两条不等式在整数域验证别只在浮点上验证。极点是否在单位圆内要以量化后的系数为准而不是以设计稿为准。5. 实战踩坑与排障实录5.1 案例一零输入下输出冒出的 ±1 LSB“鬼信号”某次用整数滤波器处理压力传感器数据把传感器断开、输入置零示波器上看输出仍有幅度 1-2 个 LSB 的小幅振荡频率还和滤波器谐振频率相关。排查半天是典型的极限环现象。极限环的成因在反馈路径的舍入误差。零输入时输出状态全靠上一拍的历史值通过 c1、c2 反馈回来。如果舍入方向总是让这个反馈量“偏向”同一边状态就不会彻底衰减到零而是停在某个极小幅度上自我维持。DF1 结构配合四舍五入极限环的幅度通常不超过 1-2 个 LSB所以在 16 位数据链路里常常被当成噪声忽略。处理方式我一般按场景选传感器场景1 个 LSB 的振荡远小于传感器本身噪声直接接受不加代码音频场景会变成可闻的本底噪声建议在反馈路径改用“朝零截断”magnitude truncation即 acc 为正向下取整、为负向上取整用牺牲微小直流偏差换极限环消失实在敏感的场合输出端加一个死区幅度小于 2 LSB 直接清零。5.2 案例二采样间隔抖动导致截止频率漂移另一个真实教训滤波代码放在主循环里用HAL_GetTick()判断“差不多到 1ms 了就跑一次”。功能看起来正常但实测截止频率和设计值对不上低通变成了“低通带一点锯齿”。原因是主循环被其他任务打断实际采样间隔在 0.9ms 到 1.3ms 之间随机抖动。IIR 滤波器所有系数都建立在固定采样率假设上间隔一变截止频率跟着变。正确做法是把滤波器挂在定时器中断或者 ADC 完成中断里保证调用间隔固定。如果主循环里干很多活就用 DMA 定时器触发 ADC转换完成中断里直接读数据、滤波、输出一气呵成。滤波器的时钟就是采样节拍节拍乱了一切白搭。5.3 案例三高 Q 窄带滤波器量化后不稳定做某个音频均衡功能时浮点仿真完全正常切到 Q14 整数实现后输出越来越大直到饱和。问题出在窄带高 Q 滤波器的系数本来就贴着稳定性边界。浮点设计时 a2 是 0.9971量化后变成 0.9976少了一位有效数字极点就挪出了单位圆。这类问题的排查方法很直接用 3.2 节的两条不等式去验整数系数一旦第二条不满足基本就是量化后的极点越界。解决路径有三条换 Q13 或 Q12牺牲一点精度换取系数表示范围把窄带滤波器拆成两个宽带二阶节级联或者改用 DF2T 结构加整节后缩放。我在这个案例里的最终方案是改用 Q13 再把 Q 值从 12 降到 9输出波形和浮点设计偏差小于 0.1dB稳定跑完整个老化测试。5.4 常见问题速查表现象原因处理输出满幅饱和或发散累加器溢出换 int64 累加器或压低输入幅值零输入仍有小振荡极限环反馈路径朝零截断或加死区截止频率不对采样间隔抖动滤波放到定时器/ADC 中断里窄带滤波器量化后发散极点越界换 Q13/Q12降低 Q 值或拆级联直流增益明显偏大/偏小系数截断误差改用四舍五入手动微调一两个 LSB上电开头出现大毛刺状态变量未清零结构体初始化时全部填 06. 把整数滤波器做成一个可复用的模块6.1 接口设计的两个原则实际项目里滤波器不会只有一个ADC 去噪一个、电压环一个、诊断通道一个。把 biquad 封装成独立模块接口设计我坚持两个原则第一状态和数据分离。系数字段在运行期不该被修改状态字段每拍都在变。初始化时用结构体赋值一次性写入系数处理函数只碰状态。这样多个滤波器实例共用同一组系数时可以在内存里只存一份系数表状态分开。第二处理函数不做任何 I/O。滤波函数只负责“输入一个样点、输出一个样点”调用者决定数据从哪来、送到哪去。这样做的好处是滤波函数天然可重入、可单元测试后续无论接 ADC 还是接通信接口都不用改。一个顺手的小设计是把处理函数做成内联的因为每个采样周期都要调用函数调用开销在极端场景下也值得省。6.2 级联二阶节四阶及以上滤波的实现方式需要更陡的衰减曲线时直接写高阶差分方程容易数值不稳定工程标准做法是把高阶滤波器分解成多个二阶节级联。以四阶 Butterworth 低通为例scipy 输出两行 SOS 数据每行生成一个 biquad_q14_t 实例按顺序处理biquad_q14_t sec1 { /* 第一节系数 */ }; biquad_q14_t sec2 { /* 第二节系数 */ }; int16_t process_4th(int16_t x) { x biquad_df1_process(sec1, x); x biquad_df1_process(sec2, x); return x; }级联时每一级的输出就是下一级的输入。中间信号可能超过 16 位范围所以中间传递我建议用 int32_t 或者把各级输入限制在 14 位以内。级联顺序也有讲究一般把 Q 值高、动态范围大的环节放在后面降低中间饱和概率。这个规律在仿真里可以验证。6.3 选型整数 IIR 和滑动平均到底怎么选常有人问我传感器去噪是不是用滑动平均就够了答案是看需求。滑动平均FIR 类没有反馈绝对不会不稳定线性相位也无失真但它衰减斜率固定对采样率附近的噪声抑制能力弱而且需要 FIFO 缓冲区。IIR biquad 用几个状态变量换来陡峭的滚降适合需要干净频带的场景代价是相位非线性和潜在稳定性问题。我的习惯是这个只压白噪声、不关心过渡带形状用滑动平均或一阶 FIR 省事需要在某个频点以下保留信号、之上彻底压掉比如 50Hz 工频之外的振动噪声就上 biquad。纯整数 biquad 这个组件我前后用了三四年踩过溢出、踩过极限环、踩过采样抖动最后沉淀下来的就是这套“公式算系数 Q14 定标 仿真验证 中断里跑”的流程。它最大的价值不在于省那几微秒而在于让滤波器的行为变得可预测同样的输入永远在同样的周期数里给出同样的输出这在控制环和实时系统里比任何精度提升都重要。最后分享一个我一直保留的土办法每一版系数合入固件前把第 4 节的 Python 脚本跑出的整数系数连同浮点参考值一起贴到注释里。三个月后回来改代码的人通常就是你自己看到那两行注释能省掉一整个下午的重新推导时间。