# Working backwards from PSD to possible signal

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.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:

-
Given a known signal you can calculate it's Fourier transform, then estimate its power spectral density-- I am looking to go in the reverse direction. I have the PSD in an analytic form and am looking to obtain any signal that fits that PSD. –  Mr. Squig Jan 30 '13 at 16:41
You need to know the phase to reconstruct the signal. Random phase isn't going to cut it. –  Jim Clay Jan 30 '13 at 17:32
Is S(k) a sampled version of a continuous power spectrum? If so, do you need to window your S(k) based samples (i.e. the numbers you end up with after procesing) to avoid Gibbs phenominon when applying the IFFT? –  user2718 Jan 30 '13 at 17:35
Note that any given PSD corresponds to an infinite number of time-domain signals, because the phase component is discarded. You can take the same PSD, and assuming different phase relationships, you can end up with wildly different time-domain signals. As an example: a sample of white noise and an impulse have the same power spectral density, but they don't look anything alike in the time domain. –  Jason R Jan 30 '13 at 18:21
The OP has a paper that explains his process. Indeed the PSD has lost its phase info, so effectively this sets all phases to 0. The IFFT is fully qualified for this. You get one answer, the answer that has all 0 phase. Now the OP wants to simulate deviations from the 0 phase case by adding in some low level random phase deviations to the spectral data before applying the IFFT. I think the issue is in how the phase data is combined with the base spectrum –  user2718 Jan 30 '13 at 22:32

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