Generate a time series from power spectral density

Question

I have a model function for a power spectral density (PSD). I want to simulate a time series with a power PSD equal to my model function. After some browsing about how to do this, I came up with the following method:

1. Discretize my continuous PSD

2. Convert the discrete PSD into amplitude points given by $\sqrt{N F_s S}$, where $N$ is the number of points in the discrete PSD, $F_s$ is the sampling rate (two times the max frequency of the discrete PSD), and $S$ is a point from the discrete PSD.

3. Give each amplitude point a randomly chosen phase. To do this, I multiplied each point by $\exp(-j \phi)$, where $\phi$ is a random value from $0$ to $2\pi$ (chosen from a uniform distribution).

I read some forums where it was advocated that one can now perform an inverse fast Fourier transform and arrive at the appropriate time series. However, the resulting time series is complex, so I added this step (as suggested in a forum) to deal with this:

4. I now have a vector of complex numbers as a function of positive frequency values, running from $0$ to, say, $f_{max}$. I then constructed a new vector with positive and negative frequencies (running from $-f_{max}$ to $f_{max}$) with the following structure. If $A$ is an amplitude point at frequency $f_1$, then the corresponding amplitude point at $-f_1$ is $A^{*}$ (the complex conjugate of $A$).

5. I now perform an inverse fast Fourier transform of this. I still get an imaginary part but its much smaller than the real part so I discard it. I then treat this as my time series. The length is $2 N - 1$ points.

6. To confirm this approach, I compute the PSD (i.e., the periodogram) of the time series and verify that it agrees perfectly with my original continuous PSD model.

I have some questions about this approach:

1. Is this a reasonable method? Can anyone confirm whether it makes sense or whether it is flawed?

2. Is there a reference for this that you can share? If you know of a better approach, would you kindly let me know?

Answer 1

Is this a reasonable method? Can anyone confirm whether it makes sense or whether it is flawed?

Yes this makes sense. The residual imaginary part is probably the Nyquist frequency and/or DC. The spectrum phase at DC and Nyquist must be real, i.e. the phase must be either 0 or $\pi$

Is there a reference for this that you can share?

This method has been around for at least 40 years. It's pretty straight forward, so there aren't a lot of articles about it. This maybe ? https://asa.scitation.org/doi/10.1121/1.397552

If you know of a better approach, would you kindly let me know?

An alternative approach is to design a filter the transfer function magnitude of which matches the square root of your desired PSD and then simply apply it to white noise.

Then has the advantage that it works for arbitrarily lengths of the signal and that every section of the signal will have the same PSD which isn't necessarily guaranteed for the FFT based method.

The downside is that you need to design the filter, which can be tricky depending on the shape of your PSD.

Here is an example for pink noise generation: https://ccrma.stanford.edu/~jos/sasp/Example_Synthesis_1_F_Noise.html