MATLAB Answers


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

Asked by Tom
on 13 Jul 2012

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




No products are associated with this question.

2 Answers

Answer by Wayne King
on 13 Jul 2012
Edited by Wayne King
on 13 Jul 2012
 Accepted 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

  1 Comment

on 15 Jul 2012

Thank you Wayne!

Answer by Wayne King
on 13 Jul 2012
Edited by Wayne King
on 13 Jul 2012

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;
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);

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

  1 Comment

on 13 Jul 2012

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?

Discover MakerZone

MATLAB and Simulink resources for Arduino, LEGO, and Raspberry Pi

Learn more

Discover what MATLAB® can do for your career.

Opportunities for recent engineering grads.

Apply Today

MATLAB Academy

New to MATLAB?

Learn MATLAB today!