Thursday 26 December 2013

Power Spectral Density Estimator

The estimation of the Power Spectral Density (PSD) of a signal is of fundamental importance in Electrical Engineering, in particular Digital Communications.

This blog post implements a Power Spectral Density (PSD) Estimator using Welch's modification of the averaged periodogram estimate method. Some testing has been done, and the results are in agreement with those produced by GNU Octave (a free Matlab-clone). The averaged periodogram estimate method is also known as Bartlett's method.

Using Welch's method, the complex time domain data input samples are broken into contiguous (possibly overlapping) frames of size $N$, which has to be an integer power of 2 (as required for implementation of the Fast Fourier Transform (FFT)). Each of these frames is then weighted by the appropriate windowing function (see drop down menu below), and the $N$-point FFT is applied to the windowed signals within a frame. The magnitude squared of each resulting FFT bin is calculated for each frame, and time-averaged over all the frames to yield the PSD estimate.

Recall that for discrete time-domain input samples $x[n]$ for $n={0,1,2,..,N-1}$ the FFT (at bin $k$ for $k={0,1,2,..,N-1}$) is defined by equation

$X[k]=\sum_{n=0}^{N-1}x[n]exp(-j2\pi\frac{nk}{N})$

where $N=2^g$ for integer $g$.

Please enter the numbers in the text areas below - one number per line, for each of the Real and Imaginary input textareas (the textareas have already been filled in with some numbers for illustration purposes). There must be no new line after the last number.

Note that if the input is real only, the imaginary input textarea can be left empty (rather than having to fill it with the same number of zeros as there are real inputs, which can be a bit more cumbersome). Conversely, if the input is imaginary only, the real input textarea can be left empty (rather than having to fill it with the same number of zeros as there are imaginary inputs).

Alternatively you can choose to load a CSV file, which must be either a single column of numbers (for a real only input) or two comma-separated columns of numbers - the first line can be a comment line, starting with the character #.

To set the number of FFT bins for spectral analysis, fill in the number of bins field. It is possible to have overlapping frames - to set the number of overlapping samples between adjacent frames, fill in the relevant field below. A value of 0 means that there is no overlap between adjacent frames.

There is the option of specifying a sampling frequency below - setting this to a non-zero value will result in the x-axes of the PSD plots being set to the frequency as opposed to the frequency bin index. Leaving this field blank will result in the x-axes of the plots being set to the FFT bin index.

To select the appropriate windowing function, use the drop down menu below. Furthermore, you can select between linear units and dB units, using the checkbox below. For each bin, the dB is calculated as $10\times log_{10}(|FFT|^2)$ where $|FFT|^2$ is the magnitude squared value of the FFT for that bin.

There is a second checkbox which allows you to center the zero frequency (i.e. DC component) of the spectrum. If the DC component is centred then the frequency indices that exceed half the sampling frequency will be set to negative values.

To perform the PSD, press the button labelled "Calculate PSD Estimate" below - the results will populate the textarea below labelled "Power Spectral Density", with two comma separated columns - the first column being the frequency (in bin index or Hz) and the second column being the PSD estimate (in dB or linear units). The first line of this textarea will be a comment line starting with character #, and describing the two columns. One can copy and paste the results from this textarea to a text file if so desired. In addition, graphs of the output PSD will be plotted.

Once you have pressed the "Calculate PSD Estimate" button to display the PSD, you can zoom in onto the plot. Note that there are two graphs representing the same data - the lower smaller graph is an overview graph representing the data in its entirety, while the larger graph is the main graph which will represent the zoomed in subsection. You can zoom by selecting sections from either the overview graph or the main graph. To completely undo the zooming operations, simply press the "Calculate PSD Estimate" button again.

If you change inputs to a smaller number of samples, press the calculate button twice for the results to take effect. Alternatively, you can simply reload the page, then fill in the input textareas.



Real Input
Imaginary Input



Number of data samples is ..


Enter number of FFT bins for PSD:-


Enter number of overlap samples (must be between 0 and half the number of FFT bins inclusive):-


Enter sampling frequency (optional):-





Power Spectral Density:-
PSD plot

Linear

FFT bin

Linear

FFT bin



Extra information..

1 comment: