PSD - Is this formula correct?

Trying to understand whether or not I am calculating the PSD correctly. Below is the equation that I am using:

$$PSD_{rms}(f_{m} = m \cdot f_{res}) = \frac{PS_{rms} }{ENBW} = \frac{2 \cdot \left | y_{m} \right |^{2}}{f_{s}\cdot S_{2}} m = 0....N/2$$

This equation is applied using the following:

1) Calculate the squared magnitude of the FFT bin:

$mag = (re*re+im*im)$

2) For each of these magnitude, I calculate the following:

$PSD(0...NFFT/2) = 2 \cdot (mag[0... NFFT/2] / Fs \cdot S2$

So it looks like the following in code:

PSD[i] = 2 * (result[i]) / (12000 * S2);


I get the following result:

This looks like the result is being scaled correctly. However, the removal of the noise...

I have been told to use a band-pass filter, but, looking at the source code of Matlab in their spectrogram algorithm they do not pass the results through any kind of filtering algorithm.

My question is: Should the output be this, and, if so, should I therefore apply a filtering algorithm even though matlab does not?

Look at matlab result:

• I think I misunderstood your previous question. I thought you wanted to keep signal content in a narrow range. That's when you would use a band pass filter. If you just want to display a PSD, you're basically done. If you want to remove noise, you have to specify what you consider noise. Perhaps another image would be helpful. – Phonon Jan 24 '14 at 17:54
• @Phonon - Please look at my update. Notice in the plot (using the same signal, same NFFT and noverlap) I get that plot from matplotlib.. There is no noise.. Why am I getting such different results? – Phorce Jan 24 '14 at 17:57
• What do the different color curves represent? – Phonon Jan 24 '14 at 17:59
• @Phonon Where do you mean, sorry? This is a plot of a spectrogram, specifically returning the PSD of the results.. So I believe they represent the PSD – Phorce Jan 24 '14 at 18:00
• @Phonon - Sorry, I just can't seem to get my head around it. I've re-implemented this over and over again, read so many papers also but nothing seems to make sense why the two results are different – Phorce Jan 24 '14 at 18:02

My question is: Should the output be this, and, if so, should I therefore apply a filtering algorithm even though matlab does not?

Without seeing your code I can't determine what's the reason for a difference in the two plots. There's multiple things that come into play like window type, overlap, and FFT size.

I can't say for sure if you should filter your data. It depends on the application. If you're just going to plot it, it may not be necessary. If the out of band energy or noise is really high and leaking into the signal band, then maybe. If this is going to be the input to some other process and you think the out of band energy is degrading some results, then sure, go ahead and filter. Without knowing your application it's tough to say what's right or wrong.

Anyway, in order to filter your signal to remove out of band noise, there's two approaches, but the easiest would be to do it in the time-domain.

Filter the data in the time-domain with either an FIR or IIR filter. Have a look at functions fir1 and fir2 in MATLAB. Once you design your filter, use freqz to quickly plot the frequency response to make sure it's correct. Once you have your filter 'b', compute:

% design filter 'b' using fir1 or fir2
y = filter(b, 1, x); % Filter input x with filter b
% compute PSD of y


Depending on your filter you should see suppression of energy in the stop bands.

• @proten - Hey, I've since sorted this :) The problem was that I was reading the .wav file incorrectly and now feel like a complete idiot!! Thanks anyway – Phorce Jan 26 '14 at 13:57