Integration of pwelch output and comparing it with the variance of a signal
15 次查看(过去 30 天)
显示 更早的评论
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 个评论
Youssef Khmou
2013-10-10
编辑:Youssef Khmou
2013-10-10
whats is the estimated variance from you velocity signal? and how is the shape of its PSD?
回答(4 个)
Youssef Khmou
2013-10-10
Your method works with these signals :
x=chirp(t,100,1,300);
y=randn(size(t));
waiting to answer the comment above..
Matthew Crema
2013-10-10
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 个评论
Matthew Crema
2013-10-10
编辑:Matthew Crema
2013-10-10
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
2014-12-1
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 个评论
另请参阅
类别
在 Help Center 和 File Exchange 中查找有关 Parametric Spectral Estimation 的更多信息
产品
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!