使用Python/Numpy计算2D FFT无法复现解析解问题求助

你的代码主要存在FFT移位逻辑错误、奇点处理不当、傅里叶变换与Hankel变换的缩放关系混淆以及语法错误这几个问题，以下是具体修复说明：

1. FFT移位逻辑错误 #

FFT默认要求信号原点位于数组左上角，正确的处理流程是：

• 先对中心在数组中间的信号执行ifftshift，将原点移到左上角；
• 执行fft2计算二维傅里叶变换；
• 最后对结果执行fftshift，将频谱原点移回数组中心。

你之前错误地用fftshift处理输入信号，破坏了FFT的输入对称性，直接导致了非零虚部和频谱变形。

2. r=0奇点处理不当 #

exp(-r)/r在r→0时呈发散特性（≈1/r），你用横向线性插值填充r=0处的值的方式破坏了径向对称性，进而引入虚部。正确做法是用极小值替代r=0，既避免除以零，又保留函数的发散特性：

r[r == 0] = 1e-10
test = np.exp(-r)/r

3. 傅里叶变换与Hankel变换的缩放关系 #

径向对称函数的2D傅里叶变换与零阶Hankel变换的关系为：
$$F(k) = 2\pi \cdot H_0f(r)$$
同时离散FFT与连续傅里叶变换之间需要考虑采样间隔的缩放因子，因此解析解需要调整为：

delta_x = x_min[1] - x_min[0]
analytic = 2 * np.pi * np.power(np.square(kx) + 1, -0.5) * delta_x * delta_x

4. 代码语法错误 #

• max_real定义时缺少右括号；
• max_imag未定义，需补充计算。

修复后的完整代码 #

import numpy as np
import matplotlib.pyplot as pl


def min_pot(x: np.ndarray, y: np.ndarray):
    xbig, ybig = np.meshgrid(x, y)
    r = np.sqrt(np.square(xbig) + np.square(ybig))
    # 处理r=0奇点
    r[r == 0] = 1e-10
    test = np.exp(-r) / r
    return test


x_min = np.linspace(-20, 20, num=1001)
y_min = np.linspace(-20, 20, num=1001)
delta_x = x_min[1] - x_min[0]

real_a = min_pot(x_min, y_min)
# 正确的FFT移位流程
fourier_a = np.fft.fftshift(np.fft.fft2(np.fft.ifftshift(real_a), norm="forward"))
kx = np.fft.fftshift(np.fft.fftfreq(len(x_min), d=delta_x))

# 修正后的解析解
analytic = 2 * np.pi * np.power(np.square(kx) + 1, -0.5) * delta_x * delta_x

# 提取x轴切片并归一化
fft_real_slice = fourier_a[fourier_a.shape[0]//2, :].real
fft_imag_slice = fourier_a[fourier_a.shape[0]//2, :].imag
max_real = np.amax(fft_real_slice)
max_imag = np.amax(np.abs(fft_imag_slice))

fig, ax = pl.subplots(1, 1)
ax.plot(kx, fft_real_slice / max_real, label='Real fft')
ax.plot(kx, fft_imag_slice / max_imag, label='Imag. fft')
ax.plot(kx, analytic / np.amax(analytic), label='Analytic')
pl.legend()
pl.show()

修复后，频谱虚部会降至数值误差级别，实部曲线形状将与解析解完全匹配。

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