我想使用指定的自相关函数在1D或2D空间中生成随机势,并且根据一些数学推导(包括Wiener-Khinchin定理和傅里叶变换的性质),可以使用以下方程式完成此操作:

其中phi(k)在间隔[0,1)中均匀分布。并且该函数满足,这是为了确保产生的电势始终是真实的。
自相关函数不应该影响我在这里所做的事情,而是采用简单的高斯分布。

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

  • 相项必须具有1的模数(根据Wiener-Khinchin定理,即函数自相关的傅立叶变换等于该函数的傅立叶变换的模数);
  • 实函数的傅立叶变换必须满足(通过直接检查整数形式的傅立叶变换的定义)。
  • 生成的电势和自相关都是真实的。

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

    有关相关数学,您可以引用以下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()
    

    并且该图类似于如下所示:


    然而,所产生的电势被证明是复杂的,并且虚部与实部具有相同的数量级,这是不期望的。我已经多次检查过数学,但是找不到任何问题。所以我在考虑它是否与实现问题有关,例如数据点是否足够密集以进行快速傅立叶变换等。

    最佳答案

    您对fft(更正确地说是DFT)的工作方式有一些误解。
    首先请注意,DFT假设序列的样本被索引为0, 1, ..., N-1,其中N是样本数。而是生成与索引-10000, ..., 10000对应的序列。其次,请注意,实际序列的DFT将为与0N/2对应的“频率”生成实际值。您似乎也没有考虑到这一点。

    我将不做进一步的详细介绍,因为这不在此stackexchange网站的范围之内。

    仅出于完整性检查的目的,下面的代码会生成一个序列,该序列具有对实值序列的DFT(FFT)期望的属性:

  • 正负频率的共轭对称性
  • 对应于频率0N/2
  • 实值元素
  • 假定
  • 序列对应于从0N-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()
    

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

    关于python - 使用快速傅立叶变换生成相关的随机势,我们在Stack Overflow上找到一个类似的问题:https://stackoverflow.com/questions/57662698/

    10-11 08:38