André Monin

Salut Gérard

Exact Arma Lattice Predictors From Autocorrelation Functions

Keywords: Linear filtering, linear predictor, reduced filters, lattice form

This paper derives an optimal linear -predictor of ARMA type in lattice form of arbitrarily fixed dimension for a process whose autocorrelation function is known. The algorithm preserves exact optimality at each step, as opposed to asymptotic convergence of more usual algorithms, at the expense of hereditary computation. Only the discrete time case is examined. It is shown how the unnormalized (respectively normalized) lattice form may be reduced to only 4n-2 parameters (respectively 2n+1) for a n-th order projection on the past. The normalization algorithm for the forward and backward residuals uses only scalar square root computations. Some examples are given which show the accuracy of this technique compared to those using the classical ARMA form for the predictor.

Introduction.

Objectives.

The paper concerns the problem of computing the exact minimum variance ARMA predictor of n-th order, as a time varying lattice, from the autocorrelation function. Approximate procedures have been proposed in the past [START_REF] Benveniste | AR and ARMA identification algorithm of Levinson type : an innovations approach[END_REF], [START_REF] Benveniste | Algorithme en treillis pour les modèles stationnaires[END_REF], [START_REF] Morf | A classification of algorithms for ARMA models and ladder realizations[END_REF], [START_REF] Morf | Recursive square root ladder estimation algorithm[END_REF], [START_REF] Lee | Recursive ladder algorithms for ARMA modelling[END_REF] to achieve the same goal, namely by trivially imbedding the past innovations into a 2-channel standard lattice problem. However, it has been shown by [START_REF] Benveniste | AR and ARMA identification algorithm of Levinson type : an innovations approach[END_REF], as well as noticed in [START_REF] Lee | Recursive ladder algorithms for ARMA modelling[END_REF], that this is equivalent to the extended leastsquares approximation [START_REF] Soderstrom | A theorical analysis of recursive identification methods[END_REF], therefore with the same pitfalls (computed transient coefficients are not the true projection coefficients, and do not guarantee convergence to the optimal, in the general case. For the same reason, the predictor is always biased in the short-term transient.). More recent papers have used stochastic gradient approximations of higher order [START_REF] Glentis | Fast adaptive algorithm for multichannel filtering and system identification[END_REF], [START_REF] Shynk | Adaptive IRR Filtering[END_REF], [START_REF] Michaut | Méthodes adaptatives pour le signal[END_REF], [START_REF] Benveniste | Algorithmes adaptatifs et approximations stochastiques[END_REF]. Except the fact that they use data records while we use autocorrelations, the problem is indeed similar to ours. However, the aim is to asymptotically determine the stationary optimal predictor parameters, not to obtain the optimal time-varying parameters. The techniques developed are essentially based on extensions of the well-known stochastic gradient search. They are recursive but do not look for the optimum at each step. Moreover, some convergence problem may occur when, for example, the dimension of the unknown system is greater that the model [START_REF] Asttröm | Uniqueness of the maximum likelihood estimates of the parameters of an ARMA model[END_REF]. It is obviously the case when the system is not finite dimensional, which is the case in the low-pass filter treated an illustration in 5.

To derive the exact ARMA lattice predictor, thus avoiding such pitfalls, is a more difficult task since the feedback parameters of the lattice lead to a non-linear implicit dependence or, equivalently, a growing memory for the computation of the time-varying optimal predictor pa-

L 2
rameters. This is equivalent to the exact stochastic gradient (i.e. an infinite order approximation). This difficulty was clearly identified in [START_REF] Benveniste | AR and ARMA identification algorithm of Levinson type : an innovations approach[END_REF], [START_REF] Benveniste | Algorithme en treillis pour les modèles stationnaires[END_REF] but was never fully investigated, to the author's knowledge. We examine here this problem for a known autocorrelation function.

Our results, in this respect, are too-fold:

-A general algorithm gives the exact (hereditary) computation of the optimal n-th order ARMA predictor, in time-varying lattice form with zero initial state, for an arbitrary autocorrelation function. The linear dependence in the vector process (infinite Hankel rank) is being exploited to obtain the minimum number of parameters.

-In the finite Hankel rank case, the computations are shown to be finitely recursive, and completely solve the so-called realization problem under lattice form.

Contents.

Let

, be a discrete-time scalar sequence of (mean square random variables) with zero mean. Assume that the only data available is the autocorrelation function of , not necessarily with a finite Hankel rank. One can derive, as in [START_REF] Monin | The optimal linear predictor in reduced dimension[END_REF], the reduced optimal n-th order non-stationary linear predictor of (in the minimum variance sense), under its observable canonical state-space form or, equivalently, its input-output ARMA form:

(1.1)
The optimal parameters are obtained by minimizing the prediction variance error ( ) with no stationary assumption but explicit state-dimension constraint on the predictor, as shown in [START_REF] Monin | The optimal linear predictor in reduced dimension[END_REF]. Otherwise, the state variable has to be reconstructed for the whole time horizon, at each step, leading to an infinite-dimensional state-space as in ARMA identification [START_REF] Benveniste | AR and ARMA identification algorithm of Levinson type : an innovations approach[END_REF].

In such a predictor, computing the values of requires the inversion of a dimensional system at each step. Despite the possibility of recursive inversion, as in [START_REF] Benveniste | AR and ARMA identification algorithm of Levinson type : an innovations approach[END_REF], the canonical state-space form is generally not tractable for higher dimensions. It is in fact well known that numerical computations are very sensitive with respect to the coefficients of the characteristic polynomial.

Recall that the optimization of leads to the orthogonality equations as:

(1.2)
Clearly, that equations (1.2) represent the projection of on the random variables and with as inner product. We write this projection:

(1.3)
The purpose of this paper is to exhibit an orthogonal basis of the linear space spanned by the random variables in order to simplify the projection formula and to improve its numerical accuracy. This orthogonal basis induces a lattice type filter but without the independent modular structure, which is well known to be impossible in the ARMA case. Note that this problem has been studied in [START_REF] Benveniste | AR and ARMA identification algorithm of Levinson type : an innovations approach[END_REF], [START_REF] Benveniste | Algorithme en treillis pour les modèles stationnaires[END_REF] in the context of the ARMA identification, but with the extended least squares approximate. y y , ( )

y t t N ∈ , { } L p 2 Ω ( )
y t z t y ˆt t 1 - = y t z t a t i z t i - i 1 = n ∑ b t i y t i - + = J t E y t z t - ( ) 2 [ ] = a t i b t i , { } i 1 = n 2n 2n × J t t ∀ 0 i ∀ , ≥ 1…n 1 -E y t z t - ( )y t i - [ ] 0 = E y t z t - ( )z t i - [ ] 0 = , = y t y t i - { } i 1 = n z t i - { } i 1 = n x y ( ) E xy [ ] = z t P y t y t 1 -… y t n -z t 1 -… z t n - , , , , , [ ] = y t i -z t i - , { } i 1 = n
The main contribution here is as follows:

-We first exhibit the optimal n-th order predictor (reduced dimension) under the unnormalized lattice form using 4n-2 parameters ( -dimensional state-space) and a complete algorithm is provided which uses hereditary computation on the autocorrelation function . -We then introduce the normalized form using 2n+1 parameters with its algorithm attached. The problem of extracting square roots of the normalization matrices is solved using only scalar square root computation because of its particular form.

-We also show that such a projection technique (i.e. onto a linear functional of the past with internal memory constraint) yields an associated innovation whose "whiteness" has the same restricted meaning: orthogonality to an n-dimensional sub-space of the past. It yields a linear model of with "n-th order white noise", whose approximation properties are illustrated through applications.

The applications concern the most useful case, and deal with autocorrelation function of infinite Hankel rank (low-pass and fractional order filters), one of which is well known for its sensitivity problems. By converting the role of and , it must be noticed that the algorithm gives the exact n-th order best -approximate with minimum phase, an open problem in the literature on optimum systems approximations [START_REF] Olivi | Approximation en norme des systèmes dynamiques[END_REF].

The orthogonalisation formulas.

Let with inner and outer products written as usual and .

The linear projection of a vector X of on the linear product space spanned by vectors of is obtained by solving the following dimensional linear system as:

(2.1)
where are matrices. The projection takes then the following form:

(2.2)
The orthogonal projection of a vector X of on a scalar random variable of leads to:

, i.e. . ( 2

.3)

In the same way, one may compute the orthogonal projection of a scalar of on a vector of . One finds:

(2.4)

We now recall some properties that will be used further to reduce the number of parameters of the predictor in lattice form.

2n y t y τ y y y L p 2 X Y L P 2 ( ) n ∈ , E X T Y [ ] E X Y T [ ] L P 2 ( ) n Y 1 … Y p , , L P 2 ( ) n n p × i ∀ 1…p E X Y i T [ ] , K j E Y j Y i T [ ] j 1 = p ∑ = = K i n n × P X Y 1 … Y p , , [ ] K i Y i i 1 = p ∑ = L P 2 ( ) n z L P 2 E X P X z [ ] - ( ) z [ ] 0 = P X z [ ] E Xz [ ] E z 2 [ ] ---------------z = z L P 2 Y L P 2 ( ) n P z Y [ ] E zY T [ ] E Y Y T [ ] (
) . Then the projection of a vector X of on the space spanned by is: .❏

1 - Y = Proposition 1. Let X

Unormalized lattice form.

Recall that our purpose is to compute the n-th order linear predictor as a projection of the observation signal on the n-th order memory of the past as in (1.3). This will be achieved in two steps. First we exhibit an orthogonal basis of the past using the Levinson formulas. We then reduce the dimension of the involved system using the degeneracy property of the projection space (clearly, is only n-dimensional). We then show how the computation of the predictor's parameters is achieved knowing the autocorrelation function of .

Forward and backward residuals formulas.

Define

, where denotes the innovation in the sense of (1.1)-(1.3). It is clear from (2.4) that the projection of on the 2n scalars can be equivalently computed as the projection of on the n vectors . That is:

(3.1)
The aim of this paragraph is to recall an orthogonalisation of the space spanned by [START_REF] Benveniste | AR and ARMA identification algorithm of Levinson type : an innovations approach[END_REF], using the Graham-Shmidt procedure. Let:

(3.2)
be the backward residuals of the family and

