能否使用C语言实现scipy.lfilter？实现方法与原理

实现可行性结论 #

完全可以在C语言中复现和scipy.lfilter输出完全一致的滤波效果，不存在技术障碍。scipy.lfilter本身不是纯Python实现的函数，其核心递推逻辑本身就是用编译型语言（C/Fortran）编写的，没有依赖任何Python专属的黑盒逻辑，只要严格对齐参数预处理规则、滤波结构、递推顺序、状态更新逻辑，就能做到双精度下逐点误差在1e-15级别，和Python原生调用结果无感知差异。

scipy.lfilter核心实现原理 #

scipy.lfilter是通用的一维线性时不变滤波器实现，同时支持FIR、IIR滤波，底层采用直接II型转置结构完成递推计算，相比直接I型结构内存占用更低、缓存访问效率更高，核心逻辑如下：

1. 系数归一化：拿到输入的分子系数向量b、分母系数向量a后，首先取a[0]作为归一化因子，将所有b、a元素除以a[0]，保证归一化后a[0] = 1。这一步是内部默认执行的，很多自行实现的版本因为漏掉这一步导致结果偏差。
2. 差分方程递推：归一化后，输入序列x到输出序列y的计算遵循以下规则：

   维护长度为max(len(b), len(a)) - 1的延迟线状态数组z，存储每一步递推的中间状态。
   对每个采样点按顺序计算：

   1. 当前输出点 y[n] = b[0] * x[n] + z[0]
   2. 依次更新前state_len-1个状态：z[i] = b[i+1] * x[n] - a[i+1] * y[n] + z[i+1]
   3. 更新最后一个状态：根据b、a长度的差异，补全剩余的b乘x项、a乘y项，无对应系数的项按0计算。

3. 初始状态处理：如果用户没有传入初始状态zi，默认状态数组全0，对应零初始条件（即输入序列开始前所有输入、输出值均为0）；如果传入zi则直接用传入值初始化状态数组，支持分段滤波的连续拼接。

C语言复现步骤及参考代码 #

按照以下步骤实现即可完全对齐scipy.lfilter的行为：

1. 系数预处理：对传入的a、b数组做归一化，注意不要修改用户传入的原始系数数组，避免影响外部逻辑。
2. 状态初始化：根据系数长度计算状态数组长度，无传入初始状态时全置0。
3. 逐点递推：严格按照直接II型转置的顺序计算输出、更新状态，不要调换计算顺序，否则会引入状态更新的时序错误。

参考可直接运行的C实现如下：

#include <string.h>
#include <stdlib.h>

/**
 * @brief 与scipy.signal.lfilter行为完全一致的C实现，双精度浮点
 * 
 * @param b 滤波器分子系数数组
 * @param a 滤波器分母系数数组
 * @param nb b数组长度
 * @param na a数组长度
 * @param x 输入信号数组
 * @param y 输出信号数组，需提前分配与输入等长的内存
 * @param sig_len 输入信号长度
 * @param zi 初始状态数组，长度为max(nb, na)-1，传NULL则默认零初始状态
 * @param zf 滤波结束后的最终状态输出，传NULL则不返回，需提前分配max(nb, na)-1长度的内存
 */
void scipy_lfilter(const double *b, const double *a, int nb, int na,
                   const double *x, double *y, int sig_len,
                   const double *zi, double *zf) {
    int filt_order = nb > na ? nb : na;
    int state_len = filt_order - 1;
    // 申请临时存储归一化系数、状态的内存
    double *b_norm = (double*)calloc(filt_order, sizeof(double));
    double *a_norm = (double*)calloc(filt_order, sizeof(double));
    double *state = (double*)calloc(state_len, sizeof(double));

    // 系数归一化
    double a0 = a[0];
    for (int i = 0; i < nb; i++) b_norm[i] = b[i] / a0;
    for (int i = 0; i < na; i++) a_norm[i] = a[i] / a0;

    // 初始化状态
    if (zi != NULL) {
        memcpy(state, zi, state_len * sizeof(double));
    }

    // 逐点递推计算
    for (int n = 0; n < sig_len; n++) {
        y[n] = b_norm[0] * x[n] + state[0];
        // 更新前state_len-1个状态
        for (int i = 0; i < state_len - 1; i++) {
            state[i] = b_norm[i+1] * x[n] - a_norm[i+1] * y[n] + state[i+1];
        }
        // 更新最后一个状态
        double last_b = (nb > state_len) ? b_norm[state_len] : 0.0;
        double last_a = (na > state_len) ? a_norm[state_len] : 0.0;
        state[state_len - 1] = last_b * x[n] - last_a * y[n];
    }

    // 按需返回最终状态
    if (zf != NULL) {
        memcpy(zf, state, state_len * sizeof(double));
    }

    // 释放临时内存
    free(b_norm);
    free(a_norm);
    free(state);
}

验证注意点 #

• 浮点类型对齐：如果Python侧使用单精度float32计算，将代码中double替换为float即可，否则会出现超出浮点精度的误差。
• 简单用例验证：可以先使用简单系数测试，比如b={0.2,0.2,0.2,0.2,0.2}, a={1}的5点滑动平均滤波，对比C实现和scipy.lfilter的输出，双精度下逐点差值应该小于1e-12。
• 分段滤波适配：分段滤波场景下，只需要把上一段滤波返回的zf作为下一段滤波的zi传入，就能得到和整段滤波完全一致的结果，和scipy的zi/zf接口行为完全对齐。

内容的提问来源于stack exchange，提问作者Sushindren S