★ FEATURED ARTICLE 简介filtfilt.rar面向数字信号处理初学者以及需要把MATLAB零相位滤波迁移到C环境的开发者核心是filtfilt双向滤波算法的C语言实现参考。压缩包共1个文件为单个c源码整体约2KB体量轻便便于直接阅读滤波系数定义、输入输出组织与前向—反向滤波流程。已有382人学习说明其在滤波实现入门场景中具备一定参考热度。代码围绕先正向滤波、再反向滤波以抵消相位失真的思路展开涉及滤波器类型与截止频率设定、IIR/FIR系数应用、反向滤波时的系数顺序处理、结果累加以及零填充或循环缓冲等边界条件处理可帮助读者理解零相位响应为何适合音频、图像与生物医学信号分析。读者可据此梳理从滤波器设计到C实现的完整链路并结合MATLAB原生filtfilt做结果对照排查相位偏移与边界不连续问题适合作为课程实验、项目移植和算法验证的轻量参考资料。1. 一个整体后移的波形把 filtfilt 推到台前某实验室做振动信号的到时差测量两路传感器采同一个敲击事件各自过一遍 IIR 低通再比对峰值位置算时间差。结果每次滤完波形都整体往右挪了几十个采样点而且两路挪的距离还不一样——因为挪多少取决于这段信号里各频率成分的占比信号变了延迟就变。测出来的时间差自然是错的排查了两天才把矛头指向滤波器本身的相位响应。换用 filtfilt 的零相位做法后峰值位置纹丝不动。这份filtfilt.rar就是把这套逻辑用 C 语言实现出来的源码包不依赖任何数值计算库拿过去就能编进采集固件或上位机程序里。适合手上有 IIR 系数、又不想让波形发生时间位移的开发者尤其是做生物电、振动、音频对齐这类对相位敏感场景的人。2. filtfilt 在 C 里该怎么落地零相位语义与函数接口设计2.1 前向-反向两次滤波为什么能抹掉相位普通 IIR 滤波是一个因果系统频率响应写成 H(ω)它的相位是频率的函数所以不同频率成分被延迟的量不同宽带信号滤完形状和时间位置都会变。filtfilt 的做法很朴素先用 H(ω) 正向滤一遍把输出序列在时间轴上翻转再用同一个 H(ω) 滤第二遍最后翻回来。第二次滤波作用在翻转后的序列上对上一步的输出相当于乘上 H(ω) 的共轭 H*(ω)。两次串起来总响应是 H(ω)·H*(ω) |H(ω)|²。这是一个实数相位恒为零群延迟为 0所有频率成分的延迟量一致——等于说没有任何一个成分被相对移动。代价也得说清楚这是很多人上手第一个翻车点幅度响应变成了原来的平方。原本 -3 dB 的截止点平方之后变成 -6 dB也就是通带边缘被压得更狠实际有效的截止频率会往低频方向挪同时阻带衰减翻倍40 dB 的滤波器串两遍就是 80 dB。设计系数的时候要是按单次滤波的指标来滤完会发现带宽比预期窄了一截。我在工程上的习惯是先用单次滤波把通带边缘放到略高于目标值的位置再跑 filtfilt最后用扫频确认一遍实际的 -3 dB 点。还有一点绕不开——整个过程需要整段数据。第二次滤波要从序列末尾往回走意味着它必须看到未来的样本所以 filtfilt 天生是非因果的、离线的。实时管线里想用只能做分块加重叠这个后面第 6 章再展开。2.2 系数、边界与初值参考实现的语义拆成三块可移植逻辑把 MATLAB、SciPy 里同名函数的行为拆开C 实现其实只需要复刻三件事剩下都是体力活。第一件是系数归一化。参考实现接收的分子系数b和分母系数a约定a[0]是首项但不要求它等于 1。内部会先把两组系数都除以a[0]再做滤波。C 代码里如果跳过这一步当a[0]是 0.5 或者 2 的时候输出会整体放大或缩小而且因为状态递推也跟着错波形不只是缩放还会出现奇怪的漂移。这是我见过的最高频的低级错误。第二件是滤波器阶数的定义。nfilt max(nb, na)状态数order nfilt - 1。这个order决定了两件事一是状态数组长度二是边界延拓的长度。参考实现用的延拓长度是3 * order前后各补这么多点。第三件是初始条件。如果直接拿全零状态开始滤第一段样本会经历一个从零状态爬到稳态的瞬态过程这段瞬态在延拓区里衰减掉一部分但残余会渗进有效数据的前若干个点。参考实现的做法是算出稳态初始状态zi让滤波器一上电就在稳态上瞬态基本消失。这两条路我都写过零状态加长延拓在多数场景够用但对边界特别敏感的应用比如截断的短信号做积分还是老实算zi。延拓方式本身也有讲究。参考实现用的是奇对称延拓左端用2*x[0] - x[pad-i]镜像翻折右端同理。这样延拓出来的序列在边界处一阶连续不会引入一个阶跃避免滤波后在首尾憋出一个大尖峰。如果偷懒用零填充或者直接复制端点值首尾会出现明显的塌陷或鼓包做频谱分析时尤其明显。2.3 资源包的文件构成与依赖取向从包名和关键词看这是一个纯 C 的滤波实现集合重点在filt和filtfilt两个接口。这类包常见的组织方式是这样下列是典型构成实际以解压后为准文件作用filtfilt.h对外接口声明、错误码、结构体定义filtfilt.c单步滤波、稳态初值计算、filtfilt 主流程demo.c生成测试信号并调用滤波的示例Makefile/ 编译说明编译目标与依赖说明依赖取向很关键这类包一般只依赖 C 标准库和libm不引任何第三方数值库。libm的用途主要在两个地方——计算稳态初值时如果涉及除法归一化还有 demo 里生成正弦测试信号用的sin()。真正跑滤波的循环里没有任何超越函数调用全是乘加这一点对嵌入式移植很友好。接口设计上常见做法是把系数、输入、输出、长度作为参数直接传返回值用int表示错误码负数代表失败。缓冲区的分配策略分两派一派在函数内部malloc调用方省事但有堆碎片风险另一派要求调用方传入一块长度为nx 6*order的工作区零动态分配。如果这套代码是要进 MCU 的选后一种只是上位机跑批处理前一种写起来更省心。提示解压后先看头文件里的函数签名和错误码定义再决定接哪一派。别急着改实现接口一旦定了后面的示例和测试都是围绕它写的。3. 核心代码拆解一次 filtfilt 调用从延拓到裁剪的完整链路3.1 用直接 II 转置型承载滤波器状态单步滤波我选直接 II 转置型transposed direct form II。理由有三状态变量个数最少等于阶数数值动态范围比直接 I 型好定点化时溢出风险低。二阶节级联的写法在高阶滤波器上更稳但对这个包来说直接把通用形式写出来更通用。/* 滤波器实例系数指针 状态缓存 */ typedef struct { int order; /* 状态个数 max(nb, na) - 1 */ const double *b; /* 已按 a[0] 归一化的分子系数 */ const double *a; /* 已按 a[0] 归一化的分母系数 */ int nb, na; /* 两组系数的有效长度 */ double *z; /* 状态数组长度 order */ } iir_t; /* 单步推进输入 x返回 y内部更新状态 */ static double iir_step(iir_t *f, double x) { double y; int k, n f-order; if (n 0) /* 退化成一阶增益没有状态 */ return f-b[0] * x; y f-b[0] * x f-z[0]; /* z[k] 对应第 k1 个状态逐个向后移位 */ for (k 0; k n - 1; k) { double bx (k 1 f-nb) ? f-b[k 1] * x : 0.0; double ay (k 1 f-na) ? f-a[k 1] * y : 0.0; f-z[k] f-z[k 1] bx - ay; } /* 最后一个状态没有 z[k1] 可借单独收尾 */ f-z[n - 1] ((n f-nb) ? f-b[n] * x : 0.0) - ((n f-na) ? f-a[n] * y : 0.0); return y; }这段循环是整份代码里被调用次数最多的部分参数含义要抠清楚。f-b[0]是分子首项f-z[0]是第一个状态两者相加就是当前输出。后面每个状态z[k]在下一次迭代里要变成z[k1]加上新的输入贡献、减去新的输出反馈这就是转置结构里状态跟着数据流走的写法。边界判断k 1 f-nb和k 1 f-na不能省。分子和分母长度往往不相等比如一个 4 阶巴特沃斯低通nb na 5但如果系数经过某种化简导致尾部是零被截掉长度就不一致了。少了这两个判断读到数组外面去结果就是随机数——这种 bug 在 debug 模式下不加边界检查根本看不出来。3.2 稳态初始条件怎么算出来上一步说了零状态起步会有瞬态。稳态初始状态的思路是假设滤波器从无穷远以前就一直吃着一个恒定输入x0那么输出必然收敛到一个恒定值y0此时各个状态也固定不变。把转置结构在稳态下的递推式反解出来就能得到一组z。/* 计算稳态初始状态。 * x0: 起始处的输入值y0: 对应稳态输出 dc_gain * x0 * z: 输出长度 order */ static void calc_steady_zi(const double *b, int nb, const double *a, int na, double x0, double y0, double *z, int order) { int k; /* 最后一个状态只依赖当前输入输出 */ z[order - 1] ((order nb) ? b[order] * x0 : 0.0) - ((order na) ? a[order] * y0 : 0.0); /* 从后往前逐个反推 */ for (k order - 2; k 0; --k) { double bx ((k 1) nb) ? b[k 1] * x0 : 0.0; double ay ((k 1) na) ? a[k 1] * y0 : 0.0; z[k] z[k 1] bx - ay; } }这里的推导是第 3.1 节那组递推式的逆运算。稳态下z[k]不再随时间变化把z[k] z[k1] b[k1]x0 - a[k1]y0这个关系从末尾往前套一遍就能把整组状态解出来。y0的取法是直流增益乘以x0直流增益等于sum(b) / sum(a)——系数已经按a[0]归一化过所以直接用求和比就行。参数里最容易搞错的是order和nb、na的关系。order是状态数比系数长度小 1当order nb不成立时说明分子比分母短那几项就得按零处理。这类下标错一位的问题在小阶数滤波器上可能完全看不出异常换成 8 阶椭圆滤波器立刻发散。算完左端的zi正向滤波从延拓区的第一个点开始状态直接灌进去。反向那一遍同理用右端点的输入值算一组新的zi从序列末尾往前走。两遍用的zi必须各自独立算共用一组是错的——我就干过这事滤完发现波形右半边有点歪查了半天才发现是反向时沿用了正向的状态。3.3 奇延拓与前后向滤波主循环主函数把前面几块拼起来。顺序是归一化系数 → 分配缓冲 → 奇延拓 → 正向滤 → 反向滤 → 裁掉延拓区。int filtfilt(const double *b_in, int nb, const double *a_in, int na, const double *x, int nx, double *y) { int order, pad, n, i, ret 0; double *bn NULL, *an NULL, *buf NULL, *z NULL; double a0, gain, x0, x1; iir_t f; if (!b_in || !a_in || !x || !y || nb 1 || na 1 || nx 1) return -1; if (a_in[0] 0.0) /* a[0] 不能为零 */ return -3; order (nb na ? nb : na) - 1; pad 3 * order; if (nx pad) return -2; /* 数据太短延拓会吃掉全部样本 */ /* 1. 按 a[0] 归一化系数后面所有计算都用归一化版本 */ a0 a_in[0]; bn (double *)malloc(sizeof(double) * nb); an (double *)malloc(sizeof(double) * na); for (i 0; i nb; i) bn[i] b_in[i] / a0; for (i 0; i na; i) an[i] a_in[i] / a0; n nx 2 * pad; buf (double *)malloc(sizeof(double) * n); z (double *)calloc(order 0 ? order : 1, sizeof(double)); if (!bn || !an || !buf || !z) { ret -4; goto done; } /* 2. 奇对称延拓左端 2*x[0] - x[pad..1]右端同理 */ for (i 0; i pad; i) buf[i] 2.0 * x[0] - x[pad - i]; memcpy(buf pad, x, sizeof(double) * nx); for (i 0; i pad; i) buf[pad nx i] 2.0 * x[nx - 1] - x[nx - 2 - i]; /* 3. 正向滤波起始状态取稳态值 */ f.order order; f.b bn; f.a an; f.nb nb; f.na na; f.z z; gain 0.0; for (i 0; i nb; i) gain bn[i]; { double den 0.0; for (i 0; i na; i) den an[i]; gain (den ! 0.0) ? gain / den : 0.0; } x0 buf[0]; calc_steady_zi(bn, nb, an, na, x0, gain * x0, z, order); for (i 0; i n; i) buf[i] iir_step(f, buf[i]); /* 4. 反向滤波序列翻转后重跑一遍状态用右端的稳态值 */ x1 buf[n - 1]; calc_steady_zi(bn, nb, an, na, x1, gain * x1, z, order); for (i n - 1; i 0; --i) buf[i] iir_step(f, buf[i]); /* 5. 裁掉两端延拓区写回调用方缓冲 */ memcpy(y, buf pad, sizeof(double) * nx); done: free(bn); free(an); free(buf); free(z); return ret; }几个点值得单独说。系数归一化放在最前面bn和an是堆上拷贝不改调用方传进来的原始数组这个习惯要守住否则调用方复用系数时会发现值被改过。pad取3 * order是跟参考实现对齐的阶数越高需要的延拓越长不足以让瞬态完全衰减。第 3 步里算gain的写法有点绕其实两行就能算完拆开来是为了把分子求和与分母求和分开方便调试时打印中间值。gain的含义是直流增益用在calc_steady_zi里推y0。第 4 步的反向循环直接从n-1开始倒着走iir_step本身不用改——滤波器是时不变的正向反向对它没有区别翻转动作完全体现在遍历方向上。这一步是我觉得整个算法最漂亮的地方一个函数两种用法。第 5 步裁剪用memcpy而不是for循环编译器能识别成一条memmove。裁剪的位置是buf pad长度nx跟延拓时memcpy的偏移一致。这两个偏移量一旦对不上输出就会整体错位几个点症状是滤完波形平移了——跟一开始要解决的问题一模一样很容易误判成算法没生效。3.4 一个能编译能跑的调用示例给一个最小可运行例子。设计一个 4 阶巴特沃斯低通截止频率归一化到 0.1对应采样率的五分之一系数用双线性变换手算太麻烦实际工程里一般由设计工具生成后粘贴进来这里直接写死一组。#include stdio.h #include math.h #include filtfilt.h int main(void) { /* 4 阶低通截止 0.1*fs系数已归一化到 a[0]1 */ static const double b[5] { 0.000416, 0.001665, 0.002497, 0.001665, 0.000416 }; static const double a[5] { 1.0, -3.180638, 3.861195, -2.112001, 0.438273 }; enum { N 512 }; double x[N], y[N]; int i; /* 构造双频信号低频分量 一个阶跃便于观察相位是否位移 */ for (i 0; i N; i) { double t (double)i; x[i] sin(2.0 * M_PI * 0.02 * t) 0.3 * sin(2.0 * M_PI * 0.30 * t); } x[N / 2] 1.0; int rc filtfilt(b, 5, a, 5, x, N, y); if (rc ! 0) { fprintf(stderr, filtfilt failed, rc%d\n, rc); return 1; } /* 打印峰值位置对比零相位滤波下二者应完全重合 */ int px 0, py 0; for (i 1; i N; i) { if (x[i] x[px]) px i; if (y[i] y[py]) py i; } printf(peak x%d peak y%d shift%d\n, px, py, py - px); return 0; }N选 512远大于pad 3*4 12不会触发长度不足的错误。信号里那个x[N/2] 1.0是人为注入的尖峰用来看滤波前后的位置关系两个正弦分量一个在通带内、一个在阻带附近能同时观察通带保形和阻带衰减。filtfilt的五个参数依次是分子系数、分子长度、分母长度、分母长度、输入、输入长度、输出。返回值rc必须检查-2是数据太短-3是a[0]为零这两个是最常见的调用错误。跑出来shift应该是 0。如果打印出来是负数或者正数几个点先别怀疑算法检查裁剪偏移是不是写成了buf nx之类如果shift恰好等于pad那就是裁剪完全没做输出的是延拓后的整段。这两种症状我看过太多次都是复制粘贴时改漏一个变量。4. 编译与验证把 rar 解开后怎么确认它真的零相位4.1 编译命令与目录约定解压之后先进源码目录看一眼有没有构建脚本。有Makefile就直接make没有的话手动编也不复杂整个包不依赖第三方库。# 假设源码在 src/ 下头文件也在同目录 gcc -O2 -stdc99 -Wall -Wextra -pedantic \ -I./src \ ./src/filtfilt.c ./src/demo.c \ -o demo -lm-stdc99是因为代码里用了//注释和混合声明别指望它能在 C89 下编过。-lm必须加在源文件后面链接器对库的顺序敏感写在前面会报未定义符号——这个坑每年都有人踩。-Wall -Wextra打开让编译器把隐式截断、符号比较这些隐患吐出来-pedantic是习惯跨平台交付前跑一遍能提前挡掉编译器扩展语法。编译过了先别急着跑自己的数据用包内自带的示例跑一遍确认输出是个合理的波形而不是一堆nan。如果 demo 一跑就段错误八成是缓冲区长度算错了回头看 3.3 节里n nx 2*pad这个式子。4.2 用单位冲激和单频正弦验证零相位验证零相位不能靠肉眼看波形好像没动得有可量化的判据。我一般做两组测试。第一组是单位冲激。构造一个全零序列中间放一个 1滤完看输出。对称滤波器的冲激响应应该以冲激位置为中心左右对称如果往一边偏说明反向那一遍没生效或者状态用错了。这组测试最快一眼就能看出问题。第二组是单频正弦测相位偏移。取一个在通带内的频率生成一段足够长的正弦滤完跟原信号算互相关峰值滞后应该是零。/* 用互相关峰值位置判断相位偏移理想结果为 0 */ static int measure_shift(const double *x, const double *y, int n) { int lag, best_lag 0; double best -1e300; /* 只搜 -32 个点足够覆盖任何合理的相位偏差 */ for (lag -32; lag 32; lag) { double sum 0.0; int i; for (i 64; i n - 64; i) { int j i lag; if (j 0 || j n) continue; sum x[i] * y[j]; } if (sum best) { best sum; best_lag lag; } } return best_lag; }lag的搜索范围取 ±32 是经验值。如果滤波器阶数很高、或者截止频率压得很低真实延迟可能超出这个范围那就把窗口放大到 ±order 的两倍。互相关开头结尾各留 64 个点是为了避开边界效应——边界处理得再好首尾几十个点的能量跟中间还是有差异把它们算进去会污染峰值位置。返回 0 说明相位偏移为零。返回非零先别慌检查测试信号是不是用了通带边缘的频率那里幅度衰减大信噪比低互相关峰值可能被噪声带偏。换到通带中心再测一次如果变成 0那滤波器本身是好的。4.3 和参考实现对齐时的容差与坑最终验收我一般会拿同一组数据、同一组系数跟参考实现MATLAB 和 SciPy 里的同名函数各跑一遍逐点比对。容差设多少合适double精度的实现相对误差在1e-10量级是正常的因为两边浮点运算顺序不同累加次序带来的舍入差异会放大到几十个 ulp。我一般设1e-9的相对容差超过这个量级说明实现里某处逻辑有差异不是舍入能解释的。对不上时按这个顺序排查效率最高先看前 100 个点是不是差得最多——如果是问题在初始条件zi没算或者算错了再看中间段是不是稳定差一个常数倍——那是a[0]归一化漏了最后看末尾若干个点——反向滤波的初始状态用错了。这三个位置基本能定位九成的不一致。还有个容易被忽略的参考实现的延拓长度和你的可能不同。3*order是常见取值但不同版本之间会调整。如果只在边界附近有微小差异、中间完全一致那就不是 bug是延拓策略不同功能上等价可以放过。5. 避坑与排查五条踩出来的记录5.1 边界和长度类的问题现象一滤完波形首尾往里塌末端几十个点直接逼近零。原因是没有做延拓或者延拓长度不足。滤波器从全零状态起步前几个点的输出被拉低到零附近这个瞬态持续的时间跟滤波器的时间常数成正比阶数越高、截止频率越低拖得越长。零填充相当于在信号两端人为制造了一个阶跃滤波后自然鼓包或塌陷。解决办法是把延拓长度提到3*order并且用奇对称延拓而不是补零。如果对边界要求极高比如做积分或者后续要算能量那就把稳态初始条件也启用两者叠加基本可以把边界误差压到可忽略。现象二函数返回 -2或者程序直接挂掉。nx 3*order时延拓区比数据本身还长越界访问是必然的。这个检查必须放在函数入口别等到分配完缓冲区才发现。工程上的处理是数据短于阈值时降级为零相位不可用直接返回错误让上层决定而不是悄悄缩小平移量——那样行为不可预测比报错更麻烦。另外注意nx为 0 的边界memcpy传 0 长度虽然安全但后面的x[nx-1]会读到x[-1]。入口处统一判定nx 1返回。5.2 系数和幅度类的问题现象三滤完幅度整体变大或变小波形形状对了但量纲不对。十有八九是a[0]没归一化。参考实现内部会除以a[0]手写实现时忘了这一步输出就是原结果的a[0]倍。检查方法很简单看输出峰值跟输入峰值的比值是不是恰好等于a[0]。修法是在函数开头统一做一次归一化拷到堆上全程用归一化后的系数。现象四滤波效果看着不错但通带比预期的窄截止频率偏低。这不是 bug是零相位滤波的固有特性——幅度响应被平方了。原本 -3 dB 的点变成 -6 dB视觉上截止频率就提前了。要么在设计系数时把指标放宽要么接受这个结果并重新标定。做定量分析前务必用扫频或已知频率的正弦测一遍实际的 -3 dB 点别照搬单次滤波的设计值。顺带一提通带内如果有纹波纹波的分贝数也会翻倍。1 dB 的纹波跑完 filtfilt 变成 2 dB对某些等波纹设计来说这个变化足够让指标不达标。5.3 数值和移植类的问题现象五在 PC 上跑得好好的移植到目标板后结果发散或者变成nan。先查数据类型。PC 上double是 64 位某些嵌入式平台配置下double会被降级成 32 位精度一掉高阶滤波器的极点靠近单位圆状态递推就会累积误差直到发散。这个问题很玄学同一组系数在 x86 上平稳换到目标板几百个点后开始抖。验证方式是打印中间状态的最大绝对值跟 PC 上的逐点对比。修法要么强制double精度加编译选项确认平台支持要么老老实实做二阶节级联——把高阶滤波器拆成若干个二阶节每节单独算状态数值稳定性比直接型高一个量级。真要是定点平台还得做系数缩放和饱和处理这块工作量不小得单独排期。注意做定点化之前先用浮点版本在同一平台上跑通并保存参考输出定点版本每一步都跟它比。没有参考输出的定点化就是在盲改。6. 二次开发把离线 filtfilt 改成能进实时管线的模块离线算法直接塞进实时管线是行不通的——它要看到整段数据才能开始反向那一遍。真正落地得做流式改造思路是重叠保留。把输入切成固定长度L的块每次处理时往前多取2*pad个历史样本拼成L 2*pad长度的窗口跑一遍 filtfilt然后只取中间L个点输出。这样每块的两端都有足够的延拓数据块与块之间的接缝不会出现跳变。pad取3*order是下限实践中我一般取4*order留一点余量。/* 流式包装每次喂 L 个点返回 L 个点 */ typedef struct { double *history; /* 长度 2*pad保存上一块尾部 */ int pad; int filled; /* 已积累的历史样本数 */ } stream_t; int stream_process(stream_t *s, const double *in, int L, double *out, const coeff_t *c) { int need 2 * s-pad; int total need L; double *win (double *)malloc(sizeof(double) * total); double *tmp (double *)malloc(sizeof(double) * total); if (!win || !tmp) { free(win); free(tmp); return -1; } memcpy(win, s-history, sizeof(double) * need); memcpy(win need, in, sizeof(double) * L); if (filtfilt(c-b, c-nb, c-a, c-na, win, total, tmp) ! 0) { free(win); free(tmp); return -2; } memcpy(out, tmp need, sizeof(double) * L); /* 滚动更新历史只保留最后 need 个样本 */ memmove(s-history, win L, sizeof(double) * need); free(win); free(tmp); return 0; }s-history的初始状态需要用第一块数据的左端点做一次奇延拓填充否则最开始那几块的输出还是会有瞬态。filled字段是用来标记历史是否已经攒够的没攒够的时候可以走一个简化分支用更短的窗口处理。这个结构有个必须接受的代价延迟。每块的输出要等到它后面的pad个样本都到齐了才能算出来所以实时性上会引入2*pad个样本的固有延迟。pad越大边界越干净延迟也越大。做闭环控制的时候这个延迟得算进环路补偿里不然会震荡。第二个改造方向是精度和性能。double改float能省一半内存带宽在 ARM Cortex-M7 这类带单精度 FPU 的核上速度提升明显但要注意极点靠近单位圆时精度损失。我的做法是保留double版本作为黄金参考float版本每次改完都跑一遍同样的测试信号逐点比相对误差超过1e-5就回头查是哪一节的问题。第三个方向是接口层面。如果这个包是要给别人用的把系数设计也一并封进去会好很多——传截止频率和阶数进去内部生成系数再滤波调用方不用关心双线性变换的细节。不过这会把包和具体的滤波器设计方法绑死看你的场景取舍。最后说个习惯。这套代码我前后改过三版第一版就是不检查a[0]、不验证边界交付后在现场数据上翻车波形首尾塌了一大截客户拿着图问是不是传感器坏了。从那以后我每次动滤波器代码都强制走一遍冲激响应看对称性、单频正弦看相位、和参考实现逐点比对这三步一步都不能少。希望帮到你。本文还有配套的精品资源点击获取
阅读完成 · 觉得有帮助？
