use modified covariance (MC) method to estimate frquency?

let x(t)=10cos(200*pi*t+1.2) which is a continuous-time sinusoid. The x(t) is sampled every Ts=0.001 sec. to obtain a sequence x[n] where x[n]=x(nTs) for n between 0 and 100 (including 0 and 100)

how can I use modified covariance (MC) method, given the following code, to determine the frquency estimate for x[n].

 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);

btw, What is purpose of setting r to 1 and -1 when its value is larger than 1 and smaller than -1,respectively?

thank u in advance.

Answer 1

[Wayne King]

I'm not familiar with this algorithm for frequency estimation. When I use the term "modified covariance method", it's in the context of autoregressive spectral estimation, but having said that it seems to work here.

 t = 0:0.001:0.1-0.001;
 x = 10*cos(200*pi*t+1.2); %although this signal is not corrupted by noise
 w = mc(x);

returns 0.6283 radians/sample, that is equivalent to 628.3 radians/second which is 100 Hz (the frequency of the input sinusoid in cycles/second).

To answer your final question, the real-valued cosine function takes values between [-1,1]. If you take the inverse cosine of values in the interval [-1,1], you will get a real-valued output. However, if you take the inverse cosine of a value larger than 1, or smaller than -1, you will get a complex-valued output. The complex-valued cosine is not bounded by 1.

The code is ensuring that you do not input any values to acos() outside the interval [-1,1], so that the output is real-valued.

To use the function, mc.m, save the file in a folder on the MATLAB path and you can run it from the command line. For example, if you have the file in a folder called c:\mfiles, use

>>addpath 'c:\mfiles'

to add that folder. or you can use

>>pathtool

[modified covariance]

thx a lot for ur answer, but i am just wondering how did u manage to get the answer? you simply typed the following?

 t = 0:0.001:0.1-0.001;
 x = 10*cos(200*pi*t+1.2); %although this signal is not corrupted by noise
 w = mc(x);

----------------------------------------------------------------------- another question is why the following code is provided since w=mc(x) seems already able to return the answer?

 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);

[Wayne King]

because having all those commands collected in a function is a convenient way to call it. You don't have to remember to copy and paste all those lines of code, or type them one by one. You just provide your input to the function and you're done! It's the same reason you bundle any algorithm into a function.

Yes, you can create your signal in the workspace, call it x, and copy and paste all that code and you get the answer, but using the function is more convenient. You don't have to call your input x for example if you use the function.

 t = 0:0.001:0.1-0.001;
 y = 10*cos(200*pi*t+1.2); 
 w = mc(y);

or

 t = 0:0.001:0.1-0.001;
 sig = 10*cos(200*pi*t+1.2); 
 w = mc(sig);

works.

[modified covariance]

i was unable to generate the output,somehow. but anyway thank u for ur kind replys in detail.

could u plz do me a favor to generate the frequency if now x(t)=10cos(200*pi*t+1.2) is changed to x(t)=10cos(267*pi*t+1.2)?

[Wayne King]

it returns 0.8388 as expected, you can just copy and paste the following to see this:

 t = 0:0.001:0.1-0.001;
 x = 10*cos(267*pi*t+1.2);
 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)

[modified covariance]

thank god!!! it finally works, such an idiot i am.

once again, thank u!