pwelch power correction factor

[Henik Vargas]

Hi all,

I am doing some spectrum analysis and I am thinking of using the pwelch function with a hann window. Now, I understand that application of a window reduces the power of the signal, so usually a correction or scaling factor has to be applied.

My question is, does anyone know if the pwelch function applies this correction factor? I implemented my own welch function to compare and the result is identical, except that MATLAB's pwelch function seems to ignore the correction factor (or I am doing something wrong)

Thank you,

[Rune Allnor]

What makes you think you need a correction factor?
Instead of replicating what matlab does, read up on
the method you play with.

Rune

[Henik Vargas]

> What makes you think you need a correction factor?
> Instead of replicating what matlab does, read up on
> the method you play with.
>
> Rune

Rune,

The motivation for the correction factor comes from the fact that windowing reduces the power of the signal, therefore, the power of the psd is much lower than that of the unwindowed signal. Another way of looking at it is that the correction factor makes the power spectrum estimate asymptotically unbiased (See Hayes, Statistical Digital Signal Processing and Modeling for a reference).

Also, I read the documentation for the function, and I think my question is a fair one. I could just blindly use the function as is, but then I lose a learning opportunity.

[Rune Allnor]

My point is that in the applications where spectrum estimation
methods like Welch's method are used, the main focus is to get
an *impression* of the spectrum shape. The exact numbers are
irrelevant.

As I said, read up on the methods.

Rune

[Henik Vargas]

>
> My point is that in the applications where spectrum estimation
> methods like Welch's method are used, the main focus is to get
> an *impression* of the spectrum shape. The exact numbers are
> irrelevant.
>
> As I said, read up on the methods.
>
> Rune

Thank you Rune. I think your point about getting the "impression" of the spectrum is a valid one. However, you keep implying that I am not reading the documentation. But I am and it says:

"Pxx is the distribution of power per unit frequency. For real signals, pwelch returns the one-sided PSD by default; ... Note that a one-sided PSD contains the total power of the input signal."

So once again, I am questioning whether or not the PSD does in fact contain the total power of the input signal. Maybe I am missing something fundamental, but again, I think it is a fair question.

[Henik Vargas]

I found an answer to my question. I am posting it for reference in case anyone stumbles upon this posting:

http://www.mathworks.com/support/solutions/en/data/1-17M0Q/index.html?product=SG&solution=1-17M0Q

[Jeremy]

I realize this is really old, but I thought I would answer this for any future readers...

pwelch appears to apply a broadband/random window correction factor. This means that the amplitude of a discrete sinusoidal will not be correct, but it can be corrected for by the ratio of the discrete signal correction factor over the random broadband/random factor. This is not discussed in the documentation anywhere, but that is what comes out of pwelch.

If you input a window function besides the default (hamming) it still applies the correction factor so it must calculate the factor from the weighted function the goes in as the "window" input, without knowing the name of the function being used. You could make up a function and I think it would still apply the correction factor.

For example, the broadband correction factor for a hanning window is 8/3. for a discrete signal the correction factor is 4.0 (2.0 on the amplitude spectrum).
So multiplying the PSD by 4/(8/3) = 1.5 will give you the correct amplitude of a sinusoidal.

This can be demonstrated by the script below:

fs=1000;
t=[1/fs:1/fs:100]';
disp(['input signal (rms):' num2str(20/sqrt(2))]);
x=randn(length(t),1)*10+20*cos(2.*pi.*56*t);
nfft=1000;
p=pwelch(x,hanning(nfft),nfft/2,nfft,fs);
disp(['calculated amplitude (rms):' num2str(sqrt(p(57)*1.5))]);
% or can just sum the 56 Hz bin with the adjacent bins where the window spilled into:
disp(['calculated amplitude (rms):' num2str(sqrt(sum(p(56:58))))]);

[olga]

Very helpful post! Thank you. I did some testing myself as well and came to the same conclusion. I experimented with a rectangular window in addition to Hanning; the idea was to see if scaling factor is needed when no window is used. The answer is no. So we can drop 1.5 (or any other factor), when rectwin is used in place of a window in pwelch.

Also, to get one to one correspondence between signal amplitude and Pxx amplitude, simply use sort(2*Pxx*scaling factor) correction factor (see code below). Scaling factor is 1 when rectwin is used or 1.5 when Hanning is used.

%
% PWELCH SCALING TEST
%
% FINDINGS: x1 and x2 are of different durations. x2 is longer than x1.
% Pxx2 will be proportially larger than Pxx1 (if x2 is twice as long as x1,
% Pxx2 will be twice as large as Pxx1). However, if window length in
% pwelch is set to the same value for both x1 and x2, Pxx1=Pxx2. However,
% in both cases, Pxx1 and Pxx2 are not equal to signal power (amplitude of
% x1 or x2 squared). Amplitude normalization is required. When hanning
% window is used, the correction is 1.5*Pxx. (1.5=4/(8/3). 8/3 for window
% and 4 for discrete signal). pwelch appears to apply a broadband/random
% window correction factor. This means that the amplitude of a discrete
% sinusoidal will not be correct, but it can be corrected for by the ratio
% of the discrete signal correction factor over the random broadband/random
% factor.
%
% PWELCH Scaling Factor with Hanning Window: sqrt(2*Pxx*1.5). It
% correspondsto signal amplitude. If no window (rectwin), then no 1.5
% factor is needed. Just use sqrt(2*Pxx).
%
Fs=128;
t1=0:1/Fs:2-1/Fs;A1=1;
t2=0:1/Fs:4-1/Fs;A2=1;
x1=A1*cos(2*pi*20*t1);RMS1=A1/sqrt(2); %RMS voltage
fprintf(1,'Vrms1=%f\n',RMS1);
x2=A2*cos(2*pi*20*t2);RMS2=A2/sqrt(2); %RMS voltage
fprintf(1,'Vrms2=%f\n',RMS2);
%
N1=length(x1);
N2=length(x2);
%
[Pxx1,F1] = pwelch(x1,hanning(N1/2),0,N1/2,Fs);Prms1=sqrt(max(Pxx1)*1.5);
[Pxx2,F2] = pwelch(x2,hanning(N1/2),0,N1/2,Fs);Prms2=sqrt(max(Pxx2)*1.5);
%
fprintf(1,'Prms1=%f\n',Prms1);
fprintf(1,'Prms2=%f\n',Prms2);
%
figure(1);
% subplot(2,1,1);plot(F1,sqrt(Pxx1*1.5));grid on;%ylim([0 1]); %RMS power
% subplot(2,1,2);plot(F2,sqrt(Pxx2*1.5));grid on;%ylim([0 1]); %RMS power
%
subplot(2,1,1);plot(F1,sqrt(2*Pxx1*1.5));grid on;%ylim([0 1]); %Voltage Amplitude
subplot(2,1,2);plot(F2,sqrt(2*Pxx2*1.5));grid on;%ylim([0 1]); %Voltage Amplitude