Time series from Power Spectral Density (PSD)

Hi, I would like to compute different time series (time realizations) of a given PSD. I've read there are several methods using random phase and amplitude or going through an IFFT. Is there any method already built in Matlab or how should this be done? I paste here below an example of a PSD. Thanks!

f = [6.5e-4 5.8e-2 5];
PSD = [1e2 1e-2 1e-6];

coef = polyfit(log10(f),log10(PSD),1);
f2 = linspace(0.8*1e-4,20,1000);
PSD2 = 10.^(polyval(coef,log10(f2)));

loglog(f2,PSD2);grid on;
xlabel('Hz');ylabel('m^2/Hz')

Answer 1

[William Rose]

@Albert Zurita, you wrote:

I would like to compute different time series (time realizations) of a given PSD. I've read there are several methods using random phase and amplitude or going through an IFFT.

You are partly correct. If you know the PSD you can compute the amplitude spectrum. You can use phases that are all zero, or that are random, or whatever you choose. Then you can compute th inverse FFT.

I assume the time sequence is to be a real sequence, not complex. If so, then the phase angle ust be zero or pi/2 at the Nyquist frequency (fNyq) (assuming an even number of samples; if the number is odd, then the Nyquist frequency will not be sampled...), and all complex components of hte spectrum above fNyq must be the complex conjugates of their "mirror images" below fNyq.

If you follow these guidelines then the ifft of the complex spectrum will yield a real sequence, as you probably desire.

[William Rose]

@Albert Zurita,

f = [6.5e-4 5.8e-2 5];
PSD = [1e2 1e-2 1e-6];
coef = polyfit(log10(f),log10(PSD),1);

Your frequency vector, f2, is not a frequency vector that could result from doing the FFT on a signal, because it does not start at zero. It should, because the inverse discrete Fourier transform implicitly assumes that the components of the complex spectrum are at equally spaced frequencies starting with f2=0. I realize that you chose your frequency vector, f2, so that you could do log10(f2) without getting an error. Let us make a frequency vector that is compatible with the forward and inverse discrete Fourier transform, in other words, a vector that satisfies th implicit assumptions of ifft(): f2 has 1001 values from f=0 to f=20, rather than 1000 values from 0.00008 to 20.

f2 = linspace(0,20,1001);

We will use your formula for the PSD, starting with the first non-zero frequency. We will assume the p.s.d.=0 at f2=0 (i.e. assume a time domain signal with zero mean value).

PSD2 = [0, 10.^(polyval(coef,log10(f2(2:end))))];

To convert the power spectral density to the equivalent power spectrum, multiply each element of PSD2 by , the frequency spacing. Also adjust by a factor of 2 (due to 1 sided versus 2 sided spectrum) and by N^2, where N=length of time domain signal, to get the normalization to come out right.

N=2*(length(f2)-1);
df=f2(3)-f2(2);
Xsq=PSD2*df*N^2/2;        % power spectrum, with factor

Now compute the amplitude spectrum:

Xabs=Xsq.^0.5;      % amplitude spectrum

Now we are ready to add phases. The phase should be 0 or pi/2 at f=0 and at the Nyquist frequency, f=fNyq. If they are not, then the time domain signal obtained by ifft() is likely to be complex. At all the other frequencies, let's make the phase be a random number between 0 and 2*pi:

Xphs=[0,2*pi*rand(1,length(f2)-2),0];

Compute the complex spectrum from the magntitude and phase:

X=Xabs.*exp(Xphs*1i);  % X is a vector of complex numbers

Now we need to add the "top half" of the spectrum, i.e. the values above the Nyquist frequency. This is done slightly differently depending on whether the spectrum so far has an even or odd number of elements.

if mod(length(f2),2)
    X=[X,conj(fliplr(X(2:end-1)))]; % do this if length(f) is odd
else                    
    X=[X,conj(fliplr(X(2:end)))];   % do this if length(f) is even
end

Now (finally!) we are ready to do the inverse transform:

y=ifft(X);

