It sounds as though you want to apply an inverse fast Fourier transform (iFFT) to your power spectral density PSD function P(f), in this way transforming from the PSD in frequency domain to time series in time domain, where that time series is random.

To this end, I assume:

(1)    Your given PDS function is sampled at M positive frequencies f(k) = k X Df, k = 1,2,3…M, where Df is the sample spacing of the PSD function in the frequency domain, and M is the number of samples in the PSD;

(2)    The number of samples in your given PDS is M = N/2 + 1 where N is the number of samples in the fast Fourier transform (FFT), N = 256, or 1024, or 2048, … or any other integer power of 2, as usual for the FFT; and

(3)    The maximum frequency in your given PSD is Fmax.

If (1) and (2) are not true, then you must re-sample your PSD such that (1) and (2) are true. If (1) and (2) are true, then you are almost ready to apply an inverse FFT (iFFT) algorithm, which you presumably already have on hand ready to use, in order to generate a corresponding time series.

All you have to do BEFORE applying the inverse FFT (iFFT) is:

(4)    Convert power PSD(f(k)) to amplitude A(f(k)),  A(f(k)) = SQRT[2 X PSD(f(k))];

(5)    Give each spectral component a random phase, PHI(f(k)) = random number, uniformly distributed between 0 and 360 degrees (or equivalently, between 0 and 2Pi radians);

(6)    Construct a frequency domain signal Z(f(k)) = A(f(k)) X e^(i X PHI(f(k))), where i is the imaginary number, i = SQRT(-1)

Use a uniform random number generator to generate the random phases.

In this way, you have gone from PSD(f(k)) to a corresponding signal in frequency domain Z(f(k)), where Z(f(k)) is ready for input into the inverse FFT (iFFT) algorithm for positive valued (one-sided) functions in frequency f(k) > 0.

Now, finally, apply the inverse FFT (iFFT) to Z(f(k)).

The output of the iFFT will be a random time series on the finite (bounded, limited) time interval t = 0 to tmax = (N-1) X Dt, where Dt = 1 / (2 X Fmax).

The PSD of this time series will be the same as your given PSD.  You should check your result, of course, by taking the FFT of the resulting random time series in order to confirm that it is giving the same PSD as you given PSD.

Note that the time series will repeat itself (no longer random), at times greater than tmax, and at times less than 0. So you cannot generate ever longer time series this way unless you increase N in the sampling of your PSD and in the inverse FFT accordingly.

You can always generate a new random time series by generating a new set of random phases PHI(f(k)) = random number, uniformly distributed between 0 and 360 degrees (or equivalently, between 0 and 2Pi radians), and then applying the iFFT once again. This assumes that the seed of the random number generator is different each time you undertake to generate a new random time series from the given PSD using the steps above.

I hope that helps.

Ronald

It sounds as though you want to apply an inverse fast Fourier transform (iFFT) to your power spectral density PSD function P(f), in this way transforming from the PSD in frequency domain to time series in time domain, where that time series is random.

To this end, I assume:

(1)    Your given PDS function is sampled at M positive frequencies f(k) = k X Df, k = 1,2,3…M, where Df is the sample spacing of the PSD function in the frequency domain, and M is the number of samples in the PSD;

(2)    The number of samples in your given PDS is M = N/2 + 1 where N is the number of samples in the fast Fourier transform (FFT), N = 256, or 1024, or 2048, … or any other integer power of 2, as usual for the FFT; and

(3)    The maximum frequency in your given PSD is Fmax.

If (1) and (2) are not true, then you must re-sample your PSD such that (1) and (2) are true. If (1) and (2) are true, then you are almost ready to apply an inverse FFT (iFFT) algorithm, which you presumably already have on hand ready to use, in order to generate a corresponding time series.

All you have to do BEFORE applying the inverse FFT (iFFT) is:

(4)    Convert power PSD(f(k)) to amplitude A(f(k)),  A(f(k)) = SQRT[2 X PSD(f(k))];

(5)    Give each spectral component a random phase, PHI(f(k)) = random number, uniformly distributed between 0 and 360 degrees (or equivalently, between 0 and 2Pi radians);

(6)    Construct a frequency domain signal Z(f(k)) = A(f(k)) X e^(i X PHI(f(k))), where i is the imaginary number, i = SQRT(-1)

Use a uniform random number generator to generate the random phases.

In this way, you have gone from PSD(f(k)) to a corresponding signal in frequency domain Z(f(k)), where Z(f(k)) is ready for input into the inverse FFT (iFFT) algorithm for positive valued (one-sided) functions in frequency f(k) > 0.

