# How to create a wide-sense stationary time series with a frequency of 40 Hz?

I want to create a time series in MATLAB which has a peak frequency of 40 Hz but is also a wide-sense stationary random process. I then want to use power spectral density estimation to recover the spectral content (with a peak at 40 Hz).

Any ideas on how to to do this?

I have tried doing a variety of things including:

1) Create a time series and add random noise. Both the Fourier transform and the PSD recover the true spectrum and the Fourier transform recovers sharper peak. So why use PSD?

2) Create a time series and add a random phase. Fourier transform recovers sharper peak at 40 Hz. Why use PSD?

3) Create a time series with random noise and then randomly sample the time series by pulling data samples from the time series with replacement. This results in random spectra with no peaks using either Fourier transform or PSD.

4) Create a time series with random noise and then re-sample the time series at a lower sample rate. This results in both Fourier transform and PSD recovering the peak at 40 Hz but, once again, the Fourier transform has sharper peaks.

I am trying to create a synthetic time series where PSD estimation is necessary and useful to recover the correct spectral information. But so far I can only create a time series where either a) Fourier transform is superior and recovers sharper peaks; or b) Neither Fourier transform or PSD recover useful spectral information.

Below is an example of a 1 s noisy time series composed of sinusoids with 200, 300 and 400 Hz. The Fourier transform recovers very sharp peaks (e.g. approximately delta functions) whereas the PSD smooths out the peaks due to windowing and averaging.

• Is there some specific issue you are having? This sounds like a do-my-homework question. Jan 31 '20 at 20:14
• Not homework. I'm a geophysicist who often works with time series but has never really been trained with statistics. Our time series are always considered "wide sense stationary" and I'm trying to get a handle on what that actually means in a practical way. I figure that if I can create a synthetic "wide sense stationary" time series then I will have a better understanding of what this type of time series actually is. I've added some more details to the question. Feb 3 '20 at 17:11
• Can you be more specific or post an example of what you mean by sharper peaks? The PSD is often estimated as just the Fourier transform multiplied by its own conjugate. There are other methods for estimating it that can give finer results under certain conditions, but in general they are not that different. And wide-sense stationary just means that the statistics are not a function of time, so, for instance, estimating the noise power over different time intervals should yield similar values. Feb 3 '20 at 17:35
• I added an example showing how the Fourier transform recovers very sharp peaks whereas the PSD smooths out the peaks. Why would someone ever use the PSD when the Fourier transform gives greater accuracy? Feb 3 '20 at 18:21
• How are you estimating the PSD? It looks like it's estimated using a much shorter time span. Can you paste in your matlab code? Feb 3 '20 at 18:29

I'm not an expert in this area but I can think of

• Add random noise to a 40-Hz sine wave: ts = std*randn(1,N) + sin(2*pi*40*t)

• Add a random phase: phi = std*randn(1,N); ts = sin(2*pi*40*t + phi)

• Neither of these work. Feb 3 '20 at 17:09
• Can you give me more details? I've tried both and I know they work.
– MBaz
Feb 3 '20 at 18:41

You do not say what distribution you wish your resulting process to have.

The canonical way to do this is to generate white noise and run it through a filter that has the frequency vs. amplitude characteristic that you desire. If you want a strictly Gaussian output then use Gaussian white noise. If you want a close-to-Gaussian output, then with enough filtering (implying a narrow filter with respect to the sample rate) you could just feed in a random binary sequence and count on averaging to make it more Gaussian.

Assuming that you're doing this on a computer, that you're sampling at appreciably above 40Hz, and that you have a decent random number generator, just generate a long vector of random floats with a normal distribution, and filter it with a bandpass filter centered on 40Hz.

(If you're doing it with analog circuitry then you'd start playing with zener diodes and amplifiers, and discovering the joys of filtering out your local line frequency from your sensitive circuits. I assume this is not going to be an issue with you.)

I've had really good luck so far just using the random number generators that come with whatever math package I'm using -- but be aware that if you want really random data this is not trivial, and as soon as people get pseudo-random generators good enough for all known practical purposes, someone comes along with a practical application that needs even more randomness.