时间戳

首页 > 解决方案 > numpy.fft 的有限精度低于某个值

python - numpy.fft 的有限精度低于某个值

问题描述

我正在尝试计算 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()

标签： pythonnumpyscipyfft

解决方案

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

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

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

推荐阅读