Equivalent prelter classes for least-squares estimates in prediction error identication

Equivalent prelter classes for least-squares estimates in prediction error identication

P. Carrette

Department of Electrical Engineering Linkping University, S-581 83 Linkping, Sweden

WWW:

Email:

April 9, 1998

REGLERTEKNIK

AUTOMATIC CONTROL LINKÖPING

Report no.: LiTH-ISY-R-2023

Submitted to IEEE Trans. Signal Processing and ECC-99

Technical reports from the Automatic Control group in Linkping are available by anony-

mous ftp at the address

. This report is contained in the com-

pressed postscript le

.

Equivalent prelter classes for least-squares estimates in prediction error identication

P. Carrette - carrette@isy.liu.se

Department of Electrical Engineering, Linkoping University S-58183 Linkoping, Sweden

Phone: ++46-13-281311 and Fax: ++46-13-282622

Abstract

We present results about classes of prelters that may result in similar model estimates by use of the classical least-squares prediction error identication methods.

It will be seen that prelters belonging to a particular class have identical transfer function moduli. We also address the issue of a possibly best prelter within such equivalent classes. The result we obtain actually advocates for prelters having the shortest extended impulse response (dened even for IIR prelters).

Keywords

: prelters, least-squares estimates, prediction identication

1 Introduction

Black-box system identication deals with the estimation of parametric models for dy- namical systems based on input-to-output (I/O) data measured on these systems during particular operating conditions.

In linear system identication, one often postulates a particular relation between the sys- tem input and output measurements in terms of linear operators and output corrupting perturbations, i.e.

y ( t ) = G

( q ) u ( t ) + H

( q ) e ( t ) (1) where u ( t ) and y ( t ) respectively denote the input and output of the system while e ( t ) stands for perturbation samples, e.g. out of a white noise process with bounded moments.

The dynamics of the system lies in the linear operators G

( q ) and H

( q ) that are called the I/O and noise-to-output (N/O) system dynamics. Here, we only consider discrete time dynamics so that these operators are digital lters 1]. As no closed loop system conguration will be considered, we will assume these lters to be stable. Furthermore, note that H

( q ) is often postulated monic and inversely stable.

This work was supported by the Swedish Research Council for Engineering Sciences.

1

When performing the identication of the system in (1), it is classical to look for a parametric dynamical model whose output prediction 2] will be made as close as possible to the system output measurements. The model dynamics is chosen as

y ( t ) = G ( q ) u ( t ) + H ( q ) " ( t )

where G ( q ) and H ( q ) are parameterized linear dynamics while " ( t ) stands for the model prediction error. This prediction error, denoted " ( t ), is identical to the deviation between the system output and the prediction of the model output ^ y ( t ) that is written as y ^ ( t ) = H

( q ) G ( q ) u ( t ) + ( H ( q )

1) y ( t )] (2) Note also that the N/O model dynamics, i.e. H ( q ), is often restricted to be monic while it has to be inversely stable for leading to a bounded model output prediction.

Now, given a set of I/O data measurements, i.e.

Z

^N =

( y ( t ) u ( t ))  t

= 1 N ]



and a particular choice of the model dynamics complexity, the prediction error methods (see 2] as a reference book) intend to estimate the model parameters within the vector



ⁿ so that a particular measure of the prediction error sequence is minimized.

Note that the parameter vector subset

is usually related to stable predictor dynamics.

Then, let us introduce the zero-extension of this original measurement set as

Z ^N =

( y ( t ) u ( t ))  t

and (0  0)  otherwise

(3) Thus, the estimated parameter vector associated to the time set

is written as follows :

 ^ (

) = argmin



V ( Z ^N 

) (4)

for which the far most often used identication cost function is V ( Z ^N 

) = 1 N

X

t

"

_m ( t )

where " m ( t ) stands for the preltered prediction error \signal" based on the data set Z ^N , i.e. " m ( t ) = M ( q ) " ( t ) with M ( q ) denoting a stable prelter.

Note that, in case

=

, the estimated parameter vector is the classical least square (LS) prediction error estimate (with prelter M ( q )).

In this paper, we deal with the prediction error preltering issue while analyzing the impact of particular choices of the time set

