具有高精度数字的 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，您会得到以下意想不到的结果...实部完全错误。 ..振幅偏离了，它不对称，没有任何意义。

我不知道是什么原因造成的。有人有什么想法吗？

Answer 1

这完全是意料之中的行为。计算机使用浮点计算，本质上是不精确的。

注意您的真实结果的 y 轴。如果不存在数值误差，则实数部分将完全为 0。根据您的 "higher precision" 结果，实数部分几乎与 0 相同（1e-14 非常接近双精度浮点数的精度）。精度较低时，实部变得更大（尽管仍然比虚部小得多）。由于数字越大，结构也越多（即误差不是由舍入误差给出的，而是由输入数据的实际特征给出的，周期略短于理想值）。

Answer 2

正如@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()

具有高精度数字的 numpy.fft.fft 的意外行为

Unexpected behavior of numpy.fft.fft with high precision numbers

python

numpy

fft