numpy.fft 的有限精度低于特定值

Question

我正在尝试计算 Python 中某些信号的傅里叶变换。我希望通过快速傅里叶变换计算的结果与根据定义计算的结果一致。然而，使用numpy.fft计算的结果与预期值存在偏差。

信号没有达到低于特定数字的值。在下图中，它大约是 10^-16。对于其他信号，这些是可比较的值（从 10^-9 到 10^-30）。在我的应用程序中，我需要更高的精度。

为了确保我也测试了 scipy.fftpack。那里出现了同样的错误，尽管错误计算的值略有不同。该问题与信号参数（长度、采样频率等）无关。

这个限制的原因是什么？如果它是 Python/Numpy 准确度，我该如何改进它？

# Fourier Transform

import numpy as np
import scipy.fftpack as fp

def gaussian_distribution(x,mean=0,variance=1):
    return (1 / (np.sqrt(2*np.pi)*variance) ) * np.exp( -((x-mean)**2) / (2 * variance**2) )

def gaussian_fourier_transform(omega, mean=0, variance=1):
    # http://mathworld.wolfram.com/FourierTransformGaussian.html
    return np.exp(-2 * np.pi**2 * variance**2 * omega**2) * np.exp(-(2j)*np.pi*mean*omega)

## signal generation
signal_range = [-2**4, 2**4]
N = 2**8
x = np.linspace(signal_range[0],signal_range[1],N, dtype='float64')
y = gaussian_distribution(x)

## calculating result
framerate = N / (signal_range[1] - signal_range[0])
frequency_axis = np.linspace(0,framerate,N)

numpy_v = np.abs( np.fft.fft(y) )
numpy_v = numpy_v / numpy_v[0] # normalization

scipy_v = np.abs( fp.fft(y) )
scipy_v = scipy_v / scipy_v[0]

symbolical_v = gaussian_fourier_transform(frequency_axis)

# ploting

import matplotlib.lines as mlines
import matplotlib.pyplot as plt

fig = plt.figure()
ax1 = fig.add_subplot()

ax1.plot(frequency_axis[0: N//2], scipy_v[0: N//2], '.r')
ax1.plot(frequency_axis[0: N//2], numpy_v[0: N//2], '.b')
ax1.plot(frequency_axis[0: N//2], symbolical_v[0: N//2], 'g')
ax1.set_yscale('log')
ax1.grid(True)

blue_line = mlines.Line2D([], [], color='blue', marker='.', markersize=15, label='result calculated by numpy.fft')
red_line = mlines.Line2D([], [], color='red', marker='.', markersize=15, label='result calculated by scipy.fftpack')
green_line = mlines.Line2D([], [], color='green', marker='', markersize=15, label='result calculated by definition')
ax1.legend(handles=[blue_line, red_line, green_line])

fig.show()

Answer 1

这是一个numerical issue。都是NumPy and SciPy use different variants of PocketFFT which is based upon FFTPack，属于exact FFT的范畴，其误差取决于机器误差ε和样本数N。我不确定这些库的确切依赖性是什么，但您可以尝试 pyFFTW which are Python bindings to FFTW and those might have 对机器错误的依赖性稍微好一些。

Answer 2

IEEE 双精度浮点数（您的计算机 CPU 可能支持的硬件）具有大约 15 位十进制数字的精度。这是因为只有 53 位尾数（或尾数）。 FFT 算法通过 O(N*Log(N)) 增加此误差范围（或量化噪声），其中 N 是 FFT 长度。

因此，为了获得更高的精度（更低的本底噪声），您可能必须找到或编写自己的内部使用四精度或任意精度算术的 FFT，并以该格式获取和输入数据.

例如，您可以尝试使用 python 的 mpmath package 对 FFT 进行编码，然后选择精度。

Answer 3

（和别人说的差不多）函数的离散傅里叶变换的每个点都依赖于初始函数的每个点，这意味着相对误差与输入数据的最大幅度值成比例并且在某种意义上是全局的。最常见的和您观察到的具有长尾和小尾巴的傅立叶变换在尾部有很大的误差。 如果你的精度是固定的并且你不能进行（部分）分析计算，则没有办法解决这个问题，但你可以使用任意精度库来提高你的精度。

但是如果进行卷积是您的最终目标，您可能会使用快速多极方法实现高达机器精度的局部相对误差。这是否以及如何工作取决于您的内核函数。

numpy.fft 的有限精度低于特定值

Limitated accuracy of numpy.fft below a certain value

python

numpy

fft

scipy