5

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.

user129412
  • 295
  • 1
  • 4
  • 16
  • 3
    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 01 '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 01 '16 at 17:49
  • 1
    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 04 '16 at 03:45
  • @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" ? – denis Oct 25 '22 at 15:43

4 Answers4

7

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.

SleuthEye
  • 1,792
  • 1
  • 10
  • 13
  • 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 02 '16 at 12:38
  • 1
    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 02 '16 at 14:02
  • I wrote some comments before but they are not really relevant to this question (made a separate question here http://mathematica.stackexchange.com/questions/108881/constructing-the-spectral-density-of-white-noise in case you are interested). One question that is related however is for ListConvolve. What type of overhang would you use? With the default settings you get just a single value as output. Also, on the topic of the Nyquist frequency, what exactly do you mean with $f_{max}$ corresponding to half the sampling rate? I thought the sampling rate was $\Delta f$ – user129412 Mar 02 '16 at 20:08
  • 1
    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 02 '16 at 20:40
  • 1
    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 02 '16 at 20:46
  • Ah, great, I understand I think. So you mean in the sense of $\Delta t = \frac{1}{2f_{max}}$, for the sampling rate. It is my inexperience with the terminology that lead to the confusion. As for the tutorial, I don't really see them doing much with the overhang part, but then again they are doing various things I don't immanently find obvious so I'll spend some time looking at it. – user129412 Mar 02 '16 at 20:51
  • Yes for $\Delta t = \frac{1}{2f_{max}}$. For the overhang you can use zero padding. Note however that you'd typically choose and then fix $K$ to design your filter, then apply it to a sequence of $N>>K$ time-domain pseudo-random samples as you generate your signal (so in that case even with the default no overhang you'd still get $N-K+1 > 1$ output sample). – SleuthEye Mar 02 '16 at 21:07
  • I see, that's interesting. So you make your filter amplitudes and apply it to some long time series ($N$ points, $N > K$) of unitary spectral density white noise. But then how do you know the time that each result corresponds to? Not sure if that makes sense, but every point $n$ in the time series corresponds to some time $t_n$ right, and after the filter you'd still want to ascribe some time value to the output data points – user129412 Mar 02 '16 at 21:16
  • To add to the above, would you just pad the end of the filter amplitudes with zeros until it has the same length as the time series? That's what I'm seemingly understanding from what I have been reading in the sources you provided – user129412 Mar 02 '16 at 22:53
  • 1
    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 03 '16 at 00:14
3

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.

fibonatic
  • 984
  • 5
  • 8
  • 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 04 '16 at 11:17
  • @user129412 No. There is no filtering involved here. Just do the IFFT and done. – tobalt Jul 29 '23 at 04:24
  • 1
    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: https://adsabs.harvard.edu/full/1995A%26A...300..707T.

     – Jean-Baptiste Oct 22 '23 at 15:12
  • Agreed. One degree of freedom per Fourier coefficient (i.e. only phase) will result in a deficient noise representation. It will have a precisely matching PSD in all windows... but that is not a characteristic of real noise. The PSD of a true stationary random process should 'dance' from window to window, but maintain a consistent overall shape when many windows are averaged. – Nathan Wiebe Mar 06 '24 at 22:15
1

I think there are some limitations in the existing answers, so I'll share my two cents. Hopefully there are some thoughts that are helpful.

The FIR-filtered white noise method will work nicely for generating arbitrary length sequences, but if you are trying to simulate a noise source with non-trivial structure in its PSD the spectrum will be difficult to emulate using a FIR filter.

The inverse Fourier method as described by other answers could also use a little further discussion. Randomizing the phase only (and not magnitude) of the Fourier coefficients would result in an infinite sequence with precisely the correct long-term average PSD (good), but also a 100% repeatable PSD in each window of a PSD estimation (bad). The PSD of each window should 'dance' a little over time for a true stationary noise source. So both the magnitude and phase should be randomized. For magnitude, you could square a Gaussian random number with variance of 1 and multiply by the component of the PSD at that frequency, and for phase you could use a uniform random from 0 to 2pi. Or... you could also generate the real and imaginary component independently using a Gaussian random with variance of 1 (magnitude has a variance of 2 at this point), and scale each of them by sqrt(PSD/2) to get to the final magnitude variance equal to the PSD at the each frequency. Note that all of this assumes you have a stable, long-term average of the PSD. If not, it might be worth considering some smoothing.

There are also some potential issues with generating an arbitrary-length sequence of noise from finite-length windows consisting of the inverse Fourier transform on a finite-length PSD. If you use the same input to each successive inverse FFT, you will have a similar problem to the above constant magnitude issue - that each 'random' sinusoidal component in the time domain will continue forever at the same magnitude (and phase in this case). This would cause an autocorrelation to spike at the lag of 1 FFT window, which is clear sign of bad noise.

On the other hand, if you create each successive frequency domain window from a different randomized magnitude/phase number scaled by the PSD value at that frequency, then there are potential discontinuities caused by stitching these together in the time domain. For example, if the PSD has a large spectral contribution from low frequencies (e.g. DC at wavenumber 0, or other low wavenumbers), then there will be an artifact at the boundary between successive windows in the time domain where the magnitude and/or phase abruptly changes from one random value to the next. I'm not sure what the textbook way to mitigate this is, but what comes to mind would be something like Welch's method in reverse: you could generate twice as many noise windows as needed, taper each one by a Hann window, and then stack them with an overlap of half the window length to maintain smooth transitions of phase and magnitude between windows. The multiplication by a Hann window in the time domain should have no smearing effect on the PSD, as it would essentially be equivalent to a convolution by a 1 wavenumber wide Dirac delta in the frequency domain.

For arbitrary-length output, you could probably also just up-sample/interpolate the PSD up to the nearest power of 2 larger than the required noise sequence length, and inverse FFT the whole monstrous thing... but that would scale terribly. Anyways, that's my 2 cents. Very open to thoughts, discussion, and criticism.

Nathan Wiebe
  • 111
  • 1
0

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.

Chris Kaplan
  • 1