Power Spectrum estimate using pwelch function

I am looking into using the pwelch function to estimate the power spectrum. As I understood, the power spectrum is simply the PSD multiplied by the frequency bandwidth (sampling rate divided by the block size). So I decided to check and compare the pwelch estimate with the PSD*BW approach. I found that for a rectangular window (no window), the methods match exactly. But if I use a window function, then the methods do not match. The ratio that they are off isn’t the window correction factor, but maybe it is related: for example with a hann window the estimated power spectra differ by a factor of 1.5x and for a hamming window it is a factor of 1.36x. Has anyone encountered anything like this or can it be explained? Attached is an example overlaying the two methods for the rectangular window and hamming window.

%% Inputs
fs=10000;           % sampling rate (sps)
t1=1;               % specified start time (s)
t2=9;               % specified end time (s)
BS=1024;            % block size (pts)
nfft=BS;            % DFT block size (pts)
OL=0.5;             % percent overlap (%/100)

%% Example Data
TimeData=0:1/fs:10;                 % time vector (s)
Data=rand(1,length(TimeData));      % data vector (EU)
tStart=TimeData(1);                 % data start time

%% PSD/PS Calculation
BW=fs/BS;               % bandwidth  (Hz)
df=fs/nfft;             % frequency resolution (Hz)
ptsOL=floor(OL*BS);     % number of points for overlap (pts)
PSDStartIndex=round((t1-tStart)*fs)+1;   % find index closest to specified PSD start time
PSDEndIndex=round((t2-tStart)*fs)+1;     % find index closest to specified PSD end time
[pxx_rect1,f_rect1]=pwelch(Data(PSDStartIndex:PSDEndIndex),ones(1,BS),ptsOL,nfft,fs);
[pxx_rect2,f_rect2]=pwelch(Data(PSDStartIndex:PSDEndIndex),ones(1,BS),ptsOL,nfft,fs,'power');
[pxx_hamming1,f_hamming1]=pwelch(Data(PSDStartIndex:PSDEndIndex),hamming(BS),ptsOL,nfft,fs);
[pxx_hamming2,f_hamming2]=pwelch(Data(PSDStartIndex:PSDEndIndex),hamming(BS),ptsOL,nfft,fs,'power');



%% Plotting
% Power Spectrum Plot - Rectangular Window
figure(1)
semilogy(f_rect1,pxx_rect1*BW)
hold on
    semilogy(f_rect2,pxx_rect2)
hold off
grid on
set(gca,'xlim',[0 fs/2])
set(gca,'ylim',[10^-4 10^-3])
xlabel(gca, 'Frequency (Hz)')
ylabel(gca, 'Amplitude (EU-rms)^2')
legend({'PSD*BW','Power'}, 'Location','NorthEast')
title(gca,{['\fontsize{14}' 'Power Spectrum - Rectangular Window' '\fontsize{10}']},'fontweight','bold')

% Power Spectrum Plot - Hamming Window    
figure(2)
semilogy(f_hamming1,pxx_hamming1*BW)
hold on
    semilogy(f_hamming2,pxx_hamming2)
hold off
grid on
set(gca,'xlim',[0 fs/2])
set(gca,'ylim',[10^-4 10^-3])
xlabel(gca, 'Frequency (Hz)')
ylabel(gca, 'Amplitude (EU-rms)^2')
legend({'PSD*BW','Power'}, 'Location','NorthEast')
title(gca,{['\fontsize{14}' 'Power Spectrum - Hanning Window' '\fontsize{10}']},'fontweight','bold')

[Sudarsanan A K]

Hello Matthew,

I understand that you are using the "pwelch" function with different window functions (rectangular and Hamming) and comparing the results with the PSD multiplied by the frequency bandwidth. You found that for a rectangular window (no window), the results match exactly. But when you use a window function, then the results do not match. You also noted that the difference in the results is not just by the correlation factor, and it also varies with the window function.

The difference you observed between the pwelch estimate and the PSD*BW approach when using a window function is expected. The reason for this discrepancy is the effect of the window function on the spectral leakage and the effective bandwidth of the signal.

The pwelch function already considers the effect of the window function on the PSD estimate. It applies a correction factor to account for the reduction in the effective bandwidth caused by windowing. Therefore, when you compare the pwelch estimate with the PSD*BW approach, you should expect some differences due to the window correction factor.

The specific factors you mentioned, such as 1.5x for a Hann window and 1.36x for a Hamming window, are likely due to the characteristics of these window functions and their impact on the effective bandwidth.