(We could have done x=ifft(X), but I don't like having x and X refer to two different things.)

Compute the corresponding time vector, and sampling rate:

fs=20*2;
dt=1/(N*df);    
t=(0:N-1)*dt;   % vector of times

Find PSD of x(t) by perioogram method:

[pyy,fp]=periodogram(y,[],N,fs);

Plot results:

figure; subplot(211); loglog(f2,PSD2,'ro',fp,pyy,'bx'); 
xlabel('Frequency (Hz)'); ylabel('(m^2/Hz)'); 
title('P.S.D.(x)'); grid on; 
legend('X=PSD','psd(y(t))')
subplot(212); plot(t,y,'-r.');
xlabel('Time (s)'); ylabel('x (m)'); 
title('Time Domain: x(t)'); grid on

If you want to try all phase angles equal to zero, the replace

Xphs=[0,2*pi*rand(1,length(f)-2),0];

with

Xphs=[0,2*pi*rand(1,length(f)-2),0];

I predict that this will result in a time domain signal with very large initial value, since every frequency sinusoid will be maximal at t=0. The time domain signal signal will also be very large at the end, since the inverse FFT results in a periodic signal, so the end is approximately equal to the beginning.

[Paul]

@William Rose

A few comments on this interesting example:

"The phase should be 0 or pi/2 at f=0 and at the Nyquist frequency"

I think that pi/2 should should really be pi.

It seems like the code is presupposing that f2(end) is the Nyquist frequency and that f2 is the "lower half" of the frequency vector. In this case, f2 must have an odd number of elements, which is true in this line

f2 = linspace(0,20,1001);

and N must be even, which is true in this line

N=2*(length(f2)-1);

But, then it seems the "if" statement here is not needed, as numel(X) must be odd (and the comments I'm sure meant to say length(f2)).

if mod(length(f2),2)
    X=[X,conj(fliplr(X(2:end-1)))]; % do this if length(f) is odd
else                    
    X=[X,conj(fliplr(X(2:end)))];   % do this if length(f) is even
end

If you want to try all phase angles equal to zero, the replace

Xphs=[0,2*pi*rand(1,length(f)-2),0];

with

Xphs=[0,2*pi*rand(1,length(f)-2),0];

Should f be replaced with f2, and should the second Xphs be multiplied by 0 (or otherwise made all zeros)?

[Albert Zurita]

Thanks for all the comments. I am quite flexible. I only have a PSD profile that I replicated. I don't have the 0 frequency in the profile because I'm not interested in the (time domain) bias terms. But I understand I need to create the full spectrum from -fs/2 to fs/2 going through 0, uniformly spaced. When I do the ifft I get NaNs in fact. I think this is due to the fact that the PSD is very steep...?

[William Rose]

[Edit: The psd of an even and of an odd length time dmain sequence can both have an odd nymber of elements. Therefore it is not possible to determine if the original time-domain sequence had an even or odd number of elements , if all you know is the PSD Therefore one can simply assume the original time domain sequence had an even number of elements. I should correct the code I posted ealrier, but din;t have time now. Sorry for the error.]

@Paul, thanks for the excellent comments and suggestions.

@Paul, you are right that the phase at 0 and the Nyquist frequency should be 0 or pi, not 0 or pi/2.

@Paul, you are right that I should have said f2, not f, in the comments.

@Paul, you are right that I did not implement all the if statements which I could have implemented for even versus odd number of elements. I made that choice for the sake of narrative simplicity.

@Albert Zurita, the power in your PSD is so small at 20 Hz (at and next to the Nyquist frequency), compared to the power at higher frequencies, that there will be no significant difference associated with how the last frequency in the PSD is treated, even versus odd. This is another reason I did not put in all the "if even, do this; if odd, do that" statements that I could have put in.

@Paul, since I specify Xphs before I compute the folded-over portion of the spectrum, it is OK to specify zero phase for the lower half of the spectrum only. The folded over part will then have zero phase too, as required.

@Albert Zurita, I am not sure why you get NaNs when you do ifft. I do not. You are right that the PSD rises sharply as frequency approaches zero. The spacing of your frequency samples is about 0.02 Hz, which implies a signal duration of 50 seconds. Therefore the lowest non-zero frequency in the PSD is 0.02 Hz. The power at this frequency is not a cause for NaNs.

[Albert Zurita]

First of all @William Rose and @Paul thank you for your very detailed analysis of the issue. I applied your comments and suggestions to the code.

I understand my PSD is log-like. I get your point on the 0 frequency, but at frequencies close to 0 the polyval is omitting values close to 0 because it spends a lot of points on the other frequencies. I think I should change the way I interpolate the initial PSD control points in [f, PSD] because there are missing points. For instance when doing

loglog(f,PSD,'*',f2,PSD2,'-');

There are gaps with respect to the original PSD. I think this is important because depending on the number of points selected in the linspace the approximation to the original PSD is better or worse, and this has a direct impact in the time-domain signal. I presume this is due to the fact that lower frequencies are more energetic, so changing the number of points impacts PSD2 and this directly impacts strongly the amplitude of the the time signal.

[Paul]

The case of zero-phase illustrates an interesting (at least to me) feature of the DFT, which is that a time domain signal that is real and conjugate symmetric has a transform that is also real and conjugate symmetric. And because of frequency/time domain symmetry, the opposite as true as well.

For example, using @William Rose's code to construct the zero-phase DFT

f2   = linspace(0,20,1001);
f    = [6.5e-4 5.8e-2 5];
PSD  = [1e2 1e-2 1e-6];
coef = polyfit(log10(f),log10(PSD),1);
PSD2 = [0, 10.^(polyval(coef,log10(f2(2:end))))];
N    = 2*(length(f2)-1);
df   = f2(3)-f2(2);
Xsq  = PSD2*df*N^2/2;   % power spectrum, with factor
Xabs = Xsq.^0.5;        % amplitude spectrum
X    = Xabs;            % X is a vector of real numbers

Create the conjugate symmetric DFT (X is real, so the conj part doesn't really matter here)

X    = [X , conj(fliplr(X(2:end-1)))];
y    = ifft(X);

As expected, y is real

max(abs(imag(y)))

ans = 0

And its plot shows the expected symmetry around the middle of the time vector

fs  = 20*2;
dt  = 1/(N*df);    
t   = (0:N-1)*dt;
figure
plot(t,y,'o')

This signal is symmetric around t = 0 after applying the appropriate shift

plot(t - N/2*dt,fftshift(y),'o')

[William Rose]

@Albert Zurita, I will address your quesiton in a separate comment.

@Paul, thank you for that insightful example.

The attached code generates the PSD of a sawtooth wave. Then it makes two signals with the same PSD: one with random phases and one with zero phase. The three time domain signals and the three PSDs are plotted - see below.

[William Rose]

@Albert Zurita,

You generate a vector of frequencies for your PSD by using linspace.

You then generate the corresponding PSD by fitting a straight line (a straight line on a log-log plot) to three points. The PSD you generate are points along this fitted straight line.

You write in your latest comment: "There are gaps with respect to the original PSD. I think this is important because depending on the number of points selected in the linspace the approximation to the original PSD is better or worse, and this has a direct impact in the time-domain signal."

You are right that the number of points selected wih linspace has a direct impact on the time domain signal. You probably know the following, but just in case you don't: The time domain signal duraiton is the reciprocal of the spacing of the frequencies. In your example, the frequency spacing was 0.02 Hz, so the time domain signal is 50 seconds long. The frequency spacing is also the lowest frequency that can exist in the reconstructed signal. So in your example, the lowest possible frequency that the reocnstructed signal can have is 0.02 Hz. If you are interested in lower freequencies, then you need finer frequency spacing in your PSD. Also, the highest frequency in the (one-sided) PSD is half the samping rate. So, in your example, the max frequency of the PSD is 20 Hz. This implies that the sampling rate is 40 Hz.

What exactly are you trying to do? Why do you fit a straight line to three points, instead of using the actual PSD? Do you have the actual PSD? How did you obtain it? I am asking because the answers might lead to a different approach that accomplishes your goals more effectively.

[William Rose]

@Albert Zurita,

You generate a vector of frequencies for your PSD by using linspace.

You then generate the corresponding PSD by fitting a straight line (a straight line on a log-log plot) to three points. The PSD you generate are points along this fitted straight line.

You write in your latest comment: "There are gaps with respect to the original PSD. I think this is important because depending on the number of points selected in the linspace the approximation to the original PSD is better or worse, and this has a direct impact in the time-domain signal."

You are right that the number of points selected wih linspace has a direct impact on the time domain signal. You probably know the following, but just in case you don't: The time domain signal duraiton is the reciprocal of the spacing of the frequencies. In your example, the frequency spacing was 0.02 Hz, so the time domain signal is 50 seconds long. The frequency spacing is also the lowest frequency that can exist in the reconstructed signal. So in your example, the lowest possible frequency that the reocnstructed signal can have is 0.02 Hz. If you are interested in lower freequencies, then you need finer frequency spacing in your PSD. Also, the highest frequency in the (one-sided) PSD is half the samping rate. So, in your example, the max frequency of the PSD is 20 Hz. This implies that the sampling rate is 40 Hz.

What exactly are you trying to do? Why do you fit a straight line to three points, instead of using the actual PSD? Do you have the actual PSD? How did you obtain it? I am asking because the answers might lead to a different approach that accomplishes your goals more effectively.

[Albert Zurita]

@William Rose thanks, I understand better the relationship. However, my original PSD is from now on a kind of 'mask PSD', and I only have an image file (png file) where I tried to obtain some points from (roughly). I now included a bit more points. But I find difficulties in doing the interpolation to have an uniformly sampled PSD to use the code. I think that's the reason behind the NaNs, it's just the interpolation problem (I changed now from polyfit to interp1 because of slightly different shape now at lower frequencies). With gaps I refer to the interpolation problem as you see in the plot below, not the fact that the frequencies are uniformly distributed. thanks!

f = [6e-5 1.1e-4 1.75e-4 7.5e-4 6.8e-2 6];
PSD = [0.85e3 1.8e3 0.8e3 1e2 1e-2 1e-6];
loglog(f,PSD,'b*');grid on;axis square;hold on;
f2 = linspace(0,20,10001);
PSD2 = 10.^(interp1(log10(f),log10(PSD),log10(f2)));
loglog(f2,PSD2);grid on;

[William Rose]

@Albert Zurita, your interpolated points look good. I see you have increased the number of points in linspace() from 1001 to 10001. Therefore the lowest non-zero frequency in PSD2 is now 0.002 Hz. I hope you have cured the NaN problem.

By the way, if you need to digitize more points from figures, I have written a Labview program to digitize points from figures. If you have Labview and if you are interested in the program, contact me by clickingon my name above, then click on the envelope in the pop-up window, to send me a private message. Or search the Matlab File Exchange, because I bet someone else has written a similar Matlab program.

[Albert Zurita]

Thanks for the offering. @William Rosev what I mean with the interpolation is that you'll have noticed that the interpolation is not good for the lower part of the spectrum in my last comment/plot, due to the log interpolation I think. This raises errors/NaNs.

[William Rose]

@Albert Zurita,

You get NaNs from interp1() because the "interpolation" is really extrapolation: you request values out to 20 Hz, but the input pairs only go as high as 6 Hz. The use of logs is not the issue. Consider the following example, in which the input frequencies range from 0 to 6, and the requested frequencies go from 0 to 20:

f=[0 6]; PSD=[10 1]; f2=linspace(0,20,6)

f2 = 1×6

0 4 8 12 16 20

PSD2=interp1(f,PSD,f2)

PSD2 = 1×6

10.0000 2.8000 NaN NaN NaN NaN

To eliminate the NaNs in PSD2, reduce the top frequency in f2 to 6 Hz, or increase the highest frequency in f to 20 Hz (and add a corresponding value to PSD).

[Albert Zurita]

However, the use of logs does a non-uniform interpolation, I should dig more on ways to interpolate straight lines in log scale...anyway thanks for the well explained answer.

[Bobby Anderson]

@William Rose

In this equation, where does the N^2 come from?

Xsq=PSD2*df*N^2/2; % power spectrum, with factor

Thanks!

[William Rose]

@Bobby Anderson,

Short answer: You need the factor N^2 due to the way the discrete Fourier transform (DFT) and the power spectral density (PSD) are normalized.

Long answer:

There is more than one way to normalize the DFT. Matlab and many other packages use a factor of 1 for the DFT and 1/N for the inverse DFT (iDFT). But you could use 1/N for the DFT and 1 for the iDFT. Or you could do 1/sqrt(N) for both the DFT and the iDFT.

Let x(n) (n=0..N-1) be the real time domain signal.

Let X(m) (m=0..N-1) be the DFT of x(n).

Let PSD(k) (k=0..N/2) be the one-sided* power spectral density of x(n), where N is even. (If N is odd, then k=0..(N-1)/2.)

Then

where , and where the factor 2 in the equaiton above becomes a 1, when k=0 and when k=N/2. It follows from the equations above that

which is why my code includes the line

Xsq=PSD2*df*N^2/2;        % power spectrum

See formulas at More About Periodogram and at More About FFT.

* Matlab's periodogram returns the one-sided PSD by default, if x(n) is real.

[Paul]

@William Rose

Regarding the one-sided vs. two-sided ....

As you and periodogram state, the factor of two is not applied for k = 1 and k corresponding to Nyquist "so that the total power is conserved" (according to the doc page).

Consider:

rng(100);
n = 0:319;
x = cos(pi/4*n)+randn(size(n)) + 20; % signal chosen with large mean relative to the rest of the signal
[pxx1,w1] = periodogram(x);
[pxx2,w2] = periodogram(x,'centered');
figure
plot(w1,pxx1,'-x',w2,pxx2,'-o'),grid,xlim([-0.1 0.1]);

Because periodogram does not muliply pxx1[0] by 2, the area under pxx1 between pxx1[0] and pxx1[1] is not twice the area under pxx2 between pxx2[0] and pxx2[1]. To wit

[sum(abs(x).^2)/numel(x), trapz(w1,pxx1), trapz(w2,pxx2)]

ans = 1×3

401.0130 276.1640 401.0113

<mw-icon class=""></mw-icon>

<mw-icon class=""></mw-icon>

Note that w2 doesn't start exactly at -pi, so we can sharpen its power estimate

w2(1)

ans = -3.1293

trapz([-pi;w2],[pxx2(end);pxx2])

ans = 401.0130

Similarly, we have to be careful if using 'twosided' to deal with the last point

[pxx2,w2] = periodogram(x,'twosided');
trapz(w2,pxx2)

ans = 248.3623

trapz([w2;2*pi],[pxx2;pxx2(1)])

ans = 401.0130

Disclaimer: I'm not as familiar with PSD stuff as I'd like to be, so I may be misinterpreting the output of periodogram.

[William Rose]

@Paul,

Thank you for those illustrations and for the reminder that the applicaiton of Parseval's theorem is a tricky business, especially when comparing one- and two-sided spectra. I always learn from your posts.

[Paul]

Note that if we use rectangular integration instead of trapezoidal, then we don't need to worry about those endpoints, a feature of discrete time processing about which I've always been amazed.

rng(100);
n = 0:319; N = numel(n);
x = cos(pi/4*n)+randn(size(n)) + 20; % signal chosen with large mean relative to the rest of the signal
[pxx1,w1] = periodogram(x);
[pxx2,w2] = periodogram(x,'centered');
[pxx2a,w2a] = periodogram(x,'twosided');
N1 = numel(w1);
N2 = numel(w2);
[sum(abs(x).^2)/N, sum(abs(fft(x)).^2)/N^2, 2*pi*N2*sum(pxx2)/N2^2, 2*pi*N2*sum(pxx2a)/N2^2, 2*pi*N2*sum(pxx1)/N2^2]

ans = 1×5

401.0130 401.0130 401.0130 401.0130 401.0130

<mw-icon class=""></mw-icon>

<mw-icon class=""></mw-icon>

[William Rose]

@Paul, nice. I didn't anticipate the factors of 2*pi in the three rightmost formulas for the total power, in your preceding post. That is because I'm not used to calling periodogram without an argument for the sampling rate, fs. The factor of 2*pi is consistent with the formulas here, for the case when "the frequencies are in radians/sample" (i.e. when fs is not supplied).

[Paul]

@William Rose

The reason why I'm amazed by Parseval's relationship for the DFT is based on a visualization of the quantities inolved.

Define a finite duration signal of length 20

rng(101);
N = 20;
n = 0:N-1;
x = rand(N,1); x = x - mean(x)*0.6;

Its DFT

Xdft = fftshift(fft(x));
wdft = (-N/2:(N/2-1))/N*2*pi;

Its DTFT

XDTFT = @(w) freqz(x,1,w);
wDTFT = linspace(-1,1,2048)*pi;

Overlay the abs(DFT)^2 on the abs(DTFT)^2

figure
plot(wDTFT,abs(XDTFT(wDTFT)).^2,'LineWidth',2);
hold on
plot(wdft,abs(Xdft).^2,'o');

Add the rectangles of height abs(DFT)^2 and width dw = 2*pi/N (rad/sample)

bar(wdft,abs(Xdft).^2,1,'FaceAlpha',0.5);
xlabel('\omega  (rad/sample)')

The amazing thing to me is that the area under the DTFT curve is the same as the area under the DFT rectangles. To me, at least, that is a very unexpected result by just looking at the plot. And there is some funny business at the edges where the left most rectangle extends to less than -pi and the right most rectangle stops before pi, and becasue N is even we have he unpaired point on the far left. Nevertheless, we have

[integral(@(w) abs(XDTFT(w)).^2,-pi,pi), 2*pi/N*sum(abs(Xdft).^2)]

ans = 1×2

14.7415 14.7415

<mw-icon class=""></mw-icon>

<mw-icon class=""></mw-icon>

which is, of course, 2*pi*E(x[n]), where E(x[n]) is the energy in x[n], exactly as Parseval tells us.

sum(abs(x).^2)*2*pi

ans = 14.7415

Intuitively, I can see why the DTFT works the way it does wrt to E(x[n]), but the DFT is another story for me. Normally, the rectangular method is an approximation to an integral, but here it is exact, despite apparent visual appearances to the contrary.

[William Rose]

@Paul, Well said and very nicely illustrated!