Generating a timeseries with an arbitrary power spectrum 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&lt;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?
You are reversing your phases okay.  You should be able to use an even value of N if you set the phase for N/2 to zero.

The problem you have is that your PSD spans the entire N.  Unless $ P_k = P_{N-k} $ you are not generating a DFT with complex conjugate symmetry which is what you need to produce a real valued signal.

Also, you should be using the "ifft" instead of the "fft" to do an inverse FFT.

Hope this helps,

Ced

Edit:  I didn't see your edit before I posted.  It seems you figured this out.

======================

Followup:

For even N:


def gen_phase(N):
    phase = np.zeros(N)
    randoms = np.random.random(size=N/2-1)*2*np.pi
    phase[1:N/2] = randoms
    phase[N/2+1:] = np.flipud(-randoms)
    return phase