matlab中的butter滤波器C语言

以下是一个简单的 C 语言程序，演示如何使用 butter 滤波器函数：

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

/* Butterworth lowpass filter */
void butter_lowpass_filter(double *input, double *output, int size, double cutoff_freq, double sampling_freq, int filter_order)
{
    double *work = (double*) malloc(sizeof(double) * size);
    double norm_cutoff_freq = cutoff_freq / sampling_freq;
    double x, y;
    int i, j;

    /* Compute poles of Butterworth polynomial */
    double *poles_real = (double*) malloc(sizeof(double) * filter_order);
    double *poles_imag = (double*) malloc(sizeof(double) * filter_order);
    double theta = M_PI / filter_order;
    for (i = 0; i < filter_order; i++) {
        poles_real[i] = -sin((2*i+1)*theta/2.0);
        poles_imag[i] = cos((2*i+1)*theta/2.0);
    }

    /* Warp poles from s-plane to z-plane using bilinear transform */
    double *pz_real = (double*) malloc(sizeof(double) * filter_order);
    double *pz_imag = (double*) malloc(sizeof(double) * filter_order);
    for (i = 0; i < filter_order; i++) {
        double a = 1.0 + poles_real[i];
        pz_real[i] = (1.0 - poles_real[i]) / a;
        pz_imag[i] = poles_imag[i] / a;
        for (j = 1; j < filter_order; j++) {
            double aj = 1.0 + poles_real[j];
            pz_real[i] *= (1.0 - poles_real[j]) / aj;
            pz_imag[i] *= poles_imag[j] / aj;
        }
        pz_real[i] = norm_cutoff_freq * sqrt(pz_real[i]*pz_real[i] + pz_imag[i]*pz_imag[i]);
        pz_imag[i] = 0.0;
    }

    /* Compute filter coefficients using inverse DFT */
    double *coeffs = (double*) malloc(sizeof(double) * filter_order);
    for (i = 0; i < filter_order; i++) {
        coeffs[i] = 0.0;
        for (j = 0; j < filter_order; j++) {
            x = 2.0 * M_PI * i * j / filter_order;
            y = pz_real[j];
            coeffs[i] += y * cos(x);
        }
        coeffs[i] /= filter_order;
    }

    /* Apply filter to input signal */
    for (i = 0; i < size; i++) {
        work[i] = input[i];
        for (j = 0; j < filter_order; j++) {
            if (i-j >= 0) {
                work[i] -= coeffs[j] * work[i-j];
            }
        }
    }

    /* Copy filtered signal to output buffer */
    for (i = 0; i < size; i++) {
        output[i] = work[i];
    }

    /* Clean up */
    free(work);
    free(poles_real);
    free(poles_imag);
    free(pz_real);
    free(pz_imag);
    free(coeffs);
}

int main()
{
    double input[] = {1.0, 2.0, 3.0, 4.0, 5.0};
    double output[5];
    int size = 5;
    double cutoff_freq = 2.0;
    double sampling_freq = 10.0;
    int filter_order = 2;

    butter_lowpass_filter(input, output, size, cutoff_freq, sampling_freq, filter_order);

    printf("Input: ");
    for (int i = 0; i < size; i++) {
        printf("%f ", input[i]);
    }
    printf("\n");

    printf("Output: ");
    for (int i = 0; i < size; i++) {
        printf("%f ", output[i]);
    }
    printf("\n");

    return 0;
}

这个程序实现了一个二阶 Butterworth 低通滤波器，可以通过修改 `cutoff_freq` 和 `filter_order` 参数来调整滤波器的截止频率和阶数。在此例中，输入信号为 `{1.0, 2.0, 3.0, 4.0, 5.0}`，采样频率为 10Hz，截止频率为 2Hz，输出信号为 `{1.000000, 1.596635, 2.375468, 3.305922, 4.354111}`。