- Email: [email protected]
Performance analysis of a correlation based single tone frequency estimator
SIGNAL
PROCESSING EISEVIER
Signal Processing 44 (1995) 223-231
Performance analysis of a correlation based single tone frequency estimator Peter Hhdel*,
Anders Eriksson, Torbjiirn Wigren
Systemsand Control Group, Uppsala University Box 27, S-751 03 Uppsala. Sweden Received 31 October 1994; revised 12 January 1995
Abstract This paper analyzes the frequency error variance of a low complexity single tone frequency estimator based on sample correlations of the input data. In the high SNR scenario it is analytically shown that the accuracy of a properly tuned algorithm is nearly optimal, i.e. nearly attains the Cramer-Rao lower bound. For low SNR the statistical efficiency of the algorithm is degraded, but it is analytically proven that for a large number of samples the error variance attains the lower bound for this class of estimators. Zusammenfassung Dieser Beitrag untersucht die Varianz des Frequenzfehlers eines einfachen Einfrequenz-Schatzers auf der Basis von Korrelationswerten des Eingangssignals. Fur grossen Rauschabstand wird analytisch gezeigt, dass die Genauigkeit eines korrekt eingestellen Verfahrens nahezu optimal ist: die Cramer-Rao Schranke wird nahezu erreicht. Fur geringen Rauschabstand verschlechtert sich die statistische Effizienz des Verfahrens. Es wird jedoch analytisch bewiesen, dass die Fehlervarianz bei einer grossen Anzahl von Abtastwerten fur diese Klasse von Schltzern ebenfalls die untere Schranke erreicht. RCumi: Nous analysons dans cet article la variance de l’erreur sur la frequence dun estimateur peu complexe de frequence de raie unique base sur les valeurs de correlation des donnees dent&. Dans le cas dun SNR tleve il est montrt analytiquement que la precision dun algorithme bien regli est presque optimale: elle presque atteint la borne inferieure de Cramer-Rao. Pour un SNR faible l’efficacite statistique de l’algorithme se degrade, mais il est prouve analytiquement que pour un nombre d’echantillons ClevC la variance de l’erreur atteint la borne inferieure pour cette classe d’estimateurs. Keywords: Frequency estimation;
*Corresponding
Spectral line analysis; Statistical analysis
author: Fax: + 46-18-503 611. E-mail: [email protected].
0165-1684/95/$9.50 0 1995 Elsevier Science B.V. All rights reserved SSDJ 0165-1684(95)00026-7
224
P.
Hiindel
et al. / Signal Processing 44 (1995) 223-231
1. Introduction Many papers have been published on the estimation of the frequencies of noisy sinusoidal signals, see, for example, [7] and the references therein. For the single frequency case, it is well known that the maximum likelihood (ML) estimator is given by the spectral position of the highest peak of the periodogram. The ML method can be effectively implemented using the FFT algorithm in combination with a numerical search procedure [lo]. In the high signal-to-noise (SNR) scenario other estimators are considered in [6,13], see also [l, 83. If the noisy signal at time-instant k is denoted by xk, then the method in [ 131 is based on properties of the phase angle L [&I, In order to avoid phase unwrapping the methods derived in [6] make use of L [xk$_ 1] (where * denotes complex conjugate). The frequency estimate is then computed as a (normalized or weighted) sum of phase angles. In order to reduce the number of phase angle calculations, the operations of taking the angle and summation can be interchanged, leading to a class of methods that make use of the phase angle of the sample correlation function. Such estimators have been studied in [S, 91, and more recently in [4]. Common elements of the estimators in [4,5,9] are that they are based on the set of estimated covariance elements {i1, . . . , tJ} computed from N data samples (where 1 < J < N - 1). The estimator considered in [S] is known as the linear prediction estimator and it is based on r*1alone. The estimator in [9] is based on the Jth element alone. For J = 1 the estimator coincides with the linear prediction estimator, and for J = 2N/3 the frequency error variance is minimized. The frequency error variance of the minimum variance estimator in [9] is approximately 1.125 times larger than the Cram&-Rao lower bound (CRB) at high SNR. Thus, the estimators above are not statistically efficient, meaning that the frequency error variance is (much) larger than the CRB. The recently proposed estimator in [4] makes use of the whole sequence {+k}, k = 1, . . . , J, resulting in an estimator with numerical complexity between the complexity of the estimators in [S, 91 and the ML estimator. For J = 1, the estimator in [4] coincides with the linear prediction estimator.
Explicitly, consider a single complex-valued sine wave (or cisoid) in noise, Xk=Aeiwk+t&,
k=l,...,
N,
(1)
where the complex-valued amplitude A and the angular frequency w E ( - 7c,7~) are deterministic but unknown constants. The noise okis a zero mean complex-valued white Gaussian noise with variance c*. Then, the problem of interest is the estimation of w. For this estimation problem the CRB for an unbiased estimator is given by [2, lo] CRB[O] =
6 SNRN(N’
- 1)’
(2)
where (;, denotes the estimated angular frequency and where the SNR is defined by SNR = JAl*/a*. In this paper, the estimator in [4] is studied in detail. In addition, the performance of the estimators in [S, 91 follows as special cases of the presented analysis. In [4], it is claimed that the performance of the estimator proposed there is quite near the CRB for small values of J (of order J = l-3), but this claim was supported by simulations only. In this paper, an expression for the error variance is derived that is shown to be valid for a wide range of scenarios. Especially, it is shown that the method, if properly tuned (that is J N 17N/20 which is considerably larger than the rule-of-thumb given in [4]) is nearly optimal for high SNR. For small values of J the statistical efficiency is proportional to N, and thus the loss in performance might be substantial for large values of N. On the other hand, for low SNR the performance is degraded (known as the SNR threshold effect) and the analytical expression for the error variance is shown to be close to the lower bound derived in [ 111. In [l 11, it is shown, among other things, that the error variance of an unbiased single tone estimator based on (fk}, k = 1, . . . , JaN, is bounded from below by min(var[&])
=
SNR* NJ(J + 1)(25 + 1)’
(3)
Clearly, from (3) it is seen that for low SNR (such that SNR/N<< 1) the expression (3) can be substantially larger than (2).
P. Hiindel et al. / Signal Processing 44 (1995) 223-231
2. Performance analysis
has
The considered estimator is given by, see [4], .I 6 1 1 m&m), o = J(J + 1)(25 + 1) nt=l
&P(m) = Im
6,(m) = L [d(m)],
225
6R (m) R(m) 1
[
(8)
and (4)
6R(m) =
;
(Aeiokvc-,,, + A*e-io(k-m)vk
k=m+l
B(m) =
5
x~x~_,,,,
m = l,...,
J,
+
(9)
ok&,).
k=m+l
where L [ .] denotes the phase angle of the complex-valued expression between the brackets. The integer J is a user chosen integer in the interval 1 < J < N - 1 that determines the trade-off between accuracy and numerical complexity, where, roughly speaking, an increased value of J results in an improved frequency estimate in terms of error variance. The influence of J on the performance of the estimator (4) is studied in detail in the sequel. In many situations J<< N, but this assumption is not a prerequisite for the analysis to hold. In order to analyze the method for a variety of possible scenarios, that is different values of N and SNR, a performance analysis technique similar to the ones of [3, 121 is used. The main assumption is that the estimation error is small, meaning that there is no ambiguity in the phase angle d(m). The equilibrium state of the estimator is related to a given set of signal parameters (here, A and 0). At the equilibrium state, the algorithm is entirely characterized by its set of state space parameters. The state space parameters, as defined below, are obtained for processing of the noise free signal. For the estimator (4), the state space parameters equals a.(m) = L [R(m)] = Im[logR(m)]
Inserting (9) into (8) gives (10) where gk is the real-valued colored noise gk
h
=
ok A*e-io(k-m)
k=m+
+ Aeiok
vkvk+-m +
IA12eiwm
1 ’
l,...,N.
(11)
Now, the frequency error becomes, see (4) dw =
6
i m&P(m) J(J + 1)(25 + 1) m=1
6 = J(J + 1)(2J +
(12)
In (12), the frequency error is given by a linear combination of the zero mean noise {gk}. Thus, the estimate is consistent since 60 -P 0, and the error variance var [&I = var [So] follows from var[&0] =
36 32(5 + 1)2(25 + 1)2
(5)
xkEJ+ ,=I+ECgkgfl.
and
1
(6) k=m+
1
Summing up (6) gives R(m) = (N - m)IA1’eiom.
A
(13)
1
in
The covariance sequence of gk is derived Appendix A. The result reads 1
&%.g,] = 2%
1 @k-n+.
-
8k-m,l
-
dk,l-.
+
dk,,)
(7)
In order to study the influence of the noise {ok> on the frequency estimation error 60 = ci, - o, one
1 1 2 +ZSNR
6k,ldk-,,I-“,
(14)
P. H&de1 et al. / Signal Processing
226
where &, = 1 for k = 1 and &,, = 0 for k # 1. Inserting (14) into (13) makes it possible to evaluate var [So] for a given scenario and any value of J. In order to derive a more convenient formula for var [SW], rewrite (13), using (14), as 36 u-1 + 1)2(25 + 1)2
varL-603 = P(J
+
(15)
T2),
3. Discussion Based on the results above the following remarks are in order.
The case J = 1
For J = 1 the estimator (4) coincides with the linear prediction estimator [S, 93 see also [6]. For this case, the frequency error variance is given by
where T, =-
44 (1995) 223-231
:3i&& x
.iI iv&S%
5
5 1 I=n+
k=m+
bLm,,-”
11 1 1 var [SW] = SNR(N - 1)2 +ZSNR2N-
- &-m.l
1
-
dk,I-n
+
(16)
8k.d
and
x
;
f
k=m+l
(17)
~k,dk-ml-n.
I=n+l
1 @k-m,l-n
-
bk-a,l
(21)
which coincides with the result derived in [9]. In particular, for high SNR such that SNR/N >>1, expression (21) reduces to var[&] z l/(SNR(N - 1)2) which equals the approximate variance expression derived in [6].
The case J = 2
For J = 2 the frequency error variance is given
In order to evaluate (16), one can show that (for m = l,..., N-l,n=l,...,N-1) ,=$+ 1 ,;g
1
-
dk,l-”
+
=2min(m,n,N-m,N-n).
bY 1 1 13N2-32N +20 var [SW] = - 25 SNR (N - 1)2(N - 2)2
8k.l)
(18)
5N-6 1 1 + % SNR’ (N - l)(N - 2)’
(22)
Inserting (18) into (16) gives T, =
s&i
j,
“iI ?&s-i
xmin(m,n,N-m,N-n). The sum T2 follows from a straightforward lation:
(19) calcu-
(20) The error variance is now given by (15) in combination with (19) and (20).
Comparing (22) with (21) one can note that for a moderate number of data (N >>1) the introduced averaging significantly reduces the error variance. The ratio between the first term in (22) and the first term in (21) (which are the dominating terms at high SNR) approximately equals to 0.5. The corresponding ratio for the second terms equals 0.2. This means that an increased value of J leads to a reduced error variance, where the reduction in error variance is especially pronounced for low SNR scenarios where the second term is prominent. The error variance (15) is, however, not monotonically decreasing with increasing values of J due to the reduced averaging in B(m) for m close to N.
P. Hiindel et al. / Signal Processing 44 (1995) 223-231
Efficiency at high SNR The statistical efficiency for SNR + co is given by 6N(N* - 1) = 52(5 + 1)2(23 + 1)2 xi
m=l
In (24) 3(J) is defined in (4). An expression for the error variance follows from calculations similar to those presented in the previous section. The result reads 1 min(N - J, J) var [So] = SNR J*(N - J)* 1 1 1 +ZSNR*J*(N-J).
t23j
i
n=l
The measure (23) depends on J and N. The minimization of (23) with respect to J (such that J is an integer in the interval 1 < J < N - 1) seems not to be analytically tractable. Numerical minimization shows that the minimum value of (23) is close to 1 for all values of N (for N > 10 the efficiency is nearly constant and close to 1.01) and that the minimizing argument J is well described by J ‘Y 17N/20. It is relevant to compare the high SNR efficiency of (4) with the efficiency of the estimator in [9] where the estimator is based on the Jth sample covariance only, ci, = f 8,(J).
(24)
0
5
10
221
(25)
The error variance (25) equals the result derived in [9] (see (70) there). The error variance is minimized for J = 2/3 N [9] and asymptotically (for SNR + K$ the statistical efficiency, cf. (23), equals 1.125.
Performance versus J Thousand realizations, each of length N = 24, are generated for the model (1) with A = ei4 where 4 is uniformly distributed in [0,27c] and o = 2xf with f = 0.05. The SNR is given by SNR = 3 dB. The empirical mean square-error (MSE) is shown in Fig. 1. Also shown are the corresponding theoretical variance and the CRB. It is seen from this
15
20
25
Fig. 1. The frequency error variance as a function of J for the case N = 24 and SNR = 3 dB. The diagram shows the theoretical error variance (solid line) and the empirical MSE ( + ). The CRB is indicated by the dashed line.
228
P. H&de1 et al. / Signal Processing
figure that there is an excellent agreement between the empirical results and the results predicted by the theory. Especially, note that there are two minima of the theoretical variance curve, where the minimum corresponding to the lower value of J results in a slightly inferior performance. The low complexity of the considered algorithm is illustrated as follows. For J = 20 the Matlab implementation of (4) in combination with the Matlab built-in phase unwrapping algorithm requires 2.4 x lo3 ‘flops’ per simulation run. On the other hand, finding the location of the highest peak of the periodogram (carefully implemented using the Matlab built-in fast radix-2 FFT) requires 2.8 x lo4 ‘flops’. Thus, in this particular scenario the numerical complexity of the FFT-based method is more than ten times larger than the complexity of the correlation-based method. For the FFT method, the data is zero padded to N = 512 samples in order to ensure that the distance between the FFT bins is ‘small enough’ (using the fact that the bin separation equals the square root of the CRB for N = 2x/$%@) = 426). For scenarios with higher SNR, the difference in numerical complexity between the two methods becomes more pronounced since the number of FFT bins has to be increased. In order to give some guidance on the user choice J, the minimum value of J is studied, such that var [So] < y CRB [&I,
(26)
where.y > 1 is a threshold value and the CRB is given by (2). The variance var[60] is computed from (15), (19) and (20) for y = 1.04 and y = 1.125, respectively. The results are illustrated in Fig. 2. Note that, when the SNR is reduced there is an abrupt rise in the minimum J for a fixed N. The reason for this is that there may be several minima of (15) as a function of J, which can then result in jumps when the SNR is changed. For high SNR and relatively low values of N, the required J is low and increases linearly with N. However, when the SNR is reduced, there is a region where the required J increases sharply. Thus the conclusion is that the method offers a low complexity unless the SNR is too low. Then, the complexity of the method increases significantly, in order to obtain a specified accuracy of the frequency estimate.
44 (1995) 223-231
SNR threshold When N/SNR>> 1 it follows that T2>>T1 and, thus, var [So] N
36 -J’(J + 1)‘(25 + 1)’
18 11J --Em2 = c12(3 + 1)2(2J + 1)2 SNR’ N ,,,El 1 = SN$
3 NJ(J + 1)(25 + 1)’
(27)
where N denotes an approximate equality in which only the dominant terms are retained; the higher-order terms are omitted using the assumption that N is large with respect to J. The last expression in (27) is the lower bound on the large sample variance of any unbiased single tone frequency estimator based on J _sample covariances (where & = l/N*R(k) and !_L = ;z)*), V 1, . ..) see [l 11. For J = 1, the expression (27) tends to the second term in (21). Thus this class of methods gradually degrades in performance with decreasing SNR contrary to the methods in [6, 131 where the threshold is sharp. The expression (27) is valid for J<
P. Hiindel et al. J Signal Processing 44 (1995) 223-231
229
1.04CRB
50
N (# samples)
(a)
SNR (da)
1.125 CRB
25-
20-
7 .E
15-
50S-50
03L
N (# samples) (b)
SNR (da)
Fig. 2. The minimum J giving an accuracy of (a) 1.04 CRB [G] and (b) 1.125 CRB[&], as a function of N and SNR. The region where J = 0 in the plot represents situations where no J that meets the accuracy specification can be found.
P. H&de1 et al. / Signal Processing 44 (199s) 223-231
230
1510-5
0
SNR(dB)
5
10
Fig.3. The frequency error variance as a function of SNR for N = 24. The diagram shows the theoretical error variances for J = 1 (dotted line), J = 20 (solid line), and for the estimator (24) with J = 16 (dashed dotted line). The empitical MSEs are depicted by 0, + and x , respectively. The CRB is indicated by the dashed line.
ML estimator, see [6]. In addition, the minimum SNR threshold value for the estimators studied in [6] is given by SNR = 6 dB, obtained for the, so called, weighted phase averager. The empirical performance of (24) is close to the theoretical values for SNR larger than the SNR threshold, that is given by SNR = 1 dB. The conclusion from this experiment is that among the different methods, estimator (4) appears to be the method of choice since it may outperform the other ones in terms of a lower SNR threshold. However, if a low complexity estimator is sought and a slight degradation in performance can be accepted, the estimator (24) is the method of choice. Note, however, that since (24) is based on rj only, the phase unwrapAping may cause problems, that is to ensure that @(J) x Jw belongs to the correct branch of the unit circle.
4. Conclusions
data is provided. The error variance depends on the scenario under study and on the user chosen parameter J that determines the trade-off between the accuracy of the estimate and the numerical complexity of the algorithm. With the user choice J N 17N/20 the performance is nearly optimal. Especially, for this particular choice of J the method is shown to outperform, due to its lower SNR threshold, the other related methods previously studied in the literature.
Appendix A. The covariance E [gkgl] Let zk = xk + iyk be a complex-valued zero mean noise with independent real and imaginary parts with zero means and variances 1/(2SNR), respectively. Let gk = Im[z,*_, + zk + &$_,,,]. Then Ebkgll
An explicit expression for the frequency error variance of the estimator in [4] for noisy sinusoidal
=
E[(
-
Yk-m
(-YF”
+
Yk
+
Ykh-m
-
XkYk-m.)
+ Yl + YI%-n - -%Y,-“)I
.
P. Hiindel et al. / Signal Processing 44 (1995) 223-231
=
& +
(ak-m,t-n-
dk-m,t
&(~k.,6k-m.l-.
-
-
8k,I-nh,k-m
-
8k,t-n
+ 8k,t)
dk,l-dl,k-m
+
8k.l~k-m,l-n),
from which (14) directly follows.
References Cl1 P.M. Baggenstoss and S.M. Kay, “On estimating the angle parameters of an exponential signal at high SNR”, IEEE Trans. Signal Process., Vol. 39, No. 5, May 1991, pp. 1203-1205. PI L.E. Brennan, “Angular accuracy of a phased array radar”, IRE Trans. Antennas and Propagation, May 1961, pp. 268-275. c31 A. Eriksson, Algorithms and performance analysis for spatial and spectral parameter estimation, Ph.D. thesis, Uppsala University, Uppsala, Sweden, December 1993. c41 M.P. Fitz, “Further results in the fast estimation of a single Trans. Communications, Vol. 42, frequency”, IEEE Nos. 21314,February/March/April 1994, pp. 862-864.
231
c51 L.B. Jackson and D.W. Tufts, “Frequency estimation by linear prediction”, Proc. 1978 IEEE Internat. Conf: Acoust. Speech Signal Process., Tulsa, OK, 1978, pp. 352-356. C61 S. Kay, “A fast and accurate single frequency estimator”, IEEE Truns. Acoust. Speech Signal Process., Vol. 37, No. 12, December 1989, pp. 1987-1990. c71 S.M. Kay, Modern Spectral Estimation, Theory and Application, Prentice Hall, Englewood Cliffs, NJ, 1987. 181 S.W. Lang and B.R. Musicus, “Frequency estimation from phase differences”, Proc. 1989 IEEE Internat. Conf: Acoust. Speech Signal Process., Glasgow, UK, 1989, pp. 2140-2143. c91 G.W. Lank, I.S. Reed and G.E. Pollon, “A semicoherent detection and Doppler estimation statistic”, IEEE Trans. Aerospace and Electronic Systems, Vol. AES-9, March 1973, pp. 151-165. cw DC. Rife and R.R. Boorstyn, “Single tone parameter estimation from discrete-time observations”, IEEE Trans. Inform. Theory, Vol. IT-20, No. 5, September 1974, pp. 591-598. Cl11 P. Stoica, P. Handel and T. Soderstrom, “Approximate maximum likelihood frequency estimation”, Automatica, Vol. 30, No. 1, January 1994, pp. 131-145. Cl21 P. Tichavsky and P. Handel, “Efficient tracking of multiple sinusoids with slowly varying parameters”, Proc. I993 IEEE Internat. Cona Acoust. Speech Signal Process., Minneapolis, MN, Pt. III, April 1993, pp. 368-371. Cl31 S.A. Tretter, “Estimating the frequency of a noisy sinusoid by linear regression”, IEEE Trans. Inform. Theory, Vol. IT-31, No. 6, November 1985, pp. 832-835.
PROCESSING EISEVIER
Signal Processing 44 (1995) 223-231
Performance analysis of a correlation based single tone frequency estimator Peter Hhdel*,
Anders Eriksson, Torbjiirn Wigren
Systemsand Control Group, Uppsala University Box 27, S-751 03 Uppsala. Sweden Received 31 October 1994; revised 12 January 1995
Abstract This paper analyzes the frequency error variance of a low complexity single tone frequency estimator based on sample correlations of the input data. In the high SNR scenario it is analytically shown that the accuracy of a properly tuned algorithm is nearly optimal, i.e. nearly attains the Cramer-Rao lower bound. For low SNR the statistical efficiency of the algorithm is degraded, but it is analytically proven that for a large number of samples the error variance attains the lower bound for this class of estimators. Zusammenfassung Dieser Beitrag untersucht die Varianz des Frequenzfehlers eines einfachen Einfrequenz-Schatzers auf der Basis von Korrelationswerten des Eingangssignals. Fur grossen Rauschabstand wird analytisch gezeigt, dass die Genauigkeit eines korrekt eingestellen Verfahrens nahezu optimal ist: die Cramer-Rao Schranke wird nahezu erreicht. Fur geringen Rauschabstand verschlechtert sich die statistische Effizienz des Verfahrens. Es wird jedoch analytisch bewiesen, dass die Fehlervarianz bei einer grossen Anzahl von Abtastwerten fur diese Klasse von Schltzern ebenfalls die untere Schranke erreicht. RCumi: Nous analysons dans cet article la variance de l’erreur sur la frequence dun estimateur peu complexe de frequence de raie unique base sur les valeurs de correlation des donnees dent&. Dans le cas dun SNR tleve il est montrt analytiquement que la precision dun algorithme bien regli est presque optimale: elle presque atteint la borne inferieure de Cramer-Rao. Pour un SNR faible l’efficacite statistique de l’algorithme se degrade, mais il est prouve analytiquement que pour un nombre d’echantillons ClevC la variance de l’erreur atteint la borne inferieure pour cette classe d’estimateurs. Keywords: Frequency estimation;
*Corresponding
Spectral line analysis; Statistical analysis
author: Fax: + 46-18-503 611. E-mail: [email protected].
0165-1684/95/$9.50 0 1995 Elsevier Science B.V. All rights reserved SSDJ 0165-1684(95)00026-7
224
P.
Hiindel
et al. / Signal Processing 44 (1995) 223-231
1. Introduction Many papers have been published on the estimation of the frequencies of noisy sinusoidal signals, see, for example, [7] and the references therein. For the single frequency case, it is well known that the maximum likelihood (ML) estimator is given by the spectral position of the highest peak of the periodogram. The ML method can be effectively implemented using the FFT algorithm in combination with a numerical search procedure [lo]. In the high signal-to-noise (SNR) scenario other estimators are considered in [6,13], see also [l, 83. If the noisy signal at time-instant k is denoted by xk, then the method in [ 131 is based on properties of the phase angle L [&I, In order to avoid phase unwrapping the methods derived in [6] make use of L [xk$_ 1] (where * denotes complex conjugate). The frequency estimate is then computed as a (normalized or weighted) sum of phase angles. In order to reduce the number of phase angle calculations, the operations of taking the angle and summation can be interchanged, leading to a class of methods that make use of the phase angle of the sample correlation function. Such estimators have been studied in [S, 91, and more recently in [4]. Common elements of the estimators in [4,5,9] are that they are based on the set of estimated covariance elements {i1, . . . , tJ} computed from N data samples (where 1 < J < N - 1). The estimator considered in [S] is known as the linear prediction estimator and it is based on r*1alone. The estimator in [9] is based on the Jth element alone. For J = 1 the estimator coincides with the linear prediction estimator, and for J = 2N/3 the frequency error variance is minimized. The frequency error variance of the minimum variance estimator in [9] is approximately 1.125 times larger than the Cram&-Rao lower bound (CRB) at high SNR. Thus, the estimators above are not statistically efficient, meaning that the frequency error variance is (much) larger than the CRB. The recently proposed estimator in [4] makes use of the whole sequence {+k}, k = 1, . . . , J, resulting in an estimator with numerical complexity between the complexity of the estimators in [S, 91 and the ML estimator. For J = 1, the estimator in [4] coincides with the linear prediction estimator.
Explicitly, consider a single complex-valued sine wave (or cisoid) in noise, Xk=Aeiwk+t&,
k=l,...,
N,
(1)
where the complex-valued amplitude A and the angular frequency w E ( - 7c,7~) are deterministic but unknown constants. The noise okis a zero mean complex-valued white Gaussian noise with variance c*. Then, the problem of interest is the estimation of w. For this estimation problem the CRB for an unbiased estimator is given by [2, lo] CRB[O] =
6 SNRN(N’
- 1)’
(2)
where (;, denotes the estimated angular frequency and where the SNR is defined by SNR = JAl*/a*. In this paper, the estimator in [4] is studied in detail. In addition, the performance of the estimators in [S, 91 follows as special cases of the presented analysis. In [4], it is claimed that the performance of the estimator proposed there is quite near the CRB for small values of J (of order J = l-3), but this claim was supported by simulations only. In this paper, an expression for the error variance is derived that is shown to be valid for a wide range of scenarios. Especially, it is shown that the method, if properly tuned (that is J N 17N/20 which is considerably larger than the rule-of-thumb given in [4]) is nearly optimal for high SNR. For small values of J the statistical efficiency is proportional to N, and thus the loss in performance might be substantial for large values of N. On the other hand, for low SNR the performance is degraded (known as the SNR threshold effect) and the analytical expression for the error variance is shown to be close to the lower bound derived in [ 111. In [l 11, it is shown, among other things, that the error variance of an unbiased single tone estimator based on (fk}, k = 1, . . . , JaN, is bounded from below by min(var[&])
=
SNR* NJ(J + 1)(25 + 1)’
(3)
Clearly, from (3) it is seen that for low SNR (such that SNR/N<< 1) the expression (3) can be substantially larger than (2).
P. Hiindel et al. / Signal Processing 44 (1995) 223-231
2. Performance analysis
has
The considered estimator is given by, see [4], .I 6 1 1 m&m), o = J(J + 1)(25 + 1) nt=l
&P(m) = Im
6,(m) = L [d(m)],
225
6R (m) R(m) 1
[
(8)
and (4)
6R(m) =
;
(Aeiokvc-,,, + A*e-io(k-m)vk
k=m+l
B(m) =
5
x~x~_,,,,
m = l,...,
J,
+
(9)
ok&,).
k=m+l
where L [ .] denotes the phase angle of the complex-valued expression between the brackets. The integer J is a user chosen integer in the interval 1 < J < N - 1 that determines the trade-off between accuracy and numerical complexity, where, roughly speaking, an increased value of J results in an improved frequency estimate in terms of error variance. The influence of J on the performance of the estimator (4) is studied in detail in the sequel. In many situations J<< N, but this assumption is not a prerequisite for the analysis to hold. In order to analyze the method for a variety of possible scenarios, that is different values of N and SNR, a performance analysis technique similar to the ones of [3, 121 is used. The main assumption is that the estimation error is small, meaning that there is no ambiguity in the phase angle d(m). The equilibrium state of the estimator is related to a given set of signal parameters (here, A and 0). At the equilibrium state, the algorithm is entirely characterized by its set of state space parameters. The state space parameters, as defined below, are obtained for processing of the noise free signal. For the estimator (4), the state space parameters equals a.(m) = L [R(m)] = Im[logR(m)]
Inserting (9) into (8) gives (10) where gk is the real-valued colored noise gk
h
=
ok A*e-io(k-m)
k=m+
+ Aeiok
vkvk+-m +
IA12eiwm
1 ’
l,...,N.
(11)
Now, the frequency error becomes, see (4) dw =
6
i m&P(m) J(J + 1)(25 + 1) m=1
6 = J(J + 1)(2J +
(12)
In (12), the frequency error is given by a linear combination of the zero mean noise {gk}. Thus, the estimate is consistent since 60 -P 0, and the error variance var [&I = var [So] follows from var[&0] =
36 32(5 + 1)2(25 + 1)2
(5)
xkEJ+ ,=I+ECgkgfl.
and
1
(6) k=m+
1
Summing up (6) gives R(m) = (N - m)IA1’eiom.
A
(13)
1
in
The covariance sequence of gk is derived Appendix A. The result reads 1
&%.g,] = 2%
1 @k-n+.
-
8k-m,l
-
dk,l-.
+
dk,,)
(7)
In order to study the influence of the noise {ok> on the frequency estimation error 60 = ci, - o, one
1 1 2 +ZSNR
6k,ldk-,,I-“,
(14)
P. H&de1 et al. / Signal Processing
226
where &, = 1 for k = 1 and &,, = 0 for k # 1. Inserting (14) into (13) makes it possible to evaluate var [So] for a given scenario and any value of J. In order to derive a more convenient formula for var [SW], rewrite (13), using (14), as 36 u-1 + 1)2(25 + 1)2
varL-603 = P(J
+
(15)
T2),
3. Discussion Based on the results above the following remarks are in order.
The case J = 1
For J = 1 the estimator (4) coincides with the linear prediction estimator [S, 93 see also [6]. For this case, the frequency error variance is given by
where T, =-
44 (1995) 223-231
:3i&& x
.iI iv&S%
5
5 1 I=n+
k=m+
bLm,,-”
11 1 1 var [SW] = SNR(N - 1)2 +ZSNR2N-
- &-m.l
1
-
dk,I-n
+
(16)
8k.d
and
x
;
f
k=m+l
(17)
~k,dk-ml-n.
I=n+l
1 @k-m,l-n
-
bk-a,l
(21)
which coincides with the result derived in [9]. In particular, for high SNR such that SNR/N >>1, expression (21) reduces to var[&] z l/(SNR(N - 1)2) which equals the approximate variance expression derived in [6].
The case J = 2
For J = 2 the frequency error variance is given
In order to evaluate (16), one can show that (for m = l,..., N-l,n=l,...,N-1) ,=$+ 1 ,;g
1
-
dk,l-”
+
=2min(m,n,N-m,N-n).
bY 1 1 13N2-32N +20 var [SW] = - 25 SNR (N - 1)2(N - 2)2
8k.l)
(18)
5N-6 1 1 + % SNR’ (N - l)(N - 2)’
(22)
Inserting (18) into (16) gives T, =
s&i
j,
“iI ?&s-i
xmin(m,n,N-m,N-n). The sum T2 follows from a straightforward lation:
(19) calcu-
(20) The error variance is now given by (15) in combination with (19) and (20).
Comparing (22) with (21) one can note that for a moderate number of data (N >>1) the introduced averaging significantly reduces the error variance. The ratio between the first term in (22) and the first term in (21) (which are the dominating terms at high SNR) approximately equals to 0.5. The corresponding ratio for the second terms equals 0.2. This means that an increased value of J leads to a reduced error variance, where the reduction in error variance is especially pronounced for low SNR scenarios where the second term is prominent. The error variance (15) is, however, not monotonically decreasing with increasing values of J due to the reduced averaging in B(m) for m close to N.
P. Hiindel et al. / Signal Processing 44 (1995) 223-231
Efficiency at high SNR The statistical efficiency for SNR + co is given by 6N(N* - 1) = 52(5 + 1)2(23 + 1)2 xi
m=l
In (24) 3(J) is defined in (4). An expression for the error variance follows from calculations similar to those presented in the previous section. The result reads 1 min(N - J, J) var [So] = SNR J*(N - J)* 1 1 1 +ZSNR*J*(N-J).
t23j
i
n=l
The measure (23) depends on J and N. The minimization of (23) with respect to J (such that J is an integer in the interval 1 < J < N - 1) seems not to be analytically tractable. Numerical minimization shows that the minimum value of (23) is close to 1 for all values of N (for N > 10 the efficiency is nearly constant and close to 1.01) and that the minimizing argument J is well described by J ‘Y 17N/20. It is relevant to compare the high SNR efficiency of (4) with the efficiency of the estimator in [9] where the estimator is based on the Jth sample covariance only, ci, = f 8,(J).
(24)
0
5
10
221
(25)
The error variance (25) equals the result derived in [9] (see (70) there). The error variance is minimized for J = 2/3 N [9] and asymptotically (for SNR + K$ the statistical efficiency, cf. (23), equals 1.125.
Performance versus J Thousand realizations, each of length N = 24, are generated for the model (1) with A = ei4 where 4 is uniformly distributed in [0,27c] and o = 2xf with f = 0.05. The SNR is given by SNR = 3 dB. The empirical mean square-error (MSE) is shown in Fig. 1. Also shown are the corresponding theoretical variance and the CRB. It is seen from this
15
20
25
Fig. 1. The frequency error variance as a function of J for the case N = 24 and SNR = 3 dB. The diagram shows the theoretical error variance (solid line) and the empirical MSE ( + ). The CRB is indicated by the dashed line.
228
P. H&de1 et al. / Signal Processing
figure that there is an excellent agreement between the empirical results and the results predicted by the theory. Especially, note that there are two minima of the theoretical variance curve, where the minimum corresponding to the lower value of J results in a slightly inferior performance. The low complexity of the considered algorithm is illustrated as follows. For J = 20 the Matlab implementation of (4) in combination with the Matlab built-in phase unwrapping algorithm requires 2.4 x lo3 ‘flops’ per simulation run. On the other hand, finding the location of the highest peak of the periodogram (carefully implemented using the Matlab built-in fast radix-2 FFT) requires 2.8 x lo4 ‘flops’. Thus, in this particular scenario the numerical complexity of the FFT-based method is more than ten times larger than the complexity of the correlation-based method. For the FFT method, the data is zero padded to N = 512 samples in order to ensure that the distance between the FFT bins is ‘small enough’ (using the fact that the bin separation equals the square root of the CRB for N = 2x/$%@) = 426). For scenarios with higher SNR, the difference in numerical complexity between the two methods becomes more pronounced since the number of FFT bins has to be increased. In order to give some guidance on the user choice J, the minimum value of J is studied, such that var [So] < y CRB [&I,
(26)
where.y > 1 is a threshold value and the CRB is given by (2). The variance var[60] is computed from (15), (19) and (20) for y = 1.04 and y = 1.125, respectively. The results are illustrated in Fig. 2. Note that, when the SNR is reduced there is an abrupt rise in the minimum J for a fixed N. The reason for this is that there may be several minima of (15) as a function of J, which can then result in jumps when the SNR is changed. For high SNR and relatively low values of N, the required J is low and increases linearly with N. However, when the SNR is reduced, there is a region where the required J increases sharply. Thus the conclusion is that the method offers a low complexity unless the SNR is too low. Then, the complexity of the method increases significantly, in order to obtain a specified accuracy of the frequency estimate.
44 (1995) 223-231
SNR threshold When N/SNR>> 1 it follows that T2>>T1 and, thus, var [So] N
36 -J’(J + 1)‘(25 + 1)’
18 11J --Em2 = c12(3 + 1)2(2J + 1)2 SNR’ N ,,,El 1 = SN$
3 NJ(J + 1)(25 + 1)’
(27)
where N denotes an approximate equality in which only the dominant terms are retained; the higher-order terms are omitted using the assumption that N is large with respect to J. The last expression in (27) is the lower bound on the large sample variance of any unbiased single tone frequency estimator based on J _sample covariances (where & = l/N*R(k) and !_L = ;z)*), V 1, . ..) see [l 11. For J = 1, the expression (27) tends to the second term in (21). Thus this class of methods gradually degrades in performance with decreasing SNR contrary to the methods in [6, 131 where the threshold is sharp. The expression (27) is valid for J<
P. Hiindel et al. J Signal Processing 44 (1995) 223-231
229
1.04CRB
50
N (# samples)
(a)
SNR (da)
1.125 CRB
25-
20-
7 .E
15-
50S-50
03L
N (# samples) (b)
SNR (da)
Fig. 2. The minimum J giving an accuracy of (a) 1.04 CRB [G] and (b) 1.125 CRB[&], as a function of N and SNR. The region where J = 0 in the plot represents situations where no J that meets the accuracy specification can be found.
P. H&de1 et al. / Signal Processing 44 (199s) 223-231
230
1510-5
0
SNR(dB)
5
10
Fig.3. The frequency error variance as a function of SNR for N = 24. The diagram shows the theoretical error variances for J = 1 (dotted line), J = 20 (solid line), and for the estimator (24) with J = 16 (dashed dotted line). The empitical MSEs are depicted by 0, + and x , respectively. The CRB is indicated by the dashed line.
ML estimator, see [6]. In addition, the minimum SNR threshold value for the estimators studied in [6] is given by SNR = 6 dB, obtained for the, so called, weighted phase averager. The empirical performance of (24) is close to the theoretical values for SNR larger than the SNR threshold, that is given by SNR = 1 dB. The conclusion from this experiment is that among the different methods, estimator (4) appears to be the method of choice since it may outperform the other ones in terms of a lower SNR threshold. However, if a low complexity estimator is sought and a slight degradation in performance can be accepted, the estimator (24) is the method of choice. Note, however, that since (24) is based on rj only, the phase unwrapAping may cause problems, that is to ensure that @(J) x Jw belongs to the correct branch of the unit circle.
4. Conclusions
data is provided. The error variance depends on the scenario under study and on the user chosen parameter J that determines the trade-off between the accuracy of the estimate and the numerical complexity of the algorithm. With the user choice J N 17N/20 the performance is nearly optimal. Especially, for this particular choice of J the method is shown to outperform, due to its lower SNR threshold, the other related methods previously studied in the literature.
Appendix A. The covariance E [gkgl] Let zk = xk + iyk be a complex-valued zero mean noise with independent real and imaginary parts with zero means and variances 1/(2SNR), respectively. Let gk = Im[z,*_, + zk + &$_,,,]. Then Ebkgll
An explicit expression for the frequency error variance of the estimator in [4] for noisy sinusoidal
=
E[(
-
Yk-m
(-YF”
+
Yk
+
Ykh-m
-
XkYk-m.)
+ Yl + YI%-n - -%Y,-“)I
.
P. Hiindel et al. / Signal Processing 44 (1995) 223-231
=
& +
(ak-m,t-n-
dk-m,t
&(~k.,6k-m.l-.
-
-
8k,I-nh,k-m
-
8k,t-n
+ 8k,t)
dk,l-dl,k-m
+
8k.l~k-m,l-n),
from which (14) directly follows.
References Cl1 P.M. Baggenstoss and S.M. Kay, “On estimating the angle parameters of an exponential signal at high SNR”, IEEE Trans. Signal Process., Vol. 39, No. 5, May 1991, pp. 1203-1205. PI L.E. Brennan, “Angular accuracy of a phased array radar”, IRE Trans. Antennas and Propagation, May 1961, pp. 268-275. c31 A. Eriksson, Algorithms and performance analysis for spatial and spectral parameter estimation, Ph.D. thesis, Uppsala University, Uppsala, Sweden, December 1993. c41 M.P. Fitz, “Further results in the fast estimation of a single Trans. Communications, Vol. 42, frequency”, IEEE Nos. 21314,February/March/April 1994, pp. 862-864.
231
c51 L.B. Jackson and D.W. Tufts, “Frequency estimation by linear prediction”, Proc. 1978 IEEE Internat. Conf: Acoust. Speech Signal Process., Tulsa, OK, 1978, pp. 352-356. C61 S. Kay, “A fast and accurate single frequency estimator”, IEEE Truns. Acoust. Speech Signal Process., Vol. 37, No. 12, December 1989, pp. 1987-1990. c71 S.M. Kay, Modern Spectral Estimation, Theory and Application, Prentice Hall, Englewood Cliffs, NJ, 1987. 181 S.W. Lang and B.R. Musicus, “Frequency estimation from phase differences”, Proc. 1989 IEEE Internat. Conf: Acoust. Speech Signal Process., Glasgow, UK, 1989, pp. 2140-2143. c91 G.W. Lank, I.S. Reed and G.E. Pollon, “A semicoherent detection and Doppler estimation statistic”, IEEE Trans. Aerospace and Electronic Systems, Vol. AES-9, March 1973, pp. 151-165. cw DC. Rife and R.R. Boorstyn, “Single tone parameter estimation from discrete-time observations”, IEEE Trans. Inform. Theory, Vol. IT-20, No. 5, September 1974, pp. 591-598. Cl11 P. Stoica, P. Handel and T. Soderstrom, “Approximate maximum likelihood frequency estimation”, Automatica, Vol. 30, No. 1, January 1994, pp. 131-145. Cl21 P. Tichavsky and P. Handel, “Efficient tracking of multiple sinusoids with slowly varying parameters”, Proc. I993 IEEE Internat. Cona Acoust. Speech Signal Process., Minneapolis, MN, Pt. III, April 1993, pp. 368-371. Cl31 S.A. Tretter, “Estimating the frequency of a noisy sinusoid by linear regression”, IEEE Trans. Inform. Theory, Vol. IT-31, No. 6, November 1985, pp. 832-835.