L p 2 ( ) 2 X x 1 x 2 = Y y 1 y 2 = E X y 1 [ ] 0 = E y 1 y 2 [ ] 0 = P X Y [ ] P X y 2 [ ] = Z 1 …Z p , L P 2 ( ) n i ∀ j , 1…n = i j E Z i Z j T [ ] ⇒ ≠ 0 = L P 2 ( ) n Z i { } i 1 = n P X Z 1 … Z p , , [ ] P X Z 1 [ ] … P X Z p [ ] + + = z t y t y t i -z t i - , { } i 1 = n y t i -z t i - , { } i 1 = n y t Z t y ˜t z t = y ˜t y t z t - y t y t i -z t i - , { } i 1 = n y t Z t i - { } i 1 = n z t P y t Z t 1 -… Z t n - , , [ ] = Z t i - { } i 1 = n X t 0 Z t = X t 1 Z t 1 - P Z t 1 -Z t [ ] - = … X t n Z t n - P Z t n -Z t n - 1 + … Z t , , [ ] - =          Z t … Z t n - 1 + , , { } (3.3)
be the forward ones.

We have then the following properties:

Property 1.
is orthogonal, that is .

Property 2. is orthogonal.

Property 3.

. From the definition of the projection, the linear predictor becomes:

(3.9)

The implicit aspect of this formula can be eliminated by using the orthogonality of the innovation with the past . It leads to an explicit expression of the predictor, gi-ven the statistics of the observation signal:

Z t 0 Z t = Z t 1 Z t P Z t Z t 1 - [ ] - = … Z t n Z t P Z t Z t n -… Z t 1 - , , [ ] - =          X t 0 … X t n 1 - , , { } i ∀ j E X t i X t j ( ) T [ ] , ≠ 0 = Z t n - 1 + 0 … Z t n 1 - , , { } i ∀ 1…n 1 - Z t … Z t i - , , { } sp , X t o … X t i , , { } sp = = i ∀ 1…n 1 - Z t … Z t i - , , { } sp , Z t i - o … Z t i , , { } sp = = X t 0 Z t = X t 1 X t 1 - 0 P X t 1 - 0 Z t 0 [ ] - = … X t n X t 1 - n 1 - P X t 1 - n 1 -Z t n 1 - [ ] - =          Z t 0 Z t = Z t 1 Z t 0 P Z t 0 X t 1 - o [ ] - = … Z t n Z t n 1 - P Z t n 1 -X t 1 - n 1 - [ ] - =          Z t n y ˜t 1 0 Z t n y ˜t P y ˜t Z t 1 -… Z t n - , , [ ] - = P y ˜t Z t 1 -… Z t n - , , [ ] 0 = 0 1 Z t n z t P z t Z t 1 -… Z t n - , , [ ] - = P z t Z t 1 -… Z t n - , , [ ] P y t Z t 1 -… Z t n - , , [ ] z t = = z t P z t X t 1 - 0 [ ] … P z t X t 1 - n 1 - [ ] + + = y ˜t X t 1 - 0 … X t 1 - n 1 - , , (3.10)
Remark: As was noticed by Benveniste in [START_REF] Benveniste | Algorithme en treillis pour les modèles stationnaires[END_REF], the vector is degenerated because it only depends on the scalar sequence . If we examine the relations (3.4), it appears at first that the lattice predictor needs 2n matrices as parameters, that is 8n parameters. In fact, because of the degeneracy of vector , the system can be reduced to only 4n-2 parameters as will be shown in the next section. (3.12)

But according to Property 4, one has:

(3.13)
Applying then the Proposition 2, we have:

.

(3.14)

The orthogonality of the family implies that is orthogonal to , so that . As a consequence, it is clear that

(3.15)
On the other hand:

.

(3.16)

By linear combination of equations (1.2) using the equality (1.1), one can easily derive that . But is also orthogonal to , according to (1.2). It leads to:

z t P y t X t 1 - 0 [ ] … P y t X t 1 - n 1 - [ ] + + = Z t y t 2 2 × Z t i ∀ 0…n 1 - = P X t 1 - i Z t i [ ] P X t 1 - i z t i [ ] = z t i Z t i Z t i Z t P Z t Z t 1 -… Z t i - , , [ ] - = 1 0 Z t i y ˜t P y ˜t Z t 1 -… Z t i - , , [ ] - y ˜t = = E X t 1 - i Z t i ( ) T [ ] 1 0 E y ˜tZ t 1 - i - [ ] E y ˜tP Z t 1 - i -Z t i -… Z t 1 - , , [ ] [ ] - = P Z t 1 - i -Z t i -… Z t 1 - , , [ ] P Z t 1 - i -Z t i - 0 … Z t 1 - i 1 - , , [ ] = P Z t 1 - i -Z t i -… Z t 1 - , , [ ] P Z t 1 - i -Z t i - 0 [ ] … P Z t 1 - i -Z t 1 - i 1 - [ ] + + = Z t i Z t 1 - i 1 -… Z t i - 0 , , , { } Z t i Z t 1 - i 1 -… Z t i - 0 , , { } E y ˜tZ t i - 0 [ ] … E y ˜tZ t 1 - i 1 - [ ] 0 = = = i ∀ 0…n 1 - E X t 1 - i Z t i ( ) T [ ] 1 0 , 0 = = 1 0 E Z t i Z t i ( ) T [ ] 0 1 E y ˜t z t P z t Z t 1 -… Z t i - , , [ ] - ( ) [ ] = E y ˜t z t [ ] 0 = y ˜t Z t i - { } i 1 = n (3.17)
Applying then Proposition 1, the result follows. ■ A useful consequence of this proposition is that the backward and forward residuals of system (3.4) need only the second component of vectors . The residuals may therefore be computed using only 4n-2 parameters by the following system.

and (3.18)

Recall that according to (3.9), . It leads to the following lattice structure of computation:

and (3.19)

where the output of the predictor is computed as:

(3.20)
According to (4) and ( 5), the parameters and are determined as follows:

(3.21)

The figure 1 illustrates the structure obtained.

i ∀ 0…n 1 - 1 0 E Z t i Z t i ( ) T [ ] 0 1 , 0 = = Z t i { } i 0 = n 1 - z t 0 z t = z t 1 z t 0 P y t X t 1 - 0 [ ] - = … z t n z t n 1 - P y t X t 1 - n 1 - [ ] - =          X t 0 Z t = X t 1 X t 1 - 0 P X t 1 - 0 z t 0 [ ] - = … X t n 1 - X t 1 - n 2 - P X t 1 - n 2 -z t n 2 - [ ] - =          z t n 0 = z t 0 z t = z t 1 z t 0 k t 0 ( ) T X t 1 - o - = … z t n 1 - z t n 2 - k t n 2 - ( ) T X t 1 - n 2 - - = 0 z t n 1 - k t n 1 - ( ) T X t 1 - n 1 - - =            X t 0 Z t y ˜t z t = = X t 1 X t 1 - 0 q t 0 z t 0 - = … X t n 1 - X t 1 - n 2 - q t n 2 -z t n 2 - - =            z t k t i ( ) T X t 1 - i i 0 = n 1 - ∑ = k t i { } i 0 = n 1 - q t i { } i 0 = n 2 - i ∀ 0…n 1 - k t i E X t 1 - i X t 1 - i ( ) T [ ] ( ) 1 - E X t 1 - i y t [ ] = q t i E X t 1 - i y t [ ] E z t i ( ) 2 [ ] --------------------------- =        , =

Off-line computation of parameters.

We now show how the above parameters are recursively computed, in the general case, (i.e. infinite Hankel rank for the process ) by a growing dimension system (as in [START_REF] Monin | The optimal linear predictor in reduced dimension[END_REF]). In such a case, the algorithm needs no particular dimensional caution, but only a stationary convergence test. It should be noted that, in case of singular initialization, the first steps call for a particular formula. This will be developed further.

When the process has a finite recursion model of order (hankel rank), but , it is shown how the corresponding (reduced) predictor's coefficients are themselves generated by finite recursion formulas.

The general algorithm.

First, let us examine the quantities to be computed to derive the lattice parameters defined in (3.21). For simplicity, we write the expectation of X. Define .

The first system of (3.19) permits us to evaluate as:

(3.22)
As a consequence, according to the Property 1, . One then simplify the formulas (3.21) as follows:

(3.23)
It appears then that there are two kinds of quantities to be computed: the finitely realizable one, and the hereditary one as . Note also and are finitely realizable.

• Finitely recursive computations.

Recall that, at each step , and are assumed to be known, so are .

+ - - + - y t z t X 1 z 1 q 0 k 0 k 0 + + - -X 2 z 2 q 1 k 1 k 1 z -1 + z n 2 - -X n 1 - k n 2 - z -1 + q n 2 - k n 1 - z -1 + + + + + + z 0 X 0 z -1 Figure 1 -Unnormalized Lattice Structure. y y m m n > X E X [ ] = M t i X t i X t i ( ) T = z t i i ∀ 0…n 1 -z t i , k t j ( ) T X t 1 - j j i = n 1 - ∑ = = X t 1 - i z t i M t 1 - i k t i = i ∀ 0…n 1 -k t i , M t 1 - i ( ) 1 - X t 1 - i y t = = i ∀ 0…n 1 -q t i , X t 1 - i y t z t i ( ) 2 ⁄ = =      X t 1 - i y t M t i z t i ( ) 2 t M t i X t 1 - i y t k t i q t i , { } i 0 = n 1 -
The next value of is computed from system (3.21) as follows:

(3.24)

But according to (3.22) and because of the orthogonality of the family , one has:

(3.25)

The system above may then be transform in:

(3.26)
In a similarly way, using system (3.19), the matrices may be recursively computed as:

(3.27)
Using (3.25) a second time, one has:

(3.28)
Note that is well defined at step because is assumed to be given data in the problem and is computed from (3.24).

Hereditary computations (for infinite Hankel rank of ).

We use such an expression in the usual sense of non-finitely recursive computation (as a function, here, of the time horizon of the autocorrelation function ).

Although are recursive quantities, the expectation is not so. In fact, its computation needs all the parameters from the beginning of the algorithm. System (3.18) allows us to derive the following formula:

