使用FFTW并评估生成的傅立叶级数

Using FFTW and evaluating the resulting Fourier series

本文关键字：傅立叶 FFTW 评估 使用      更新时间：2023-10-16

我是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;
  }

我得到的结果是假的：我只是在绘制系列时得到振荡波。

有人能给我一个线索，告诉我如何正确地评估得到的傅立叶级数吗？

您非常接近预期结果！获得正确输出的两点：

  • 第一项CCD_ 1是输入信号的平均值。所以是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+%gin",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提供专用于实数输入的函数，用于实数到复数的Output0和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);
相关文章：