Frequency estimation of real-valued single-tone in colored noise using multiple autocorrelation lags

ARTICLE IN PRESS Signal Processing 90 (2010) 2303–2307

Contents lists available at ScienceDirect

Signal Processing journal homepage: www.elsevier.com/locate/sigpro

Fast communication

Frequency estimation of real-valued single-tone in colored noise using multiple autocorrelation lags$ Rim Elasmi-Ksibi a,, Hichem Besbes a, Roberto Lo´pez-Valcarce b, Sofiane Cherif a a b

Unite´ de Recherche TECHTRA, Sup’Com, Tunisia Department of Signal Theory and Communications, University of Vigo, Spain

a r t i c l e in fo

abstract

Article history: Received 18 May 2009 Received in revised form 14 January 2010 Accepted 27 January 2010 Available online 4 February 2010

The problem of frequency estimation of a real-valued sinusoid in colored noise is addressed. A new estimator using higher-order lags of the sample autocorrelation is proposed, and a statistical analysis is provided. Computer simulations are included to validate the theoretical development and to demonstrate that the proposed approach outperforms several existing frequency estimators in terms of estimation accuracy. & 2010 Elsevier B.V. All rights reserved.

Keywords: Frequency estimation Single sinusoid Autocorrelation

1. Introduction Detection of sinusoidal components and estimation of their frequencies in the presence of broadband noise are common problems in signal processing with a broad range of areas of application; numerous techniques have been developed for the white noise case [1,2]. The maximum likelihood (ML) method is statistically efficient, achieving the Cramer–Rao lower bound (CRLB) asymptotically, but it is computationally demanding [3]. Simpler, suboptimal frequency estimators can be obtained using the linear prediction properties of sinusoidal signals, such as the Pisarenko harmonic decomposer (PHD) [4], reformed PHD [5] and modified covariance (MC) [6] $ This work was supported by the Spanish Agency for International Cooperation (AECI) and the Tunisian Government under Project A/5405/ 06, and by the Spanish Ministry of Science and Innovation under Projects SPROACTIVE (ref. TEC2007- 68094-C02-01/TCM) and COMONSENS (CONSOLIDER-INGENIO 2010 CSD2008-00010).  Corresponding author. E-mail addresses: [email protected] (R. Elasmi-Ksibi), [email protected] (H. Besbes), [email protected] (R. Lo´pez-Valcarce), sofi[email protected] (S. Cherif).

0165-1684/$ - see front matter & 2010 Elsevier B.V. All rights reserved. doi:10.1016/j.sigpro.2010.01.025

methods; a number of statistical analysis have shown their inefficiency [7,8]. For the single sinusoid case, these simple estimators make use of just two coefficients of the sample autocorrelation. For complex tone frequency estimation, it is shown in [9–11] that improved frequency estimates can be obtained using multiple autocorrelation lags. For applications where real-time estimation is required, the use of multiple autocorrelation lags has been considered in [12] (the so-called P-estimator) and in [13] (the socalled modified PHD, MPHD). Although these different approaches result in performance improvement, they require phase unwrapping to resolve the frequency ambiguities that then appear. The phase unwrapping stage adds extra complexity to the procedure, which is clearly undesirable in applications requiring rapid frequency estimation. We focus on the real-valued case, and present a novel estimator that exploits higher lags of the sample autocorrelation. It enjoys low computational complexity since no phase unwrapping is required. The new method compares favorably to its forerunners, and it is applicable to environments with colored noise of finite memory (e.g. moving average, MA, noise models).

ARTICLE IN PRESS 2304

R. Elasmi-Ksibi et al. / Signal Processing 90 (2010) 2303–2307

This letter is organized as follows. Section 2 states the problem and presents the new estimator. A statistical analysis is given in Section 3. Section 4 shows some numerical examples, and conclusions are drawn in Section 5. 2. Novel frequency estimator Consider the problem of estimating the unknown frequency o0 2 ½0; p of a real-valued sine wave sn immersed in colored noise un. The N available samples of the observed signal, yn, are given by yn ¼ sn þun ¼ a sinðo0 n þ jÞ þ un ;

1r n r N;

ð1Þ

