Closed-form unbiased frequency estimation of a noisy ... - IEEE Xplore
IEEE TRANSACTIONS ON AUTOMATIC CONTROL, VOL. 48, NO. 7, JULY 2003
small parameter " > 0, has been established. This proposition is based on the assumption of the stabilizability of the boundary-layer system. It was also shown that this connection is only one-direction valid, i.e., the controllability of the reduced-order and boundary-layer systems always yields the controllability of the original system, but not vice versa. The criterion of the impulse-free E0 -controllability of the reduced-order system is derived in the terms of an auxiliary gain matrix K (t). The invariance of this criterion to K (t) is shown. Due to the duality, similar results can be obtained for the Euclidean space observability of singularly perturbed linear time-dependent systems with multiple small delay. In this case, the assumption of the stabilizability of the boundary-layer system has to be replaced by the assumption of the detectability of this system.
REFERENCES [1] P. V. Kokotovic and A. H. Haddad, “Controllability and time-optimal control of systems with slow and fast modes,” IEEE Trans. Automat. Contr., vol. AC-20, pp. 111–113, 1975. [2] P. V. Kokotovic, H. K. Khalil, and J. O’Reilly, Singular Perturbation Methods in Control: Analysis and Design. London, U.K.: Academic, 1986. [3] H. K. Khalil, “Feedback control of nonstandard singularly perturbed systems,” IEEE Trans. Automat. Contr., vol. 34, pp. 1052–1060, Oct. 1989. [4] Y. Y. Wang, P. M. Frank, and N. E. Wu, “Near-optimal control of nonstandard singularly perturbed systems,” Automatica, vol. 30, pp. 277–292, 1994. [5] H. Krishnan and N. H. McClamroch, “On the connection between nonlinear differential-algebraic equations and singularly perturbed control systems in nonstandard form,” IEEE Trans. Automat. Contr., vol. 39, pp. 1079–1084, May 1994. [6] V. Kecman and Z. Gajic, “Optimal control and filtering for nonstandard singularly perturbed linear systems,” J. Guid. Control Dyna., vol. 22, pp. 362–365, 1999. -optimal con[7] H. Xu and K. Mizukami, “Nonstandard extension of trol for singularly perturbed systems,” in Advances in Dynamic Games and Applications. Boston, MA: Birkhauser, 2000, vol. 5, pp. 81–94. [8] E. Fridman, “A descriptor system approach to nonlinear singularly perturbed optimal control problem,” Automatica, vol. 37, pp. 543–549, 2001. [9] V. Y. Glizer and E. Fridman, “ control of linear singularly perturbed systems with small state delay,” J. Math. Anal. Applicat., vol. 250, pp. 49–85, 2000. [10] T. B. Kopeikina, “Controllability of singularly perturbed linear systems with time-lag,” Diff. Equat., vol. 25, pp. 1055–1064, 1989. [11] V. Y. Glizer, “Euclidean space controllability of singularly perturbed linear systems with state delay,” Syst. Control Lett., vol. 43, pp. 181–191, 2001. , “Controllability of singularly perturbed linear time-dependent [12] systems with small state delay,” Dyna. Control, vol. 11, pp. 261–281, 2001. [13] R. B. Vinter and R. H. Kwong, “The infinite time quadratic control problem for linear systems with state and control delays: an evolution equation approach,” SIAM J. Control Optim., vol. 19, pp. 139–153, 1981. [14] M. C. Delfour, C. McCalla, and S. K. Mitter, “Stability and the infinite-time quadratic cost problem for linear hereditary differential systems,” SIAM J. Control, vol. 13, pp. 48–88, 1975. [15] V. Y. Glizer, “Asymptotic solution of a singularly perturbed set of functional-differential equations of Riccati type encountered in the optimal control theory,” Nonlinear Diff. Equat. Applicat., vol. 5, pp. 491–515, 1998. [16] R. B. Zmood, “The Euclidean space controllability of control systems with delay,” SIAM J. Control, vol. 12, pp. 609–623, 1974. [17] A. Halanay, Differential Equations: Stability, Oscillations, Time Lags. New York: Academic, 1966. [18] I. M. Cherevko, “An estimate for the fundamental matrix of singularly perturbed differential-functional equations and some applications,” Diff. Equat., vol. 33, pp. 281–284, 1997. [19] M. C. Delfour and S. K. Mitter, “Controllability, observability and optimal feedback control of affine hereditary differential systems,” SIAM J. Control, vol. 10, pp. 298–328, 1972.
H
H
1285
[20] P. V. Kokotovic, “Applications of singular perturbation techniques to control problems,” SIAM Rev., vol. 26, pp. 501–550, 1984. [21] H. K. Khalil, “Feedback control of implicit singularly perturbed systems,” in Proc. 23rd Conf. Decision Control, Las Vegas, NV, 1984, pp. 1219–1223.
Closed-Form Unbiased Frequency Estimation of a Noisy Sinusoid Using Notch Filters Sergio M. Savaresi, S. Bittanti, and H. C. So Abstract—In this note, the problem of the frequency estimation of a sinusoid embedded in white noise is considered. The approach used herein is the minimization of the sample variance of the output of constrained notch filters fed by the noisy sinusoid. In particular, this note focuses on closed-form expressions of the frequency estimate, which can be obtained using notch filters having an all-zeros finite-inpulse response (FIR) structure. The results presented in this note are as follows. 1) It is shown that the FIR notch filters obtained from standard second-order infinite-impulse response (IIR) filters are inadequate. 2) A new second-order IIR notch filter is proposed, which provides an unbiased estimate of the frequency. 3) The FIR filter obtained from the new IIR filter provides a closed-form unbiased frequency estimate. 4) The closed-form frequency estimate obtained using the new FIR notch filter asymptotically converges toward the Pisarenko Harmonic Decomposition estimator and the Yule–Walker estimator. Index Terms—Frequency estimation, harmonic analysis, notch filters, unbiased parameter identification.
I. INTRODUCTION AND PROBLEM STATEMENT This note deals with the problem of estimating the frequency 0 of a harmonic signal s(t) = A cos( 0 t + '), given its noisy measurement y(t) = s(t)+ n(t), t = 1; 2; . . . ; N , where n(t) is a zero-mean white Gaussian noise (n W GN (0; 2 )). This problem is frequently encountered in real-world applications, especially in the fields of adaptive control and signal processing, and numerous techniques have been developed for its treatment (see, e.g., [4]–[7], [10]–[13], [15], [18], [21]). This note focuses on the class of estimation methods based on constrained notch filters (see, e.g., [8] and the references cited therein). The basic idea underlying notch-filters-based estimation techniques is the minimization, with respect to , of the loss function
J ( ) =
N t=1
"(t; )2
(1)
where "(t; ) = G(z 01 ; )y (t) is the output of a notch filter with transfer function G(z 01 ; ), fed by the measured signal y (t). The notch of G(z 01 ; ) is centered around the frequency . In general, the dependence of J ( ) on is nonlinear and nonconvex; hence, iterative quasi-Newton minimization methods must be used. Manuscript received November 5, 2002; revised January 31, 2002. Recommended by Associate Editor A. Garulli. This work was supported by MIUR project “New Methods for Identification and Adaptive Control for Industrial Systems,” and by the EU project “Nonlinear and Adaptive Control.” S. M. Savaresi and S. Bittanti are with the Dipartimento di Elettronica e Informazione, Politecnico di Milano, 20133 Milan, Italy (e-mail: [email protected]). H. C. So is with the Department of Computer Engineering and Information Technology, City University of Hong Kong, Hong Kong. Digital Object Identifier 10.1109/TAC.2003.814278
0018-9286/03$17.00 © 2003 IEEE
1286
IEEE TRANSACTIONS ON AUTOMATIC CONTROL, VOL. 48, NO. 7, JULY 2003
G1 z 01 ;
( ; 0) is that it provides a rigorously unbiased estimate of 0 . • Section IV: It is shown that the closed-form frequency estimate provided by the FIR notch filter G3 (z 01 ; ; 0), if the number N of data snapshots is large, tends to the frequency estimators provided by the Pisarenko Harmonic Decomposition (PHD) approach, and by the Yule–Walker (YW) approach.
Obviously, the most crucial design choice in a notch-based estimation technique is the selection of the structure and of the parameterization of the filter G(z 01 ; ). Usually, second-order IIR filter with a strongly constrained parameterization are used. Starting from second-order filters, simple finite-impulse response (FIR) filters or more sophisticated higher order infinite-impulse response (IIR) filters have been developed and proposed ([2], [9], [17]). Two slightly different second-order IIR notch filters are typically used in practice. They have the following expressions:
1 0 2 cos( )z01 + z02 01 G1 (z ; ; ) = 1 0 2 cos( )z01 + 2 z02 1 0 2 cos( )z01 + z02 : 01 G2 (z ; ; ) = 1 0 (1 + 2 ) cos( )z01 + 2 z02
II. NEW UNBIASING SECOND-ORDER CONSTRAINED NOTCH FILTER (2) (3)
In (2) and (3), the parameter (0 < 1) is known as the de-biasing parameter or the poles-contraction factor (note that only affects the position of the poles). In the literature, filters of this type are also known as constrained notch filter, where the term constrained refers to the fact that their structure is strongly under parameterized: the five parameters of a fully-parameterized second-order digital IIR filter are reduced to one parameter only. As a matter of fact, since is regarded as a design parameter, the only unknown parameter of (2) and (3) is the angular frequency . The main difference between (2) and (3) is that (3) provides a rigorously unbiased estimation of the frequency of a pure tone embedded in white noise, whereas (2) provides a biased estimate. It is easy to see that such bias is negligible if 1; the problem of the bias becomes severe if 1. The properties of such filters have been discussed and analyzed in a large number of works (see, e.g., [3], [14], [25], [26], and the references cited therein). Note that if the unknown frequency 0 is time-varying, and the minimization of (1) is made recursively, the estimation algorithm usually is called frequency tracker (see e.g. [22]). Notch filters are frequently used for real-time recursive frequency estimation: in the literature this problem is referred to as adaptive notch filtering (ANF). This work does not focus on ANF but the results presented herein can be straightforwardly extended to ANF as well. The goal of this note is to develop closed-form frequency estimators based on notch filters. The starting point of this work can be summarized in the following simple observations. First, note that closed-form expressions of the frequency estimator cannot be obtained if the notch filter has a IIR structure, due to the autoregressive part of the filter. Moreover, observe that a constrained FIR notch filter can be easily obtained from G1 (z 01 ; ; ) by setting = 0; unfortunately, this is not possible using G2 (z 01 ; ; ) (note that G2 (z 01 ; ; 0) is not a FIR). Finally, note that the closed-form frequency estimate obtained from G1 (z 01 ; ; 0) is severely affected by a bias error. Starting from these observations, the main results and original contributions of this note are the following. • Section II: A new second-order IIR unbiasing constrained notch filter G3 (z 01 ; ; ) is developed and analyzed. • Section III: It is shown that a closed-form frequency estimate can be obtained using the FIR filters G1 (z 01 ; ; 0) and G3 (z 01 ; ; 0); the major advantage of G3 (z 01 ; ; 0) over
As already remarked in Section I, one of the major drawbacks of the notch filter (2) (the most widely used in practice) is that it provides a biased estimation of 0 . This bias is particularly severe when 1. Starting from the cost function (1), a new unbiasing second-order IIR notch filter can be obtained as follows. • Consider the long-run (asymptotic) version of the cost function (1), namely N J
where
1 ( ) = Nlim !1 N
t=1
( )2
" t;
( ) = G(z01 ; )y(t): It is easy to see that J ( ) can be given the following expression " t;
(see [3]):
1 J ( ) = 2
+
0
(
j! ; ) 2 S (! )d! y
G e
where Sy (! ) is the power spectrum of y (t), which can be split into the power spectra of s(t) and n(t), namely
( ) = Ss (!) + Sn (!)
Sy !
and
( ) = 2
Sn !
1 (! + ) + 1 (! 0 ) : 0 0 2 2 Compute the asymptotic cost function J 1 ( ) associated with the notch filter G1 (z 01 ; ; ), by plugging in (4) the expression of the notch filter (2) and the expressions of Ss (! ) and Sn (! ) (s) (n) J 1 ( ) = J 1 ( ) + J 1 ( ); + 2 1 (n) j! 2 J 1 ( ) = G1 (e ; ; ) d! 2 0 1 A2 G (ej! ; ; ) : (s) J 1 ( ) = 0 2 2 1 (n) For the computation of J 1 ( ) (the contribution to J 1 ( ) 2
( ) = A2
Ss !
•
due to the noise) we have resorted to the Rugizka algorithm (this algorithm is based on theory of residues; see [2]). The calculus (s) of J 1 ( ) calls for cumbersome but easier computations. The (s) (n) expressions obtained for J 1 ( ) and J 1 ( ) are shown in (5a) and (b) at the bottom of the page. It is interesting to note that an
2 (cos( ) 0 cos( 0 ))2 ( ) = 1 + 4 + 42 cos2 ( ) + 42 cos2A2 ( ) 0 43 cos( ) cos( 0 ) 0 4 cos( ) cos( 0 ) 0 22 0 2 3 2 2 + 0 6 cos ( ) + + 2 cos2 ( ) + 1 (n) J 1 ( ) = : (1 + ) (2 + 2 cos( ) + 1) (2 0 2 cos( ) + 1) (s)
J1
(4)
(5a) (5b)
IEEE TRANSACTIONS ON AUTOMATIC CONTROL, VOL. 48, NO. 7, JULY 2003
Fig. 1. Shape of the function ( ; ) in the ranges
1287
2 [0; ] and 2 [0; 1].
expression very similar to (5b) is proposed in [25] (10). The two expressions (both correctly derived) do not coincide since they have a slightly different meaning. • Observe that the bias in the frequency estimate obtained using (n) G1 (z 01 ; ; ) is due to the fact that J 1 ( ) is a function of . (n) This dependence of J 1 ( ) on has the effect of moving the minimum of J 1 ( ) away from 0 (note, instead, that 0 is the (s) minimum of J 1 ( )). (s) Now, observe that the minimum of J 1 ( ) does not change (s) if J 1 ( ) is multiplied by a strictly positive function of and , say ( ; ) (obviously for 2 [0; ) and 2 [0; 1)). This is (s) due by the presence of the factor (cos( ) 0 cos( 0 )) in J 1 ( ), which is null if = 0 . Consider then the following function ( ; ):
( ; )
(1+ ) (2 +2 cos( )+1)(2 0 2 cos( )+1) : = (1+ 2 ) (3 + 2 0 6 cos2 ( )+ +1+2cos2 ( ))
(6)
Note that such function is the square-root of the inverse of (n) J 1 ( ) (but for the coefficient 2 = ), multiplied by (1 + 2 )). A new unbiasing filter G3 (z 01 ; ; ) can be obtained from G1 (z 01 ; ; ) and (6) as follows:
G3 (z 01 ; ; ) = ( ; )G1 (z 01 ; ; )
1 0 2 cos( )z01 + z02 : = ( ; ) 1 0 2 cos( )z01 + 2 z02
(7)
Due to the fact that G3 (z 01 ; ; ) is simply obtained by multiplying G1 (z 01 ; ; ) by ( ; ), some remarks on the shape of ( ; ) are due (see Fig. 1 where ( ; ) is plotted in the ranges 2 [0; ] and 2 [0; 1]). • ( ; ) is not null in the ranges 2 [0; ] and 2 [0; 1]; this can be easily seen from (6). This guarantees the well posedness of the optimization problem based on the cost function (1). • Note that ( ; 1) = 1, whereas ( ; 0) strongly differs from 1; this is expected since ( ; ) is a sort of “de-biasing factor” of G1 (z 01 ; ; ). Therefore, ( ; ) leaves G1 (z 01 ; ; ) almost unchanged if is close to 1, whereas ( ; ) provides a strong correction to G1 (z 01 ; ; ) for small values of . • Note that (=2; ) = 1 8 2 [0; 1], and that ( ; ) is symmetric with respect to = =2 in the range 2 [0; ]. This is consistent with a peculiar feature of G1 (z 01 ; ; ): it provides an unbiased estimate 2 [0; 1] if and only if
0 = =2 ([3]). In order to get a complete understanding of the differences between the three second-order constrained notch filters G1 (z 01 ; ; ), G2 (z 01 ; ; ), and G3 (z 01 ; ; ), it is interesting to compare the corresponding asymptotic cost functions J 1 ( ), J 2 ( ), and J 3 ( ), respectively. The closed-form expressions of J 1 ( ) has already been computed (s) (n) in (5). Following the same procedure, J 2 ( ) = J 2 ( ) + J 2 ( ) (s) (n) and J 3 ( ) = J 3 ( ) + J 3 ( ) can be obtained as (8a), (8b), (9a), and (9b), shown at the bottom of the next page By comparing J 1 ( ), J 2 ( ), and J 3 ( ) it is apparent that J 2 ( ) and J 3 ( ) have the minimum exactly at 0 (unbiased estimate), since
1288
IEEE TRANSACTIONS ON AUTOMATIC CONTROL, VOL. 48, NO. 7, JULY 2003
Fig. 2. Shape of the cost functions J
( ), J ( ), and J ( ) (A = 1, = 1, = =6, = 0:5).
(n) (n) J 2 ( ) and J 3 ( ) do not depend on . On the contrary, as already (n) observed, J 1 ( ) provides a biased estimate since J 1 ( ) is dependent. The shapes of these three different cost functions can be better appreciated from Fig. 2, where J 1 ( ), J 2 ( ), and J 3 ( ) are displayed for A = 1, 2 = 1 (signal-to-noise ratio (SNR)= 0:5), 0 = =6, and = 0:5. To conclude this section, it is worth remarking that the new filter G3 (z 01 ; ; ) merges the two main appealing features of G1 (z 01 ; ; ) and G2 (z 01 ; ; ): similarly to G1 (z 01 ; ; ), an FIR filter can be obtained from G3 (z 01 ; ; ) by simply using = 0; similarly to G2 (z 01 ; ; ), G3 (z 01 ; ; ) provides an unbiased estimate of 0 8 2 [0; 1). These features will be fully exploited in the following section, in order to obtain closed-form frequency estimates based on FIR notch filters.
III. CLOSED-FORM FREQUENCY ESTIMATION VIA NOTCH FILTERS A closed-form notch-based frequency estimate cannot be obtained if the filter has a IIR structure. Consider the FIR filters obtained by simply setting = 0 in (2) and in (7) (it has been already observed that setting = 0 in G2 (z 01 ; ; ) does not yield a FIR filter), namely
G1 (z 01 ; ; 0) = 1 0 2 cos( )z 01 + z 02 G3 (z 01 ; ; 0) =
1 01 02 2 cos2 ( )+1 1 0 2 cos( )z + z
:
Using such filters, closed-form frequency estimators from the data can be obtained as follows.
(s) J 2 ( ) =
2A2 (cos( ) 0 cos( 0 ))2 1 + cos ( ) 0 2 cos( ) cos( 0 ) 0 2 + + cos2 ( ) + 22 cos2 ( ) 0 42 cos( ) cos( 0 ) 0 2 cos( ) cos( 0 ) + 42 cos2 ( 0 ) 4
2
4
2
4
(8a) (n) J 2 ( ) = (s) J 3 ( ) = (n) J 3 ( ) =
2 (1 + 2 )
(8b)
2A2 (cos( ) 0 cos( 0 ))2 2 ( ; ) 1 + + 4 cos ( 0 ) + 42 cos2 ( ) 0 43 cos( ) cos( 0 ) 0 4 cos( ) cos( 0 ) 0 22 4
2 (1 + 2 ) :
2
2
(9a) (9b)
IEEE TRANSACTIONS ON AUTOMATIC CONTROL, VOL. 48, NO. 7, JULY 2003
Fig. 3.
Asymptotic bias error of (11) (A
= 1,
1289
= 1).
A. Closed-Form Frequency Estimator Obtained Using G1 (z 01 ; ; 0)
is obtained: 2
^ 1 = N!1 ! arccos AA 2 + 2 cos( 0 ) :
Consider the following cost function, obtained by plugging in (1) the FIR notch filter G1 (z 01 ; ; 0):
J1 ( ) =
N
t=1
(y(t) 0 2 cos ( )y(t 0 1) + y(t 0 2))2
and differentiate J1 ( ) with respect to :
dJ1 ( ) d
=2
N
t=1
[(y(t) 0 2 cos( )y(t 0 1) + y(t 0 2))
2 (2 sin( )y(t 0 1))] :
(10)
Note that (10) is quadratic with respect to cos( ); hence, by solving dJ1 ( )=d = 0 with respect to cos( ), it is easy to see that the following holds:
N
cos( ) = t=1 y(t 0N1) (y(t) + 2y(t 0 2)) : t=1 y(t 0 1) N y(t 0 1) (y(t) + y(t 0 2)) t=1 : N y(t 0 1)2 t=1
From (12), it is apparent that the frequency estimate is affected by a severe bias; note that the bias is null in the (trivial and unrealistic) case of zero noise ( 2 = 0); it grows as the SNR decreases. To get a quantitative idea of this bias error, in Fig. 3 the asymptotic bias error in the case of SNR= 0:5 (A = 1, 2 = 1) is displayed. As expected, the bias is null if 0 = =2; it is maximum for 0 = 0 or 0 = . Note that, in average, the bias error is huge; hence, the frequency estimate (11) is of no use in practice. B. Closed-Form Frequency Estimator Obtained Using G3 (z 01 ; ; 0) Consider the following cost function, obtained by plugging in (1) the FIR notch filter G3 (z 01 ; ; 0):
N
(y(t) 0 2 cos( )y(t 0 1) + y(t 0 2))2 2 cos2 ( ) + 1 t=1 and differentiate J3 ( ) with respect to ; (13), as shown at the bottom J3 ( ) =
The closed-form frequency estimator, therefore, is given by
^ 1 = arccos
(12)
(11)
^ 1 tends to the minimum of the As the number N of data grows,
asymptotic cost function (5) (in the special case of = 0). After some cumbersome computation, the following asymptotic expression of (11)
of the next page holds. Consider now the problem of solving dJ3 ( )=d = 0 with respect to . After some manipulation the following expression is obtained:
N
t=1
(sin( )(y(t) 0 2 cos( )y(t 0 1) + y(t 0 2))
2 (y(t 0 1) + cos( )y(t) + cos( )y(t 0 2))) = 0:
(14)
1290
IEEE TRANSACTIONS ON AUTOMATIC CONTROL, VOL. 48, NO. 7, JULY 2003
Equation (14) admits a trivial solution: sin( ) = 0. Assuming that
0 6= f0; g, the following quadratic form (with respect to cos( )) can be obtained from (14):
2 N t=1
N t=1
y(t 0 1) (y(t)+ y(t 0 2))
cos2 ( )+
t=1
Given a signal y (t) = s(t) + n(t), where s(t) = A cos( 0 t + '), W GN (0; 2 ), the autocorrelation coefficients of order 1 and 2, say r1 and r2 , respectively, are given by
n
r1 = E [y(t)y(t 0 1)] = A2 cos( 0 ) : (17) r2 = E [y(t)y(t 0 2)] = A2 cos(2 0 ) By eliminating the parameter A in (17), the following equation is
2y(t 0 1)2 0 (y(t)+ y(t 0 2))2 cos( )0 N
A. YW Approach
obtained:
y(t 0 1)(y(t)+ y(t 0 2))
= 0:
(15)
From (15), a closed-form frequency estimator can be computed. It has the expression shown in (16) at the bottom of the page. ^ 3 tends to the minimum of the As the number N of data grows,
asymptotic cost function (8) (in the special case of = 0)
^ 3 N! !1 0 :
Thus, the new filter G3 (z 01 ; ; 0) provides a simple closed-form unbiased estimate of 0 . Interestingly, (16) is closely related to the method given in [19] and [20], even if it the derivation of this result is completely different. We conclude this section by briefly discussing the problem of using a finite number of data. This note mainly deals with asymptotic results, even if (16) is used for N finite. To get a rough and preliminary indication on its behavior when N is small, in Fig. 4 the average frequency estimate and error variance obtained using (16) for N = 100, N = 1:000, N = 10:000, and N = 100:000 are displayed. The results in Fig. 4 have been obtained using 100 different uncorrelated (obviously for the noisy part only) realizations of the signal y . It is apparent that both the bias and the variance errors rapidly decrease when N gets large. A comparison—for finite (and low) values of N —between this estimation algorithm and other estimation algorithms goes out of the scope of this note and might be the subject of future work.
2r1 cos2 ( 0 ) 0 r2 cos( 0 ) 0 r1 = 0: Its solution with respect to 0 provides the YW frequency estimator,
given by
2 2
^ YW = arccos r2 + 4rr2 + 8r1 : 1 B. PHD Approach
Given a zero-mean stationary signal y (t), its autocorrelation matrix of order three is given by
r0 r1 r2 R = r1 r0 r1 r0 = E y(t)2 r2 r1 r0 r1 = E [y(t)y(t 0 1)] r2 = E [y(t)y(t 0 2)] : The eigenvector associated with the smallest eigenvalue of R has the following form:
p 1 0 r + 2rr +8r
T
1 :
In the literature, other closed-form frequency estimators for harmonic signals in white noise have been proposed and analyzed. Two celebrated estimators are the YW estimator, and the “ PHD estimator (see, e.g., [1], [7], [16], [19], [20], [23], and [24]). In this section, they will be briefly recalled and compared with the asymptotic version of the notch-based estimator (16).
(19)
Pisarenko (see [16] and the analysis proposed in [7]) has proven that, if y (t) = s(t)+ n(t) (s(t) = A cos( 0 t + '), n W GN (0; 2 ), the smallest eigenvalue of R must have the following simple expression:
[ 1 02 cos( 0 ) 1 ]T :
IV. RELATED FREQUENCY-ESTIMATION METHODS
(18)
(20)
By comparing (19) and (20), the PHD frequency estimator is obtained
2 2
^ PHD = arccos r2 + 4rr2 + 8r1 : 1
(21)
dJ3 ( ) = d
N
2 (y(t) 0 2 cos( )y(t 0 1) + y(t 0 2))((2 + cos(2 ))2 sin( )y(t 0 1) + sin(2 )(y(t) 0 2 cos( )y(t 0 1) + y(t 0 2))) : (2 + cos(2 ))2 t=1 (13)
^ 3 = arccos
1
0
N 2y(t 0 1)2 0 (y(t)+ y(t 0 2))2 t=1
+ 4
N 2y(t 0 1)2 0 (y(t)+ y(t 0 2))2 2 +8 t=1 N y(t 0 1)(y(t)+ y(t 0 2)) t=1
N y(t 0 1)(y(t)+ y(t 0 2)) 2 t=1
:
(16)
IEEE TRANSACTIONS ON AUTOMATIC CONTROL, VOL. 48, NO. 7, JULY 2003
Fig. 4.
1291
Average frequency estimate and error variance for different values of N (A
Interestingly enough, the YW and PHD approaches provide exactly the same results. This has been recently proven and discussed in [23] and [24]. Consider now the notch-based closed-form estimator (16). The following result holds. Proposition 1: Given a signal y (t) = s(t) + n(t), where s(t) = A cos( 0 t + '), n W GN (0; 2 ), the notch-filter based estimator ^ YW and ^ PHD , namely (16) asymptotically converges toward
^
Nlim !1 3 = arccos
r2 +
r22 + 8r12 4r1
Proof: If N is large, the following hold:
2 2 ^ 3 = arccos r2 + r2 + 8r1 lim
N !1 4r1 r1 = E [y(t)y(t 0 1)] r2 = E [y(t)y(t 0 2)] :
^ YW , From a theoretical point of view the fact that (asymptotically)
^ PHD and ^ 3 are exactly the same, is particularly interesting: it shows
V. CONCLUSION AND FUTURE WORK
1 N y(t 0 1) (y(t) + y(t 0 2)) lim N !1 N t=1
1 N 2y(t)y(t 0 1) = Nlim !1 N t=1
(22a)
N !1 N t=1
= 02E [y(t)y(t 0 2)] :
By plugging in (16) the asymptotic expressions (22), it is easy to see that
the equivalence of three classical approaches which have been independently conceived and developed following three completely different paths.
r1 = E [y(t)y(t 0 1)] r2 = E [y(t)y(t 0 2)] :
= 2E [y(t)y(t 0 1)] ; N lim 1 2y(t 0 1)2 0 (y(t) + y(t 0 2))2
= 1, = 1, = =6).
(22b)
In this work, a notch FIR filter which provides and unbiased closed-form frequency estimate of harmonic signals in white noise has been proposed. This estimator has been proven to converge asymptotically to the well-known YW and PHD estimators. These three equivalent estimators are very appealing since they admit a very simple closed-form expression starting from a set of measured data. However, their main flaw is the sensitivity to the noise: they guarantee an unbiased estimate if the noise is white; when the noise is colored, the bias error can be large. In order to try to overcome this pitfall, probably the most interesting rationale is that proposed by Quinn and Fernandes in [15]: the idea is
1292
IEEE TRANSACTIONS ON AUTOMATIC CONTROL, VOL. 48, NO. 7, JULY 2003
to prefilter the data with a frequency enhancer, in an iterative fashion; at each step the frequency enhancer is centered around the frequency estimated at the previous step. This approach seems to fit perfectly to simple closed-form estimators. The reduction of the noise-sensitivity of FIR-based frequency estimators using frequency enhancers is currently the subject of further investigation.
[25] Y. Xiao, Y. Takeshita, and K. Shida, “Steady-state analysis of a plain gradient algorithm for a second-order adaptive IIR notch filter with constrained poles and zeros,” IEEE Trans. Circuits Systems, vol. 48, pp. 733–740, July 2001. , “Tracking properties of a gradient-based second-order adaptive [26] IIR notch filter with constrained poles and zeros,” IEEE Trans. Signal Processing, vol. 50, pp. 878–888, Apr. 2002.
REFERENCES [1] E. Anarim and Y. Istefanopulos, “Statistical analysis of Pisarenko type tone frequency estimator,” Signal Processing, vol. 24, pp. 291–298, 1991. [2] K. J. Åström, Introduction to Stochastic Control Theory. New York: Academic, 1970. [3] S. Bittanti, M. Campi, and S. M. Savaresi, “Unbiased estimation of a sinusoid in noise via adapted notch filters,” Automatica, vol. 33, no. 2, pp. 209–215, 1997. [4] S. Bittanti and G. Picci, Eds., Identification, Adaptation, Learning—The Science of Learning Models from Data. Berlin, Germany: SpringerVerlag, 1996, Computer and Systems Sciences. [5] S. Bittanti and S. M. Savaresi, “On the parametrization and design of an extended Kalman filter frequency tracker,” IEEE Trans. Automat. Contr., vol. 45, pp. 1718–1724, Sept. 2000. [6] B. Boashash, “Estimating and interpreting the instantaneous frequency of a signal,” Proc. IEEE, vol. 80, pp. 520–568, Apr. 1992. [7] A. Eriksson and P. Stoica, “On statistical analysis of Pisarenko tone frequency estimator,” Signal Processing, vol. 31, pp. 349–353, 1993. [8] P. Händel and A. Nehorai, “Tracking analysis of an adaptive notch filter with constrained poles and zeros,” IEEE Trans. Signal Processing, vol. 42, pp. 281–291, Feb. 1994. [9] P. Händel, P. Tichavsky, and S. M. Savaresi, “Large error recovery for a class of frequency tracking algorithms,” Int. J. Adapt. Control Signal Processing, vol. 12, pp. 417–436, 1998. [10] L. Hsu, R. Ortega, and G. Damm, “A globally convergent frequency estimator,” IEEE Trans. Automat. Contr., vol. 44, pp. 698–713, Apr. 1999. [11] S. M. Kay, Modern Spectral Estimation: Theory and Applications. Upper Saddle River, NJ: Prentice-Hall, 1988. [12] B. La Scala and R. Bitmead, “Design of an extended Kalman filter frequency tracker,” IEEE Trans. Signal Processing, vol. 44, pp. 739–742, Mar. 1996. [13] B. La Scala, R. Bitmead, and B. G. Quinn, “An extended Kalman filter frequency tracker for high-noise environments,” IEEE Trans. Signal Processing, vol. 44, pp. 431–434, Feb. 1996. [14] G. Li, “A stable and efficient adaptive notch filter for direct frequency estimation,” IEEE Trans. Signal Processing, vol. 45, pp. 2001–2009, Aug. 1997. [15] B. G. Quinn and J. M. Fernandes, “A fast efficient technique for the estimation of frequency,” Biometrika, vol. 78, pp. 489–497, 1991. [16] V. F. Pisarenko, “The retrieval of harmonics from a covariance function,” J. Roy. Astr., vol. 33, pp. 374–376, 1973. [17] S. M. Savaresi, “Funnel filters: a new class of filters for frequency estimation of harmonic signals,” Automatica, vol. 33, no. 9, pp. 1711–1718, 1997. [18] S. M. Savaresi, R. Bitmead, and W. Dunstan, “Nonlinear system identification using closed-loop data with no external excitation: the case of a lean combustion process,” Int. J. Control, vol. 74, no. 18, pp. 1796–1806, 2001. [19] H. C. So, “A closed form frequency estimator for a noisy sinusoid,” in 45th IEEE Midwest Symp. Circuits Systems, Tulsa, OK, 2002. [20] H. C. So and S. K. Ip, “A novel frequency estimator and its comparative performances for short record lengths,” in 11th Eur. Signal Processing Conf., Toulouse, France, 2002. [21] P. Stoica, “List of references on spectral line analysis,” Signal Processing, vol. 31, pp. 329–340, 1992. [22] P. Strobach, “Single section least squares adaptive notch filter,” IEEE Trans. Signal Processing, vol. 43, pp. 2007–2010, Aug. 1995. [23] Y. Xiao and Y. Tadokoro, “On Pisarenko and constrained Yule-Walker estimators of tone frequency,” IEICE Trans., vol. E77-A, no. 8, pp. 1404–1406, 1994. [24] , “Statistical analysis of a simple constrained high-order Yule-Walker tone frequency estimator,” IEICE Trans., vol. E78-A, no. 10, pp. 1415–1418, 1995.
Comments on “A Robust State Observer Scheme” M. Boutayeb and M. Darouach Abstract—In this note, we point out that the sufficient conditions given in a recent paper to assure convergence of the proposed robust state observer are incomplete. Index Terms—Robust observer design, stability analysis, uncertain systems.
I. INTRODUCTION In [2], a robust state observer scheme was proposed for uncertain linear systems. The main result and theorem may be summarized as follows. Consider the linear map
A0 + 1A x + Bu 61 yx_ == Cx + Du
1 0
where A represents the model uncertainty, assumed to be bounded by > with k A k2 < . The authors propose the following observer dynamics:
1
Bu + H (y 0 y^) + 62 xy^^_ == CA0x^x^++Du: We are given the uncertain system 61 and the state ob-
Theorem: server 2 . If we set
6
2 T = x^T x^ P 01 C T r; 2r r
where P
r = y 0 y^
> 0 satisfies the following matrix inequality:
(A0 0 HC )T P + P (A0 0 HC ) + 2P 2 + 2 I < 0 and H , such that (A0 0 HC ) is stable, then the state error vector e = < x 0 x^ converges to zero. Manuscript received June 7, 2002; revised September 27, 2002, January 10, 2003, and January 17, 2003. Recommended by Associate Editor M. E. Valcher. M. Boutayeb is with the LSIIT-CNRS UMR 7005, University of Louis Pasteur, 67400 Illkirch, France (e-mail: [email protected]). M. Darouach is with the CRAN-CNRS UMR 7039, University of Henri Poincaré, IUT of Longwy 54400, France (e-mail: [email protected]). Digital Object Identifier 10.1109/TAC.2003.812795
0018-9286/03$17.00 © 2003 IEEE
small parameter " > 0, has been established. This proposition is based on the assumption of the stabilizability of the boundary-layer system. It was also shown that this connection is only one-direction valid, i.e., the controllability of the reduced-order and boundary-layer systems always yields the controllability of the original system, but not vice versa. The criterion of the impulse-free E0 -controllability of the reduced-order system is derived in the terms of an auxiliary gain matrix K (t). The invariance of this criterion to K (t) is shown. Due to the duality, similar results can be obtained for the Euclidean space observability of singularly perturbed linear time-dependent systems with multiple small delay. In this case, the assumption of the stabilizability of the boundary-layer system has to be replaced by the assumption of the detectability of this system.
REFERENCES [1] P. V. Kokotovic and A. H. Haddad, “Controllability and time-optimal control of systems with slow and fast modes,” IEEE Trans. Automat. Contr., vol. AC-20, pp. 111–113, 1975. [2] P. V. Kokotovic, H. K. Khalil, and J. O’Reilly, Singular Perturbation Methods in Control: Analysis and Design. London, U.K.: Academic, 1986. [3] H. K. Khalil, “Feedback control of nonstandard singularly perturbed systems,” IEEE Trans. Automat. Contr., vol. 34, pp. 1052–1060, Oct. 1989. [4] Y. Y. Wang, P. M. Frank, and N. E. Wu, “Near-optimal control of nonstandard singularly perturbed systems,” Automatica, vol. 30, pp. 277–292, 1994. [5] H. Krishnan and N. H. McClamroch, “On the connection between nonlinear differential-algebraic equations and singularly perturbed control systems in nonstandard form,” IEEE Trans. Automat. Contr., vol. 39, pp. 1079–1084, May 1994. [6] V. Kecman and Z. Gajic, “Optimal control and filtering for nonstandard singularly perturbed linear systems,” J. Guid. Control Dyna., vol. 22, pp. 362–365, 1999. -optimal con[7] H. Xu and K. Mizukami, “Nonstandard extension of trol for singularly perturbed systems,” in Advances in Dynamic Games and Applications. Boston, MA: Birkhauser, 2000, vol. 5, pp. 81–94. [8] E. Fridman, “A descriptor system approach to nonlinear singularly perturbed optimal control problem,” Automatica, vol. 37, pp. 543–549, 2001. [9] V. Y. Glizer and E. Fridman, “ control of linear singularly perturbed systems with small state delay,” J. Math. Anal. Applicat., vol. 250, pp. 49–85, 2000. [10] T. B. Kopeikina, “Controllability of singularly perturbed linear systems with time-lag,” Diff. Equat., vol. 25, pp. 1055–1064, 1989. [11] V. Y. Glizer, “Euclidean space controllability of singularly perturbed linear systems with state delay,” Syst. Control Lett., vol. 43, pp. 181–191, 2001. , “Controllability of singularly perturbed linear time-dependent [12] systems with small state delay,” Dyna. Control, vol. 11, pp. 261–281, 2001. [13] R. B. Vinter and R. H. Kwong, “The infinite time quadratic control problem for linear systems with state and control delays: an evolution equation approach,” SIAM J. Control Optim., vol. 19, pp. 139–153, 1981. [14] M. C. Delfour, C. McCalla, and S. K. Mitter, “Stability and the infinite-time quadratic cost problem for linear hereditary differential systems,” SIAM J. Control, vol. 13, pp. 48–88, 1975. [15] V. Y. Glizer, “Asymptotic solution of a singularly perturbed set of functional-differential equations of Riccati type encountered in the optimal control theory,” Nonlinear Diff. Equat. Applicat., vol. 5, pp. 491–515, 1998. [16] R. B. Zmood, “The Euclidean space controllability of control systems with delay,” SIAM J. Control, vol. 12, pp. 609–623, 1974. [17] A. Halanay, Differential Equations: Stability, Oscillations, Time Lags. New York: Academic, 1966. [18] I. M. Cherevko, “An estimate for the fundamental matrix of singularly perturbed differential-functional equations and some applications,” Diff. Equat., vol. 33, pp. 281–284, 1997. [19] M. C. Delfour and S. K. Mitter, “Controllability, observability and optimal feedback control of affine hereditary differential systems,” SIAM J. Control, vol. 10, pp. 298–328, 1972.
H
H
1285
[20] P. V. Kokotovic, “Applications of singular perturbation techniques to control problems,” SIAM Rev., vol. 26, pp. 501–550, 1984. [21] H. K. Khalil, “Feedback control of implicit singularly perturbed systems,” in Proc. 23rd Conf. Decision Control, Las Vegas, NV, 1984, pp. 1219–1223.
Closed-Form Unbiased Frequency Estimation of a Noisy Sinusoid Using Notch Filters Sergio M. Savaresi, S. Bittanti, and H. C. So Abstract—In this note, the problem of the frequency estimation of a sinusoid embedded in white noise is considered. The approach used herein is the minimization of the sample variance of the output of constrained notch filters fed by the noisy sinusoid. In particular, this note focuses on closed-form expressions of the frequency estimate, which can be obtained using notch filters having an all-zeros finite-inpulse response (FIR) structure. The results presented in this note are as follows. 1) It is shown that the FIR notch filters obtained from standard second-order infinite-impulse response (IIR) filters are inadequate. 2) A new second-order IIR notch filter is proposed, which provides an unbiased estimate of the frequency. 3) The FIR filter obtained from the new IIR filter provides a closed-form unbiased frequency estimate. 4) The closed-form frequency estimate obtained using the new FIR notch filter asymptotically converges toward the Pisarenko Harmonic Decomposition estimator and the Yule–Walker estimator. Index Terms—Frequency estimation, harmonic analysis, notch filters, unbiased parameter identification.
I. INTRODUCTION AND PROBLEM STATEMENT This note deals with the problem of estimating the frequency 0 of a harmonic signal s(t) = A cos( 0 t + '), given its noisy measurement y(t) = s(t)+ n(t), t = 1; 2; . . . ; N , where n(t) is a zero-mean white Gaussian noise (n W GN (0; 2 )). This problem is frequently encountered in real-world applications, especially in the fields of adaptive control and signal processing, and numerous techniques have been developed for its treatment (see, e.g., [4]–[7], [10]–[13], [15], [18], [21]). This note focuses on the class of estimation methods based on constrained notch filters (see, e.g., [8] and the references cited therein). The basic idea underlying notch-filters-based estimation techniques is the minimization, with respect to , of the loss function
J ( ) =
N t=1
"(t; )2
(1)
where "(t; ) = G(z 01 ; )y (t) is the output of a notch filter with transfer function G(z 01 ; ), fed by the measured signal y (t). The notch of G(z 01 ; ) is centered around the frequency . In general, the dependence of J ( ) on is nonlinear and nonconvex; hence, iterative quasi-Newton minimization methods must be used. Manuscript received November 5, 2002; revised January 31, 2002. Recommended by Associate Editor A. Garulli. This work was supported by MIUR project “New Methods for Identification and Adaptive Control for Industrial Systems,” and by the EU project “Nonlinear and Adaptive Control.” S. M. Savaresi and S. Bittanti are with the Dipartimento di Elettronica e Informazione, Politecnico di Milano, 20133 Milan, Italy (e-mail: [email protected]). H. C. So is with the Department of Computer Engineering and Information Technology, City University of Hong Kong, Hong Kong. Digital Object Identifier 10.1109/TAC.2003.814278
0018-9286/03$17.00 © 2003 IEEE
1286
IEEE TRANSACTIONS ON AUTOMATIC CONTROL, VOL. 48, NO. 7, JULY 2003
G1 z 01 ;
( ; 0) is that it provides a rigorously unbiased estimate of 0 . • Section IV: It is shown that the closed-form frequency estimate provided by the FIR notch filter G3 (z 01 ; ; 0), if the number N of data snapshots is large, tends to the frequency estimators provided by the Pisarenko Harmonic Decomposition (PHD) approach, and by the Yule–Walker (YW) approach.
Obviously, the most crucial design choice in a notch-based estimation technique is the selection of the structure and of the parameterization of the filter G(z 01 ; ). Usually, second-order IIR filter with a strongly constrained parameterization are used. Starting from second-order filters, simple finite-impulse response (FIR) filters or more sophisticated higher order infinite-impulse response (IIR) filters have been developed and proposed ([2], [9], [17]). Two slightly different second-order IIR notch filters are typically used in practice. They have the following expressions:
1 0 2 cos( )z01 + z02 01 G1 (z ; ; ) = 1 0 2 cos( )z01 + 2 z02 1 0 2 cos( )z01 + z02 : 01 G2 (z ; ; ) = 1 0 (1 + 2 ) cos( )z01 + 2 z02
II. NEW UNBIASING SECOND-ORDER CONSTRAINED NOTCH FILTER (2) (3)
In (2) and (3), the parameter (0 < 1) is known as the de-biasing parameter or the poles-contraction factor (note that only affects the position of the poles). In the literature, filters of this type are also known as constrained notch filter, where the term constrained refers to the fact that their structure is strongly under parameterized: the five parameters of a fully-parameterized second-order digital IIR filter are reduced to one parameter only. As a matter of fact, since is regarded as a design parameter, the only unknown parameter of (2) and (3) is the angular frequency . The main difference between (2) and (3) is that (3) provides a rigorously unbiased estimation of the frequency of a pure tone embedded in white noise, whereas (2) provides a biased estimate. It is easy to see that such bias is negligible if 1; the problem of the bias becomes severe if 1. The properties of such filters have been discussed and analyzed in a large number of works (see, e.g., [3], [14], [25], [26], and the references cited therein). Note that if the unknown frequency 0 is time-varying, and the minimization of (1) is made recursively, the estimation algorithm usually is called frequency tracker (see e.g. [22]). Notch filters are frequently used for real-time recursive frequency estimation: in the literature this problem is referred to as adaptive notch filtering (ANF). This work does not focus on ANF but the results presented herein can be straightforwardly extended to ANF as well. The goal of this note is to develop closed-form frequency estimators based on notch filters. The starting point of this work can be summarized in the following simple observations. First, note that closed-form expressions of the frequency estimator cannot be obtained if the notch filter has a IIR structure, due to the autoregressive part of the filter. Moreover, observe that a constrained FIR notch filter can be easily obtained from G1 (z 01 ; ; ) by setting = 0; unfortunately, this is not possible using G2 (z 01 ; ; ) (note that G2 (z 01 ; ; 0) is not a FIR). Finally, note that the closed-form frequency estimate obtained from G1 (z 01 ; ; 0) is severely affected by a bias error. Starting from these observations, the main results and original contributions of this note are the following. • Section II: A new second-order IIR unbiasing constrained notch filter G3 (z 01 ; ; ) is developed and analyzed. • Section III: It is shown that a closed-form frequency estimate can be obtained using the FIR filters G1 (z 01 ; ; 0) and G3 (z 01 ; ; 0); the major advantage of G3 (z 01 ; ; 0) over
As already remarked in Section I, one of the major drawbacks of the notch filter (2) (the most widely used in practice) is that it provides a biased estimation of 0 . This bias is particularly severe when 1. Starting from the cost function (1), a new unbiasing second-order IIR notch filter can be obtained as follows. • Consider the long-run (asymptotic) version of the cost function (1), namely N J
where
1 ( ) = Nlim !1 N
t=1
( )2
" t;
( ) = G(z01 ; )y(t): It is easy to see that J ( ) can be given the following expression " t;
(see [3]):
1 J ( ) = 2
+
0
(
j! ; ) 2 S (! )d! y
G e
where Sy (! ) is the power spectrum of y (t), which can be split into the power spectra of s(t) and n(t), namely
( ) = Ss (!) + Sn (!)
Sy !
and
( ) = 2
Sn !
1 (! + ) + 1 (! 0 ) : 0 0 2 2 Compute the asymptotic cost function J 1 ( ) associated with the notch filter G1 (z 01 ; ; ), by plugging in (4) the expression of the notch filter (2) and the expressions of Ss (! ) and Sn (! ) (s) (n) J 1 ( ) = J 1 ( ) + J 1 ( ); + 2 1 (n) j! 2 J 1 ( ) = G1 (e ; ; ) d! 2 0 1 A2 G (ej! ; ; ) : (s) J 1 ( ) = 0 2 2 1 (n) For the computation of J 1 ( ) (the contribution to J 1 ( ) 2
( ) = A2
Ss !
•
due to the noise) we have resorted to the Rugizka algorithm (this algorithm is based on theory of residues; see [2]). The calculus (s) of J 1 ( ) calls for cumbersome but easier computations. The (s) (n) expressions obtained for J 1 ( ) and J 1 ( ) are shown in (5a) and (b) at the bottom of the page. It is interesting to note that an
2 (cos( ) 0 cos( 0 ))2 ( ) = 1 + 4 + 42 cos2 ( ) + 42 cos2A2 ( ) 0 43 cos( ) cos( 0 ) 0 4 cos( ) cos( 0 ) 0 22 0 2 3 2 2 + 0 6 cos ( ) + + 2 cos2 ( ) + 1 (n) J 1 ( ) = : (1 + ) (2 + 2 cos( ) + 1) (2 0 2 cos( ) + 1) (s)
J1
(4)
(5a) (5b)
IEEE TRANSACTIONS ON AUTOMATIC CONTROL, VOL. 48, NO. 7, JULY 2003
Fig. 1. Shape of the function ( ; ) in the ranges
1287
2 [0; ] and 2 [0; 1].
expression very similar to (5b) is proposed in [25] (10). The two expressions (both correctly derived) do not coincide since they have a slightly different meaning. • Observe that the bias in the frequency estimate obtained using (n) G1 (z 01 ; ; ) is due to the fact that J 1 ( ) is a function of . (n) This dependence of J 1 ( ) on has the effect of moving the minimum of J 1 ( ) away from 0 (note, instead, that 0 is the (s) minimum of J 1 ( )). (s) Now, observe that the minimum of J 1 ( ) does not change (s) if J 1 ( ) is multiplied by a strictly positive function of and , say ( ; ) (obviously for 2 [0; ) and 2 [0; 1)). This is (s) due by the presence of the factor (cos( ) 0 cos( 0 )) in J 1 ( ), which is null if = 0 . Consider then the following function ( ; ):
( ; )
(1+ ) (2 +2 cos( )+1)(2 0 2 cos( )+1) : = (1+ 2 ) (3 + 2 0 6 cos2 ( )+ +1+2cos2 ( ))
(6)
Note that such function is the square-root of the inverse of (n) J 1 ( ) (but for the coefficient 2 = ), multiplied by (1 + 2 )). A new unbiasing filter G3 (z 01 ; ; ) can be obtained from G1 (z 01 ; ; ) and (6) as follows:
G3 (z 01 ; ; ) = ( ; )G1 (z 01 ; ; )
1 0 2 cos( )z01 + z02 : = ( ; ) 1 0 2 cos( )z01 + 2 z02
(7)
Due to the fact that G3 (z 01 ; ; ) is simply obtained by multiplying G1 (z 01 ; ; ) by ( ; ), some remarks on the shape of ( ; ) are due (see Fig. 1 where ( ; ) is plotted in the ranges 2 [0; ] and 2 [0; 1]). • ( ; ) is not null in the ranges 2 [0; ] and 2 [0; 1]; this can be easily seen from (6). This guarantees the well posedness of the optimization problem based on the cost function (1). • Note that ( ; 1) = 1, whereas ( ; 0) strongly differs from 1; this is expected since ( ; ) is a sort of “de-biasing factor” of G1 (z 01 ; ; ). Therefore, ( ; ) leaves G1 (z 01 ; ; ) almost unchanged if is close to 1, whereas ( ; ) provides a strong correction to G1 (z 01 ; ; ) for small values of . • Note that (=2; ) = 1 8 2 [0; 1], and that ( ; ) is symmetric with respect to = =2 in the range 2 [0; ]. This is consistent with a peculiar feature of G1 (z 01 ; ; ): it provides an unbiased estimate 2 [0; 1] if and only if
0 = =2 ([3]). In order to get a complete understanding of the differences between the three second-order constrained notch filters G1 (z 01 ; ; ), G2 (z 01 ; ; ), and G3 (z 01 ; ; ), it is interesting to compare the corresponding asymptotic cost functions J 1 ( ), J 2 ( ), and J 3 ( ), respectively. The closed-form expressions of J 1 ( ) has already been computed (s) (n) in (5). Following the same procedure, J 2 ( ) = J 2 ( ) + J 2 ( ) (s) (n) and J 3 ( ) = J 3 ( ) + J 3 ( ) can be obtained as (8a), (8b), (9a), and (9b), shown at the bottom of the next page By comparing J 1 ( ), J 2 ( ), and J 3 ( ) it is apparent that J 2 ( ) and J 3 ( ) have the minimum exactly at 0 (unbiased estimate), since
1288
IEEE TRANSACTIONS ON AUTOMATIC CONTROL, VOL. 48, NO. 7, JULY 2003
Fig. 2. Shape of the cost functions J
( ), J ( ), and J ( ) (A = 1, = 1, = =6, = 0:5).
(n) (n) J 2 ( ) and J 3 ( ) do not depend on . On the contrary, as already (n) observed, J 1 ( ) provides a biased estimate since J 1 ( ) is dependent. The shapes of these three different cost functions can be better appreciated from Fig. 2, where J 1 ( ), J 2 ( ), and J 3 ( ) are displayed for A = 1, 2 = 1 (signal-to-noise ratio (SNR)= 0:5), 0 = =6, and = 0:5. To conclude this section, it is worth remarking that the new filter G3 (z 01 ; ; ) merges the two main appealing features of G1 (z 01 ; ; ) and G2 (z 01 ; ; ): similarly to G1 (z 01 ; ; ), an FIR filter can be obtained from G3 (z 01 ; ; ) by simply using = 0; similarly to G2 (z 01 ; ; ), G3 (z 01 ; ; ) provides an unbiased estimate of 0 8 2 [0; 1). These features will be fully exploited in the following section, in order to obtain closed-form frequency estimates based on FIR notch filters.
III. CLOSED-FORM FREQUENCY ESTIMATION VIA NOTCH FILTERS A closed-form notch-based frequency estimate cannot be obtained if the filter has a IIR structure. Consider the FIR filters obtained by simply setting = 0 in (2) and in (7) (it has been already observed that setting = 0 in G2 (z 01 ; ; ) does not yield a FIR filter), namely
G1 (z 01 ; ; 0) = 1 0 2 cos( )z 01 + z 02 G3 (z 01 ; ; 0) =
1 01 02 2 cos2 ( )+1 1 0 2 cos( )z + z
:
Using such filters, closed-form frequency estimators from the data can be obtained as follows.
(s) J 2 ( ) =
2A2 (cos( ) 0 cos( 0 ))2 1 + cos ( ) 0 2 cos( ) cos( 0 ) 0 2 + + cos2 ( ) + 22 cos2 ( ) 0 42 cos( ) cos( 0 ) 0 2 cos( ) cos( 0 ) + 42 cos2 ( 0 ) 4
2
4
2
4
(8a) (n) J 2 ( ) = (s) J 3 ( ) = (n) J 3 ( ) =
2 (1 + 2 )
(8b)
2A2 (cos( ) 0 cos( 0 ))2 2 ( ; ) 1 + + 4 cos ( 0 ) + 42 cos2 ( ) 0 43 cos( ) cos( 0 ) 0 4 cos( ) cos( 0 ) 0 22 4
2 (1 + 2 ) :
2
2
(9a) (9b)
IEEE TRANSACTIONS ON AUTOMATIC CONTROL, VOL. 48, NO. 7, JULY 2003
Fig. 3.
Asymptotic bias error of (11) (A
= 1,
1289
= 1).
A. Closed-Form Frequency Estimator Obtained Using G1 (z 01 ; ; 0)
is obtained: 2
^ 1 = N!1 ! arccos AA 2 + 2 cos( 0 ) :
Consider the following cost function, obtained by plugging in (1) the FIR notch filter G1 (z 01 ; ; 0):
J1 ( ) =
N
t=1
(y(t) 0 2 cos ( )y(t 0 1) + y(t 0 2))2
and differentiate J1 ( ) with respect to :
dJ1 ( ) d
=2
N
t=1
[(y(t) 0 2 cos( )y(t 0 1) + y(t 0 2))
2 (2 sin( )y(t 0 1))] :
(10)
Note that (10) is quadratic with respect to cos( ); hence, by solving dJ1 ( )=d = 0 with respect to cos( ), it is easy to see that the following holds:
N
cos( ) = t=1 y(t 0N1) (y(t) + 2y(t 0 2)) : t=1 y(t 0 1) N y(t 0 1) (y(t) + y(t 0 2)) t=1 : N y(t 0 1)2 t=1
From (12), it is apparent that the frequency estimate is affected by a severe bias; note that the bias is null in the (trivial and unrealistic) case of zero noise ( 2 = 0); it grows as the SNR decreases. To get a quantitative idea of this bias error, in Fig. 3 the asymptotic bias error in the case of SNR= 0:5 (A = 1, 2 = 1) is displayed. As expected, the bias is null if 0 = =2; it is maximum for 0 = 0 or 0 = . Note that, in average, the bias error is huge; hence, the frequency estimate (11) is of no use in practice. B. Closed-Form Frequency Estimator Obtained Using G3 (z 01 ; ; 0) Consider the following cost function, obtained by plugging in (1) the FIR notch filter G3 (z 01 ; ; 0):
N
(y(t) 0 2 cos( )y(t 0 1) + y(t 0 2))2 2 cos2 ( ) + 1 t=1 and differentiate J3 ( ) with respect to ; (13), as shown at the bottom J3 ( ) =
The closed-form frequency estimator, therefore, is given by
^ 1 = arccos
(12)
(11)
^ 1 tends to the minimum of the As the number N of data grows,
asymptotic cost function (5) (in the special case of = 0). After some cumbersome computation, the following asymptotic expression of (11)
of the next page holds. Consider now the problem of solving dJ3 ( )=d = 0 with respect to . After some manipulation the following expression is obtained:
N
t=1
(sin( )(y(t) 0 2 cos( )y(t 0 1) + y(t 0 2))
2 (y(t 0 1) + cos( )y(t) + cos( )y(t 0 2))) = 0:
(14)
1290
IEEE TRANSACTIONS ON AUTOMATIC CONTROL, VOL. 48, NO. 7, JULY 2003
Equation (14) admits a trivial solution: sin( ) = 0. Assuming that
0 6= f0; g, the following quadratic form (with respect to cos( )) can be obtained from (14):
2 N t=1
N t=1
y(t 0 1) (y(t)+ y(t 0 2))
cos2 ( )+
t=1
Given a signal y (t) = s(t) + n(t), where s(t) = A cos( 0 t + '), W GN (0; 2 ), the autocorrelation coefficients of order 1 and 2, say r1 and r2 , respectively, are given by
n
r1 = E [y(t)y(t 0 1)] = A2 cos( 0 ) : (17) r2 = E [y(t)y(t 0 2)] = A2 cos(2 0 ) By eliminating the parameter A in (17), the following equation is
2y(t 0 1)2 0 (y(t)+ y(t 0 2))2 cos( )0 N
A. YW Approach
obtained:
y(t 0 1)(y(t)+ y(t 0 2))
= 0:
(15)
From (15), a closed-form frequency estimator can be computed. It has the expression shown in (16) at the bottom of the page. ^ 3 tends to the minimum of the As the number N of data grows,
asymptotic cost function (8) (in the special case of = 0)
^ 3 N! !1 0 :
Thus, the new filter G3 (z 01 ; ; 0) provides a simple closed-form unbiased estimate of 0 . Interestingly, (16) is closely related to the method given in [19] and [20], even if it the derivation of this result is completely different. We conclude this section by briefly discussing the problem of using a finite number of data. This note mainly deals with asymptotic results, even if (16) is used for N finite. To get a rough and preliminary indication on its behavior when N is small, in Fig. 4 the average frequency estimate and error variance obtained using (16) for N = 100, N = 1:000, N = 10:000, and N = 100:000 are displayed. The results in Fig. 4 have been obtained using 100 different uncorrelated (obviously for the noisy part only) realizations of the signal y . It is apparent that both the bias and the variance errors rapidly decrease when N gets large. A comparison—for finite (and low) values of N —between this estimation algorithm and other estimation algorithms goes out of the scope of this note and might be the subject of future work.
2r1 cos2 ( 0 ) 0 r2 cos( 0 ) 0 r1 = 0: Its solution with respect to 0 provides the YW frequency estimator,
given by
2 2
^ YW = arccos r2 + 4rr2 + 8r1 : 1 B. PHD Approach
Given a zero-mean stationary signal y (t), its autocorrelation matrix of order three is given by
r0 r1 r2 R = r1 r0 r1 r0 = E y(t)2 r2 r1 r0 r1 = E [y(t)y(t 0 1)] r2 = E [y(t)y(t 0 2)] : The eigenvector associated with the smallest eigenvalue of R has the following form:
p 1 0 r + 2rr +8r
T
1 :
In the literature, other closed-form frequency estimators for harmonic signals in white noise have been proposed and analyzed. Two celebrated estimators are the YW estimator, and the “ PHD estimator (see, e.g., [1], [7], [16], [19], [20], [23], and [24]). In this section, they will be briefly recalled and compared with the asymptotic version of the notch-based estimator (16).
(19)
Pisarenko (see [16] and the analysis proposed in [7]) has proven that, if y (t) = s(t)+ n(t) (s(t) = A cos( 0 t + '), n W GN (0; 2 ), the smallest eigenvalue of R must have the following simple expression:
[ 1 02 cos( 0 ) 1 ]T :
IV. RELATED FREQUENCY-ESTIMATION METHODS
(18)
(20)
By comparing (19) and (20), the PHD frequency estimator is obtained
2 2
^ PHD = arccos r2 + 4rr2 + 8r1 : 1
(21)
dJ3 ( ) = d
N
2 (y(t) 0 2 cos( )y(t 0 1) + y(t 0 2))((2 + cos(2 ))2 sin( )y(t 0 1) + sin(2 )(y(t) 0 2 cos( )y(t 0 1) + y(t 0 2))) : (2 + cos(2 ))2 t=1 (13)
^ 3 = arccos
1
0
N 2y(t 0 1)2 0 (y(t)+ y(t 0 2))2 t=1
+ 4
N 2y(t 0 1)2 0 (y(t)+ y(t 0 2))2 2 +8 t=1 N y(t 0 1)(y(t)+ y(t 0 2)) t=1
N y(t 0 1)(y(t)+ y(t 0 2)) 2 t=1
:
(16)
IEEE TRANSACTIONS ON AUTOMATIC CONTROL, VOL. 48, NO. 7, JULY 2003
Fig. 4.
1291
Average frequency estimate and error variance for different values of N (A
Interestingly enough, the YW and PHD approaches provide exactly the same results. This has been recently proven and discussed in [23] and [24]. Consider now the notch-based closed-form estimator (16). The following result holds. Proposition 1: Given a signal y (t) = s(t) + n(t), where s(t) = A cos( 0 t + '), n W GN (0; 2 ), the notch-filter based estimator ^ YW and ^ PHD , namely (16) asymptotically converges toward
^
Nlim !1 3 = arccos
r2 +
r22 + 8r12 4r1
Proof: If N is large, the following hold:
2 2 ^ 3 = arccos r2 + r2 + 8r1 lim
N !1 4r1 r1 = E [y(t)y(t 0 1)] r2 = E [y(t)y(t 0 2)] :
^ YW , From a theoretical point of view the fact that (asymptotically)
^ PHD and ^ 3 are exactly the same, is particularly interesting: it shows
V. CONCLUSION AND FUTURE WORK
1 N y(t 0 1) (y(t) + y(t 0 2)) lim N !1 N t=1
1 N 2y(t)y(t 0 1) = Nlim !1 N t=1
(22a)
N !1 N t=1
= 02E [y(t)y(t 0 2)] :
By plugging in (16) the asymptotic expressions (22), it is easy to see that
the equivalence of three classical approaches which have been independently conceived and developed following three completely different paths.
r1 = E [y(t)y(t 0 1)] r2 = E [y(t)y(t 0 2)] :
= 2E [y(t)y(t 0 1)] ; N lim 1 2y(t 0 1)2 0 (y(t) + y(t 0 2))2
= 1, = 1, = =6).
(22b)
In this work, a notch FIR filter which provides and unbiased closed-form frequency estimate of harmonic signals in white noise has been proposed. This estimator has been proven to converge asymptotically to the well-known YW and PHD estimators. These three equivalent estimators are very appealing since they admit a very simple closed-form expression starting from a set of measured data. However, their main flaw is the sensitivity to the noise: they guarantee an unbiased estimate if the noise is white; when the noise is colored, the bias error can be large. In order to try to overcome this pitfall, probably the most interesting rationale is that proposed by Quinn and Fernandes in [15]: the idea is
1292
IEEE TRANSACTIONS ON AUTOMATIC CONTROL, VOL. 48, NO. 7, JULY 2003
to prefilter the data with a frequency enhancer, in an iterative fashion; at each step the frequency enhancer is centered around the frequency estimated at the previous step. This approach seems to fit perfectly to simple closed-form estimators. The reduction of the noise-sensitivity of FIR-based frequency estimators using frequency enhancers is currently the subject of further investigation.
[25] Y. Xiao, Y. Takeshita, and K. Shida, “Steady-state analysis of a plain gradient algorithm for a second-order adaptive IIR notch filter with constrained poles and zeros,” IEEE Trans. Circuits Systems, vol. 48, pp. 733–740, July 2001. , “Tracking properties of a gradient-based second-order adaptive [26] IIR notch filter with constrained poles and zeros,” IEEE Trans. Signal Processing, vol. 50, pp. 878–888, Apr. 2002.
REFERENCES [1] E. Anarim and Y. Istefanopulos, “Statistical analysis of Pisarenko type tone frequency estimator,” Signal Processing, vol. 24, pp. 291–298, 1991. [2] K. J. Åström, Introduction to Stochastic Control Theory. New York: Academic, 1970. [3] S. Bittanti, M. Campi, and S. M. Savaresi, “Unbiased estimation of a sinusoid in noise via adapted notch filters,” Automatica, vol. 33, no. 2, pp. 209–215, 1997. [4] S. Bittanti and G. Picci, Eds., Identification, Adaptation, Learning—The Science of Learning Models from Data. Berlin, Germany: SpringerVerlag, 1996, Computer and Systems Sciences. [5] S. Bittanti and S. M. Savaresi, “On the parametrization and design of an extended Kalman filter frequency tracker,” IEEE Trans. Automat. Contr., vol. 45, pp. 1718–1724, Sept. 2000. [6] B. Boashash, “Estimating and interpreting the instantaneous frequency of a signal,” Proc. IEEE, vol. 80, pp. 520–568, Apr. 1992. [7] A. Eriksson and P. Stoica, “On statistical analysis of Pisarenko tone frequency estimator,” Signal Processing, vol. 31, pp. 349–353, 1993. [8] P. Händel and A. Nehorai, “Tracking analysis of an adaptive notch filter with constrained poles and zeros,” IEEE Trans. Signal Processing, vol. 42, pp. 281–291, Feb. 1994. [9] P. Händel, P. Tichavsky, and S. M. Savaresi, “Large error recovery for a class of frequency tracking algorithms,” Int. J. Adapt. Control Signal Processing, vol. 12, pp. 417–436, 1998. [10] L. Hsu, R. Ortega, and G. Damm, “A globally convergent frequency estimator,” IEEE Trans. Automat. Contr., vol. 44, pp. 698–713, Apr. 1999. [11] S. M. Kay, Modern Spectral Estimation: Theory and Applications. Upper Saddle River, NJ: Prentice-Hall, 1988. [12] B. La Scala and R. Bitmead, “Design of an extended Kalman filter frequency tracker,” IEEE Trans. Signal Processing, vol. 44, pp. 739–742, Mar. 1996. [13] B. La Scala, R. Bitmead, and B. G. Quinn, “An extended Kalman filter frequency tracker for high-noise environments,” IEEE Trans. Signal Processing, vol. 44, pp. 431–434, Feb. 1996. [14] G. Li, “A stable and efficient adaptive notch filter for direct frequency estimation,” IEEE Trans. Signal Processing, vol. 45, pp. 2001–2009, Aug. 1997. [15] B. G. Quinn and J. M. Fernandes, “A fast efficient technique for the estimation of frequency,” Biometrika, vol. 78, pp. 489–497, 1991. [16] V. F. Pisarenko, “The retrieval of harmonics from a covariance function,” J. Roy. Astr., vol. 33, pp. 374–376, 1973. [17] S. M. Savaresi, “Funnel filters: a new class of filters for frequency estimation of harmonic signals,” Automatica, vol. 33, no. 9, pp. 1711–1718, 1997. [18] S. M. Savaresi, R. Bitmead, and W. Dunstan, “Nonlinear system identification using closed-loop data with no external excitation: the case of a lean combustion process,” Int. J. Control, vol. 74, no. 18, pp. 1796–1806, 2001. [19] H. C. So, “A closed form frequency estimator for a noisy sinusoid,” in 45th IEEE Midwest Symp. Circuits Systems, Tulsa, OK, 2002. [20] H. C. So and S. K. Ip, “A novel frequency estimator and its comparative performances for short record lengths,” in 11th Eur. Signal Processing Conf., Toulouse, France, 2002. [21] P. Stoica, “List of references on spectral line analysis,” Signal Processing, vol. 31, pp. 329–340, 1992. [22] P. Strobach, “Single section least squares adaptive notch filter,” IEEE Trans. Signal Processing, vol. 43, pp. 2007–2010, Aug. 1995. [23] Y. Xiao and Y. Tadokoro, “On Pisarenko and constrained Yule-Walker estimators of tone frequency,” IEICE Trans., vol. E77-A, no. 8, pp. 1404–1406, 1994. [24] , “Statistical analysis of a simple constrained high-order Yule-Walker tone frequency estimator,” IEICE Trans., vol. E78-A, no. 10, pp. 1415–1418, 1995.
Comments on “A Robust State Observer Scheme” M. Boutayeb and M. Darouach Abstract—In this note, we point out that the sufficient conditions given in a recent paper to assure convergence of the proposed robust state observer are incomplete. Index Terms—Robust observer design, stability analysis, uncertain systems.
I. INTRODUCTION In [2], a robust state observer scheme was proposed for uncertain linear systems. The main result and theorem may be summarized as follows. Consider the linear map
A0 + 1A x + Bu 61 yx_ == Cx + Du
1 0
where A represents the model uncertainty, assumed to be bounded by > with k A k2 < . The authors propose the following observer dynamics:
1
Bu + H (y 0 y^) + 62 xy^^_ == CA0x^x^++Du: We are given the uncertain system 61 and the state ob-
Theorem: server 2 . If we set
6
2 T = x^T x^ P 01 C T r; 2r r
where P
r = y 0 y^
> 0 satisfies the following matrix inequality:
(A0 0 HC )T P + P (A0 0 HC ) + 2P 2 + 2 I < 0 and H , such that (A0 0 HC ) is stable, then the state error vector e = < x 0 x^ converges to zero. Manuscript received June 7, 2002; revised September 27, 2002, January 10, 2003, and January 17, 2003. Recommended by Associate Editor M. E. Valcher. M. Boutayeb is with the LSIIT-CNRS UMR 7005, University of Louis Pasteur, 67400 Illkirch, France (e-mail: [email protected]). M. Darouach is with the CRAN-CNRS UMR 7039, University of Henri Poincaré, IUT of Longwy 54400, France (e-mail: [email protected]). Digital Object Identifier 10.1109/TAC.2003.812795
0018-9286/03$17.00 © 2003 IEEE
