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.
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
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.