简介面向数字信号处理学习者与嵌入式开发者的轻量参考代码包把 MATLAB 中 filtfilt 的双向滤波思路移植到 C 语言帮助在缺少 MATLAB 环境的工程里实现零相位、无延迟的滤波输出。压缩包共 1 个文件为单个 filt.c 源码体积约 2KB代码给出滤波器系数定义、输入输出数据处理以及前向与反向两次滤波的完整流程并涉及反向系数顺序应用、结果相乘累加、零填充与循环缓冲等边界处理技巧。读者可据此理解 IIR、FIR 滤波器由 MATLAB 到 C 的转换要点对照 MATLAB 原生 filtfilt 校验结果一致性排查相位失真与边界不连续问题适用于音频处理、生物医学信号分析等需要保持原始时序关系的场景。目前已有 382 人学习代码体量小便于直接嵌入或二次修改。1. 用 C 做零相位滤波到底难在哪filtfilt 这个名字做信号处理的人都不陌生它就是零相位滤波的经典实现先正向滤一遍把信号倒过来再滤一遍用两次相位抵消换回一条不落后的波形。麻烦在于MATLAB 和 Python 里一行调用的事落到 C 语言实现上得自己管二阶节级联、稳态初值、边缘延拓和缓冲区反转。这份资源就是干这件事的——一套纯 C、无第三方依赖的 filtfilt 实现附带示例和对比工具。它适合两类人一类是在嵌入式或上位机里做传感器数据处理需要实时或离线滤掉基线漂移、又不能让波形整体后移的工程师另一类是学生或转行开发者想在 C 里复现出和参考实现逐点对得上的滤波结果。压缩包里该有的都有滤波器设计、延拓、前后向滤波、CSV 输入输出示例编一下就能跑。2. 工程骨架一次 filtfilt 调用经过哪些文件解包之后看到的不是一个大杂烩 .c 文件而是按职责切开的几个模块。这种拆法在信号处理类的小型 C 工程里很常见好处是后面换成定点、换成 FIR、塞进 MCU都只动其中一两个文件。我一般拿到这类包先看头文件头文件看懂了整个项目的设计意图基本就清楚了。2.1 目录结构与模块职责典型组织方式是这样的filtfilt/ ├── include/ │ ├── filtfilt.h /* 对外唯一入口filtfilt() 与 PadType */ │ ├── biquad.h /* 单节二阶滤波系数结构与单步运算 */ │ └── sos_design.h /* 滤波器设计Butterworth 拆节 */ ├── src/ │ ├── filtfilt.c /* 延拓 正反两遍 裁剪主流程 */ │ ├── biquad.c /* DF2T 单步、稳态初值 */ │ ├── pad.c /* 奇/偶/常数三种边缘延拓 */ │ └── sos_design.c /* 双线性变换 预畸变生成系数 */ ├── examples/ │ ├── demo_sine.c /* 造一段含噪正弦滤完写 CSV */ │ └── compare_csv.c /* 读 CSV、滤波、再写 CSV供外部比对 */ ├── test/ │ └── test_filtfilt.c /* 脉冲响应、直流增益等断言 */ ├── data/ │ └── sample.csv /* 示例输入 */ └── Makefile这个结构里有两条线值得注意。第一条是biquad.c 和 filtfilt.c 的分离单节滤波只管一个样本进、一个样本出状态用两个 double 传filtfilt.c 只管信号层面的延拓、反转、裁剪完全不知道滤波器内部长什么样。第二条是sos_design.c 独立存在它属于「离线设计」跑一次拿到系数就固定下来运行时不需要三角函数。很多嵌入式项目恰恰卡在这里——把设计代码和运行代码混在一起结果链接时把 libm 也带进去了还多花几百字节 Flash。2.2 对外接口与数据流头文件是契约。看懂下面这三段剩下的代码基本是填空/* 单个二阶节a0 已归一化为 1所以只存 a1/a2 */ typedef struct { double b0, b1, b2; double a1, a2; } Biquad; /* 滤波器若干二阶节级联 一个总增益 */ typedef struct { const Biquad *sos; size_t nsections; double gain; } SosFilter; typedef enum { PAD_ODD 0, /* 奇对称延拓默认斜率连续 */ PAD_EVEN, /* 偶对称延拓值连续 */ PAD_CONSTANT, /* 常数延拓最简单也最容易在首尾起振 */ PAD_NONE } PadType; /* 返回 0 成功负数表示参数错误或内存不足 */ int filtfilt(const SosFilter *f, const double *x, size_t n, double *y, PadType padtype, size_t padlen);参数含义与改法如下参数作用常见取值改错了会怎样f-sos二阶节数组2 节对应 4 阶节数写多了读越界结果随机f-nsections节数15少写一节等于只滤了一半f-gain总增益一般 1.0输出整体缩放波形对但幅值错x / n输入与长度n ≥ 4×padlenn 太小直接返回错误y输出缓冲与 x 等长允许 xy 原地滤波内部已开临时区padtype延拓方式PAD_ODD改常数延拓后首尾容易出现台阶padlen延拓点数3 × 滤波器阶数给小了首尾振铃给大了短信号被吃光数据流是一条直线调用方把原始采样点交给filtfilt它先按padlen把两端撑开得到长度n 2*padlen的工作缓冲然后正向滤一遍写入临时缓冲接着把临时缓冲整体反转、用同一组系数再滤一遍最后反转回来掐掉两端的延拓段把中间n个点拷给调用方。这里有个容易看漏的细节第二遍滤波的初值不是零。按参考实现的做法第二遍也就是反向那一遍的初始状态应该由正向滤波结束时的末状态推出来或者干脆用正向输出的最后一个样本去算稳态初值。前者更贴近参考结果后者实现更简单、效果也够。资源里的实现走的是后一条我后面第 3 章会把两种写法的差异点出来。/* filtfilt.c 主流程骨架 */ int filtfilt(const SosFilter *f, const double *x, size_t n, double *y, PadType padtype, size_t padlen) { if (!f || !x || !y || n 0) return -1; if (padtype PAD_NONE) padlen 0; if (padlen * 2 n) padlen n / 4; /* 延拓不能吃掉整段信号 */ size_t m n 2 * padlen; double *buf (double *)malloc(m * sizeof(double)); double *tmp (double *)malloc(m * sizeof(double)); if (!buf || !tmp) { free(buf); free(tmp); return -2; } pad_edges(x, n, padlen, padtype, buf); /* 两端延拓 */ run_sections(f, buf, tmp, m); /* 正向 */ reverse_inplace(tmp, m); /* 翻转 */ run_sections(f, tmp, buf, m); /* 反向 */ reverse_inplace(buf, m); /* 翻回来 */ memcpy(y, buf padlen, n * sizeof(double)); /* 掐掉延拓段 */ free(buf); free(tmp); return 0; }m是工作长度padlen是单侧延拓点数。reverse_inplace就是首尾对调的交换循环没有技巧但要保证长度算对否则中间那个元素会被多换一次——这在奇数长度上不会出错在偶数长度上才显形属于典型的「跑小样例没事、跑真实数据翻车」的一类问题。3. 核心算法拆解二阶节级联与前后向两遍滤波这一章是全篇最该抄的部分。零相位滤波的结果对不对八成取决于三个东西滤波器用什么结构落地、初值怎么给、反向那一遍怎么处理。下面按调用顺序一个一个拆。3.1 二阶节为什么要用 DF2T 结构先把差分方程摆出来二阶节的标准形式是y[n] b0*x[n] b1*x[n-1] b2*x[n-2] - a1*y[n-1] - a2*y[n-2]直接照这个式子写代码叫直接 I 型需要保存 2 个历史输入和 2 个历史输出。转置一下把历史量合并成两个中间状态就得到直接 II 转置型DF2T只需要 2 个状态变量而且数值特性更好/* 单节 DF2Tz1/z2 由调用方持有跨样本保留 */ static inline double biquad_step(const Biquad *s, double xn, double *z1, double *z2) { double yn s-b0 * xn *z1; /* 当前输出只依赖 b0 和状态 */ *z1 s-b1 * xn - s-a1 * yn *z2; /* 更新状态 1 */ *z2 s-b2 * xn - s-a2 * yn; /* 更新状态 2 */ return yn; }xn是当前输入样本z1、z2是指向状态变量的指针。三个系数b0/b1/b2和两个反馈系数a1/a2都来自Biquad结构。注意yn先算出来再拿去更新状态顺序反了结果完全不同。选 DF2T 而不是直接 I 型理由很实际状态变量个数少一半且中间结果的动态范围更小。做高通或者窄带滤波器时直接 I 型的历史输出可能比输入大好几个数量级落到 float 上就是精度崩掉。这一点在低频段尤其明显属于不看不知道、一看吓一跳的坑。级联的时候把上一节的输出直接喂给下一节状态各自独立/* 把所有节串起来跑完一整段信号 */ static void run_sections(const SosFilter *f, const double *in, double *out, size_t n) { memcpy(out, in, n * sizeof(double)); for (size_t s 0; s f-nsections; s) { double z1, z2; /* 关键初值用首样本推稳态而不是清 0 */ biquad_steady_state(f-sos[s], out[0], z1, z2); for (size_t i 0; i n; i) out[i] biquad_step(f-sos[s], out[i], z1, z2); } if (f-gain ! 1.0) for (size_t i 0; i n; i) out[i] * f-gain; }in和out允许是同一块内存因为一开始就memcpy到out了后面是原地覆盖。f-gain单独留在最后乘是为了让每一节的系数都归一化到相近的量级避免某一节系数过大把精度吃掉。3.2 稳态初值不清零反而更准零初值是最直觉的做法也是最容易在首尾留下一段「爬坡」的做法。原因很简单滤波器默认自己从「静默」状态启动而真实信号在第一个采样点就已经有一个非零电平这个突变会被滤波器当成阶跃响应慢慢衰减掉衰减的过程就是你看到的那段翘起来的头。正确做法是让滤波器的内部状态和输入首样本「对上」。假定输入恒为常数 c输出也恒为 c·H(1)那么状态变量应该停在一个固定值上这个值就是稳态初值/* 输入恒为 c 时本节应处的内部状态 */ static void biquad_steady_state(const Biquad *s, double c, double *z1, double *z2) { /* H(1)直流增益等于 (b0b1b2)/(1a1a2) */ double h1 (s-b0 s-b1 s-b2) / (1.0 s-a1 s-a2); double y h1 * c; /* 输出稳态值 */ *z2 s-b2 * c - s-a2 * y; /* 由状态方程反解 */ *z1 s-b1 * c - s-a1 * y *z2; return; }推导过程就是令z1[n]z1[n-1]、z2[n]z2[n-1]、x[n]c、y[n]h1*c代回 DF2T 的两个状态方程解出来。c一般取当前段的第一个样本。有个边界要提一句分母1a1a2在高通或极低频截止时会接近 0除法会炸。稳妥的写法是加个判断分母绝对值小于 1e-12 时退化成清零初值宁可首尾有点翘也别输出 NaN——NaN 一旦进了状态变量后面整段数据全废而且这种污染在反向那一遍会被放大排查起来非常费劲。3.3 反向那一遍反转信号别反转系数前向跑完把结果整体首尾对调再跑一遍同样的滤波再翻回来。听起来像魔术原理其实一句话正向滤波给每个频率分量加了 φ(ω) 的相位延迟反向跑一遍相当于再加一个 -φ(ω)两者抵消总相位为 0。代价是幅频响应被平方这个代价后面第 5 章单独说。/* 原地反转长度 m */ static void reverse_inplace(double *v, size_t m) { for (size_t i 0; i m / 2; i) { double t v[i]; v[i] v[m - 1 - i]; v[m - 1 - i] t; } }m/2这个上界保证了偶数长度时中间两个元素正常互换、奇数长度时中间那个不动。写反成(m1)/2就会多换一次结果看起来「好像也对」只是中间样本错了肉眼几乎发现不了。这里最容易翻车的地方是有人觉得反向滤波要把 a1 取负、或者把系数数组倒过来。恰恰相反反向那一遍用的还是同一组b0/b1/b2/a1/a2唯一变化的是信号顺序和初始状态。系数一动滤波器就变成了另一个东西输出必然发散或者完全跑偏。我见过把a1取反导致后半段指数增长的例子波形图上看就像信号「自激」了。第二遍的初值资源里用的是run_sections内部的稳态初值逻辑也就是拿反转后缓冲区的第一个样本原来信号的最后一个样本去推。更贴近参考实现的做法是把正向结束时的状态直接接过来/* 更贴近参考实现的写法正向末状态作为反向初状态 */ static void run_sections_chained(const SosFilter *f, const double *in, double *out, size_t n, double *z1_io, double *z2_io, int forward) { memcpy(out, in, n * sizeof(double)); for (size_t s 0; s f-nsections; s) { double z1, z2; if (forward) { biquad_steady_state(f-sos[s], out[0], z1, z2); } else { z1 z1_io[s] ; z2 z2_io[s]; /* 接手正向的末状态 */ } for (size_t i 0; i n; i) out[i] biquad_step(f-sos[s], out[i], z1, z2); if (forward) { z1_io[s] z1; z2_io[s] z2; } } }z1_io/z2_io是按节数分配的数组长度等于nsections。两种写法在长信号上差异极小在只有几十个点的短信号上能差出可见的头部偏差。做对拍验证时如果发现只有头尾几十个点对不上基本就是这里的选择不同。3.4 系数从哪来双线性变换与预畸变资源里带了sos_design.c负责把「4 阶 Butterworth、截止 30 Hz、采样 500 Hz」这类人话翻译成系数。核心是二阶低通节的双线性变换加一步频率预畸变/* 二阶低通节fc 截止频率fs 采样率q 品质因数 */ void lpf_biquad(double fc, double fs, double q, Biquad *s) { double k tan(M_PI * fc / fs); /* 预畸变把模拟频率映射到数字域 */ double k2 k * k; double norm 1.0 / (1.0 k / q k2); s-b0 k2 * norm; s-b1 2.0 * k2 * norm; s-b2 k2 * norm; s-a1 2.0 * (k2 - 1.0) * norm; s-a2 (1.0 - k / q k2) * norm; } /* N 阶 Butterworth 拆成 N/2 个二阶节Q 值按极点角度分配 */ void butter_sos(double fc, double fs, int order, Biquad *sos) { int nsec order / 2; for (int i 0; i nsec; i) { double theta (2.0 * i 1.0) * M_PI / (2.0 * order); double q 1.0 / (2.0 * cos(theta)); /* 0.5412, 1.3066, ... */ lpf_biquad(fc, fs, q, sos[i]); } }tan(M_PI * fc / fs)这一步不能省。省掉之后截止频率会随fc/fs增大而明显偏移fc/fs 0.1时误差还能忍到0.3就已经偏得很离谱了。q是每个节的品质因数4 阶时是 0.5412 和 1.3066q越大那个节在截止点附近越「尖」级联后才是平坦的 Butterworth 响应。order只支持偶数奇数阶要额外补一个一阶节。资源里按偶数写够用要改奇数阶得在Biquad之外再加一个一阶结构或者把一阶节塞进二阶节并把b2/a2置零——后者省事但要注意稳态初值公式在一阶退化时仍然成立不用改。4. 编译与结果验证把曲线对上去算法看着对不代表跑出来对。这一章是把包里的东西变成可执行文件并且和参考结果做逐点比对。我一般的习惯是先跑示例看形状再对拍看数值两步都过才算真的跑通。4.1 用 Makefile 一次编出两个可执行文件包里给的 Makefile 是标准的分离编译写法-Iinclude指定头文件路径-lm链接数学库tan、cos需要CC ? gcc CFLAGS ? -O2 -stdc11 -Wall -Wextra -Iinclude LDFLAGS ? -lm SRC : src/filtfilt.c src/biquad.c src/pad.c src/sos_design.c OBJ : $(SRC:.c.o) all: demo compare demo: $(OBJ) examples/demo_sine.o $(CC) $^ -o $ $(LDFLAGS) compare: $(OBJ) examples/compare_csv.o $(CC) $^ -o $ $(LDFLAGS) test: $(OBJ) test/test_filtfilt.o $(CC) $^ -o $ $(LDFLAGS) ./test clean: rm -f $(OBJ) examples/*.o test/*.o demo compare test直接三条命令跑起来make clean make # 编出 demo 和 compare ./demo # 造一段含噪正弦输出 out.csv python3 check.py # 与参考实现逐点比对make test会跑一组断言主要检查三件事单位脉冲输入下输出是否近似偶对称零相位的直接证据、直流输入下输出是否稳定在同一个值稳态初值有没有生效、全零输入下输出是否恒为零状态有没有残留。CFLAGS里-Wall -Wextra别省。这类代码里最常见的告警是有符号/无符号比较size_t和int混用减出负数后变成巨大的正数循环直接跑飞。开了告警能省下大量排查时间。4.2 示例程序从 CSV 到 CSVcompare_csv.c的职责很单一读入一列浮点、调filtfilt、按原顺序写出一列浮点。这样外部任何工具都能拿来做对拍。#include stdio.h #include stdlib.h #include filtfilt.h #include sos_design.h int main(int argc, char **argv) { if (argc 4) { fprintf(stderr, usage: %s in.csv out.csv fc\n, argv[0]); return 1; } double fc atof(argv[3]), fs 500.0; /* 4 阶 Butterworth 低通拆成 2 个二阶节 */ Biquad sos[2]; butter_sos(fc, fs, 4, sos); SosFilter f { sos, 2, 1.0 }; /* 读数据一个文件里一列逐行扫 */ double *x NULL; size_t n 0, cap 0; FILE *fp fopen(argv[1], r); if (!fp) return 1; double v; while (fscanf(fp, %lf, v) 1) { if (n cap) { cap cap ? cap * 2 : 1024; x (double *)realloc(x, cap * sizeof(double)); } x[n] v; } fclose(fp); double *y (double *)malloc(n * sizeof(double)); size_t padlen 3 * 4; /* 3 倍阶数够用且不至于吃掉短信号 */ int rc filtfilt(f, x, n, y, PAD_ODD, padlen); fp fopen(argv[2], w); for (size_t i 0; i n; i) fprintf(fp, %.12g\n, y[i]); fclose(fp); free(x); free(y); return rc 0 ? 0 : 2; }padlen 3 * order是经验值阶数越高、过渡带越陡需要的延拓就越长。%.12g保留 12 位有效数字是为了让对拍时能看出 1e-10 量级的差异——用%f打印会把所有误差抹平看起来「完全一致」实际早就偏了。4.3 与参考实现对拍对拍是最省事的验证手段。手边有 Python 的话几行就够import numpy as np from scipy.signal import butter, sosfiltfilt x np.loadtxt(data/sample.csv, delimiter,) sos butter(4, 30, btypelow, fs500, outputsos) y_ref sosfiltfilt(sos, x, padlen12) y_c np.loadtxt(out.csv) err np.abs(y_c - y_ref) print(max abs err %.3e at %d % (err.max(), err.argmax())) print(rms err %.3e % np.sqrt((err ** 2).mean()))判读标准我通常这么定max abs err小于信号峰值的 1e-6rms err小于 1e-8就算通过。如果误差集中在首尾各padlen个点那基本是延拓方式不同PAD_EVEN对odd或者反向初值策略不同如果误差沿全程均匀分布那是系数精度或者节的排列顺序不同如果中段某个位置突然冒出一个大尖峰八成是缓冲区长度算错了。提示对拍时务必用同一组系数。两边各自设计滤波器即使参数写的一样浮点三角函数在最后一位上的差异也可能被级联放大导致误差看着有 1e-9 但怎么改都对不齐。还有一种不看数值的验证法扫频看群延迟。喂一段频率缓慢上升的正弦把输入输出叠在一张图上理想情况下两条曲线应该完全重合没有任何水平位移。如果输出整体右移说明相位没抵消干净如果幅度在某段被压掉一大截说明截止频率被迫平方效应拉低了那就是下一章要讲的事。5. 避坑与排查零相位滤波里最容易翻车的五个点这一章是我自己踩过和帮别人看过的坑按出现频率排序。每条都按现象、原因、解决三段写遇到问题可以先对着现象找。5.1 现象高通滤完后半段越来越大直到溢出原因几乎都是用了直接型b/a而不是二阶节级联。高阶滤波器的极点密集分布在单位圆附近直接型的系数对舍入误差极其敏感极点稍微一偏就跑到圆外系统变成不稳定。低频截止越低越明显因为极点越贴近 1。解决办法是老老实实拆成二阶节每个节单独跑节间用double传值。两个二阶节串联的数值稳定性远好于一个四阶直接型代价只是多一次循环。已经在用二阶节还出现这个现象的检查节的顺序——把 Q 值大的节放在后面通常更稳因为前面的节先把大幅度成分压下去了。5.2 现象首尾各有一小段明显翘起中间完全正常padlen给小了。前后向滤波的瞬态在两端各出现一次延拓长度不够滤波器还没「忘掉」启动瞬间的突变就进入正文了于是那一段翘起被保留在输出里。把padlen提到滤波器阶数的 3 倍以上。4 阶滤波器给 12 点起步8 阶给 24 点。同时也别给太大padlen * 2 n时整个信号都被延拓段挤没了代码里那个n / 4的兜底就是防这个。做短帧处理比如每帧 64 点时宁可把帧拼长一点再滤也不要把padlen压到 4 以下。顺带一提延拓方式也会影响首尾形状。奇对称延拓保持斜率连续绝大多数情况最自然偶对称延拓保持值连续在信号两端本来就有明显趋势时更合适常数延拓最省事但首尾容易出现台阶一般只在调试时用。5.3 现象对拍时只有前 20 个点对不上正向那一遍没问题问题出在反向的初值。零初值会让反向那一遍从「静默」启动在反转后的开头也就是原信号的尾部产生一段瞬态。如果你的实现里run_sections是统一用稳态初值那这一段会小很多但不为零。两个选择一是反向也用稳态初值拿反转缓冲区的首样本去推二是把正向结束时的状态接过来。前者实现简单后者更贴近参考实现。做产品时我更倾向第一种因为状态数组不用跨函数传代码干净做对拍工具时用第二种能压到 1e-10 以内。5.4 现象滤波器参数明明按 -3dB 设计的实际通带窄了一圈这是零相位滤波的固有代价不是 bug。前后向各滤一次幅频响应被平方原来 -3dB 的点变成了 -6dB阻带衰减翻倍通带纹波也翻倍。视觉上就是「滤得更狠了」。要补偿的话把设计截止频率往上抬。n 阶 Butterworth 的等效 -3dB 点约在单次截止频率的0.4142^(1/(2n))倍处反过来要把fc乘以0.4142^(-1/(2n))阶数等效截止/设计截止补偿系数20.801.2540.901.1160.931.0780.951.05想要等效 -3dB 落在 30 Hz4 阶就得设计在 33 Hz 左右。这个补偿不做也行但心里要有数否则和别人的结果对比时永远差那么一点。5.5 现象换成 float 后输出出现缓慢漂移或台阶状态变量用 float 累积误差。二阶节的状态是递归的误差会一直在里面打转double下看不出来float跑几万个点就开始有可见偏移尤其在低频截止很低的滤波器上。能上double就上Cortex-M4/M7 带 FPU 的话double只是慢一点精度收益很大。实在要用float两条路一是把信号分段每段重新用稳态初值启动把累积误差截断二是改用定点实现但定点要自己做缩放每个节的增益都得重新分配工作量不小。上 MCU 前先拿长数据在 PC 上跑一遍float版本看误差有没有超预算比烧进去再调省事得多。6. 二次开发换滤波器、改精度、搬上 MCU拿到这份代码之后真正花时间的往往不是读懂它而是把它改成自己项目要的样子。下面三条改造路线是我实际用过的按改动量从小到大排。第一条换成带通或者高通。sos_design.c里只给了低通加高通最简单——把lpf_biquad的分子系数做个频域翻转就行b0 b2 b1 b0这一类推导容易出错更稳的做法是按标准双线性变换重写一个hpf_biquad高通的原型是s^2/(s^2 s/q 1)代入s (1-z⁻¹)/(1z⁻¹)后分子变成(1-z⁻¹)²系数是b0 norm、b1 -2*norm、b2 norm分母不变。带通就是低通和高通级联但要注意两级的截止频率要拉开否则中间那段通带会被压得很薄。第二条改成块处理为上板做准备。现在的filtfilt一次性吃完整段信号内存是2*(n2*padlen)个 double。做离线分析没问题做实时就得改成「攒够一帧再滤」。这里有个取舍帧越长首尾延拓浪费的比例越小但延迟越大。我一般按帧长 8 * padlen来配4 阶滤波器padlen12帧长至少 96 点取 128 或 256 比较舒服。改的时候把filtfilt拆成filtfilt_begin / process / end三段中间那段复用同一块静态缓冲区避免每次 malloc。/* 分块版本的核心保留上一帧尾部作为本帧的延拓素材 */ typedef struct { const SosFilter *f; double tail[32]; /* 上一帧最后 padlen 个样本 */ size_t padlen, tail_n; } FiltfiltCtx; int filtfilt_block(FiltfiltCtx *ctx, const double *in, double *out, size_t n) { size_t m ctx-tail_n n; /* 拼上历史尾巴 */ double *buf (double *)malloc(m * sizeof(double)); if (!buf) return -2; memcpy(buf, ctx-tail, ctx-tail_n * sizeof(double)); memcpy(buf ctx-tail_n, in, n * sizeof(double)); /* 用拼好的段做一次完整 filtfilt取后 n 个点输出 */ double *y (double *)malloc(m * sizeof(double)); int rc filtfilt(ctx-f, buf, m, y, PAD_ODD, ctx-padlen); if (rc 0) memcpy(out, y ctx-tail_n, n * sizeof(double)); /* 更新尾巴保留本段最后 padlen 个样本 */ ctx-tail_n ctx-padlen m ? ctx-padlen : m; memcpy(ctx-tail, buf m - ctx-tail_n, ctx-tail_n * sizeof(double)); free(buf); free(y); return rc; }tail的长度取padlen就够了因为延拓只需要这么长。tail_n单独存是为了处理第一帧——那时候还没有历史长度是 0。这个写法不追求理论上的严格分段等价追求的是实现简单、边界不炸实测在 256 点帧长、4 阶滤波器下和整段滤波的差异在 1e-6 量级够用。第三条精度改造。如果目标平台只有单精度 FPU把double全换成float之后先在 PC 上跑一遍长数据对比。我的判据是max abs err / peak小于 1e-4 可以接受超过 1e-3 就得考虑分段重初始化。另外sos_design.c可以整个从固件里剔掉——系数在 PC 上算好直接写成常量数组编进去省下tan、cos和整个 libm 的开销。系数数组加上static const修饰还能进 Flash 而不是 RAM。最后说一个我自己的习惯。从那以后我每次改完滤波器系数或者延拓参数都会强制走一遍「脉冲 直流 扫频」三件套脉冲看对称性直流看稳态扫频看群延迟和幅频。三张图都正常了再谈上板。这三个用例加起来不到五十行代码但救回来的调试时间远超写它的成本。压缩包里的源码是可以直接编的示例数据也在改之前先在本地跑通一遍对拍对后面的改动心里才有底。希望帮到你。本文还有配套的精品资源点击获取