(accordingly to the prelter M ( q ) and, possibly, to  ) on the results of the system identication as presented above. More pre- cisely, we show that, when this time set is extended from the original set, i.e.

, to a set for which the preltered prediction errors are no more corrupted by ltering edge eects,

2

the identication cost function becomes invariant over the class of prelters having the same transfer function modulus. This invariance property obviously propagates to the associated model parameter estimate. Furthermore, while restricting the time set

from the original time set, it is possible to distinguish best prelters out of a given equivalent prelter class. They are dened as the prelters for which the set of preltered prediction errors that are free from ltering edge eects is the most extended.

Therefore, we proceed as follows. In Section 2, we recall some properties of linear lter- ing as well as of ltered signal energy. In Section 3, we apply these ltering properties to our prediction identication framework. Then we derive the characterization of the equivalent class of prelters (associated to appropriate times sets) over which the related identication cost function is invariant. In Section 4, we evaluate the deviation between the dierent identication cost functions within such an equivalent prelter class. Then, given prelters out of this class, we derive time sets that lead to \edge eect"-free predic- tion error sets. The accuracy of the corresponding parameter estimates gives rise to the notion of best prelters out of that class. Finally, some simulations are given in Section 5 in order to illustrate the obtained results.

2 Some discrete time ltering properties

In this section, we intend to recall properties of discrete time ltering while putting into light problems occuring with nite length signals 1].

For the purpose of clarity, let us introduce our building blocks. There are three entities we will be dealing with.

Firstly, the (innite length) sequence s , i.e. s =

s ( t )  t

with

the set of integers, constitutes the zero-extension of a nite length sequence s

, i.e. s

=

s

( t )  t

. Its samples satisfy

s ( t ) =



s

( t ) if t

0 otherwize

Secondly, the sequence s (

) is the zero-extension of the restriction of the sequence s to the time interval

. Thus, its samples satisfy s ( t

) = s ( t ) for t

and zero otherwise.

Obviously, for

, we have s (

) = s .

In the sequel, such (innite length) sequence is also called a \signal" with reference to the time argument, e.g. the signal s ( t ) for the sequence s .

Thirdly, the lter F ( q ) has a L f -length impulse response with real bounded coecients, i.e.

F ( q ) = q ^{d f L}

^f

i

f i q

ⁱ (5)

where d f stands for a nite delay corresponding to a d f -forward time-shift operator.

Remark 1 This lter formulation may also hold in case of a stable innite impulse re- sponse (IIR) lter provided its \eective" length is less than L f with appropriate denition

3

of the delay d f , i.e. F ( q ) = q ^{d f}

_i f i q

ⁱ and L f

f j f j

k

X

0

i<L f f _i

i f _i

for any integer j provided

k

L f . The stability of the lter implies that the most right

hand side of these inequalities is nite.

Now, let us evaluate the ltered signals conformably to classical signal processing theory.

For simplicity, we rst consider the ltering of the signal s ( t ). We have s f ( t ) = F ( q ) s ( t ) (with abuse of notation)

= ^{L f}

;1

X

i

f i s ( t + d f

i )

Obviously, the sequence s f is still the zero-extension of a nite sequence because s f ( k ) = 0 for

k

large enough. In fact, we can state that

T T

f = 1

d f N + ( L f

1)

d f ]

s f (

) = s f

Thus, s f (

) = s f only in the case d f = 0 and L f = 1 that corresponds to a trivial dynamics for the lter, i.e. F ( q ) = f

.

In order to be complete, it is interesting to put into light the shortest time interval

for which the ltering of the signal s ( t

) by F ( q ) is identical to the corresponding signal s f ( t

). As s (

) = s and s f (

f ) = s f , the set

is simply obtained by the union of these two time intervals. This is stated in the following lemma.

Lemma 1 Let the signal s ( t ) be nonzero only for t

and let the lter F ( q ) have a (forward) delay d f and a length L f . Then,

T T



)

s f ( t

) = F ( q ) s ( t

) 

t

where

= 1

 ( d f ) N +  (( L f

1)

d f )] with  ( x ) = max(0 x ).

In Figure 1, we have drawn the dierent time sets that appear in the derivations.

Note that the time set

=

( d f L f N ) has the smallest number of elements when d f