z t i ( ) 2 z t n 1 - ( ) 2 k t n 1 - ( ) T X t 1 - n 1 -X t 1 - n 1 - ( ) T k t n 1 - = z t i ( ) 2 z t i 1 + ( ) 2 2 k t i ( ) T X t 1 - i z t i k t i ( ) T X t 1 - i X t 1 - i ( ) T k t i - + i ∀ , 0…n 2 - = =      X t 1 - 0 … X t 1 - n 1 - , , { } X t 1 - i z t i X t 1 - i X t 1 - i ( ) T k t i = z t n 1 - ( ) 2 k t n 1 - ( ) T M t 1 - n 1 -k t n 1 - = z t i ( ) 2 z t i 1 + ( ) 2 k t i ( ) T M t 1 - i k t i + i ∀ , 0…n 2 - = =      M t i i ∀ 1 X t 1 - i X t 1 - i ( ) T , ≥ X t 1 - i 1 -X t 1 - i 1 - ( ) T q t i X t 1 - i 1 - ( ) T z t i 1 - - X t 1 - i 1 -z t i 1 -q t i 1 - ( ) T - = q t i 1 -z t i 1 - ( ) 2 q t i 1 - ( ) T + M t i M t 1 - i 1 - q t i 1 -k t i 1 - ( ) T M t 1 - i 1 - - i ∀ 1 ≥ , = M t 0 y t ( ) 2 z t 0 ( ) 2 - 0 0 z t 0 ( ) 2 =          M t 0 t y t 2 z t 0 ( ) 2 y y t y τ X t i { } i 0 = n 1 - X t 1 - i y t k τ i q τ i τ t ≤ , , { } i 0 = n 1 - (3.29)
It appears then that a second hereditary quantity is to be computed, that is . A similar computation leads to:

(3.30)
where the autocorrelation function of appears as the input of the hereditary recursion (3.29)-(3.30). The figure 2 presents the sequence of computation.

Additional formulas for singular initialization.

Without a priori information on the past, the initial estimator at is the a priori expectation of the stochastic process to be estimated, so that . For a zero mean process, this value is zero.

At the first step( ), the only relevant data is since no information on past values of the predictor is available. As a consequence, the projection of on the past is and only the first components of and are determined, since is singular.

At the second step ( ), one may then project on . This yields the matrix non-singular but not . Therefore, and are well defined as is regular. At step , and are well defined but the first components of and are de-

0 ∀ τ t 1 - ≤ ≤ X τ 0 y t y τ y t z τ 0 y t - z τ 0 y t = X τ i y t X τ 1 - i 1 -y t q τ i 1 -z τ i 1 -y t - i ∀ , 1…n 1 - = =          , z τ i 1 -y t 0 ∀ τ t 1 - ≤ ≤ z τ n 1 -y t k τ n 1 - ( ) T X τ 1 - n 1 -y t = z τ i y t z τ i 1 + y t k τ i ( ) T X τ 1 - i y t + i ∀ , 0…n 2 - = =      , y t Initialization Computation of k t o … k t n 1 - , , Computation of z t o ( ) 2 … z t n 1 - ( ) 2 , , (20)
Computation of q t o … q t n 2 - , , (17)
Computation of M t 1 - o … M t 1 - n 1 - , , (21)
Computation of X t 1 - ), all parameters are well determined.❏

o y t … X t 1 - n 1 - y t , , (22,23)
t 0 = (17) t t 1 + → Figure 2 -Unnormalized Lattice Algorithm. t 0 = z 0 y 0 = t 1 = y 0 y 1 z 1 P y 1 y 0 [ ] = k 1 0 q 1 0 M 0 0 t 2 = y 2 y 1 z 1 , { } M 1 0 M 1 1 k 2 0 q 2 0 M 1 0 t 3 = k 3 0 q 3 0 k 3 1 q 3 1 termined.

Completely finite recursive computation for finite Hankel rank of y.

Suppose that the state-space model of is a finite-dimensional one. For example:

(3.31)
where , and are independent white noises of variances and respectively.

The computations of and are obviously still finitely recursive and are computed by the same formulas that in 3.3.1. The only difference lies in the computation of . It is achieved through:

(3.32)
Computation of is finitely recursive through the matrices as follows: since only depends on , one has:

(3.33) i n ≤ z 2i 2 - P y 2i 2 -X 2i 3 - o … X 2i 3 - i 2 - , , [ ] = k 2i 2 - j q 2i 2 - j , { } j 0 = i 2 - M 2i 2 - j { } j 0 = i 2 - z 2i 1 - P y 2i 1 -X 2i 2 - o … X 2i 2 - i 2 - 1 0 X 2i 2 - i 1 - , , , [ ] = k 2i 1 - j q 2i 1 - j , { } j 0 = i 2 - 1 0 k 2i 1 - i 1 - 1 0 q 2i 1 - i 1 - ,       ∪ M 2i 1 - j { } j 0 = i 2 - M 2i 1 - i 1 - i n = t 2n = k t q t , y x t 1 + F x t Gw t 1 + + = y t H x t v t + =    x t R m ∈ m n ≥ , w t v t Q R M t i z t i ( ) 2 y t 2 y t 2 H x t x t T H T R + = x t x t ( ) T F x t 1 -x t 1 - T F T GQG T + =      X t 1 - i y t R 2 m × X t 1 - i x t T       i 0 = n 1 - X t 1 - i { } i 0 = n 1 - w τ v τ , { } τ t 1 - ≤ X t 1 - i x t T X t 1 - i x t 1 - T F T i ∀ , 0…n

Normalized lattice form.

Recall that the variances of the forward residuals are decreasing with the orthogonalization process. As a consequence, the values of the parameters grow and their computation degenerates. Therefore, we now replace these residuals by normalized one. Normalization of forward residuals is quite easy. One just considers the variables which are colinear to the previous one and have a unit variance. These quantities being scalar, there is a unique way to achieve it, that is:

(4.1)
Recall that the backward residuals are vectors of . So, the normalization procedure is not unique. If one looks for a matrix of such that:

(4.2)
are normalized ( ), the solution is undetermined. In fact, any product of a solution with an orthogonal matrix will succeed. However, in such a special situation, there is a natural way to make this normalization.

Backward residuals normalization.

Define .

It is clear that is then normalized. Let us examine the recursive equation of for . According to system (3.19), one has: But if and are normalized, according to the definition of , one has:

X t 1 - 0 x t 1 - T H x t 1 -x t 1 - T z t 1 - 0 x t 1 - T - z t 1 - 0 x t 1 - T = X t 1 - i x t 1 - T X t 2 - i 1 -x t 1 - T q t 1 - i 1 -z t 1 - i 1 -x t 1 - T - i ∀ , 1…n 1 - = =          z t 1 - i x t 1 - T k t 1 - j ( ) T X t 2 - j x t 1 - T j i = n 1 - ∑ i ∀ , 0…n 1 - = = X t 1 - i x t T       k t i q t i , { } i 0 = n 1 - q t i { } i 0 = n 2 - i ∀ 0 … n 1 -u t i , , , z t i z t i ( ) 2 ⁄ = = X t i { } i 0 = n 1 - L p 2 ( ) 2 Q t i { } i 0 = n 1 - R 2 i ∀ 0 … n 1 -V t i , , , Q t i X t i = = V t i V t i ( ) T I = Q t 0 M t 0 ( ) 1 2 ⁄ - y t 2 z t 2 - ( )
1 2 ⁄ - 0 0 z t 2 ( ) 1 2 ⁄ - = = V t 0 Q t 0 X t 0 = V t i i 0 >
(4.8)
Substituting above the value of , the result follows.■

Forward residuals normalization.

One has to compute the variances of each residual. System (3.19) leads to the expression of the normalized forward residuals as: (4.9) But according to (3.21) and the definition of , one has . Equation (4.9) takes then the following form:

(4.10) V t i Q t i X t 1 - i q t i 1 -u t i 1 -z t i 1 - - ( ) = z z 2 ( ) 1 2 ⁄ = z q t i 1 -z t i 1 - X t 1 - i 1 -y t z t i 1 - ⁄ = Q t 1 - i 1 - ( ) 1 - V t 1 - i 1 -y t ( ) z t 1 - i 1 - ⁄ = p t i V t 1 - i u t i V t 1 - i y t z t i ⁄ = = V t i Q t i Q t 1 - i 1 - ( ) 1 - V t 1 - i 1 - p t i 1 -u t i 1 - - ( ) = i ∀ 1 … n 1 - , , = Q t i I p t i 1 - p t i 1 - ( ) T - ( ) 1 2 ⁄ - Q t 1 - i 1 - = V t i Q t i X t i = { } i 0 = n 1 - i 0 = V t i V t i V t i ( ) T Q t i Q t 1 - i 1 - ( ) 1 - V t 1 - i 1 -V t 1 - i 1 - ( ) T p t i 1 -u t i 1 - ( ) 2 p t i 1 - ( ) T +   = V t 1 - i 1 -u t i 1 - p t i 1 - ( ) T - p t i 1 -u t i 1 -V t i 1 - ( ) T -   Q t 1 - i 1 - ( ) T - Q t i ( ) T V t 1 - i 1 - u t i 1 - p t i V t i V t i ( ) T Q t i Q t 1 - i 1 - ( ) 1 - I p t i 1 - p t i 1 - ( ) T - ( )Q t 1 - i 1 - ( ) T - Q t i ( ) T = Q t i u t i u t i 1 -z t i 1 - k t i 1 - ( ) T Q t 1 - i 1 - ( ) 1 - V t 1 - i 1 - - ( )z t i ⁄ = p t i k t i Q t 1 - i ( ) T p t i z t i = u t i z t i 1 - z t i ---------------u t i 1 - p t i 1 - ( ) T V t 1 - i 1 - - ( ) = Equation (3.22) leads to
. As a consequence, using the fact that , one may derive: Proposition 6. The normalized lattice form is therefore:

with . ❑
The figure 3 shows the lattice structure generated. The coefficients and denote the normalization factors.

Remarks:

-Note that, in this case, only 2n+1 parameters are used to derive the lattice form predictor. That is and . It results from the fact that the forward and backward residuals are both computed from the same set of parameters.

-Square roots of matrices may be computed here from scalar square roots, using the following algebraic identity: , which derives from a straightforward calculation. ------------------------------------------

z t i 2 k t j ( ) T M t 1 - j k t j j i = n 1 - ∑ = M t i Q t i ( ) 1 - Q t i ( ) T - = z t i 1 - z t i --------------- 1 p t i 1 - ( ) T p t i 1 - - ( ) 1 2 ⁄ - = i ∀ 1 … n 1 - u t i 1 p t i 1 - ( ) T p t i - ( ) 1 2 ⁄ - u t i 1 - p t i 1 - ( ) T V t 1 - i 1 - - ( ) = V t i I p t i 1 - p t i 1 - ( ) T - ( ) 1 2 ⁄ - V t 1 - i 1 - p t i 1 -u t i 1 - - ( ) =      , , , = u t o z t z t ⁄ = V t o y t 2 z t 2 - ( )
1 2 ⁄ - 0 0 z t 1 - y ˜t z t =        q i Q i + - - + - y t z t V 1 u 1 p 0 p 0 p 0 + + - -X 2 p 1 p 1 p 1 + u n 2 - -V n 1 - p n 2 - + p n 2 - p n 1 - + + + + + + Q 0 z 1 - u 0 z 1 - Q 1 -z 1 - 1 z 0 ⁄ q 0 q 1 V 0 q n 2 - Q n 2 -z 1 - q 0 Figure 3 -Normalized Lattice Structure. p t i { } i 0 = n 1 - z t 0 I pp T - ( ) 1 2 ⁄ - I pp T 1 p T p - ( )
1 2 ⁄ - 1 - p T p
-       + =

