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

Question

我想用指定的自相关函数在一维或二维空间中生成一个随机势，根据一些数学推导，包括维纳-辛钦定理和傅里叶变换的性质，结果证明这可以是使用以下等式完成：其中 phi(k) 均匀分布在区间 [0, 1) 中。而这个函数满足，就是保证产生的势总是真实的。自相关函数应该不会影响我在这里做的事情，我取一个简单的高斯分布。

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

相位项的模必须为1（根据维纳-欣钦定理，即函数自相关的傅里叶变换等于该函数傅里叶变换的模）；
一个实函数的傅里叶变换必须满足（直接看傅里叶变换的积分形式的定义）。
生成的势和自相关都是真实的。

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

相关数学可以参考以下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]))

第一个元素是偏移量（零频率模式）。 "Negative" 频率在中点之后。

这给了我

Answer 2

您对fft（更准确地说，DFT）的运作方式有一些误解。首先请注意，DFT 假设序列的样本索引为 0, 1, ..., N-1，其中 N 是样本数。相反，您生成一个对应于索引 -10000, ..., 10000 的序列。其次，请注意，实数序列的 DFT 将为对应于 0 和 N/2 的 "frequencies" 生成实数值。您似乎也没有考虑到这一点。

我不会详细介绍，因为这超出了这个 stackexchange 站点的范围。

只是为了完整性检查，下面的代码生成了一个序列，该序列具有实值序列的 DFT (FFT) 预期的属性：

正负频率共轭对称，
对应频率0和N/2
假定序列对应于索引 0 到 N-1

可以看到，这个序列的ifft确实生成了一个实数序列

from scipy.fftpack import ifft

N = 32 # number of samples
n_range = np.arange(N) # indices over which the sequence is defined
n_range_positive = np.arange(int(N/2)+1) # the "positive frequencies" sample indices
n_range_negative = np.arange(int(N/2)+1, N) # the "negative frequencies" sample indices

# generate a complex-valued sequence with the properties expected for the DFT of a real-valued sequence
abs_FFT_positive = np.exp(-n_range_positive**2/100)
phase_FFT_positive =  np.r_[0, np.random.uniform(0, 2*np.pi, int(N/2)-1), 0] # note last frequency has zero phase
FFT_positive = abs_FFT_positive * np.exp(1j * phase_FFT_positive)
FFT_negative = np.conj(np.flip(FFT_positive[1:-1]))
FFT = np.r_[FFT_positive, FFT_negative] # this is the final FFT sequence

# compute the IFFT of the above sequence
IFFT = ifft(FFT)

#plot the results

plt.plot(np.abs(FFT), '-o', label = 'FFT sequence (abs. value)')
plt.plot(np.real(IFFT), '-s', label = 'IFFT (real part)')
plt.plot(np.imag(IFFT), '-x', label = 'IFFT (imag. part)')
plt.legend()

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

Generating correlated random potential using fast Fourier transform

python

arrays

numpy

fft