Integration of pwelch output and comparing it with the variance of a signal

Hi,

On a sinusoidal signal with noise (an example given in matlab help for fft), I am using a rectangular window, with zero overlap and using the entire data length as one segment in pwelch. I obtained the output of pwelch, integrated the output and checked if it is equal to the variance of the signal. The example is shown below:

Fs = 1000; % Sampling frequency

T = 1/Fs; % Sample time

L = 1000; % Length of signal

t = (0:L-1)*T;% Time vector

% Sum of a 50 Hz sinusoid and a 120 Hz sinusoid

x = 0.7*sin(2*pi*50*t) + sin(2*pi*120*t);

y = x + 2*randn(size(t));

yvar=var(y);

na=1;

wr=rectwin(floor(L/na));

[pyr,fr]=pwelch(y,wr,0,NFFT,Fs);

diffpyvar=abs(trapz(pyr*(fr(2)-fr(1)))-yvar)*100/yvar

I found that diffpyvar was only 0.15%.

However if I repeat the same procedure on a real signal (velocity output from my computations) whose length is chosen as NFFT = 2^18 and keeping the rest of the parameters same, the relative difference between the pwelch output and the variance of the signal is over 1000%.

I cannot understand why the differenc is huge and I am grateful to those who can answer this riddle !

[Youssef Khmou]

whats is the estimated variance from you velocity signal? and how is the shape of its PSD?

Answer 1

[Youssef Khmou]

Your method works with these signals :

x=chirp(t,100,1,300);
y=randn(size(t));

waiting to answer the comment above..

[Venkata]

Hi Youssef, Thanks for the reply. I don't think that the pwelch method works for chirp signals alone. As per the theory, the integration of power spectral density of any signal should be equal to the variance and matlab pwelch gives the power spectral density.

Answer 2

[Matthew Crema]

It seems to me that the equality you are assuming is var(y) == int(pyr*deltaf)

And I agree, this should be correct. However I find that by modifying your test data that the error gets bigger if I make L large, and change some of the parameters to the pwelch statement. For example, try L to 100000.

[pyr,fr]=pwelch(y,wr,0,L/4,Fs);

Instead, try replacing your pwelch statement with

[pyr,fr]=pwelch(y,[],[],[],Fs);

I find that the above code still works with your example data, and the error remains small with my test data. But does it work with your actual data?

[Matthew Crema]

Another thing you might try as a test is is changing two lines:

[pyr, fr] = periodogram(y, rectwin(L), L, Fs, 'onesided', 'ms');

diffpyvar=abs(trapz(pyr)-yvar)*100/yvar

For some reason I feel that although this may give you a poorer estimate of the power spectrum, it may give you a better estimate of the total power.

[Venkata]

Thank you for the reply Matthew. I followed your suggestion of using periodogram on a data of length 2^18. Then the percentage relative difference between the integration of the power spectral density and variance has come down to 99% from 826% that was found with pwelch. Still a large difference ! May I know what is 'ms' in the periodogram function ?

[Matthew Crema]

I should also think 99% is too much.

'ms' tells MATLAB to return the specturm in units of Mean-squared power and not in units of spectral density (mean-squared power per Hz). There may be more to it than this, but basically I think it does not divide by the sample rate when computing pyr. That is why you also have to change the line of code that computes the total power. Note that you should compute the "sum" of the components instead of the area under the density (do not multiply by delta f).

What about the first test, where you let MATLAB choose the window length, noverlap and nfft by passing [] to pwelch for those inputs?

[pyr,fr]=pwelch(y,[],[],[],Fs);

Generally speaking, I don't think it's a good idea to have a window length equal to the length of the signal (because you effectively aren't windowing at all). I'd recommend starting with the above statement and then tailoring the inputs based for your particular signal.

[Venkata]

Thanks for the clarification on 'ms' used in periodogram.

I have intentionally avoided windowing in pwelch as I would like to see if the integration of output from pwelch gives the variance of a signal used. I used the default arguments in pwelch, after which the integration of its output did not yield the variance! That's the reason I was using a single rectangular window to make further checks.

Answer 3

[Marek]

Had the same problem,

found a better solution for me, so I thought sharing might help other people.

[pyr,fr]=pwelch(y,[],[],[],Fs);

was a good hint, but the variance was still far off.

Use this function instead of integrating with trapz which gives me great results: variance = bandpower(pyr,fr,'psd')

See: http://se.mathworks.com/help/signal/ref/bandpower.html#btrdcz_-3

Answer 4

[Venkata]

Hi Marek, Thanks for your reply. I found the following also answers my question. [pyr,fr]=pwelch(detrend(y,'constant'),[],[],[],Fs); Then the difference between variance of y and trapz(pyr,dfr), where dfr is the frequency bin, is found to be small.