MATLAB Central Newsreader - Generate a time series from a PSD

Generate a time series from a PSD

Rainer — Tue, 16 Sep 2008 10:03:02 -0400

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

Re: Generate a time series from a PSD

David — Tue, 16 Sep 2008 11:22:02 -0400

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

Re: Generate a time series from a PSD

NZTideMan — Tue, 16 Sep 2008 20:07:23 -0400

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.

Re: Generate a time series from a PSD

SmartEngineer — Wed, 24 Sep 2008 13:05:05 -0400

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.

Re: Generate a time series from a PSD

Greg Heath — Sat, 27 Sep 2008 08:16:45 -0400

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

Re: Generate a time series from a PSD

Matt — Sat, 27 Sep 2008 13:30:05 -0400