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.

• I'm nowhere near knowledgable in this area, but I would start by generating white noise and filtering it with a filter whose frequency response $H(f)$ satisfies $|H(f)|^2=S_{xx}(f)$. – MBaz Mar 1 '16 at 17:49
• This previous answer may help. You can get whatever spectral density you want by generating white noise to begin with, then filtering it to the desired shape. – Jason R Mar 1 '16 at 17:49
• The spectral density says nothing about the distribution of values in the time series. A sequence of independent random variables that are uniformly distributed on $[-1,+1]$ has the same spectral density as a sequence of independent $N(0,\sigma^2)$ random variables. – Dilip Sarwate Mar 4 '16 at 3:45

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.

• This seems like a really great answer, I'll start working on it immediately. I'm wondering however if you could perhaps point me towards some theory on how this works. I'm a physicist by training and this whole field is still a bit new to me, and Google is letting me down a bit. In any case definitely something to go from. – user129412 Mar 2 '16 at 12:38
• Some digital signal processing background could help. I can think of dspguide as a free starter resource. Otherwise the Digital filter wiki, FIR wiki and Spectral density wiki are good starting points. Finally, you can find some more notes and links in this other post of mine about the specific $S_y(f) = |H(f)|^2 S_x(f)$ relationship. – SleuthEye Mar 2 '16 at 14:02
• See this wolfram tutorial for a sample usage of ListConvolve for signal filtering. According to your post $f_{max}$ is the maximum frequency you want to model. For a discrete time signal, Nyquist indicates that you must sample your signal in time at least twice as fast. Now, $\Delta f$ in your post correspond to some modeling frequency component (with $\Delta f = f_{max}/K$). This would be the frequency steps between your sampled $S_{xx}(f)$ function for the Frequency sampling-based FIR filter design. – SleuthEye Mar 2 '16 at 20:40
• Sorry for the confusion arising from 'sampling' being done in the time-domain (due to the discrete time signal representation) and in the frequency-domain as part of the FIR filter design method. The "sampling-rate" usually refers to the discrete time signal sampling. – SleuthEye Mar 2 '16 at 20:46
• For time mapping, in general filtering delay of a linear phase FIR filter with $K$ coefficients is $(K-1)/2$. That said, in your case you care about the generated filtered sequence, but not the unfiltered input sequence. So no need to worry about their time correspondence. In that same spirit, to generate $n=N-K+1$ filtered output samples, I'd generate $N=n+K-1$ pseudo-random inputs, then filter them with ListConvolve (and avoid padding). Otherwise, pad the end of the $N$ input samples with $K-1$ zeros to get $(N+K-1)-K+1=N$ outputs. – SleuthEye Mar 3 '16 at 0:14

The spectral density $S_{xx}(f)$ basically defines the magnitude of the signal in the frequency domain. This leaves the phase information to be chosen at will. For this, as mentioned by others, you could start of with white noise and filter it, but you can also start in the frequency domain and use a uniform random distribution of the phase from $0$ to $2\pi$ (or $360^{\circ}$ if you prefer degrees) and take the inverse fft to get a time series.

The first option will probably be easier to implement and more dynamic in length, however the implementation of the filter might not give the exact spectral density you wanted if you use something like a FIR or IIR filter. So the second method might give you more control over it and you can also play around a bit with different phase distributions.

• That's interesting, so you'd create a series in frequency space to begin with and convert it to the time domain, and filter that. It is definitely worth giving a try. – user129412 Mar 4 '16 at 11:17