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

Question

我想要在1D或2D空间中生成具有指定自相关函数的随机势能。经过一些数学推导，包括维纳-金钦定理和傅里叶变换的属性，发现可以使用以下方程进行操作： 其中， phi(k) 在区间[0,1)内是均匀分布的。这个函数需要满足条件 ，以确保生成的势能始终为实数。 自相关函数不应影响我在此处执行的操作，我选取了简单的高斯分布。

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

1. 相位项必须具有模1 (根据维纳-金钦定理，即函数的自相关的傅里叶变换等于该函数的傅里叶变换的模);

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

以下是相似的情节：

然而，所生成的电位是复杂的，虚部与实部的数量级相同，这是意料之外的。我已经多次检查了数学计算，但没有发现任何问题。因此，我在考虑它是否与实现问题有关，例如数据点是否足够密集进行快速傅里叶变换等。

- Sato

Answer 1

在连接字符串时需要更加小心：

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

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

这给了我

- Paul Panzer

---

区别在于我在0处进行镜像：phi1、phi2、phi3 -> -phi3，-phi2，-phi1，0，phi1，phi2，phi3等（它绕过零点），而你在两个网格点之间进行镜像：phi1、phi2、phi3 -> -phi3，-phi2，-phi1，phi1，phi2，phi3，这不是fft中频率排列的方式。 - Paul Panzer

我明白了。那左边的第一个[0]呢？我知道如果没有最左边的[0]，维度就不对了，但为什么它必须是零呢？ - Sato

我所说的零是指左侧的零。按照惯例，左边的所有内容都放在末尾。如果fft的顺序是偶数，则会出现另一个零。它必须为零，因为它是phi和-phi的自身模式（位于n/2处，因此你将得到类似于exp（2pi k（n/2）/n）的内容，这恰好与exp（2pi k（-n/2）/n）相同）。 - Paul Panzer

谢谢。这应该如何调整为二维情况？我尝试以类似的方式连接生成的矩阵：统一生成矩阵 matrix1 和 matrix2，然后按正确的维度沿 axis=0 连接 np.zeros, matrix1, np.zeros, matrix2，然后定义 matrix3 = np.flipur(flipld(matrix2)) 和 matrix4 = np.flipur(flipld(matrix1))，并沿 axis=1 连接上述两个连接的矩阵，并在最上面加上一个 np.zeros，中间再加上一个 np.zeros。但是那似乎效果不好。 - Sato

我将这些矩阵都翻转了两次，因为根据数学原理，phi(-k) = -phi(k)，而当k是多维的时候，“另一侧的元素”应该对应于关于“原点”k = (0,0)的镜像。 - Sato

Answer 2

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

我不会深入讨论此问题，因为这超出了此站点的范围。

只是为了进行健全性检查，下面的代码生成具有实值序列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()

- Stelios

---

感谢您的回答。这应该如何调整为二维情况？我尝试类比地将均匀生成的矩阵按以下方式连接起来：均匀生成矩阵matrix1和matrix2，然后沿着axis=0连接np.zeros、matrix1、np.zeros、matrix2，在正确的维度上进行连接，然后定义matrix3=np.flipur(flipld(matrix2))和matrix4=np.flipur(flipld(matrix1))，并沿着axis=1连接上述两个连接好的矩阵，其中最上面是一个np.zeros，中间是一个np.zeros。但是这似乎效果不佳。 - Sato

我猜，然而，不必将每个样本对应到一个整数索引上吧？我索引的方式也像Paul Panzer在另一个答案中展示的那样有效。 - Sato

@Sato，将每个样本映射到整数索引0和N-1是必要的。这是fft/ifft所期望的。我在另一个答案中看到了非零虚部，因此我不明白你为什么说“它起作用了”。我强烈建议您先熟悉一维情况，然后再进行二维情况。一旦您完成了这个过程，做二维情况就会变得简单（尽管繁琐）。 - Stelios

当我调整了连接数组的方式并采取 x_min，x_max，interval_x = -10，10，10000 后，实部的数量级为1e-1，而虚部为1e-19，因此我猜测通过在一个较窄的范围内添加更多样本来使其工作？ - Sato

谢谢。还有一个问题，输入零的数组是（如您的代码中所示）phase_FFT_positive，而不是用于执行ifft的数组。在您的代码中，FFT类似于[1，...，1 ...]，而不是[0，...，0，...]，但这似乎不正确？因为据我所知，应该直接输入进行fft或ifft的数组应具有零填充。但是在您的代码中，相位具有零填充。这是个问题吗？ - Sato