Python中FFT/IFFT传播实函数时偶数采样的虚部误差问题

问题背景 #

通过傅里叶变换（FFT）与逆傅里叶变换（IFFT）求解实函数的自由传播，对应对流方程：
$$\frac{\partial \psi(z,t)}{\partial t} - v \frac{\partial \psi(z,t)}{\partial z} = 0$$
傅里叶变换后频域方程为：
$$\frac{\partial \psi(k,t)}{\partial t} + ikv \psi(k,t) = 0$$
其解为：
$$\psi(k,t) = e^{-ikvt}\psi(k,0)$$
逆傅里叶变换后得到传播后的实函数$\psi(z-vt,0)$。

现象复现 #

在Python实现时发现：

• 偶数采样点数下，FFT/IFFT处理后结果会出现明显微小虚部，仅取实部时范数偏差较大（仅2-3位小数精度），但复数整体范数守恒
• 奇数采样点数下，虚部误差为浮点级别，实部范数精度可达16位小数

复现代码如下：

import numpy as np

def free_prop(vectr,Nd,vel=0.92,dt=1):
    dz=1/Nd
    psi_k = np.fft.fft(vectr)
    k_vals = 2.0*np.pi*np.fft.fftfreq(Nd, d=dz)
    return np.fft.ifft(np.exp(-1j*dt*k_vals*vel)*psi_k)

vectr_even=np.random.rand(10) ## Even case
vectr_odd=np.random.rand(11) ## Odd case

print('Norm of input array (with even sampling):',np.linalg.norm(vectr_even))
print('Norm of output complex array (with even sampling):',np.linalg.norm(free_prop(vectr_even,10)))
print('Norm of real part of output array (with even sampling):',np.linalg.norm(np.real(free_prop(vectr_even,10))))
print()
print('Norm of input array (with odd sampling):', np.linalg.norm(vectr_odd))
print('Norm of output complex array (with odd sampling):',np.linalg.norm(free_prop(vectr_odd,11)))
print('Norm of real part of complex array (with odd sampling):',np.linalg.norm(np.real(free_prop(vectr_odd,11))))

原因分析 #

实函数的FFT满足共轭对称性：$\psi(k) = \psi^*(-k)$，即正频率分量与负频率分量互为共轭，保证IFFT结果为实数。

• 奇数采样：所有频率分量$k$都存在唯一的$-k \neq k$（在FFT的频率索引映射中），乘以相位因子$e^{{-ikvt}$后，$\psi(k)e}{-ikvt}$的共轭为$\psi^*(-k)e{ikvt} = \psi(-k)e^{{ikvt}$，而$-k$对应的相位因子正好是$e}{-i(-k)vt}=e^{ikvt}$，因此共轭对称性得以保持，IFFT结果为纯实数（仅浮点误差级虚部）。
• 偶数采样：存在一个特殊的频率分量——$k = \pi/dz$（对应FFT索引中的$N/2$位置），该分量的$-k$与自身在FFT频率映射中等价（即$-k \equiv k \mod 2\pi/dz$），因此原FFT结果中该分量必须为实数。但乘以相位因子$e^{-ikvt}$后，该分量变为复数（当$kvt$不是$\pi$的整数倍时），直接破坏了共轭对称性，导致IFFT结果出现非浮点级别的虚部，取实部时丢失部分能量，范数偏差明显。

解决方法 #

方法1：使用实值FFT（rfft/irfft） #

利用NumPy的np.fft.rfft和np.fft.irfft，这两个函数专门处理实输入信号：

• rfft仅计算正频率及零频率分量（一半的采样点数），自动维护共轭对称性
• irfft根据对称性质重构实值输出，避免手动维护对称性的麻烦

修改后的代码：

import numpy as np

def free_prop_real(vectr, Nd, vel=0.92, dt=1):
    dz = 1/Nd
    psi_k = np.fft.rfft(vectr)
    k_vals = 2.0*np.pi*np.fft.rfftfreq(Nd, d=dz)
    return np.fft.irfft(np.exp(-1j*dt*k_vals*vel)*psi_k)

vectr_even=np.random.rand(10) ## Even case
vectr_odd=np.random.rand(11) ## Odd case

print('Norm of input array (with even sampling):',np.linalg.norm(vectr_even))
print('Norm of output array (with even sampling):',np.linalg.norm(free_prop_real(vectr_even,10)))
print('Max imaginary part (even sampling):', np.max(np.abs(np.imag(free_prop_real(vectr_even,10)))))
print()
print('Norm of input array (with odd sampling):', np.linalg.norm(vectr_odd))
print('Norm of output array (with odd sampling):',np.linalg.norm(free_prop_real(vectr_odd,11)))
print('Max imaginary part (odd sampling):', np.max(np.abs(np.imag(free_prop_real(vectr_odd,11)))))

运行后，偶数采样下的虚部误差会降至浮点级别（约$10^{-16}$），实部范数精度与奇数采样一致。

方法2：手动维护频域共轭对称性 #

如果必须使用普通FFT/IFFT，处理后需手动修正频域分量的对称性：

1. 对偶数采样的中间分量（索引$N/2$），将其强制为实数（取实部或乘以共轭相位）
2. 确保$\psi(k)$与$\psi(N-k)$互为共轭（$k=1$到$N/2-1$）

这种方法需要额外的代码逻辑，不如rfft/irfft简洁高效。

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