Matlab fft 和 Scipy fft 的 FFT 结果略有不同

Posted 2023-03-12

【中文标题】Matlab fft 和 Scipy fft 的 FFT 结果略有不同

【英文标题】：Slightly different FFT results from Matlab fft and Scipy fft

【发布时间】：2015-11-16 15:54:11

【问题描述】：

我一直在制作一个使用 NumPy/Scipy 测量两个光谱之间相位差的例程。

我已经有了用Matlab编写的例程，所以我基本上用NumPy重新实现了函数和相应的单元测试。但是，我发现单元测试失败了，因为scipy.fftpack.fft 引入了一些小的数值错误：

import numpy as np
import scipy.fftpack.fft
x = np.array([0.0, 1.0, 2.0, 3.0, 4.0, 3.0, 2.0, 1.0])
X = scipy.fftpack.fft(x)

在这种情况下，由于时域信号是对称的，所以期望的输出是

[16.0000   -6.8284         0   -1.1716         0   -1.1716         0   -6.8284]

如下Matlab代码所示：

>> x = [0.0, 1.0, 2.0, 3.0, 4.0, 3.0, 2.0, 1.0];
>> X = fft(x)

X =

   16.0000   -6.8284         0   -1.1716         0   -1.1716         0   -6.8284

根据 DSP 理论，结果不应包含 任何 虚构成分。但是，scipy的结果如下：

array([ 16.00000000 +0.00000000e+00j,  -6.82842712 -2.22044605e-16j,
         0.00000000 -0.00000000e+00j,  -1.17157288 -2.22044605e-16j,
         0.00000000 +0.00000000e+00j,  -1.17157288 +2.22044605e-16j,
         0.00000000 +0.00000000e+00j,  -6.82842712 +2.22044605e-16j])

为什么scipy.fftpack.fft 会引入小的虚构组件？我真的很想避免这个问题。谁能给我一个建议？

【问题讨论】：

在处理有限精度计算时，您应该始终准备好预期一定程度的数值误差。在这种情况下，它比信号电平低 16 个数量级，非常合理。

我明白，但是从DSP理论来看，这种情况下，FFT结果应该是真实的。那么这是否意味着 Scipy FFT 库比 Matlab FFT 库差？使用 Matlab 库，FFT 结果是实数，不包含任何复数。

【参考方案1】：

一方面，保证scipy.fftpack.fft 总是返回一个复杂的结果，而 MATLAB 的 fft 函数的结果有时是真实的，有时是复杂的，这取决于是否存在非零虚部。但是，这并不能解释为什么 scipy.fftpack.fft 的结果实际上包含非零虚部，而 MATLAB 的 fft 函数的结果却没有。

我怀疑造成这种差异的根本原因与 MATLAB 的 fft 函数显然是 based 上的 FFTW，而 scipy 和 numpy 由于许可限制使用 FFTPACK 这一事实有关。

pyfftw 但是，确实提供了与 FFTW 的 Python 绑定。如果我们比较 FFTPACK 和 FFTW 结果的虚部：

from pyfftw.interfaces import scipy_fftpack as fftw

Fx1 = fftpack.fft(x)
print(Fx1.imag)
# [  0.00000000e+00  -2.22044605e-16  -0.00000000e+00  -2.22044605e-16
#    0.00000000e+00   2.22044605e-16   0.00000000e+00   2.22044605e-16]
print(Fx1.imag == 0)
# [ True False  True False  True False  True False]

Fx2 = fftw.fft(x)
print(Fx2.imag)
# [ 0.  0.  0.  0.  0.  0.  0.  0.]
print(Fx2.imag == 0)
# [ True  True  True  True  True  True  True  True]

我们看到 FFTW 结果的虚部比较完全等于零，而 FFTPACK 有少量浮点舍入误差。

除此之外，我不知道为什么 FFTW 的实现比 FFTPACK 受舍入误差的影响要小，但无论如何，重要的是要注意这些舍入误差足够小，通常不会引起问题（你知道你应该'不是要测试浮点值之间的精确相等性，对吗？）。

通常你会简单地取结果的真实部分，例如：

scipy.fftpack.fft(x).real

如果这些错误是问题，那么您可以切换到使用pyfftw 而不是 numpy/scipy，但如果您的代码 对舍入错误很敏感，那么它可能意味着你做错了什么。

【讨论】：

您可能理解这一点，但要明确一点：scipy.fftpack.rfft 返回一个“非复杂结果”，但这只是因为它将复杂组件表示为真实数组的不同元素。 scipy.fftpack.rfft([0, 1, 2, 3, 4, 3, 2, 1]) 仍然给出包含非零（但很小）虚部的结果。

通常在处理浮点数时，您应该使用误差范围进行测试。大约几个粒度步骤就足够了。 en.wikibooks.org/wiki/Floating_Point/Epsilon