Algebraic parameter estimation of a biased sinusoidal waveform signal from noisy data

16th IFAC Symposium on System Identification The International Federation of Automatic Control Brussels, Belgium. July 11-13, 2012

Algebraic parameter estimation of a biased sinusoidal waveform signal from noisy data Rosane Ushirobira ∗ Wilfrid Perruquetti ∗∗ Mamadou Mboup ∗∗∗ Michel Fliess ∗∗∗∗ ∗

Institut de Math´ematiques de Bourgogne (CNRS), Universit´e de Bourgogne & Non-A Inria Lille - Nord Europe, France (e-mail: [email protected]) ∗∗ ´ Universit´e Lille Nord de France & Ecole Centrale de Lille & Non-A Inria Lille - Nord Europe & LAGIS (CNRS), France (e-mail: [email protected]) ∗∗∗ CReSTIC, Universit´e de Reims Champagne Ardenne, France & Non-A Inria Lille - Nord Europe (e-mail: [email protected]) ∗∗∗∗ ´ LIX (CNRS), Ecole polytechnique, 91128 Palaiseau (e-mail: [email protected]) & AL.I.E.N., 24-30 rue Lionnois, BP 60120, 54003 Nancy, France, www.alien-sas.com Abstract: The amplitude, frequency and phase of a biased and noisy sum of two complex exponential sinusoidal signals are estimated via new algebraic techniques providing a robust estimation within a fraction of the signal period. The methods that are popular today do not seem able to achieve such performances. The efficiency of our approach is illustrated by several computer simulations. Keywords: Parameter identification, differential algebra, sinusoidal waveform. 1. INTRODUCTION The parameter estimation of a biased sinusoidal signal in a noisy environment is an important issue occurring in many practical engineering problems, e.g. the signal demodulation in communications, the regulation of electronic converters power, the circadian rhythm of biological cells and the modal identification for flexible structures (see Trapero et al. [2007b]). Many different resolution methods have been developed, such as linear or nonlinear regression, subspace methods (Haykin [1991], Roy et al. [1989], Kahn et al. [1992]), the extended Kalman filter (Bittanti et al. [2000]), the notches filter (Regalia et al. [1995]), or alternatively, the use of techniques borrowed from adaptive nonlinear control (Hsu et al. [1999], Mojiri et al. [2004]). However, the robust parameter estimation in a fraction of the time signal, in the presence of noise and of an unknown constant bias, is not yet fully solved. This paper draws its inspiration from the algebraic analysis of Fliess et al. [2003, 2010, 2008], Fliess [2008], Mboup [2009]. In addition to numerical simulations found in these papers, we refer to Neves et al. [2006], Trapero et al. [2007a,b, 2008b] for more very encouraging results in concrete examples. In his 1795 seminal paper, Riche de Prony studies the parameter estimation of a finite sum of sinusoidal functions (see Riche de Prony [1795], Kahn et al. [1992], Osborne et al. [1995]). In this paper, we are interested in Prony’s problem for a two-terms sum of complex sinusoidal func978-3-902823-06-9/12/$20.00 © 2012 IFAC

tions, meaning that our aim is to estimate the parameters of the signal (see also Neves et al. [2007] for a quite related study) x(t) = α1 exp i(ω1 t + φ1 ) + α2 exp i(ω2 t + φ2 ) from the biased and noisy output measure y(t) = x(t) + β + $

(1)

where β is an unknown constant bias and $ is a noise 1 . A linear parametric estimation problem may often be formalized as finding a good approximation of some vector Θ on the basis of an observed signal that is a linear functional of the true signal depending on a set of parameters and a noise corrupting the observation. Here the signal z(t) = x(t) + β and Θ are linearly differentially algebraic. Indeed, we have the following differential equation: z¨(t) − i (ω1 + ω2 ) z(t) ˙ − ω1 ω2 (z(t) − β) = 0. Notice that the signal frequencies appear as roots of the characteristic equation of the above relation which is also behind Prony’s method and other techniques such as adaptive notch filters: take this equation as the numerator of some transfer function, then the filter zeros align with the signal frequencies. The further applied algebraic operator will take advantage of this remark. In the operational domain, we obtain 1