In summary, the differences you observed between the "pwelch" estimate and the PSD*BW approach when using a window function are expected because of the window function on spectral leakage and the effective bandwidth. The pwelch function already accounts for this effect through a correction factor, which explains the discrepancy between the two methods.

To gain a deeper understanding of the "pwelch" function and its nuances, you can additionally refer to the MathWorks documentation, which provides comprehensive information and examples on how to use this function effectively.

https://mathworks.com/help/signal/ref/pwelch.html#btuf68p-2

Hope this helps.

[Matthew Casiano]

Thank you for the comprehensive review. It’s quite clear and descriptive.

However, I’m still uncertain about which method yields the most representative power spectrum. I might argue that before calculating power, the ‘power’ option in the pwelch function should first remove the window correction factor.

In the example I provided, I believe the calculations using the rectangular window are accurate. They yield identical spectra whether we use the PSD*BW method or the pwelch 'power' option method, as seen in the first plot. There’s no window correction factor in this case. Moreover, if we overlay the PSD*BW method using the Hamming window with these, it also matches. The only outlier in this example is the pwelch 'power' option method using the Hamming window.

While this is by no means a rigourous review, I suspect that the pwelch 'power' option method should first remove the window correction factor. The PSD spectrum should have the correct overall RMS before power calculation, which is only achievable without window corrections. {deleted last sentence because it was confusing}

Thanks!

Answer 1

[Sudarsanan A K]

Hello Matthew,

I understand that you are using the "pwelch" function with different window functions (rectangular and Hamming) and comparing the results with the PSD multiplied by the frequency bandwidth. You found that for a rectangular window (no window), the results match exactly. But when you use a window function, then the results do not match. You also noted that the difference in the results is not just by the correlation factor, and it also varies with the window function.

The difference you observed between the "pwelch" estimate and the PSD*BW approach when using a window function is expected. The reason for this discrepancy is the effect of the window function on the spectral leakage and the effective bandwidth of the signal.

The "pwelch" function already considers the effect of the window function on the PSD estimate. It applies a correction factor to account for the reduction in the effective bandwidth caused by windowing. Therefore, when you compare the "pwelch" estimate with the PSD*BW approach, you should expect some differences due to the window correction factor.

The specific factors you mentioned, such as 1.5x for a Hanning window and 1.36x for a Hamming window, are likely due to the characteristics of these window functions and their impact on the effective bandwidth.

In summary, the differences you observed between the "pwelch" estimate and the PSD*BW approach when using a window function are expected because of the window function on spectral leakage and the effective bandwidth. The "pwelch" function already accounts for this effect through a correction factor, which explains the discrepancy between the two methods.

To gain a deeper understanding of the "pwelch" function and its nuances, you can additionally refer to the MathWorks documentation, which provides comprehensive information and examples on how to use this function effectively.

https://mathworks.com/help/signal/ref/pwelch.html#btuf68p-2

To the best of my knowledge, if you need the correct overall RMS in the power spectrum, the 'PSD*BW' method or using the rectangular window is appropriate as they do not include the window correction factor.

However, if you want accurate power estimation within windowed segments and need to compensate for spectral leakage, the 'power' option in "pwelch" with window correction is useful. It provides accurate power estimation within the windowed segments, but the overall RMS may differ.

Hope this clarifies your query.

Answer 2

[David Goodmanson]

Hi Matthew,

I don't think there is an issue here, exactly. pxx_hamming1 is the default psd for pwelch. If you take pxx_hamming1*BW, the fig. 2 blue plot, and compare it with the rectangular window curves in plot 1, they agree pretty well. So pwelch is doing the window correction factor correctly.

(Incidentally

Data=rand(1,length(TimeData)) -1/2;      % data vector (EU)

is probably a better and more realistic choice since that data has mean close to 0 and for the psd makes the large peak at zero frequency go away).

The only question is what is going on with the 'power' option in pxx_hamming2, the fig. 2 red plot. It differs by a constant factor of

Q = mean(pxx_hamming2./(BW*pxx_hamming1))
  Q = 1.3638

at every frequency, not just the mean.

From help pwelch: " 'power' - scales each estimate of the PSD by the equivalent noise bandwidth of the window (in hertz). Use this option to obtain an estimate of the power at each frequency. "

Well, I guess. pwelch is doing something funky with the 'power' option scaling and after trying several examples I have not been able to figure it out. Maybe someone else will weigh in.

These days you can't use type pwelch (or edit pwelch) to see what's happening. It's a black box. In this case I am going with the philosophy, if you don't understand it, don't use it. So I don't use that option and for several other reasons wrote my own version of pwelch to have some control over what is going on.

[Matthew Casiano]

Thank you, you're comment provides helpful information.