0 L f ) for which we have ] (

) = N + ( L f

1).

With this result, we are ready to turn to energy considerations. Recall that the energy of any (innite length) sequence is given by its 2-norm square as

s

=

_t

s

( t ) where the samples s ( t ) are indeed dened over the whole time domain.

Here, it is easily seen that

s (

)

=

s

for any

while

s f (

)

=

s f

for any

f . Hence, by the denition of

, we have

k

s (

)

=

s

and

s f (

)

=

s f

(6)

4

-

r r

1

d f N + ( L f

1)

d f

r r

1 N

 r

 r

- r

d f d f

L f

1

Figure 1: The time set

( d f L f N ) (thick line) in the case d f

0 L f ).

The last step of our derivations is now to establish the link between the energy of the

ltered signals and the characteristics of the lter F ( q ) itself. Therefore, we use Parseval's identity relating the energy of any (innite length) sequence expressed in terms of its time representation and of its z -transform representation (see e.g. 1, Chap. 4]). It is written as

k

s f

= 12 

Z



;



F ( e ^j! )

S ( e ^j! )

d!

where S ( z ) =

_t

s ( t ) z

^t stands for the z -transform of s ( t ) with z = e ^j! . Then, by use of (6), we obtain

k

s (

)

= 12 

Z



;



F ( e ^j! )

S ( e ^j!

)

d! (7) where S ( z

) denotes the z -transform of s ( t

) (identical to S ( z )).

In fact, this last relation constitutes the basis for the results that will be derived in the next sections. Let us comment it in order to put into light its characteristics:

This relation shows up two dierent ways of expressing the energy of a signal, i.e.

s f ( t

), for which the well-known edge eects due to ltering has been removed by the appropriate choice of the time set

, i.e.

=

( d f L f N ).

In the right hand side of the Parseval's identity (7), we recognize the square modulus of the transfer function of the lter F ( q ). This means that this lter character is the basic issue for this relation to hold.

As the modulus of the lter transfer function, i.e.

F ( z )

, is invariant under delay operation, it is obviously seen that the lter delay d f only aects the time set extension of any

( d f L f N ) (see Lemma 1).

As a partial conclusion, it is interesting to note that if two lters, i.e. F i ( q ) (for i = 1  2) with delay d i and length L f i , have the same transfer function modulus then the energies of the corresponding ltered signals, i.e. s f i ( t

), are identical provided an appropriate denition of the time set

, i.e.

with

_i

=

( d f i L f i N ).

Finally, it is interesting to evaluate the impact of the edge eects due to ltering within

5

the original time set

. More precisely, we intend to answer the following question: in which set

lie the time instants for which the nontrivial samples of the signal s f are found in the signal corresponding to the ltering of the original (nite length) signal s

( t ) by F ( q ) disregarding its initial conditions? In fact, this time set

is the largest for which the samples s f ( t ) are not inuenced by the zero-extension we have imposed.

Therefore, the following result holds.

Lemma 2 Let the signal s ( t ) be nonzero only for t

and let the lter F ( q ) have a (forward) delay d f and a length L f . Then,

T T )

s f ( t ) = 0  t

where

= 1 +  (( L f

1)

d f ) N

 ( d f )].

In Figure 2, we have schematically drawn the time set

=

( d f L f N ).

-

 r

1 + ( L f

1)

d f N

d f

r r

1 N

 r

 r

- r

d f

d f

L f

1

Figure 2: The time set

( d f L f N ) (thick line) in the case d f

0 L f ).

Note that it can be called the \edge eect"-free time set. Furthermore, it has the largest number of elements when d f

0 L f ) for which we have ] (

) = N

( L f

1). It obviously depends on the length of the lter.

3 Prediction prelter equivalence class

In this section, we apply the result obtained in expression (7) to the prediction error identication cost function presented in the introduction. This will allow us to dene a time set

_p

for which this identication cost function, i.e. V ( Z ^N 

_p

), becomes invari- ant under prelters, i.e. generically denoted M ( q ), belonging to an appropriately dened prelter class that we intend to characterize.

In the framework of least squares prediction error identication, it is seen that the iden- tication cost function is nothing else but the mean energy of the preltered prediction error samples, i.e. " m ( t ) = M ( q ) " ( t ), evaluated over a subset of the zero-extended data set Z ^N , i.e.