We use here the framework developed in Fliess [2006, 2008]. It is independent of any probabilistic modeling. In this point of view, the noise should be viewed as fast oscillations (see Fliess [2006]).

167

10.3182/20120711-3-BE-2027.00393

16th IFAC Symposium on System Identification Brussels, Belgium. July 11-13, 2012

s(s − iω1 )(s − iω2 )Z(s) = s(s − i(ω1 + ω2 ))z(0) + sz(0) ˙ − βω1 ω2 . (2) Among the unknown parameters, we wish to estimate Θest := {θ1 , θ2 , θ3 , θ4 }, but not the bias Θest = {θ5 }: θ1 = −ω1 ω2 ,

θ2 = −i(ω1 + ω2 ),

θ3 = −x(0) = β − z(0),

θ4 = −z(0) ˙ = −x(0), ˙

(4) (5)

Notice that the parameters α1 , ω1 , φ1 , α2 , ω2 , φ2 can be easily deduced from Θest , see Lemma 15 in the Appendix. Using the notation (3)-(5), the equation (2) reads as 2


s + θ2 s + θ1 s Z(s) + (s2 + θ2 s)θ3 + θ4 s +
s2 + θ2 s + θ1 θ5 = 0. (6) From (6) we would like to obtain a system of equations on Θest and independent of Θest . We consider the algebraic extensions CΘ := C(Θ), CΘest := C(Θest ) and CΘ := C(Θest ). Moreover, CΘ [s] (respectively CΘest [s], est CΘ [s]) denotes the polynomial ring in the variable s with est coefficients in CΘest (respectively in CΘ ). We obtain the est relation: R (s, Z(s), Θest , Θest ) := P (s) Z(s) + Q(s) + Q(s) = 0 (7) with P (s) = s T (s) and Q(s) = s2 θ3 + s(θ4 + θ2 θ3 ) ∈ CΘest [s] and Q(s) = T (s)θ5 ∈ CΘ

est

[s],

where we set T (s) = s2 + θ2 s + θ1 ∈ CΘest [s]. Now, we proceed in three steps: (1) Algebraic elimination of all terms in Θest : all differential operators that annihilate be generated by a  dQ  can 2 single operator in CΘest (s) ds , called a Q-minimal annihilator. These annihilators will be rewritten in a canonical form. (2) Obtaining a system of equations on Θest : we apply on R several differential operators annihilating Q. Using their canonical forms, choices will be made to obtain a system of equations with good numerical properties (once back in the time domain). (3) Resolution of the obtained system: we use the inverse Laplace transform Z 1   p p m+p −1

L

1 d Z(s) sm dsp

(−1) t = (m − 1)!

w

m−1,p

(τ )z(tτ )dτ

0

(8)

with wm,p (t) = (1−t)m tp , ∀ p, m ∈ N, m ≥ 1 to bring the equations back in the time domain. The integers m, p will be chosen as small as possible so that the resulting estimation is as least as possible sensitive to the noise. The first point emphasizes the central role played by Qminimal annihilators. They are defined in Section 2 where we describe the algebraic structure behind them and the algebraic elimination technique to eliminate Θest (subsection 2.1). Minimal annihilators can be rewritten in a 2

The polynomial ring in

d ds

2. AN ALGEBRAIC FRAMEWORK FOR OBTAINING ANNIHILATORS

(3)

θ5 = −β

3

canonical form, detailed in subsection 2.1. The parameter estimation is provided in section 3. Section 4 contains convincing numerical experiments that illustrate our techniques and are easily implementable. These experiments are compared to the well-known modified Prony’s method.

with coefficients in CΘest [s]

The algebraic framework is borrowed from Fliess et al. [2010, 2003, 2008], Fliess [2008], Mboup [2009] 3 . More details about the algebraic notions can be found in Dixmier et al. [1974] and McConnell et al. [2000]. Recall that we wish to annihilate
Q = s2 + θ2 s + θ1 θ5 ∈ CΘ

est

[s].

