c++实现巴特沃斯带通滤波器

巴特沃斯带通滤波器是一种常用的滤波器，可以在一定频率范围内滤除不需要的信号，而保留需要的信号。以下是一个用C语言实现的巴特沃斯带通滤波器的示例代码：

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

#define PI 3.14159265358979323846

// 计算阻带衰减系数A
double calc_A(int n, double epsilon) {
    double a = log10(pow(10.0, epsilon) - 1.0);
    double b = log10(pow(10.0, epsilon) + 1.0);
    return n * (a + b) / 2.0;
}

// 计算极点的实部和虚部
void calc_pole(int n, int k, double epsilon, double* re, double* im) {
    double A = calc_A(n, epsilon);
    double theta = PI * (2 * k + 1) / (2 * n);
    *re = -sin(theta) * sinh(A);
    *im = cos(theta) * cosh(A);
}

// 计算零点的实部和虚部
void calc_zero(int n, int k, double* re, double* im) {
    double theta = PI * (2 * k + 1) / (2 * n);
    *re = -sin(theta);
    *im = cos(theta);
}

// 计算巴特沃斯带通滤波器的系数
void calc_coeffs(int n, double epsilon, double f0, double Q, double* b, double* a) {
    double re_pole, im_pole, re_zero, im_zero;
    double K = 1.0;
    int i;

    for (i = 0; i < n / 2; i++) {
        calc_pole(n, i, epsilon, &re_pole, &im_pole);
        calc_zero(n, i, &re_zero, &im_zero);

        double den = 1.0 / (1.0 + Q * re_pole + pow(re_pole, 2.0) + pow(im_pole, 2.0));
        double num_b = K * (Q * re_zero + im_zero);
        double num_a = K * (Q * re_pole + im_pole);

        b[i * 2] = num_b * den;
        b[i * 2 + 1] = -2.0 * num_b * den;
        b[i * 2 + 2] = num_b * den;
        a[i * 2] = 1.0 * den;
        a[i * 2 + 1] = 2.0 * (pow(re_pole, 2.0) - Q * re_pole + pow(im_pole, 2.0)) * den;
        a[i * 2 + 2] = -(1.0 - Q * re_pole + pow(re_pole, 2.0) + pow(im_pole, 2.0)) * den;
    }

    if (n % 2 == 1) {
        calc_pole(n, n / 2, epsilon, &re_pole, &im_pole);
        calc_zero(n, n / 2, &re_zero, &im_zero);

        double den = 1.0 / (1.0 + Q * re_pole + pow(re_pole, 2.0) + pow(im_pole, 2.0));
        double num_b = K * (Q * re_zero + im_zero);
        double num_a = K * (Q * re_pole + im_pole);

        b[n - 1] = num_b * den;
        a[n - 2] = 2.0 * (pow(re_pole, 2.0) - Q * re_pole + pow(im_pole, 2.0)) * den;
        a[n - 1] = -(1.0 - Q * re_pole + pow(re_pole, 2.0) + pow(im_pole, 2.0)) * den;
    }

    double W0 = 2.0 * PI * f0;
    double BW = W0 / Q;
    double gain = pow(W0, n) / (pow(W0, n) + pow(BW, n));
    for (i = 0; i < n; i++) {
        b[i] *= gain;
    }
}

// 应用滤波器
void apply_filter(double* x, int n, double* b, int nb, double* a, int na, double* y) {
    int i, j;

    for (i = 0; i < n; i++) {
        y[i] = 0.0;
        for (j = 0; j < nb; j++) {
            if (i - j >= 0) {
                y[i] += b[j] * x[i - j];
            }
        }

        for (j = 1; j < na; j++) {
            if (i - j >= 0) {
                y[i] -= a[j] * y[i - j];
            }
        }
    }
}

int main() {
    int n = 4;  // 阶数
    double epsilon = 0.2;  // 阻带衰减量
    double f0 = 1000.0;  // 中心频率
    double Q = 1.0 / sqrt(2.0);  // 带宽因子

    double b[n + 1], a[n];
    calc_coeffs(n, epsilon, f0, Q, b, a);

    int i;
    printf("b = [");
    for (i = 0; i <= n; i++) {
        printf("%.15g, ", b[i]);
    }
    printf("]\n");

    printf("a = [");
    for (i = 0; i < n; i++) {
        printf("%.15g, ", a[i]);
    }
    printf("]\n");

    int len = 1000;
    double x[len], y[len];

    // 生成测试信号，包含100Hz和1000Hz两个频率分量
    for (i = 0; i < len; i++) {
        x[i] = sin(2.0 * PI * 100.0 * i / len) + sin(2.0 * PI * 1000.0 * i / len);
    }

    apply_filter(x, len, b, n + 1, a, n, y);

    // 输出滤波后的信号
    for (i = 0; i < len; i++) {
        printf("%g\n", y[i]);
    }

    return 0;
}