数字滤波器的c语言实现

核心概念：数字滤波器是什么？

数字滤波器是一种对输入的数字序列（信号）进行运算，从而改变其频率特性，得到期望输出数字序列的算法或设备，它的核心思想是差分方程。

最常见的一种是无限脉冲响应和有限脉冲响应滤波器,它们的实现方式截然不同。

• FIR (Finite Impulse Response) - 有限脉冲响应滤波器

  • 特点：输出仅取决于当前的输入和过去的有限个输入样本。
  • 结构：没有反馈回路，像一条“流水线”。
  • 优点：总是稳定的,易于实现线性相位。
  • 缺点：通常需要更多的系数（阶数）才能达到与IIR相同的频率选择性。

• IIR (Infinite Impulse Response) - 无限脉冲响应滤波器

  • 特点：输出不仅取决于输入和过去的输入，还取决于过去的输出（有反馈）。
  • 结构：包含反馈回路。
  • 优点：可以用较少的系数（阶数）实现很好的频率选择性。
  • 缺点：可能不稳定,且难以实现线性相位。

FIR滤波器的C语言实现

FIR滤波器的实现非常直观,其差分方程为：

y[n] = b₀x[n] + b₁x[n-1] + b₂x[n-2] + ... + bₙx[n-N]

• y[n] 是当前输出。
• x[n] 是当前输入。
• x[n-1], x[n-2] 是过去1个、2个...的输入。
• b₀, b₁, ..., bₙ 是滤波器的系数。
• N 是滤波器的阶数。

实现方法：环形缓冲区

为了高效地存储过去的输入样本，使用环形缓冲区 是最佳选择，它避免了在每次新样本到来时移动整个数据数组,只需要更新一个指针即可。

完整代码示例：FIR低通滤波器

#include <stdio.h>
#include <stdint.h>
// 定义滤波器参数
#define FIR_ORDER 32  // 滤波器阶数 (必须是2的幂-1，方便取模，但非必须)
#define SAMPLE_RATE 44100 // 采样率 44.1kHz
#define CUTOFF_FREQ 5000  // 截止频率 5kHz
// FIR滤波器系数 (这里为了演示，使用一个简单的窗函数法生成的系数)
// 在实际应用中，应该使用MATLAB, Python (scipy)或滤波器设计工具生成
const float fir_coeffs[FIR_ORDER + 1] = {
    // ... (这里应该是FIR_ORDER+1个系数，示例代码中省略，实际使用时请填充)
    // 一个简单的低通滤波器系数可能如下（非最优，仅作示例）
    0.0021, 0.0035, 0.0060, 0.0094, 0.0138, 0.0190, 0.0250, 0.0316,
    0.0386, 0.0457, 0.0528, 0.0596, 0.0659, 0.0714, 0.0760, 0.0795,
    0.0820, 0.0835, 0.0840, 0.0835, 0.0820, 0.0795, 0.0760, 0.0714,
    0.0659, 0.0596, 0.0528, 0.0457, 0.0386, 0.0316, 0.0250, 0.0190,
    0.0138, 0.0094, 0.0060, 0.0035, 0.0021
};
// FIR滤波器结构体
typedef struct {
    float buffer[FIR_ORDER + 1]; // 存储输入样本的环形缓冲区
    uint16_t index;             // 当前写入位置的索引
} FIRFilter;
// 初始化FIR滤波器
void fir_init(FIRFilter* fir) {
    for (int i = 0; i < FIR_ORDER + 1; i++) {
        fir->buffer[i] = 0.0f;
    }
    fir->index = 0;
}
// 执行FIR滤波
float fir_filter(FIRFilter* fir, float input_sample) {
    float output = 0.0f;
    // 1. 将新样本存入环形缓冲区
    fir->buffer[fir->index] = input_sample;
    // 2. 计算卷积 (加权求和)
    // 从当前索引开始，反向遍历历史样本
    for (int i = 0; i < FIR_ORDER + 1; i++) {
        // 计算历史样本的索引 (环形)
        int history_index = (fir->index - i + (FIR_ORDER + 1)) % (FIR_ORDER + 1);
        output += fir_coeffs[i] * fir->buffer[history_index];
    }
    // 3. 更新索引，为下一次输入做准备
    fir->index = (fir->index + 1) % (FIR_ORDER + 1);
    return output;
}
// --- 主函数演示 ---
int main() {
    FIRFilter my_filter;
    fir_init(&my_filter);
    // 模拟一个输入信号 (一个10kHz的正弦波)
    float input_signal;
    // 假设我们处理1024个样本
    for (int i = 0; i < 1024; i++) {
        input_signal = sin(2 * 3.14159 * 10000 * i / SAMPLE_RATE);
        // 滤波
        float filtered_output = fir_filter(&my_filter, input_signal);
        // 打印部分结果观察
        if (i % 100 == 0) {
            printf("Input: %f, Output: %f\n", input_signal, filtered_output);
        }
    }
    return 0;
}

IIR滤波器的C语言实现

IIR滤波器的差分方程为：

y[n] = b₀x[n] + b₁x[n-1] + ... - a₁y[n-1] - a₂y[n-2] - ...

