I have PSD data as a 1D array, function of frequency(rad/sec). It is like powerspectrum for damped oscillator.

PSD formula

Considering 3 modes, PSD is shown here.

I want to generate noise time series which should have PSD exactly like this. Power of the signal should be similar with the power calculated using PSD.

PSD for first 3 modes Here the frequency range is from 10^3 to 10^7 with steps of 10^3.

If some filter will be applied on white noise than it will change the power. One way I thought is to get band limited white noise for band of 1000 rad/sec and in this way I cover 10^3 to 10^7 range and then I will sum up all the noises. But it will be very very time consuming (10^4 such steps) and I also don't know how to get band limited white noise.

I am very new regarding this field. Can someone please tell me, how I can I do this in Mathematica? It will be very helpful to me.

  • $\begingroup$ I don't see a problem filtering white noise to the desired frequency response. You can use your PSD data and InverseFourier[] to transform this back to a time series. Either by starting with a white noise time series, transforming it to frequency space, multiplying by your desired response and transforming back, or by directly starting in frequency space with your PSD data and multiplying rnadom phases and transforming back to a time series. You can also use a FIR filter like described here. $\endgroup$ Commented Nov 3, 2018 at 23:31
  • 1
    $\begingroup$ Also, can you upload an example of your PSD data, e.g. to pastebin.com? This usually makes it easier for people to give you a practical answer to your problem. $\endgroup$ Commented Nov 3, 2018 at 23:34

1 Answer 1


One approach is to take your desired magnitude spectrum and invert it. If you choose the phases randomly, then each time you invert, you will get a different realization, and the realization is guaranteed to have the desired spectrum.

To proceed, first make up a spectrum that looks vaguely like the plot you have shown. I do this by convolving a small blip function with a "spikes" list that places the spikes at locations 100, 400 and 600. This specifies a desired "spectrum" which you will no doubt want to replace with one of your own.

blip = {0, 1, 2, 3, 5, 9, 6, 3, 2, 1};
spikes = RandomReal[{0.001, 0.002}, 1000];
spikes[[{100, 400, 600}]] = {1, 0.5, 2};
spec = ListConvolve[blip, spikes];
ListLogPlot[spec, Joined -> True, PlotRange -> All]

enter image description here

The inverse transform can be calculated as:

randPhase = RandomReal[{0, 2 Pi}, Length[spec]];
compL = spec Exp[I randPhase];
compL2 = Flatten[{compL, Reverse[Conjugate[compL[[Range[2, Length[spec]]]]]]}]; 
outWave = Re[InverseFourier[compL2, FourierParameters -> {1, -1}]];

enter image description here

outWave is your random process. To verify that this has the desired Fourier transform, just take the Fourier transform and verify that it is the same as what you wanted:

      FourierParameters -> {1, -1}]][[1 ;; 1000]], Joined -> True]

You will see that it is identical to the first plot above. Moreover, each time you rerun it, the outWave is different, due to the different randomized phases.

  • $\begingroup$ Thanks @bill. It is working, but the frequencies in the outWave are the element numbers of Spec, like 1,2,......1000. That is because Spec is as usual array of 1 to 1000 elements. If I use Spectrum for frequencies, {1,2,3,4,5........,10^5}, I get perfect result. But, I have spectrum array for the frequencies, in the step of 500 rad/sec : {500,1000, .... 10^5....8 * 10^6}. When I insert this, Mathematica as usual takes up 1 to 16000 numbers for frequency. How I can fed the information of the respective frequency also with all the spectrum value? or if some any other method can be used. $\endgroup$
    – Hari
    Commented Nov 13, 2018 at 17:08

Your Answer

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

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