FFTW实现正反FFT后数组首元素出现NaN值的问题排查与解决

咱们一步步拆解你的问题——这根本不是数据类型导致的，而是代码里的几个关键逻辑错误直接引发了NaN值和结果异常。下面我来逐个指出问题并给出修复方案：

1. 反向FFT函数的核心逻辑颠倒 #

你的FFT2complex函数犯了一个低级但致命的顺序错误：

• 你先创建了反向变换计划，但此时in是刚用fftwl_malloc分配的未初始化内存，里面全是垃圾值；
• 接着你直接把这些垃圾值赋值给了U；
• 最后才执行fftwl_execute(plan_b)——这时候变换出来的正确结果根本没被用到！

除此之外，你还搞反了FFT计划的输入输出：反向变换应该把频域数据（也就是U里的内容）作为输入，变换后输出到时域。但你的代码错误地把U的数据放到了out变量，而计划是从in到out的变换，完全错位了。

2. 错误重复调用fftwl_cleanup() #

你在complex2FFT和FFT2complex两个函数里都调用了fftwl_cleanup()，这个函数会清除FFTW的所有内部状态（包括已分配的Wisdom缓存等）。如果在后续还要调用FFTW操作（比如你先跑正向再跑反向），这会导致后续的FFT出现不可预料的错误，包括内存问题或计算异常。正确的做法是只在程序结束前调用一次。

3. 归一化因子设置不合理 #

FFTW的正向和反向变换本身不带归一化，行业通用的做法是：

• 正向变换后除以N（这里N=Ypt）；
• 反向变换不做额外归一化；
  或者反过来。你在complex2FFT里用了2.0/Ypt，这会导致结果幅度异常，虽然不是NaN的直接原因，但会影响最终结果的正确性。

修正后的完整代码 #

下面是修复了所有问题的版本，我标注了关键修改点：

#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= 1.0/Ypt; // 修改为标准的1/N归一化
    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);
    
    // 将U的时域数据拷贝到FFT输入in
    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);
    
    // 将频域结果赋值回U并应用归一化
    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);
    
    // 将U的频域数据拷贝到反向FFT输入in
    for (int i = 0; i < Ypt; ++i){
        in[i][0] = creal(U[i]);
        in[i][1] = cimag(U[i]);
    }
    
    // 创建反向变换计划：in(频域) → out(时域)
    plan_b = fftwl_plan_dft_1d(Ypt, in, out, FFTW_BACKWARD, FFTW_ESTIMATE);
    fftwl_execute(plan_b); // 先执行变换，再读取结果
    
    // 将时域结果赋值回U
    for (int i = 0; i < Ypt; ++i){
        U[i] = normalizing_factor * out[i][0] + normalizing_factor * out[i][1] * I;
    }
    
    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){
        // 避免不必要的double强制转换，直接用long double计算保持精度
        long double x = 2.0 * PI * i * dy;
        U[i] = sin(x);
        V[i] = sin(x);
    }
    
    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 %Lf %Lf \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 %Lf %Lf \n", (long double)i*dy, creal(U[i]), cimag(U[i]) );
    }
    
    fclose(stream); // 别忘了关闭文件，确保内容写入磁盘
    free(U);
    free(V);
    fftwl_cleanup(); // 只在程序结束时调用一次
}

关键修改点总结 #

• 修正反向FFT逻辑顺序：先拷贝输入数据，执行变换，再读取结果赋值回U；
• 移除重复的fftwl_cleanup()，仅在程序末尾调用一次；
• 调整归一化因子：采用FFT标准的归一化方式，保证结果幅度正确；
• 优化数值计算：避免不必要的类型转换，全程用long double计算，保持精度；
• 添加文件关闭操作：修复原代码未关闭文件的潜在问题。

运行修正后的代码，你会发现变换后的结果和原始正弦值几乎完全一致，不会再出现NaN值，第一个元素也能正常显示。

内容的提问来源于stack exchange，提问作者Aziz