I need to generate random numbers with a given PSD. To do so, I found this recipe. In case it's behind a paywall, here's now it works:
Given your power spectrum $P_k$ with $0\leq k<N$ points generated using a timeseries with timestep $h$, you define
$$A_k=\sqrt{\frac{P_k}{2Nh}}\exp{(ia_k)}$$
where $a_k$ are uniformly randomly chosen in $[0,2\pi)$. The desired timeseries is then the Fourier transform of $A_k$. The idea here is that since the PSD throws away the phase information, we simply have to generate random phases for the transform to get a possible realization of the timeseries that generated our PSD.
However, to ensure that all of the resulting time points are real, the author gives the condition that $A_{N-k}=A^*_k$, or equivalently $a_{N-k}=-a_k$.
When I implement this, however, I don't get a real signal, so I seem to be doing something wrong. I see one problem right away: for $k=0$, the FFT assumes a periodic signal, so the only way to satisfy the condition is if $a_0=0$. Then you can generate the rest of the phases as follows (in python):
def gen_phase(N):
phase = np.zeros(N)
randoms = np.random.random(size=N/2)*2*np.pi
phase[1:N/2+1] = randoms
phase[N/2+1:] = np.flipud(-randoms)
return phase
The other problem is that it is not possible to satisfy the given condition is $N$ is even, as far as I can tell. This is a pain practically, because the fft is very slow for odd numbers (particularly ones of the form $2^n+1$...)
Can someone suggest what I might be doing wrong? The signal I get looks more or less right ion terms of the amplitude, but it has roughly equivalent real and imaginary parts, so something is wrong. The rest of the code is too simple to be the culprit, I think:
def main():
psd = read_psd()
N = len(psd)
phase = gen_phase(N)
h = 0.24e-6
A = np.sqrt(psd/(2*N*h)) * np.exp(1j*phase)
timeseries = np.fft.fft(A)
pl.plot(np.real(timeseries))
pl.show()
if __name__=='__main__':
main()
EDIT: I misinterepreted. Because $A_k$ depends on $P_k$, I don't have the freedom to choose $A_{N-k}=A^*_k$ in the first place, because I don't have control over the amplitude. Any ideas?