matlab - 来自 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 的结果实际上包含非零虚部，而 MATLAB 的 fft 函数的结果却没有。

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

pyfftw但是，确实为 FFTW 提供了 Python 绑定(bind)。如果我们比较 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，但是如果您的代码那个敏感到舍入误差那么它可能意味着你做错了什么。