Path: news.mathworks.com!newsfeed-00.mathworks.com!newsfeed2.dallas1.level3.net!news.level3.com!postnews.google.com!w24g2000prd.googlegroups.com!not-for-mail
From: NZTideMan <mulgor@gmail.com>
Newsgroups: comp.soft-sys.matlab
Subject: Re: Generate a time series from a PSD
Date: Tue, 16 Sep 2008 13:07:23 -0700 (PDT)
Organization: http://groups.google.com
Lines: 97
Message-ID: <bc4f9947-a80e-41a6-889d-296f82f447eb@w24g2000prd.googlegroups.com>
References: <gao08m$ad3$1@fred.mathworks.com>
NNTP-Posting-Host: 202.78.152.105
Mime-Version: 1.0
Content-Type: text/plain; charset=ISO-8859-1
Content-Transfer-Encoding: quoted-printable
X-Trace: posting.google.com 1221595643 23925 127.0.0.1 (16 Sep 2008 20:07:23 GMT)
X-Complaints-To: groups-abuse@google.com
NNTP-Posting-Date: Tue, 16 Sep 2008 20:07:23 +0000 (UTC)
Complaints-To: groups-abuse@google.com
Injection-Info: w24g2000prd.googlegroups.com; posting-host=202.78.152.105; 
	posting-account=qPexFwkAAABOl8VUndE6Jm-9Z5z_fSpR
User-Agent: G2/1.0
X-HTTP-UserAgent: Mozilla/4.0 (compatible; MSIE 7.0; Windows NT 5.1; Maxthon; 
	Mozilla/4.0 (compatible; MSIE 6.0; Windows NT 5.1; SV1) ; .NET CLR 
	1.1.4322),gzip(gfe),gzip(gfe)
X-HTTP-Via: 1.1 bc3
Xref: news.mathworks.com comp.soft-sys.matlab:490521


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.