Clear Filters
Clear Filters

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

8 views (last 30 days)
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 !
  1 Comment
Youssef  Khmou
Youssef Khmou on 10 Oct 2013
Edited: Youssef Khmou on 10 Oct 2013
whats is the estimated variance from you velocity signal? and how is the shape of its PSD?

Sign in to comment.

Answers (4)

Youssef  Khmou
Youssef Khmou on 10 Oct 2013
Your method works with these signals :
x=chirp(t,100,1,300);
y=randn(size(t));
waiting to answer the comment above..
  1 Comment
Venkata
Venkata on 10 Oct 2013
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.

Sign in to comment.


Matthew Crema
Matthew Crema on 10 Oct 2013
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?
  4 Comments
Matthew Crema
Matthew Crema on 10 Oct 2013
Edited: Matthew Crema on 10 Oct 2013
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
Venkata on 21 Oct 2013
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.

Sign in to comment.


Marek
Marek on 1 Dec 2014
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')

Venkata
Venkata on 4 Dec 2014
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.

Products

Community Treasure Hunt

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

Start Hunting!