c++ - 使用FFTW并评估所得的傅里叶级数

我是FFTW的新手。我想将一个函数分解为一个傅立叶级数。到目前为止，我还没有做到。我的代码如下:

  // 1) Create discretizations for my function 'my_function'
  int N = 100; // number of discretizations
  float x_step = (x_end - x_start) / ((float)(N - 1));
  fftw_complex *Input, *Output;
  fftw_plan my_plan;
  Input = (fftw_complex*)fftw_malloc(sizeof(fftw_complex) * N);
  Output = (fftw_complex*)fftw_malloc(sizeof(fftw_complex) * N);
  float x = x_start;
  ForIndex(i, N) {
    Input[i][0] = my_function(x);
    cout << "Input[" << i << "]=" << Input[i][0] << endl;
    Input[i][1] = 0;
    x += x_step;
  }

  my_plan = fftw_plan_dft_1d(N, Input, Output, FFTW_FORWARD, FFTW_ESTIMATE);

  fftw_execute(my_plan);

  // 3) Evaluation - this is the part I am confused with
  // I should get something very close to 'my_function' when I plot, shouldn't I?
  ForIndex(i, N) {
    float sum = 0;
    float x = (float)i;
    sum = Output[0][0] / 2.0f;
    for (int k = 1; k < N; k++) {
      // Fourier series
      // http://en.wikipedia.org/wiki/Fourier_series
      float s = 2.0f*(float)M_PI * (float)k * x / (float)N;
      sum += Output[k][0] * std::cos(s) + Output[k][1] * std::sin(s);
      // I also tried
      // sum += Output[k][0] * std::sin(s + Output[k][1]);
      // to no avail
    }
    cout << "For x=" << x << ", y=" << (sum / (float)N) << endl;
  }

我得到的结果是虚假的:在绘制序列时，我只会得到振荡波。

有人可以给我一个线索，以了解如何正确评估所得的傅立叶级数吗？

最佳答案

您非常接近预期的结果!获得正确的输出有两点:

第一项Output是输入信号的平均值。所以它是sum = Output[0][0];，而不是sum = Output[0][0] / 2.0f;

i是i*i==-1的复数，因此，当虚部相乘时，必须在结果中添加减号。因此:

sum += Output[k][0] * std::cos(s) - Output[k][1] * std::sin(s);

这是由g++ main.cpp -o main -lfftw3 -lm编译的更正的代码:

#include <fftw3.h>
#include <iostream>
#include <cmath>
#include <stdlib.h>

using namespace std;

float my_function(float x){
    return x;
}
int main(int argc, char* argv[]){

    // 1) Create discretizations for my function 'my_function'
    int N = 100; // number of discretizations
    float x_end=1;
    float x_start=0;
    int i;
    float x_step = (x_end - x_start) / ((float)(N - 1));
    fftw_complex *Input, *Output;
    fftw_plan my_plan;
    Input = (fftw_complex*)fftw_malloc(sizeof(fftw_complex) * N);
    Output = (fftw_complex*)fftw_malloc(sizeof(fftw_complex) * N);
    float x = x_start;
    for(i=0;i<N;i++){
        Input[i][0] = my_function(x);
        cout << "Input[" << i << "]=" << Input[i][0] << endl;
        Input[i][1] = 0;
        x += x_step;
    }

    my_plan = fftw_plan_dft_1d(N, Input, Output, FFTW_FORWARD, FFTW_ESTIMATE);

    fftw_execute(my_plan);

    for(i=0;i<N;i++){
      printf("%d %g+%gi\n",i,Output[i][0],Output[i][1]);
    }

    // 3) Evaluation - this is the part I am confused with
    // I should get something very close to 'my_function' when I plot, shouldn't I?
    for(i=0;i<N;i++) {
        float sum = 0;
        float x = (float)i;
        sum = Output[0][0];
        for (int k = 1; k < N; k++) {
            // Fourier series
            // http://en.wikipedia.org/wiki/Fourier_series
            float s = 2.0f*(float)M_PI * (float)k * x / (float)N;
            sum += Output[k][0] * std::cos(s) - Output[k][1] * std::sin(s);
            // I also tried
            // sum += Output[k][0] * std::sin(s + Output[k][1]);
            // to no avail
        }
        cout << "For i=" << i <<" x="<<x_start+x_step*i<< " f(x)="<<my_function(x_start+x_step*i)<<" y=" << (sum / (float)N) << endl;
    }

}

我还添加了一些东西来打印傅里叶变换的系数。由于输入是真实信号，因此频率k和N-k的系数是共轭的。由于N=100，请查看频率49和51或48和52。

因此，库fftw提供functions dedicated to real inputs，fftw_plan fftw_plan_dft_r2c_1d()(从实到复数)和fftw_plan fftw_plan_dft_c2r_1d()。存储并计算一半的频率(实际上是(N+1)/2)

要回到现实世界，您可以使用 FFTW_BACKWARD 使用第二个计划:

my_plan2 = fftw_plan_dft_1d(N, Output, input, FFTW_BACKWARD, FFTW_ESTIMATE);

关于c++ - 使用FFTW并评估所得的傅里叶级数，我们在Stack Overflow上找到一个类似的问题：https://stackoverflow.com/questions/28630810/