简介面向数字信号处理学习者与嵌入式开发者的双向滤波 C 语言实现示例重点把 MATLAB 中 filtfilt 的零相位滤波思路迁移到 C 环境解决单向滤波引起的相位失真与时序偏移适合音频处理、生物医学信号分析等对波形时间关系敏感的场景。压缩包总共 1 个文件主体为一份 C 源码整体约 2KB体量轻巧便于逐行阅读和移植内容围绕滤波器类型与系数设定、前向滤波、反向滤波、样本乘积累加以及边界条件处理展开并配有与 MATLAB 原生函数结果对照验证的思路可帮助读者理解 IIR、FIR 滤波在代码层面的差异。已有 382 人学习适合作为课程实验、算法复现或工程移植的参考底稿读者能据此梳理零相位滤波的完整实现脉络掌握从理论公式到 C 代码落地的关键方法。1. 一段对不齐的振动波形把我拉回了 filtfilt 的 C 实现前年帮某实验室处理一批设备振动数据加速度计采回来的曲线上叠着高频毛刺我随手写了个二阶 IIR 低通滤了一遍波形是好看了可峰值时刻相对原始信号整体往后挪了十几个采样点。两个通道做互相关算时延算出来的数里混进了滤波器本身的群延迟标定怎么都对不上。问题不在参数在于单向 IIR 滤波天生带相位滞后。filtfilt 就是冲着这件事来的正向滤一遍、反向再滤一遍相位响应一乘就抵消幅频响应变成原来的平方。代价是它非因果只能离线或分块跑。这套 C 实现适合手里只有 C 工具链、又不想为一次滤波把整个数值计算环境拖进来的场景——上位机离线分析、采集端回传后的二次清洗、嵌入式设备上的批量数据处理源码拉进来就能编。2. 解压后先看代码结构filtfilt 的两次滤波为什么要带状态拿到压缩包别急着编译先确认它把「单次滤波」「初值计算」「两次串联」拆成了哪几层。常见实现里这三件事各占一个文件混在一起的版本通常意味着初值被写死成 0那种版本在端点附近的表现会明显差一截后文会讲怎么判断。2.1 一份典型 C 实现的文件构成不同人打包的命名习惯不一样但函数职责跑不出下面这几类。看代码结构比看文件名可靠打开头文件扫一眼函数签名基本就能定位主流程在哪。文件职责关键函数filtfilt.c主流程延拓、正反两次滤波、去延拓filtfilt()iir.c单次 IIR 滤波通常用直接 II 型转置结构iir_df2t()coeff.c/h滤波器系数 b/a 与预计算好的初值 zi常量数组demo.c读文本数据、调用、写结果方便对拍main()Makefile/CMakeLists.txt构建脚本一般只依赖 libm—这里最容易出问题的是coeff.c。系数是谁设计的、用什么结构设计的、有没有做过归一化源文件里往往没有注释。如果a[0]不是 1后面所有公式都得重推第一次跑之前务必打印确认。2.2 相位抵消的数学前提幅频被平方了单次滤波的频响记作H(e^{jω}) |H|·e^{jφ(ω)}。时间反转在频域等价于取共轭所以反向再滤一遍得到|H|·e^{-jφ}两次串起来是|H|²·e^{j0}相位项彻底消掉这才是「零相位」的真正含义。它不是把延迟补回来而是让延迟结构上不存在。代价有两个都容易被忽略。第一幅频响应被平方了原来截止频率处 -3dB 的点滤完之后变成 -6dB实际通带比设计值窄想要同样的有效带宽得把设计截止频率往上抬一点。第二滤波器变成非因果输出第 n 点依赖第 n 点之后的数据实时流里用不了只能整段处理或者分块加重叠。选型上也因此有个基本判断IIR 阶数低、算力省但相位非线性严重是 filtfilt 收益最大的场景FIR 本身可以做成线性相位群延迟恒定直接时移补偿就够了硬套 filtfilt 只是把群延迟也消掉代价是阶数高、内存翻倍。振动、生理信号、地震波这类关心「峰值落在哪一时刻」的离线数据IIR filtfilt 是性价比最高的组合。2.3 零初值为什么会毁掉开头几十个点单次filter从状态全零起步等价于假设输入在 t0 一直是 0。如果x[0]本身不是 0这个跳变就是一次阶跃激励输出的前若干点会叠上一段按极点衰减的瞬态。极点越靠近单位圆衰减越慢有时候几百个点都散不掉。解决办法有两条成熟实现是两条一起用。一条是延拓在数据两端造出一段「比真实数据更早/更晚」的假数据让瞬态发生在假数据上滤完再切掉另一条是初值直接算出「如果这段信号从无穷远就一直存在此刻的状态该是多少」让滤波器从稳态起步。前者解决数据边界后者解决状态边界缺一个端点都会有残差。2.4 单次滤波的骨架直接 II 型转置不管上层怎么组织底层那一次滤波大概率是直接 II 型转置DF2T状态长度等于阶数数值上比直接 I 型稳中间变量也少。/* iir.c —— 单次 IIR 滤波直接 II 型转置结构 */ #include iir.h /* b、a 需已归一化(a[0]1)并补齐到同一长度 nz 长度 n-1原地更新 */ void iir_df2t(const double *b, const double *a, int n, const double *x, double *y, int N, double *z) { if (n 2) { /* 纯增益无状态可维护 */ for (int i 0; i N; i) y[i] b[0] * x[i]; return; } for (int i 0; i N; i) { double xi x[i]; double yi b[0] * xi z[0]; /* 输出只依赖旧状态 */ for (int k 0; k n - 2; k) z[k] b[k 1] * xi - a[k 1] * yi z[k 1]; /* 必须顺序更新 */ z[n - 2] b[n - 1] * xi - a[n - 1] * yi; y[i] yi; } }参数上要盯三件事。n是补零对齐之后的长度max(len(b), len(a))短的那个尾部补 0z的长度是n-1调用方负责分配和初始化函数内部不碰它的初值z是原地改写的所以同一份状态不能同时喂给两个线程。那个n2的分支不是凑数一阶增益环节没有状态漏掉它会读到长度为 0 的数组。内层循环不能写成z[k] ...的并行形式因为z[k]依赖本轮的z[k1]从后往前算就错了。这是我在移植到向量指令时踩过的坑编译器自动向量化打开后结果会漂加#pragma或者改回标量最稳。3. 系数与初值实现从导出 b/a 到 C 数组的三步主流程写对只是及格线结果能不能和参考实现对得上八成取决于系数和初值这两个输入。这两样东西都不是在 C 里现算的而是在设计环境里算好、抄进来、再在 C 里补初值。3.1 系数怎么来的参数表比代码更值得先确认设计阶段要定的就四个量任何一个改了都得重新导出参数含义常见取值与影响fs采样率Hz由采集硬件决定改不了fc截止频率Hz取关心的最高频率的 1.21.5 倍留出平方衰减的余量order滤波器阶数26 阶够用超过 8 阶直接型系数开始敏感rp/rs通带/阻带指标切比雪夫、椭圆型才需要巴特沃斯不用导出时确认三件事a[0]是不是 1不是就整体除掉b、a 长度是否相等不等就补 0 对齐系数是双精度还是单精度贴进 C 的时候别截断。我一般会把导出的 b、a 连同fs、fc一起写进注释半年后回来改参数没有这行注释基本等于重做。阶数选择上有个血泪经验为了压住某个窄带干扰很容易一路把阶数堆到 10 阶以上直接型在 double 下还能撑换成 float 跑长数据就开始飘。真要高阶用二阶节级联SOS每个节单独滤初值逐节传稳定性好得多代价是代码量翻几倍。3.2 初值的两条路硬编码和现场迭代初值的物理含义是「输入恒为x[0]时滤波器稳态下的内部状态」。参考实现的做法是解一个线性方程组C 里没必要照抄。路线一在设计环境里算好直接硬编码。数据长度固定、系数不常改的项目这是最省事的也最不容易错# 在设计环境里执行一次把输出抄进 coeff.h # b, a 设计函数(阶数, 截止频率/奈奎斯特) # zi 初值函数(b, a) # 长度 max(len(b),len(a))-1 # print(, .join(f{v:.17g} for v in zi)) # 17 位有效数字别省打印出来带 17 位有效数字是有讲究的双精度十进制往返要 17 位才无损抄成 6 位小数端点稳态就会有可见偏差这种偏差在小信号上尤其明显。路线二在 C 里现场迭代。原理是把单位阶跃喂进去反复更新状态直到收敛收敛后的状态就是稳态初值/* 以单位阶跃为激励做不动点迭代求稳态初值b、a 已归一化且等长 */ static int lfilter_zi(const double *b, const double *a, int n, double *zi) { double z[IIR_MAX_ORDER] {0}; for (int it 0; it 5000; it) { double y b[0] z[0]; /* 输入恒为 1 */ for (int k 0; k n - 2; k) z[k] b[k 1] - a[k 1] * y z[k 1]; z[n - 2] b[n - 1] - a[n - 1] * y; } memcpy(zi, z, sizeof(double) * (size_t)(n - 1)); return 0; }收敛速度取决于极点离单位圆有多近。截止频率很低、阶数又高的滤波器5000 次可能不够稳妥的做法是加一个增量判据两次迭代的z差值的无穷范数小于1e-14就提前退出超了上限还没收敛说明系数本身有问题不稳定这时候该报错而不是硬算。拿到归一化的zi之后实际初值是zi[k] * x[0]这一步不能漏否则不同幅值的信号会得到不同的端点行为。3.3 主流程延拓、正向、反转、反向、去延拓把上面几块拼起来就是完整流程。下面这份实现可以直接对着改返回码区分了「数据太短」和「内存不足」两种失败便于上层判断。/* filtfilt.c —— 零相位滤波主流程 */ #include stdlib.h #include string.h #include iir.h #include filtfilt.h int filtfilt(const double *b, const double *a, int n, const double *x, int N, double *y, int padlen) { if (n 2) return -3; /* 无状态环节不需要零相位处理 */ if (N padlen) return -1; /* 数据比延拓段还短判错 */ int M N 2 * padlen; double *ext (double *)malloc(sizeof(double) * (size_t)M); double *tmp (double *)malloc(sizeof(double) * (size_t)M); if (!ext || !tmp) { free(ext); free(tmp); return -2; } /* 奇对称延拓2*x[0]-x[k] 保证首尾一阶差分连续 */ for (int i 0; i padlen; i) { ext[i] 2.0 * x[0] - x[padlen - i]; ext[M - 1 - i] 2.0 * x[N - 1] - x[N - 2 - i]; } memcpy(ext padlen, x, sizeof(double) * (size_t)N); double zi[IIR_MAX_ORDER], z[IIR_MAX_ORDER]; lfilter_zi(b, a, n, zi); /* 第一遍正向初值按首样本缩放 */ for (int k 0; k n - 1; k) z[k] zi[k] * ext[0]; iir_df2t(b, a, n, ext, tmp, M, z); /* 时间反转 */ for (int i 0; i M / 2; i) { double t tmp[i]; tmp[i] tmp[M - 1 - i]; tmp[M - 1 - i] t; } /* 第二遍对反转序列再滤一次初值按反转后的首样本缩放 */ for (int k 0; k n - 1; k) z[k] zi[k] * tmp[0]; iir_df2t(b, a, n, tmp, ext, M, z); /* 再反转一次并切掉两端延拓段 */ for (int i 0; i N; i) y[i] ext[M - 1 - padlen - i]; free(ext); free(tmp); return 0; }几个参数要交代清楚。padlen是两端各延拓多少点建议取3 * (n - 1)第 4 章展开讲这个 3 是怎么来的。ext和tmp都是M长度内存峰值大约是输入的三倍长数据要留够堆空间。最后一步的下标M - 1 - padlen - i是两次反转叠加的结果写错成padlen i不会崩只会得到一段莫名其妙的结果这类下标建议在纸上画一遍再落笔。zi只在函数开头算一次两次滤波复用同一份归一化初值分别乘各自的起点值这是这套流程里最容易写错的地方——第二次滤波乘的是tmp[0]也就是原始数据的最后一个点不是x[0]。4. 端点瞬态的根因padlen、奇延拓与那个 3 倍阶数延拓长度和延拓方式是两个独立参数很多人把它们混成一件事结果调来调去总有一头不对。搞清各自管什么参数就不用试了。4.1 瞬态从哪来两个边界各管一头正向滤波的瞬态发生在起点反向滤波的瞬态发生在终点两次叠加之后原始数据的头和尾各有一段不可信区间。这个区间的长度和滤波器的时间常数成正比大致是极点衰减到可忽略所需的点数。延拓的作用是把这个区间挪到假数据上去。头部延拓padlen点正向滤波的瞬态就落在假的头段里尾部同理反向滤波的瞬态落在假的尾段里。切掉之后真实数据的首尾就都是「已经稳下来」的结果。初值的作用是让起点本来就不产生瞬态。两者叠加才是完整方案初值把残余瞬态压到很小延拓把这个很小的残余挪出真实数据区。只做初值不做延拓尾部还是会有残留只做延拓不做初值padlen得取很长才够白白浪费算力。4.2 零延拓、镜像延拓、奇延拓的差别三种常见做法我都试过差别在端点附近几十个点上肉眼不一定看得出来但做互相关或者算峰值时刻时影响很直接。零延拓最简单两端补 0。问题是如果信号本身有直流偏置补 0 相当于在边界造了一个大阶跃瞬态反而更猛端点会出现明显的上冲或下坠这是最容易翻车的一种。镜像延拓补x[padlen..1]的翻转连续性好一些但一阶导数在接缝处反号对有趋势的信号会在端点造出一个假的拐点。奇对称延拓补2*x[0] - x[k]本质是把信号关于端点做中心对称接缝处值和一阶差分都连续对线性趋势的信号几乎不引入偏差。成熟实现默认用这个代价是实现上多一次减法索引也更容易写错——2*x[0]-x[padlen-i]里的padlen-i不是i1这点对不上下标的话端点会整体错位一个点。4.3 padlen 取 3 倍阶数的由来瞬态的衰减包络大致按r^n走r是主导极点的模。对常见的巴特沃斯低通r和阶数的关系决定了3*(n-1)这个经验值足够把残余压到千分之一量级。低于 2 倍端点还能看出抖动高于 5 倍收益已经不明显只是多花内存。有个硬约束padlen必须小于数据长度否则奇延拓取不到x[padlen]直接越界。短数据比如长度 20 点、4 阶滤波器要延拓 9 点就得把padlen压到N/2以下此时端点精度必然打折这是物理限制不是代码问题。提示改完padlen一定要重跑端点对比只看中间段的波形是看不出问题的。截取前padlen个点和参考实现逐点做差误差应该在 1e-10 量级双精度、系数一致的前提下。5. 避坑排查五个把结果搞错却不报错的细节这套代码最麻烦的地方是所有错误都是静默的编译过、跑得动、波形也像那么回事但数值就是不对。下面五条是我和身边人实际踩过的按出现频率排。5.1 数值与系数类问题现象一首尾各几十个点明显上翘或下坠中间完全正常。原因基本是没做延拓或者padlen给小了瞬态没被挪出真实数据区。也可能是延拓写了但下标错位比如把x[padlen-i]写成x[i1]边界值看着连续实际差了一个点。解决先把padlen临时调到5*(n-1)看端点是否恢复正常恢复就说明是长度问题没恢复就打印延拓段的前 10 个值和手算的2*x[0]-x[k]逐个对。现象二整体幅值偏大或偏小形状对得上。大概率是a[0]没有归一化成 1。设计工具导出的系数有时带一个公共缩放因子直接用会让输出乘以这个因子。另一种情况是把二阶节SOS的系数矩阵当成单个 b/a 传进来了n算出来是 3 但语义完全不对。解决进函数第一件事断言fabs(a[0]-1.0) 1e-12SOS 就老老实实逐节调用别想着拼成高阶。现象三float 版本跑长数据后段开始发散或者出现规则振荡。高阶直接型在单精度下系数敏感极点稍微偏出单位圆就完蛋。双精度跑同样的系数没事换 float 就出事说明是精度问题不是逻辑问题。解决内部状态和中间变量强制用 double只让输入输出用 float或者降到 4 阶以下、改用级联结构。这类问题最迷惑人因为它和数据长度相关短样本测不出来。5.2 流程与使用方式问题现象四和参考实现逐点对比中间段能对上端点差 1e-3 量级。不是延拓的问题是初值没生效。常见原因有两个zi算完忘了乘起点值或者第二次滤波沿用了第一次滤波结束时的状态而不是重新用zi * tmp[0]。还有一种是zi硬编码时截断到 6 位小数这个精度对双精度来说不够。解决把zi和参考实现打印出来的值逐位比小数点后对不上就重新导出别用 printf 的默认精度。现象五把这套代码搬到实时流处理里输出整体错位。filtfilt 是非因果的第 n 点的输出依赖第 n 点之后的数据实时场景根本用不了。有人会把它套在滑动窗口上窗口内做 filtfilt、只取中间一段输出这叫重叠保留能用但窗口长度必须大于2*padlen 输出段长度否则边界还是脏的。解决实时场景老老实实用单向滤波加固定的相位补偿IIR 群延迟不是常数只能近似补偿或者改用线性相位 FIR 加时移。零相位这件事在实时里没有免费午餐。注意这五条里有四条不会触发任何断言验证手段只有一条——和参考实现逐点对比残差。别靠肉眼看波形波形的宽容度远超你的想象。6. 验证到位再交付扫频对相位、脉冲看对称、分段处理长数据代码能跑通不等于结果可信。我现在交付任何滤波实现之前固定跑三组验证加起来不到二十行脚本能挡掉九成以上的隐性错误。第一组脉冲响应看对称性。构造一个长度 512 的单位脉冲第 0 点为 1其余为 0走一遍 filtfilt输出应该关于中点对称且左右两侧的差在 1e-12 以内。不对称就说明两次滤波的某一步方向搞反了或者反转下标写错。这组测试最灵敏下标错一位它立刻现形。第二组线性扫频看峰值对齐。生成一段从 1Hz 扫到 40Hz、时长 10 秒的 chirp采样率 200Hz过一遍滤波然后逐段找输入输出的过零点位置。零相位滤波的正确表现是过零点时刻基本重合偏差应当小于一个采样周期。如果输出系统性地滞后若干个点说明只滤了一遍如果滞后量随频率变化说明相位没消干净回去查第二次滤波是不是真的作用在反转序列上。这组测试比看幅频曲线直观得多用眼睛就能判。第三组残差比对定量。同一组系数、同一段数据和参考实现的结果逐点做差取最大绝对误差和均方根误差。双精度、系数一致的前提下端点误差应该在 1e-10 量级中间段更小。误差在 1e-3 量级通常是初值问题1e-6 量级多半是系数截断。把这两个数写进单元测试以后改padlen或者换系数跑一遍就知道有没有退化。长数据还有个工程技巧值得单独说。一次性 filtfilt 需要三倍内存几百兆的采集文件直接爆。稳妥做法是分段处理每段长度L前后各多取2*padlen点上下文段内做 filtfilt只保留中间L点段与段之间按L步进重叠区直接丢弃。分段长度建议不小于50*padlen太短了每段的瞬态占比过高累积误差会显出来。我一般会先用整段跑一次做基准再用分段跑一次比较两者的最大差值确认在 1e-9 以内才敢用到正式流程里。最后留一个习惯任何一次改完系数我都会把fs、fc、order、padlen四个值和残差指标一起打印到日志里。有次排查一个三个月前的异常数据翻日志才发现当时用的是另一个截止频率参数写在代码里但没人记得改过。从那以后我每次交付滤波模块都强制走一遍「脉冲对称 扫频对齐 残差定量」这三步参数和指标一并落日志省下的排查时间远比写测试的时间多。希望帮到你。本文还有配套的精品资源点击获取