Time series from Power Spectral Density (PSD)

Question

Albert Zurita 2023-6-2

0
链接

此问题的直接链接

https://ww2.mathworks.cn/matlabcentral/answers/1977509-time-series-from-power-spectral-density-psd

评论： William Rose 2024-10-30，5:18

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')

0 个评论
显示 -2更早的评论隐藏 -2更早的评论

请先登录，再进行评论。

请先登录，再回答此问题。

Answer 1

William Rose 2023-6-3

1
链接

此回答的直接链接

https://ww2.mathworks.cn/matlabcentral/answers/1977509-time-series-from-power-spectral-density-psd#answer_1249464

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

22 个评论
显示 20更早的评论隐藏 20更早的评论

William Rose 2023-6-3

在 MATLAB Online 中打开

@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 2023-6-3

在 MATLAB Online 中打开

@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)?

William Rose 2023-6-4

编辑：William Rose 2023-6-4

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

Paul 2023-6-4

在 MATLAB Online 中打开

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 2023-6-5

@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 2023-6-5

@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 2024-10-29，5:15

编辑：William Rose 2024-10-29，5:17

在 MATLAB Online 中打开

@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 2024-10-29，12:25

在 MATLAB Online 中打开

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

Paul 2024-10-30，0:49

编辑：Paul 2024-10-30，2:51

在 MATLAB Online 中打开

@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 2024-10-30，5:18

@Paul, Well said and very nicely illustrated!

请先登录，再进行评论。

Time series from Power Spectral Density (PSD)

0 个评论
显示 -2更早的评论隐藏 -2更早的评论

采纳的回答

22 个评论
显示 20更早的评论隐藏 20更早的评论

更多回答（0 个）

另请参阅

类别

标签

产品

版本

Community Treasure Hunt

Time series from Power Spectral Density (PSD)

0 个评论 显示 -2更早的评论隐藏 -2更早的评论

采纳的回答

22 个评论 显示 20更早的评论隐藏 20更早的评论

更多回答（0 个）

另请参阅

类别

标签

产品

版本

Community Treasure Hunt

0 个评论
显示 -2更早的评论隐藏 -2更早的评论

22 个评论
显示 20更早的评论隐藏 20更早的评论