巴特沃斯滤波器 c++代码

以下是一个简单的巴特沃斯滤波器C代码示例：

#include <stdio.h>
#include <math.h>

#define PI 3.14159265359

// 定义巴特沃斯滤波器的阶数和截止频率
#define ORDER 2
#define CUTOFF_FREQ 1000.0

// 计算阻带增益
double calc_stopband_gain(double freq, int order) {
    double omega_c = 2.0 * PI * CUTOFF_FREQ;
    double omega = 2.0 * PI * freq;
    double s = sin(omega);
    double c = cos(omega);
    double s2 = s * s;
    double c2 = c * c;
    double a = 1.0 + pow(omega / omega_c, 2.0 * order);
    double b = 2.0 * pow(omega / omega_c, 2.0 * order) * c;
    double c1 = pow(omega / omega_c, 2.0 * order) * (2.0 * order * s2 - c2 + 1.0);
    double d = 1.0 + pow(omega / omega_c, 2.0 * order) * (2.0 * order * s2 + c2 - 1.0);
    double stopband_gain = sqrt(a * d / (b * b + c1 * d));
    return stopband_gain;
}

// 巴特沃斯滤波器函数
void butterworth_filter(double *input, double *output, int num_samples) {
    double omega_c = 2.0 * PI * CUTOFF_FREQ;
    double a[ORDER+1], b[ORDER+1];
    double gain = calc_stopband_gain(2.0 * CUTOFF_FREQ, ORDER);
    // 计算巴特沃斯滤波器的系数
    a[0] = 1.0;
    a[1] = 2.0 * exp(-PI * omega_c / sqrt(gain)) * cos(sqrt(gain) * PI / 2.0 + PI * omega_c);
    a[2] = exp(-2.0 * PI * omega_c / sqrt(gain));
    b[0] = 1.0 / (1.0 + gain);
    b[1] = 2.0 * b[0];
    b[2] = b[0];
    // 巴特沃斯滤波器的差分方程
    double x[ORDER+1], y[ORDER+1];
    for (int i = 0; i <= ORDER; i++) {
        x[i] = 0.0;
        y[i] = 0.0;
    }
    for (int n = 0; n < num_samples; n++) {
        x[0] = input[n];
        y[0] = b[0] * x[0];
        for (int i = 1; i <= ORDER; i++) {
            y[0] += b[i] * x[i] - a[i] * y[i];
        }
        for (int i = ORDER; i >= 1; i--) {
            x[i] = x[i-1];
            y[i] = y[i-1];
        }
        output[n] = y[0];
    }
}

int main() {
    // 模拟输入信号
    double input[1000];
    for (int i = 0; i < 1000; i++) {
        input[i] = sin(2.0 * PI * i * 100.0 / 1000.0) + sin(2.0 * PI * i * 200.0 / 1000.0);
    }
    // 应用巴特沃斯滤波器
    double output[1000];
    butterworth_filter(input, output, 1000);
    // 输出结果
    for (int i = 0; i < 1000; i++) {
        printf("%f\n", output[i]);
    }
    return 0;
}