V ( Z ^N 

) = 1 N

X

t

"

_m ( t )  (8)

6

Before elaborating on this expression, let us analyze the properties of the lter that links these preltered prediction error \signal" to the samples taken from this data set.

First, as it appears from the expression (2), the prediction error \signal" originates from the ltering of the zero-extended input and output signals by the predictor dynamics.

More precisely, we have

" ( t ) = H

( q ) y ( t )

H

( q ) G ( q ) u ( t )

From this expression, let d p (  ) and L p (  ) respectively denote the delay and the length characteristics of the class of lters (see the formal denition in (5)) to which these two prediction dynamics belong. Note that the causality of the model dynamics implies that the delay is nonpositive, i.e. d p (  )

0. In fact, due to the monotonicity of the model noise dynamics, this delay d p (  ) is zero whenever this dynamics is present.

Furthermore, the quantity L p (  ) is identical to the length of the most extended impulse response between that of \ H

( q )" and of \ H

( q ) G ( q )" at  .

Note that it is valuable to extend these two prediction lter characters to a subset of the parameter vector  set, i.e.

ⁿ . This leads to

d p =

min _

2D

j

d p (  )

and L p = max _

2D

L p (  ) (9)

For example, in the case of ARX models, the delay d p cancels out while the length L p

only depends on the model complexity. Indeed, if A ( q ) and B ( q ) denote the ARX polynomials, i.e. H

( q ) = A ( q ) and  H

( q ) G ( q )] = B ( q ) with B ( q ) having a (negative) delay d b , then it is easily found that

( L p (  ) =) L p = max( L a L b +

d b

1)

where L a and L b are the numbers of parameters in A ( q ) and B ( q ), respectively.

Secondly, the preltered prediction error \signal" " m ( t ) is obtained by ltering the above prediction error \signal" " ( t ) by the (generic) prelter M ( q ) whose characteristics are a delay d and a length L , say. In other words, the preltered prediction errors are related to the (zero-extended) input and output signals by global lters F ( q ) that can be described with a delay d f and a length L f respectively identical to

d f = d + d p and L f = L + ( L p

1) for any 

.

Now, we are ready to apply the results derived in the expression (7) to the identication cost function in (8). This leads to the following theorem.

Theorem 3 Let the zero-extended data set Z ^N and the prediction lter characters d p

and L p (see expression (9)) be given. Then, provided 

and

_p

( L ) with

T

p

( L ) = 1

 ( L

d p

) N + ( L + L p

2) +

d p

] 

7

the quantity V ( Z ^N 

) is invariant under prelters M ( q ) belonging to the equivalent class

( L

M

) dened as

M

( LX ) =

(

M ( q ) = q ^{d L}

i

m i q

ⁱ











d

0 L ) 

M ( z )

= X ( z ) for

z

= 1

)

where X ( z ) is a positive function over

z

= 1 that represents the modulus of the transfer function of any prelter belonging to this class.

Proof: First, the prelters belonging to the equivalent class

( L

M

) have characters d and L with d

0 L ).

Then, with F

( q ) = M ( q ) H

( q ) and F

( q ) = M ( q )  H

( q ) G ( q )], we have that

T

p

( L )

( d p + Ld f

(  ) L f

(  ) N )

( d f

(  ) L f

(  ) N )

as d p

0, d p

d f i (  )

L + d p and L f i

L + ( L p

1) (with i = 1  2) for any 

. This implies that

" m ( t

) = F

( q ) y ( t

)

F

( q ) u ( t

) for

_p

( L ).

Hence, the result in expression (7) can be applied. It straightforwardly yields N 1

X

t

"

_m ( t ) = 1 2  N

Z



;



j

M ( e ^j! )

j

H ( e ^j!  )





Y ( e ^j!

)

G ( e ^j!  ) U ( e ^j!

)

d!

where U ( z

) (resp. Y ( z

)) denotes the z -transform of the (innite length) sequence u (

) (resp. y (

)). The left hand side of this expression is identical to the identication cost function on the considered time set, i.e. V ( Z ^N 

), while its right hand side is obviously invariant under prelters having the same transfer function modulus with d

