how to use a MATLAB function called MC estimator to determine frequency?

2 次查看（过去 30 天）
modified covariance
modified covariance 2012-10-25
the following is the matlab function for the MC estimator
function w = mc(x);
% w = mc(x) is used to estimate frequency
% the MC method from the vector x
% x is supposed to be a noisy single-tone sequence
% w is the estimated frequency in radian
N=max(size(x));
t1=0;
t2=0;
for n=3:N
t1=t1+x(n-1)*(x(n)+x(n-2));
t2=t2+2*(x(n-1))^2;
end
r = t1/t2;
if (r>1)
r=1;
end
if (r<-1)
r=-1;
end
w=acos(r);
attempted to determine the frequency estimate of the data, which can be found at the following link,using the MC function.
http://solarscience.msfc.nasa.gov/greenwch/spot_num.txt
(please note that the first row and rows after July 2012 are removed)
this is my attempt where the data of thrid column are used only (observing the data, sampling interval is 1 month, that is Fs=1)
load spot_num.txt
ssn = spot_num(:,3);
ssn = ssn-mean(ssn);
x=ssn;
N=max(size(x));
t1=0;
t2=0;
for n=3:N
t1=t1+x(n-1)*(x(n)+x(n-2));
t2=t2+2*(x(n-1))^2;
end
r = t1/t2;
if (r>1)
r=1;
end
if (r<-1)
r=-1;
end
w=acos(r);
w should be close to 0.0076 according to model answer, but somehow i just couldnt get the right answer? could anyone help, thx!

请先登录，再回答此问题。

回答（1 个）

Wayne King
Wayne King 2012-11-2
编辑：Wayne King 2012-11-2
It may be that this method is just not very robust for noise-corrupted data.
For example:
t = 0:159;
x = cos(pi/4*t);
w = mc(x);
I get 0.7854, which is pi/4 and the method does a perfect job, but if I add some noise
xnew = x+0.5*randn(size(t));
plot(xnew)
w = mc(xnew);
then the estimate is way off.
Why not use the periodogram, which will get the frequency correct for your sunspot data?
[Pxx,F] = periodogram(ssn,rectwin(length(ssn)),length(ssn),1);
[~,idx] = max(Pxx);
F(idx)
  1 个评论
modified covariance
modified covariance 2012-11-2
编辑：modified covariance 2012-11-2
thank you Wayne King, the periodogram u suggested certainly gets the answer which is 0.0076.
unfortunately, i must have to use the following code
load spot_num.txt
ssn = spot_num(:,3);
ssn = ssn-mean(ssn);
together with MC estimator function to get the frequency "As Required by Assignment Guidelines"...
but anyway, thx a lot!

请先登录，再进行评论。

类别

在 Help Center 和 File Exchange 中查找有关 Logical 的更多信息

标签

Community Treasure Hunt

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

Start Hunting!

Translated by