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.
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)