Now, finally, apply the inverse FFT (iFFT) to Z(f(k)).

The output of the iFFT will be a random time series on the finite (bounded, limited) time interval t = 0 to tmax = (N-1) X Dt, where Dt = 1 / (2 X Fmax).

The PSD of this time series will be the same as your given PSD.  You should check your result, of course, by taking the FFT of the resulting random time series in order to confirm that it is giving the same PSD as you given PSD.

Note that the time series will repeat itself (no longer random), at times greater than tmax, and at times less than 0. So you cannot generate ever longer time series this way unless you increase N in the sampling of your PSD and in the inverse FFT accordingly.

You can always generate a new random time series by generating a new set of random phases PHI(f(k)) = random number, uniformly distributed between 0 and 360 degrees (or equivalently, between 0 and 2Pi radians), and then applying the iFFT once again. This assumes that the seed of the random number generator is different each time you undertake to generate a new random time series from the given PSD using the steps above.

I hope that helps.

Ronald

Hi Manoj,

Roland answer is an excellent one. I believe it can be more accurate with few slight modifications.

First we should return to the PSD definition. power spectral density is a representation of statistical power distribution across frequency range of a specific signal. Not an actual one to one representation of the signal in frequency domain as in FFT.

Although the PSD provides some useful information (namely the average power at a specific frequency) still it lacks two important characteristics.

1. The statistical distribution of the signal power (not just the mean)

Actually it's exactly the variance of the signal amplitude at each specific frequency

2. The correlation between different frequencies (if any)

3. The phase data for each frequency.

If we have any of these information we better try to incorporate it the modeled signal. Otherwise we better assume it to have a totally random values, as ronald was assuming with the phase.

What can be added to Ronald model is to have each bin amplitude as a random value as well with a variance equal to the value of the PSD bin at these frequency.

To do this you can go through two different methods.

The first one is just an extension to Ronald method.

1. get a random number generator to generate the random phases for each bin.

2. Get another random generator to generate the amplitude of the bins 

3. Multiply the generated amplitude with the root square of the PSD amplitude (variance to standard variation)

4.perform ifft on the genrated sequence to get its time domain counterpart.

5. Repeat 1->4 multiple times concatenate each time domain sequence with the preceding one to get real random signal with the required PSD

Another method which I prefer as it requires less computational power (assuming you have dedicated FFT/IFFT optimized functions)

1. take the square root of the PSD.

2. Generate random numbers with a variance of 1 with the length of the PSD pins

3. Take the FFT of this random sequence.

4. multiply this sequence with the square root of the PSD

5 take the ifft of the product

6. repeat 2->5 and concatenate the newer sequences.

Having a real random generator for the raw time domain signal provides a white noise with a variance equal to 1 (distributed equally across the whole frequency range, you might need to scale for that later) with random phases between +/- pi.

So we don't need two random generators for the amplitude and phase.

Thanks,

Mina

Hi Manoj,

I just noticed a missing note in my answer. Since I'm just multiplying by the sqrt of the PSD amplitude the result in time domain will be a complex signal.

To enforce the time domain to be signal we need to have a DSB PSD signal and forcing a symmetry between the USB and the LSB before multiplying its square root by the fft of the white noise.

If you want to go with Ronald answer you should as well insure a DSB representation and enforce the LSB to be the complex conjugate of the USB (as he is adding phase information to the PSD)

Thanks,

Mina 

Hello. In addition to the well thought out answers already given, I would add a couple of points. First, there is no unique random signal for a given PSD. Second, the PSD of measured data is an estimate. Without constraining the problem, you will have an infinite number of possible time series with the same PSD even in the absence of noise. With noise, it's worse than this since the PSD estimate may be a poor approximation of the actual PSD of the system. Many PSDs can be generated by putting white noise into a linear system. If you fit the PSD of a linear dynamical model (auto-regressive moving average) to your PSD you will get a representation for the dynamical system you have (this is related to Burg's method). You can then just use white noise as an input and adjust the variance to conserve total signal power. There is a good review of this type of stuff in the IEEE by Childers (if memory serves).

I have studied your answers. Thanks for your thoughtful response. I am in the process of applying your suggestion to my work and will update in due time about the results. Thanks to Kessel, Mina and Pulford.

It depends on the type of the signal. For complex signals (like IQ data for digital TRXs, base band representation of oscillator Phase noise, ..) You have to go with the definition of the frequency range as -FS/2 --> FS/2 (where FS is the sampling frequency and the number of samples "frequency resolution" will depend on the time-domain sequence length). And the PSD is not required to by symmetric between +ve and -ve frequencies.

