As previously discussed, finding the frequency of a complex sinusoid embedded in noise is a classical problem in Signal Processing. The problem is compounded by the fact that number of samples available is usually quite small. So far, we have discussed Zero Crossing, FFT, MUSIC and ESPRIT methods of frequency estimation. Zero Crossing method is simplest of the above four but it can detect only one sinusoid at a time. Advantage of Zero Crossing method is that it is computationally not that complex. It does not require complex matrix manipulations as some of the other methods do.
Now we present another single frequency estimator that is quite easy to implement. It is known as Kay’s Estimator [1] after its founder Steven M. Kay. Although quite simple to implement it achieves Cramer Rao Lower Bound (CRLB) for moderate to high SNRs (SNR>10dB). The basic concept of Kay’s estimator is quite simple. It finds the phase advanced from one sample to the next and this gives us the frequency. Averaging is performed over the time window to get more accurate results. Averaging can be performed using uniform weighting or a weighting function defined by Kay in [1]. In the latter case CRLB is achieved with equality, so there is no need to look for another estimator.
Mathematically CRLB is given as:
Mean Squared Error (MSE) in Angular Frequency=6/(SNR*N*(N^2-1))
Where N is the number of samples (set to 100 in this case) and SNR is the Signal to Noise Ratio (set to 10dB). MATLAB/Octave code for both the methods (with different weighting functions) is given below. For the example given CRLB is calculated as:
MSE=6/(10*100*(100^2-1))=-62.22dB
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% FREQUENCY ESTIMATION USING
% KAY's ESTIMATOR
% (SINGLE FREQUENCY ONLY)
% www.raymaps.com
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
clear all
close all
N=100;
fc=1.0e9;
fs=20*fc;
w=2*pi*fc/fs;
n=0:N-1;
s=sqrt(1.0)*exp(i*w*n);
wn=sqrt(0.1/2)*(randn(1,N)+i*randn(1,N));
x=s+wn;
% KAY ESTIMATOR-1 (uniform weighting)
w_est=0;
M=N;
for m=1:M-1
w_est=w_est+(1/(M-1))*arg((x(m))'*x(m+1));
end
error_1=w-w_est;
% KAY ESTIMATOR-2 (using weighting function)
w_est=0;
M=N;
for m=1:M-1
p(m)=6*(m)*(M-m)/(M*(M^2-1));
w_est=w_est+p(m)*arg((x(m))'*x(m+1));
end
error_2=w-w_est;
Note:
- If you want the complete code for comparing the MSE of the two methods with CRLB please leave a comment below and I will get back with you.
- Please note that as the number of samples N is increased the CRLB and Mean Squared Error (MSE) both go down. But after a time the law of diminishing returns sets in and MSE does not go down as rapidly as the CRLB.
- The MSE also depends upon number of samples per cycle and having a very small number does not help in reducing the MSE. A few complete cycles work the best (in our code we use five complete cycles to do the estimation).
Reference:
[1] S. M. Kay, “A Fast and Accurate Single Frequency Estimator,” IEEE Transactions on Acoustics, Speech, and Signal Processing, vol. 37, no. 12, Dec 1989.
A concern always with methods like this is phase wrapping. Steven Kay claims in his paper to have solved this problem (which he correctly notes was an issue with Tretter’s equivalent method – which just involved computing a linear regression of the phases).
However, I don’t think the problem has really been solved. Any (weighted) sum of phases relies entirely on those phases not straddling the “wrap phase” (typically chosen to be either -π/π or 0/2π). That means our choice of where the phase wraps can have a major impact on estimation performance.
For example, wrapping at -π/π will tend to give good results for |ω| close to zero (but terrible results for ω close to ±π)… and vice versa if we choose to wrap our phases at 0/2π. If we’re lucky, we’ll perhaps only be interested in some narrow sub-band of frequencies, in which case one choice of “wrap phase” may work just fine. However, more generally, it would perhaps be desirable to pre-rotate our N phases onto some optimal center phase, to minimize the risk of wrapping (for those N specific phases). I would be interested to hear your thoughts on this.
It’s worth noting that Kay’s simplified method described in the Appendix of his paper (where the order of summation and angle calculation are interchanged) does completely solve the phase-wrapping problem. However, it only achieves really good performance (CRLB) at much higher SNRs.
Thank you very much dear John for your response. It give the following error:
Undefined function ‘arg’ for input arguments of type ‘double’.
Error in KAYSingleFrequencyEstimator (line 24)
w_est=w_est+(1/(M-1))*arg((x(m))’*x(m+1));
If arg(Z) does not work try angle(Z). I would also recommend that if you do not have licensed MATLAB than use Octave, its free.
Thank you very much dear John for your guidance. Your 1st point in above note is as:
1- If you want the complete code for comparing the MSE of the two methods with CRLB please leave a comment below and I will get back with you.
So can you give the complete code here. Thanks.
CRLB(SNR_index)=6*(1/SNR_linear)/(N*((N^2)-1))
I will let you figure out the rest. Its not that difficult!
Thank you very much dear John for your response. I didn’t get you. What should I do with the command:
CRLB(SNR_index)=6*(1/SNR_linear)/(N*((N^2)-1))
I put it at the end of the given program but it gave me the following error:
Undefined function or variable ‘SNR_linear’.
Error in KAYSingleFrequencyEstimator1 (line 37)
CRLB(SNR_index)=6*(1/SNR_linear)/(N*((N^2)-1))
You have to set SNR_linear value (like SNR_linear=10) before using this line of code.
Dear John this code gives error in Matlab. Also the plot has three graphs but the code seems to have two of them? Can you provide the complete code?
Paste the error here and I will look into it.