搜索
C语言

C语言如何高效实现IIR滤波器？

99ANYc3cd6
99ANYc3cd6 / / 0 评论 / 3 阅读
作者有点忙，还没来得及写简介......
预计阅读时长 31 分钟
位置： 首页 › C语言 › 正文

IIR滤波器的理论基础

IIR滤波器的差分方程可以表示为：

y[n] = Σ (b_k * x[n-k]) - Σ (a_k * y[n-k])

  • y[n] 是当前的输出样本。
  • x[n] 是当前的输入样本。
  • b_k 是前向（Feedforward）系数，也称为“零点”系数。
  • a_k 是后向（Feedback）系数，也称为“极点”系数。
  • N 是滤波器的阶数（通常指分母的阶数）。

这个公式可以这样理解：

  1. 前向部分（Feedforward）: 对输入信号 x 进行加权求和，权重由 b_k 决定，这部分是FIR滤波器的核心。
  2. 后向部分（Feedback）: 对过去的输出信号 y 进行加权求和，权重由 a_k 决定，然后从结果中减去这部分，这是IIR滤波器的关键。

注意: 在很多实现中，分母的 a_0 会被归一化为1。a_0 不为1，我们需要先将所有系数除以 a_0。

实现方法：直接I型（Direct Form I）

直接I型结构是直接根据上述差分方程实现的,它需要两组延迟线（缓存）：

  • 一组用于存储过去的输入值 x[n-1], x[n-2], ...。
  • 一组用于存储过去的输出值 y[n-1], y[n-2], ...。

这是最直观、最容易理解的结构，但对于高阶滤波器，它对系数的量化误差比较敏感。

C语言实现步骤

步骤 1: 定义滤波器结构体

为了封装滤波器的状态和系数,我们使用一个结构体 IIRFilter。

#include <stdio.h>
#include <stdlib.h>
#include <string.h>
// 定义IIR滤波器结构体
typedef struct {
    int order;             // 滤波器阶数
    double *b_coeffs;     // 前向系数数组 (b0, b1, b2, ...)
    double *a_coeffs;     // 后向系数数组 (a0, a1, a2, ...), 通常a0=1
    double *x_buffer;      // 输入延迟线 (存储过去的输入值)
    double *y_buffer;      // 输出延迟线 (存储过去的输出值)
    int index;             // 当前延迟线的写入/读取索引
} IIRFilter;

步骤 2: 初始化滤波器

这个函数负责分配内存,并设置初始状态。

// 初始化IIR滤波器
IIRFilter* iir_filter_init(int order, const double *b_coeffs, const double *a_coeffs) {
    IIRFilter *filter = (IIRFilter*)malloc(sizeof(IIRFilter));
    if (!filter) return NULL;
    filter->order = order;
    // 分配并复制系数
    filter->b_coeffs = (double*)malloc((order + 1) * sizeof(double));
    filter->a_coeffs = (double*)malloc((order + 1) * sizeof(double));
    memcpy(filter->b_coeffs, b_coeffs, (order + 1) * sizeof(double));
    memcpy(filter->a_coeffs, a_coeffs, (order + 1) * sizeof(double));
    // 分配延迟线缓冲区，并初始化为0
    filter->x_buffer = (double*)calloc(order + 1, sizeof(double));
    filter->y_buffer = (double*)calloc(order + 1, sizeof(double));
    filter->index = 0;
    return filter;
}

步骤 3: 核心滤波函数

这是处理每个输入样本并生成输出的函数,它实现了直接I型的算法。

// 处理单个样本
double iir_filter_process(IIRFilter *filter, double input) {
    double output = 0.0;
    int i;
    // 1. 将当前输入存入延迟线
    filter->x_buffer[filter->index] = input;
    // 2. 计算输出 y[n] = Σ(b_k * x[n-k]) - Σ(a_k * y[n-k])
    //    注意：a_0 通常为1，所以公式中不包含 a_0 * y[n]
    for (i = 0; i <= filter->order; i++) {
        // 获取延迟线的索引（循环缓冲）
        int x_idx = (filter->index - i + filter->order + 1) % (filter->order + 1);
        int y_idx = (filter->index - i + filter->order + 1) % (filter->order + 1);
        output += filter->b_coeffs[i] * filter->x_buffer[x_idx];
        // 对于反馈部分，跳过 k=0 的情况 (a_0 * y[n])
        if (i > 0) {
            output -= filter->a_coeffs[i] * filter->y_buffer[y_idx];
        }
    }
    // 3. 将当前输出存入延迟线
    filter->y_buffer[filter->index] = output;
    // 4. 更新延迟线索引
    filter->index = (filter->index + 1) % (filter->order + 1);
    return output;
}

