Could you please let me know how to write code to draw the power spectral density (PSD) vs frequency for my data? Thanks.

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

Answer 1

[Vilnis Liepins]

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.

[Ray Huang]

I think that it is the excellent answer. Let me try the first way. However, I'm also interested in the second way. Thank you very much for help.

[Vilnis Liepins]

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)

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.