where a is the sinusoid amplitude and j is a random phase uniformly distributed in ½p; p½. The noise process {un} is zero-mean wide-sense stationary and Gaussian, and independent of {sn}. It is modeled as an MA(M) process, generated by passing a white Gaussian process with variance s2 through an order-M FIR filter. The normalized noise autocorrelation is denoted rk ¼ Efun unk g=s2 ; thus, rk ¼ 0 for k 4M. The power P jko . spectral density of {un} is Su ðejo Þ ¼ s2 M k ¼ M rk e The signal to noise ratio is defined as SNR 9a2 =ð2s2 Þ. The autocorrelation sequence of the observations is rk 9Efyn ynk g ¼

a

2

2

cosðko0 Þ þ s2 rk :

ð2Þ

for all k4 M þ 1:

rk ðrk1 þ rk þ 1 Þ ¼ 2 coso0

k¼p

q X

rk2 :

1

a2 Dðo0 Þ

½cosðpo0 Þ cosððp1Þo0 Þ 0    0

cosððq þ1Þo0 Þ cosðqo0 ÞT ; P with DðoÞ9 qk ¼ p cos2 ðkoÞ. It is shown in the Appendix that C

1 4 ½s B þ 2a2 Su ðejo0 ÞwwT ; N

ð9Þ

ð10Þ

with B a Toeplitz matrix with elements ½Bi;j ¼ bij , P where bl 9 M k ¼ M rk rkl , and w a vector with elements ½wi ¼ cosððp þ i2Þo0 Þ. It turns out that vT w ¼ 0, and thus, using (9) and (10) ^ 0 o0 Þ2 g  Efða^ 0 a0 Þ2 g=sin2 ðo0 Þ, one obtains and Efðo ^ 0 o0 Þ2 g  Efðo where GðoÞ ¼ 1; qp þ2g and

Gðo0 Þ ; 4NSNR D2 ðo0 Þsin2 ðo0 Þ 2

P

i2I bi Gi ð

ð11Þ

oÞ, with I ¼ f0; 1; qp; qp þ

G0 ðoÞ9cos2 ððp1ÞoÞ þ cos2 ðpoÞ þ cos2 ðqoÞ þ cos2 ððq þ 1ÞoÞ; ð12Þ

G1 ðoÞ92½cosððp1ÞoÞcosðpoÞ þ cosðqoÞcosððqþ 1ÞoÞ; ð13Þ

Gqp ðoÞ92 cosððp1ÞoÞcosððq þ1ÞoÞ;

ð14Þ

ð15Þ

ð3Þ

Gqp þ 2 ðoÞ92 cosðpoÞcosðqoÞ:

Hence, we observe that for any integers q 4 p 4M þ1, q X

v¼

Gqp þ 1 ðoÞ92½cosðpoÞcosððq þ 1ÞoÞ þ cosððp1ÞoÞcosðqoÞ;

Thus, using trigonometric relations one can show that rk1 þ rk þ 1 ¼ 2rk cosðo0 Þ;

^ One has with C the covariance matrix of r.

ð4Þ

k¼p

ð16Þ

Note that if q, p are chosen such that qp 42M, then

bqp ¼ bqp þ 1 ¼ bqp þ 2 ¼ 0. For the case of white Gaussian noise, bl ¼ dl and thus the mean square error (MSE) of the

Eq. (4) suggests a means to obtain an estimate ^ 0 Þ as a^ 0 9cosðo Pq þ r^ k þ 1 Þ r^ ðr^ k ¼ p k k1 ; ð5Þ a^ 0 ¼ Pq 2 2 k ¼ p r^ k

estimator is

where r^ k is the unbiased autocorrelation estimate

The influence of the design parameters p, q in the MSE is via the factor Gðo0 Þ=D2 ðo0 Þ. Although this term is highly oscillatory, the amplitude of the oscillations decreases as q increases. It is instructive to consider (17) for o0 ¼ p=2, for which Gðo0 Þ ¼ 2 and Dðo0 Þ  ðqp þ 1Þ=2; thus, the MSE behaves as 1/(q p +1)2. Roughly speaking, by doubling the number of lags included in the estimation process, the MSE is decreased by 6 dB. Of course, this reduction of the MSE cannot be achieved indefinitely by increasing q, since the finite sample size limits the number of autocorrelation estimates available; at the same time, the approximations used when deriving (11) become less accurate with larger q. In the following, we compare the computational complexity of the proposed approach to that of the interpolated discrete Fourier transform (DFT) [14], RPHD [5], P-estimator [12] and modified PHD (MPHD) [13]. In P-estimator and MPHD, P designs the autocorrelation sample lag used by the algorithms. Similarly, in MPHD, Pmax ¼ bðN1Þ=5c and P* defines the optimal choice of the

