手动实现二维逆FFT积分时的偏移与幅值异常问题求助

问题背景 #

构建4个随时间依次平移1单位的高斯密度分布二维数组，执行FFT2后，通过ω-k域积分手动实现IFFT以还原原始密度，结果出现：

• 密度整体向右偏移
• 不同时间点的密度幅值异常
  尝试多种偏移修正后，仅高斯分布右侧出现变形，未解决根本问题。

原始代码 #

from os import path
import numpy as np
import matplotlib.pyplot as plt
from operator import add
from operator import sub
import gc
from pathlib import Path
from pylab import *
import cmath

x = []
indataxx = []

mmn = int(10)
dd = Path("""/home/maciek/test""")
outpath = str(dd/'test')


filename = dd/'nb.dat'
#filename = '/media/maciej/Nowy/trap_per_mu_0.6_Nb_300/gbf_' + str(i) + '/energy.dat'
indata = np.genfromtxt(filename ,usecols=(0))

mx = 300
nx = int(mx/2)
dx = 0.1

dl = len(indata)/(int(mx))
dl = int(dl)
print("ilosc zrzutow ", dl)


for i in range(0,mx):
    x.append(dx * (i-nx))

l2=[]

for j in range(0,dl):
    col = []
    for i in range (0,mx):
        col.append(indata[j*mx + i])
    l2.append(col)


l2 = np.array(l2)


# Apply 2D FFT
l2_shift = np.fft.ifftshift(l2)
l2_shift_fft = np.fft.fft2(l2_shift)
l2_shift_fft = np.fft.fftshift(l2_shift_fft)

magnitude_spectrum = np.abs(l2_shift)**2
print("Shape of magnitude_spectrum:", magnitude_spectrum.shape)


dt = mmn
Freqomega = np.linspace(-np.pi/(dt),np.pi/(dt),l2_shift_fft.shape[0],endpoint = False)
Freq_k_wec = np.linspace(-np.pi/(dx),np.pi/(dx),l2_shift_fft.shape[1],endpoint = False)

print("Shape of Freqomega:", Freqomega.shape)       #omegi
print("Shape of Freq_k_wec:", Freq_k_wec.shape)    #wektor k


S,D = np.meshgrid(Freq_k_wec, Freqomega)


with open(dd/'test_freqx.txt', "w") as f:
    np.savetxt(f, Freq_k_wec)

with open(dd/'test_freqy.txt', "w") as f:
    np.savetxt(f, Freqomega)

filename1 = dd/'test_freqy.txt'
indata1 = np.loadtxt(filename1, usecols=(0))

filename11 = dd/'test_freqx.txt'
indata11 = np.loadtxt(filename11, usecols=(0))


evo = []
evo1 = []
time = 4


print(" dk step ", Freq_k_wec[1]-Freq_k_wec[0])
print(" dw step ", Freqomega[1]-Freqomega[0])

ss = mx*dx/2
print(" przesuniecie po fft ", ss)


for yy in range(time):
    print("time " , yy)
    for u in range(len(indata1)):
        print("omega ", indata1[u])

        for j in range(mx):

            tmp = 0.

            for i in range(l2_shift_fft .shape[1]):
                tmp = tmp + l2_shift_fft[u , i] * np.exp( 1j * indata11[i] * x[j] )             #integral over k

            tmp = tmp * np.exp( 1j * indata1[u] * yy)
            tmp = tmp * np.abs(Freq_k_wec[1]-Freq_k_wec[0]) * np.abs(Freqomega[1]-Freqomega[0])  #dk*dw
            #tmp = tmp * np.exp( 1j * np.pi/dx ) * np.exp( 1j * np.pi/dt ) #shift
            evo.append(tmp)




for u in range(0,time,1):
    for i in range(mx):
        with open(dd/'xx_den.txt', "ab") as f:
            np.savetxt(f, np.c_[x[i],evo[u*mx + i].real/(l2_shift_fft.shape[0] * l2_shift_fft.shape[1])])

for yy in range(time):
    print(yy)
    for j in range(mx):
        atmp = 0.
        for u in range(len(indata1)):
            atmp = atmp + evo[yy*mx + u*(mx) + j]
        evo1.append(atmp.real/(l2_shift_fft.shape[0] * l2_shift_fft.shape[1]))



for u in range(time):
    for i in range(0,mx):
        with open(dd/'yy_test.txt', "ab") as f:
            np.savetxt(f, np.c_[x[i],evo1[u*mx + i]])

错误排查与修正方案 #

1. FFT移位的相位补偿缺失 #

• 错误原因：代码中对时域数据执行ifftshift后，又对频域结果执行fftshift，两次移位会给频域分量引入额外相位因子，但手动积分时未抵消该相位差，导致空间偏移。
• 修正：要么取消频域的fftshift，要么在积分时加入相位补偿项：

  # 在k域积分的exp项中加入补偿
  tmp += l2_shift_fft[u,i] * np.exp(1j * (indata11[i]*x[j] + np.pi*(i + u)))
  

2. 循环索引与积分逻辑混乱 #

• 错误原因：嵌套循环顺序错误（时间→频率→空间），导致evo数组的索引映射完全混乱，幅值叠加错误。
• 修正：调整循环顺序为空间→频率→时间，并使用数组存储结果避免索引错误：

  evo = np.zeros((time, mx), dtype=np.complex128)
  dk = Freq_k_wec[1] - Freq_k_wec[0]
  dw = Freqomega[1] - Freqomega[0]
  for yy in range(time):
      for j in range(mx):
          total = 0.0+0.0j
          for u in range(l2_shift_fft.shape[0]):
              omega = Freqomega[u]
              k_sum = 0.0+0.0j
              for i in range(l2_shift_fft.shape[1]):
                  k = Freq_k_wec[i]
                  k_sum += l2_shift_fft[u,i] * np.exp(1j * k * x[j])
              total += k_sum * np.exp(1j * omega * yy)
          evo[yy,j] = total * dk * dw
  

3. 归一化因子错误 #

• 错误原因：使用数组尺寸乘积作为归一化因子，不符合连续FFT的积分归一化规则（连续逆FFT需要除以$(2\pi)^2$）。
• 修正：将归一化因子改为$(2\pi)^2$：

  norm_factor = (2 * np.pi)**2
  # 后续幅值计算时使用
  evo[yy,j].real / norm_factor
  

4. 频率采样与FFT输出不匹配 #

• 错误原因：手动生成的频率数组与NumPyfft2的输出频率顺序不完全匹配，导致相位偏移。
• 修正：改用np.fft.fftfreq生成标准FFT频率，再移位：

  # 生成ω频率（时间轴）
  Freqomega = np.fft.fftfreq(l2.shape[0], d=dt) * 2 * np.pi
  Freqomega = np.fft.fftshift(Freqomega)
  # 生成k频率（空间轴）
  Freq_k_wec = np.fft.fftfreq(l2.shape[1], d=dx) * 2 * np.pi
  Freq_k_wec = np.fft.fftshift(Freq_k_wec)
  

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