c++ 写一个滤波程序

A：

滤波程序一般由两个部分构成：滤波器设计和滤波器应用。下面以实现一种简单的FIR(有限脉冲响应)低通滤波器为例，展示如何编写一个滤波程序。

1. 滤波器设计

考虑设计一个截止频率为$f_c$的低通FIR滤波器。采用窗函数法，可得到其传递函数为：

$$
H(z)=\frac{1}{N}\cdot\sum_{n=0}^{N-1}h[n]z^{-n}=\frac{1}{N}\cdot\sum_{n=0}^{N-1}h[n]\cdot e^{-j2\pi fn}
$$

其中$h[n]$为滤波器的冲激响应，$N$为滤波器的阶数，$f=\frac{f_c}{f_s}$为归一化截止频率，$f_s$为采样频率。这里给出一个简单的FIR滤波器的代码实现：

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

#define PI 3.1415926

//定义滤波器的阶数和截止频率
#define N 31
#define fc 0.4

//FIR低通滤波器设计
void fir_design(double h[], int n, double fc)
{
    int i, m;
    double wn, sum = 0;

    for(i = 0; i < n; i++) {
        m = i - n / 2;
        wn = 2 * fc * (m==0 ? 1 : sin(2 * PI * fc * m) / (2 * PI * fc * m));
        h[i] = wn * (0.54 - 0.46 * cos(2 * PI * i / (n - 1))); // Hamming 窗
        sum += h[i];
    }

    for(i = 0; i < n; i++) {
        h[i] /= sum;
        printf("h[%d] = %f\n", i, h[i]);
    }
}

int main()
{
    int i, j;
    double h[N];
    double x[N], y[N];
    double y_tmp;

    fir_design(h, N, fc);

    //发现了一个错误（原文是y[i] = 0）
    for(i = 0; i < N; i++) {
        x[i] = y[i] = 0;
    }

    //模拟输入信号，这里取三角波
    for(i = 0; i < 1000; i++) {
        x[i % N] = i / N % 2 == 0 ? i % N : N - i % N;
        y_tmp = 0;
        for(j = 0; j < N; j++) {
            y_tmp += h[j] * x[(i-j+N)%N];
        }
        y[i % N] = y_tmp;
        printf("%f\n", y_tmp);
    }

    return 0;
}

2. 滤波器应用

将上面的滤波器设计程序中的$h$作为已知量，可以通过卷积实现对输入信号的滤波。将滤波器应用程序集成到主程序中，即可完成滤波程序的编写，并使用所需的输入。在本例中，我们对一个三角波进行FIR低通滤波操作：将三角波信号通过设计好的滤波器滤波以后，输出的信号即为经过滤波器处理后的信号。注意到这里是采用的 overlap-add 的形式。