MATLAB Answers

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

5 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!

Translated by