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

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

Answer 1

[Wayne King]

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

[Tom]

Thank you Wayne!

Answer 2

[Wayne King]

The periodogram is:

Fs = 500; % Sampling frequency 
T = 1/Fs; % Sample time 
L = 4000; % Length of signal 
t = (0:L-1)*T; % Time vector
x = cos(2*pi*100*t)+randn(size(t));
xdft = fft(x);
Pxx = 1/(L*Fs)*abs(xdft(1:length(x)/2+1)).^2;
freq = 0:Fs/L:Fs/2;
plot(freq,10*log10(Pxx));
xlabel('Hz'); ylabel('dB/Hz');

I don't think you need to zero pad and I've created a test signal, x, in the above so you can copy and paste the code to generate a PSD estimate.

Now, MATLAB scales the single-sided PSD estimate by two for all frequencies except 0 (DC) and Fs/2 (the Nyquist), but this is a convention aimed at conserving total power and is NOT part of the definition of the periodogram.

So, if you have the Signal Processing Toolbox and you want to get perfect agreement with MATLAB's periodogram.m, you can do:

    Pxx(2:end-1) = 2*Pxx(2:end-1);
    plot(freq,10*log10(Pxx));

And yes, the units are in N^2/Hz, i.e.

[Tom]

Thank you that's great. Could you explain to me why I should multiply the absolute value of the transform by 1/(L*Fs)? I'm confused as to why this normalises the signal to Hz?

[Shiva]

@Wayne King, @Tom i did the same thing as you explained, but when i calculated the gRMS of the PSD curve it is not equal or near the gRMS of the original input curve(x).

We need the PSD to reflect the rms of the original curve right? If wrong please tell me.

And how can we convert the PSD curve into Octave band like 1/3,1/12 etc..

Thanks in advance