Since we have a polynomial in s, the natural idea is to d look for operator in ds . It is clear that operators of order 4 d d greater than 3 will annihilate Q, e.g. Π1 = (s ds −2)◦(s ds −
3 d d 1) ◦ s ds and Π2 = ds3 . One question would be whether these annihilators are the same. Another question would be whether there exists a lower order annihilator. Answers are provided  d  by the algebraic structure of the Weyl algebra CΘ (s) ds . 2.1 Weyl algebra We consider K a field of characteristic zero (here K = C). Definition 1. Let k ∈ N \ {0}. The Weyl algebra Ak (K) is the K-algebra generated by p1 , q1 , . . . , pk , qk satisfying the relations [pi , qj ] = δij , [pi , pj ] = [qi , qj ] = 0, ∀ 1 ≤ i, j ≤ k where [·, ·] is the commutator defined by [u, v] := uv − vu, ∀ u, v ∈ Ak (K). We will simply write Ak instead of Ak (K) when we do not need to make explicit the base field. A well-known fact is that Ak can be realized as the algebra of polynomial differential operators on K[s1 , . . . , sk ] with ∂ pi = and qi = si × · , ∀ 1 ≤ i ≤ k. (9) ∂si We can also write h i Ak = K[q1 , . . . , qk ][p1 , . . . , pk ] = K[s1 , . . . , sk ]

∂ ∂ ,..., ∂s1 ∂sk

(notice that we use the same notation for the variable si and for the operator multiplication by si here). There is a closely related algebra to Ak : it is defined as the set of differential operators on K[s1 , . . . , sk ] with coefficients in the rational functions field K(s1 , . . . , sk ). We denote it by Bk (K) and write i h Bk (K) := K(q1 , . . . , qk )[p1 , . . . , pk ] = K(s1 , . . . , sk )

∂ ∂ ,..., ∂s1 ∂sk

Sometimes we will simply write Bk instead of Bk (K). In the case k = 1 for instance, we have     d d A1 = hp, q | pq−qp = 1i = K[s] and B1 = K(s) ds ds 3

Similar tools were also used for numerical differentiation of noisy signal (see Mboup et al. [2009], Liu et al. [2011]) and change-point detection (see Fliess et al. [2010]).   4 The order of an operator Π ∈ C (s) d is its degree as a Θ ds polynomial in the variable

168

d . ds

16th IFAC Symposium on System Identification Brussels, Belgium. July 11-13, 2012

 Proposition 2. A basis for Ak is given by q I pJ | I, J ∈ Nk where q I := q1i1 . . . qkik and pJ := pj11 . . . pjkk . An element F ∈ Ak can be written in its canonical form, X F = λIJ q I pJ with λIJ ∈ K. I,J

Example 3. In this paper the following identities are useful: m n

n m

m    X m n

k!(−1)k pn−k q m−k k k k=1 n    X n m pn q m = q m pn + i!q m−i pn−i i i

q p =p q

+

Remark 9. Remark 8 still holds if we replace B by BΘest in Definition 7. That means that AnnBΘ (R) is generated by est a unique generator Πmin ∈ BΘest (up to multiplication by a polynomial in CΘest (s)) called a minimal R-annihilator w.r.t. BΘest . We have the following lemmas: Lemma 10. Let Qn = sn , n ∈ N. A minimal Qn annihilator is d Πn = s − n. ds (unique up to a multiplication by a polynomial in C(s)). Let us note that for m, n ∈ N, the operators Πm and Πn commute. Thus one can take advantage of the following lemma Lemma 11. Let P1 , P2 ∈ CΘ [s]. Let Πi be a Pi est annihilator (i = 1, 2) such that Π1 Π2 = Π2 Π1 . Then Π1 Π2 is a (µP1 + ηP2 )-annihilator for all µ, η ∈ CΘ .

k=1

Similarly, an element F ∈ Bk can be written as X F = gI (s)pI with gI (s) ∈ K(s1 , . . . , sk ). I

est

