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

巴特沃斯带通滤波器是一种数字信号处理滤波器，可以滤除低于或高于一定频率的信号，同时保留其它频率范围内的信号。以下是一个C语言实现的巴特沃斯带通滤波器的示例代码：

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

#define PI 3.14159265358979323846

//函数原型
void butterworth_bandpass_filter(double *input, double *output, int N, double fs, double f1, double f2, int order);

int main()
{
    //输入信号
    double input[100] = {0.0};
    //输出信号
    double output[100] = {0.0};
    //采样率
    double fs = 1000.0;
    //滤波器通带频率范围
    double f1 = 50.0;
    double f2 = 150.0;
    //滤波器阶数
    int order = 4;
    //信号长度
    int N = 100;

    //生成信号，这里生成了一个10Hz和100Hz的正弦波
    for(int i = 0; i < N; i++)
    {
        input[i] = sin(2 * PI * 10 * i / fs) + sin(2 * PI * 100 * i / fs);
    }

    //调用巴特沃斯带通滤波器函数
    butterworth_bandpass_filter(input, output, N, fs, f1, f2, order);

    //打印输出信号
    for(int i = 0; i < N; i++)
    {
        printf("%lf\n", output[i]);
    }

    return 0;
}

//巴特沃斯带通滤波器函数
void butterworth_bandpass_filter(double *input, double *output, int N, double fs, double f1, double f2, int order)
{
    //定义中间变量
    double w1 = 2 * PI * f1 / fs;
    double w2 = 2 * PI * f2 / fs;
    double bw = w2 - w1;
    double wc = sqrt(w1 * w2);
    double Q = wc / bw;
    double alpha = sin(bw / 2) / (2 * Q);
    double beta = sqrt(1 - alpha * alpha);
    double a0, a1, a2, b1, b2;
    double w0 = wc;

    //计算巴特沃斯滤波器的系数
    switch(order)
    {
        case 1:
            a0 = 1 + alpha;
            a1 = -2 * cos(w0);
            a2 = 1 - alpha;
            b1 = -2 * cos(w0);
            b2 = 1 - alpha;
            break;
        case 2:
            a0 = 1 + 2 * alpha + alpha * alpha;
            a1 = -2 * cos(w0) * (1 + alpha);
            a2 = 1 - 2 * alpha + alpha * alpha;
            b1 = -2 * cos(w0);
            b2 = 1 - alpha;
            break;
        case 3:
            a0 = 1 + 3 * alpha + 3 * alpha * alpha + alpha * alpha * alpha;
            a1 = -2 * cos(w0) * (1 + alpha + alpha * alpha);
            a2 = 1 - 3 * alpha + 3 * alpha * alpha - alpha * alpha * alpha;
            b1 = -2 * cos(w0) * (1 + alpha);
            b2 = 1 - 2 * alpha + alpha * alpha;
            break;
        case 4:
            a0 = 1 + 4 * alpha + 6 * alpha * alpha + 4 * alpha * alpha * alpha + alpha * alpha * alpha * alpha;
            a1 = -2 * cos(w0) * (1 + 2 * alpha + 2 * alpha * alpha + alpha * alpha * alpha);
            a2 = 1 - 4 * alpha + 6 * alpha * alpha - 4 * alpha * alpha * alpha + alpha * alpha * alpha * alpha;
            b1 = -2 * cos(w0) * (1 + alpha + alpha * alpha);
            b2 = 1 - 2 * alpha + alpha * alpha;
            break;
        default:
            printf("The order of filter is not correct!\n");
            return;
    }

    //滤波器初始化
    double x1 = 0.0;
    double x2 = 0.0;
    double y1 = 0.0;
    double y2 = 0.0;

    //滤波器计算
    for(int i = 0; i < N; i++)
    {
        output[i] = (input[i] * a0 + x1 * a1 + x2 * a2 - y1 * b1 - y2 * b2) / a0;
        x2 = x1;
        x1 = input[i];
        y2 = y1;
        y1 = output[i];
    }
}

在上面的代码中，`butterworth_bandpass_filter()`函数是实现巴特沃斯带通滤波器的函数，其中`input`是输入信号，`output`是输出信号，`N`是信号长度，`fs`是采样率，`f1`和`f2`是滤波器通带频率范围，`order`是滤波器阶数。在主函数中，我们生成了一个10Hz和100Hz的正弦波作为输入信号，然后调用`butterworth_bandpass_filter()`函数进行滤波，并打印输出信号。

需要注意的是，巴特沃斯滤波器的阶数越高，滤波器的通带截止频率和阻带截止频率的陡度就越大，但是也会带来计算量的增加和相位失真的问题。在实际应用中需要根据具体需求选择适当的滤波器阶数。