使用C语言实现MATLAB同款的巴特沃斯高通滤波器

非常感谢您的提问。以下是使用C语言实现MATLAB同款的巴特沃斯高通滤波器的代码：

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

#define PI 3.14159265358979323846

typedef struct {
    double real;
    double imag;
} complex;

void butterworth_highpass(int n, double fc, double fs, double *b, double *a) {
    int i;
    double w0 = 2 * PI * fc / fs;
    double bw = w0 / 10;
    double alpha = sin(w0) / (2 * sqrt(2) * bw);
    double a0 = 1 + alpha;
    double a1 = -2 * cos(w0);
    double a2 = 1 - alpha;
    double b0 = (1 + cos(w0)) / 2 / a0;
    double b1 = -(1 + cos(w0)) / a0;
    double b2 = (1 + cos(w0)) / 2 / a0;
    a[0] = a0;
    a[1] = a1;
    a[2] = a2;
    b[0] = b0;
    b[1] = b1;
    b[2] = b2;
    for (i = 0; i < n; i++) {
        b[i] /= a0;
        a[i] /= a0;
    }
}

void complex_multiply(complex *a, complex *b, complex *c) {
    c->real = a->real * b->real - a->imag * b->imag;
    c->imag = a->real * b->imag + a->imag * b->real;
}

void complex_add(complex *a, complex *b, complex *c) {
    c->real = a->real + b->real;
    c->imag = a->imag + b->imag;
}

void complex_subtract(complex *a, complex *b, complex *c) {
    c->real = a->real - b->real;
    c->imag = a->imag - b->imag;
}

void complex_divide(complex *a, complex *b, complex *c) {
    double denominator = b->real * b->real + b->imag * b->imag;
    c->real = (a->real * b->real + a->imag * b->imag) / denominator;
    c->imag = (a->imag * b->real - a->real * b->imag) / denominator;
}

void complex_exp(complex *a, complex *b) {
    b->real = exp(a->real) * cos(a->imag);
    b->imag = exp(a->real) * sin(a->imag);
}

void fft(complex *x, int n) {
    int i, j, k;
    complex t, u;
    for (i = 0, j = 0; i < n; i++) {
        if (j > i) {
            t = x[j];
            x[j] = x[i];
            x[i] = t;
        }
        k = n;
        while (j >= k) {
            j -= k;
            k /= 2;
        }
        j += k;
    }
    for (k = 1; k < n; k *= 2) {
        for (j = 0; j < k; j++) {
            t.real = cos(PI * j / k);
            t.imag = -sin(PI * j / k);
            for (i = j; i < n; i += 2 * k) {
                complex_multiply(&t, &x[i + k], &u);
                complex_add(&x[i], &u, &x[i + k]);
                complex_subtract(&x[i], &u, &x[i]);
            }
        }
    }
}

void ifft(complex *x, int n) {
    int i;
    for (i = 0; i < n; i++) {
        x[i].imag = -x[i].imag;
    }
    fft(x, n);
    for (i = 0; i < n; i++) {
        x[i].real /= n;
        x[i].imag /= -n;
    }
}

void filter(double *x, int n, double *b, double *a, double *y) {
    int i, j;
    complex *X = (complex *) malloc(n * sizeof(complex));
    complex *H = (complex *) malloc(n * sizeof(complex));
    complex *Y = (complex *) malloc(n * sizeof(complex));
    for (i = 0; i < n; i++) {
        X[i].real = x[i];
        X[i].imag = 0;
    }
    fft(X, n);
    for (i = 0; i < n; i++) {
        H[i].real = b[0];
        H[i].imag = 0;
        for (j = 1; j <= i && j < 3; j++) {
            complex t;
            complex_multiply(&H[i - j], &X[i], &t);
            complex_add(&t, &H[i], &H[i]);
            complex_multiply(&H[i - j], &X[i - j], &t);
            complex_subtract(&H[i], &t, &H[i]);
        }
        for (j = 1; j < 3 && j < n - i; j++) {
            complex t;
            complex_multiply(&H[i + j], &X[i], &t);
            complex_add(&t, &H[i], &H[i]);
            complex_multiply(&H[i + j], &X[i + j], &t);
            complex_subtract(&H[i], &t, &H[i]);
        }
        complex_divide(&H[i], &X[i], &H[i]);
        complex_exp(&H[i], &H[i]);
    }
    for (i = 0; i < n; i++) {
        complex_multiply(&H[i], &X[i], &Y[i]);
    }
    ifft(Y, n);
    for (i = 0; i < n; i++) {
        y[i] = Y[i].real * a[0];
        for (j = 1; j <= i && j < 3; j++) {
            y[i] -= y[i - j] * a[j];
        }
        for (j = 1; j < 3 && j < n - i; j++) {
            y[i] -= y[i + j] * a[j];
        }
    }
    free(X);
    free(H);
    free(Y);
}

int main() {
    int i, n = 1000;
    double *x = (double *) malloc(n * sizeof(double));
    double *y = (double *) malloc(n * sizeof(double));
    double *b = (double *) malloc(3 * sizeof(double));
    double *a = (double *) malloc(3 * sizeof(double));
    double fc = 50, fs = 1000;
    for (i = 0; i < n; i++) {
        x[i] = sin(2 * PI * i * fc / fs);
    }
    butterworth_highpass(2, fc, fs, b, a);
    filter(x, n, b, a, y);
    for (i = 0; i < n; i++) {
        printf("%f\n", y[i]);
    }
    free(x);
    free(y);
    free(b);
    free(a);
    return 0;
}

这段代码实现了一个二阶巴特沃斯高通滤波器，可以通过调整参数来实现不同的滤波效果。