MATLAB: How to create power spectral density from fft (fourier transform)

fftfourierpower spectral density

Hi All.
Apologies if this is a basic post! – I am by no means a mathematician (my background is in biomechanics).
I would like to use MATLAB to plot power spectral density of force platforms traces from various impacts. Using the fft function, so far I have this (where x is my signal):
Fs = 500; % Sampling frequency T = 1/Fs; % Sample time L = 4000; % Length of signal t = (0:L-1)*T; % Time vector
NFFT = 2^nextpow2(L); % Next power of 2 from length of y X = fft(x,NFFT)/L; f = Fs/2*linspace(0,1,NFFT/2+1); AMP = 2*abs(X(1:NFFT/2+1));
This gives me the absolute value of my transform. As I understand it is ‘per unit bin’, so could be plotted against bin number on the x axis.
To get power spectral density do I simply need to square this and normalise to frequency? Thus:
%PSDx = AMP.^2/L; % Power Spectral Density
And total power of the spectrum is:
Tend = T*L; Pxf = Tend*sum(PSDx); % Total spectral power
I’m not sure of these last 2 sections. If they are correct, and my initial signal was in N, does that mean my power spectrum is in N^2/Hz?
Thanks in advance
Tom

Best Answer

The 1/L comes from the fact that you are using a "biased" estimate of the autocorrelation function to produce the PSD estimate. Think of taking the sample mean, you divide by the number of elements. While 1/L would be unbiased in the case of the mean, it turns out to be biased in the estimation of the autocorrelation, but there are good reasons to use the biased estimate. Remember, the autocorrelation is an average.
The 1/Fs, which is equal to the sampling interval comes from the Fourier transform of an infinite discrete-time sequence with a given sampling interval (the elements of the sequence are spaced by 1/Fs), so discrete-time/continuous frequency. Again, the sequence would be the autocorrelation sequence.
Hope that helps, Wayne
Related Question