0 L ), i.e. M ( q )

( L

M

).

In Figure 3, we have depicted the time set

_p

( L ) =

_p

( Ld p L p N ) for which the delay d p

-

r r

1

( L

d p

) N + ( L + L d

2) +

d p

r r

1 N

- r

 r

- r

- r

j

d p

L

j

d p

L + L p

2

Figure 3: The time set

_p

( Ld p L p N ) (thick line) in the case

d p

L .

is nonpositive. As mentioned above, this delay is zero whenever a model noise dynamics

8

is present.

It is interesting to note that any linear phase non-causal prelter M ( q ) belongs to an equivalent class

( L

M

) that at least contains a causal prelter, provided the impulse response length L is chosen appropriately. For example, if we consider M ( q ) = R ( q ) R ( q

) to be a zero-phase non-causal prelter with R ( q ) a L r -length causal lter then M

( q ) = R ( q ) R ( q ) is a causal lter having a (2 L r

1)-length impulse response

. Thus, we obviously have that these two prelters belong to

(2 L r

1 

R

).

As an immediate consequence of Theorem 3, the model parameter estimate resulting from the minimization of the related identication cost function, i.e.

 ^ (

) = argmin _

2D

V ( Z ^N 

) for

_p



inherits the invariance property over the equivalent class

( L

M

) of prediction error prelters.

Finally, let us consider the situation where the time set on which the identication cost is evaluated is identical to the original interval

. Due to the fact that L > 0, this means that the invariance of this cost value with respect to prelters M ( q ) belonging to the class

M

( L

M

) does not hold anymore. Nevertheless, the time set

_p

( L ) can still be used as an intermediate for inferring the inuences of such prelters onto the corresponding identication cost values.

It is actually seen that the possible dierence between two such values becomes marginal for increasing numbers of data measurements. Indeed, it is easy to derive that

V ( Z ^N 

) = V ( Z ^N 

_p

( L ))

1 N

X

t

_p

L

"

_m ( t )

Obviously, the contribution of the second term in the right hand side has a relatively decreasing impact on the cost value for increasing number of data because it remains generated by a constant number of preltered prediction error samples that is

#(

_p

( L )

)

2 L + L p

2

This then implies that the identication cost function is asymptotically invariant under the application of prelters belonging to the same equivalent class. So will be the corre- sponding model parameter estimates, i.e. ^  (

).

As a conclusion, we can mentioned that this is a well-known result from the frequency do- main interpretation of the model parameter estimation while considering the asymptotic situation of an innite number of data measurements 2, chap. 8].

1

These two prelters are easily implemented by use of \

" and by the double use of \

"

(in Matlab) with the lter

(

), respectively.

9

4 Best prelter within an equivalent class

The last part of the preltering issue deals with the deviation between dierent least- squares estimates ^  (

) (from the minimization in (4)) obtained while using prelters belonging to the same equivalent class, i.e. M ( q )

( LX ), say.

This issue is related to the question of the existence of a possibly best prelter given an equivalent class. Actually, this question is rather tricky to answer because we must state what best means under this preltering issue. Here, we will consider the following assumption as a starting point for the comparison between dierent parameter estimates.

Assumption 1 There exists a monic (causal) prelter M

( q ) dening the equivalent class

M

( L

M

) so that

G ( q

) = G

( q ) as well as H ( q

) = H

( q ) M

( q )

for 

ⁿ .

This assumption refers to the fact that the system dynamics belongs to the preltered model set. In that case, the prelter M

( q ) can be seen as a priori known contributions to the N/O system dynamics, e.g. pole locations.

Then, we may look for the prelter M ( q ) belonging to the equivalent class

( L

M

) such that the corresponding model parameter estimate is the closest to 

in some sense.

First of all, it is important to note that some edge eects remain hidden in the origi- nal identication cost V ( Z ^N 

). By that we mean that we cannot pretend that the restriction of the sequence " m ( 

) to the original time set, i.e.

" m ( t

)  t

, is a part of the sequence that would have been generated by an innitely long original data set (except in case initial conditions are known when applying the prelter as well as of the prediction lters to the data measurements). In fact, only a restriction of that subsequence is!