Computation of parameters.

The general algorithm.

It is clear that: . (

According to definition (4.5) and the first remark above, the computation of parameters requires only the quantities . As it was noticed in the unnormalized case, these quantities are also hereditary ones and call for . In fact, suppose that and are known, according to (4.11), their computation is achieved by the following algorithm, with τ going from 1 to :

• First, assuming that is well defined, one computes:

• Then, one computes:

(4.13)

Additional formulas for singular initialization.

As in 3.3.2, the first steps of computation of the parameters lead to a singular problem if one chooses to initialize the predictor in the absence of a priori information on the past. For a zero mean process, one may consider that at , one has: Note that after these first singular steps, stands normalized. The figure 4 illustrates the algorithm obtained:

V t 1 - i y t       i 0 = n 1 - V τ i y t u τ i y t τ t 1 - ≤ , ,       i 0 = n 1 - p τ i τ t 1 - ≤ , { } i 0 = n 1 - z τ 0 τ t 1 - ≤ , { } t 1 - V τ 1 - i y t       i 0 = n 1 - u τ n 1 -y t p τ n 1 - ( ) T V τ 1 - n 1 -y t = u τ i y t p τ i ( ) T V τ 1 - i y t 1 p τ i ( ) T p τ i - u τ i 1 + y t + i ∀ , 0…n 2 - = =      V τ 1 - o y t y t y τ u τ 0 y t z τ 0 - ( ) y τ y τ z τ 0 2 -     ⁄ u τ 0 y t = V τ i y t I p τ i 1 - p τ i 1 - ( ) T - ( ) 1 2 ⁄ - V τ 1 - i 1 -y t p τ i 1 -u τ i 1 -y t - ( ) i ∀ , 1…n 1 - = =          2n t 0 = (4.
u 0 i 0 i ∀ , 0…n 1 - = = V 0 i 0 i ∀ , 1…n 1 - = = V 0 0 y 0 y 0 ⁄ 0 =          t 1 = V 0 0 u 1 0 p 1 0 ( ) T V 0 0 = p 1 0 V o u 0 p 1 0 1 = p 1 0 1 0 = 1 ∀ i n 1 - ≤ ≤ V 1 i , 0 = t 2 = p 2 0 i ∀ 1…n 1 -p 2 i , 0 = = p 2 0 1 = t 3 = p 3 0 p 2 0 1 = I p 2 0 p 2 0 ( ) T - p 3 1 1 = p 3 1 1 0 = i n ≤ 2i 1 - p 2i 1 - j { } j 0 = i 2 - p 2i 1 - i 1 - 1 0 = j ∀ i…n 1 -p 2i 1 - j , 0 = = 2i p 2 j { } j 0 = i 1 - p 2i i 1 - 1 = j ∀ i 1…n 1 - + p 2i j , 0 = = p n 1 -
Obviously, the finite Hankel rank case leads to a finitely computation of parameters with the normalized lattice form, as in 3.3.3 for the unnormalized lattice form.

Examples of applications.

As is well known the optimal linear predictor of a finite Hankel rank signal process may be viewed as a particular ARMA model (with zero inside the unit circle), the innovation process being a white noise input since the model and its whitening filter are inverse of one another. Otherwise, the linear predictor in reduced dimension makes the innovation process a "pseudowhite noise" in the sense that its orthogonalization property is limited to delays. Its variance is defined here by the asymptotic value (in the stationary case) of . In this sense, the ARMA predictor yields the best n-th order approximate model of the signal in the metric. This will be exploited to illustrate the accuracy of the method presented here because one can compare easily the approximation of the model with the model itself.

The algorithm being particularly designed to the prediction of a signal which can only be represented by an infinite dimensional system, we deal here with the two following cases:

-The low-pass filter.

-A fractional system.

Test on finite dimensional linear systems.

An obvious preliminary step before any application is to test the algorithm in a situation analog to Kalman filtering (n-dimensional linear realization of the autocorrelation). The compared impulse responses gave identical values up to the tenth significant number.

The figures 5 and 6 illustrate the convergence accuracy of the algorithm for a 4-dimensional linear system. The module and phase of the model thus identified are plotted for the 20 first steps of computation. Consider a signal resulting of the sum of a random frequency uniformly distributed in the interval and a white noise of variance . The autocorrelation of this signal is not realizable by a finite dimension system. Its value is obtained by the inverse Fourier transform of the function , i.e.:

t t 1 + → (32,33) V t 1 - 0 y t … V t 1 - n 1 -y t , , z t 0 … z t n 1 - , , p t 0 … p t n 1 - , , t 1 =
(5.1)

where is the Kronecker discrete function. Taking successively different dimensions of the predictor (n = 2,4 and 12 here), we obtained the optimal linear reduced predictor with its gain and phase characteristics of the figures 7 and 8 (with , and , ratio signal-noise).

t 0 f , [ ] R 1 f -f , [ ] y t y t τ + 4πfτ ( ) sin 4πfτ -------------------------Rδ 0 τ ( ) + = δ 0 τ ( ) n f 0.1 sampling-frequence × = R

Fractional system.

It can be shown that the use of fractional systems in control is of great interest with regard to the robustness of the control involved [START_REF] Oustaloup | Systèmes asservis linéaires d'ordre fractionnaire[END_REF].

Consider an irrational impulse response of the form . It is clear that this system can't be realized by a finite system. The autocorrelation of the signal is the inverse z-transform of this spectrum. A series development of the square root leads to the following autocorrelation: with .

The figure 17 illustrates the approximation modelling yielded with (near the continuous time situation). The power spectrum is plotted in the "loglog" coordinates to show the singular attenuation of by octave (or by decade) which characterizes this fractional system. ------------------------

Conclusion.

Let us recall as a conclusion that our problem deals with the computation of the exact time-varying ARMA predictor of order for any square-integrable signal , in the basis of its correlation function . It is well known that such a general problem has no finite dimensional solution, in term of a dynamical system with as the input. It is however an essential objective to design the optimal predictor in a recursive fashion which:

• preserves exact optimality at each step (minimum variance of predictor error)

• uses a minimum amount of new computations when the number of points of the autocorrelation function increases.

Our paper derives an algorithm which fulfills these goals, at the expense of the hereditary computation of the coefficients, in terms of the autocorrelation function. The hereditary aspect of the algorithm makes it essentially different from other algorithms. It is crucial throughout the paper and accounts for:

• exact optimality, for each , as a time-varying predictor of fixed order.

• convergence behavior, as a (unique) minimum-phase synthesis of approximations. The usefulness of the algorithm derived is that it yields the exact solution of best -approximation in infinite-dimensional situations such as ideal filters, distributed parameters, etc. The lattice form involved by the orthogonalization process is of interest concerning numerical calculus and allows to deal with greater dimensions than with classical ARMA forms.

A final point concerns the extension of these ideas when applied to raw data. This is a quite different problem, since one has to replace theoretical minimum variance in by effective least square residuals which comes a maximum likelihood formulation (with a necessarily gaussian assumption for the problem to make sense). Contrarily to the kind of separation principle often invoked, the optimal result does not consist in replacing the autocorrelation terms by the experimental average of corresponding correlation products, since this is suboptimal in the adaptive transient. This is the object of a forthcoming paper which extends to real-time applications the approach introduced here.

4 )

 4 According to the last equation of system(3.3), the first component of is equal to while the second one is equal to zero. In fact, a straightforward calculation leads to: .

3. 2 .

 2 Reduction of the number of parameters. Recall that, according to (3.3), one has: Let us examine the first component of these vectors. One has . (3.11) because of the orthogonality relations (1.2). Let us then check the hypothesis of Proposition 1:

  recursion requires only the storage at each step of the finite dimensional matrix . The parameters are obtained with (3.21).

  Recursively, the proposition is true for . Let us compute the covariance matrix of . According to (4.3), one has:

Figure 4

 4 Figure 4 -Normalized Lattice Algorithm.

Figure 5 -Figure 6 -

 56 Figure 5 -Module of a 4-Dimensional System Identified.

  y

0. 1 =Figure 7 -Figure 8 -FrequencyFigure 12 -Figure 13 -FrequencyFigure 14 -

 178121314 Figure 8 -Phase of the Optimal Low-pass Predictor in Reduced Dimension

Figure 15 -Figure 16 -

 1516 Figure 15 -Evolution of the module of Optimal Filter.

Figure 17 -

 17 Figure 17 -Power Spectrum of the Optimal Approximation of a Fractional System.

  And so on... Continuing this procedure, one finds that:

	Proposition 4. For all	, after 2i-2 iterations, one may just define:
	and		are well determined because the matrices	are
	regular. After 2i-1 iterations, one may define:
	a n d			a re w e l l d e t e r m i n e d s i n c e
		are regular but not	.
	From	(i.e.

14)

	At the first step (	), one can only project the observation on the first component of
	. As a consequence, according to Proposition 6, one has	and the second
	component of	is zero. But	and	being normalized, one has	and so
	. Therefore,			.
	At the second step (	),	is defined and	. Remark that at
	this step,	.
	At third step (	), only	is well defined but	implies that the matrix
	is singular. As a consequence, because	, one has	.
	And so on... Continuing this procedure, one find that:
	Proposition 7. For all	, computation of the predictor's parameters using definition (27)
	proceeds as follows:
	-after	iterations,		are well determined,	and
			.
	-a f t e r	i t e ra t i o n s ,		a re w e l l d e t e r m i n e d ,	, a n d
			.
	-after 2n iterations, all parameters are well determined.q


André Monin
email: monin@laas.laas.fr

Salut Gérard

ARMA Lattice Identification: A New Hereditary Algorithm

This paper derives an identification solution of the optimal linear predictor of ARMA type, as a time varying lattice of arbitrarily fixed dimension, for a process whose output signal only is known. The projection technique introduced here leads to an hereditary algorithm which is the adaptive extension to raw data of previous results of the authors on lattice realization from given autocorrelation functions ([1]). It produces a minimum phase linear model of the signal whose n-th order "whiteness" of the associated innovation has the following restricted meaning: orthogonality to an n-dimensional subspace memory of the past in a suitable Hilbert sequence space. The metric of that sequence space leads to a least-squares Identification algorithm which possesses a "certainty equivalence principle" with respect to the corresponding Realization algorithm (i.e., sample correlation products replace true correlation terms). Due to the detailed state-space time-varying computations, this is possible here while avoiding the well known "side errors" from missing correlation products which usually occur in a blunt replacement of the output autocorrelation by averaged sample products. Application examples show the superiority of the hereditary algorithm over classical recursive and non-recursive algorithms, in terms of accuracy, adaptativity and order reduction capabilities.

Introduction.

The ARMA stochastic identification problem has been extensively investigated for many years and numerous algorithms have been proposed. Because of its numerical robustness, linear identification under lattice form is of special interest. Moreover, it is often desirable to obtain o recursive algorithms on account of their ability to cope with adaptive filtering situations. Recall that the aim of ARMA identification is to fit a linear stochastic realization to an experimental set of measurements of the output (possibly of the input/output ) of some unknown system. For synthesis results on this subject, the reader is referred to our closely related paper [START_REF] Monin | Exact ARMA Lattice Predictor From Autocorrelation Functions[END_REF], on optimal rational approximations, which is the "realization" counterpart of the present "identification" algorithm.

The common approach consists in deriving a one-step-ahead predictor of the output and then evaluating the prediction performance through the experimental mean square error, i. e., in the stationary case, , in order to optimize the predictor parameters. This kind of formulation is useful to avoid the gaussian hypothesis under which the searched solution coincides with the maximum likelihood predictor.

Many different techniques have been proposed to achieve this goal and principally differ according to the class of predictors chosen as well as approximations made. One may first dis-

L 2 y τ τ , 1…t = { } u τ y τ τ , , 1…t = { } L P 2 y ˆτ J t 1 t ⁄ ( ) y τ y ˆτ - ( ) 2 τ 1 = t ∑ =
tinguish them according to the existence or not of feedback in the dynamic of the filter and then bring out two main classes of systems:

Transversal filters or Finite Impulse Response (FIR) filters.

The one-step-ahead predictor is represented as a moving average (MA) over the past output and takes the following canonical form:

(1.1)
It is the most popular class of predictors because the computation of the optimal parameters can be straightfully achieved as a linear quadratic optimum in a one-step recursive form.

Noting

and substituting to , it leads to an autoregressive (AR) model of the output. Many algorithms have been developed according to this formulation, particularly under various lattice forms ([2]- [START_REF] Carlson | Least Squares ARMA Modeling of Linear Time-Varying Systems: Lattice Filter Structures and Fast RLS Algorithm[END_REF] for instance). Unfortunately, it is well known that if the observation signals is corrupted by noise , which is most often the case in practice, such models are biased [START_REF] Shynk | Adaptive IRR Filtering[END_REF].

Rational filters or Infinite Impulse Response (IIR) filters.

The one-step-ahead predictor is represented as an ARMA over the past output and takes the following canonical forms:

(1.2)
This type of model is well known to be the most efficient linear form, with respect to the number of parameters, as it may represent IIR systems with only a finite set of parameters. Substituting to leads to an ARMA model:

(1.3)
whose general predictor is again of the type (1.2). This is why the general ARMA model is closed under additive measurement noises. Computation of the optimal parameters is however a more difficult task since the criterium is no longer quadratic with respect to the parameters. Moreover, the non-linearity involved is not explicitly tractable.

Finitely Recursive algorithms:

Many algorithms have been developed using essentially two kinds of approximations:

• The Extended Least Square (ELS) algorithm ignores the dependance of and with respect to the parameters when performing the minimization of the criterium [START_REF] Benveniste | AR and ARMA Algorithms of Levinson Type: An Innovations Approach[END_REF] [START_REF] Benveniste | Algorithmes Adaptatifs Et Approximations Stochastiques[END_REF]. The minimum is then achieved via classical linear regression algorithms or gradient type algorithms (Recursive Extended Least Square algorithm [START_REF] Ljung | System Identification: Theory for the User[END_REF] for

y ˆt a i y t i - i 1 = n ∑ = a i { } y ˜t y t y ˆt - = y t y ˆt y ˜t + y t v t y ˆt a i y ˆt i - i 1 = n ∑ b j y ˜t j - j 1 = q ∑ + = y t y ˜t y ˆt + y t a i y t i - i 1 = n ∑ a i y ˜t i - i 1 = n ∑ - b j y ˜t j - j 1 = q ∑ y + + = y ˆt i - y ˜t i - a i b j , { } J t instance).
• The innovation approach or prediction error methods leads to a prefiltered state vector which takes into account the implicit dependance neglected above. The optimal parameters are generally computed using local search methods such as stochastic gradient or Gauss-Newton algorithms [11] [12] [14] [START_REF] Michaut | Méthodes Adaptatives pour le Signal[END_REF].

(Recursive Maximum Likelihood algorithm [START_REF] Ljung | System Identification: Theory for the User[END_REF] for instance).

Hereditary algorithms:

In none of these approaches transient behavior is optimized. We are concerned in this paper with a time varying approach of IIR identification preserving a transient convergence similar to that of [START_REF] Monin | Exact ARMA Lattice Predictor From Autocorrelation Functions[END_REF] for known auto correlations, which has proved to be useful in terms of convergence rate as well as final accuracy. The price to be paid by such algorithms is their hereditary nature (growing memory) in terms of time , due to updating of the overall least-square criterium as data accumulate. We shall deal in the sequel with:

• Global hereditary criterium (stationary case):

(

• Locally limited memory criterium (adaptive case):

(1.5)

We develop the stationary case and will just underline the difference in the second case. We deal here with ARMA problems where the AR and MA orders are chosen equal for convenience. (Extension to the general case is easily obtained using representations such as in [START_REF] Benveniste | AR and ARMA Algorithms of Levinson Type: An Innovations Approach[END_REF]).

Problem setting: Suppose that the predictor which minimizes the criterium is available. Its parameters have been obtained from the knowledge of the data sequence . We will write it . When, at time , the new data is available, the n-th order linear predictor of is computed as the projection: ,

In this time varying situation, the criterium is quadratic with respect to the parameters since it only depends on the past values of the previously computed parameters, that is . On the other hand, it will lead to an hereditary computation. Note that although the representations (1.2) and (1.6) are different, they are asymptotically equiva-

y t 1 -… y t n - y ˜t 1 -… y ˜t n - , , , , { } t E t y • y ˆ• - ( ) 2 [ ] 1 t --- y τ y ˆτ - ( ) 2 τ 1 = t ∑ = E T t y • y ˆ• - ( ) 2 [ ] 1 T --- y τ y ˆτ - ( ) 2 τ t T - 1 + = t ∑ = E t 1 - y • y ˆ• - ( ) 2 [ ] a i t 1 -b i t 1 - , { } y 1 … y t 1 - , , { } y ˆτ t 1 - t y t y y ˆτ t a i t y ˆτ i - t i - b i t y ˜τ i - t i - + i 1 = n ∑ = τ ∀ 1…t 1 + = E t y • y ˆ• t - ( ) 2 [ ] a i t b i t , { } a i τ b i τ τ t < , , { }
lent in the steady state. Let us define the projection in that context.

The sequence space: Let be the Hilbert space spanned by sequences of length of random vector realizations of (with finite variance). The quadratic norm of a sequence of is defined as:

(1.7)
where is the euclidean norm of real vectors.

Orthogonal projections:

The projection at time of an element of over another element is defined as the linear transformation which minimizes the norm . The projector is the optimal matrix of appropriate dimension. Minimizing the criterium with respect to leads to the following well known projection equation:

(1.8)
When the process is not degenerated (non-singular autocorrelation), the orthogonal projection takes the following form:

(1.9)

Notations:

We will write for the above projection and for the extended inner product of sequences of . Although this notation may seems a bit heavy, it is necessary because of the non stationarity of the extended inner product. It means precisely the projection of the trajectory over the trajectory , the optimization being made over the interval . With this definition, the optimal predictor of (1.6) can be rewritten as:

(1.10) Note that the above formulation is not complete as and are not defined for . For instance, the projection of cannot be computed from the formula (1.10) as

are not yet defined. Remember, however, that an unknown variable is equivalent to a degenerate random variable (infinite variance) so that the linear projection of any other variable over it is zero. In other words, the projection over an unknown variable is the same as the projection over zero. With this adding, the projector is then well defined and leads to consider the usual windowed case, that is . A useful consequence of this assumption is:

(1.11)
This is a direct consequence of the definition of the projector (1.9).

The Levinson formulas for ARMA process.

Our aim now is to recursively compute an orthogonal basis, according to the inner prod-

H t X • t L P 2 X • X τ τ 1…t = , { } = H t X • t 2 1 t --- X τ 2 τ 1 = t ∑ = X X T X = t X • H t Y • X ˆτ t M t Y τ = X • X ˆt• - t 2 M t M t E t X • M t Y • - ( ) Y • T [ ] 0 = Y X ˆτ t E t X • Y • T [ ] E t Y • Y • T [ ] ( ) 1 - Y τ τ ∀ , 1…t = = X ˆτ t P τ t X • Y • [ ] = X Y 〈 | 〉 t E t X • Y • ( ) T [ ] = H t τ X τ → τ Y τ → 1 t , [ ] y ˆτ t P τ t y • y ˜• 1 - t 1 - y ˆ• 1 - t 1 - , … y ˜• n - t n - y ˆ• n - t n - , , , [ ] τ ∀ , 1…t = = y ˆτ y ˜τ τ 0 ≤ y ˆ1 1 y ˆ0 1 y 0 , y τ y ˆτ 0 τ 0 ≤ ∀ , = = Y 0 0 = ( ) P τ t X • 1 -Y • 1 - [ ] P τ 1 - t 1 - X • Y • [ ] = ⇒ uct
, of the space over which the projection is made. This is commonly achieved using the classical Gram-Schmidt procedure. Note that strict application of this procedure, which run in a 2n-dimensional space , would lead to a rather involved computation [START_REF] Carlson | Least Squares ARMA Modeling of Linear Time-Varying Systems: Lattice Filter Structures and Fast RLS Algorithm[END_REF]. It is convenient in our case to work with a 2-channel structure in order to preserve the recursive structure of the set with respect to time [START_REF] Benveniste | AR and ARMA Algorithms of Levinson Type: An Innovations Approach[END_REF] [START_REF] Monin | Exact ARMA Lattice Predictor From Autocorrelation Functions[END_REF]. In fact, we orthogonalize this set by 2-dimensional block corresponding to one step backward. The derived basis is therefore not really orthogonal (which would lead to a diagonal covariance matrix) but only block orthogonal (a 2-dimensional block diagonal covariance matrix). The space is here decomposed as . It is defined below.

The unormalized form.

Definition 1 : Let us define first the 2-dimensional vectors . We define the projection space as:

(2.1)

The projection of at time can then be rewritten as:

(2.2)
Our goal is first to find a set of 2-dimensional vectors of which spans the same space, and such that . At time , it is equivalent, according to (1.11), to make sure that .

Definition 2 : Backward residuals obey the following defining relations:

(2.3)
With this definition, it is clear that:

(2.4)
and the one step ahead predictor can be computed as the sum of 2-dimensional projections as:

(2.5)

These are the 2-dimensional backward residuals of the Levinson formulas. Note that the above direct sum is related to the extended inner product of , It would be convenient to make the residuals computation recursive in the same way as in the classical case. To this aim, let us set the following:

〈 | 〉 t 1 - y ˜• 1 - t 1 - y ˆ• 1 - t 1 - , … y ˜• n - t n - y ˆ• n - t n - , , , { } sp H t 1 - 2n H t 2n H t 2 ( ) n Z τ t y ˜τ t y ˆτ t , [ ] T τ , ≡ 1…t 1 - = H t 1 - 2n H t 1 - 2n y ˜• 1 - t 1 - y ˆ• 1 - t 1 - , … y ˜• n - t n - y ˆ• n - t n - , , , { } sp Z • 1 - t 1 -… Z • n - t n - , , { } sp ≡ ≡ y t y ˆτ t P τ t y • Z • 1 - t 1 -… Z • n - t n - , , [ ] = X • 1 - i t 1 - , { } i 0…n 1 - = H t 1 - 2 i j ≠ ∀ X • 1 - i t 1 - , X • 1 - j t 1 - , 〈 | 〉 t , 0 = t X • i t , X • j t , 〈 | 〉 t 0 = X τ 0 t , Z τ t = X τ i t , Z τ i - t i - P τ t Z • i - t i -Z • i - 1 + t i - 1 + … Z • t , , [ ] - i ∀ , 1…n 1 - = =      y ˜• 1 - t 1 - y ˆ• 1 - t 1 - , … y ˜• n - t n - y ˆ• n - t n - , , , { } sp X • 1 - 0 t 1 - , { } sp … X • 1 - n 1 -t 1 - , { } sp ⊕ ⊕ = n y ˆτ t P τ t y • X • 1 - 0 t 1 - , [ ] … P τ t y • X • 1 - n 1 -t 1 - , [ ] + + = H t 1 - 2
Definition 3 : Forward residuals obey the following defining relations:

(2.6)
With this definition, for all , we have the following decomposition:

(2.7)

the orthogonality being related one more time to the inner product . With this decomposition, the backward residuals can be written as follows:

(2.8)

Using then that, at the previous time , we had:

(2.9)

we may rewrite (2.3), using formula (1.11), under the following recursive form:

(2.10)

Carrying this expression into (2.3) and using the orthogonality of within the space , we find that:

(2.11)

Note that the set is degenerated since, according to (2.3), the first component of is equal to , the first component of . Therefore, because of the orthogonality , the second component is also orthogonal to . As a consequence, the projection of over being zero, one may simplify the projection term as:

(2.12)

with the following definition:

Definition 4 : second component of .
The systems (2.3) and (2.6) are then finally equivalent to:

Z τ 0 t , Z τ t = Z τ i t , Z τ t P τ t Z • t Z • 1 - t 1 -…Z • i - t i - , [ ] - i ∀ , 1…n 1 - = =      i 1…n 1 - = Z • t … Z • i - t i - , , { } sp Z • i t , { } sp … Z • i - 0 t i - , { } sp ⊕ ⊕ = 〈 | 〉 t X τ i t , Z τ i - t i - P τ t Z • i - t i -Z • 1 - i 2 -t 1 - , …Z • i - 1 + 0 t i - 1 + , , [ ] - P τ t Z • i - t i -Z • i 1 -t , [ ] - = t 1 - X τ 1 - i 1 -t 1 - , Z τ i - t i - P τ 1 - t 1 - Z • i - 1 + t i - Z • i 2 -t 1 - , …Z • i - 2 + 0 t i - 1 + , , [ ] - = X τ i t , X τ i 1 -t 1 - , P τ t Z • i - t i -Z • i 1 -t , [ ] - = Z τ i 1 -t , Z • 1 - i 2 -t 1 - , { } sp … Z • i - 1 + 0 t i - 1 + , { } sp ⊕ ⊕ X τ i t , X τ 1 - i 1 -t 1 - , P τ t X • 1 - i 1 -t 1 - , Z • i 1 -t , [ ] - = Z • i t , … Z • i - 0 t i - , , , { } sp Z τ i t , y ˜τ t Z τ t E t y ˜t t y ˆτ t [ ] 0 = y ˜τ t X τ i t 1 - , y ˜τ t P τ t X • 1 - i t 1 - , Z • i t , [ ] P τ t X • 1 - i t 1 - , z • i t , [ ] = z • i t , ≡ Z • i t , (2.13)
(2.14)

According to the orthogonal projection formula (1.9), if we define the projection coefficients as: , the systems (2.13) and (2.14) take the following form:

(2.15) (2.

16)

Remark: Note that the projection coefficients may seem to be implicitly defined as they depend on the forward residuals which are not available at the step . But, according to (2.16) and to the orthogonality properties (2.4), it is clear that:

(2.17)

Moreover, as is orthogonal to , we have:

(2.18)
and the predictor parameters are well defined.

The normalized form.

The above orthogonalization implies that the magnitude of the forward residuals decreases quickly with . As a consequence, it is interesting, from a numerical point of view, to normalize the forward and backward residuals in view to deal with quantities of similar magnitude. Concerning the forward residuals , there is a unique way to make the normalization, that is to define the variables as:

(2.19)

X τ 0 t , Z τ t = X τ i t , X τ 1 - i 1 -t 1 - , P τ t X • 1 - i 1 -t 1 - , z • i 1 -t , [ ] - i ∀ , 1…n 1 - = =      z τ 0 t , y ˆτ t = z τ i t , z τ i 1 -t , P τ t z • i 1 -t , X • 1 - i 1 -t 1 - , [ ] - i ∀ , 1…n 1 - = =      K t i E t X • 1 - i t 1 - , z • i t , [ ] E t z • i t , ( ) 2 [ ] ⁄ = Q t i E t z • i t , X • 1 - i t 1 - , ( ) T [ ]E t X • i t 1 - , X • i t 1 - , ( ) T [ ] ( ) 1 - = X τ 0 t , y ˜τ t y ˆτ t , [ ] T = X τ i t , X τ 1 - i 1 -t 1 - , K t i 1 -z τ i 1 -t , - i ∀ , 1…n 1 - = =      z τ 0 t , y ˆτ t = z τ i t , z τ i 1 -t , Q t i 1 -X τ i 1 -t , - i ∀ , 1…n 1 - = =      z • i t , { } i 0…n 1 - = t E t z • i t , X • 1 - i t 1 - , [ ] E t z • 0 t , X • 1 - i t 1 - , [ ] = E t y ˆτ t X • 1 - i t 1 - , [ ] = y ˜τ t y τ y ˆτ t - = X τ 1 - i t 1 - , { } i 0 = n 1 - E t z • i t , X • 1 - i t 1 - , [ ] E t y • X • 1 - i t 1 - , [ ] = z σ i t , i z τ i t , n u τ i t , u τ i t , z τ i t , z • i t , t ⁄ = with
as a definition.

On the other hand, the normalization of the backward residuals is not uniquely defined as they are 2-dimensional vectors, but only up to any non singular orthogonal matrix.There is however a natural way to achieve this normalization:

Consider the 2-dimensional normalized vectors of such that with , the identity matrix. At time , it leads to make sure that which is equivalent to make sure that . We proceed to compute the normalization matrices recursively on the order . Let us develop the normalized backward vectors. Equations (2.13) and (2.19) lead to:

(2.

20)

Let us then compute the autocorrelation matrix of :

(2.21)

At the previous step, one has: (2.24)

As a consequence, it is natural to choose the normalization matrices such that:

(2.25)

A similar procedure, applied to the normalization coefficients of the forward residuals leads to the relation:

(2.26)
Finally, the orthonormalized form writes:

z • i t , t E t z • i t , ( ) 2 [ ] = n H t 2 V • i t , { } i 0…n 1 - = V τ i t , R t i X τ i t , i ∀ , 0…n 1 - = = E t V • 1 - i t 1 - , V • 1 - i t 1 - , ( ) T [ ] I = t E t 1 + V • 1 - i t , V • 1 - i t , ( ) T [ ] I = t t 1 + -----------E t V • i t , V • i t , ( ) T [ ] I = i 0…n = V τ i t , R t i R t 1 - i 1 - ( ) 1 - V τ 1 - i 1 -t 1 - , E t V • 1 - i 1 -t 1 - , u • i 1 -t , [ ] u τ i 1 -t , - ( ) = V τ i t , E t V • i t , V • i t , ( ) T [ ] R t i R t 1 - i 1 - ( ) 1 - E t V • 1 - i 1 -t 1 - , V • 1 - i 1 -t 1 - , (
)

T [ ] ( = E - t V • 1 - i 1 -t 1 - , u • i 1 -t , [ ] E t u • i 1 -t , ( ) 2 [ ]E t V • 1 - i 1 -t 1 - , ( ) T u • i 1 -t , [ ] ) R t i R t 1 - i 1 - ( ) 1 - ( ) × T E t V • 1 - i 1 -t 1 - , V • 1 - i 1 -t 1 - , (
) 3. The algorithm.

T [ ] I = S t i E t V • 1 - i t 1 - , u • i t , [ ] = I R t i R t 1 - i 1 - ( )
1 - t t 1 + -----------     I S t i 1 -S t i 1 - ( ) T - ( )R t i R t 1 - i 1 - ( ) 1 - ( ) T = R t i R t i R t 1 - i 1 - ( )
1 - t t 1 + -----------I S t i 1 -S t i 1 - ( ) T - ( )     1 2 ⁄ - = z • i t , t z • i 1 -t , t -------------------- 1 S t i 1 - ( ) T S t i 1 - - ( ) 1

The recursive steps.

The aim of the algorithm is to compute the parameters of the predictor, that is . As it has been set in our problem statement, they depend on the whole trajectories and and are not actually finitely recursive. Let us examine more precisely what this means. According to (2.23), (2.18) and (2.19), the lattice parameters are defined as follows:

(3.1)
Recall that the forward residuals obey the following recursion:

(3.2)
So, we have:

(3.3)
and it appears that the parameters are directly linked to the quantities . In fact, one has to remember the whole past trajectory in order to compute the new trajectory which allows to extract the expression computed as the sum over the interval and finally the parameters with

V τ 0 t , E t 1 + y • 1 - y ˆ• 1 - t - ( ) 2 [ ] 0 0 E t 1 + y ˆ• 1 - t ( ) 2 [ ] 1 2 ⁄ - y τ y ˆτ t - ( )
y ˆτ t = V τ i t , t t 1 + -----------     I S t i 1 -S t i 1 - ( ) T - ( ) 1 2 ⁄ - V τ 1 - i 1 -t 1 - , S t i 1 -u τ i 1 -t , - ( )i , 1…n 1 - = = u τ 0 t , y ˆτ t z • 0 t , t ⁄ = u τ i t , 1 S t i 1 - ( ) T S t i 1 - - ( ) 1 2 ⁄ - u τ i 1 -t , ( S t i 1 - ( ) T - V τ 1 - i 1 -t 1 - , ) i ∀ 1…n 1 - = = N i I S i 1 -S i 1 - ( ) T - (
)

1 2 ⁄ - ≡ n i 1 S i 1 - ( ) T S i 1 - - ( )
1 2 ⁄ ≡ S t i y τ τ , 1…t = { } y ˆτ t 1 -τ , 1…t 1 - = { } S t i E t V • 1 - i t 1 - , y • [ ] z • i t , t ⁄ = z τ i t , z τ i 1 t , + P τ t y • X • 1 - i t 1 - , [ ] + = z τ i 1 t , + E t y • V • 1 - i t 1 - , [ ]V τ 1 - i t 1 - , + = z • i t , t E t y • V • 1 - k t 1 - , [ ] 2 k i = n 1 - ∑ i ∀ , 0…n 1 - = = E t V • 1 - i t 1 - , y • [ ] V τ i t 1 - , τ , 1…t 1 - = { } i 0…n 1 - = V τ i t , τ , 1…t = { } i 0…n 1 - = E t V • 1 - i t 1 - , y • [ ] 1 t , [ ] (3
.1) and (3.3). It should be noted that square roots of the matrices can be computed using the algebraic following identity:

(3.4)

The algorithm is described by Figure 2.

The initial steps.

As was noticed in 1., the data prior to time are assumed to be completely unknown, so that we cannot make any projection over them. As a consequence, the coefficients related to these variables are undetermined (or zero) and the algorithm established above doesn't run. Recall that the predictor is defined by (1.10). Because we deal with the one-step-ahead predictor, it is clear that the first estimator to be non singular is . The only variable over which the projection can be made is indeed . So, at this step (time ), the predictor is defined by , for and the criterium to minimize is . At the next step (time ), being available, one might project over the set but it is easy to check that this set is singular due to the fact that the predictor is proportional to . In other words, it is equivalent to project over the spaces or minimizing the criterium . We have chosen the former and therefore define . In fact, as soon as a new observation is available, only one new parameter is determined. So, during the first steps of the algorithm ( being the system order), one has to distinguish two kinds of projections:

• , then

• then

Clearly, the odd steps are quite identical to the recursive steps. For the even steps, one has to modify the general procedure. Note first that, for odd steps , one has, . As a consequence, implies that . The normalization matrix of the backward residuals is then singular according to (3.4) as noticed above. So, define the unormalized backwards residuals as follows:

(3.5)
where Definition 5 : first component of .

Clearly we have the orthogonal decomposition:

with respect to the inner product .

Let us now normalize these backward residuals. Formulas concerning the first residuals are obviously the same as for the recursive steps. Concerning the last one, that is

I S t i S t i ( ) T - S R n ∈ S 1 ≠ ⁄ ( ) I SS T - ( ) 1 2 ⁄ - ⇒ I SS T - ( ) 1 S T S - ( ) 1 2 ⁄ - 1 - ( )S T S ⁄ = 0 y ˆτ 2 τ , 2 3 , = { } y 1 t 2 = y ˆτ 2 P τ 2 y • y • 1 - [ ] = τ 2 3 , = y 2 y ˆ2 2 - ( ) 2 t 3 = y ˆ2 2 y 1 y 2 y ˆ2 2 , , { } y ˆ2 2 y 1 y 2 y ˆ2 2 , { } y 1 y 2 , { } y 2 y ˆ2 3 - ( ) 2 y 3 y ˆ3 3 - ( ) 2 + y ˆτ 3 P τ 3 y • y ˜• 1 - 2 y ˆ1 - 2 , [ ] = 2n n t 2k = y ˆτ t P τ t y • Z • 1 - t 1 -… Z • k - 1 + t k - 1 + y ˜• k - t k - , , , [ ] = t 2k 1 + = y ˆτ t P τ t y • Z • 1 - t 1 -… Z • k - t k - , , [ ] = t 2k 1 - = u τ k 2 -t , S t k 2 -V τ 1 - k 2 -t 1 - , = E t u • k 2 -t , ( ) 2 [ ] 1 = S t k 2 -2 1 = X • k 1 -t , X τ 0 t 1 - , Z τ t 1 - = X τ i t 1 - , X τ 1 - i 1 -t 2 - , P τ t 1 -X • 1 - i 1 -t 2 - , z • i 1 -t 1 - , [ ] - i ∀ , 1…k 2 - = = x τ k 1 -t 1 - , x τ 1 - k 2 -t 2 - , P τ t 1 -x • 1 - k 2 -t 2 - , z • k 2 -t 1 - , [ ] - =        x τ k 2 -t 2 - , ≡ X τ k 2 -t 2 - , Z • 1 - t 1 -… Z • k - 1 + t k - 1 + y ˜• k - t k - , , , { } sp X • 1 - 0 t 1 - , { } sp … X • 1 - k 2 -t 1 - , { } sp x • 1 - k 1 -t 1 - , { } sp ⊕ ⊕ ⊕ = 〈 | 〉 t 1 - k 1 -
, a particular step is necessary. Let be the normalized residual. The last equation of (3.5) allows to compute the normalization coefficient as:

(3.6)
Setting and assuming that lead to:

(3.7)
Let us note the component of as . using the fact that:

(3.8)
the matrix may be written:

(3.9)
So, one natural solution of equation (3.7) is:

(3.10)
Finally, we have:

(3.11)

Adaptativity.

The hereditary nature of our algorithm, which preserves all correlation terms, proves to be essential for accuracy. However, for adaptive purposes, a limit must be imposed upon heredity in order to follow-up possible drifts with time. This is achieved by imposing a finite past time-horizon upon the hereditary past. Note this horizon length. In this case, the hereditary part of the algorithm is modified as follows:

• While , the algorithm is the same that the one described in Figure 2.

For , the hereditary loop is computed form to and the formulas of the algorithm have to be modified as follows:

(3.12)

x τ k 1 -t 1 - , v τ k 1 -t 1 - , r t 1 - k 1 -x τ k 1 -t 1 - , = r t 1 - k 1 - v τ 1 - k 1 -t 1 - , r t 1 - k 1 - 1 0 R t 1 - k 2 - ( ) 1 - V τ 2 - k 2 -t 2 - , S t 2 - k 2 -u τ 1 - k 2 -t 1 - , - ( ) = ρ r t 1 - k 1 - 1 0 R t 1 - k 2 - ( ) 1 - = E t x τ 1 - k 1 -t 1 - , ( ) 2 [ ] 1 = ρ t 1 - t ----------I S t 1 - k 2 -S t 1 - k 2 - ( ) T - ( ) ρ T 1 = S S S ( ) 1 S ( ) 2 , [ ] T = S t 1 - k 2 -2 S t 1 - k 2 - ( ) 1 2 S t 1 - k 2 - ( ) 2
2 In a previous paper [START_REF] Monin | Exact ARMA Lattice Predictor From Autocorrelation Functions[END_REF], the authors have established a realization algorithm of the exact ARMA lattice predictor of order , from the autocorrelation function of the process. This is valid only when a precise model of the process exists (ideal spectra, fractional or distributed parameters systems,...).When only experimental data are available, constructing the ARMA lattice predictor of order from those data solves in particular the identification problem of an -th order linear system in ARMA form, when that is the real case.

+ 1 = = I SS T - I SS T - 1 S ( ) 1 2 - S ( ) -1 S ( ) 2 S ( ) -1 S ( ) 2 1 S ( ) 2 2 - S ( ) 2 2 S ( ) -1 S ( ) 2 S ( ) -1 S ( ) 2 S ( ) 1 2 = = S ( ) 2 - S ( ) 1 S ( ) 2 - S ( ) 1 = ρ t t 1 - ----------S t 2 - k 2 - - ( ) 2 S t 2 - k 2 - ( ) 1 , = v τ 1 - k 1 -t 1 - , t t 1 - ----------S t 2 - k 2 - - ( ) 2 S t 2 - k 2 - ( ) 1 , V τ 2 - k 2 -t 2 - , S t 2 - k 2 -u τ 1 - k 2 -t 1 - , - ( ) = N t N < t N ≥ τ t N - 1 + = t 1 - V τ 0 t , y τ y ˆτ t - ( ) E t y • 2 [ ] z • 0 t , t - ( ) ⁄ y ˆτ t z • 0 t , ⁄ = V τ i t , I S t i 1 -S t i 1 - ( ) T - [ ] 1 2 ⁄ - V τ 1 - i 1 -t 1 - , S t i 1 -u τ i 1 -t , - ( )i , 1…n
To show that, its suffices to take a note of:

The certainty equivalence property.

The lattice realization algorithm of [START_REF] Monin | Exact ARMA Lattice Predictor From Autocorrelation Functions[END_REF], and the lattice identification algorithm of the present paper are identical, up to the substitution:

(4.1)
Proof: It suffices to compare Figure 2 of the present paper with Figure 4 of [START_REF] Monin | Exact ARMA Lattice Predictor From Autocorrelation Functions[END_REF] (together with its formula 4.11,4.12 and 4.13).

Asymptotic convergence.

Proposition 1 . Whenever is an ergodic process in the sense, the hereditary algorithm converges in probability to the unique stationary optimal predictor.

Proof: When one compares the algorithms appearing in each of the papers mentioned above, with property 4.1 in mind, one immediately concludes by simple inspection that they both represent the same (state+parameters) hereditary transition map, with simply a different input. Indeed, we may write, on the one hand:

(4.2)
for the realization algorithm with , , and predictor as:

(4.3) E t y • 2 [ ] 1 N ---- y σ 2 σ t N - 1 + = t ∑ = E t 1 + y • V • 1 - i t , [ ] 1 N ---- y σ V σ 1 - i t , σ t N - 2 + = t 1 + ∑ = n n n E y t V t 1 - i t , [ ] E t 1 + y • V • 1 - i t , [ ] → y t { } L 2 v t p t z t 0 r v t 1 - p t 1 - z t 1 - 0 E y t y t τ - [ ] { } 0 τ t < ≤ ,             = v R 2n t × ∈ p R 2n ∈ z 0 R ∈ y ˆt y ˆt l v t p t z t 0 , , ( ) =
and, on the other hand:

(4.4)
for the identification algorithm with , , and predictor as:

(4.5)
In both cases, and are the same maps, from which it results that the input-output maps of both algorithms are the same. They are rational (containing only the four arithmetical operations) and continuous (all inversions are non-singular by definition). Note that the symbols , and have been used, instead of , and in both cases, to stress the fact that they are not numerically identical for finite .

From [START_REF] Monin | Exact ARMA Lattice Predictor From Autocorrelation Functions[END_REF], system (4.2) is the exact time-varying solution, for all initial states, of the optimal predictor. Since ergodicity implies stationarity, the asymptotic solution of (4.2) is, by definition, the unique value of the stationary optimal predictor. System (4.4)-(4.5) is identical to system (4.2)-(4.3), except for its input signal instead of . One has to prove that the asymptotic output of (4.4)-(4.5) tends to the asymptotic output of (4.2)-(4.3) in a suitable topology.

Let us choose the -ergodic topology for the input, i. e.: iff .

Let be the common input-output map of both algorithms. By continuity of , we have the following diagram:

As -erg convergence implies convergence in probability, the proof is completed.

Examples and comparisons.

As was announced in the abstract, the algorithm proposed here is shown to be more accurate than classical solutions, in term of spectral resolution, especially in the transient phase.

Spectral resolution

. V t S t Z t 0 r V t 1 - S t 1 - Z t 1 - 0 E t y • y • τ - [ ] { } 0 τ t < ≤ ,             = V R 2n t × ∈ S R 2n ∈ Z 0 R ∈ Y ˆt Y ˆt l V t S t Z t 0 , , ( ) = r l Y V Z y ˆv z t E t y • y • τ - [ ] E y t y t τ - [ ] L 2 L 2 -erg E t y • y • τ - [ ] → E y t y t τ - [ ] 1 T ---E t y • y • τ - [ ] E y t y t τ - [ ] - 2 τ d 0 T ∫ 0 → Φ Φ Φ E y t y t τ - [ ] → y ˆt { } ↑L 2 -erg ↑L 2 -erg Φ E t y • y • τ - [ ] → Y ˆt { } L 2
To illustrate this property, we have first chosen a simple 4-dimensional linear stochastic system whose poles and zeros are represented in Figure 7 (a). The module of the corresponding transfer function (square root of the spectral density) is represented from 0 to the half of the sampling frequency in Figure 7 (b). Figure 4 shows the transient spectrum obtained with finitely recursive algorithms such as the "Extended Least Square" and "Recursive Maximum Likelihood" algorithms [START_REF] Ljung | System Identification: Theory for the User[END_REF]- [START_REF] Michaut | Méthodes Adaptatives pour le Signal[END_REF]. It appears clearly that these algorithms do not separate the two distinct poles of the system, even after 2000 iterations. Note that, for such low-order examples as above, the lattice character of our algorithm is of little concern. The differences mainly lie in frequency resolution. In Figure 9, we exhibit the power spectra obtained with the classical Gauss-Newton Block Data algorithm (GNBD) [START_REF] Ljung | System Identification: Theory for the User[END_REF], from the MATLAB identification tool-box, compared with the one obtained using the Hereditary Least Square algorithm (HLS). We have run the GNBD algorithm for successive data sequences of length 50, 150, 200... Although the GNBD algorithm converges to the real power spectrum of the signal (after a number of steps greater than 750), it clearly appears that the hereditary algorithm detects the two poles much faster (after 150 steps only). In order to underline this fact, we represent in Figure 6 all the power spectra (more precisely, their square roots) obtained after 150 iterations with all the methods we tested.

As is well known, it is the more difficult to distinguish two resonance modes as they ar nearer to each other. We therefore have chosen, as a second example, the linear stochastic system described in Figure 7 through its pole/zero representation. We see in Figure 9 that the GNBD algorithm does not separate the two frequencies for a data length lower than 2000, which is not the case with our algorithm which early detects them. (Note that the frequency scale has been expanded to make sure a better reading).

Concerning the finitely recursive algorithms (ELS and RML), it appears clearly in Figure 9 that they do not distinguish between the two frequencies. The lack of resolution on the part of ELS and RML algorithms is due to the approximations which are inherent to the usual finitely recursive algorithms, whose consequence is to skip higher interval correlation terms.

Adaptativity.

The signal to be identified consists of four frequencies with a drift along time which contains information. This is the case in the SONAR passive estimation, the drift being due to the Doppler effect. A simulation is presented in Figure 10 and shows the accuracy of our algorithm to match the emitted spectrum. The dotted line represents the spectrum to be identified and the continuous line the results of our HLS algorithm.

Due to the superior transient capabilities of the HLS algorithm, a local time-horizon of 500 appears to be sufficient to follow-up the drift while preserving accuracy on modes which are closed to each other. It should be noticed, in the above example, that the lattice approach is necessary from a numerical point of view.

Order reduction.

Finally, we exhibit the performance of our algorithm when order reduction is required. As is well known, a good test for spectral parameters identification is to choose a signal spectrum with abrupt slopes. This implies higher order whose low order approximations are difficult to get in an optimal way [START_REF] Monin | Exact ARMA Lattice Predictor From Autocorrelation Functions[END_REF]. We simulated a low-pass filter of dimension 20 (synthesized with the algorithm of [START_REF] Monin | Exact ARMA Lattice Predictor From Autocorrelation Functions[END_REF]) and use it to generate the signal to be identified with a 4-dimensional system. The results appear in Figure 11, Figure 12 and Figure 13 and confirm again the accuracy of the hereditary algorithm.

Conclusion.

We have derived in this work a lattice solution to the identification problem of the optimal linear predictor of ARMA type, which leads to a recursive algorithm with hereditary terms. These terms, which are inherent to the statistical optimization, are due to the increasing memory structure of higher interval correlation computations which are essential for an accurate spectral resolution. Heredity consists in preserving such terms which implies a growing memory algorithm as opposed to more usual approximations in view of finite recursiveness. The projection technique introduced produces a minimum phase linear model of the signal whose n-th order "whiteness" of the associated innovation has the following restricted meaning: orthogonality to an n-dimensional subspace memory of the past in a suitable Hilbert sequence space. This is suitable for stationary signals, in which case we have shown the optimality of the hereditary algorithm for fixed samples of data as well as its asymptotic convergence, for unlimited records. In non-stationary situations, the algorithm exhibits also superior adaptive properties, when a limited window over the past is considered. Finally, we give an example of the better parametric properties of the hereditary algorithm, when trying to fit a low order predictor to higher dimensional spectrum such as the ideal low-pass. For t 1 =

Computation of the backward residuals

V τ i t , τ , 1…t 1 - = V τ i t , t t 1 + -----------     I S t i 1 -S t i 1 - ( ) T - ( ) 1 2 ⁄ - V τ 1 - i 1 -t 1 - , S t i 1 -u τ i 1 -t , - ( ) = V τ 0 t , y τ u τ 0 t , z • 0 t , t - ( ) t t 1 + -----------E t y • 2 [ ] z • 0 t , t - ( )     ⁄ u τ 0 t , t t 1 + ----------- ⁄ = u τ i t , 1 S t i ( ) T S t i - ( ) 1 2 ⁄ u τ i 1 + t , S t i ( ) T V τ 1 - i t 1 - , + = For i 1 = to to T For τ 1 = to t For i n 2 - = down to 0 n 1 - u τ n 1 -t , S t n 1 - ( ) T V τ 1 - n 1 -t 1 - , = Computation of the parameters S t 1 + i S t 1 + i E t 1 + V • 1 - i t , y • [ ] z • i t 1 + , t 1 + ⁄ = Computation of the means E t 1 + y • V • 1 - i t , [ ]
For i 0 = to n 1 - For i 0 = to n 1 - Hereditary terms E t 1 + y • V • 1 - i t , [ ] 1 t 1 + ----------- y σ V σ 1 - i t , σ 1 = t 1 + ∑ = z • i t 1 + , t 1 + E t 1 + y • V • 1 - k t , [ ] 2 k i = n 1 - ∑ =

(2. 22 )

 22 If one defines the projection coefficients as: , (2.23) the relation (2.21) takes the following form:

Figure 1

 1 Figure 1 illustrates the involved lattice structure with and

Figure 2 -

 2 Figure 2 -The hereditary algorithm.

Figure 3 -

 3 Figure 3 -Poles (x) and zeros (o) location (a) and power spectrum (b) (example 1).

Figure 4 -Figure 5 -Figure 6 -

 456 Figure 4 -Spectral density obtained with classical finitely recursive identification algorithms (example 1)

Figure 7 -

 7 Figure 7 -Poles (x) and zeros (o) location (a) and power spectrum (b) (example 2).

Figure 8 -

 8 Figure 8 -Transient spectral densities obtained with Gauss-Newton and hereditary algorithms (example 2).

Figure 9 -

 9 Figure 9 -Transient spectral densities obtained with finitely recursive algorithms (example 2).

Figure 10 -Figure 11 -Figure 12 -

 101112 Figure 10 -Non stationary situation (example 3).

Figure 13 -Figure 1 -Figure 2 -

 1312 Figure 13 -Final results after 2000 steps.

Figure 3 -

 3 Figure 3 -Poles (x) and zeros (o) location (a) and power spectrum (b) (example 1).

Figure 4 - 1 )Figure 5 -

 415 Figure 4 -Spectral density obtained with classical finitely recursive identification algorithms (example 1) Figure 5 -Transient spectral densities obtained with Gauss-Newton and hereditary algorithms (example 1).

Figure 6 -

 6 Figure 6 -Power spectra obtained for a 150 long sequence (example 1).

Figure 7 -

 7 Figure 7 -Poles (x) and zeros (o) location (a) and power spectrum (b) (example 2).

Figure 8 -

 8 Figure 8 -Transient spectral densities obtained with Gauss-Newton and hereditary algorithms (example 2).

Figure 9 -

 9 Figure 9 -Transient spectral densities obtained with finitely recursive algorithms (example 2).

Figure 10 -

 10 Figure 10 -Non stationary situation (example 3).

Figure 11 -Figure 12 -

 1112 Figure 11 -The low-pass filter approximations with ELS and RML Figure 12 -The low-pass filter approximations with HLS and GNBD.

Figure 13 -

 13 Figure 13 -Final results after 2000 steps.

Figure and table