Python/NumPy - 二维FFT无法得到解析解

编程语言 2026-07-10

我正在编写一个代码，作为第一步计算一个函数的二维FFT。我用函数exp(-r)/r进行测试，其中r = sqrt(x²+y²)，它的解析Hankel变换是1/sqrt(r²+1)。零阶Hankel变换与二维傅里叶变换之间通过一个重新缩放相关联。然而，我无法恢复这个解析解。我的变换比它更窄，且有非零的虚部。

最小可复现示例：

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))
    r0 = np.unravel_index(np.argmin(r, axis=None), r.shape)
    test = np.exp(-r)*np.reciprocal(r)
    test[r0] = 2*test[r0[0]+1, r0[0]] - test[r0[0]+2, r0[0]]  #this suppresses the r=0 singularity
    return test


x_min = np.linspace(-20, 20, num=1001)
y_min = np.linspace(-20, 20, num=1001)


real_a = min_pot(x_min, y_min)
shift_a = np.fft.ifftshift(real_a)
fourier_a = np.fft.fftshift(np.fft.fft2(np.fft.fftshift(real_a), norm="forward"))
kx = np.fft.fftshift(np.fft.fftfreq(len(x_min), d=(x_min[1]-x_min[0])))

#we represent a cut in the x axis

analytic = np.power(np.square(kx) + 1, -.5)

max_real = np.amax(fourier_a[fourier_a.shape[0]//2, :].real  # normalization to compare curve shapes better
fig, ax = pl.subplots(1, 1)
ax.plot(kx, fourier_a[fourier_a.shape[0]//2, :].real/max_real), label='Real fft')
ax.plot(kx, fourier_a[fourier_a.shape[0]//2, :].imag/max_imag), label='Imag. fft')
ax.plot(kx, analytic, label='Analytic')
pl.legend()
pl.show()

结果图是：数值结果与解析结果对比

有人知道我到底哪里做错了吗？

解决方案

这是一个不完整的回答，因为我无法准确看到你到底在做错什么，也不使用Python。凭直觉从你的图看，FFT的结果在宽度上似乎被放大了 1/2pi 倍，所以我在Excel做了一个快速探索性模型，以你的图为背景。解析函数普适化到 f(x,a) = a/sqrt(a^2+x^2)，因此它可以很容易生成一个归一化到 1 的图，且具有变化的全宽半最大值（FWHM）。

你的解析形式对于 a=1 来说是完全正确的。我的粗略模型在X 值较小时能很好拟合FFT结果，但当X 增大时发散，因为该函数理应渐近趋向于 a/x，但并非如此。