Time series from Power Spectral Density (PSD)
86 次查看(过去 30 天)
显示 更早的评论
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')
![](https://www.mathworks.com/matlabcentral/answers/uploaded_files/1401604/image.png)
采纳的回答
William Rose
2023-6-3
@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 个评论
William Rose
2023-6-3
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.
![](https://www.mathworks.com/matlabcentral/answers/uploaded_files/1401774/image.png)
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
![](https://www.mathworks.com/matlabcentral/answers/uploaded_files/1401779/image.png)
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
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
2023-6-3
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
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.
Albert Zurita
2023-6-4
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
2023-6-4
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.
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')
![](https://www.mathworks.com/matlabcentral/answers/uploaded_files/1402614/image.png)
This signal is symmetric around t = 0 after applying the appropriate shift
plot(t - N/2*dt,fftshift(y),'o')
![](https://www.mathworks.com/matlabcentral/answers/uploaded_files/1402619/image.png)
William Rose
2023-6-5
@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.
![](https://www.mathworks.com/matlabcentral/answers/uploaded_files/1402829/image.jpeg)
William Rose
2023-6-5
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
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
2023-6-5
@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;
![](https://www.mathworks.com/matlabcentral/answers/uploaded_files/1402924/image.png)
William Rose
2023-6-5
@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
2023-6-5
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
2023-6-6
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.
更多回答(0 个)
另请参阅
类别
在 Help Center 和 File Exchange 中查找有关 Parametric Spectral Estimation 的更多信息
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!发生错误
由于页面发生更改,无法完成操作。请重新加载页面以查看其更新后的状态。
您也可以从以下列表中选择网站:
如何获得最佳网站性能
选择中国网站(中文或英文)以获得最佳网站性能。其他 MathWorks 国家/地区网站并未针对您所在位置的访问进行优化。
美洲
- América Latina (Español)
- Canada (English)
- United States (English)
欧洲
- Belgium (English)
- Denmark (English)
- Deutschland (Deutsch)
- España (Español)
- Finland (English)
- France (Français)
- Ireland (English)
- Italia (Italiano)
- Luxembourg (English)
- Netherlands (English)
- Norway (English)
- Österreich (Deutsch)
- Portugal (English)
- Sweden (English)
- Switzerland
- United Kingdom (English)
亚太
- Australia (English)
- India (English)
- New Zealand (English)
- 中国
- 日本Japanese (日本語)
- 한국Korean (한국어)