With the help of Lemma 2, we can state the following result whose proof is left to the reader.

Theorem 4 Let the zero-extended data set Z ^N and the prediction lter characters d p

and L p (see expression (9)) be given. Let the prelter M ( q ) have delay d and length L with d

0 L ). Then, provided 

,

T T

p ( dL )

" m ( t ) = M ( q ) H

( q ) y ( t )

G ( q ) u ( t )]









i:c:  t

where

p ( dL ) = 1 + ( L + L p

2)

( d

d p

) N

 ( d

d p

)] and \

i:c: " means that correct initial conditions are used in the prediction lters and prelter application.

10

-

 r

1 + ( L + L p

2)

( d

d p

) N

( d

d p

)

r r

1 N

- r

 r

-

r

- r

 r

j

d p

d p

d L + L p

2 d

Figure 4: The time set

p ( dLd p L p N ) (thick line) in the case

d p

d

0 L ).

In Figure 4, we have represented the time set

p ( dL ) =

p ( dLd p L p N ) in the case

j

d p

d

0 L ). Note that in this situation, the number of elements of this set is iden- tical to ] (

p ( dL )) = N

( L + L p

2). Hence, the smaller the length of the prelter, the larger that number. Surprisingly enough, it is independent of the delays within the prediction lters and the prelters.

Now, let us evaluate the expectation (over the characteristics of the system white noise) of the identication cost function for this particular time set in situations where Assump- tion 1 holds. It is shown that, provided the prelter is taken in the equivalent class

M

( L

M

) with L larger than the length of the prelter itself, this expectation cost is minimum at 

and its value only depends on the length of the time set corresponding to the chosen prelter.

More precisely, the following results holds.

Theorem 5 Let the original data set Z

^N be given. Let Assumption 1 holds and the pre-

lter M ( q ) (having a delay d m and a length L m ) belongs to

( L

M

). Then, provided



,

E ( V ( Z ^N 

pm ))

E ( V ( 

Z ^N 

pm )) = #(

^pm )

N

where

pm =

p ( d m L m ) (for short) and

is the variance of the system white-noise.

Proof: ^{for }

, the expectation of the identication cost function is written as E ( V ( Z ^N 

pm )) = 1 N

X

t

pm

E ( "

_m ( t ))

= 1 N

X

t

pm

n



M ( q ) H

( q )( G

( q )

G ( q )) u ( t )

+ E

M ( q ) H

( q ) H ( q

) M

( q ) e ( t )

because the system input and noise sequences are decorrelated, i.e. E ( u ( t ) e ( t + k )) = 0 for all k . The rst term in the last expression is minimum at 

for which it cancels out. The

11

second term is more tricky to analyze but it is also seen to be minimum at 

while using the ergodic property of the system white-noise (in order to make the time set innitely large) as well as the derivation in 2, p. 226] after getting to frequency domain.

As a consequence, the quantity E ( V ( Z ^N 

pm )) is minimum at 

for which it is iden- tical to

E ( V ( 

Z ^N 

pm )) = #(

^pm )

N

This completes the proof.

For what concerns the model parameter estimate, we can mention that the classical asymp- totic (in N ) characterization of the statistical moments of the minimizing solution of prediction identication cost function 2, chap. 8-9] holds. That is to say that

N lim

;

V ( Z ^N 

pm )

E ( V ( Z ^N 

pm ))

= 0

so that

#(

pm )

⁼

 (

pm )



AsN (0 

) as N

(10) where

is a covariance matrix that is independent of N and invariant under the equivalent class

( L

M

) for any L

L m .

With this result and Theorem 4 in mind, it is possible to answer the question relative to the dierent performance of prelters belonging to a particular prelter in

( L

M

) with respect to the model parameter estimation.

This actually lies in the number of elements in the time set

p ( d m L m ). As this number increases for smaller L m (provided

d p

d m

0 L m ), no additional delay in M ( q )), we can state that prelters M ( q ) having shorter impulse response (eective) extents lead to more accurate model parameter estimates in view of expression (10). Then, we would say that

the best prelter belonging to a given equivalent class is the one whose impulse response has the shortest extent.

