巴特沃斯滤波器 C/C++代码

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

// 巴特沃斯滤波器
// 输入：x（输入信号数组）
// 输出：y（输出信号数组）
// N（信号长度）
// fc（截止频率）
// fs（采样频率）
// order（滤波器阶数）

void butterworth_filter(double *x, double *y, int N, double fc, double fs, int order) {
    int i, j;
    double pi = 3.14159265358979323846;
    double wc = 2.0 * pi * fc / fs;
    double *b = (double*)malloc((order+1)*sizeof(double));
    double *a = (double*)malloc((order+1)*sizeof(double));
    double *xn = (double*)malloc((order+1)*sizeof(double));
    double *yn = (double*)malloc((order+1)*sizeof(double));
    double *tmp;

    // 计算巴特沃斯滤波器系数
    double c = 1.0 / tan(wc/2.0);
    double d = pow(c, 2.0*order);
    for (i = 0; i <= order; i++) {
        double num = 1.0;
        double den = 1.0;
        for (j = 0; j < order; j++) {
            if (i == j) {
                num *= 1.0;
                den *= 1.0;
            } else {
                num *= 2.0*c*cos((pi/2.0) + ((pi*i) + (pi*j))/(2.0*order));
                den *= pow(c, 2.0) - 2.0*c*cos((pi*i - pi*j)/(2.0*order)) + 1.0;
            }
        }
        b[i] = num / den;
        a[i] = d / den;
    }

    // 零相延迟滤波器
    for (i = 0; i < N; i++) {
        xn[0] = x[i];
        yn[0] = b[0] * xn[0];
        for (j = 1; j <= order; j++) {
            xn[j] = x[i-j+1];
            yn[0] += b[j] * xn[j] - a[j] * yn[j];
            yn[j] = yn[j-1];
        }
        y[i] = yn[0];
    }

    // 释放内存
    free(b);
    free(a);
    free(xn);
    free(yn);
}

这段代码实现了一个零相延迟的巴特沃斯滤波器。输入为信号数组x、信号长度N、截止频率fc、采样频率fs和滤波器阶数order。输出为滤波后的信号数组y。在函数体内，先根据截止频率和采样频率计算出滤波器系数b和a，然后使用零相延迟的方式进行滤波。最后释放内存。