5
$\begingroup$

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.

Improve this question
$\endgroup$
4
  • 3
    $\begingroup$ 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)$. $\endgroup$
    – MBaz
    Mar 1, 2016 at 17:49
  • $\begingroup$ 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. $\endgroup$
    – Jason R
    Mar 1, 2016 at 17:49
  • 1
    $\begingroup$ 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. $\endgroup$
    – Dilip Sarwate
    Mar 4, 2016 at 3:45
  • $\begingroup$ @Dilip Sarwate, how would one go about generating a time series with a given distribution (e.g. Weibull) and spectral density ? Is that even possible, other than "set up an objective function, run a big optimizer" ? $\endgroup$
    – denis
    Oct 25, 2022 at 15:43

3 Answers 3

Reset to default
7
$\begingroup$

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.

Improve this answer
$\endgroup$
10
  • $\begingroup$ 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. $\endgroup$
    – user129412
    Mar 2, 2016 at 12:38
  • 1
    $\begingroup$ 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. $\endgroup$
    – SleuthEye
    Mar 2, 2016 at 14:02
  • 1
    $\begingroup$ 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. $\endgroup$
    – SleuthEye
    Mar 2, 2016 at 20:40
  • 1
    $\begingroup$ 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. $\endgroup$
    – SleuthEye
    Mar 2, 2016 at 20:46
  • 1
    $\begingroup$ 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. $\endgroup$
    – SleuthEye
    Mar 3, 2016 at 0:14
3
$\begingroup$

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.

Improve this answer
$\endgroup$
3
  • $\begingroup$ 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. $\endgroup$
    – user129412
    Mar 4, 2016 at 11:17
  • $\begingroup$ @user129412 No. There is no filtering involved here. Just do the IFFT and done. $\endgroup$
    – tobalt
    Jul 29 at 4:24
  • $\begingroup$ I'm afraid that this is actually not entirely correct. The Fourier coefficients of a signal $x(t)$ whose underlying process has a given power spectral density (PSD) $S_{xx}(f)$ are not given directly by $S_{xx}(f)$, but are themselves random, and the mean of the distribution is given by the PSD (they are complex, so the real and imaginary parts are). Only randomizing the phase is not enough (you won't sample over the entire space). This is explained here: adsabs.harvard.edu/full/1995A%26A...300..707T. $\endgroup$
    – Jean-Baptiste
    Oct 22 at 15:12
0
$\begingroup$

There's a fairly simple method I've used in the past:

  1. Generate your desired power spectral density in the frequency domain (for example Gaussian, 1/f, etc).
  2. Randomize phase.
  3. FFT back to the time domain.

Since the phase doesn't impact the PSD, you now have a random time domain signal defined by the initial PSD.

Improve this answer
$\endgroup$
1
  • 2
    $\begingroup$ Welcome to SE.SP! How is this different from @fibonatic's answer from Mar 4, 2016? $\endgroup$
    – Peter K. ♦
    Jul 28 at 22:52

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge that you have read and understand our privacy policy and code of conduct.

Not the answer you're looking for? Browse other questions tagged or ask your own question.