Working backwards from PSD to possible signal

Question

I have been trying to reconstruct a random signal from its PSD and am running into trouble. I know that many different signals in the time or spatial domains can result in the same PSD-- I am interested in any of these signals. I am not the greatest at signal processing, so I've gotten a bit confused. I am reading this paper. They state:

Therefore, PSD is approximated as $$ S(k)=\frac{2\sigma^2L_c}{(1+k^2L_c^2)^{0.5+\alpha}} $$ where $$ k=2\pi f, f=i/(N\Delta Z), 0\lt i \lt N/2 $$ Lastly, N is the number of points along the line.

...

With the magnitude information provided by S(k), we can reconstruct random line edges by applying a random phase to each frequency component of the PSD to form a unique signal in the frequency domain. A line edge with roughness can be simulated by doing an inverse Fourier transform of this signal. Random lines are distinguished through applied random phases.

I have computed values for S(k), per the equation given. I then tried taking the square root of those values, multiplied them by the number of values, and subtracted their mean. Finally, I generated random phase information before doing an inverse FFT. Even after subtracting the mean I still see a very large spike at the beginning and end of the signal.

Does anybody have any suggestions as to what I'm doing wrong?

For those that speak Python, here is the code I've written so far:

import numpy as np
from matplotlib import pyplot as plt

if __name__ == '__main__':
    SIGMA = 4.0
    ALPHA = 0.2
    Lc = 25

    L = 2000.0
    N = 1000

    dz = L/N

    i = np.linspace(0,N/2,N)
    f = i/(N*dz)
    k = 2*np.pi*f
    PSD = lambda x: (2*3*SIGMA**2*Lc)/(1+(x*Lc)**2)**(0.5+ALPHA)

    window = 25.0/46-21.0/46*np.cos((2*np.pi)/(N-1)*np.array(range(N)))
    magnitude = N*np.sqrt(PSD(k))
    phase = 1j*np.random.randint(0,6,N)+np.random.randn(N)
    FFT = magnitude-np.mean(magnitude)+phase
    FFT = window * FFT
    FFT = np.concatenate((FFT[::-1][0:-1],FFT))

    fig = plt.figure()
    ax = fig.add_subplot(111)

    ax.plot(np.fft.ifft(FFT))

    plt.show()

Edit: I forgot about Gibbs phenomenon! As per Bruce Zenone's suggestion, I applied a Hamming window and mirrored my data. I added plots and updated my Python code. I am still seeing some effects of Gibbs phenomenon. Truncating the affected parts, the signal does not seem to be varying with sigma, alpha, and Lc.

Here is my resulting signal:

And a figure from the paper showing the effect of the different parameters:

Answer 1

There are several issues with your code:

1. The "Gibbs" issue is a non-issue, so you shouldn't need to window.

2. You are better off generating the complete $-\pi$ to $+\pi$ spectrum that you want, rather than relying on getting the symmetry correct.

3. The code

   FFT = magnitude-np.mean(magnitude)+phase
   

   is changing the phase of the PSD magnitude by addition. It should be using complex multiplication:

   FFT = magnitude * np.exp(1j*phase)
   

The code below is a modified version that seems to do the trick. The top graph shows the PSD and the (absolute value of the) FFT, showing that they overlap exactly. The bottom graph shows the generated time domain signal for this particular realization of the phase noise.

---

UPDATE: I just realized that the original code probably shifted the PSD by $\pi$ because the peak was in the center, rather than at the extremities.

---

import math
import numpy as np
from matplotlib import pyplot as plt

    if __name__ == '__main__':
    SIGMA = 4.0
    ALPHA = 0.2
    Lc = 25

    L = 2000.0
    N = 1000

    dz = L/N

    i = np.concatenate([np.linspace(0,N/2,N/2+1), np.linspace(-N/2+1,-1,N/2-1)])

    f = i/(N*dz)
    k = 2*np.pi*f
    PSD = lambda x: (2*3*SIGMA**2*Lc)/(1+(x*Lc)**2)**(0.5+ALPHA)

    magnitude = N*np.sqrt(PSD(k))

    phase = 2*np.pi*np.random.randn(N)
    FFT = magnitude * np.exp(1j*phase)

    fig = plt.figure()


    ax1 = fig.add_subplot(211)
    ax1.plot(np.abs(FFT))
    ax1.plot(magnitude)

    ax2 = fig.add_subplot(212)
    ax2.plot(np.fft.ifft(FFT))

    plt.show()