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

Question

我想在具有指定自相关函数的 1D 或 2D 空间中生成随机势，并且根据一些数学推导，包括 Wiener-Khinchin 定理和傅立叶变换的性质，结果证明这可以使用以下等式完成:

哪里phi(k)在区间 [0, 1) 内均匀分布。而这个函数满足 ，这是为了确保产生的势总是真实的。
自相关函数应该不会影响我在这里做的事情，我取了一个简单的高斯分布 .
phi(k)相项的选择和条件基于以下属性

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

实函数的傅立叶变换必须满足 (通过直接检查积分形式的傅立叶变换的定义)。

生成的电位和自相关都是真实的。

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

相关的数学可以引用以下pdf的p.16:
https://d-nb.info/1007346671/34

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

一维和二维情况我都试过，下面只显示一维情况。

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()

该图类似于如下所示:
.

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

Answer 1

连接时需要更加小心:

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

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

这给了我