Generating correlated random potential using fast Fourier transform

Question

I would like to generate a random potential in 1D or 2D spaces with a specified autocorrelation function, and according to some mathematical derivations including the Wiener-Khinchin theorem and properties of the Fourier transforms, it turns out that this can be done using the following equation: where phi(k) is uniformly distributed in interval [0, 1). And this function satisfies , which is to ensure that the potential generated is always real. The autocorrelation function should not affect what I am doing here, and I take a simple Gaussian distribution .

The choice of the phase term and the condition of phi(k) is based on the following properties

1. The phase term must have a modulus of 1 (by Wiener-Khinchin theorem, i.e. the Fourier transform of the autocorrelation of a function equals the modulus of the Fourier transform of that function);

2. The Fourier transform of a real function must satisfy (by directly inspecting the definition of Fourier transform in integral form).

3. Both the generated potential and the autocorrelation are real.

By combining these three properties, this term can only take the form as stated above.

For the relevant mathematics, you may refer to p.16 of the following pdf: https://d-nb.info/1007346671/34

I randomly generated a numpy array using uniform distribution and concatenated the negative of the array with the original array, such that it satisfies the condition of phi(k) stated above. And then I performed the numpy (inverse) fast Fourier transform.

I have tried both 1D and 2D cases, and only the 1D case is shown below.

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

And the plot is similar to what is shown as follows: .

However, the generated potential turns out to be complex, and the imaginary parts are of the same order of magnitude as that of the real parts, which is not expected. I have checked the math many times, but I couldn't spot any problems. So I am thinking whether it's related to the implementation problems, for example whether the data points are dense enough for Fast Fourier Transform, etc.

Answer 1

More care needs to be taken when concatenating:

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

The first element is the offset (zero frequency mode). "Negative" frequencies come after the midpoint.

This gives me

Answer 2

You have a few misunderstandings about how fft (more correctly, DFT) operates. First note that DFT assumes that the samples of the sequence are indexed as 0, 1, ..., N-1, where N are the number of samples. Instead, you generate a sequence corresponding to indices -10000, ..., 10000. Second, note that the DFT of a real sequence will generate real values for the "frequencies" corresponding to 0 and N/2. You also seem to not take this into account.

I won't go into further details as this is out of the scope of this stackexchange site.

Just for a sanity check, the code below generates a sequence that has the properties expected for the DFT (FFT) of a real-valued sequence:

• conjugate symmetry of positive and negative frequencies,
• real-valued elements corresponding to frequencies 0 and N/2
• sequence assumed to correspond to indices 0 to N-1

As you can see, the ifft of this sequence indeed generates a real-valued sequence

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