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

Question

我一直在做一个程序，它使用 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引入小的虚部？我真的很想避免这个问题。谁能给我一个建议？

Answer 1

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

我怀疑差异的根本原因与 MATLAB 的fft函数显然基于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 结果的虚部比较正好等于 0，而 FFTPACK 有少量的浮点舍入误差。

除此之外，我不知道为什么 FFTW 的实现比 FFTPACK 受到的舍入误差更少，但在任何情况下，重要的是要注意这些舍入误差足够小，它们通常不会引起问题（你知道你不应该测试浮点值之间的精确相等，对吗？）。

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

scipy.fftpack.fft(x).real

如果这些错误是一个问题，那么您可以切换到 usingpyfftw而不是 numpy/scipy，但是如果您的代码对舍入错误如此敏感，那么这可能意味着您无论如何都做错了。