P The order of an element F ∈ Bk , F = I gI (s)pI is defined as ord(F ) := max{[I| | gI (s) 6= 0}. Notice that the same definition is valid for the Weyl algebra Ak since Ak ⊂ Bk . Proposition 4. Ak is a domain. Moreover, Ak is simple and Noetherian. However, Ak is neither a principal right domain, nor a principal left domain. Proposition 5. For the algebra Bk , one has: (1) Bk is a domain. Moreover, Bk is simple and Noetherian. (2) B1 admits a left division algorithm, that is, if F , G ∈ B1 , then there exists q, r ∈ B1 such that F = q G + r and ord(r) < ord(G). As a consequence, B1 is a principal left domain. d Since ds is a derivation operator we have: Proposition 6. (Derivation). For F ∈ C[s] and X(s) the Laplace transform of a signal x(t), we have (Leibniz rule): n   k X n d F dn−k X dn (F X) = dsn k dsk dsn−k k=0

2.2 Annihilator In the sequel, to distinguish whether an annihilator de- d pends on Θest or not, we work with B := B1 = C(s) ds d and BΘest := CΘest (s) ds . Definition 7. Let R ∈ CΘ [s]. Consider the left ideal: est

AnnB (R) = {F ∈ B | F (R) = 0} . A R-annihilator w.r.t. B is an element Π ∈ AnnB (R). Remark 8. Let us note that AnnB (R) is a left principal ideal. Thus it is generated by a single generator Πmin ∈ B called a minimal R-annihilator w.r.t. B. So AnnB (R) = B Πmin . Remark that AnnB (R) contains annihilators in finite inte gral form, i.e. operators with coefficients in C 1s .

Using this Lemma and considering Q in (7), we obtain a minimal
Q-annihilator w.r.t. B where Q = s2 + θ2 s + θ1 θ5 ∈ CΘ [s]: est       d d d Πmin = s − 2 ◦ s − 1 ◦ s . (10) ds ds ds The identities in Example 3 give: d3 . ds3 Lemma 12. Let R ∈ CΘest [s]. Then a minimal Rannihilator w.r.t BΘest is Πmin = s3

Πmin = R

d dR − . ds ds

In our example, by Lemma 12 we obtain that the Qannihilator w.r.t. BΘest is:
d Θest − (2s + θ2 ) . Πmin = s2 + θ2 s + θ1 ds 3. PARAMETER ESTIMATION Our estimation problem is equivalent to find a family d r of annihilators (Πi )i=1 in C(s) ds that applied on (7) provides a set of equations enabling the computation of Θest in the time domain. This family will be generated by a minimal Q-annihilator. Formula (8) justifies the use of finite-integral form annihilators. Moreover, these operators should be of minimal d degree in ds to minimize noise sensitivity. Lastly, the obtained system of equations should be well-balanced to provide good numerical estimation. Recall that the minimal Q-annihilator Πmin is a generator of the ideal AnnB (Q), so a general Q-annihilator will be of the form: ! ` X di Π= gi (s) i ◦ Πmin (11) ds i=0 where gi (s) ∈ C(s), ∀ i = 0, . . . , `.

169

16th IFAC Symposium on System Identification Brussels, Belgium. July 11-13, 2012

3.2 Estimation of θ3 and θ4

3.1 Estimation of θ1 and θ2 3

d Applying Πmin = s3 ds 3 , the minimal Q-annihilator w.r.t. B, on relation (7) gives

d3 Z(s) d2 Z(s) + P (s) 2 ds3 ds2 dZ(s) +P1 (s) + P0 (s)Z(s) ds where P0 (s) = 6s3 , P2 (s) = 9s5 + 6s4 θ2 + 3s3 θ1 , 6 5 4 P1 (s) = 18s4 + 6s3 θ2 ,
P3 (s) = s + s θ2 + s θ1 , and Πmin (Q(s)) = Πmin Q(s) = 0. That provides the following algebraic relation Πmin (P (s) Z(s)) = P3 (s)

P3 (s)

d2 Z(s) dZ(s) d3 Z(s) + P2 (s) + P1 (s) + P0 (s)Z(s) = 0. 3 ds ds2 ds

We obtain a single equation in θ1 and θ2 . To linearly identify these two parameters, we need two independent equations. However, the following result show that this is not possible in the operational domain (see the appendix for a proof): Theorem 13. There do not exist two Q-annihilators w.r.t B leading to two independent equations in θ1 and θ2 . Remark 14. Let us note that for a similar parameter identification problem of a single sinusoid, it is indeed possible to find two independent equations in the operational domain (see Ushirobira et al. [2011]). Therefore we will use such a construction in the time domain. For this, since Q-annihilators are of the form (11), to get two equations we select ` = 1. That leads to the following 4th -order annihilator written in the canonical form: d3 d4 Π = g0 (s) 3 + g1 (s) 4 , (12) ds ds where g0 (s), g1 (s) ∈ C(s). The choices of g0 (s) = 1, g1 (s) = 0 and then g0 (s) = 0, g1 (s) = 1, give two equations in the operational domain leading to the following system in the time domain:   1    1  − t I1 −I2  θ1 = − t I5  16 6

where R1 I1 =

0

I2 =

R1

I3 =

R1

I4 =

R1

I5 =

R1

I6 =

R1

0 0 0 0 0

t2 I3 −tI4

θ2

I6




6w2,3 (τ ) − 3w2,2 (τ ) z(tτ )dτ

Using annihilators generated by the minimal Q-annihilator  d w.r.t. B := C(s) ds given by (10), we could linearly identify the parameters θ1 and θ2 . These annihilators do not depend on the parameters to be found. Now, one can show that it is not possible to identify linearly the remaining parameters θ3 and θ4 , so we will use nonlinear equations in θ1 and θ2 . So, let us consider the Q-annihilator w.r.t. BΘest :
d Θest − (2s + θ2 ) . Πmin = s2 + θ2 s + θ1 ds Applying it on (7) gives:   dZ(s) Θest Πmin (P (s) Z(s)) = T (s)2 s + Z(s) ds Θ

est (Q(s)) = −θ s2 + 2θ θ s + θ (θ + θ θ ) Πmin 4 1 3 1 4 2 3 Θ

est (Q(s)) = 0 Πmin

where T (s) = s2 +θ2 s+θ1 ∈ C[s]. That gives the following algebraic relation:   dZ(s) 2 T (s) Z(s) + s + (θ1 − s2 )θ4 + θ1 (2s + θ2 )θ3 = 0. ds We obtain a single equation in θ3 and θ4 . As in the previous subsection, we use (11) with ` = 1, where Πmin is replaced Θest by Πmin . That leads to the following 4th -order annihilator written in the canonical form:   d2 Θest (13) Π = g0 (s)Πmin + g1 (s) T (s) 2 − 2 , ds where g1 (s), g2 (s) ∈ C(s). Selecting g0 (s) = 1, g1 (s) = 0 and then g0 (s) = 0, g1 (s) = 1, we obtain two equations:      B3 θ3 θ1 (2s + θ2 ) (θ1 − s2 ) , (14) =− B4 θ4 2θ1 −2s with   dZ(s) B3 = T (s)2 Z(s) + s (15) ds d2 Z(s) dZ(s) B4 = sT (s)2 + T (s)(6s2 + 4sθ2 + 2θ1 ) ds2 ds +T (s)(4s + 2θ2 )Z(s) (16) The above matrix can be easily written in a diagonal form and the expressions for θ3 and θ4 are thus obtained. Notice that they depend on θ1 and θ2 .




4. SIMULATIONS

5w1,3 (τ ) − 5w1,2 (τ ) + w1,1 (τ ) z(tτ )dτ




7w2,4 (τ ) − 4w2,3 (τ ) z(tτ )dτ




−7w1,4 (τ ) + 8w1,3 (τ ) − 2w1,2 (τ ) z(tτ )dτ




w3,0 (τ ) − 9w2,1 (τ ) + 9w1,2 (τ ) − w0,3 (τ ) z(tτ )dτ




−4w3,1 (τ ) + 18w2,2 (τ ) − 12w1,3 (τ ) + w0,4 (τ ) z(tτ )dτ

The expressions for θ1 and θ2 are thus obtained: 6 (−I4 I5 + I6 I2 ) θ1 = − 2 t I4 I1 + I3 I2 1 I1 I6 − I3 I5 θ2 = t I4 I1 + I2 I3

Figure 1 shows the simulation results for the estimation of the parameters θ1 and θ2 vs the estimation time. The modified Prony’s method (PM) is used as a reference. Each point is obtained by averaging the results over 200 trials. The first plot displays the real part of a sample realization of the noisy signal (1), with N = 1024 samples and SNR = 20 dB). The corresponding value of the constant bias is β = 7.9. The results in the second and third plots show that the estimation is unbiased and insensitive to the constant bias β for the presented method (solid line curve). For the same experiment, the results obtained with the modified Prony’s method display a bias in the estimates (dot-dashed line curve). Note that even by

