Copyright © IFAC Identification and System Parameter Estimation 1985, York , UK, 1985
P ARAMETER ESTIMATION AND IDENTIFICATION OF NON MINIMUM PHASE LINEAR STOCHASTIC SYSTEMS Chu Jia-jin and Wang Ying Department oJ Automation , T singhua Universil)" Beijing, The People'S Republic oJ China
Abstract.
In this paper, we introduce the problem of identification on
nonminimum phase linear stochastic systems (NPLSS), and discuss the problem about uniqueness that is whether we can decompose the output of li near stochastic system into convolution of two signals (one of them is the impulse response of the system, the other is a random signal) uniquely. Then the identifiability conditions of NPLSS are obtained and simulation result is given. Keywords. Identification; Parameter estimation; Stochastic system; Signal processing; Nonminimum phase system. ~(q-l) == 1 + '~' 1 ..." -1 + •.• + (-n, q- 0 -
INTRODUCTION A system is some nhy sical object, and its behavior can be described by some equations called its mathematical model. We usually indicate a system by a block for conveniec e , as shown in Fig. 1,
A(q-1) has no factors common wi th both B(q-l) and C(q-1) ( c ontrolable and observable). In practice, we often encounter t~e sys tem in which the r e is no control input u(t), such as, in signa l processing, In this c ase, the Eq(l .1) is changed into
zct)
u(t)
A(q-1)Y(t)
Fig.
Define
where y(t) is output of the system, u(t) is control inp~t, p(t) is unmeasurable process noise, z(t) is measurement value of output and n(t) is measurment noise.
A(z)
0<1
13(%1) =
1S VOL I-GG
•
.L: b
In this paper, we only discuss the system described by Eq(1.3). If A(z) and C(z) are both stabl e . the system is called minimum phase linear stoc has tic sys t em (MP LSS) . When some zeros of C(z) a re located outside the ~ nit circte in Z-plane, the syst em is called nonminimum phase line 8r stochastic system (NPLSS). At present, there are a
t) ( 1.1 )
z(t) =y(t) +n(t) where A(q-1 )
( 1.4)
where z is Z-transform operator. The Z-tran s form of a discrete t im e Si gnal b(O),b(l), is defined by
In this paper, we only discuss discrete linear time-invariant systems. In this case, the mathematical model of the system can be deseribed as A ( 0 -1 ) Y ( t) = B ( q -1 ) u ( t) + C ( q -1 ))(
(1. CS )
(1.2)
991
Chu Jia-jin and Wang Ying
992
lot of methods discussing estimation of MPLSS, but there are a few methods discussing identification of NPLSS. The problem of identification of NPLSS exists in practice, such as, in seismic signnl processing, in communication for blind adjustment of a linea r e!uazier, in astronomical data processing, as in Kormylo and Mendel (1983), Benveniste(1980)and Scargle (1981). In e80pnysical signal pro cess ing, when a impulse gener~ted by explosive transmits in subsurface, the phase of the impulse changes when it goes t ~rough different rocks. ~he acknowledge of phase is us eful to understand the property of rocks. Because the identification of phase is valuable for practice, more and more reserc~ers have begun to pay attention to NPLSS. lJefine H(z) = C(z)/A (z)
( 1.5)
OC>
I/.c.
UNIQUE:NESS ABOUT THE DECOMPOSITI ON OF' NPLSS In thi s section, we shall discuss the uniquehess aLout t ~e decomposition of NPLSS. we shall give two t heo rems. In this section, we use h(t) to express the impulse response of the system which is stable and causal and use i(t) to express its inV8rse satisfying h( t)*i( t) =~( t)
H(z) is a rational function, its inve rs e Z- tra nsfo r m h(t) is called wavelet in seismic exploration (i.e. t he impulse rssponse of the system described by Eq(1.3». USing convolution model, we express the system as
Y(t)=Zh(k~t-k)=b(t)
theorem in wold ('954). In NPLSS if we know nothing about the unmea s urable noise JIJ.. t) and system's impulse response h( t) , we can not estimate system's parameters uniquely, even though we have known that YJ- t) is w'li t e noise. So, it is nece ssary to discuss t he unioueness for parameter estimation in NPLSS. In the later section, we shall discuss t'lis problem.
*jt(t)
whe re n * n denots convolution operator. As we know, nonminimum phase system can be expressed as
C>Q
•
Z i "( t) c::: ""'" "t;=-...
Theorem I Assume t f)'1 t J.I-( t) 0 if t L. 0 and fit) = e(t) i f t~O, where eft) is w'lite noise of zero-average, ~(e(t)e(j»= Ca(t-j) (C is constant g reater tha n 0). If y(t) i s decomposed into
and
y(t)
h ,(t)*J(,(t)
(2.' )
y(t)
h ( t)* )1..2 (t) 2
(2.2)
whe re where Hmin(z) is the transfer f unction of minimum phase system, Hap (z) is the all pass filter. Assume tha t ~(t) is the inverse Z-transform of H(z) and hmin(t) is tne inverse Z-transform of Hmin(z), the relation
E(.At (t ~ (j»
a,S(t-j)
E ~( t )Jl:?( j»
a),(t-j)
E( JI.;>(t»
=0,
Jl. 2 (t) = 0, w'len C is known, a, = a 2 = C, Prom t'le above, we have conclusion
i f t.:. O,)J.{t) = 0 nnd
holds for all t ~O. Thus, the :nl"bu;n pha se system is also called minimum energy delay system. In MPLSS, we can estimate parameters in Eq(1.3) using the method of spectral analysis, aC80rding to Wold's d8composition
,.
h,(t) =
~h2(t)
when C is known
~
h (t) =th,(t)(a,/a 2 ) ~ when C is unknown 2 Proof: From Eq(2.'), we obtain
Nonminimum Phase Linear Stochastic Systems
i,(t)*y(t)
993
where
i,(t)*h2(t)*~(t)
a,b(t-j)
E (Jt, ( t)Jt, ( j ) )
Define E()l.2( t then
E(Jt,(t)) =0
=
-
E(){2(t»
;=0 r( t ,k)
,.
E(A (t)Jl, (t+k»
then r(t,k) = a 2
~
J:.
= 0,
)t,(t) and~2(t) are both i.i.d random sequence and of the same distribution as J(t). From the above, we obtain
c?:vM2( j ) f( t-j )
Define
f(t-j)f(t+k-j)
From the known condition, we have
so
)~( j) )
r( t,k) = a~k), t~ 0 ...... 2 . r(-' ,0) = a 2 Z f (-J) = 0
2.
h 2 (t) = ah,(t) when~(t) is non-Gaussian distribution h 2 (t) = a(g(t)*h,(t», when~(t) is Gaussian distribution.
Where g(t) is an all pass filter and
:J:f
f( j) = 0, j£O Let k = 0, t varies from zero to infinite large, we obtain f 2 (0) + f2( 1) f 2 (0) (a,/a 2 ) f 2 (0) + f2( 1) + f2(2) Hence
Proof: From theorem I, we obtain J(,(t) = f(t)*~(t) =
.2 JA
2(.il f( t-j)
j="'"
j>O
f( j) = 0,
therefore when k
!
0 00
a so
2
~ f(j)f(j+k) = 0
;=-00
when k = 0, When e is known, a, = a 2 = e, so Q.E.D
Define (f(t)/a)
From this t~eorem, we see that decomposition is unique whether the system is minimum phase or not, provided that ~(t) is one-sided w~ite noise. Theorem II
g(t)
then .,.... 2
s,,&g and when k
!
(j)
= ,
0
Zg(j)g(j+k) =
Assumej.l(t) is identically and independently distributed (i.i.d) random sequence of zero-average ()l(t) is not equal to zero, i f t.::O). Ify(t) is decomposed into
j=-
0
y(t)
h,(t)*fl,(t)
thus g(t) is an all pass filter. Define o(t) =JA.,(t)/a and q(O) =.Z)t2(j)g(-j)
y(t)
h 2 (t)*)S(t)
Let fq(x) and f~(x) express the characteristic function of q(t) and ~(t). Then,
-
3=-
994
Chu Jia-jin and Wang Ying
from the known conditions, we have fq(x) = f,}4(x) So, q(O) and)t2(t) are both of Gaussian distribution from the theoreM in Kagan, Linik and Rao (1973) if g( t) t- b ( t). In the first part of the theorem, we have known the~(t) is of non-Gaussian distribution, so, the above result contradicts with the known condition. It follows that g(t)=~(t). In the second part of the theorem, the result that q(t) is of Gaussian distribution does not contradict with the known condition, so get) is a all pass filter with ·zero phase shift or non-zero phase shift. Q.E.D
THE IDENTIFIABILITY AND IDENTIFICATION METHODS ABOUT NPLSS In this paper, assume that the structure of the system is already known, so we only discuss the problem on identification of sy s tem's parameters. There are several definitions on identifiability of parameters. The definitions which are often used are: Definition I If p(limg = eT) = 1, we say that the paN rameters eT is identifiable. Defini tion Il I N(9 T )) = 1, we say that If P(lim IN(i N) the pa rameters eT is identifiable. Definition III
e
N - OTt~f » = 0, where E is given small positive number arbitrarily, we say that the parameters eT is identifiable.
I f lim(p(\
Definition IV I f lim(P( IJN(gN) -
IN(Or)l~f))
= 0, we say
that parameters8T is identifiable.
eN'
In the above, 6TE@(® is parameters space of Rn ), eT is true parameters, " is estimates for Q T' I N(9) is an object function, N is the number of sample data,
eN
n is the number of parameters. Which definitions can be used depends on particular problem. We do not want to discuss this problem further. From theorem 11, we see that the NPLSS is not identifiable only using the second order moment, so we need higher order moment to identify the system. Rosenblatt(1980) and Lii (1982) have proved that het) is identifiable using higher order moment, ~(t) being i.i.d non-Gaussian white noise. The advantange of this method is that we need not know the probability distribution ofjJ..(t).
Fig. 2 When we know~(t) is i.i.d non-Gaussian white noise, e(t) (as shown in Fig. ~) is i.i.d white noise with the same distributi on as)t(t), if and only if "' (t)*i(t) = cb'(t), where C is constant. It follows that we can estimate het) by adjusting i(t) to make e(t) i.i.d white noise with the same distribution asY-(t), where the system which is ch ~ racterized by i(t) is noncausal stable system. When we have known the distribution of we can use maximum-likelihood method to estimate het). This method can usually be converted into the problem in which we make some functional rea c h maximum or minimum value. The methods of Benveniste (1980) and Kormylo (1979) are particular examples for this method. Which functional can be used depends on the distribution of }let). Benveniste (1980) has proved that het) is identifiable using his method when JJ.i.. t) is sub or super-Ganssian distribution. Kormylo (1979) has also proved t hat het) is identifiabl e using his method when .At) is Gaussian-Bernoulli distribution. The disadvantage of this method is that we must know the probability distribution of}t(t).
~(t),
From theorem I, we sided white noise, h(t)*i(t) = CS(t), white noise, where this principle, we
see that e(t) is oneif and only if when~(t) is one-sided C is constant. From can establish a func-
Nonminimum Phase Linear Stochastic Systems
tional J k (8)
=
le
~
j~
E{e(t)e{t+j»
In practice. tie ensemble average about e{t) can not be obtained. only the sample data of y{t) can be used. so we need use the sample average instead of the ensemble average. BecausevM(t) is a nonstationary random seouence defined in (-OO.bO) we must use multiple sample paths Y.i(t) (j .2 ..... M) to estimate autocov~riance of €(t). Let D(i) E(e(t)e(t + i» Then
995
Ci( t) is Gaussian distribution. fJ{ t) is Bernoulli sequence. We aSSume that this sequence is known. We use Kormylo's method to estimate h{t). The simulation result is given in Fig. '5.
0-+
=,
DU)
=
(Z ej(t)ej(t+i»!M j=o
-.·t ·Iot
(4.1)
-D-'
A new functional k " ..... ,.,
J~(e) =( Z
ZZ
\·~I;=I
F"3·1
e j ( t)e j( t+i»/M/N
(4.2)
...
is feasible in practice. The inverse of i(t) which makes the above functional (Eq(4.2» reach minimum value is the estimate of h(t). We are researching this method. the result will present in the future paper.
,. ,
o'f
,·t
Because the number of parameters in h(t) is large. we use the ARMA model to describe the system.
Fig. '3
In Fig. '3 G(q-' ) F(q-')
fo + f,q-' + •.. + fnq-n
We adjust parameters (g, instead of adjusting i(t).
~
fn ... fo)
SIMULATION RESULT We use the waveform as shown in Fig. 4 as impulse response of system and u~e pulse sequence as shown in Fig. 6 as process noise to generat the output of system as shown in Fig. 7. The pro"ess noise.)J..( t) is Gaussian-Bernoulli distribution, the measurement noise n(t) is Gaussian distribution. )t{t) = a(t)~{t)
I f ,
..
I'
I
I
I~
~ I11
Ill,
11
I
11111, III
ri,
I
I
JI
11)1
"
R;. ,
-I~~~~ F~.
7
996
Chu Jia-jin and Wang Ying
CONCLUSION According to the above, we see that NPLSS is identifiable in some special ca~es. From the papers we have seen and the work we have done, the following results are concluded. In the condition that there is nothing known about the impulse response of the system h(t), the system should be identifiable in the condition (1))(( t) is stationary random sequence of non-Gaussian distribution (2)~(t~ is nonstationary in order two random seouence of Gaussian distribution (3)~t) is nonstationary in order two random seauence of non-Gaussian distribution.
In the above, ~(t) is defined in (-~,~). In this field, there is a lot of work to do. This is a significant subject.
REFERENCES Benveniste, A., M. Goursat and G. Ruget (1980). Robust identification of nonminimum phase system: blind adjustment of a linear equalizer in data communication. IEEE Trans. on Automatic Control, Vol. AC-25, No. '3, pp. 285-298. Kagan, A.M., YU.V. Linik and C.R. Rao (1973). Characterization Problems in Mathematical Statistics. Chapter five. Jo~n Wiley & Sons, Inc .• Kormylo, J. (1979). Maximum-likelihood seismiC deconvolution. Ph. D. dissertation, University of Southern California, Los Angeles, Ca .. Kormylo, J. and J. Mendel (1983). Maximum-likelihood seismic deconvolution. IEEE Trans. on Geosci. Electron., Vol. GE-21 , No. 1, pp. 72-81. Lii, K.S. and M. Rosenblatt (1982). ~ volution and estimation of transfer function phase and coefficients for non-Gaussian linear process. The Annals of Statistics, Vol. 10, No. 4, pp. 1195-1208. Rosenblatt, M. (1980). Linear processes and bispectra. J. Appl. Prob. Vol. 17, pp. 265-270.
Scargle, J.D. (1981). Phase-sensitive deconvolution to model random process?~ with special reference to astronomical data. Applied Time Series AnalYSiS, edited by D.F. Findley, Academic Press pp. 549-564. Wold, H. (1954). A Study in the Analysis of Stationary Time Series. 2nd ed. Upplsala, Sweden.