how to improve the resolution of the PSD using Matlab

Question

I have and audio signal, which I read with Matlab, and use pwelch to get its PSD, here ist the code that I'm using

[x,Fs] = audioread('audioFile.wav');
x= x(:,1) %  mono
[xPSD,f] = pwelch(x,hamming(512),16,1024,Fs);
plot(f,xPSD);

since the FS=96000 and I'm only interrested in Frequencies bellow 5khz, I would like to calculate the PSD only for the area, and also being able to adjust the resolution of the PSD ! any idea hwo to do that !

Answer 1

When calculating PSDs with pwelch, there is always a trade-off between spectral resolution, number of averages and the amount of data you need. My preferred way to use it is:

[psd_x, freq] = pwelch(x, hann(nfft), nfft/2, nfft, fsample);

A few differences with your code:

1. I prefer to use hann windows, since I have bad experiences with hamming windows, they are not very good if your signal contains for example a large DC component. See this comparison, which shows that the roll-off of hann is much better, at the only cost of a slightly higher first sidelobe.

2. I use windows that overlap 50% (by using noverlap = nfft/2), in this way you 'get the most out of your data'. In your case, there is only 16/512 = 3% overlap between windows, and since the window function is close to zero at its edges, the data points at the edges do not contribute as much as points in the middle of a window. With half-overlapping windows this effect is much smaller. Making the overlap bigger than 50% is useless, you will get more averages, but since you will use the same points many times, this does not add any extra information. Just stick to 50%.

3. I usually make the length of the fft (fourth argument to pwelch) the same as the window length. The only case you want to have this different is when you use zero-padding, which has limited use.

There are a few simple formulas, which you should memorize when working with pwelch and similar functions:

• The spectral resolution is given by the window length only: df = 1 / t_window

• The length of a single window is t_window = nfft / f_sample.

• With half-overlapping windows, the total amount of needed data is t_total = t_window * (n_average + 1) / 2

• For one-sided spectra, the number of spectral bins of the PSD is nfft / 2.

• Nyquist: f_max = f_sample / 2

To get a reasonably smooth spectrum, I would typically use on the order of 20 averages. Combining the equations above and filling in the required spectral resolution then gives you the total length of data you need. Or the other way around, if you only have a limited amount of data available, you could calculate the frequency resolution you can obtain.