Path: news.mathworks.com!newsfeed-00.mathworks.com!newsfeed2.dallas1.level3.net!news.level3.com!postnews.google.com!e53g2000hsa.googlegroups.com!not-for-mail
From: Greg Heath <heath@alumni.brown.edu>
Newsgroups: comp.soft-sys.matlab
Subject: Re: Generate a time series from a PSD
Date: Sat, 27 Sep 2008 01:16:45 -0700 (PDT)
Organization: http://groups.google.com
Lines: 91
Message-ID: <36d2a9f2-dcc6-49ca-9aeb-16d90fc5fbc1@e53g2000hsa.googlegroups.com>
References: <gao08m$ad3$1@fred.mathworks.com>
NNTP-Posting-Host: 172.132.80.159
Mime-Version: 1.0
Content-Type: text/plain; charset=ISO-8859-1
Content-Transfer-Encoding: quoted-printable
X-Trace: posting.google.com 1222503405 17997 127.0.0.1 (27 Sep 2008 08:16:45 GMT)
X-Complaints-To: groups-abuse@google.com
NNTP-Posting-Date: Sat, 27 Sep 2008 08:16:45 +0000 (UTC)
Complaints-To: groups-abuse@google.com
Injection-Info: e53g2000hsa.googlegroups.com; posting-host=172.132.80.159; 
	posting-account=mUealwkAAACvQrLWvunjg50tRAnsNtJR
User-Agent: G2/1.0
X-HTTP-UserAgent: Mozilla/4.0 (compatible; MSIE 7.0; Windows NT 
	5.1),gzip(gfe),gzip(gfe)
Xref: news.mathworks.com comp.soft-sys.matlab:492402


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