我正在尝试使用 C 中的 FFTW 从 Python 中的一些已知工作代码重新创建结果。我发现结果中有一些小错误。
scipy.fft
我的输入数据是真实的 3d,尺寸 = (294, 294, 294)。
我的 scipy.fft 调用如下所示:
我的 fftw 代码如下所示这个:
complex_data_out = scipy.fft.fftn(real_data_in,s=dimensions)
注意
fftw_plan plan = fftw_plan_dft(3, dimensions, complex_data_in, complex_data_out, FFTW_FORWARD, FFTW_ESTIMATE);
fftw_execute(plan);
所有虚值都为0。
complex_data_in == real_data_in
我已经验证了
两个程序之间的匹配。
real_data_in
无论如何,我得到的结果略有不同。如果我从两个complex_data_out数组中抓取几个任意点,我会得到这些值:
有关缩放因子的一些调试:
SciPy:
pt 1: 187904.8906250 + 132361.1718750 i
pt 2: 112964.4062500 + 237751.6250000 i
FTTW:
pt 1: 187709.8070491 + 132244.8899144 i
pt 2: 113043.6615987 + 238079.0344479 i
这些差异非常小,以至于它们可能由某些归一化因子引起是没有意义的。而且尺度并不一致。
scipy val one r / fftw val one r = 1.001039283
scipy val one i / fftw val one i = 1.000879293
scipy val two r / fftw val two r = 0.999298896
scipy val two r / fftw val two r = 0.998624787
是否有人有过使用 FFTW 与 SciPy、MATLAB 或 NumPy 等得到稍微不同的结果的经验?我看到人们在互联网上讨论很大的差异(由于标准化因素),但不是小的不一致差异。
我的第一个预感是,也许我的 C 代码在数据分配或数据初始化的某个地方出现了一个错误,但是我打猎了又打猎,但我不认为就是这样。
My first hunch is that maybe my C code has an off by one error somewhere in data assignment or data initialization but I've hunted and hunted and I don't think that is it.
这些微小的差异很可能是由浮点数精度和 C 和 Python 实现的 FFTW 所使用的特定算法之间的差异造成的。以下是一些需要考虑的因素:
-
浮点数精度: 浮点数运算本身就不精确,并且执行操作的顺序会影响结果。由于 SciPy 和 FFTW 使用不同的库(可能是不同优化的库)来执行浮点数运算,因此它们可能会在结果中引入微小的差异。
-
算法差异: FFTW 使用各种算法来计算 FFT,而 SciPy 使用的具体实现可能会使用与的 C 代码略有不同的算法或算法变体。即使是算法中的细微差异也会导致不同的舍入误差和最终结果的细微差异。
-
数据类型: 请确保 C 代码和 Python 代码中使用的数据类型一致。例如,确保使用的是
fftw_complex
或double[2]
来表示 C 代码中的复数,并且的 Python 代码使用相应的数据类型(通常是numpy.complex128
)。 -
FFTW 规划标志: 正在使用
FFTW_ESTIMATE
规划标志,这可能会导致 FFTW 选择并非所有情况都最精确的算法。尝试使用FFTW_MEASURE
或FFTW_PATIENT
来查看是否有助于提高精度,但请注意,这些标志会导致计划时间更长。 -
数据初始化: 提到已经检查了数据,但请仔细检查是否已在 C 代码和 Python 代码中以完全相同的方式初始化了输入数据。即使是很小的差异也可能传播并导致结果不同。
调试技巧:
-
减少问题范围: 尝试使用更小的数据集,看看差异是否仍然存在。如果差异出现在较大的数据集上,则可能表明数值误差累积。
-
比较中间结果: 如果可能,请比较 C 代码和 Python 代码中的算法的中间步骤。这可以帮助查明差异出现的位置。
-
使用其他库: 尝试使用另一个 Python FFT 库(如 NumPy 的
numpy.fft
)来查看是否获得了与 SciPy 相似的结果。
总的来说,看到的结果差异很小,并且可能在数值计算中预料到的范围内。但是,如果需要更高的精度,则可能需要研究使用不同的 FFT 库、探索不同的算法或使用更高精度的浮点数据类型。
标签:python,c,scipy,fft,fftw From: 78795322