企业官网建设流程全解析

简介：面向数字信号处理学习者与嵌入式开发者的双向滤波 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 实现的文件构成

不同人打包的命名习惯不一样，但函数职责跑不出下面这几类。看代码结构比看文件名可靠，打开头文件扫一眼函数签名，基本就能定位主流程在哪。

这里最容易出问题的是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从状态全零起步，等价于假设输入在 t<0 一直是 0。如果x[0]本身不是 0，这个跳变就是一次阶跃激励，输出的前若干点会叠上一段按极点衰减的瞬态。极点越靠近单位圆，衰减越慢，有时候几百个点都散不掉。

解决办法有两条，成熟实现是两条一起用。一条是延拓，在数据两端造出一段「比真实数据更早/更晚」的假数据，让瞬态发生在假数据上，滤完再切掉；另一条是初值，直接算出「如果这段信号从无穷远就一直存在，此刻的状态该是多少」，让滤波器从稳态起步。前者解决数据边界，后者解决状态边界，缺一个端点都会有残差。

2.4 单次滤波的骨架：直接 II 型转置

不管上层怎么组织，底层那一次滤波大概率是直接 II 型转置（DF2T），状态长度等于阶数，数值上比直接 I 型稳，中间变量也少。

/* iir.c —— 单次 IIR 滤波，直接 II 型转置结构 */ #include "iir.h" /* b、a 需已归一化(a[0]=1)并补齐到同一长度 n；z 长度 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))，短的那个尾部补 0；z的长度是n-1，调用方负责分配和初始化，函数内部不碰它的初值；z是原地改写的，所以同一份状态不能同时喂给两个线程。那个n<2的分支不是凑数，一阶增益环节没有状态，漏掉它会读到长度为 0 的数组。

内层循环不能写成z[k] = ...的并行形式，因为z[k]依赖本轮的z[k+1]，从后往前算就错了。这是我在移植到向量指令时踩过的坑，编译器自动向量化打开后结果会漂，加#pragma或者改回标量最稳。

3. 系数与初值实现：从导出 b/a 到 C 数组的三步

主流程写对只是及格线，结果能不能和参考实现对得上，八成取决于系数和初值这两个输入。这两样东西都不是在 C 里现算的，而是在设计环境里算好、抄进来、再在 C 里补初值。

3.1 系数怎么来的：参数表比代码更值得先确认

设计阶段要定的就四个量，任何一个改了都得重新导出：

导出时确认三件事：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不是i+1，这点对不上下标的话端点会整体错位一个点。

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[i+1]，边界值看着连续，实际差了一个点。解决：先把padlen临时调到5*(n-1)看端点是否恢复正常，恢复就说明是长度问题，没恢复就打印延拓段的前 10 个值，和手算的2*x[0]-x[k]逐个对。

现象二：整体幅值偏大或偏小，形状对得上。大概率是a[0]没有归一化成 1。设计工具导出的系数有时带一个公共缩放因子，直接用会让输出乘以这个因子。另一种情况是把二阶节（SOS）的系数矩阵当成单个 b/a 传进来了，n算出来是 3 但语义完全不对。解决：进函数第一件事断言fabs(a[0]-1.0) < 1e-12，SOS 就老老实实逐节调用，别想着拼成高阶。

现象三：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四个值和残差指标一起打印到日志里。有次排查一个三个月前的异常数据，翻日志才发现当时用的是另一个截止频率，参数写在代码里但没人记得改过。从那以后我每次交付滤波模块，都强制走一遍「脉冲对称 + 扫频对齐 + 残差定量」这三步，参数和指标一并落日志，省下的排查时间远比写测试的时间多。希望帮到你。

本文还有配套的精品资源，点击获取