r^ k 9

N X 1 yn ynk : Nk n ¼ k þ 1

ð6Þ

3. Performance analysis ^ r^ p1 r^ p    r^ q þ 1 T . We can express the novel Let r9½ ^ where f is defined by (5). estimator as a^ 0 ¼ f ðrÞ, It is important to note that the following analysis holds for sufficiently large N and for high SNR level. A first-order ^ around r9½rp1 rp    rq þ 1 T Taylor expansion of f ðrÞ yields ^ a^ 0  a0 þ vT ðrrÞ;

ð7Þ

^ r^ ¼ r . with v9rf ðrÞj Thus, Efða^ 0 a0 Þ2 g  vT Cv;

ð8Þ

^ 0 o0 Þ2 g  Efðo

cos2 ððp1Þo0 Þ þ cos2 ðpo0 Þþ cos2 ðqo0 Þ þ cos2 ððq þ 1Þo0 Þ : P 2 q 4NSNR2 cos2 ðko0 Þ sin2 ðo0 Þ k¼p

ð17Þ

ARTICLE IN PRESS R. Elasmi-Ksibi et al. / Signal Processing 90 (2010) 2303–2307

2305

Table 1 Computational complexity of frequency estimators. Method

Addition

Multiplication

Division

Square root

cos1 ð Þ

tan1 ð Þ

RPHD MPHD (P = 1) MPHD (P = P*) P-estimator DFT interpolator Proposed

3N  3 4N  16 4Pmax(N + 1)  10Pmax(Pmax + 1) + P* + 4 2N  3P+ 2 2N log2(N + 54) (q  p + 1)(N + 3  (p + q)/2)  2

2N +5 2N  4 Pmax(2N + 7)  5Pmax(Pmax + 1) + 12 2N  3P + 6 2N log2(N +31) (q  p + 1)(N + 3  (p + q)/2)

1 1 Pmax 1 0 1

1 1 Pmax 1 0 0

1 1 P* +1 1 0 1

0 0 0 0 1 0

0

0 Proposed, q=7 (simulation) Proposed, q=7 (theory) Proposed, q=20 (simulation) Proposed, q=20 (theory) RPHD P−estimator, P=20 MPHD CRLB

MSE, dB

−20 −30

Proposed, q=20 (simulation) Proposed, q=20 (theory) RPHD P−estimator, P=20 MPHD CRLB

−10 −20

MSE, dB

−10

−40

−30 −40

−50

−50

−60

−60

−70

−70

0

0.2

0.4

0.6

0.8

1

w0 /π

0

0.2

0.4

0.6

0.8

1

w0 /π

Fig. 1. MSE vs. o0 . White noise case, SNR = 10 dB, N = 100, p = 2.

Fig. 2. MSE vs. o0 . Colored noise case, SNR = 10 dB, N =100, p = 6.

sample autocorrelation lag. It is seen from Table 1 that the proposed method with q = p +1, MPHD (with P =1), P-estimator and RPHD estimators involve similar operation requirements. With increased value of q, there is an increase in complexity of approximately q  p+ 1 times compared with that of q = p + 1. Compared to MPHD estimator with optimal choice of P=P*, the proposed estimator is computationally simpler since no additional computational load is required to avoid the phase wrapping problem. The comparison with the interpolated DFT estimator in terms of computational complexity depends on whether q p + 1 is smaller than 2 log2(N).

