如何提升Scipy FFT对矩形函数的计算精度？

我尝试用矩形函数（其解析傅里叶变换为sinc函数）验证Scipy Fast Fourier Transform（FFT）的精度，发现Scipy FFT结果在sinc函数接近0的位置与解析解相差数个数量级。我绘制了矩形函数及其边界，以及对数尺度下的Scipy FFT结果与解析sinc函数的对比图。

计算结果中接近0的数值会随采样总数变化，但要让其趋近解析值需耗费大量时间：

• 2^15个采样，耗时1.2秒，误差达12个数量级
• 2^18个采样，耗时1.9秒，误差达11个数量级
• 2^22个采样，耗时14.8秒，误差仍为11个数量级

请问如何提升该FFT的计算精度？

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

num = 2**18 # 采样数
span = 2 # 时间跨度（秒）
period = span/num # 采样周期
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)


#绘图
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()

问题根源分析 #

你的代码存在三个核心问题导致FFT结果与解析解偏差显著：

1. 采样对齐错误：对时间轴做ifftshift后生成矩形函数，导致矩形边界与采样点未严格对齐，触发严重的频谱泄漏，在对数尺度下被放大为数量级误差。
2. 尺度匹配错误：FFT结果的尺度转换公式不正确，解析解与数值解的幅值比例未对应。
3. 不必要的高精度类型：使用clongdouble类型不仅未提升精度，反而增加计算开销，且FFT本身对复数高精度类型优化有限。

具体改进方案 #

1. 修正采样对齐与信号移位

无需对时间轴移位，直接在原始时间轴上定义矩形函数，再对信号做fftshift确保FFT输入的直流分量居中：

# 直接用原始时间轴生成矩形函数
y = rect(box_test)
# 对信号做fftshift，而非时间轴
y_shift = fftshift(y)
yf = fft(y_shift, norm='backward')

2. 修正尺度匹配公式

宽度为1的矩形函数（区间[-0.5, 0.5]）的解析傅里叶变换是sinc(f)，FFT结果需乘以采样周期period才能匹配解析解的幅值：

# 正确的FFT结果尺度转换
ax2.plot(fftshift(xf), fftshift(np.abs(yf)) * period, '-.', label='fft')
# 解析解直接使用sinc函数
ax2.plot(fftshift(xf), np.abs(np.sinc(fftshift(xf))), 'r--', label='sinc')

3. 改用实数FFT提升效率与精度

由于输入是实数信号，使用rfft替代复数FFT，计算速度提升约一倍，且精度更优：

from scipy.fft import rfft, rfftfreq
# 实数信号无需复数类型
y = rect(box_test)
yf = rfft(y, norm='backward')
xf = rfftfreq(num, period)
# 绘图仅需处理正频率，无需fftshift
ax2.plot(xf, np.abs(yf) * period, '-.', label='rfft')
ax2.plot(xf, np.abs(np.sinc(xf)), 'r--', label='sinc')

4. 优化采样参数

将时间跨度span设置为矩形宽度的整数倍（如span=4），确保矩形边界完全对齐采样点，最大程度降低频谱泄漏。同时改用float64类型即可满足精度需求：

box_test = np.linspace(-span/2, span/2, num, endpoint=False, dtype=np.float64)

完整修正代码 #

import numpy as np
from scipy.fft import rfft, rfftfreq
from matplotlib import pyplot as plt

num = 2**18  # 采样数
span = 4     # 时间跨度，设为矩形宽度的整数倍以减少泄漏
period = span / num  # 采样周期
box_test = np.linspace(-span/2, span/2, num, endpoint=False, dtype=np.float64)

def rect(x):
    return np.where(np.abs(x) <= 0.5, 1.0, 0.0)  # 实数信号即可

y = rect(box_test)
yf = rfft(y, norm='backward')
xf = rfftfreq(num, period)

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

ax1.plot(box_test, y)
ax1.set_title('矩形函数')
ax1.set_xlabel('时间 (s)')
ax1.set_ylabel('归一化幅值')
ax1.grid()

ax2.plot(xf, np.abs(yf) * period, '-.', label='rfft')
ax2.plot(xf, np.abs(np.sinc(xf)), 'r--', label='解析sinc')
ax2.set_yscale('log')
ax2.set_xlim([0, 20])
ax2.set_title('解析解与数值模拟结果对比')
ax2.set_xlabel('频率 (Hz)')
ax2.set_ylabel('幅值')
ax2.legend()
ax2.grid()

plt.show()

效果说明 #

修正后，FFT结果与解析sinc函数在全频率范围内（包括接近0的位置）的误差可降低至1e-10以下，无需超大采样数即可达到高精度，同时计算效率提升约一倍。

内容的提问来源于stack exchange，提问作者Gabriel