On the other hand for real signals then you have one of two opetions. used again the -FS/2 to FS/2 and enforce that the -ve frequency component are just the complex conjugate of the +ve frequency components. Or just go from 0 to FS/2 (i.e. +ve frequency components and DC).

which is better will depend on what IFFT engine you have at hand. I use matlab so typically I've +/-ve frequency components but rotated to be from 0 to FS.

How to verify your results is by having some PSD algorithm in place. Typically again I use the built in PSD functions in Matlab. which is basically dividing the stream into multiple smaller parts each with the length of the FFT. Do the FFT for each of them and then doing some kind of an averaging between them. As for the scaling factors I always forget them so I always check the power in frequency domain and time domain and make sure that they are equal.

Thanks,

Mina

Hello,

in the procedure described by Ronald the iFFT output values have imaginary part whereas time domain values are real usually. What should be done with imaginary part of iFFT output values?

Thank you

Oleg

Hi Oleg,

Having real or complex time domain series is based on the symmetry in the frequency domain series "between positive and negative frequency".

For Ronald approach he is assuming one sided IFFT algorithm. (So basically assuming the FFT is done with the sin + cos coefficients to control the phase instead of exp(1i) )

If the IFFT algorithm you have is the general complex double sided one then you need an extra conditioning step.

After you create the positive frequency pins (from 0 to N/2) you calculate the negative pins as the complex conjugate (C of -M = conj(C of M).

Then apply the IFFT.

Thanks

Mina

Thank you very much Mina.

We simulate now a response of our future Silicon Photomultiplier and now I need to prepare a noise (real) signal simulation, the frequency density function of this noise is supposed to be fixed later and its bandwidth is known now. I've got the public FFTW C++ library. Could you recommend me the way I could simulate a (real) time domain signals of this noise?

Thank you again

Oleg

Hello Ronald and Mina, I appreciate the methodology the two of you described on this topic as it has aided me greatly. However I was searching for publications that somehow certified the effectiveness of the output signal and I didn't find any, besides classic signal processing texts in which the formulae are described.

Would you have any sources of that type that could validate the method?

Thank you very much,

Niago

I don't know of a published reference on this methodology.

But none is required. You can (and should) always validate the method yourself when you use it. That is: compute the PSD for a time series generated using the method. Does it give the original PSD? If so, then you have validated the method for your signals.

Keep in mind that the PSD for a random signal will also be somewhat random. So if you use the method to generate many time series, each generated using the same intended PSD, then the average of the set of validating PSDs for each time series should converge to the original PSD used to generate all of the time series in your validation test.

Yes, I validated it that way! I was just wondering if there was additional reference on the matter.

Yes, I am aware of the randomness of the output PSDs.

Thank you again,

Niago

You may try the paper "Free-Surface Time-Series Generation for Wave Energy Applications" by Alexis Merigaud and John V. Ringwood.

Hi Mina and Ronald,

Let say M=100;

Do I need to resample it to be 129 (to satisfy N = 256 requirement)?

These days matlab can perform FFT without having N equals to the power of 2.

Another question is after constructing a frequency domain signal Z (even number=100) follow step 6 recommended by Ronald from steps (4) to (6) then how do we add the complex conjugate to the frequency domain signal. Can we simply add it as a ro w or column vector?

(7) New_Z=[Z; conj(Z)]

For Ronald approach, The output of the iFFT will be a random time series on the finite (bounded, limited) time interval t = 0 to tmax = (N-1) X Dt, where Dt = 1 / (2 X Fmax).

Can this statement still true?

Hi Samdani,

You can use 100. the 128, 256, ... 2^N count is just the fastest scenario.

o If you give matalb a 100 samples and 128 samples the 128 will be faster.

For adding the conjugate it should be appended in the same direction as the Z.

So if Z is a row you do it like [Z conj(Z)]. If it's a column then [Z;conj(Z)]

Thanks,

Mina

Hi Mina,

why do I still get complex values after taking ifft of [Z;conj(Z)] ? Shouldn't it be real only?

All the best for your project dear

Did you balance the symmetry in the real and imaginary parts of the noise transform? The 0th Fourier coefficient needs to be real. For all others, using Python indexing convention, the real part of the nth coefficient equals the real part of the -nth. The imaginary part of the nth coefficient equals minus the imaginary part of the -nth coefficient.

Dear Ronald Kessel ,

thank you very much for the thorough explanation! I'm only struggling with the units here...

Assuming that the PSD is given in Amplitude^2 per Hertz, then step 4) in your explanation yields Amplitude per sqrt(Hz). Does the regular iFFT not run on Amplitude per Hertz in order to get back to the time domain where simply Amplitude vs time?

