Could you please let me know how to write code to draw the power spectral density (PSD) vs frequency for my data? Thanks.
3 次查看(过去 30 天)
显示 更早的评论
Time(s) Amplitude 0.707 0.707 1.401 0.694 2.104 0.703 2.803 0.699 3.497 0.694 4.183 0.686 4.869 0.686 5.575 0.706 6.288 0.713 7.007 0.719 7.711 0.704 8.401 0.69 9.098 0.697 9.8 0.702 10.505 0.705 11.192 0.687 11.876 0.684
0 个评论
采纳的回答
Vilnis Liepins
2018-5-23
Hi Ray,
Your data are:
amplitude=[0.707 0.694 0.703 0.699 0.694 0.686 0.686 0.706 0.713 0.719 0.704 0.69 0.697 0.702 0.705 0.687 0.684];
t=[0.707 1.401 2.104 2.803 3.497 4.183 4.869 5.575 6.288 7.007 7.711 8.401 9.098 9.8 10.505 11.192 11.876];
The vector t values are non-uniformly spaced in time. I see two options how to calculate PSD:
1) You can use the Matlab function 'interp1' to interpolate data on uniform time grid and then apply the function 'fft' or any other Matlab function that calculates PSD. Code example:
T_mean=(t(length(t))-t(1))/(length(t)-1); %Calculate mean sampling period
amplitude_i=interp1(t,amplitude,t(1):T_mean:t(length(t))); %Interpolate data on uniform time grid
Psd_i=abs(fft(amplitude_i)).^2/(length(t)*T_mean); %Calculate PSD of interpolated data
semilogy(Psd_i); %Plot PSD on log scale as PSD has a large peak on zero frequency
2) Another option is to use the function 'nedft' downloadable on fileexchange:
https://se.mathworks.com/matlabcentral/fileexchange/11020-extended-dft
This function allows you to calculate the PSD directly from the non-uniform data (amplitude and t vectors) without interpolating them to uniform grid. Let me known if you are interested in applying it.
2 个评论
Vilnis Liepins
2018-5-26
编辑:Vilnis Liepins
2018-5-28
Please see the code bellow where the FFT of interpolated uniform data is compared with the NEDFT of non uniform data.
amplitude=[0.707 0.694 0.703 0.699 0.694 0.686 0.686 0.706 0.713 0.719 0.704 0.69 0.697 0.702 0.705 0.687 0.684];
t=[0.707 1.401 2.104 2.803 3.497 4.183 4.869 5.575 6.288 7.007 7.711 8.401 9.098 9.8 10.505 11.192 11.876];
NDFT=length(amplitude); %The length of the DFT equal to length of the data
T_mean=(t(NDFT)-t(1))/(NDFT-1); %Calculate mean sampling period
t_u=t(1):T_mean:t(NDFT); %Uniform time grid, sampling period T_mean
amplitude_i=interp1(t,amplitude,t_u,'spline'); %Interpolate the data onto uniform time grid by 'spine' method
Psd_i=abs(fft(amplitude_i)).^2/(NDFT*T_mean); %Calculate the PSD of interpolated data by FFT algorithm
fr=[-ceil((NDFT-1)/2):floor((NDFT-1)/2)]/(NDFT*T_mean); %Calculate frequencies vector in range ]-1/(2T_mean),1/(2T_mean)]
Psd_a=abs(nedft(amplitude,t,fr)).^2/(NDFT*T_mean); %Calculate the PSD of nonuniform data (amplitude, t) by NEDFT algorithm
semilogy(fr,fftshift(Psd_i),'-b', fr, Psd_a, '-g'); %Plot the PSD by frequency calculated by FFT (blue) and NEDFT (green)
figure
NDFT=1000; %Increase the length of the DFT
Psd_i=abs(fft(amplitude_i,NDFT)).^2/(length(t)*T_mean); %Calculate PSD of interpolated data by FFT algorithm
fr=[-ceil((NDFT-1)/2):floor((NDFT-1)/2)]/(NDFT*T_mean); %Calculate frequencies vector in range ]-1/(2T_mean),1/(2T_mean)]
Psd_a=abs(nedft(amplitude,t,fr)).^2/(NDFT*T_mean); %Calculate PSD of nonuniform data (amplitude, t) by NEDFT algorithm
semilogy(fr,fftshift(Psd_i),'-b', fr, Psd_a, '-g'); %Plot the PSD by frequency calculated by FFT (blue) and NEDFT (green)
![](https://www.mathworks.com/matlabcentral/answers/uploaded_files/188997/image.png)
![](https://www.mathworks.com/matlabcentral/answers/uploaded_files/188999/image.png)
As you see from the plots if the length of FFT/NEDFT is the same as data (17 samples) then PSD estimates almost coincide (Figure 1). This is because of your data are very close to being uniformly sampled - check the difference (t-t_u).
The increasing of the length of the DFT up to 1000 allows NEDFT algorithm refine the PSD estimate (Figure 2, green plot) while the FFT in this case (Figure 2, blue plot) shows even worse picture - the high power of zero frequency masking the week data components (spectral leakage).
This, in short, on the application of NEDFT algorithm for your data.
更多回答(0 个)
另请参阅
类别
在 Help Center 和 File Exchange 中查找有关 Fourier Analysis and Filtering 的更多信息
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!