做滤波器落到C代码这件事，我前前后后帮朋友调过好几轮，最大的感触是：MATLAB里点两下出来的曲线很漂亮，但把系数搬进C语言那一刻，才是真正的开始。尤其是巴特沃斯滤波器，很多人以为butter(n, Wn)拿到b和a就完事了，结果往嵌入式里一跑，波形不对、相位乱跳、甚至直接发散。这篇内容我就拿一个完整的例子，把“设计指标 -> MATLAB求系数 -> 二阶节拆分 -> C代码实现 -> 验证排查”这条链路讲透，重点说说那些文档里不会写、但实操里一定会遇到的事情。适合刚接触数字滤波的嵌入式工程师、自动化方向的学生，以及任何想把MATLAB算法真正变成可运行C代码的人。
1. 滤波方案选型：为什么IIR偏偏用巴特沃斯
1.1 巴特沃斯 vs 切比雪夫/贝塞尔：工程选择的底层逻辑
数字滤波器家族里，IIR滤波器其实有不少分支：巴特沃斯、切比雪夫I型/II型、椭圆滤波器、贝塞尔滤波器。为什么工程里默认首选巴特沃斯？核心原因在于它的通带内幅度响应最平坦，没有纹波。
通带内有纹波意味着什么？比如你在采集一个传感器信号，原始信号里有一个50Hz的有效分量，切比雪夫滤波器在通带边缘可能会让你这个分量幅度忽大忽小，误差就直接带进测量结果。巴特沃斯不会这样，它在通带内是一条非常平缓的曲线，代价是过渡带比切比雪夫宽一些。如果你的应用是采集、控制和一般的数据分析，没有极其苛刻的频率选择性要求，巴特沃斯是最省心的选择。
我做过一个对比测试，同样是4阶低通滤波，巴特沃斯的-3dB点在100Hz，切比雪夫I型纹波1dB时，相同的过渡带宽度可以做到大约70Hz开始衰减，看起来切比雪夫更“陡”。但换来的结果是通带内幅度有正负0.5dB左右的起伏。对于闭环控制系统，这种起伏会在特定频率上变成增益误差，严重时影响稳定性。所以我的原则是：追求“平坦”和“可靠”，选巴特沃斯；追求“极限陡峭”，才去考虑切比雪夫或椭圆，而且要提前想清楚纹波代价。
1.2 指标分解：通带、阻带、过渡带，怎么从需求逆推阶数
滤波器的设计不是直接拍脑袋定一个“几阶”就完了，而是先有指标，再由指标反推阶数。经典指标有三样：通带上限频率、阻带下限频率、阻带最小衰减。
举个具体例子。一个数据采集系统，采样率fs = 1000Hz，需要提取100Hz以下的有效信号，同时希望200Hz以上的干扰至少被压掉40dB。那么：
- 通带频率：
fp = 100Hz - 阻带频率：
fr = 200Hz - 通带最大允许衰减：
Ap = 1dB（这个可以理解成我允许100Hz以内的信号幅度稍微有点损失，但控制在1dB以内，约等于12%的幅度偏差） - 阻带最小衰减：
As = 40dB（200Hz以上至少要衰减到1%）
在MATLAB里，这段设计就是两行代码的事：
fs = 1000; % 采样率 fp = 100; % 通带频率 fr = 200; % 阻带频率 Ap = 1; % 通带最大衰减(dB) As = 40; % 阻带最小衰减(dB) % 根据指标反推最小阶数 n 和实际截止频率 Wn [n, Wn] = buttord(fp/(fs/2), fr/(fs/2), Ap, As); [b, a] = butter(n, Wn, 'low');这里buttord返回的n就是满足指标的最小阶数。在这个例子里，算出来通常是5阶。注意：MATLAB里fp/(fs/2)和fr/(fs/2)使用的都是归一化频率，分母是奈奎斯特频率，也就是采样率的一半。这是因为数字滤波器的频率响应是周期性的，设计时只看0到fs/2这一段，把这段映射成0到1的归一化值，MATLAB内部才好统一处理。
不少人在这一步容易踩坑：直接传了100和200进去，没除以fs/2。这样出来的滤波器截止频率完全不是100Hz，因为MATLAB会把输入当作归一化值0到1范围来处理。归一化频率这个概念，是整个设计流程里第一个需要刻进脑子里的东西。
2. MATLAB里把系数和前因后果搞明白
2.1 归一化频率：最容易翻车的第一个地方
再单独说一次归一化，因为它太容易错了。数字滤波器工作时，每个采样点之间的时间间隔是1/fs，对应的最高可表示频率就是fs/2。比如采样率1000Hz，你最多能表示0到500Hz的信号，超过500Hz的信号会混叠成别的频率，这在信号采集阶段就该处理掉。
因此MATLAB里的butter函数要求截止频率用“相对于奈奎斯特频率的比例”来表示。就是：
Wn = fc / (fs/2);假如采样率1000Hz，想要100Hz截止，那Wn = 100/500 = 0.2。如果你的采样率是2000Hz，同样想要100Hz截止，Wn = 100/1000 = 0.1，完全不一样。
我实际遇到过这样一个问题：某个开发者从别人的工程里复制了一段滤波代码，参数还是按采样率1000Hz算的Wn = 0.2，但他的系统实际采样率是500Hz。结果Wn = 0.2对应的实际截止频率就成了0.2 * 250 = 50Hz，设计目标100Hz的信号被削掉了一截。所以无论你是自己写还是抄代码，第一步先确认采样率对不对。这个错误最隐蔽，因为它不会让程序报错，曲线看起来也正常，就是结果不对。
2.2 直接型b/a与SOS到底选哪个
butter(n, Wn, 'low')返回的b和a是直接型传递函数系数：
H(z) = (b0 + b1*z^-1 + b2*z^-2 + ...) / (1 + a1*z^-1 + a2*z^-2 + ...)5阶的时候还行，但如果阶数上到10阶、15阶，直接用这个a向量在C代码里写差分方程，数值敏感性会急剧上升。因为高阶IIR的直接型实现里，滤波器极点对系数误差非常敏感，你在C语言里用float类型存系数，小数点后几位一截断，极点可能就飘出单位圆了——滤波器直接不稳定，输出乱飞。
工程上更稳的做法是拆成二阶节级联（SOS，Second-Order Sections）。把一个高阶滤波器分解成好几个二阶滤波器串联，每一级的系数都小得多，量化误差的影响被限制在局部，不会因为高阶传递函数的极点集中而放大。
MATLAB里拿到SOS很简单：
[sos, g] = tf2sos(b, a);或者一步到位：
z = zpk(butter(n, Wn, 'low')); % 这个方法也行，但不如直接点更干脆的是直接用butter的SOS输出：
[sos, g] = butter(n, Wn, 'low', 'sos');butter本身就支持直接生成二阶节格式，不需要先转到b/a再转回来。这样生成的sos是一个L x 6的矩阵，每一行对应一个二阶节，格式是：
[b0 b1 b2 a0 a1 a2]a0通常等于1，C代码里可以不存，也可以存下来当校验用。g是整个滤波器的增益，通常在所有二阶节的前面统一乘一次。
比如前面那个5阶滤波器，拆出来后通常是3个二阶节（5阶需要3级，最后一级是一个一阶节加一个二阶节的等效拆分）。sos矩阵里每一行代表一个独立的二阶滤波器：
y[n] = b0*x[n] + b1*x[n-1] + b2*x[n-2] - a1*y[n-1] - a2*y[n-2]这个形式写C代码非常直观，后面我会给完整实现。
2.3 可视化验证：别急着写C
拿到系数后千万别立刻去写C，先在MATLAB里把频率响应确认一遍。我习惯用freqz先看基本幅度谱，再用fvtool打开图形化界面，GUI里可以做更详细的检查：
freqz(sos, 2048, fs);freqz第一参数可以直接传SOS矩阵，第二参数是计算频率响应的点数，点数越多越平滑，第三参数传采样率，这样横轴显示的就是Hz而不是归一化弧度。看几个关键点：
- 100Hz处幅度是不是大约在-3dB位置；
- 通带内（0到100Hz）幅度是不是平直，没有明显纹波；
- 200Hz处幅度是不是低于-40dB。
如果这几个点都对，再继续下一步。如果不对，先回设计阶段检查buttord的指标参数，不要在C代码里徒劳折腾。
fvtool还能看零极点图、群延迟、阶跃响应。零极点图值得看一项：所有极点是否都严格在单位圆内。有些系统边界情况下极点在单位圆上，理论上稳定，但量化后很容易越界，看到这种情况直接加一点阻尼，或者换更高精度存储。
3. 从MATLAB到C：最实用的代码化路线
3.1 手写差分方程：二阶节级联的C实现
拿到SOS系数后，我推荐手写C代码来做实时滤波，因为生成的代码简洁可控，没有一堆MATLAB Coder自动生成的中间变量。二阶节级联的C实现大概是这样的：
typedef struct { float b0, b1, b2; // 前向系数 float a1, a2; // 反馈系数（a0永远=1，不存） float z1, z2; // 状态变量 } Biquad; float Biquad_Process(Biquad* f, float x) { float y = f->b0 * x + f->z1; f->z1 = f->b1 * x - f->a1 * y + f->z2; f->z2 = f->b2 * x - f->a2 * y; return y; }这个实现对应的是转置直接II型（Transposed Direct Form II）。为什么选这个结构？因为它用两个状态变量z1和z2就完成了滤波，不需要额外数组缓存过去的输入和输出。在嵌入式DSP里，这种结构计算效率高，状态更新自然，而且是业内默认做法。
再串联所有二阶节，做一个完整的滤波器：
#define NUM_SECTIONS 3 // 5阶=3个二阶节 typedef struct { Biquad sections[NUM_SECTIONS]; float gain; // 总增益，MATLAB里的g } BiquadFilter; float BiquadFilter_Process(BiquadFilter* filter, float x) { float y = x * filter->gain; // 先乘总增益 for (int i = 0; i < NUM_SECTIONS; i++) { y = Biquad_Process(&filter->sections[i], y); } return y; }这里有个容易忽略的点：MATLAB的butter(..., 'sos')返回的sos矩阵里，g是总增益，可能在分解时已经乘进某些系数里，也可能没有。取决于MATLAB版本和调用方式，有的返回g，有的会把g合并到第一个二阶节。稳妥的用法是显式在外面乘一次，像我代码里写的那样。如果你不放心，可以先用MATLAB的sosfilt和一个正弦波验证一下，直接用系数整体滤波，再看输出幅度对不对。
初始化状态变量也很简单，就是把所有z1、z2清零。一般来说，全局变量在嵌入式C里默认清零，但为了可移植性和明确性，我还是建议写一个BiquadFilter_Init：
void BiquadFilter_Init(BiquadFilter* filter) { filter->gain = 1.0f; for (int i = 0; i < NUM_SECTIONS; i++) { filter->sections[i].b0 = 0.0f; filter->sections[i].b1 = 0.0f; filter->sections[i].b2 = 0.0f; filter->sections[i].a1 = 0.0f; filter->sections[i].a2 = 0.0f; filter->sections[i].z1 = 0.0f; filter->sections[i].z2 = 0.0f; } }然后从MATLAB里把系数抄到C数组里，建议用一个const数组统一管理：
static const float FILTER_SOS[NUM_SECTIONS][6] = { {1.000000f, 2.000000f, 1.000000f, 1.000000f, -0.400000f, 0.050000f}, // 示例值 // 真实值以MATLAB生成结果为准 };注意MATLAB的SOS顺序是[b0 b1 b2 a0 a1 a2]，C数组里存的时候a0一般默认1，可以不写；如果你非要保留6个元素对齐，那么代码里要跳过a0。
3.2 利用MATLAB Coder半自动生成
如果不想手写，或者滤波器阶数很高、结构复杂，可以用MATLAB Coder自动生成C代码。思路是写一个入口函数，然后在Coder App里配置输入输出类型。
先写一个可生成的函数，比如：
function y = filter_100hz(x) % x 是输入信号数组，y 是滤波后的输出信号 persistent sos g if isempty(sos) [sos, g] = butter(5, 0.2, 'low', 'sos'); end y = sosfilt(sos, x) * g; % 注意如果g已并入sos就别乘两次 end然后打开MATLAB Coder，选择这个源文件作为入口函数，指定输入x的类型（比如double(1x1000)或single(1x1000)），生成代码。生成出来的代码会自动包含二阶节状态管理等中间结构，虽然不如手写那么精简，但胜在不用自己抠细节，适合快速原型验证。
做这类生成时有两个经验：
persistent变量会被Coder静态化为函数内静态变量，上电后保留上次状态。如果你的应用希望每次调用都从清零状态开始，需要额外写一个reset函数，或者在C接口里手动清结构体。- 输入信号类型建议直接用
single生成，这样生成出的C代码就是float运算，在嵌入式上更友好。如果你在MATLAB里用的是double，生成代码会是double，在某些M0/M3内核上会带来额外计算开销。
手写和自动生成的适用场景，我做了个简单对比：
| 对比项 | 手写C代码 | MATLAB Coder生成 |
|---|---|---|
| 代码可读性 | 高，自己控制逻辑 | 中等，自动生成变量多 |
| 嵌入式适配 | 好，可精确控制内存 | 需额外配置，有时需要裁剪 |
| 开发速度 | 慢，要理解结构 | 快，写个入口就行 |
| 系数可维护性 | 纯手动填，容易抄错 | 自动嵌入，一致性高 |
| 适合场景 | 长期维护、资源受限 | 原型验证、快速交付 |
我的习惯是：先在MATLAB Coder里生成一版跑通算法，确认滤波效果正确；再手写一版干净代码做最终交付。这样两边一对比，至少能确认手写代码里的系数顺序和状态更新逻辑没有错。
3.3 把系数固化成头文件，避免手动抄错
手动抄系数最大的风险是抄错。多位小数在眼睛和手之间转一轮，出错概率非常高。我一般用一个小脚本直接生成头文件：
num_sections = size(sos, 1); fid = fopen('biquad_coeffs.h', 'w'); fprintf(fid, '#ifndef BIQUAD_COEFFS_H\n'); fprintf(fid, '#define BIQUAD_COEFFS_H\n\n'); fprintf(fid, '#define NUM_SECTIONS %d\n\n', num_sections); fprintf(fid, 'static const float FILTER_SOS[%d][6] = {\n', num_sections); for i = 1:num_sections fprintf(fid, ' {%0.8ff, %0.8ff, %0.8ff, %0.8ff, %0.8ff, %0.8ff},\n', ... sos(i,1), sos(i,2), sos(i,3), sos(i,4), sos(i,5), sos(i,6)); end fprintf(fid, '};\n\n'); fprintf(fid, '#endif // BIQUAD_COEFFS_H\n'); fclose(fid);这样生成的.h文件里系数是8位小数精度，直接能放进工程里。%0.8f这个精度对float来说其实已经超过单精度有效位数了，写入时C编译器会做一次合法截断。如果你打算用双精度，可以写到12位以上。注意头文件里的系数顺序和代码里的顺序一定要对应，最好用一种带a0列的完整6元素数组，保留a0只是为了对齐MATLAB的SOS格式，代码里跳过它就行。
生成完头文件后，你可以在C代码里加一个编译期断言，用_Static_assert或者直接在初始化函数里校验第一个系数不为0，防止编译器在数组初始化时出错没发现。这些细节在项目大后期排查时会省很多时间。
4. 常见坑与排查实录
4.1 系数量化后性能劣化，极点飘出单位圆
前面提过，高阶直接型系数容易量化敏感，但哪怕拆成二阶节，如果系数精度给得太低，一样会出问题。我做过一个实验：用MATLAB设计一个8阶巴特沃斯低通滤波器，把SOS系数保留4位小数写入C，在MATLAB里用同样的4位小数系数重新算频率响应，发现阻带衰减从设计值的60dB退化到了45dB左右，而且在某些频段出现了奇怪的隆起。
解决办法有两个方向：
一是提高系数存储精度。能用float的尽量用float，但打印系数时至少保留6到7位有效数字；如果芯片资源允许，直接上double，省心很多。现代主流MCU如果带FPU，double代价还是不小的，建议先用float验证，再用double对比。
二是在MATLAB里提前模拟量化。你可以把系数转成single类型后重新做freqz:
sos_s = single(sos); freqz(double(sos_s), 2048, fs);如果曲线和原始设计几乎重叠，说明单精度够用；如果差别明显，要么提高精度，要么换结构。这个验证成本极低，但很多人都跳过不做，等嵌入式里跑出问题才回来查，很浪费。
4.2 初始化状态量没清零导致开机瞬间乱跳
滤波器状态变量在第一次调用前必须是0。如果没清零，输出会有一段莫名其妙的瞬态。物理上可以理解成滤波器在上电瞬间“看到”了一个巨大的虚假输入。
这个坑常见于两种情况：一是用RTOS，滤波器实例是动态分配的内存，malloc出来的内存内容未知，直接当结构体用，状态变量全是随机数；二是有人在调试时手动改过状态变量，没复位就继续跑。
解决方法其实很简单：每个滤波器定义的地方，显式调用一次初始化函数。别依赖calloc或BSS段的自动清零，因为不同平台行为不一样，移植到别的平台可能就翻车。我在自己的代码里始终坚持写一个_Init函数，即使初始化为全0，也写清楚意图。
4.3 采样率写错：最隐蔽的错误
这个前面在归一化频率部分提过，但值得放在排查实录里单独说一次。症状是：滤波器在MATLAB仿真里完全正常，搬到C里以后，输出波形幅度小了，但看不出明显问题。用扫频仪（或者手动输入不同频率正弦波）检查时，发现-3dB截止点跟设计要求差了整整一倍。
八成原因是采样率不对。比如你设计时用的是fs=1000Hz，但嵌入式里ADC实际工作在fs=500Hz，那么归一化频率0.2对应的不是100Hz而是50Hz，低频截止点整个左移了。
排查方法也比较直接：在代码里打印fs和Wn的值，和MATLAB设计脚本对比一遍。如果要严谨，还可以在嵌入式端产生一个精确的模拟方波，通过计算滤波后输出波形过零点的位置来反推实际截止频率。这个方法的精度还可以，能帮你验证整个链路的一致性。
4.4 数据截断和中间变量溢出
实时滤波C代码里，每个二阶节输出对下一级来说就是输入。如果输入信号幅度接近满量程，而滤波器又有DC增益，中间节点的数值很可能超过1.0甚至更多，特别是在级联结构里，增益会逐级累积。如果你用的是16位定点处理器，这个问题会极其明显：中间变量一饱和，滤波结果就畸形了。
解决办法有二。一是在输入到滤波器之前做一次总的幅度归一化，把总增益g先乘在输入上，而不是放在每个节里。二是在每级之间做防饱和处理，比如用饱和指令__SSAT或者手动加if语句限幅。对浮点MCU来说，只要确认了总增益分摊合理，基本不会溢出；但对定点平台，这一条必须当成核心问题处理。
我在实际项目里还碰到过一个反例：一个16位定点DSP上实现8阶巴特沃斯，所有系数都用Q15格式存储，当时只检查了系数本身不超过1，忽略了中间状态可能超过1，导致输出在高频段时有时无。最后是在每个二阶节后面加了一个移位和限幅，才彻底解决问题。所以，定点平台上的滤波器，状态变量范围和中间节点范围必须单独计算，不要想当然。
4.5 用正弦扫频验证整个链条
所有的排查技巧汇总起来，最终还是要靠验证。我在交付一段滤波器C代码前，一定会建立一个“测试台”：在C代码里输入一个已知频率的正弦波，同时用MATLAB做同样的输入，然后把两者输出做对比。
具体操作分三步：
- 在MATLAB里生成一个100Hz正弦波，长度1秒，用
sosfilt滤波得到参考输出； - 把同一个100Hz正弦波的C数组导入嵌入式工程，用C代码滤波，得到待测输出；
- 把两组数据拉到同一张图里对比，或者计算误差均方根。如果误差小于千分之一，基本说明代码实现正确。
这样做最大的好处是：如果C端输出和MATLAB参考对不上，你立刻就知道不是算法问题而是实现问题。很多人在调试时纠结“滤波器是不是不稳定”，其实只要对比这两条曲线，问题一下就能定位到系数、状态初始化或数据类型上。
5. 几点个人体会
前面这些经验，大多是我在一次次的折腾里攒出来的。滤波这件事，理论上不难，但把它从MATLAB搬到真实嵌入式环境，中间隔着的不是代码，而是采样率、量化精度、状态初始化、定点溢出这些细节问题。我个人的体会是：设计参数一定要“先算后测”，不要靠肉眼调；验证一定要“MATLAB和嵌入式两条腿走路”，不要只信其中一边。
最后再分享一个小技巧：如果你做的滤波器是要长期用在某个项目里的，建议在MATLAB脚本里把设计、生成头文件、验证三件事写成一个完整的脚本保存下来。下次改采样率或者改截止频率，只需要改最上面的几个参数，一键重新生成，C代码几乎不用动。这个习惯看起来不起眼，但在项目迭代里能省下大量的重复劳动。做技术的人，最值得花的功夫就是让“可复现”这件事变成默认选项。
