Python/Numpy中2D FFT无法还原exp(-r)/r函数解析解的技术求助

阿华AIGC实验室

2026-4-27

问题分析与修复方案

我帮你找出了代码里的几个关键问题，这些问题直接导致了FFT结果和解析解不匹配，以及非零虚部的出现：

1. FFT移位顺序错误（核心问题）

你当前的FFT处理顺序完全搞反了：

fourier_a = np.fft.fftshift(np.fft.fft2(np.fft.fftshift(real_a), norm="forward"))

正确的2D FFT流程应该是：

先用np.fft.ifftshift把图像的原点（中心）移到数组左上角（FFT算法默认原点在数组角落）
执行fft2计算二维傅里叶变换
再用np.fft.fftshift把变换后的频谱原点移回中心

错误的移位顺序会让频谱的位置彻底错乱，这是结果形状和解析解差异巨大的主要原因。

2. r=0奇点的处理方式不合理

你用线性外推设置r=0处的值，但实际上当r→0时，exp(-r)/r趋近于无穷大。离散网格中r=0是单个点，直接用极小值替代r=0更合理，既能避免除以零，又能减少数值误差：

r = np.sqrt(xbig**2 + ybig**2)
# 用极小值替代r=0，避免除以零
epsilon = 1e-6
r = np.where(r < epsilon, epsilon, r)
test = np.exp(-r) / r

3. 汉克尔变换与2D FFT的缩放关系未考虑

零阶汉克尔变换和2D FFT之间存在明确的缩放因子：
对于径向对称函数f(r)，其2D傅里叶变换F(k)满足：
F(k) = 2π * H₀[f](k)
其中H₀[f](k)是你提到的解析解1/√(k²+1)，所以解析解需要乘以2π才能和FFT结果匹配。

4. 语法错误与归一化问题

你的代码里有两处语法错误：

max_real的定义缺少右括号
max_imag未定义就直接使用

另外，归一化时应该让解析解和FFT结果的最大值对齐，这样才能更直观地对比曲线形状。

修复后的完整代码

import numpy as np
import matplotlib.pyplot as plt

def min_pot(x: np.ndarray, y: np.ndarray):
    xbig, ybig = np.meshgrid(x, y)
    r = np.sqrt(xbig**2 + ybig**2)
    # 处理r=0的奇点
    epsilon = 1e-6
    r = np.where(r < epsilon, epsilon, r)
    test = np.exp(-r) / r
    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)

# 正确的FFT流程
shift_a = np.fft.ifftshift(real_a)
fourier_a = np.fft.fftshift(np.fft.fft2(shift_a, norm="forward"))

# 计算波数kx
kx = np.fft.fftshift(np.fft.fftfreq(len(x_min), d=(x_min[1]-x_min[0])))

# 解析解：加入2π缩放因子
analytic = 2 * np.pi * np.power(np.square(kx) + 1, -0.5)

# 提取中心行的FFT结果
fft_real = fourier_a[fourier_a.shape[0]//2, :].real
fft_imag = fourier_a[fourier_a.shape[0]//2, :].imag

# 归一化到最大值（方便对比形状）
max_fft_real = np.amax(fft_real)
norm_fft_real = fft_real / max_fft_real
norm_analytic = analytic / np.amax(analytic)

# 绘图
fig, ax = plt.subplots(1, 1)
ax.plot(kx, norm_fft_real, label='Real FFT (normalized)')
ax.plot(kx, fft_imag, label='Imag. FFT')
ax.plot(kx, norm_analytic, label='Analytic (normalized)')
ax.legend()
plt.xlim(-5, 5)  # 限制x轴范围，聚焦核心特征
plt.show()