Thread Subject: Generate a time series from a PSD

Hi,

I tried to generate a time series (one realization) from a Power Spectral Density (PSD). Since the PSD only contains amplitude information but no phase information, I defined the phase (at each frequency) as a random variable, uniformly distributed between 0 and 2*pi. Basically the code is:

============================================================
function [x, t] = psd2timeseries(PSD, fNyq)
% x = time series
% t = time vector
% PSD = power spectral density (one-sided, from 0 to fNyq)
% fNyq = nyquist frequency

N = length(PSD);

% Compute amplitude of frequency spectrum
SpectrumAmplitude = sqrt(PSD);

% Compute phase of frequency spectrum
SpectrumPhase = rand(size(PSD))*2.*pi;

% Compute complex spectrum
Spectrum = SpectrumAmplitude .* exp(i*SpectrumPhase);

% This is a one-sided PSD, for frequencies between 0 and
% the Nyquist frequency fNyq. We construct a two sided-PSD:

% Sampling of frequency vector
delta_f = fNyq/(N-1);

% Generate two-sided frequency vector (as it would result
% from Matlab's fft command for an even number of samples):
f_TwoSided = ...
    [...
    [0 : delta_f : fNyq] ...
    [-fNyq + delta_f : delta_f : -delta_f]...
    ];
% ... and generate the two-sided PSD:
Spectrum_TwoSided = [Spectrum(1:1:end) Spectrum(end-1:-1:2)];

% Compute inverse FFT
x= ifft(Spectrum);
x= fftshift(x);
x= real(x);

% Compute time vector
delta_t = 1./(2.*fNyq);
t= [0:delta:(N-1)*delta];
===========================================================

What I am not sure about is:
* if it is required to build a two-sided PSD prior to computing the inverse FFT with ifft.
* In which ordering the time vector of the resulting time series x (from ifft) will be computed. Is fftshift required?
* The problematic code seems to be:

x= ifft(Spectrum);
x= fftshift(x);
x= real(x);

In any case the function does not produce the results I expected. When I am trying, once x and t have been computed, to compute a periodogram using

[Pxx,freq] = periodogram(x,hanning(length(x)),length(x),2*fNyq);

the resulting periodogram should resemble to PSD but it does not (at all!).

Any ideas what went wrong?

Thanks
Rainer

"Rainer " <wilhelm.remove.this.rainer@gmx.net> wrote in message <gao08m$ad3$1@fred.mathworks.com>...
> Hi,
>
> I tried to generate a time series (one realization) from a Power Spectral Density (PSD). Since the PSD only contains amplitude information but no phase information, I defined the phase (at each frequency) as a random variable, uniformly distributed between 0 and 2*pi. Basically the code is:
>
> ============================================================
> function [x, t] = psd2timeseries(PSD, fNyq)
> % x = time series
> % t = time vector
> % PSD = power spectral density (one-sided, from 0 to fNyq)
> % fNyq = nyquist frequency
>
> N = length(PSD);
>
> % Compute amplitude of frequency spectrum
> SpectrumAmplitude = sqrt(PSD);
>
> % Compute phase of frequency spectrum
> SpectrumPhase = rand(size(PSD))*2.*pi;
>
> % Compute complex spectrum
> Spectrum = SpectrumAmplitude .* exp(i*SpectrumPhase);
>
> % This is a one-sided PSD, for frequencies between 0 and
> % the Nyquist frequency fNyq. We construct a two sided-PSD:
>
> % Sampling of frequency vector
> delta_f = fNyq/(N-1);
>
> % Generate two-sided frequency vector (as it would result
> % from Matlab's fft command for an even number of samples):
> f_TwoSided = ...
> [...
> [0 : delta_f : fNyq] ...
> [-fNyq + delta_f : delta_f : -delta_f]...
> ];
> % ... and generate the two-sided PSD:
> Spectrum_TwoSided = [Spectrum(1:1:end) Spectrum(end-1:-1:2)];
>
> % Compute inverse FFT
> x= ifft(Spectrum);
> x= fftshift(x);
> x= real(x);
>
> % Compute time vector
> delta_t = 1./(2.*fNyq);
> t= [0:delta:(N-1)*delta];
> ===========================================================
>
> What I am not sure about is:
> * if it is required to build a two-sided PSD prior to computing the inverse FFT with ifft.
> * In which ordering the time vector of the resulting time series x (from ifft) will be computed. Is fftshift required?
> * The problematic code seems to be:
>
> x= ifft(Spectrum);
> x= fftshift(x);
> x= real(x);
>
> In any case the function does not produce the results I expected. When I am trying, once x and t have been computed, to compute a periodogram using
>
> [Pxx,freq] = periodogram(x,hanning(length(x)),length(x),2*fNyq);
>
> the resulting periodogram should resemble to PSD but it does not (at all!).
>
> Any ideas what went wrong?
>
> Thanks
> Rainer

why bother with a random phase? it would seem to add complexity while you are trying to figure out the basics. I would leave it as zero until you get the basic ifft working then consider whether you really have to muck with it at all.

also, why are you doing the fftshift after the ifft? fftshift takes the results of doing an fft and swaps the vector halves so the zero frequency is in the middle... doing the fftshift after an ifft swaps the left and right halves of the time sequence which wouldn't seem to make sense. you might need to do ifftshift before the ifft if your psd that you have duplicated has the zero frequency in the middle of the vector.

On Sep 16, 10:03=A0pm, "Rainer " <wilhelm.remove.this.rai...@gmx.net>
wrote:
> Hi,
>
> I tried to generate a time series (one realization) from a Power Spectral=
 Density (PSD). Since the PSD only contains amplitude information but no ph=
ase information, I defined the phase (at each frequency) as a random variab=
le, uniformly distributed between 0 and 2*pi. Basically the code is:
>
> =3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=
=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=
=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D
> function [x, t] =3D psd2timeseries(PSD, fNyq)
> % x =3D time series
> % t =3D time vector
> % PSD =3D power spectral density (one-sided, from 0 to fNyq)
> % fNyq =3D nyquist frequency
>
> N =3D length(PSD);
>
> % Compute amplitude of frequency spectrum
> SpectrumAmplitude =3D sqrt(PSD);
>
> % Compute phase of frequency spectrum
> SpectrumPhase =3D rand(size(PSD))*2.*pi;
>
> % Compute complex spectrum
> Spectrum =3D SpectrumAmplitude .* exp(i*SpectrumPhase);
>
> % This is a one-sided PSD, for frequencies between 0 and
> % the Nyquist frequency fNyq. We construct a two sided-PSD:
>
> % Sampling of frequency vector
> delta_f =3D fNyq/(N-1);
>
> % Generate two-sided frequency vector (as it would result
> % from Matlab's fft command for an even number of samples):
> f_TwoSided =3D ...
> =A0 =A0 [...
> =A0 =A0 [0 : delta_f : fNyq] ...
> =A0 =A0 [-fNyq + delta_f : delta_f : -delta_f]...
> =A0 =A0 ];
> % ... and generate the two-sided PSD:
> Spectrum_TwoSided =3D [Spectrum(1:1:end) Spectrum(end-1:-1:2)];
>
> % Compute inverse FFT
> x=3D ifft(Spectrum);
> x=3D fftshift(x);
> x=3D real(x);
>
> % Compute time vector
> delta_t =3D 1./(2.*fNyq);
> t=3D [0:delta:(N-1)*delta];
> =3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=
=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=
=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D
>
> What I am not sure about is:
> * if it is required to build a two-sided PSD prior to computing the inver=
se FFT with ifft.
> * In which ordering the time vector of the resulting time series x (from =
ifft) will be computed. Is fftshift required?
> * The problematic code seems to be:
>
> x=3D ifft(Spectrum);
> x=3D fftshift(x);
> x=3D real(x);
>
> In any case the function does not produce the results I expected. When I =
am trying, once x and t have been computed, to compute a periodogram using
>
> [Pxx,freq] =3D periodogram(x,hanning(length(x)),length(x),2*fNyq);
>
> the resulting periodogram should resemble to PSD but it does not (at all!=
).
>
> Any ideas what went wrong?
>
> Thanks
> Rainer

Fundamental error:
x=3D ifft(Spectrum);
should be:
x=3D ifft(Spectrum_TwoSided);
You went to all the trouble of making Spectrum_TwoSided, then failed
to use it.

One thing I noticed with your algorithm is that the phase at zero
frequency in Spectrum is nonzero. It should be zero. I'm not sure
what effect this will have.
Also, you need to be careful how you butt the two parts of the
spectrum together:
Spectrum_TwoSided =3D [Spectrum(1:1:end) Spectrum(end-1:-1:2)];
This will yield an odd number of data in x.
You may like to try inserting a zero between them. This will give an
even number in x.

Hello Rainer,
             Is the code working for you with the suggestion mentioned above.Can you please the final working code.I am trying to work on a similar application and it will be great if you dont mind sharing the work you have carried out.

Thanks and Best Wishes.

On Sep 16, 6:03=A0am, "Rainer " <wilhelm.remove.this.rai...@gmx.net>
wrote:
> Hi,
>
> I tried to generate a time series (one realization) from a Power Spectral=
 Density (PSD). Since the PSD only contains amplitude information but no ph=
ase information, I defined the phase (at each frequency) as a random variab=
le, uniformly distributed between 0 and 2*pi.

Wouldn't it be better to estimate phase using log(PSD) and
a Hilbert transform?

Hope that helps.

Greg


Basically the code is:
>
> =3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=
=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=
=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D
> function [x, t] =3D psd2timeseries(PSD, fNyq)
> % x =3D time series
> % t =3D time vector
> % PSD =3D power spectral density (one-sided, from 0 to fNyq)
> % fNyq =3D nyquist frequency
>
> N =3D length(PSD);
>
> % Compute amplitude of frequency spectrum
> SpectrumAmplitude =3D sqrt(PSD);
>
> % Compute phase of frequency spectrum
> SpectrumPhase =3D rand(size(PSD))*2.*pi;
>
> % Compute complex spectrum
> Spectrum =3D SpectrumAmplitude .* exp(i*SpectrumPhase);
>
> % This is a one-sided PSD, for frequencies between 0 and
> % the Nyquist frequency fNyq. We construct a two sided-PSD:
>
> % Sampling of frequency vector
> delta_f =3D fNyq/(N-1);
>
> % Generate two-sided frequency vector (as it would result
> % from Matlab's fft command for an even number of samples):
> f_TwoSided =3D ...
> =A0 =A0 [...
> =A0 =A0 [0 : delta_f : fNyq] ...
> =A0 =A0 [-fNyq + delta_f : delta_f : -delta_f]...
> =A0 =A0 ];
> % ... and generate the two-sided PSD:
> Spectrum_TwoSided =3D [Spectrum(1:1:end) Spectrum(end-1:-1:2)];
>
> % Compute inverse FFT
> x=3D ifft(Spectrum);
> x=3D fftshift(x);
> x=3D real(x);
>
> % Compute time vector
> delta_t =3D 1./(2.*fNyq);
> t=3D [0:delta:(N-1)*delta];
> =3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=
=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=
=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D
>
> What I am not sure about is:
> * if it is required to build a two-sided PSD prior to computing the inver=
se FFT with ifft.
> * In which ordering the time vector of the resulting time series x (from =
ifft) will be computed. Is fftshift required?
> * The problematic code seems to be:
>
> x=3D ifft(Spectrum);
> x=3D fftshift(x);
> x=3D real(x);
>
> In any case the function does not produce the results I expected. When I =
am trying, once x and t have been computed, to compute a periodogram using
>
> [Pxx,freq] =3D periodogram(x,hanning(length(x)),length(x),2*fNyq);
>
> the resulting periodogram should resemble to PSD but it does not (at all!=
).
>
> Any ideas what went wrong?
>
> Thanks
> Rainer

"Rainer " <wilhelm.remove.this.rainer@gmx.net> wrote in message <gao08m$ad3$1@fred.mathworks.com>...
> Hi,
>
> I tried to generate a time series (one realization) from a Power Spectral Density (PSD). Since the PSD only contains amplitude information but no phase information, I defined the phase (at each frequency) as a random variable, uniformly distributed between 0 and 2*pi. Basically the code is:
-------------------


If the idea is to get a realization of the time series, do you really need phase information? Do you really need to take this approach at all?

It seems like what you're trying to do is express the time series as the output of some filter when fed with white noise input. You're trying to recreate that filter and are finding that the PSD only gives you its magnitude response.

But if the goal is just to get a realization of the time series, it seems to me that all you have to do is compute the autocorrelation function by doing an ifft on the PSD.

Once you have the autocorrelation function, you can reconstruct the joint probability distribution function of the time series, if you assume the time series is Gaussian.

Once you have the joint distribution, you can compute a realization.