MATLAB 1D FFT与C++ FFTW计算结果严重不匹配排查

AI 小北

2026-8-28

问题背景

我编写的MATLAB代码如下：

Nx = 10; 
Ny = 10; 
Lx = 2*pi; 
ygl = -cos(pi*(0:Ny)/Ny)'; %Gauss-Lobatto切比雪夫点 
x  = (0:Nx-1)/Nx*2*pi;
%构造网格
[X,Y]   = meshgrid(x,ygl); 
A = 2*pi / Lx;
u = sin( (2*pi / Lx) * X);
uh = fft(u)

上述代码运行输出结果如下：

0    6.4656   10.4616   10.4616    6.4656    0.0000   -6.4656  -10.4616  -10.4616   -6.4656
         0    0.0000         0         0    0.0000    0.0000   -0.0000         0         0         0
         0    0.0000         0         0    0.0000    0.0000         0         0         0   -0.0000
         0         0         0         0    0.0000    0.0000         0         0         0         0
         0    0.0000         0         0    0.0000    0.0000   -0.0000         0         0         0
         0    0.0000         0         0    0.0000    0.0000         0         0         0   -0.0000
         0    0.0000         0         0    0.0000    0.0000         0         0         0   -0.0000
         0    0.0000         0         0    0.0000    0.0000   -0.0000         0         0         0
         0         0         0         0    0.0000    0.0000         0         0         0         0
         0    0.0000         0         0    0.0000    0.0000         0         0         0   -0.0000
         0    0.0000         0         0    0.0000    0.0000   -0.0000         0         0         0

我同时编写了调用FFTW库的C++代码实现相同计算逻辑，代码如下：

static const int nx = 10;
static const int ny = 10; 
static const int nyk = ny/2 + 1;

double Lx = 2 * M_PI;
double A = (2 * M_PI)/Lx;

    double *XX;
    XX = (double*) fftw_malloc((nx*(ny+1))*sizeof(double));
    memset(XX, 42, (nx*(ny+1))* sizeof(double)); 
    

    double *YY;
    YY = (double*) fftw_malloc((nx*(ny+1))*sizeof(double));
    memset(YY, 42, (nx*(ny+1))* sizeof(double)); 

    double *u;
    u = (double*) fftw_malloc((((ny+1)*nx))*sizeof(double)); 
    

    fftw_complex *uh; 
    uh = (fftw_complex*) fftw_malloc(((ny+1)*nyk)*sizeof(fftw_complex)); 
    memset(uh, 42, ((ny+1)*nyk)* sizeof(fftw_complex)); 


    
    for(int i = 0; i< nx+1; i++){
        for(int j = 0; j< ny; j++){
            XX[i + (ny+1)*j] = (j)*2*M_PI/nx; 
            YY[i + (ny+1)*j] = -1. * cos(((i) * M_PI )/ny); 
            u[i + (ny+1)*j] = sin(A * XX[i + (ny+1)*j]);
        
        }   
    
    }
fftw_plan r2c1; 
r2c1 = fftw_plan_dft_r2c_1d(nx , &u[0], &uh[0], FFTW_ESTIMATE);
fftw_execute(r2c1);
fftw_destroy_plan(r2c1);
fftw_cleanup();
//打印uh结果

C++代码运行输出结果如下：

