numpy.fft.fft 函数在使用高精度数时表现异常

Question

我有以下代码...请注意# generate sine curve下的两行代码。其中一行使用比另一行更高精度的2pi值，尽管它们应该仍然给出接近相同的结果。

import numpy as np
import matplotlib.pyplot as plt


t1 = np.arange(0., 1., .01)

# generate sine curve
y1 = np.sin(6.28318*5.*t1)  
#y1 = np.sin(6.283185307179586*5.*t1) # equivalent to np.sin(2*np.pi*t1)

# calculate the fft (no averaging!) of the time series
ffty = np.fft.fft(y1)

fig, ax_list = plt.subplots(3,1)
ax_list[0].plot(t1,y1, '.-')

ax_list[1].plot(ffty.real, '.-', label='Real Part')
ax_list[1].legend()

ax_list[2].plot(ffty.imag, '.-', label='Imag Part')
ax_list[2].legend()


plt.show()

如果您使用较低的精度6.28318运行代码，则可以获得fft的预期结果... 然而，如果您使用更高的精度6.283185307179586（等于2.*numpy.pi）运行代码，则会得到以下意外结果...实部明显错误...幅值偏差很大，不对称，没有任何意义。 我不知道是什么原因导致了这种情况。有人有什么想法吗？

- J Webster

---

5e-14大约等于零，而且几乎所有的浮点数都是近似值。 - hotpaw2

Answer 1

如@Cris Luengo所说，您需要查看y轴的刻度以准确比较两个图形。另一种方法是将您要比较的两个内容都绘制在同一张图上，就像我下面做的那样。

使用对数刻度显示FFT的幅度，并且很明显，使用更少的pi有效数字确实会导致更低的精度结果。 当使用浮点数时，大多数值并不完全为零，这是可以预期的，但是使用更多有效数字将获得许多数量级的改进，这在分别绘制FFT时不会立即显现。 使用的代码：

import numpy as np
import matplotlib.pyplot as plt


t1 = np.arange(0., 1., .01)

values = {
'low':6.28318,
'higher':6.283185307179586,
'highest':2*numpy.pi,
}
styles = {
'low':'-',
'higher':'-',
'highest':'.-'
}

fig, ax_list = plt.subplots(3,1)
for name, tau in values.items():
    y1 = np.sin(tau*5.*t1)
    ffty = np.fft.fft(y1)

    ax_list[0].plot(t1,y1, styles[name], label=name)
    ax_list[1].plot(abs(ffty.real), styles[name],label=name)
    ax_list[2].plot(abs(ffty.imag), styles[name], label=name)

[ax.legend() for ax in ax_list]
ax_list[0].set_title('time domain')
ax_list[1].set_title('real part')
ax_list[2].set_title('imaginary part')
ax_list[1].set_yscale('log')
ax_list[2].set_yscale('log')
plt.draw()

- user2699

Answer 2

这是完全预期的行为。计算机使用浮点运算，这本身就是不精确的。

请注意实际结果的y轴标尺。如果没有数值误差存在，实部将恒等于0。使用“更高的精度”结果时，实部几乎与0相同（1e-14非常接近双精度浮点数的精度）。当精度较低时，实部变得更大（尽管仍然比虚部小得多）。由于数字更大，也有更多结构（即误差不是由四舍五入误差导致，而是由输入数据的实际特征引起的，即一个略短于理想值的周期）。

- Cris Luengo

---

你的意思是说，如果我在numpy.sin中使用numpy.pi来生成正弦波，然后将其放入numpy.fft.fft中...并且得到了实部结果的无意义输出...那么这是预期行为吗？不是。 - J Webster

@JWebster：这不是胡言乱语，而是舍入误差。你有注意到那个y轴顶部的“1e-14”吗？如果你真的期望双精度浮点计算是无误差的，你需要多读一些资料。你可以从这里开始：《计算机科学家应该了解的浮点运算知识》（David Goldberg著） - Cris Luengo

@JWebster：好的，如果你认为那就是你想要的，那就这样吧。祝你找到一个能够给你不同结果的FFT工具。像我建议的那样设置轴限制，然后再告诉我实部与零点有何不同。 - Cris Luengo

哦...也许这对我来说是一个可教的时刻。这是否意味着错误实际上是在“低精度”实部中存在5 Hz的峰值？而正确的结果是“高精度”情况，其中所有值几乎都是零？我可以理解...我想我只是期望看到5 Hz的峰值，所以当我没有看到它们时...另外，我希望结果在正负频率上镜像，但是，这只是浮点数问题。 - J Webster

@JWebster：是的，正弦波具有纯虚频谱。请注意，当虚部较大时，噪声方差也会更大--复数的幅度更大，因此舍入误差也更大。这就是为什么第二种情况中实部有一个峰值的原因。在第一种情况下，由于您的周期略小于理想值，信号不完全是奇函数，因此变换不是完全虚数，它具有与空间域的小偶分量相对应的小实分量。 - Cris Luengo