Of course, nding this \best" prelter is an optimization problem that minimizes a time domain character (its impulse response length) under a frequency domain constraint (its transfer function modulus). Furthermore, the prelter order is part of the optimization.

5 Simulations

The aim of this simulation section is twofold:

Illustrate the fact that the time set

_p

( L ) derived in Theorem 3 leads to model parameter estimates ^  (

_p

( L )) that are invariant under the use of prelters M ( q ) belonging to the same equivalent class, i.e.

( L

M

).

12

Demonstrate the expression for the accuracy of the model parameter estimate

 ^ (

pm ), i.e. minimizing solution of V ( Z ^N 

pm ) dened over the associated \edge eect"-free time interval

pm , for prelters M ( q ) having dierent impulse response (eective) length in case M

( q ) = 1 in Assumption 1.

Therefore, the system dynamics is chosen as

(1

1 : 5 q

+ 0 : 7 q

) y ( t ) = (1 + 0 : 5 q

) u ( t

1) + e ( t )

which is an ARX dynamics. The input signal is a pseudo-random binary sequence having values in the set

1  1

while the system noise process is a Gaussian white-noise with variance

= 0 : 01. The original data set Z

^N containing N = 1000 I/O measurements.

The model dynamics is chosen as a (2  2  1)-ARX structure so that the true parameter vector is identical to 

= 

1 : 5  0 : 7  1  0 : 5] ^T (in order to satisfy Assumption 1 with M

( q ) = 1). The number of parameter of this model structure is n = 4. Note that any prediction lter has a delay d p that is zero and a length L p that is identical to

L p = max( L a L b +

d b

1) = max(3  2) = 3 independently of  because of the ARX structure.

First, we consider the application of two dierent prelters belonging to the same class.

These prelters are dened as follows. Let R ( q ) be a 4-th order Butterworth lter with a cut-o frequency identical to 0 : 2 times the half sampling frequency

. Then, the two prelters are M

( q ) = R ( q ) R ( q

) (being zero-phase non-causal) and M

( q ) = R ( q ) R ( q ) (being causal), respectively. The impulse response of these two prelters are given in Fig- ure 5. The corresponding equivalent class is

( L

R

) with L chosen as an upper-bound for the (eective) extend of these two impulse responses. Under a tolerance identical to 10

, we take L = 220. For what concerns the prelter delay, the zero-phase property of M

( q ) implies d

= L= 2 while the causality of M

( q ) imposes d

= 0.

Now, the time set

_p

( L ) is evaluated as

T

p

( L ) = 1

 ( L

d p

) N + ( L + L p

2) +

d p

] = 

219  1221]

for L = 220. The estimation of the model parameter estimate, i.e. ^  (

_p

( L )), correspond- ing to the application of the prelter M i ( q ) is denoted ^  _pi

and is evaluated by use of an ARX routine that does not take care of the edge eects generated by the predictor dynamics (as these have been eliminated by the extension of the data set). We obtain the following estimates:

 ^

_pi = 

1 : 511  0 : 7056  1 : 053  0 : 4092] ^T with i = 1  2 for which we (not surprisingly) have computed

 ^ _p

 ^{^} _p

= 4 : 6074 10

!

2

This refers to the Matlab command \

"(4

0

2).

13

−20 −10 0 10 20 30 40

−0.05 0 0.05 0.1 0.15 0.2

Figure 5: Impulse responses of the pre-

lters:

(

) in (

) and

(

) in (

).

0 20 40 60 80 100 120 140

10⁻¹⁰ 10⁻⁵ 10⁰

L

Figure 6: Quantity

^

_pi

_r

for the time set

_p

(

):

= 1 (

) or 2 (

).

To get more illustrative results, we can vary the length L without changing the prelters.

For example, dene a reference model parameter estimate, denoted  r for either prelter as the vectors ^  _pi

become closer to one another for increasing L , while taking L very large. Then, for dierent values of L , we can evaluate the deviation between the model parameter estimates ^ 

_pi (according to such L ) and that  r , i.e.

 ^ _pi

 r

as a function of L . The result is depicted in Figure 6 for  r taken as the average between the two ^  _pi

for L = 500. This is a semilog-plot that reects the exponential decaying of the impulse responses of the prelter used. It is also interesting to note that the causal prelter, i.e.

