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

评论： Albert Zurita 2023-6-6

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.

14 个评论
显示 12更早的评论隐藏 12更早的评论

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

在 MATLAB Online 中打开

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

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.

请先登录，再进行评论。