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

巴特沃斯带通滤波器是一种数字滤波器，它可以通过改变滤波器的截止频率和阶数来调整其频率响应。以下是一种用C语言实现的巴特沃斯带通滤波器的示例代码：

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

#define PI 3.14159265358979323846

// 求解倒数
double reciprocal(double x) {
    return 1.0 / x;
}

// 求解复数的实部
double real(double complex x) {
    return creal(x);
}

// 求解复数的虚部
double imag(double complex x) {
    return cimag(x);
}

// 求解巴特沃斯带通滤波器的系数
void butterworth_bandpass(int order, double f_low, double f_high, double sample_rate, double *b_coef, double *a_coef) {
    double complex s[order * 2], z[order * 2];
    double Wn_low = 2 * PI * f_low / sample_rate, Wn_high = 2 * PI * f_high / sample_rate;
    double complex Hs, Hz;
    double a = 0, b = 0, c = 0, d = 0;
    int i, j;

    // 计算极点
    for (i = 1; i <= order; i++) {
        s[i - 1] = cexp(I * PI / (2 * order) * (2 * i - 1 + order));
        s[2 * order - i] = conj(s[i - 1]);
    }

    // 变换到数字域
    for (i = 0; i < order * 2; i++) {
        z[i] = (2 * PI * sample_rate + s[i]) / (2 * PI * sample_rate - s[i]);
    }

    // 计算系数
    for (i = 0; i < order * 2; i++) {
        Hs = 1;
        for (j = 0; j < order * 2; j++) {
            if (j != i) {
                Hs *= (s[i] - s[j]);
            }
        }
        Hz = 1;
        for (j = 0; j < order * 2; j++) {
            if (j != i) {
                Hz *= (z[i] - z[j]);
            }
        }
        a += real(Hs / Hz);
        b += real((Wn_high + Wn_low) * Hs / Hz);
        c += real(Wn_high * Wn_low * Hs / Hz);
        d += real(Hz / Hs);
    }
    a /= d;
    b /= d;
    c /= d;
    d = 1;

    // 归一化系数
    for (i = 0; i <= order; i++) {
        a_coef[i] = 0;
        b_coef[i] = 0;
    }
    b_coef[0] = (c + b + 1) / 2;
    b_coef[1] = -(c - b + 1) / 2;
    b_coef[2] = (c + b + 1) / 2;
    a_coef[0] = 1;
    a_coef[1] = -(2 * a - 2) / (c + b + 1);
    a_coef[2] = (a - b + c - 1) / (c + b + 1);
}

// 应用巴特沃斯带通滤波器
void apply_butterworth_bandpass(double *data, int length, int order, double f_low, double f_high, double sample_rate) {
    double b_coef[3], a_coef[3];
    double x[order + 1], y[order + 1];
    int i, j;

    // 计算滤波器系数
    butterworth_bandpass(order, f_low, f_high, sample_rate, b_coef, a_coef);

    // 初始化状态变量
    for (i = 0; i <= order; i++) {
        x[i] = 0;
        y[i] = 0;
    }

    // 应用滤波器
    for (i = 0; i < length; i++) {
        x[0] = data[i];
        y[0] = b_coef[0] * x[0] + b_coef[1] * x[1] + b_coef[2] * x[2]
            - a_coef[1] * y[1] - a_coef[2] * y[2];
        for (j = order; j >= 1; j--) {
            x[j] = x[j - 1];
            y[j] = y[j - 1];
        }
        data[i] = y[0];
    }
}

int main() {
    double data[] = {0.1, 0.2, 0.3, 0.4, 0.5, 0.6, 0.7, 0.8, 0.9, 1.0};
    int length = sizeof(data) / sizeof(double);
    int order = 2;
    double f_low = 0.2, f_high = 0.3, sample_rate = 1.0;
    int i;

    // 应用巴特沃斯带通滤波器
    apply_butterworth_bandpass(data, length, order, f_low, f_high, sample_rate);

    // 输出滤波后的数据
    for (i = 0; i < length; i++) {
        printf("%f ", data[i]);
    }
    printf("\n");

    return 0;
}

该代码实现了一个二阶巴特沃斯带通滤波器，可以通过修改`order`、`f_low`、`f_high`和`sample_rate`来调整滤波器的阶数、截止频率和采样率。在`main`函数中，我们定义了一个包含10个数据点的测试数据，并使用该滤波器对其进行滤波，最后输出滤波后的数据。