注意: index 的更新和缓冲区访问方式构成了一个高效的循环缓冲区，避免了频繁的内存移动。

步骤 4: 释放资源

使用完毕后,必须释放所有分配的内存。

// 释放滤波器资源
void iir_filter_free(IIRFilter *filter) {
    if (filter) {
        free(filter->b_coeffs);
        free(filter->a_coeffs);
        free(filter->x_buffer);
        free(filter->y_buffer);
        free(filter);
    }
}

完整示例：一个2阶低通IIR滤波器

下面是一个完整的、可运行的C程序，它创建一个2阶低通IIR滤波器，并对一个包含高频噪声的输入信号进行滤波。

1 设计滤波器系数

我们将使用一个简单的2阶Butterworth低通滤波器作为例子,其传递函数为： H(s) = 1 / (s^2 + √2*s + 1)

通过双线性变换法,我们可以将其离散化为Z域传递函数，并得到系数，假设我们设定的采样频率 fs 为1000Hz，截止频率 fc 为50Hz。

计算出的系数（归一化后，a_0=1）为：

  • b_0 = 0.0201
  • b_1 = 0.0402
  • b_2 = 0.0201
  • a_1 = -1.5610
  • a_2 = 0.6414

2 C代码 (iir_example.c)

#include <stdio.h>
#include <stdlib.h>
#include <string.h>
#include <math.h> // 用于 sin 和 cos
// 定义IIR滤波器结构体
typedef struct {
    int order;             // 滤波器阶数
    double *b_coeffs;     // 前向系数数组 (b0, b1, b2, ...)
    double *a_coeffs;     // 后向系数数组 (a0, a1, a2, ...), 通常a0=1
    double *x_buffer;      // 输入延迟线 (存储过去的输入值)
    double *y_buffer;      // 输出延迟线 (存储过去的输出值)
    int index;             // 当前延迟线的写入/读取索引
} IIRFilter;
// 初始化IIR滤波器
IIRFilter* iir_filter_init(int order, const double *b_coeffs, const double *a_coeffs) {
    IIRFilter *filter = (IIRFilter*)malloc(sizeof(IIRFilter));
    if (!filter) return NULL;
    filter->order = order;
    // 分配并复制系数
    filter->b_coeffs = (double*)malloc((order + 1) * sizeof(double));
    filter->a_coeffs = (double*)malloc((order + 1) * sizeof(double));
    memcpy(filter->b_coeffs, b_coeffs, (order + 1) * sizeof(double));
    memcpy(filter->a_coeffs, a_coeffs, (order + 1) * sizeof(double));
    // 分配延迟线缓冲区，并初始化为0
    filter->x_buffer = (double*)calloc(order + 1, sizeof(double));
    filter->y_buffer = (double*)calloc(order + 1, sizeof(double));
    filter->index = 0;
    return filter;
}
// 处理单个样本
double iir_filter_process(IIRFilter *filter, double input) {
    double output = 0.0;
    int i;
    // 1. 将当前输入存入延迟线
    filter->x_buffer[filter->index] = input;
    // 2. 计算输出 y[n] = Σ(b_k * x[n-k]) - Σ(a_k * y[n-k])
    for (i = 0; i <= filter->order; i++) {
        // 获取延迟线的索引（循环缓冲）
        int x_idx = (filter->index - i + filter->order + 1) % (filter->order + 1);
        int y_idx = (filter->index - i + filter->order + 1) % (filter->order + 1);
        output += filter->b_coeffs[i] * filter->x_buffer[x_idx];
        // 对于反馈部分，跳过 k=0 的情况 (a_0 * y[n])
        if (i > 0) {
            output -= filter->a_coeffs[i] * filter->y_buffer[y_idx];
        }
    }
    // 3. 将当前输出存入延迟线
    filter->y_buffer[filter->index] = output;
    // 4. 更新延迟线索引
    filter->index = (filter->index + 1) % (filter->order + 1);
    return output;
}
// 释放滤波器资源
void iir_filter_free(IIRFilter *filter) {
    if (filter) {
        free(filter->b_coeffs);
        free(filter->a_coeffs);
        free(filter->x_buffer);
        free(filter->y_buffer);
        free(filter);
    }
}
// 生成测试信号：一个低频正弦波 + 一个高频正弦波（噪声）
double generate_test_signal(int n, double fs) {
    double low_freq = 5.0;  // 5 Hz
    double high_freq = 50.0; // 50 Hz
    double amplitude_low = 1.0;
    double amplitude_high = 0.2;
    return amplitude_low * sin(2 * M_PI * low_freq * n / fs) +
           amplitude_high * sin(2 * M_PI * high_freq * n / fs);
}
int main() {
    // --- 1. 定义滤波器系数 ---
    // 2阶低通滤波器, fs=1000Hz, fc=50Hz
    int order = 2;
    double b_coeffs[] = {0.0201, 0.0402, 0.0201}; // b0, b1, b2
    double a_coeffs[] = {1.0, -1.5610, 0.6414};   // a0, a1, a2 (a0 is always 1)
    // --- 2. 初始化滤波器 ---
    IIRFilter *filter = iir_filter_init(order, b_coeffs, a_coeffs);
    if (!filter) {
        fprintf(stderr, "Failed to initialize filter.\n");
        return 1;
    }
    // --- 3. 准备测试 ---
    double fs = 1000.0; // 采样率 1000 Hz
    int num_samples = 200;
    printf("Sample No.\tInput Signal\tFiltered Output\n");
    printf("------------------------------------------------\n");
    // --- 4. 处理信号 ---
    for (int i = 0; i < num_samples; i++) {
        // 生成带噪声的输入信号
        double input = generate_test_signal(i, fs);
        // 使用滤波器处理信号
        double output = iir_filter_process(filter, input);
        // 打印结果
        printf("%d\t\t%.4f\t\t%.4f\n", i, input, output);
    }
    // --- 5. 清理 ---
    iir_filter_free(filter);
    return 0;
}