可以看到，输出 y[n] 的计算需要用到过去的输出值 y[n-1], y[n-2] 等,这构成了反馈回路。

实现方法：直接I型或直接II型

直接II型结构更节省内存（只需要一组延迟单元）,是实现IIR滤波器的标准方法。

完整代码示例：IIR低通滤波器

我们以一个二阶IIR滤波器（也称为双二节滤波器，Biquad）为例,它是构建高阶IIR滤波器的基本模块。

#include <stdio.h>
#include <stdint.h>
// 定义滤波器参数
#define IIR_ORDER 2 // 二阶滤波器
// IIR滤波器系数 (示例: 一个低通滤波器)
// 这些系数通常通过MATLAB的 `butter`, `cheby1` 等函数设计得到
const float b_coeffs[3] = { 0.2929, 0.5858, 0.2929 }; // b0, b1, b2
const float a_coeffs[3] = { 1.0000, -0.1716, 0.2500 }; // a0, a1, a2 (注意a0通常为1)
// IIR滤波器结构体 (直接II型实现)
typedef struct {
    float x_history[IIR_ORDER]; // 存储过去的输入 x[n-1], x[n-2]
    float y_history[IIR_ORDER]; // 存储过去的输出 y[n-1], y[n-2]
} IIRFilter;
// 初始化IIR滤波器
void iir_init(IIRFilter* iir) {
    for (int i = 0; i < IIR_ORDER; i++) {
        iir->x_history[i] = 0.0f;
        iir->y_history[i] = 0.0f;
    }
}
// 执行IIR滤波
float iir_filter(IIRFilter* iir, float input_sample) {
    float output = 0.0f;
    // 1. 计算分子部分 (前馈): b0*x[n] + b1*x[n-1] + b2*x[n-2]
    output = b_coeffs[0] * input_sample;
    for (int i = 1; i <= IIR_ORDER; i++) {
        output += b_coeffs[i] * iir->x_history[i - 1];
    }
    // 2. 计算分母部分 (反馈): - a1*y[n-1] - a2*y[n-2]
    for (int i = 1; i <= IIR_ORDER; i++) {
        output -= a_coeffs[i] * iir->y_history[i - 1];
    }
    // 3. 更新历史数据
    // 将当前输入和输出存入历史记录，为下一次计算做准备
    for (int i = IIR_ORDER - 1; i > 0; i--) {
        iir->x_history[i] = iir->x_history[i - 1];
        iir->y_history[i] = iir->y_history[i - 1];
    }
    iir->x_history[0] = input_sample;
    iir->y_history[0] = output;
    return output;
}
// --- 主函数演示 ---
int main() {
    IIRFilter my_filter;
    iir_init(&my_filter);
    // 模拟一个输入信号 (一个10kHz的正弦波)
    float input_signal;
    // 假设我们处理1024个样本
    for (int i = 0; i < 1024; i++) {
        input_signal = sin(2 * 3.14159 * 10000 * i / 44100);
        // 滤波
        float filtered_output = iir_filter(&my_filter, input_signal);
        // 打印部分结果观察
        if (i % 100 == 0) {
            printf("Input: %f, Output: %f\n", input_signal, filtered_output);
        }
    }
    return 0;
}

关键注意事项与优化

1. 定点数 vs. 浮点数

   • 浮点数：实现简单，代码清晰，但计算量大,在资源受限的微控制器上可能性能不足。
   • 定点数：使用整数来模拟小数（Q15格式：1位符号位，15位小数），计算速度快，占用资源少，但需要仔细处理溢出和量化噪声，代码更复杂，在大多数现代ARM Cortex-M系列MCU上，浮点运算单元已非常普及,浮点实现是首选。

2. 滤波器系数的生成

   • 手动设计滤波器系数非常困难且不精确，你应该使用专业的工具：
     • MATLAB / Octave: 使用 fir1, butter, cheby1 等函数，然后使用 fi 工具箱进行定点数转换。
     • Python (SciPy): 使用 scipy.signal 模块,功能与MATLAB类似。
     • 在线工具: 有很多在线的滤波器设计器。

3. 性能优化

   • 查表法：对于某些特定滤波器（如梳状滤波器），可以将计算结果预存入查找表中，运行时直接查表,速度极快。
   • SIMD指令: 如果你的CPU支持（如ARM的NEON，x86的SSE/AVX）,可以使用单指令多数据流指令集并行处理多个样本。
   • 汇编优化: 对性能要求极高的核心部分,可以手动编写汇编代码。

4. 稳定性 (仅针对IIR)

   IIR滤波器的稳定性至关重要，一个不稳定的滤波器会使输出迅速发散到无穷大，在设计IIR滤波器时，必须确保其所有极点都在单位圆内,使用专业设计工具可以保证这一点。

选择FIR还是IIR，取决于你的具体应用需求：如果对相位敏感（如音频），优先选FIR；如果对计算资源和内存要求苛刻，且可以接受非线性相位，优先选IIR,希望这些详细的解释和代码能帮助你理解并实现数字滤波器。