170

16th IFAC Symposium on System Identification Brussels, Belgium. July 11-13, 2012

11 10 9 8 7 6 5 4 3

REFERENCES

A sample realisation of the real part of the noisy signal − SNR = 20dB

0

1

2

0

1

2

0

1

2

2.0

3

4

5

6

7

8

9

10

3

4

5

6

7

8

9

10

3

4

5

6

7

8

9

10

θ1 vs estimation time: mean values

S. Bittanti, M. Campi and S. Savaresi, Unbiased estimation of a sinusoid in colored noise via adapted notch filters, Automatica, vol. 33, pp. 209–215, 1997. S. Bittanti and S. Savaresi, On the parameterization and design of an extended Kalman filter frequency tracker, IEEE Trans. Automat. Control, vol. 45, pp. 1718–1724, 2000. P. Cartier and Y. Perrin, Integration over finite sets, Nonstandard Analysis in Practice, F. & M. Diener (Eds), Springer, Berlin, pp. 195–204, 1995. M. H. Cheng and J. L. Tsai , A new IIR adaptive notch filter, Signal Process, vol. 86, pp. 1648–1655, 2006. J. Dixmier, Alg`ebres enveloppantes, Gauthier-Villars, 1974. M. Fliess, Analyse non standard du bruit, C.R. Acad. Sci. Paris, ser. I, 342, pp. 797–802, 2006. M. Fliess, Critique du rapport signal `a bruit en communications num´eriques, ARIMA, vol. 9, pp. 419–429, 2008.