M

( q ), generates less edge eects when L > 40. This is due to the fact that this prelter has a shorter impulse response (eective) length for a given tolerance.

To end with this part of the simulations, we can mention that if the two prelters do not belong to the same equivalent class then the deviation between the two model parameter estimates really depends upon the system noise by classical asymptotic formulas. For example, if we take M

( q ) = R ( q ) as a prelter then we obtain

 ^

_pi

 r

= 7 : 1101 10

for L = 220. Obviously, this result has nothing to do with the preceding deviations!

The second part of our simulations deals with the statistical moments of the model pa- rameter estimate ^  (

p ( d m L m )) while taking care of removing the edge eects that are present when evaluating the preltered prediction errors related to the original data set Z

^N .

In order to satisfy Assumption 1, let the reference prelter be M

( q ) = 1 that leads to an equivalent class of prelters identical to

( L 1) for a given L . Thus, the system dynam- ics belongs to the preltered (by M

( q )) model structure, i.e. the corresponding model parameter vector is 

. The result derived in Section 4 stipulates that if the prelters belong to that equivalent class, the statistical moments of the associated estimates are

14

0 5 10 15 20

−0.4

−0.2 0 0.2 0.4 0.6

Figure 7: Impulse response of the all-pass lter

(

).

asymptotically expressed by formula (10). To illustrate this result we will consider three causal prelters generated by a reference causal and all-pass lter denoted Q ( q ). These prelters are

M

( q ) = 1 M

( q ) = Q ( q ) M

( q ) = Q

( q ) while the lter Q ( q ) is chosen as

Q ( q ) =

0 : 512 + q

1

0 : 512 q

whose impulse response is presented in Figure 7. For a tolerance 10

, the (eective) length of this impulse response is L q = 105. The length of the dierent prelters is simply L

= 1, L

= L q and L

= 2 L q

1, respectively. These lengths will be used in the evaluation of the corresponding \edge eects"-free time sets, i.e.

T

pi = 1 + ( L i + L p

2)

( d

d p

) N

 ( d

d p

)] =  L i + 2 N ]

as d = 0 and L p = 3, over which the identication cost function V i ( Z ^N 

pi ) is dened.

Finally, the minimizing solution of this cost function is denoted ^  pi whose statistical moments are given in the expression (10), i.e. no bias with respect to 

and a variance that only depends on the number of elements in the corresponding time-domain set, i.e.

#(

pi ) = N

L i

1.

In order to compute these moments, we have performed 400 Monte-Carlo simulations of the data set Z

^N for a given input signal. For each experiment, we have computed the model parameter estimates ^  pi . The results are presented in Table 1 where we have shown the statistical moments of  pi = #(

pi ) =N ]

⁼

 ^ pi . Obviously, this table completely validates our derivations: no bias and variances inversely proportional to the number of

\edge eect"-free data in the corresponding time set, i.e. #(

pi ).

15

i 0 1 2

k

Bias(  pi )

3.3011 10

2.3403 10

4.9329 10

Var(  pi ) 2.9561 10

2.9230 10

2.9526 10

MSE(  pi ) 2.9594 10

2.9254 10

2.9576 10

Table 1: Statistical moments of the \normalized" model parame- ter estimates

pi = #(

pi )

]

⁼

^ _pi with respect to

.

6 Conclusions

In this note, we have address the issue of the use of prelters in prediction error identica- tion methods. We rst have derived equivalent classes of prelters whose characteristics are essentially to have identical transfer function moduli. That is to say that the phase of the prelter transfer functions does not play any particular role in the model estimates.

In other words, the prelter causality is not an issue for characterizing the equivalent class, i.e. given a noncausal prelter, one can always nd an equivalent causal prelter that will lead to the same parameter estimate for the identied model.

On the other hand, beyond an equivalent class, we have advocated for a best prelter which is the one whose impulse response has the shortest (eective) extend.

It thus appears that the problem of nding such a best prelter given a particular func- tion transfer modulus is an optimization problem that mixes frequency domain constrains with time domain length minimization.

References

1] M. Bellanger. Digital processing of signals. John Wiley and Sons, 1984.

2] L. Ljung. System Identication: theory for the user. Prentice Hall, 1987.

16