How to create power spectral density from fft (fourier transform)
127 views (last 30 days)
Show older comments
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
0 Comments
Accepted Answer
Wayne King
on 13 Jul 2012
Edited: Wayne King
on 13 Jul 2012
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
More Answers (1)
Wayne King
on 13 Jul 2012
Edited: 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;
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.
2 Comments
Shiva
on 29 Feb 2024
@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
See Also
Categories
Find more on Spectral Measurements in Help Center and File Exchange
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!