1.5 1.0 0.5 0.0 −0.5 3.5 3.0 2.5 2.0 1.5 1.0 0.5

θ2 vs estimation time: mean values

Fig. 1. Comparison with the modified Prony’s method: sample mean setting β = 0 with the Modified Prony’s method (dashed line curve), the results of the proposed method still show more robustness to the noise corruption. This is confirmed by the corresponding variances which are displayed below

1.0 0.9 0.8 0.7 0.6 0.5 0.4 0.3 0.2 0.1 0.0

Variance of θ1

http://hal.archives-ouvertes.fr/inria-00001059/en/

0

1

2

3

4

5

6

7

8

9

10

5

6

7

8

9

10

M. Fliess and H. Sira-Ram´ırez, An algebraic framework for linear identification. ESAIM Control Optim. Calc. Variat., vol. 9, pp. 151–168, 2003. M. Fliess, C. Join, M. Mboup and H. Sira-Ram´ırez, Compression diff´erentielle de transitoires bruit´es. C.R. Acad. Sci. Paris,, ser. I(339), pp. 821–826, 2004. M. Fliess and H. Sira-Ram´ırez, Closed-loop parametric identification for continuous-time linear systems via new algebraic techniques, in H. Garnier & L. Wang (Eds): Identification of Continuous-time Models from Sampled Data, Springer, pp. 362–391, 2008. M. Fliess, C. Join and H. Sira-Ram´ırez, Non-linear estimation is easy. Int. J. Modeling Identif. Control, vol. 4, pp. 12–27, 2008.

Variance of θ2

1.6 1.4 1.2 1.0 0.8 0.6 0.4 0.2 0.0

http://hal.archives-ouvertes.fr/inria-00311719/en/

M. Fliess, C. Join and M. Mboup, Algebraic changepoint detection, Applicable Algebra Engin. Communic. Comput., vol. 21, pp. 131–143, 2010. M. Fliess, M. Mboup, H. Mounier and H. Sira-Ram´ırez, Questioning some paradigms of signal processing via concrete examples, in H. Sira-Ram´ırez, G. Silva-Navarro (Eds.): Algebraic Methods in Flatness, Signal Processing and State Estimation, Editorial Lagares, M´exico, pp. 1– 21, 2003.

0

1

2

3

4

Fig. 2. Comparison with the modified Prony’s method: sample variance

