使用Python中的FFT和IFFT以偶数点采样时，实函数自由传播的离散化函数变得复杂

Question

我的代码涉及到使用傅里叶变换和逆傅里叶变换来传播一个真实的函数。具体来说，这个函数的演变可以用以下公式表示：

∂ψ(z,t)/∂t - v ∂ψ(z,t)/∂z =0

我通过对上面的方程进行傅里叶变换来解决这个问题，变换后的方程是：

∂ψ(k,t)/∂t + ikv ψ(k,t)=0

这个方程的解是：

ψ(k,t)=e^(-ikvt)ψ(k,0)

现在对这个结果进行逆傅里叶变换，就能得到传播后的函数：

ψ(z-vt,0)

然而，当我在Python中执行这个过程（代码见下文）时，最终得到的函数（经过快速傅里叶变换和逆变换后）似乎有一个微小的复数部分，特别是在函数采样为偶数个点时。这个误差在我的模拟过程中逐渐累积，导致结果不够准确。不过，当我将函数采样为奇数个点时，这个问题似乎就消失了。有人能帮我解释一下这是怎么回事吗？

下面是一个简单的代码示例，展示了我的观点。运行这段代码后，我们可以看到，当函数采样为10个点时，经过快速傅里叶变换和逆变换后，函数的实部的范数有时只能准确到小数点后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))))

从这里可以看到，在偶数采样的情况下（只考虑实部），范数在第二位小数上有所不同，如果这只是浮点误差，我会感到非常惊讶。

Answer 1

看看 k_vals：

当 Nd=6 时：

[0., 6.28318531, 12.56637061, -18.84955592, -12.56637061, -6.28318531]

当 Nd=5 时：

[0., 6.28318531, 12.56637061, -12.56637061, -6.28318531]

当 Nd 是偶数时，有一个负频率没有对应的正频率。我们知道，输入到逆傅里叶变换（IFFT）的数据必须是共轭对称的，这样输出才会是实数值。对于奇数长度的数组，计算出的指数函数会自动是共轭对称的。但对于偶数长度的数组，有一个频率没有正的对应项。我们需要让这个频率成分变成实数值，这样指数函数才能是共轭对称的，逆傅里叶变换的结果也才能是实数值。

例如，下面这个对你函数的修改就能实现这个目的：

def free_prop(vectr, Nd, vel=0.92, dt=1):
    dz=1/Nd
    psi_k = np.fft.fft(vectr)
    k_vals = 2 * np.pi * np.fft.fftfreq(Nd, d=dz)
    k_kernel = np.exp(-1j * dt * k_vals * vel)
    if np.mod(Nd, 2) == 0:
       k_kernel[Nd // 2] = k_kernel[Nd // 2].real
    return np.fft.ifft(k_kernel * psi_k)

至于这是否仍然符合物理规律，我就不太清楚了。也有人可能会说那个频率区间应该是0，因为根据定义它包含了混叠的数据（这是奈奎斯特频率，为了正确采样信号，在奈奎斯特频率及以上的频率必须是0能量）。