<pre>

Path: news.mathworks.com!not-for-mail
From: &lt;HIDDEN&gt;
Newsgroups: comp.soft-sys.matlab
Subject: Generate a time series from a PSD
Date: Tue, 16 Sep 2008 10:03:02 +0000 (UTC)
Organization: EADS Deutschland GmbH
Lines: 67
Message-ID: &lt;gao08m$ad3$1@fred.mathworks.com&gt;
Reply-To: &lt;HIDDEN&gt;
NNTP-Posting-Host: webapp-02-blr.mathworks.com
Content-Type: text/plain; charset=&quot;ISO-8859-1&quot;
Content-Transfer-Encoding: 8bit
X-Trace: fred.mathworks.com 1221559382 10659 172.30.248.37 (16 Sep 2008 10:03:02 GMT)
X-Complaints-To: news@mathworks.com
NNTP-Posting-Date: Tue, 16 Sep 2008 10:03:02 +0000 (UTC)
X-Newsreader: MATLAB Central Newsreader 1502580
Xref: news.mathworks.com comp.soft-sys.matlab:490441


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

<pre>