5. CONCLUSION The algebraic method we used in this paper proved to be very efficient in this parameter estimation problem of a sum of two sinusoidal waveform signal. Indeed, the two triplets (amplitude, frequency, phase) could be easily identified using the minimal annihilators proposed here. In the case of a single sinusoidal wave, the results are also convincing (see Ushirobira et al. [2011]). The perspective of a extension of this algebraic method to a sum of several sinusoidal waveform signal is challenging, but hopefully very likely to be succeeded. A positive result in the three-sinusoid case has already been obtained.

http://hal.archives-ouvertes.fr/inria-00158855/en/

S. Haykin, Adaptive Filter Theory (2nd ed.), Englewood Cliffs, 1991. L. Hsu, R. Ortega and G. Damm, A globally convergent frequency estimator, IEEE Trans. Automat. Control, vol. 44, pp. 698–713, 1999. M. Kahn, M. Mackisack, M. Osborne and G. K. Smyth, On the consistency of Prony’s method and related algorithms, J. Comput. Graph. Statist., vol. 1, pp. 329– 349, 1992. T. H. Li and B. Kedem, Strong consistency of the contraction mapping method for frequency estimation, IEEE Trans. Inform. Theory, vol. 39, pp. 989–998, 1993. D. Y. Liu, O. Gibaru and W. Perruquetti, Error analysis for a class of numerical differentiator: application to state observation, 48th IEEE Conference on Decision and Control , China, 2009. D. Y. Liu, O. Gibaru and W. Perruquetti, Differentiation by integration with Jacobi polynomials, Jour-

171

16th IFAC Symposium on System Identification Brussels, Belgium. July 11-13, 2012

nal of Computational and Applied Mathematics, 2011, doi:10.1016/j.cam.2010.12.023 D. Y. Liu, O. Gibaru and W. Perruquetti, Error analysis of a class of derivative estimators for noisy signals, Numerical Algo., 2011. M. Mboup, Parameter estimation for signals described by differential equations, Applicable Analysis, vol. 88, pp. 29–52, 2009. M. Mboup, C. Join and M. Fliess, Numerical differentiation with annihilators in noisy environment, Numerical Algo., vol. 50, pp. 439–467, 2009. J. C. McConnell and J. C. Robson, Noncommutative Noetherian Rings, Amer. Math. Soc., 2000. M. Mojiri and A. R. Bakhsahi, An adaptive notch filter for frequency estimation of a periodic signal, IEEE Trans. Automat. Control, vol. 49, pp. 314–318, 2004. A. Neves, M. Mboup and M. Fliess, An algebraic receiver for full response CPM demodulation, VI Int. Telecom. Symp. (ITS2006), Fortaleza , CE, Brazil, 2006. http://hal.archives-ouvertes.fr/inria-00086115/en/

A. Neves, M. Miranda and M. Mboup, Algebraic parameter estimation of damped exponentials, EUSIPCO’07, Poznan, Poland, 2007.

Appendix A. FINDING THE ORIGINAL PARAMETERS Lemma 15. The parameters α1 , ω1 , φ1 , α2 , ω2 , φ2 can be obtained from Θest . Proof. From θ1 and θ2 one can easily deduce ω1 and ω2 since iω1 and iω2 are the roots of the polynomial ω 2 + θ2 ω + θ1 = 0. From θ3 and θ4 one obtains x(0) = −θ3 = α1 exp iφ1 + α2 exp iφ2 and x(0) ˙ = −θ4 = i(ω1 α1 exp iφ1 + ω2 α2 exp iφ2 ) and then easily deduce the remaining parameters (α1 , φ1 , α2 , φ2 ):      −θ3 1 1 α1 exp iφ1 = −θ4 iω1 iω2 α2 exp iφ2     1 −iθ4 − ω2 θ3 α1 exp iφ1 = iθ4 + ω1 θ3 α2 exp iφ2 ω2 − ω1 Taking absolute values and arguments and ordering the frequencies ω2 > ω1 : p

http://hal.archives-ouvertes.fr/inria-00601655/en/

Y. Xiao and Y. Tadokoro, LMS-based notch filter for the estimation of sinusoidal signals in noise, Signal Processing, vol. 46 (2), pp. 223–231, 1995.

(A.1)

ω2 − ω1 <(θ4 ) + ω2 =(θ3 ) tan(φ1 ) = − =(θ4 ) − ω2 <(θ3 )

