python - 使用快速傅里叶变换生成相关随机势

Question

我想在具有指定自相关函数的 1D 或 2D 空间中生成随机势，并且根据一些数学推导，包括 Wiener-Khinchin 定理和傅里叶变换的性质，事实证明这可以使用以下等式完成： 其中phi(k)均匀分布在区间 [0, 1) 中。而这个功能满足，就是保证所产生的势永远是真实的。自相关函数应该不会影响我在这里做的事情，我采取了一个简单的高斯分布。

相位项和条件的选择phi(k)基于以下性质

1. 相位项的模数必须为 1（根据 Wiener-Khinchin 定理，即函数的自相关的傅里​​叶变换等于该函数的傅里叶变换的模数）；

2. 实函数的傅里叶变换必须满足（通过直接检查傅里叶变换的积分形式的定义）。

3. 产生的电位和自相关都是真实的。

通过结合这三个属性，该术语只能采用上述形式。

相关数学可参考以下pdf第16页： https ://d-nb.info/1007346671/34

我使用均匀分布随机生成了一个numpy数组，并将数组的负数与原始数组连接起来，使其满足上述条件phi(k)。然后我执行了 numpy（逆）快速傅里叶变换。

1D 和 2D 案例我都试过了，下面只显示 1D 案例。

import numpy as np
from numpy.fft import fft, ifft
import matplotlib.pyplot as plt

## The Gaussian autocorrelation function 
def c(x, V0, rho):
    return V0**2 * np.exp(-x**2/rho**2) 

x_min, x_max, interval_x = -10, 10, 10000
x = np.linspace(x_min, x_max, interval_x, endpoint=False)

V0 = 1
## the correlation length
rho = 1 

## (Uniformly) randomly generated array for k>0
phi1 = np.random.rand(int(interval_x)/2)
phi = np.concatenate((-1*phi1[::-1], phi1))
phase = np.exp(2j*np.pi*phi)

C = c(x, V0, rho) 
V = ifft(np.power(fft(C), 0.5)*phase)
plt.plot(x, V.real)
plt.plot(x, V.imag)
plt.show()

并且情节类似于如下所示 ：

然而，生成的势能变得很复杂，虚部与实部处于同一数量级，这是意料之中的。我已经检查了很多次数学，但我找不到任何问题。所以我在想它是否与实现问题有关，例如数据点是否足够密集以进行快速傅里叶变换等。

标签： pythonarraysnumpyfft

Answer 1

连接时需要更加小心：

phi1 = np.random.rand(int(interval_x)//2-1)
phi = np.concatenate(([0], phi1, [0], -phi1[::-1]))

第一个元素是偏移量（零频率模式）。“负”频率出现在中点之后。

这给了我