如何用C语言实现MATLAB的filtfilt滤波功能？

下面我将详细解释 filtfilt 的工作原理，并提供两种 C 语言实现方案：一种是基于 IIR 滤波器的直接实现（但仍有相位延迟），另一种是推荐使用的、真正实现零相位滤波的 FFT 方法。

MATLAB filtfilt 的工作原理

filtfilt 的核心目标是消除相位延迟，普通的 filter 函数（无论是 FIR 还是 IIR）都会在滤波过程中引入相位延迟，导致输出信号相对于输入信号在时间上发生偏移。

filtfilt 通过一种巧妙的方法实现了零相位滤波，其步骤如下：

1. 正向滤波：将输入信号 x 通过滤波器 H(z)，得到一个中间结果 y1。 y1 = filter(b, a, x);

2. 信号反转：将中间结果 y1 的时间轴反转（即从后往前播放），得到 y1_rev。 y1_rev = flip(y1);

   
   （图片来源网络，侵删）

3. 反向滤波：将 y1_rev 再次通过同一个滤波器 H(z)，得到另一个中间结果 y2。 y2 = filter(b, a, y1_rev);

4. 结果反转：将 y2 的时间轴再次反转，得到最终的输出信号 y。 y = flip(y2);

为什么这样做能消除相位延迟？ 可以这样理解：正向滤波引入的相位延迟，在信号反转后进行反向滤波时，会以相反的方式被“抵消”掉，两次反转操作保证了最终输出信号的时间顺序与输入一致，而两次滤波则保证了幅频响应的正确性。

重要前提：这种方法要求数字滤波器是零相位的，对于 IIR 滤波器（如 Butterworth, Chebyshev），它们本身不是零相位的。filtfilt 通过上述“前向-后向”滤波法，为 IIR 滤波器构造了一个等效的零相位滤波器，但对于 FIR 滤波器，如果其系数本身就是对称的（线性相位），则可以直接使用。

C 语言实现方案

在 C 语言中，我们通常不会从头实现所有滤波算法，而是依赖于成熟的数学库，最常用的是 FFTW (The Fastest Fourier Transform in the West)，它用于快速傅里叶变换。

直接模拟 filtfilt 步骤（适用于 IIR 滤波器）

这种方法直接将 MATLAB 的四个步骤翻译成 C 语言，你需要一个 IIR 滤波器的实现库，

• DSP libraries: 像 Intel's IPP (Integrated Performance Primitives) 或 ARM's CMSIS-DSP 提供了高效的 IIR 滤波函数。
• DIY implementation: 自己用直接I型或直接II型结构实现，但要注意数值稳定性。

伪代码/C 语言风格逻辑：

#include <stdlib.h> // for malloc, free
#include <string.h> // for memcpy
// 假设你有一个IIR滤波器函数
// void iir_filter(const double* b, int nb, const double* a, int na, 
//                 const double* x, int n, double* y);
// b: 分子系数, a: 分母系数, x: 输入, n: 数据点数, y: 输出
void filtfilt_like_iir(const double* b, int nb, const double* a, int na,
                       const double* x, int n, double* y) {
    double* y1 = (double*)malloc(n * sizeof(double));
    double* y1_rev = (double*)malloc(n * sizeof(double));
    double* y2 = (double*)malloc(n * sizeof(double));
    if (!y1 || !y1_rev || !y2) {
        // 处理内存分配失败
        free(y1); free(y1_rev); free(y2);
        return;
    }
    // 1. 正向滤波
    iir_filter(b, nb, a, na, x, n, y1);
    // 2. 信号反转
    for (int i = 0; i < n; i++) {
        y1_rev[i] = y1[n - 1 - i];
    }
    // 3. 反向滤波
    iir_filter(b, nb, a, na, y1_rev, n, y2);
    // 4. 结果反转
    for (int i = 0; i < n; i++) {
        y[i] = y2[n - 1 - i];
    }
    free(y1);
    free(y1_rev);
    free(y2);
}

缺点：

• 这种方法在时域进行,计算量可能很大，特别是对于长信号和高阶滤波器。
• IIR 滤波器的数值稳定性在有限精度（如 float）的 C 语言环境中可能比在 MATLAB 中更容易出现问题。

基于 FFT 的零相位滤波（推荐方案）

这是 filtfilt 在 MATLAB 中更高效的实现方式，也是 C 语言中实现零相位滤波的最佳实践，其核心思想是：在频域中，滤波就是乘法，而乘法是可交换的，没有顺序之分，因此没有相位延迟。

步骤如下：

1. 将滤波器系数转换为频域响应：

   • 计算滤波器系数 b 的频率响应 H(f)，这通常通过计算 b 的 DTFT（离散时间傅里叶变换）或 FFT 来完成。filtfilt 会自动根据输入信号的长度对滤波器进行零填充，以确保在频域卷积时循环卷积等于线性卷积。

2. 对输入信号进行 FFT：

   • 将输入信号 x 补零到合适的长度（通常是下一个2的幂次方，以优化FFT性能），然后计算其 FFT，得到频域表示 X(f)。

3. 频域相乘：

   • 在频域中,将信号的频谱 X(f) 与滤波器的频响 H(f) 逐点相乘，得到输出信号的频谱 Y(f)。
   • Y(f) = X(f) * H(f)

4. 逆 FFT (IFFT)：

   • 对 Y(f) 进行逆 FFT，将结果转换回时域，得到最终的零相位滤波输出 y。

C 语言实现（使用 FFTW 库）

FFTW 是一个用于计算一维、多维离散傅里叶变换的 C 语言库，以其速度极快而闻名。

准备工作：

1. 下载并安装 FFTW：从 FFTW 官网 下载并按照说明安装，你需要包含头文件 #include <fftw3.h> 并链接 FFTW 库。

示例代码：

#include <stdio.h>
#include <stdlib.h>
#include <math.h>
#include <fftw3.h>
// 计算滤波器的频率响应 H(f)
// b: 滤波器分子系数, nb: 阶数+1
// N_fft: FFT 的长度
void get_filter_response(const double* b, int nb, int N_fft, fftw_complex* H) {
    fftw_plan plan;
    double* b_padded = (double*)fftw_malloc(sizeof(double) * N_fft);
    fftw_complex* b_fft = (fftw_complex*)fftw_malloc(sizeof(fftw_complex) * N_fft);
    // 1. 将滤波器系数 b 补零到 N_fft 长度
    memset(b_padded, 0, sizeof(double) * N_fft);
    memcpy(b_padded, b, sizeof(double) * nb);
    // 2. 计算 b 的 FFT，得到 H(f)
    plan = fftw_plan_dft_r2c_1d(N_fft, b_padded, b_fft, FFTW_ESTIMATE);
    fftw_execute(plan);
    // 3. 归一化并存储结果到 H
    // 注意：FFW 的 r2c 输出长度为 N/2 + 1
    for (int i = 0; i <= N_fft / 2; i++) {
        H[i][0] = b_fft[i][0] / (double)N_fft; // 实部
        H[i][1] = b_fft[i][1] / (double)N_fft; // 虚部
    }
    fftw_destroy_plan(plan);
    fftw_free(b_padded);
    fftw_free(b_fft);
}
// 基于 FFT 的零相位滤波函数
// x: 输入信号, n: 信号长度
// b: 滤波器分子系数, nb: 阶数