Does pwelch function adopt already the correction factor due to the window?

Does pwelch function adopt already the correction factor due to the window?

When we use a window like hamming a correction factor for the amplitude or energy is needed. In my case the energy correction factor is needed as also explained here. Is it possible that pwelch does already a correction, or should I use a correction due to my window type. Does the overlap affect the energy?

I've seen that I am not the only one to ask this question like here, but nobody answered.

Thanks

Answer 1

[David Goodmanson]

hello sy,

At this point I believe that the documentation is misleading in that it appears to imply Power/unit_frequency but the result is in Power/radian, smaller by a factor of 1/2pi (but doubled for a one-sided spectrum). At least, that seems the most likely possibility. pwelch does not seem to be accessable code, so here is an example for which I constructed a function called pwelch1 to look at the scaling.

M = 10000;
y = 2*rand(1,M)-1;
S = rms(y)^2       % average power
figure(1); grid on
plot(y)            % not very interesting

[Z f] = pwelch1(y);
fig(2); grid on
semilogy(f,Z)
S2 = trapz(f,Z)  % agrees

Here the normalized f runs from -1/2 to 1/2. Integrating dS/df over all frequencies should give the average power S, which it does. Both are right around 1/3, the theoretical value. The function does contain the Hamming function correction factor U which is approximately 0.4, depending on the number of points.

pwelch1 always averages 8 windows with 50% overlap and there is no attempt to use ffts with length of 2^n since I think that effort is usually a waste of time.

In the Matlab code, w runs from 0 to pi, which certainly looks like normaliazed circular frequency. The code outputs positive frequencies only and multiplies the spectrum by A = 2 to make up for that. For Power/radian instead of Power/unit_frequecy, pwelch is smaller than pwlech1 by a factor of A/(2*pi) = 1/pi as you can see on the plot.

[z w] = pwelch(y);
figure(3); grid on
semilogy(w,z)
S3 = trapz(w,z)   % agrees

It would be interesting to hear what Mathworks might say on this toplic.

function [Z f] = pwelch1(y)

% similar to pwelch except that
% [1] the outputs are rows
% [2] The normailzed frequency array runs from -1/2 to 1/2
%     and the Z array is not doubled.  Z is also not divided by a factor of 2pi.
%
% function [Z f] = pwelch1(y)

N = length(y);
nwin = round(N/4.5);
h = hamming(nwin)';          % make it a row vector
U = (1/nwin)*trapz(h.^2);
indstart = zeros(1,8);
for k = 1:4                  % set up window starting points
  indstart(2*k-1)  = 1+(k-1)*nwin;
  indstart(10-2*k) = N+1-k*nwin;
end
Z = zeros(1,nwin);
for k = 1:8
  ind = indstart(k)+(0:(nwin-1));
  Z = Z + abs(fft(y(ind).*h)).^2;
end
Z = fftshift(Z)/(nwin*8*U);
f = (ceil(-nwin/2):ceil(nwin/2)-1)/nwin;  % works for odd or even
end

[SYML2nd]

I did the same, with an other code and I came to the same conclusion. Btw, I would like to have reference for this. I didn't find anything on the doc file about this correction.

The only things that can suggest a scaling in the guide is the following:

Power spectrum scaling, specified as one of 'psd' or 'power'. Omitting the spectrumtype, or specifying 'psd', returns the power spectral density. Specifying 'power' scales each estimate of the PSD by the equivalent noise bandwidth of the window. Use the 'power' option to obtain an estimate of the power at each frequency.

Unfortunately this is inconsistent with our result. Therefore omitting the spectrum type no scaling seems to be applied, but our results are consistent with the real energy

[SYML2nd]

Hi David. Have you applied the energy correction factor due to the window?

I am start thinking that there is something wrong in the documentation.

I hope that someone of the Matlab staff can answer.

[David Goodmanson]

Hello sy, see updated answer