(0.0000,0.0000),  (0.0000,0.0000),  (0.0000,0.0000),  (0.0000,0.0000),  (0.0000,0.0000),  (0.0000,0.0000),
 (0.0000,0.0000),  (0.0000,0.0000),  (0.0000,0.0000),  (0.0000,0.0000),  (0.0000,0.0000),  (0.0000,0.0000),
 (0.0000,0.0000),  (0.0000,0.0000),  (0.0000,0.0000),  (0.0000,0.0000),  (0.0000,0.0000),  (0.0000,0.0000),
 (0.0000,0.0000),  (0.0000,0.0000),  (0.0000,0.0000),  (0.0000,0.0000),  (0.0000,0.0000),  (0.0000,0.0000),
 (0.0000,0.0000),  (0.0000,0.0000),  (0.0000,0.0000),  (0.0000,0.0000),  (0.0000,0.0000),  (0.0000,0.0000),
 (0.0000,0.0000),  (0.0000,0.0000),  (0.0000,0.0000),  (0.0000,0.0000),  (0.0000,0.0000),  (0.0000,0.0000),
 (0.0000,0.0000),  (0.0000,0.0000),  (0.0000,0.0000),  (0.0000,0.0000),  (0.0000,0.0000),  (0.0000,0.0000),
 (0.0000,0.0000),  (0.0000,0.0000),  (0.0000,0.0000),  (0.0000,0.0000),  (0.0000,0.0000),  (0.0000,0.0000),
 (0.0000,0.0000),  (0.0000,0.0000),  (0.0000,0.0000),  (0.0000,0.0000),  (0.0000,0.0000),  (0.0000,0.0000),
 (0.0000,0.0000),  (0.0000,0.0000),  (0.0000,0.0000),  (0.0000,0.0000),  (0.0000,0.0000),  (0.0000,0.0000),
 (0.0000,0.0000),  (0.0000,0.0000),  (0.0000,0.0000),  (0.0000,0.0000),  (0.0000,0.0000),  (0.0000,0.0000),

为何二者的计算结果完全不匹配？我知晓MATLAB与C的FFT实现存在部分差异，但本次为简单1D FFT计算，结果偏差远超合理范围，远大于常规实现差异。
注：应评论区要求，我补充了MATLAB与C的运行结果用于对比。

结果不匹配的核心原因

C++代码存在4个关键错误，和MATLAB计算逻辑完全不一致：

FFT计算维度完全错误：MATLAB中fft(u)未指定维度时，默认对矩阵第一个维度（列方向，即y方向，长度为Ny+1=11）做FFT。构造的u每一列（固定x坐标）的所有y点值都是相同的sin(x)，属于常数序列，因此FFT结果只有直流分量（第一行）有值，其余位置全为0，和贴出的MATLAB输出完全吻合。但C++代码创建的是长度nx=10的1D r2c plan，目标是做x方向FFT，和MATLAB实际计算的维度根本不是同一个。
内存布局错误，传入FFT的输入序列全为0：代码用索引i + (ny+1)*j访问数组，其中i是y方向索引、j是x方向索引，意味着x方向的10个点在内存中不是连续存储，相邻x点的索引间隔为ny+1=11。调用fftw_plan_dft_r2c_1d时传入&u[0]，只会读取从u[0]开始的连续10个元素，这些元素对应j=0（x=0）时i从0到9的y点值，sin(0)=0，输入序列全0自然输出全0。
仅执行1次FFT，未做批量计算：MATLAB对11列y方向序列全部做了FFT，共11次1D变换，但C++代码只执行了1次FFT，就算输入正确也只能得到1行结果，和MATLAB的11行输出维度不匹配。
memset用法错误：memset是按单字节赋值的，写memset(XX, 42, ...)会把double类型数组的每个字节都设为0x2A，得到完全无意义的垃圾值。虽然循环覆盖了u数组的所有元素，但uh数组除了第一次FFT写入的前6个元素，其余位置全是垃圾值；如果循环边界写错，部分u元素也会是垃圾值。

修正方案

先明确FFT计算方向：如果要复现MATLAB当前结果，需要做沿y方向（长度11）的批量实数FFT；如果要做x方向FFT，MATLAB代码需要改成fft(u, [], 2)指定沿第二维度计算。
调整数组内存布局，把待做FFT的维度放在快维，即同一组做FFT的点在内存中连续存放，比如做x方向FFT时，用索引u[y*nx + x]存储，每个y对应的nx个x点连续存放。
批量FFT推荐用fftw_plan_many_dft_r2c接口，指定变换次数、输入输出步长和距离，无需手动循环调用单个plan。
浮点数组初始化不要用memset赋非0值，double数组请用循环赋值，仅需要填充0时可以用memset(..., 0, ...)。

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