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

Question

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.

Answer 1

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)

Answer 2

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.