Scipy FFT 精度比箱函数的分析 FT 差得多

Question

我正在尝试使用具有已定义解析表达式（sinc 函数）的框函数来检查 scipy 快速傅里叶变换函数的准确性。 scipy FFT 的结果给出的值与 sinc 函数结果接近 0 时相差几个数量级。

在我下面的图表中，我绘制了框函数和传递给它的边界以创建框，在第二个图中，我以对数刻度显示了 scipy fft 和分析等价的 sinc 函数。

“接近 0”的计算值确实会随着样本总数的变化而变化，但我很快需要大量时间才能将这些值推低至分析值。 例如：

• 2^15个样本，1.2s，相差12个数量级
• 2^18个样本，1.9s，11个数量级的差异
• 2^22个样本，14.8s，相差11个数量级

如何提高此 FFT 的准确性？

import numpy as np
from scipy.fft import fft, fftshift, fftfreq
from matlabplot import pyplot as plt

num = 2**18 # sample spacing
span = 2 # seconds
period = span/num #1.0 / 1e9
box_test = np.linspace(-period*num/2, period*num/2, num, endpoint= False, dtype='clongdouble')
shifted_box = ifftshift(box_test)
def rect(x):
    return np.where(np.abs(x)<=0.5, 1 + 0j, 0 + 0j)


y = rect(box_test)
y_shift = rect(shifted_box)
yf = fft(y_shift, norm='backward')
xf = fftfreq(num, period)


#plotting
fig, (ax1, ax2) = plt.subplots(2, 1, figsize=(20,10))

plt.grid()
ax1.plot(box_test, y)
ax2.plot(fftshift(xf),  fftshift(np.absolute(yf))*span/num, '-.', label='fft')
ax2.plot(fftshift(xf), fftshift(np.absolute(np.sinc(xf))), 'r--', label = 'sinc')
ax1.set_title('box function')
ax1.set_xlabel('time (s)')
ax1.set_ylabel('arb units')
ax2.set_yscale('log')
ax2.set_xlim([-20,20])
# ax2.set_ylim([-10,10])
ax2.set_title("analytical and numerical simulation Optical homodyne PSD")
ax2.set_xlabel("Frequency (Hz)")
ax2.set_ylabel("theoretical detector output (arb units)")
plt.legend()
plt.grid()
plt.show()

Answer 1

FFT 或 DFT 对 FFT 长度的窗函数进行循环卷积（不是线性卷积）。循环卷积将数据包裹在 FFT 周围。由于 Sinc 函数具有无限的非零“尾部”，因此它将围绕 FFT 结果缠绕其尾部，附加地，几乎无限次（达到数值噪声的极限）。较长的 FFT 意味着环绕的尾部部分离中心较远，因此幅度较低。

如果你想要一个“类似盒子”的函数，它在 FFT 后产生有限宽度的加窗 Sinc，使用 FFT 库的最佳方法是采用所需加窗 Sinc 的逆 FFT，加窗到长度等于或短于 FFT.