FFT 返回 NaN 值答案

【问题标题】：FFT returns NaN valuesFFT 返回 NaN 值
【发布时间】：2021-09-18 16:56:46
【问题描述】：

我正在尝试将 FFTW 用于我的大型项目，因此我编写了一个基本程序来检查 FFT 是否正常工作。我正在尝试将正弦值发送到 FFT 并返回。在使用 FFT 前向和 FFT 后向我得到完全相同的结果，但我在数组的第一个元素处得到一个 NaN 值。我在这里读到一些问题，问题是数据类型，我需要使用long double 和fftwl 来保持结果的准确性。哪里有问题？是否在数据类型中，我该如何解决？

#include <fftw3.h>
#include <math.h>
#include <stdio.h>
#include <complex.h> 
#include <stdlib.h>
#include <inttypes.h>
#include <assert.h>

int Ypt=128;
long double PI=3.14159265358979323846;

void complex2FFT( complex long double *U)
{
  long double normalizing_factor= 2.0/Ypt;
  fftwl_plan plan_f;
  fftwl_complex *in;
  fftwl_complex *out;

  in  = (fftwl_complex *) fftwl_malloc(sizeof(fftwl_complex) * Ypt);
  out = (fftwl_complex *) fftwl_malloc(sizeof(fftwl_complex) * Ypt);
  
  for (int i = 0; i < Ypt; ++i){
    in[i ][0]= creal(U[i]);
    in[i ][1]= cimag(U[i]);
  }

  plan_f = fftwl_plan_dft_1d(Ypt, in, out, FFTW_FORWARD , FFTW_ESTIMATE);

  fftwl_execute(plan_f);

  for (int i = 0; i < Ypt; ++i){
    U[i] = normalizing_factor*out[i][0] + normalizing_factor*out[i][1]*I;
  }

  fftwl_destroy_plan(plan_f);
  fftwl_free(in);
  fftwl_free(out);
  fftwl_cleanup();
}

void FFT2complex( complex long double *U)
{
  long double normalizing_factor= 1.0;
  fftwl_plan plan_b;
  fftwl_complex *in;
  fftwl_complex *out;

  in  = (fftwl_complex *) fftwl_malloc(sizeof(fftwl_complex) * Ypt);
  out = (fftwl_complex *) fftwl_malloc(sizeof(fftwl_complex) * Ypt);
  
  for (int i = 0; i < Ypt; ++i){
    out[i ][0]= creal(U[i]);
    out[i ][1]= cimag(U[i]);
  }

  plan_b = fftwl_plan_dft_1d(Ypt, in, out, FFTW_BACKWARD, FFTW_ESTIMATE);
  for (int i = 0; i < Ypt; ++i){
    U[i] = normalizing_factor*in[i ][0] + normalizing_factor*in[i ][1]*I;
  }

  fftwl_execute(plan_b);

  fftwl_destroy_plan(plan_b);
  fftwl_free(in);
  fftwl_free(out);
  fftwl_cleanup();
}

int main(int argc, char **argv){
  long double dy=( (long double)1)  / ( (long double)Ypt);
  complex long double *U = malloc(Ypt * sizeof(*U)); 
  complex long double *V = malloc(Ypt * sizeof(*V)); 

  for (int i = 0; i < Ypt; ++i){
    U[i] = sin( (double) (2.0*PI* (double)i * dy)) ;
    V[i] = sin( (double) (2.0*PI* (double)i * dy)) ;
  }

  char name[45];
  FILE *stream;
  sprintf(name, "V%d.txt", 0);
  stream= fopen(name,"w");

  for (int i = 0; i < Ypt; ++i){
    fprintf(stream, "%Lf %2.5f %2.5f \n", (long double)i*dy,creal(U[i]),cimag(U[i]));
  }
       
  complex2FFT(U);

  FFT2complex(U);

  for (int i = 0; i < Ypt; ++i){
    fprintf(stream, "%Lf %2.5f %2.5f \n", (long double)i*dy,creal(U[i]),cimag(U[i]) );
  }

  free(U);
  free(V);
}

【问题讨论】：

fftwl_execute(plan_b) 在您复制输出后运行。
我改变了它，仍然得到 NaN 值。问题是，如果我使用 Ypt=16，我不会得到任何 nan 值，但对于 Ypt=32,128，我会得到 nan 值
嗯，U[i] = sin( ...) ; V[i] = sin( ...) ;} --> 我希望其中一个使用cos()。
V[] 已分配、设置和释放，但从未使用过。为什么会出现在代码中？
我使用 V[] 来比较 fft_forward 和 fft_backward 前后 Sin(U) 的值。因此，我将正弦的值签名为 U 和 V，然后 fft U 并比较结果

标签： c double fftw

【解决方案1】：

@Dietrich Epp 注意到， fftwl_execute(plan_b) 在复制输出后运行，因此输出保持不变。此外，FFTW_BACKWARD 执行反向 DFT，但并不意味着将参数 in 用作输出。即使有标志FFTW_BACKWARD，fftw_plan_dft_1d() 的输入参数是in，输出是out。例如，参见前向然后后向变换的示例：FFTW forward and back ward yield in different results why?

for (int i = 0; i < Ypt; ++i){
  out[i ][0]= creal(U[i]);
  out[i ][1]= cimag(U[i]);
}

plan_b = fftwl_plan_dft_1d(Ypt, out, in, FFTW_BACKWARD, FFTW_ESTIMATE);

fftwl_execute(plan_b);
for (int i = 0; i < Ypt; ++i){
  U[i] = normalizing_factor*in[i ][0] + normalizing_factor*in[i ][1]*I;
}

【讨论】：