4阶巴特沃斯滤波器C语言代码

以下是一个4阶巴特沃斯滤波器的C语言代码实现，其中假设采样率为Fs，截止频率为Fc，通带最大衰减为Ap，阻带最小衰减为As：

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

// 差分方程系数
#define A1 1.4142135623730951
#define A2 1.0
#define A3 0.0

// 计算巴特沃斯滤波器系数
void butterworth(int N, double Ap, double As, double Fc, double Fs, double *a, double *b) {
    double Wp = 2 * M_PI * Fc / Fs;
    double Ws = 2 * M_PI * (Fc + (Fs - Fc) / 2) / Fs; // 阻带截止频率
    double Ep = sqrt(pow(10, Ap / 10) - 1);
    double Es = sqrt(pow(10, As / 10) - 1);
    double K = Wp / Ws;
    double K1 = Ep / Es;
    double K2 = pow(K1, 1.0 / N);
    double K3 = (1 - K2) / (1 + K2);
    double K4 = cos(Wp);
    double K5 = 1 / sqrt(1 + pow(Ep, 2));
    double K6 = K5 * pow(K3, N / 2) * K4;
    double K7 = pow(K3, N);
    double K8 = sqrt(1 - pow(K6, 2)) / K7;
    double K9 = (1 - K7) * K8;
    double K10 = 2 * K6;
    double K11 = pow(K5, 2) * pow(K3, N - 1);
    double K12 = sqrt(1 - pow(K11, 2));
    double K13 = pow(K3, N / 2) * sin(Wp);
    double K14 = K12 / K13;
    double K15 = pow(K14, 2) - pow(K10, 2);
    double K16 = 2 * pow(K10, 2) - pow(K14, 2);
    double K17 = pow(K14, 2) + pow(K10, 2);
    double K18 = K16 / K17;
    double K19 = K15 / K17;
    double K20 = K10 / K17;
    double K21 = 2 * K6 * K18;
    double K22 = pow(K6, 2) - pow(K19, 2) - pow(K20, 2);
    double K23 = 2 * K6 * K19;
    double K24 = pow(K6, 2) - pow(K21, 2) - pow(K20, 2);
    double K25 = 2 * K6 * K20;
    double K26 = pow(K6, 2) - pow(K21, 2) - pow(K19, 2);
    double K27 = pow(K6, 2) - pow(K21, 2) - pow(K19, 2) - pow(K20, 2);

    a[0] = 1;
    a[1] = K21;
    a[2] = K22;
    a[3] = K23;
    a[4] = K24;
    a[5] = K25;
    a[6] = K26;
    a[7] = K27;

    b[0] = K5;
    b[1] = K6 / K7;
    b[2] = K11;
    b[3] = K6 / K7;
    b[4] = K5;
}

// 差分方程计算滤波器输出
double filter(double x, double *a, double *b, double *z, int N) {
    double y = b[0] * x + z[0];
    for (int i = 1; i <= N; i++) {
        z[i - 1] = b[i] * x - a[i] * y + z[i];
    }
    return y;
}

int main(void) {
    double Fs = 1000; // 采样率
    double Fc = 100; // 截止频率
    double Ap = 0.5; // 通带最大衰减
    double As = 40; // 阻带最小衰减
    int N = 4; // 阶数
    double a[N + 1], b[N + 1], z[N]; // 系数和状态变量
    butterworth(N, Ap, As, Fc, Fs, a, b); // 计算系数

    // 输入信号
    double x[] = {1, 2, 3, 4, 5, 6, 7, 8, 9, 10};
    int len = sizeof(x) / sizeof(double);

    // 滤波器输出
    double y[len];
    for (int i = 0; i < len; i++) {
        y[i] = filter(x[i], a, b, z, N);
    }

    // 输出结果
    for (int i = 0; i < len; i++) {
        printf("%f\n", y[i]);
    }

    return 0;
}

该代码实现了一个4阶巴特沃斯滤波器，输入信号为长度为10的实数序列，输出为长度为10的实数序列。在实现中使用了butterworth函数计算巴特沃斯滤波器系数，filter函数使用差分方程计算滤波器输出。