Power spectral density(PSD)

Power spectral density function (PSD) shows the strength of the variations(energy) as a function of frequency. In other words, it shows at which frequencies variations are strong and at which frequencies variations are weak. The unit of PSD is energy per frequency(width) and you can obtain energy within a specific frequency range by integrating PSD within that frequency range. Computation of PSD is done directly by the method called FFT or computing autocorrelation function and then transforming it.






The power spectral density, PSD, describes how the power (or variance) of a time series is distributed with frequency. Mathematically, it is defined as the Fourier Transform of the autocorrelation sequence of the time series.




PSD is a very useful tool if you want to identify oscillatory signals in your time series data and want to know their amplitude. For example let assume you are operating a factory with many machines and some of them have motors inside. You detect unwanted vibrations from somewhere. You might be able to get a clue to locate offending machines by looking at PSD which would give you frequencies of vibrations. PSD is still useful even if data do not contain any purely oscillatory signals. For example, if you have sales data from an ice-cream parlor, you can get rough estimate of summer sales peak by looking at PDF of your data. We quite often compute and plot PSD to get a "feel" of data at an early stage of time series analysis. Looking at PSD is like looking at simple time series plot except that we look at time series as a function of frequency instead of time. Here, we could say that frequency is a transformation of time and looking at variations in frequency domain is just another way to look at variations of time series data. PSD tells us at which frequency ranges variations are strong and that might be quite useful for further analysis.





PSD in matlab

Sampling frequency
Fs=50000;

%# of samples in the data
datasize=size(SINE512);
numsample=datasize(1);

%Windowing
H=hann(numsample);
W=H.*(SINE512(:,2));

%Fourier Transform
FFTX=fft(W,numsample);

%Power: magnitude^2
X=FFTX(1:floor(numsample/2)).*conj(FFTX(1:floor(numsample/2)));

%Bandwidth
BW=1.5*Fs*numsample;

%PSD=magnitude^2/bandwidth
PSD=X/BW;

%Computing the corresponding frequency values
Omega=Fs*(0:numsample-1)/numsample;
Omega=Omega(1:floor(numsample/2));

%Plot PSD
H=plot(Omega,PSD);
set(H,'Color','BLACK');