Integration of pwelch output and comparing it with the variance of a signal
8 views (last 30 days)
Show older comments
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
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?
Answers (4)
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..
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
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.
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')
0 Comments
See Also
Categories
Find more on Spectral Estimation in Help Center and File Exchange
Products
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!