0 r kr P1 (P=20 in this example). For the new estimator (5), p= 2 was taken, and two different values of q (7 and 20) were considered. As expected, performance improves with increasing q, since more autocorrelation lags are then included in the estimation process. Good agreement with the approximate MSE (17) is observed. With q= 20, the performance of the novel estimator is very close to that of MPHD [13] with optimal choice of P=P*. The performance of the proposed method approaches the CRLB [15] for a wide range of frequencies. Next, we consider a colored noise case in which the noise coloring filter has transfer function B(z)= 1+ 0.1z  1 + 0.7z  2 +0.05 z  3 +0.3z  4. Since M =4, now p ¼ 6 4 M þ 1 is adopted for the novel estimator. The RPHD estimate is biased in the presence of noise coloring, resulting in the poor MSE behavior observed in Fig. 2. The MPHD estimator also suffers from the edge problem, since the optimal choice of the lag P will depend also on the coloring noise. The novel estimator still offers a means to trade off performance and complexity by increasing q. Fig. 3 shows the MSE of the estimators together with the interpolator DFT [14], as a function of the SNR, fixing o0 ¼ 0:4p. The sequence length is set to N= 100. This measurement time corresponds in fact to 20 cycles of the sinusoid component. Using (p,q)= (6,20), the proposed estimator outperforms the PHD, RPHD and MPHD estimators methods. In addition, the proposed algorithm has close performance to that of DFT method for SNR 4 10 dB in terms of computational load and accuracy, approaching the CRLB. However, when the

4. Simulation results Computer simulations have been carried out to evaluate the performance of the proposed estimator for a single pffiffiffireal sinusoid. The sinusoid amplitude is taken as a ¼ 2 and a random phase is used, whereas different SNRs were obtained by properly scaling the noise variance. All simulation results were averaged over 1000 independent runs. Fig. 1 shows the MSE of several estimators in a white noise setting as a function of o0 , for N=100 and SNR = 10 dB. The RPHD [5] scheme uses two autocorrelation coefficients only, performing worse than methods exploiting higher-order lags. Among these, the P-estimator [12] suffers from the so-called ’edge frequency’ problem: performance degrades for o0 close to kp=P, with

ARTICLE IN PRESS 2306

R. Elasmi-Ksibi et al. / Signal Processing 90 (2010) 2303–2307

20

0

−20

MSE, dB

T1 9

Proposed PHD RPHD P−estimator, P=30 MPHD Interpolated DFT CRLB

N N X 1 X Efui uj uik ujm g: 2 N i ¼ kþ1 j ¼ mþ1

ð20Þ

Assume w.l.o.g. that k Z m, and define a trapezoidal window wl as 8 > < Nkjlj; N þ k þ1 r l r0; 0 o l r km; wl ¼ Nk; ð21Þ > : Nmjlj; km o l rNm1:

−40

−60

Then, for large N, the following approximation holds:

−80

N N X 1 X Efsik sjm gEfui uj g N2 i ¼ k þ 1 j ¼ m þ 1

−100 −10

0

10

20

30

¼

X 1 Nm1 w Efs s gEfun un þ l g N2 l ¼ k þ 1N l nk n þ lm



1 1 X Efsnk sn þ lm gEfun un þ l g N l ¼ 1

¼

1 1 X a2 cosððlmþ kÞo0 Þ  s2 rl N l ¼ 1 2

¼

1 a2 cosððkmÞo0 ÞSu ðejo0 Þ: N 2

40

SNR, dB Fig. 3. MSE vs. SNR. Colored noise case, o0 ¼ 0:4p, N = 100, q= 20, p =6.

SNR drops below a ceratin threshold of about 5 dB, the performance of the proposed methods degrades, exhibiting a threshold behavior at a higher SNR than interpolated DFT (about 0 dB). In very low SNR settings, below  5 dB, the performance deterioration is more prominent for interpolated DFT.

The remaining terms in (19) are obtained analogously, yielding T0 

ð22Þ

Now, since {un} is Gaussian,

5. Conclusion A frequency estimator for real-valued sinusoids in noise exploiting higher-order lags of the sample autocorrelation has been presented. For MA noise, if an upper bound is available for the order of the noise coloring filter, the ‘corrupted’ low-order lags can be left out, and thus bias is avoided. Knowledge of the MA noise model parameters is not needed. The proposed method compares favorably to previous approaches and allows a tradeoff between performance and complexity by choosing the number of lags to be included in the estimation process.

Efui uj uik ujm g ¼ Efui uj gEfuik ujm g þEfui ujm gEfuik uj g þEfui uik g Efuj ujm g ; ð23Þ |fflfflfflfflfflffl{zfflfflfflfflfflffl} |fflfflfflfflfflffl{zfflfflfflfflfflffl} ¼ 0 for k 4 M ¼ 0 for m 4 M

and therefore T1 ¼

N N X 1 X Efui uj gEfuik ujm g þEfui ujm gEfuik uj g N2 i ¼ k þ 1 j ¼ m þ 1

¼

 ¼

