Path: news.mathworks.com!not-for-mail
From: <HIDDEN>
Newsgroups: comp.soft-sys.matlab
Subject: Re: Generate a time series from a PSD
Date: Tue, 16 Sep 2008 11:22:02 +0000 (UTC)
Organization: The MathWorks, Inc.
Lines: 72
Message-ID: <gao4sq$3bf$1@fred.mathworks.com>
References: <gao08m$ad3$1@fred.mathworks.com>
Reply-To: <HIDDEN>
NNTP-Posting-Host: webapp-02-blr.mathworks.com
Content-Type: text/plain; charset="ISO-8859-1"
Content-Transfer-Encoding: 8bit
X-Trace: fred.mathworks.com 1221564122 3439 172.30.248.37 (16 Sep 2008 11:22:02 GMT)
X-Complaints-To: news@mathworks.com
NNTP-Posting-Date: Tue, 16 Sep 2008 11:22:02 +0000 (UTC)
X-Newsreader: MATLAB Central Newsreader 791003
Xref: news.mathworks.com comp.soft-sys.matlab:490445



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