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

Hi All,

My first post on this forum – apologies if it’s a basic one! - 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 (AMP). 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