Appendix Let k, m4 M. One has Efðr^ k rk Þðr^ m rm Þg ¼ Efr^ k r^ m g rm rk , and Efr^ k r^ m g ¼

2a2 Su ðejo0 Þcosðko0 Þcosðmo0 Þ: N

N N X 1 X Efyi yj yik yjm g N2 i ¼ k þ 1 j ¼ m þ 1

s4

Nm1 X

N 2 l ¼ k þ 1N M s4 X

N

l ¼ M M 4 X

s

N

wl ðrl rlk þ m þ rl þ m rlk Þ

rl ðrlk þ m þ rlkm Þ rl rlk þ m ;

ð24Þ

l ¼ M

where the last lkm oM.

step

follows

from

the

fact

that

References

4

¼

a

4

cosðko0 Þcosðmo0 Þ þ T0 þT1 ;

ð18Þ

where T0 and T1 are defined as T0 9

N N X 1 X ½Efsik sjm gEfui uj g þEfsi sj gEfuik ujm g 2 N i ¼ kþ1 j ¼ mþ1

þ Efsi sjm gEfuik uj g þ Efsik sj gEfui ujm g þ Efsi sik gEfuj ujm g þEfsj sjm gEfui uik g ; |fflfflfflfflfflffl{zfflfflfflfflfflffl} |fflfflfflfflfflffl{zfflfflfflfflfflffl} ¼ 0 for m 4 M

¼ 0 for k 4 M

ð19Þ

[1] P. Stoica, R.L. Moses, Spectral Analysis of Signals, Prentice-Hall, Englewood Cliffs, NJ, 2005. [2] B. Quinn, E. Hannan, The Estimation and Tracking of Frequency, Cambridge University Press, Cambridge, 2001. [3] R.J. Kenefic, A.H. Nuttall, Maximum likelihood estimation of the parameters of tone using real discrete data, IEEE J. Oceanic Eng. OE-12 (1) (1987) 279–280. [4] V.F. Pisarenko, The retrieval of harmonics by linear prediction, Geophys. J. R. Astron. Soc. 33 (1973) 347–366. [5] H.C. So, K.W. Chan, Reformulation of Pisarenko harmonic decomposition method for single tone frequency estimation, IEEE Trans. Signal Process. 52 (April 2004) 1128–1135.

ARTICLE IN PRESS R. Elasmi-Ksibi et al. / Signal Processing 90 (2010) 2303–2307

[6] L.B. Fertig, J.H. McClellan, Instantaneous frequency estimation using linear prediction with comparisons to the DESAs, IEEE Signal Process. Lett. 3 (February 1996) 54–56. [7] H. Sakai, Statistical analysis of Pisarenko’s method for sinusoidal frequency estimation, IEEE Trans. Acoust. Speech Signal Process. ASSP-32 (February 1984) 95–101. [8] A. Eriksson, P. Stoica, On statistical analysis of Pisarenko tone frequency estimator, Signal Processing 31 (April 1993) 349–353. ¨ ¨ [9] B. Volcker, P. Handel, Frequency estimation from proper sets of correlations, IEEE Trans. Signal Process. 50 (April 2002) 791–802. [10] Z. Wang, S.S. Abeysekera, Frequency estimation from multiple lags of correlations in the presence of MA colored noise, in: Proceedings of the ICASSP, 2005, pp. 693–696.

2307

[11] S.S. Abeysekera, Efficient frequency estimation using the pulse-pair method at various lags, IEEE Trans. Commun. 54 (9) (September 2006) 1542–1546. [12] R. Elasmi-Ksibi, R. Lo´pez-Valcarce, H. Besbes, S. Cherif, A family of real single-tone frequency estimators using higher-order sample covariance lags, in: Proceedings of the European Signal Processing Conference (EUSIPCO), 2008, pp. 956–959. [13] K.W.K. Lui, H.C. So, Modified Pisarenko harmonic decomposition for single-tone frequency estimation, IEEE Trans. Signal Process. 56 (7) (July 2008) 3351–3356. [14] B.G. Quinn, Estimating frequency by interpolation using Fourier coefficients, IEEE Trans. Signal Process. 42 (1) (May 1994) 1264–1268. [15] D.N. Swingler, Approximate bounds on frequency estimates for short cisoids in colored noise, IEEE Trans. Signal Process. 46 (May 1998) 1456–1458.