C++ fftw3实现matlabz=filter(b,a,y);

在C++中使用fftw3库实现MATLAB中的filter函数，需要进行以下步骤：

1. 安装fftw3库，可以从官网下载对应平台的库文件；
2. 将需要进行滤波的信号进行FFT变换，得到频域数据；
3. 根据滤波器系数b和a，生成对应的滤波器；
4. 将滤波器也进行FFT变换，得到频域滤波器；
5. 将频域数据和频域滤波器进行点乘，得到滤波后的频域数据；
6. 对滤波后的频域数据进行反FFT变换，得到时域滤波后的信号。

下面是一个简单的C++代码示例，实现了和MATLAB中filter函数相似的功能：

#include <fftw3.h>
#include <iostream>
#include <vector>

using namespace std;

vector<double> filter(vector<double> x, vector<double> b, vector<double> a)
{
    int N = x.size();    // 信号长度
    double *in, *out;    // 输入和输出缓冲区
    fftw_complex *out_freq;    // 输出频域缓冲区
    fftw_plan p, q;    // FFT变换计划

    // 分配内存空间
    in = (double*) fftw_malloc(sizeof(double) * N);
    out = (double*) fftw_malloc(sizeof(double) * N);
    out_freq = (fftw_complex*) fftw_malloc(sizeof(fftw_complex) * (N / 2 + 1));

    // 初始化输入信号
    for (int i = 0; i < N; i++) {
        in[i] = x[i];
    }

    // 创建FFT变换计划
    p = fftw_plan_r2r_1d(N, in, out, FFTW_R2HC, FFTW_ESTIMATE);
    q = fftw_plan_dft_r2c_1d(N, in, out_freq, FFTW_ESTIMATE);

    // 执行FFT变换
    fftw_execute(p);
    fftw_execute(q);

    // 生成滤波器
    int nb = b.size();
    int na = a.size();
    vector<double> filter_b(nb, 0);
    vector<double> filter_a(na, 0);
    for (int i = 0; i < nb; i++) {
        filter_b[i] = b[i] / a[0];
    }
    for (int i = 0; i < na; i++) {
        filter_a[i] = a[i] / a[0];
    }

    // 将滤波器转化为频域滤波器
    fftw_complex *filter_b_freq, *filter_a_freq;
    filter_b_freq = (fftw_complex*) fftw_malloc(sizeof(fftw_complex) * (N / 2 + 1));
    filter_a_freq = (fftw_complex*) fftw_malloc(sizeof(fftw_complex) * (N / 2 + 1));
    for (int i = 0; i < N / 2 + 1; i++) {
        filter_b_freq[i][0] = filter_b[i];
        filter_b_freq[i][1] = 0;
        filter_a_freq[i][0] = filter_a[i];
        filter_a_freq[i][1] = 0;
    }
    fftw_plan fb = fftw_plan_dft_1d(N, filter_b_freq, filter_b_freq, FFTW_FORWARD, FFTW_ESTIMATE);
    fftw_plan fa = fftw_plan_dft_1d(N, filter_a_freq, filter_a_freq, FFTW_FORWARD, FFTW_ESTIMATE);
    fftw_execute(fb);
    fftw_execute(fa);

    // 将频域数据和频域滤波器进行点乘
    for (int i = 0; i < N / 2 + 1; i++) {
        double re_x = out_freq[i][0];
        double im_x = out_freq[i][1];
        double re_b = filter_b_freq[i][0];
        double im_b = filter_b_freq[i][1];
        double re_a = filter_a_freq[i][0];
        double im_a = filter_a_freq[i][1];
        out_freq[i][0] = (re_x * re_b + im_x * im_b) / N;
        out_freq[i][1] = (im_x * re_b - re_x * im_b) / N;
        if (i > 0) {
            out_freq[i][0] /= re_a;
            out_freq[i][1] /= re_a;
        }
    }

    // 执行IFFT变换，得到时域信号
    fftw_plan r = fftw_plan_dft_c2r_1d(N, out_freq, out, FFTW_ESTIMATE);
    fftw_execute(r);

    // 输出结果
    vector<double> y(N);
    for (int i = 0; i < N; i++) {
        y[i] = out[i];
    }

    // 释放内存空间
    fftw_destroy_plan(p);
    fftw_destroy_plan(q);
    fftw_destroy_plan(fb);
    fftw_destroy_plan(fa);
    fftw_destroy_plan(r);
    fftw_free(in);
    fftw_free(out);
    fftw_free(out_freq);
    fftw_free(filter_b_freq);
    fftw_free(filter_a_freq);

    return y;
}

这里的filter函数接收三个参数，分别是输入信号x、滤波器系数b和a，返回滤波后的输出信号y。在函数中，先将输入信号进行FFT变换得到频域数据，然后根据滤波器系数生成频域滤波器，将频域数据和频域滤波器进行点乘，最后对点乘后的频域数据进行IFFT变换得到时域信号。注意，在计算频域滤波器时，需要将滤波器系数b和a进行归一化，以避免在频域计算时产生误差。