Power Spectral Density estimate (PSD)

Question

I'm trying to understand PSD, and, it's functions.

So far I have carried out the following:

• Segmented a signal with NFFT = 256 and hop_size = 128
• Windowed each of these segments (Hamming)
• DFT of each Window
• Taken the first half of the resulting matrix (n/2+1)

From what I have read so far, the algorithm can be best described as the following:

We take the squared magnitude of the LENGTH NFFT, we calculate this by: mag = (bin.re * bin.re + bin.im * bin.im). This gives the magnitudes of each of the points in the DFT bin.

Q1: These are "averaged" to form Pxx Pxx is a bin of size NFFT/2+1 which I assume store the averages, but, do these store the averages for each of the bins? Because I don't quite understand:

If I have a resulting 2D matrix containing values:

a = [[1, 2, 3, 4, 5, 6, 9, 8],
     [2, 5, 8, 9, 7, 8, 9, 7],
     [4, 5, 6, 4, 7, 8, 9, 8]

I then take the "average" of these.. I.e. (1+2....., n)/NFFT/2+1. This is only going to give one value, for each of the bins?

OK, so now onto Scaling. Can we therefore scale to the Frequency? If so, do we scale by Fs * window_size

EDIT:

I don't think my last question made that much sense..

So far I have a 2D vector containing 146x128..

block1 = [1, 2, 3, ...... n] (where `n` is 128) 

block2 = [0, 5, 6, ...... n] 

.....

block146 = [1, 4, 5, .... n]

Now, if I calculate the average magnitude of each of these blocks, i.e. sum(sqrt(re*re+im*im))/128 this will give me the averages for each block (a vector containing 146 samples).

For each of my samples in the blocks do I therefore divide each of these by the average magnitude for that particular block? I.e.

block1 = [0, 1, 2, ...... n]
average = (0+1+2)/128 = 0.58 (for example)

result = [0/0.58, 1/0.58, 2/0.58.... n]

Therefor my resulting vector will still be a 2D vector.

Does this make sense?

EDIT:

1 - What type of signal are you trying to process? I.e: What type of application is it? Is it a unidimensional signal, like Audio sampled at 44.1Khz? Or is it a multidimensional signal like an image or video?

I'm trying to process a Audio signal.

What I have done so far is:

Taken a 1D input vector, this vector is then split into blocks (of size 256 with an overlap/hop of 128 this gives me a total of 146 blocks with each containing 128 samples. These samples have been multiplied by the Hamming Window and then a DFT has been passed through.

These blocks contain the resulting values after Windowing and DFT has taken place.

2 - Where are you performing these calculations? Matlab, Octave, python, or some other?

I am performing these calculations using C++

3 - When you say 2D vector (maybe could be less confusing to call this a Matrix), What does the numbers inside mean? what are the columns and what are the rows?

The columns and rows contain the DFT results of the Hanning windowing of each block at 256/2 = 128

Essentially, what I'm attempting to do is create a spectrogram. Currently, I can plot the spectrogram with magnitudes of the values, i.e. (re*re+im*im) however, instead of plotting the magnitudes, I want to plot the overall power spectrum (PSD) so I can determine where the power is mainly within the signal.

If I do this using matplotlib I can get the following results:

So in theory, instead of showing the magnitudes which is what I'm currently doing, I want to show the PSD (like in the graph). I looked at the matplotlib source code and this is what was given:

result, windowVals = apply_window(result, window, axis=0,
                                  return_window=True)
result = np.fft.fft(result, n=pad_to, axis=0)[:numFreqs, :]

result = np.conjugate(result) * result

# for PSD
result /= (np.abs(windowVals)**2).sum()

result[1:-1] *= scaling_factor

# where scaling_factor is either 1 or 2

To me, the above code looks like they are taking the squared sum of the WindowValues (Which is Hamming) and dividing each value from the result (which is the resulting DFT) and finally multiplying by the scaling factor.

Does this look correct?

Answer 1

First of all, what do you mean when you say "Taken the first half of the resulting matrix (n/2+1)"?

The FFT will result in a vector, but maybe you mean taking only the first half of the FFT samples because the others are the negative frequencies and hence redundant?

Anyway, the process you describe is actually correct, and there is no need to do anything else. In your example with the resulting matrix a, you aparently took an FFT of size 8 (Or 16 if you then choped the second half of it). But what you do is simply take the average of all the first numbers of these vectors, and that would be an estimate of the Power spectral Density at f = 1*fs/NFFT. Taking the average of all the second number will be an estimate of PSD at f = 2*fs/NFFT, etc.

In the example the PDS would be

PSD = 1/3*[(1+2+4) , (2+5+5) ,(3+8+6) , ... ];

When you are analizing in real time as you carry out the computation, if the signal is Wide sense stationary, then it should converge or hover around some value, but if it is not, it will deviate over time to adjust to the new coming statistics of the signal.

If i understood you correctly, that should be it.

Edit: oops i forgot about the important scaling factor. The scaling factor should be introduced to preserve the parseval theorem between the PSD and its transform (The autocorrelation). If you windowed your signal first, i think you should apply the scaling factor 1/sum(W*fs/NFFT), where W is your window vector. Since you are merely dividing by the integral of the window function in the frequency domain This is what you add to the estimate of the frequency component when you do the FFT by windowing first. Because the FFT alone already preserves the energy.

Edit 2:

I am even more confused with the 2-D thing. Let me try an ask you a few questions:

1 - What type of signal are you trying to process? I.e: What type of application is it? Is it a unidimensional signal, like Audio sampled at 44.1Khz? Or is it a multidimensional signal like an image or video?

2 - Where are you performing these calculations? Matlab, Octave, python, or some other?

3 - When you say 2D vector (maybe could be less confusing to call this a Matrix), What does the numbers inside mean? what are the columns and what are the rows?

4 - Did you check the link to mathworks? i think you will find what you are trying to do, is well explained there.

Also, just to make sure you understand the backgroun well enough, you could refer to some classical books that have the topic very well explained:

• Digital Signal Processing, Principles, Algorithms And Applications -Proakis & Manolakis - 3rd Ed. Ch. 12

And the Father of all books in spectral Estimation:

• Spectral Analysis of Signals - Stoica & Moses

Edit 3:

Ahh, ok, now i get what you are trying to do. So, here's the thing. I gave a quick read to the sections on FFT and PSD in the books by Stoica and Proakis, and still can't quite understand why they propose a scaling factor of simply 1/N to the FFT, while the matlab function for PSD using periodogram method scales the FFT by 2/(Fs*N), being N the length of the FFt, in your case 256.

So i don't know a lot of C++ but i would say do the following:

1 - Take first frame X, which consists of 256 samples of audio signal.

2 - apply 256 FFT to it. you obtain Xfft.

3 - To each complex value Xfft(i), for i = 0...255, do xfft(i)*conj(xfft(i)) which is equal to abs(xfft(i))^2. this you already did,so up to here nothing new.

4 - Multiply the whole vector Xfft by 2/(Fs*N).

This will give you the Power spectral density for your frame X. Now take the next frame and do the same, And when you reach a number of frames that you would like to plot, you stop and plot your vectors.

now since you are taking frames every 128 samples, your estimates at each time step will be sort of crammed up together, so a trick that youwould use is the fact that a Speech/audio signal is approximately wide sense statinary (WSS) in an interval of 20-50 msec (See Speech Enhancement - Jacob Benesty). So what you could do, is assuming you have an Fs = 44100, that gives you a period of 50msec of around 2205 samples.

With the dimensions of the frames you took you could easily average 15 or 16 overlaped frames and get a Smoother estimate of the PSD. But the choice is yours, experiment with it and have fun!