5 views (last 30 days)

Show older comments

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?

Thanks

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

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

Start Hunting!