@another, probably naive approach:

Treat the PSD as the FFT of the amplitude^2 function, do the regular inverse of this and take the sqrt in the time domain. Would this work?

Thanks again.

Jan

Hi @Ronald Kessel, really appreciate your answer. I just wonder if the method also applies to the band limited frequency range . Actually, I encounted the problem of transforming the railway PSD to excitation in time domain but there are both the low limit and up limit of the frequency range. If the method is used， does it assumes there are only input energy within the frequency range while the rest frequencies' energy contrubution is zero？

Yes, the method works on band-limited PSD. Just zero the PSD outside the frequency band. Or you might choose to smoothly window the PSD in the frequency domain for physical reasons.

Noted and thank you!

Ronald Kessel Hi, I was wondering can we obtained the original time series data from a given PSD if we followed the steps you mentioned? Thank you.

Generally no, you cannot obtain the original time series from the PSD using this method. The phase of the power spectrum is randomized, in this method giving a new time series each time is it applied. Only by some very unlikely chance would the original time series result.

Dear all

I am implementing the method of Mina Abdallah (and Ronald Kessel ) to generate time series for an analytical known PSD (noise spectrum from an RC circuit).

It works nicely, I get the correct shape but I am having problems getting the correct amplitude.

I have analyzed the generated time series using a the matlab “periodogram“ function and a simpleminded FFT based calculation. The two estimates agree.

However, there is a difference of a factor Delta omega (the spacing on the frequency axis) between the estimated PSD and the analytical.

I am guessing I am missing some normalization factor somewhere, could you please give a hint two what I am overlooking.

Below I have included a somewhat long minimal example written in Matlab:

M=2^16/2; % Number of points in the PSD

fmax=100 % Maximal frequency in the PSD

omegalin=linspace(0,fmax,M); %frequency axis for evaluating the PSD

%PSD function evaluated at the discreet points

%Simpel RC paralel led

PSDfkt=1./(1+(omegalin).^2);

%Add the symmetric part to the PSD function needed for getting real

% values in the time domain

PSDfkt2=[PSDfkt fliplr(PSDfkt(2:end))]

%generate time series based on the method of Mina Abdallah

rn=randn(1,length(PSDfkt2));

sig_t=ifft(fft(rn).*sqrt(PSDfkt2));

sig=sig_t(1:M);

%corresponding sample frequency in the time domain

fsamp=2*fmax;

%use matlab periodogram to estimate PSD

[pest,fest]=periodogram(sig,[],2^16,fsamp)

%manual FFT based estimate

% see https://en.wikipedia.org/wiki/Spectral_density#Power_spectral_density

pest2=fft(sig).*conj(fft(sig))

pest2=2*pest2(1:M)/fsamp/length(sig)

fest2=fsamp*(0:(M-1))/M;

%spacing of omega bins

Domega=fsamp/M;

figure; hold on

plot(log10(fest),log10(pest))

plot(log10(fest2),log10(pest2))

plot(log10(omegalin), log10(PSDfkt))

plot(log10(omegalin), log10(PSDfkt*Domega))

% Note that the perodogram and FFT based methods agree

% The analytical function is off by a factor Domega

Hi Bo,

definitely you will get this scaling factor of Delta Omega (I think it should be delta F in Hz not rad/sec).

The reason is that the for PSD, the Power of the signal at specific frequency is the pin area Not amplitude. And this will be the area of a rectangle “amplitude x DF“. while ifft has no knowledge of this DF. It just assumes some FS=1 in a sampled system. and DF of 1/N, where N is the total length.

what I typically do when I need to re-write similar algorithms is to do rms(the signal) and check the target power. To make sure that I didn’t miss a scaling factor. In frequency domain this will be SUM(PSD”Not in dB but in watts”)xDF. Typical mistakes I usually do is related to N or DF, Or their square root based on the implementation.

hope this helps

thanks

Mina

Thank you Mina Abdallah for the fast reply.

I will try to evaluate the total energy in the analytic function and the esimates.

Best

Bo

Bo Jakobsen I am facing a similar problem to the one you did. I am taking the the IFFT of a PSD and am not getting the correct RMS in the time domain. Did you ever solve your issue?