3 如何编译和运行

  1. 将代码保存为 iir_example.c。
  2. 使用GCC编译器进行编译：
    gcc iir_example.c -o iir_example -lm
    • -lm 是链接数学库所必需的，因为我们使用了 sin 和 M_PI。
  3. 运行生成的可执行文件：
    ./iir_example

运行后,你会看到一列输入信号和经过滤波后的输出信号，你可以将输出数据导入到Excel、Python的Matplotlib或类似工具中绘制成图表，可以直观地看到高频噪声被显著抑制，而低频信号（5Hz）被保留了下来。

重要注意事项

  1. 稳定性: IIR滤波器的稳定性至关重要，如果极点（由 a_coeffs 决定）位于Z平面的单位圆外，滤波器就不稳定，输出会发散，在设计滤波器时，必须确保其是稳定的，直接I型结构对系数的舍入误差比较敏感，对于高阶滤波器，可以考虑使用直接II型或二阶节 结构，它们更稳定且计算效率更高。

  2. 量化误差: 在使用浮点数（如double）时，这个问题不那么突出，但如果要在嵌入式系统等资源受限的环境中使用定点数（如int或fixed-point），系数和中间计算的量化会严重影响滤波器的频率响应和稳定性。

  3. 系数设计: 手动设计IIR系数非常复杂，通常使用专业的工具（如MATLAB, Python的scipy.signal库）来设计滤波器，然后导出系数，在Python中：

    from scipy import signal
import numpy as np
# 设计一个2阶低通滤波器，截止频率为50Hz，采样频率1000Hz
b, a = signal.butter(2, 50, btype='low', fs=1000)
print("b coefficients:", b)
print("a coefficients:", a)

    这会输出与我们示例中非常接近的系数。

-- 展开阅读全文 --
头像
圆周率在c语言怎么表示
« 上一篇 今天
C语言函数定义的规则和步骤是什么？
下一篇 » 今天

相关文章

取消
微信二维码
支付宝二维码

目录[+]