http://hal.inria.fr/inria-00179732/en/

M. R. Osborne and G.K. Smyth, A modified Prony algorithm for exponential function fitting, SIAM Journal on Scientific Computing, vol. 16 (1), pp. 119-138, 1995. P. A. Regalia, Adaptive IIR Filtering in Signal Processing and Control, Marcel Dekker, New York, 1995. G. M. Riche de Prony, Essai exp´erimental et analytique : sur les lois de la dilatabilit´e de fluides ´elastiques et sur celles de la force expansive de la vapeur de l’eau et de la vapeur de l’alcool ` a diff´erentes temp´eratures, Journal de l’´ecole polytechnique, vol. 1, no. 22, pp. 2476, 1795. R. Roy and T. Kailath, ESPRIT-estimation of signal parameters via rotational invariance techniques, IEEE Trans. Signal Process, vol. 37, pp. 984–995, 1989. J. R. Trapero, H. Sira-Ram´ırez and V. F. Battle, An algebraic frequency estimator for a biased and noisy sinusoidal signal, Signal Processing, vol. 87, pp. 1188– 1201, 2007. J. R. Trapero, H. Sira-Ram´ırez and V. F. Battle, A fast on-line frequency estimator of lightly damped vibrations in flexible structures, J. Sound Vibration, vol. 307, pp. 365–378, 2007. J. R. Trapero, M. Mboup, E. Pereira-Gonzalez and V. B. Feliu, On-line frequency and damping estimation in a single-link flexible manipulator based on algebraic identification, MED, 2008. J. R. Trapero, H. Sira-Ram´ırez and V. F. Batlle, On the algebraic identification of the frequencies, amplitudes and phases of two sinusoidal signals from their noisy sum, Int. J. Control, vol. 81, pp. 507–518, 2008. R. Ushirobira, W. Perruquetti, M. Mboup M. and M. Fliess M., Estimation alg´ebrique des param`etres intrins`eques d’un signal sinuso¨ıdal biais´e en environnement bruit´e, in GRETSI, Bordeaux, 2011.

(=(θ4 ) − ω2 <(θ3 ))2 + (<(θ4 ) + ω2 =(θ3 ))2

α1 =

p α2 =

(A.2)

(−=(θ4 ) + ω1 <(θ3 ))2 + (<(θ4 ) + ω1 =(θ3 ))2

ω2 − ω1 <(θ4 ) + ω1 =(θ3 ) tan(φ2 ) = − −=(θ4 ) + ω1 <(θ3 )

(A.3) (A.4)

Appendix B. PROOF OF THEOREM 13 A Q-annihilator is of the form (11), that is ` X di+3 Π= gi (s) i+3 , ds i=0

(B.1)

where gi (s) ∈ CΘest (s), ∀ 1 ≤ i ≤ `. We apply such annihilator on relation R (7). Since Π(Q(s)) = Π(Q) = 0, we compute Π(P (s) Z(s)). Using proposition 6, one obtains:  4  X i+3 di+3 P (k) (s)Z (i+3−k) (s) (P (s) Z(s)) = k dsi+3 k=0

since P ∈ C[s] is of degree 3. Moreover: P (0) (s) = s3 + θ2 s2 + θ1 s, P (1) (s) = 3s2 + 2θ2 s + θ1 , P (2) (s) = 6s + 2θ2 , and P (2) (s) = 6. Thus using notations (9), Π given by (B.1) applied on relation R (7) reads in the Weyl algebra framework as: ` X gi (q) (Ai (p, q)θ1 + Bi (p, q)θ2 + Ci (p, q)) = 0 (B.2) i=0

where Ai (p, q) = (qp + ai ) pi+2 ,
Bi (p, q) = q 2 p2 + 2ai qp + bi pi+1 ,
Ci (p, q) = q 3 p3 + 3ai q 2 p2 + 3bi qp + ci pi , with ai = (i + 3), bi = ai (i + 2), ci = bi (i + 1), ∀0 ≤ i ≤ `. Denote by Li := [Ai (p, q)Bi (p, q)Ci (p, q)], then for i > j we have by induction on (j − i) and by example 3: Lj = p(j−i) Li and that completes the proof of theorem 13.

172

(B.3)