fft2 (matlab) 和 fftw (C) 的不同结果

Posted 2023-02-16

【中文标题】fft2 (matlab) 和 fftw (C) 的不同结果

【英文标题】：Different results with fft2 (matlab) and fftw (C)

【发布时间】：2014-06-30 08:31:37

【问题描述】：

我正在尝试使用 FFTW3 库在 C 中实现 Matlab fft2() 函数。

但是，我得到了不同的结果。

考虑下一个矩阵：

Z=[ 
    0.4791    0.4765    0.4791    0.4765    0.4791    0.4765    0.4791    0.4765     
    0.4798    0.4695    0.4798    0.4695    0.4798    0.4695    0.4798    0.4695     
    0.4791    0.4765    0.4791    0.4765    0.4791    0.4765    0.4791    0.4765     
    0.4798    0.4695    0.4798    0.4695    0.4798    0.4695    0.4798    0.4695     
    0.4791    0.4765    0.4791    0.4765    0.4791    0.4765    0.4791    0.4765     
    0.4798    0.4695    0.4798    0.4695    0.4798    0.4695    0.4798    0.4695     
    0.4791    0.4765    0.4791    0.4765    0.4791    0.4765    0.4791    0.4765     
    0.4798    0.4695    0.4798    0.4695    0.4798    0.4695    0.4798    0.4695     
    .... 
] 

并使用以下代码：

Z-> Double* 

fftw_complex* fft2; 
fft2 = (fftw_complex *)fftw_malloc(sizeof(fftw_complex)*Samples*(Lines)); 

fftw_plan p1; 

p1 = fftw_plan_dft_r2c_2d(Lines,Samples, Z, fft2, FFTW_ESTIMATE); 
fftw_execute(p1); 

Matlab 的结果：

fft2= [ 
 5534,25859596829 + 0,00000000000000i     186,747610745237 - 529,515274347496i
 42,6452471730436 - 321,074636721419i    -21,4495750160608 - 190,407528614266i
-50,3875107145668 - 50,5480303619799i     30,1151029075525 + 378,240946095017i
-196,295569635431 + 228,972218925794i     35,6434356803659 - 5,46216875816971i
 36,2702126322693 - 38,5502177293316i     18,5093049539101 - 33,4608602804025i
     .... 
 ]

我的 C 代码的结果：

5534.260423 + 0.000000 i           186.731496 + -529.495788 i 
  42.655319 + -321.068356 i        -21.425010 + -190.382717 i 
 -50.277195 + -50.384210 i          29.909846 + 377.823957 i 
 -195.767224 + 228.693862 i         35.241375 + -5.315382 i 
 36.134134 + -38.527643 i           18.406395 + -33.467351 i 
    .... 
] 

我做错了什么？

【问题讨论】：

你有错误的期望。计算机计算通常不准确。

@JoachimPileborg 否。MATLAB 引擎是用 Fortran 和/或 C 编写的，使用良好的老式 IEEE754 浮点类型。

@DavidHeffernan MATLAB 使用双精度，而不是浮点数

@mch 我没有另外说。还是说 IEEE754 双精度类型不是浮点类型？

@undur_gongor 你是对的，但 OP 的差异足以引起 IMO 的关注。如果使用相同的类型和相似的算法，那么结果应该只在机器精度的数量级上偏离。您使用的 FFTW 版本可能是单精度的（我很确定 MATLAB 版本是双精度的）。但是根据我的经验，即使 MATLAB 声称使用了 FFTW，我认为他们已经以某种方式更改了代码（可能更好地优化了 codelet）。我注意到速度明显不同here。

【参考方案1】：

您在 C 实现中没有做错任何事情，但您无法确定您是否正在比较同类。

您可以在 C 中的 FFTW 和 MATLAB 中的 fft2() 之间获得不同结果的原因多种多样。

MATLAB 使用自己专门优化的 FFTW 版本，编译时支持多线程和 SSE/AVX 矢量化指令。您可以在 MATLAB 中使用version('-fftw') 找到它的详细信息。 MATLAB 在 FFTW 之上使用它拥有的抽象 (libmwmfl_fft) 来取消规划器例程并公开简单的函数，例如 fft() 和 fft2()。您不能确定它是否会选择与您的问题相同的计划程序例程。 FFTW 使用启发式方法来确定针对您指定的数据大小和类型计算 FFT 的最佳算法。这些启发式方法可能因在相同数据和同一台机器上的运行而异，尤其是FFTW_ESTIMATE，因为它使用的测试程序比FFTW_MEASURE、FFTW_PATIENT 或FFTW_EXHAUSTIVE 更不严格。 浮点运算不一定是关联的，并且可能导致计算不同的值，具体取决于 FFTW 决定用于计算 FFT 的算法。根据所选择的算法，浮点舍入误差可能会在某些计算中传播，从而导致准确性降低。

---

我查看了您在问题中发布的数据的相对误差，对于某些值，它是 1e^-5，而平均相对误差是 0.2341。

总体而言，尽管考虑到浮点计算的性质，但谁能说哪个计算 正确 值？ C 或 MATLAB 中的 FFTW？

【参考方案2】：

在比较 Matlab 的 fft2 与 FFTW 之前，您应该验证在给定相同输入的情况下，Matlab 和 FFTW 执行之间的可重现结果。因为对于相同的输入，默认的 FFTW 不会每次都产生相同的结果。我知道这很令人不安，（并且可能会在这个答案上引起很多负面评价。）差异很小。

我自己在使用 FFTW 时遇到了这个问题。我使用输出的 MD5 校验和来验证没有区别。通过更改默认的Planner Flags for FFTW，我能够重现相等的校验和（相同的输出）。此问题已在in Q3.8 of the FFTW FAQ 讨论。

确定性 FFTW 工作后，对 Matlab 执行相同操作（因为它只是包装 FFTW，并提供修改 FFTW 计划的选项）。

然后比较 FFTW 和 Matlab 的结果。

【讨论】：

搜索关键字“FFTW”+“deterministic”会在这个主题上找到很多结果。