Numerically generating noise time series from spectral density I'm trying to figure out a way to construct (numerically) time series of random variables given some spectral density $S_{xx}(f)$. I'm interested in various types of spectral densities (flat, lorentzian, 1/f) so if there is a general method that starts from an expression of $S_{xx}(f)$ and generates a suitable time series that would be ideal of course, but I am not sure if there is.

Honestly I am not so sure where to start. My guess is that you compose your random signal as a sum of different components with amplitudes related to your spectral density at a certain frequency. I understand that if you want your noise to have a complete frequency range you need infinitely many terms, so I suppose you cut off your noise at some frequency $f_{max}$, giving you some maximum amount of terms $K$. You'll also need some step size $\Delta f$ for the different frequency components, but after that I get kind of stuck. How do you incorporate your spectral density into this?

I understand that the question is perhaps a bit too general, so I'm also very much interested in advice on some sources on the subject. The ones I have been able to find are very technical, while for now I'm simply trying to create some time series in Mathematica given some spectral density as input, after which I want to see if I can reconstruct the spectral density from the output. That would get me started quite well, and I can think about what features are essential from there on.
As mentioned by @MBaz in comments, you could use the fact that filtering white noise with unitary power spectral density (for example from independent Gaussian samples zero mean and unitary standard deviation) with a filter having a frequency response $H(f)$ gives you colored noise with power spectral density $S_{xx}(f) = |H(f)|^2$. 

Then comes the question of designing such a filter. A flexible approach would be to use a Frequency sampling-based FIR filter design. Note that this is readily available as FrequencySamplingFilterKernel in Mathematica (or similarly fir2 in matlab for future readers), using an evenly spaced sampling of $|H(f)| = \sqrt{S_{xx}(f)}$ up to the Nyquist frequency (i.e. $f_{max}$ which would correspond to half your time-domain sampling rate). Once you have the FIR filter coefficients, you can use those to filter your white noise with ListConvolve.

Note that the steeper the slopes in your power spectral density, the more frequency samples (resulting in a higher number of FIR filter coefficients) you'll need to get an accurate representation of your power spectral density.