c++实现巴特沃斯带通滤波器完整代码

由于巴特沃斯带通滤波器需要指定截止频率和阶数，因此代码可能会根据具体需求有所不同。以下是一种实现巴特沃斯带通滤波器的完整代码示例：

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

#define PI 3.14159265358979323846

typedef struct {
    int order;          // 阶数
    double *a_coeffs;   // 分母系数
    double *b_coeffs;   // 分子系数
} Filter;

double *allocate_array(int n) {
    double *arr = (double*)malloc(n * sizeof(double));
    if (arr == NULL) {
        printf("Memory allocation failed!\n");
        exit(1);
    }
    return arr;
}

void free_array(double *arr) {
    free(arr);
}

Filter *allocate_filter(int order) {
    Filter *f = (Filter*)malloc(sizeof(Filter));
    if (f == NULL) {
        printf("Memory allocation failed!\n");
        exit(1);
    }
    f->order = order;
    f->a_coeffs = allocate_array(order);
    f->b_coeffs = allocate_array(order + 1);
    return f;
}

void free_filter(Filter *f) {
    free_array(f->a_coeffs);
    free_array(f->b_coeffs);
    free(f);
}

void calc_butterworth_coeffs(double *a_coeffs, double *b_coeffs, int order, double f1, double f2, double fs) {
    double bw = f2 - f1;
    double wc1 = 2.0 * PI * f1 / fs;
    double wc2 = 2.0 * PI * f2 / fs;
    double Wc1 = tan(wc1 / 2.0);
    double Wc2 = tan(wc2 / 2.0);
    double Bw = Wc2 - Wc1;
    double m = order / 2.0;
    double *q_coeffs = allocate_array(order + 1);
    double *p_coeffs = allocate_array(order + 1);
    double a = sin(PI * m);
    double b = sinh(1.0 / m * asinh(1.0 / a));
    double c = cosh(1.0 / m * asinh(1.0 / a));
    double d = sqrt(1.0 + pow(b, 2.0));
    double e = sqrt(1.0 + 1.0 / pow(b, 2.0));
    double f = pow(c, 2.0 * m);
    double g = pow(bw / Bw, 1.0 / (2.0 * m));
    double h = pow(d, 2.0 * m);
    double i = pow(e, 2.0 * m);
    double j = pow(g, 2.0);
    double k = 1.0 / (1.0 + j * i / h);
    double l = 2.0 * k;
    for (int n = 0; n <= order; n++) {
        double s = sin((2.0 * n + 1.0) * PI / (2.0 * order));
        double u = pow(s, 2.0);
        double v = pow(j, m);
        double w = 1.0 + l * v * pow(e, 2.0 * m) / h * u;
        double x = pow(bw / (2.0 * Bw) * s, 2.0);
        q_coeffs[n] = 1.0 / w;
        p_coeffs[n] = x / w;
    }
    double prod_q = 1.0;
    double prod_p = 1.0;
    for (int n = 0; n <= order; n++) {
        prod_q *= q_coeffs[n];
        prod_p *= p_coeffs[n];
    }
    double K = pow(prod_p, 1.0 / (2.0 * order)) / prod_q;
    b_coeffs[0] = K;
    for (int n = 1; n <= order; n++) {
        double sum_a = 0.0;
        double sum_b = 0.0;
        for (int k = 0; k < n; k++) {
            sum_a += a_coeffs[k] * b_coeffs[n - k - 1];
            sum_b += a_coeffs[k] * a_coeffs[n - k - 1];
        }
        b_coeffs[n] = K * sum_a / sum_b;
        a_coeffs[n - 1] = -b_coeffs[n] * sum_b / a_coeffs[0];
    }
    a_coeffs[order] = K;
    // Normalize coefficients
    double sum_sq_b = 0.0;
    for (int n = 0; n <= order; n++) {
        sum_sq_b += pow(b_coeffs[n], 2.0);
    }
    double norm_factor = 1.0 / sqrt(sum_sq_b);
    for (int n = 0; n <= order; n++) {
        b_coeffs[n] *= norm_factor;
    }
    free_array(q_coeffs);
    free_array(p_coeffs);
}

void apply_filter(double *input, double *output, int n, Filter *f) {
    for (int i = 0; i < n; i++) {
        double y = f->b_coeffs[0] * input[i];
        for (int j = 1; j <= f->order; j++) {
            if (i >= j) {
                y += f->b_coeffs[j] * input[i - j] - f->a_coeffs[j - 1] * output[i - j];
            }
        }
        output[i] = y;
    }
}

int main() {
    int order = 4;
    double f1 = 10.0;
    double f2 = 100.0;
    double fs = 1000.0;
    Filter *f = allocate_filter(order);
    calc_butterworth_coeffs(f->a_coeffs, f->b_coeffs, order, f1, f2, fs);
    // Generate test signal
    int n = 1000;
    double *input = allocate_array(n);
    double *output = allocate_array(n);
    for (int i = 0; i < n; i++) {
        input[i] = sin(2.0 * PI * 50.0 * i / fs) + sin(2.0 * PI * 200.0 * i / fs);
    }
    apply_filter(input, output, n, f);
    // Print output to file
    FILE *fp = fopen("output.txt", "w");
    if (fp == NULL) {
        printf("Cannot open file!\n");
        exit(1);
    }
    for (int i = 0; i < n; i++) {
        fprintf(fp, "%f\n", output[i]);
    }
    fclose(fp);
    free_array(input);
    free_array(output);
    free_filter(f);
    return 0;
}

此代码中，`allocate_array` 和 `free_array` 函数用于分配和释放动态数组内存；`Filter` 结构体用于存储滤波器的阶数和系数；`allocate_filter` 和 `free_filter` 函数用于分配和释放滤波器内存；`calc_butterworth_coeffs` 函数用于计算巴特沃斯带通滤波器的系数；`apply_filter` 函数用于应用滤波器到输入信号上，生成输出信号。`main` 函数中，首先指定巴特沃斯滤波器的阶数和截止频率，然后计算系数，接着生成一个包含两个正弦波的测试信号，应用滤波器并将输出信号写入到文件中。