MATLAB Answers

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

2 views (last 30 days)
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

Accepted Answer

David Goodmanson
David Goodmanson on 30 Jul 2021
Edited: David Goodmanson on 2 Aug 2021
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
  3 Comments

Sign in to comment.

More Answers (0)

Tags

Products

Community Treasure Hunt

Find the treasures in MATLAB Central and discover how the community can help you!

Start Hunting!