A general smoothing equation for signal estimation ... - Semantic Scholar
Digital Signal Processing 16 (2006) 369–388 www.elsevier.com/locate/dsp
A general smoothing equation for signal estimation using randomly delayed observations in the correlated signal-noise case S. Nakamori a,∗ , A. Hermoso-Carazo b , J. Linares-Pérez b a Department of Technology, Faculty of Education, Kagoshima University, 1-20-6, Kohrimoto,
Kagoshima 890-0065, Japan b Departamento de Estadística e Investigación Operativa, Universidad de Granada, Campus Fuentenueva, s/n,
18071 Granada, Spain Available online 16 June 2005
Abstract This paper treats the least-squares linear smoothing problem for signal estimation using measurements contaminated by additive white noise correlated with the signal, with stochastic delays. We derive a general smoothing equation which is applied to obtain specific smoothing algorithms, which are referred in the signal estimation literature as fixed-point, fixed-interval, and fixed-lag smoothing. Using an innovation approach, the general smoothing equation is derived without requiring the whole knowledge of the state-space model generating the signal, but only covariance information of the signal and the observation noise, as well as the delay probabilities. © 2005 Elsevier Inc. All rights reserved. Keywords: Fixed-point smoothing; Fixed-interval smoothing; Fixed-lag smoothing; Covariance information; Delayed observations
* Corresponding author.
E-mail addresses: [email protected] (S. Nakamori), [email protected] (A. Hermoso-Carazo), [email protected] (J. Linares-Pérez). 1051-2004/$ – see front matter © 2005 Elsevier Inc. All rights reserved. doi:10.1016/j.dsp.2005.04.013
370
S. Nakamori et al. / Digital Signal Processing 16 (2006) 369–388
1. Introduction For a long time, least-squares estimation problems of random signals from observed data have been the subject of intensive research. There are several solution methods, which depend on the process assumptions and the models that represent the possible relationships between the unknown signal and the observable variables. If the observations of the signal are without stopping in time, there are many algorithms capable of obtaining the signal estimates (see, for example, [1–4] in an uncertain measurement context). However, many real-world problems exist in which errors due to transmission delays may be present and, hence, the measurements do not arrive in time but are randomly delayed. Such situations arise frequently in engineering applications and in certain economic studies (for example, the delay between submitting a patent application and actually obtaining the patent). Systems with random delays are considered in communication theory (Nilson et al. [5], discuss modelling and analysis of systems subject to random time-delays in the communication network and they present a method for different control schemes) and also in some control applications (Kolmanovsky and Maizenberg [6] consider a finite horizon optimal control problem for a class of linear systems with a randomly varying time-delay, affecting only the state variable, which is modelled by a Markov process with a finite number of states). In the cases in which measurements are randomly delayed, it is necessary to modify of the standard algorithms for the signal estimation problem in order to incorporate the effects of the random delay. Using different approaches, the state estimation problem for system models with randomly varying time-delays has been investigated and modifications of conventional signal estimation algorithms have been proposed to accommodate the randomly delayed observations. Let us cite, among other papers, the following: Evans and Krishnamurthy [7] consider state estimation for a discrete-time hidden Markov model in which the observations are delayed by a random time and the delay process is modelled as a finite state Markov chain. In Yaz et al. [8], state estimation in linear discrete-time state-space models with stochastic parameters is treated using linear matrix inequalities, and the results are applied to the problem of state estimation with a random sensor delay. Su and Lu [9] designed an extended Kalman filtering algorithm which provides optimal estimates of interconnected network states for systems in which some or all measurements are delayed. More recently, Matveev and Savkin [10] proposed a recursive minimum variance state estimator in a linear discrete-time partially-observed system perturbed by white noises when the observations are transmitted via communication channels with random transmission times and when various measurement signals may incur independent delays. In this paper, we consider the situation in which the state-space model of the signal is not available and the signal observations may be delayed by one sampling time; hence, the observed data may not contain the actual signal (randomly delayed observations). The delay is modelled by a sequence of independent Bernoulli random variables whose values, one or zero, indicate whether the measurements arrive on time or are delayed by one sampling time, respectively. Moreover, since the assumption that the signal process and the observation noise are independent can be restrictive in many important engineering and stochastic control problems [11–13], we consider the situation in which randomly delayed observations are perturbed by a white noise which is correlated with the signal process.
S. Nakamori et al. / Digital Signal Processing 16 (2006) 369–388
371
In Nakamori et al. [14], a recursive algorithm for the least-squares linear one-stage prediction and filtering problems was obtained under the above hypotheses concerning the model; in this paper we are interested in complementing the latter study by considering linear interpolation or smoothing problems, and derive recursive algorithms satisfied by smoothed estimates. A general smoothing equation is deduced from an innovation approach, which simplifies obtention since the innovation process is a white noise; from this equation we obtain the fixed-point, fixed-interval, and fixed-lag smoothing algorithms. The proposed algorithms do not require the whole knowledge of the state-space model generating the signal, but only covariance information of the signal and the observation noise, as well as the delay probabilities. Specifically, it is assumed that both, the autocovariance functions of the signal and the crosscovariance function between the signal and the observation noise are expressed in a semi-degenerate kernel form. The precision of the least-squares estimators is measured by the estimation error covariance matrices, which provide a global measurement of the performance of the estimators. For this reason, the formulas to obtain such matrices (for the fixed-point, fixed-interval, and fixed-lag smoothers) are also included in the present study.
2. Problem formulation We consider observations of a n × 1 signal zk described by the equation y˜k = zk + vk ,
k 0,
where {vk ; k 0} is the observation noise. It is assumed that the real observation at time k can be delayed by one sampling period, with a certain probability pk , or updated, with a probability 1 − pk . Therefore, if {γk ; k 1} denotes a sequence of independent Bernoulli random variables with P [γk = 1] = pk , the received measurement at time k is given by yk = (1 − γk )y˜k + γk y˜k−1 ,
k 1.
(1)
Remark. If γk = 1, then yk = y˜k−1 and the measurement is delayed by one sampling period. Otherwise, if γk = 0, then yk = y˜k , which means that the measurement is up-todate. Therefore, the value pk represents the probability of a delay in the measurement yk . In applications of communication networks, the noise {γk ; k > 1} usually represents the random delay from sensor to controller and the assumption of one-step sensor delay is based on the reasonable supposition that the induced data latency from the sensor to the controller is restricted so as not to exceed the sampling period [14]. To treat the least-squares (LS) estimation problem of the signal from the observations given by Eq. (1), we assume the following hypotheses on the signal and noise processes: (i) The signal process {zk ; k 0} has zero mean and its autocovariance function is expressed in a semi-degenerate kernel form as T Ak BsT , 0 s k, Kz (k, s) = E zk zs = Bk ATs , 0 k s,
372
S. Nakamori et al. / Digital Signal Processing 16 (2006) 369–388
where Ak and Bk are known n × M matrix functions. (ii) The noise process {vk ; k 0} is a zero-mean white sequence with known autocovariance function 0, if s = k, E vk vsT = Rk , if s = k. (iii) The signal process {zk ; k 0} and the noise process {vk ; k 0} are correlated, being T Kzv (k, s) = E zk vsT = Ck Ds , 0 s k, 0, 0 k < s, where Ck and Dk are known n × N matrix functions. (iv) The multiplicative noise {γk ; k 1} is a sequence of independent Bernoulli random variables with known probabilities, P [γk = 1] = pk > 0, ∀k 1. (v) The process {γk ; k 1} is independent of {zk ; k 0} and {vk ; k 0}. Remark. The special form of correlation in hypothesis (iii) allows us to consider the estimation problem of signals from randomly delayed observations when the evolution of the signal is influenced by the observation noise up to the considered time. This fact is shown in the numerical example (Section 5), where the disturbances of the signal and the observation are correlated at consecutive time instants and, consequently, the signal at time k is correlated with the observation noise up to the time k. Our aim is to approach the LS linear estimation problem of the signal zk , based on the observations {y1 , . . . , yL } with L > k; that is, the time k at which it is desired to estimate the signal is less than the time L of the last measurement. More specifically, we are interested in deriving of recursive smoothing algorithms which can be applied to various practical situations; specifically, we seek to obtain the smoothers of the signal at a fixed point, k, when the number of available measurements, L, is increasing (fixed-point smoother), the smoothers of the signal at different instants, k, based on a fixed number, L, of measurements (fixed-interval smoother) and the smoothers of the signal at each time, k, based on the observations up to the time k + λ, where λ is a fixed number (fixed-lag smoother). In order to simplify obtaining the algorithms, we use an innovation approach; the innovation at time i represents the new information which provides the observation at time i and is given by νi = yi − yˆi,i−1 , where yˆi,i−1 denotes the LS linear estimator of the observation yi based on the observations {y1 , . . . , yi−1 }. It is known (Kailath [15]) that the LS linear estimator of zk based on the observations {y1 , . . . , yj }, which is denoted by zˆ k,j , is equal to the LS linear estimator based on the innovations {ν1 , . . . , νj }. Hence, the fact that the innovations constitute a white process leads us to the following expression for the estimator of the signal (Nakamori et al. [4]), zˆ k,j =
j Sk,i Πi−1 νi , i=1
where
Sk,i = E zk νiT
and Πi = E νi νiT .
(2)
S. Nakamori et al. / Digital Signal Processing 16 (2006) 369–388
373
Taking into account expression (2), it is possible to obtain the smoothers, zˆ k,L , L k + 1, from the filter, by adding the remaining terms of the sum in (2), whose coefficients, Sk,i , i k + 1, usually have a common form (Nakamori et al. [4]); however, for the model at hand, the correlation between the signal and the noise specified in hypothesis (iii), together with the delay in the observations, leads to different expressions for the coefficients Sk,k+1 and Sk,i , i k + 2. Hence, in this case, it is more appropriate to start from the single-stage smoother, zˆ k,k+1 , and express the smoothers, zˆ k,L for L > k + 1 as zˆ k,L = zˆ k,k+1 +
L
Sk,i Πi−1 νi ,
L > k + 1.
(3)
i=k+2
Therefore, we first derive an algorithm to obtain zˆ k,k+1 and, afterwards, the different smoothing algorithms (fixed-point, fixed-interval, and fixed-lag) are derived from a general formula for the second term of the right member of (3).
3. Single-stage smoother To obtain the single-stage smoother at time k, we start from (2) for j = k + 1 and, by separating the last term of the sum, it is immediately obtained that that estimator verifies −1 νk+1 . zˆ k,k+1 = zˆ k,k + Sk,k+1 Πk+1
Then, since the filter, zˆ k,k , is required as an input for each k, the single-stage smoother is performed accompanied by the former; hence, we first present the recursive filtering algorithm established in Nakamori et al. [14]. Since some expressions of the proof of the filtering algorithm are used in the deduction of the smoothers, we present a brief proof of this algorithm, in Appendix A. Theorem 1. Consider the observation equation (1) verifying hypotheses (i)–(v); then the signal filter, zˆ k,k , is given by 1 zˆ k,k = (Ak | Ck )Ok ,
k 1,
(4)
where the vector Ok is recursively calculated from Ok = Ok−1 + Jk Πk−1 νk ,
k 1;
O0 = 0,
(5)
where νk is the innovation, which satisfies νk = yk − G(Ak |Ck ) Ok−1 − Hk νk−1 ,
k 1;
ν0 = 0,
(6)
and T T Jk = G(B − rk−1 G(A − Jk−1 HkT , k |Dk ) k |Ck )
k 1;
J0 = 0,
(7)
where the matrices GYk (Yk = (Ak | Ck ), (Bk | Dk )) and Hk are given by −1 T GYk = (1 − pk )Yk + pk Yk−1 , Hk = pk (1 − pk−1 ) Dk−1 Ck−1 + Rk−1 Πk−1 . 1 (U | V ) denotes a partitioned matrix into sub-matrices U and V .
374
S. Nakamori et al. / Digital Signal Processing 16 (2006) 369–388
The matrix rk = E[Ok OkT ] is recursively calculated from rk = rk−1 + Jk Πk−1 JkT ,
k 1;
r0 = 0,
(8)
and Πk , the covariance matrix of the innovation νk , is given by Π = (1 − pk ) Ak BkT + Ck DkT + Dk CkT + Rk T T T + Ck−1 Dk−1 + Dk−1 Ck−1 + Rk−1 + pk Ak−1 Bk−1 T T T − G(Ak |Ck ) rk−1 G(A +Jk−1 HkT − Hk Jk−1 G(A + Πk−1 HkT , k |Ck ) k |Ck )
k 1;
Π0 = 0.
(9)
Proof. The proof of Theorem 1 is given in Appendix A.
2
In the following theorem, we present an algorithm for computing zˆ k,k+1 , the singlestage smoothed estimator of the signal zk , for any given k 1. Theorem 2. Consider the observation equation (1), satisfying hypotheses (i)–(v). Then, the single-stage smoother of the signal zk can be obtained from −1 νk+1 , zˆ k,k+1 = zˆ k,k + Sk,k+1 Πk+1
k 1,
(10)
where 2 T Sk,k+1 = (Bk | 0n×N ) − (Ak | Ck )rk G(A k+1 |Ck+1 ) T − (Ak | Ck )Jk Hk+1 + pk+1 Ck DkT .
(11)
The filter zˆ k,k , the matrices G(Ak |Ck ) , Hk , rk and Jk , the innovations νk and their covariances Πk are obtained from the formulas given in Theorem 1. Proof. As indicated above, Eq. (10) for the single-stage smoother is immediately obtained from (2) for j = k + 1. T ], we use formula (6) for the In order to calculate expression (11) for Sk,k+1 = E[zk νk+1 innovation νk+1 and we obtain T T T − E zk νkT Hk+1 . − E zk OkT G(A Sk,k+1 = E zk yk+1 k+1 |Ck+1 ) T ] and expression (A.5) for O , it is Then, using the hypotheses on the model for E[zk yk+1 k deduced that T Sk,k+1 = (Bk | 0n×N )G(A + pk+1 Ck DkT k+1 |Ck+1 )
−
k T T Sk,j Πj−1 JjT G(A − Sk,k Hk+1 . k+1 |Ck+1 ) j =1
Now, using (A.3) for Sk,j , 1 j k and afterwards (A.7), the formula (11) is proven. 2 0 a×b denotes the a × b zero matrix.
2
S. Nakamori et al. / Digital Signal Processing 16 (2006) 369–388
375
Error covariance matrices. The precision of the filter and the single-stage smoother is measured by the covariance matrices of the corresponding estimation errors, T Σk,j = E zk − zˆ k,j zk − zˆ k,j , j = k, k + 1. Since the error zk − zˆ k,j is orthogonal to the estimator zˆ k,j , it is immediately verified that T Σk,j = Kz (k, k) − E zˆ k,j zˆ k,j . Then, hypothesis (i) on the model and expressions (4) and (10) for the filter and the singlestage smoother, respectively, lead to Σk,k = Ak BkT − (Ak | Ck )rk (Ak | Ck )T , −1 T Σk,k+1 = Σk,k − Sk,k+1 Πk+1 Sk,k+1 ,
k 1,
(12)
k 1.
(13)
4. General expression for the linear smoother Once the algorithm to obtain the single-stage smoother, zˆ k,k+1 , has been deduced, our purpose is to derive a general expression for the remaining smoothers, zˆ k,L , L > k + 1. For
−1 this purpose, as expressed in (3), we must determine the sum L i=k+2 Sk,i Πi νi or, more T specifically, the coefficients Sk,i = [zk νi ], for i k + 2 which, from expression (6) for νi , are given by T T T T G(Ai |Ci ) − E zk νi−1 Hi . Sk,i = E zk yiT − E zk Oi−1 Taking into account the hypotheses on the model and (A.5), it is obtained that T − Sk,i = (Bk | 0n×N )G(A i |Ci )
i−1 T Sk,j Πj−1 JjT G(A − Sk,i−1 HiT , i |Ci )
i k + 2.
j =1
Then, using (A.3) and considering (A.7) for rk , these coefficients are expressed as follows: −1 T T Jk+1 G(Ai |Ci ) Sk,i = (Bk | 0n×N ) − (Ak | Ck )rk − Sk,k+1 Πk+1 −
i−1 j =k+2
T Sk,j Πj−1 JjT G(A − Sk,i−1 HiT , i |Ci )
i > k + 2;
−1 T T T Jk+1 G(Ak+2 |Ck+2 ) − Sk,k+1 Hk+2 . Sk,k+2 = (Bk | 0n×N ) − (Ak | Ck )rk − Sk,k+1 Πk+1 This latter expression for Sk,i allows us to deduce that Sk,i = Mk k,i ,
i k + 2,
(14)
where
−1 T Jk+1 − Sk,k+1 Mk = (Bk | 0n×N ) − (Ak | Ck )rk − Sk,k+1 Πk+1
and k,i verifying
(15)
376
S. Nakamori et al. / Digital Signal Processing 16 (2006) 369–388
k,i = (G(Ai |Ci ) | 0n×n )T −
i−1 j =k+2
T k,j Πj−1 JjT G(A − k,i−1 HiT , i |Ci )
i > k + 2;
T k,k+2 = G(Ak+2 |Ck+2 ) | Hk+2 .
(16)
Hence, if we define qk,L =
L
k,i Πi−1 νi ,
L > k + 1;
qk,k+1 = 0,
(17)
i=k+2
the following general expression for the smoother is obtained from (3), using (14), zˆ k,L = zˆ k,k+1 + Mk qk,L ,
L k + 1.
(18)
General expression for the smoother error covariance matrices. Taking into account that the error covariance matrices of the smoothers can be expressed as Σk,L = Kz (k, k) − T ], the general expression (18) for the estimators leads to the following general E[ˆzk,L zˆ k,L expression for such matrices Σk,L = Σk,k+1 − Mk Qk,L MTk ,
L k + 1,
(19)
where Σk,k+1 , the single-stage smoother error covariance matrices, are given in (13) and T ]. Qk,L = E[qk,L qk,L Next, different smoothing algorithms are obtained from the general expression (18); specifically, the smoothing estimators of the signal are recursively obtained in three different ways: • Fixed-point smoothing; a fixed integer k 1 is considered and the estimators zˆ k,L with L = k + 1, k + 2, . . . , are obtained (recursions for increasing L are proposed in Section 4.1). • Fixed-interval smoothing; the estimators zˆ k,L with L a fixed positive integer and k = L − 1, L − 2, . . . , 1, are derived (recursions in k are established in Section 4.2). • Fixed-lag smoothing; here, the obtained estimators are zˆ k,k+λ , where λ is a positive integer and k = 1, 2, . . . (recursions for increasing k are obtained in Section 4.3). Recursive relationships for the error covariance matrices of the smoothers in the three above cases are obtained in the respective sections mentioned above. 4.1. Fixed-point smoothing algorithm In this section we wish to determine a recursive algorithm to implement the computation of zˆ k,L for a fixed k 1, and for any L k + 1; it is not necessary to know, a priori, the time point of the final observation, since the algorithm proceeds forward recursively in time and can be stopped at will. So, taking into account the general expression (18), to solve this problem we need to obtain a forward recursive formula in L for qk,L . From definition (17) it is immediately deduced that qk,L = qk,L−1 + k,L ΠL−1 νL ,
L > k + 1;
qk,k+1 = 0.
(20)
S. Nakamori et al. / Digital Signal Processing 16 (2006) 369–388
377
Now, if we denote δk,i =
i
k,j Πj−1 JjT ,
i > k + 1;
δk,k+1 = 0,
(21)
j =k+2
expression (16) is rewritten as 3 T k,L = (IM+N | 0(M+N )×n )T − δk,L−1 G(A − k,L−1 HLT , L |CL ) k,k+2 = (G(Ak+2 |Ck+2 ) | Hk+2 ) , T
L > k + 2; (22)
and, from (21), it is clear that the function δ verifies δk,L = δk,L−1 + k,L ΠL−1 JLT ,
L > k + 1;
δk,k+1 = 0.
(23)
So, Eqs. (20), (22) and (23) provide the desired forward recursive formula in L for qk,L . Fixed-point smoother error covariance matrices. The error covariance matrices of the smoother, Σk,L , are given in (19). Now, we need to obtain a forward recursive expression T ]. in L for Qk,L = E[qk,L qk,L Since qk,L−1 and νL are uncorrelated, the following expression is immediately obtained from (20), Qk,L = Qk,L−1 + k,L ΠL−1 Tk,L ,
L > k + 1;
Qk,k+1 = 0.
(24)
Remark. From (2), we have the following alternate form of the fixed-point smoothing algorithm zˆ k,L = zˆ k,L−1 + Sk,L ΠL−1 νL ,
L > k + 1,
(25)
whose initial condition is the single-stage smoother, zˆ k,k+1 , given in Theorem 2. The matrices Sk,L can be calculated from (14); that is, Sk,L = Mk k,L ,
L > k + 1,
where Mk is given in (15) and k,L can be recursively calculated from expressions (22) and (23). Finally, from (25), we also have the following alternate form for obtaining the fixedpoint smoothing error covariance matrices, T Σk,L = Σk,L−1 − Sk,L ΠL−1 Sk,L ,
L > k + 1,
(26)
with initial condition Σk,k+1 given in (13). 4.2. Fixed-interval smoothing algorithm In this section we assume that the measurement data over the fixed interval [1, L] are available and that, for each time point k within the interval, we wish to obtain the LS linear estimator of the signal zk based on all the observations {y1 , . . . , yL }; thus, since 3 I denotes the a × a identity matrix. a
378
S. Nakamori et al. / Digital Signal Processing 16 (2006) 369–388
the estimators of zL and zL−1 are the filter, zˆ L,L , and the single-stage smoother, zˆ L−1,L , respectively, our purpose is to obtain the remaining estimators, zˆ k,L , for k < L − 1. For this purpose, we start from the general expression (18) and the fixed-interval smoothing algorithm, whose starting point is zˆ L−1,L , is obtained by deducing a backward recursive formula in k (k < L − 1) for qk,L . From (16) and (17), it is immediate that −1 qk,L = (G(Ak+2 |Ck+2 ) | Hk+2 )T Πk+2 νk+2 +
L
k,i Πi−1 νi ,
k < L − 2;
i=k+3
qL−2,L = (G(AL |CL ) | HL )T ΠL−1 νL . (1)T
Now, we express k,i = (k,i (a)
(27)
(2)T
| k,i )T and calculate, using (16), the differences
(a)
k,i − k+1,i (a = 1, 2), for i k + 3; by comparing the resultant expressions with (a)
those obtained from (16) for k+1,i , we deduce the following relations for i k + 3, (1) (1) (2) T T k,i = IM+N − G(A Π −1 J T k+1,i − G(A , k+2 |Ck+2 ) k+2 k+2 k+2 |Ck+2 ) k+1,i −1 T T T k,i = −Hk+2 Πk+2 Jk+2 k+1,i − Hk+2 k+1,i . (2)
(1)
(2)
Hence, we have k,i = Φk+2 k+1,i , where Φk+2 =
i k + 3,
(28)
T Π −1 J T IM+N − G(A k+2 |Ck+2 ) k+2 k+2
T −G(A k+2 |Ck+2 )
T Π −1 J T −Hk+2 k+2 k+2
T −Hk+2
,
and by substituting (28) in (27), we obtain the desired backward recursive formula for the vectors qk,L , −1 qk,L = Φk+2 qk+1,L + (G(Ak+2 |Ck+2 ) | Hk+2 )T Πk+2 νk+2 ,
qL−2,L = (G(AL |CL ) | HL )
T
ΠL−1 νL .
k < L − 2; (29)
Fixed-interval smoother error covariance matrices. The formula for the error covariance matrices of the fixed-interval smoother, Σk,L , is obtained from (19), using the following T ], which is obtained from (29) taking into account recursive relation for Qk,L = E[qk,L qk,L that qk+1,L and νk+2 are uncorrelated, T Qk,L = Φk+2 Qk+1,L Φk+2 + (G(Ak+2 |Ck+2 ) | Hk+2 )T −1 × Πk+2 (G(Ak+2 |Ck+2 ) | Hk+2 ),
QL−2,L = (G(AL |CL ) | HL )
T
k < L − 2;
ΠL−1 (G(AL |CL ) | HL ).
(30)
4.3. Fixed-lag smoothing algorithm We wish now to obtain an algorithm for the smoothers zˆ k,k+λ , where k = 1, 2, . . . , and λ 1 is a positive integer which represents a fixed-lag of λ units of time.
S. Nakamori et al. / Digital Signal Processing 16 (2006) 369–388
379
Note that, in the case λ = 1, the request smoothers, zˆ k,k+1 , can be obtained from the filtering algorithm as is proposed in Theorem 2. On the other hand, if we consider λ = 2, in view of general expression (18), using definition (17) for qk,k+2 and (16) for k,k+2 , we have −1 νk+2 , zˆ k,k+2 = zˆ k,k+1 + Mk (G(Ak+2 |Ck+2 ) | Hk+2 )T Πk+2
k 1.
Hence, we must only consider the cases λ 3. Like the fixed-point and fixed-interval smoothing algorithms, the fixed-lag one is again obtained from the general expression (18); hence, the smoothers zˆ k,k+λ are calculated from the following relation: zˆ k,k+λ = zˆ k,k+1 + Mk qk,k+λ ,
k 1,
(31)
using a recursive expression in k for qk,k+λ . For k 2, we can use definition (17) to calculate and compare the vectors qk,k+λ and qk−1,k−1+λ ; then, by considering the expressions (16) for k−1,k+1 and (28) to deduce that k−1,i = Φk+1 k,i , i k + 2, we obtain the following relation between the above two vectors, −1 −1 νk+1 qk−1,k−1+λ − (G(Ak+1 |Ck+1 ) | Hk+1 )T Πk+1 qk,k+λ = Φk+1 −1 + k−1,k+λ Πk+λ νk+λ , k 2. Thus, since −1 Φk+1 =
IM+N
T −G(A H −T k+1 |Ck+1 ) k+1
−1 T −Πk+1 Jk+1
−1 T −T T −[In − Πk+1 Jk+1 G(A ]Hk+1 k+1 |Ck+1 )
and, hence, −1 (G(Ak+1 |Ck+1 ) | Hk+1 )T = (0n×(M+N ) | −In )T , Φk+1
the above equation for qk,k+λ can be rewritten as −1 −1 qk−1,k−1+λ + (0n×(M+N ) | In )T Πk+1 νk+1 qk,k+λ = Φk+1 −1 + k,k+λ Πk+λ νk+λ ,
k 2.
(32)
In order to initiate this recursive relation, the vector q1,1+λ must be introduced into Eq. (32) at time k = 2; this vector can be obtained from expression (20), starting with q1,2 = 0 and processing the observations to obtain q1,3 , q1,4 , . . . , q1,1+λ . In view of (32), the fixed-lag smoothing algorithm also requires a recursive expression in k for k,k+λ , which is obtained using (21), (22), and (28), −1 −1 T δk−1,k−1+λ − (0n×(M+N ) | In )T Πk+1 Jk+1 k,k+λ = (IM+N | 0(M+N )×n )T − Φk+1 −1 T T × G(A − Φk+1 k−1,k−1+λ Hk+λ , k+λ |Ck+λ )
k 2.
(33)
The start point for (33), 1,1+λ , is obtained, starting with 1,3 = (G(A3 |C3 ) | H3 )T , from expressions (22) and (23), to obtain 1,4 , 1,5 , . . . , 1,1+λ , and δ1,1+λ is obtained from (23), starting with δ1,2 = 0. Finally, taking into account (33), it is also necessary to determine a recursive formula in k for δk,k+λ , k 1. From (17) and (A.5), since the innovation is a white process, it is clear
380
S. Nakamori et al. / Digital Signal Processing 16 (2006) 369–388
T ]. Then, by using (32) that δk,k+λ , as defined in (21), verifies that δk,k+λ = E[qk,k+λ Ok+λ for qk,k+λ , (5) for Ok+λ and taking into account that νk+λ and qk−1,k−1+λ are uncorrelated, we deduce the following relation: −1 −1 T δk−1,k−1+λ + (0n×(M+N ) | In )T Πk+1 Jk+1 δk,k+λ = Φk+1 −1 T + k,k+λ Πk+λ Jk+λ ,
k 2,
(34)
where δ1,1+λ is obtained from (23), starting with δ1,2 = 0. Equations (32)–(34) provide the desired recursive expression in k for qk,k+λ . Fixed-lag smoother error covariance matrices. From (19), the formula for the error covariance matrices of the fixed-lag smoother, zˆ k,k+λ , is Σk,k+λ = Σk,k+1 − Mk Qk,k+λ MTk ,
k 1,
(35)
T with Qk,k+λ = E[qk,k+λ qk,k+λ ]. It is clear that, for λ = 1, Qk,k+1 = 0 and, for λ = 2, −1 Qk,k+2 = (G(Ak+2 |Ck+2 ) | Hk+2 )T Πk+2 (G(Ak+2 |Ck+2 ) | Hk+2 ).
The recursive formula to obtain Qk,k+λ for λ 3 is derived from (32), taking into acT ]= T count that E[qk−1,k−1+λ νk+1 k−1,k+1 = (G(Ak+1 |Ck+1 ) | Hk+1 ) and E[qk−1,k−1+λ × T νk+λ ] = 0; so, −1 −T −1 Qk,k+λ = Φk+1 Qk−1,k−1+λ Φk+1 − (0n×(M+N ) | In )T Πk+1 (0n×(M+N ) | In ) −1 Tk,k+λ , + k,k+λ Πk+λ
k 2,
(36)
and the initial condition, Q1,1+λ , is obtained from expression (24), starting with Q1,2 = 0.
5. Numerical example In this section, we present a numerical simulation example to illustrate the application of the different smoothing algorithms proposed in the paper. For this purpose, we have performed a program in MatLab which simulates, at each iteration, the signal and the observation values and provides the smoothing estimates as well as the smoother error covariance matrices. We have considered the same example as in Nakamori et al. [14]; specifically, the signal {zk ; k 0} is a scalar process generated by the following first-order autoregressive model, zk+1 = 0.95zk + wk ,
(37)
where {wk ; k 0} is a zero-mean white Gaussian noise with Var[wk ] = 0.1, for all k. The covariance function of the signal is given by Kz (k, s) = 1.025641 × 0.95k−s ,
0 s k,
which can be expressed in a semi-degenerate kernel form, with Ak = 1.025641 × 0.95k ,
Bk = 0.95−k .
S. Nakamori et al. / Digital Signal Processing 16 (2006) 369–388
381
Fig. 1. Filter error variances, Σk,k , and smoother error variances, Σk,k+lag , with lag = 1 and for the values p = 0.2 and 0.6, when S = 0.15.
The observation model, affected by a random delay, is given by y˜k = zk + vk ,
k 0,
yk = (1 − γk )y˜k + γk y˜k−1 ,
k 1,
where the process {vk ; k 0} is a zero-mean white Gaussian noise with Var[vk ] = 9, for all k; we also assume that E[wk−1 vs ] = 0, for s = k, and E[wk−1 vk ] = S. Hence, the covariance function of the signal and the observation noise is expressed as Kzv (k, s) = S × 0.95k−s ,
0 s k,
and, consequently, Ck = 0.95k and Dk = S × 0.95−k . Finally, the noise {γk ; k 1} is a sequence of independent Bernoulli random variables with P [γk = 1] = p, for all k; that is, we assume that the probability of a delay in the measurements is constant at any time. We have performed one hundred iterations of the proposed algorithms, for different values of the delay probability p (specifically, p = 0.2 and 0.6) and different values of S, which lead to different degrees of correlation between the signal and the observation noise; the smoothing estimates (fixed-point, fixed-interval, and fixed-lag), as well as the corresponding error variances, have been calculated for these values.
382
S. Nakamori et al. / Digital Signal Processing 16 (2006) 369–388
Fig. 2. Single-stage smoother error variances, Σk,k+1 , fixed-interval smoother error variances, Σk,100 , k < 100, and Σ100,100 , for the values p = 0.2 and 0.6, when S = 0.15.
First, for the models corresponding to S = 0.15 and the above-mentioned values of the probability p, the smoother error variances were computed and compared with those of the filter. Figure 1 shows the filter error variances, Σk,k , and the smoother error variances, Σk,k+lag , with lag = 1 and 5. This figure shows that, for both values of probability p, the error variances corresponding to the smoothers are less than the filter ones and, also, that as probability p decreases, so do the error variances and, consequently, the performance of the estimators is better. From this figure, we also deduce that the accuracy of the smoother at each time k is better as the number of available observations increases. For the same models, the fixed-interval smoother error variances, Σk,100 , and the singlestage smoother ones, Σk,k+1 , for k < 100, are shown in Fig. 2 together with the error variances of the filter, Σ100,100 . Results similar to those of Fig. 1 are inferred. For the values p = 0.2 and S = 0.15, Fig. 3 shows a simulated signal together with the smoothing estimates zˆ k,k+lag , with lag = 1 and while Fig. 4 shows a simulated signal together with the single-stage smoothing estimates and the fixed-interval smoothing estimates, zˆ k,100 . According to the results given in Figs. 1 and 2, respectively, these figures show that the smoothing estimates follow the signal evolution better as the number of the observations used for estimation increases.
S. Nakamori et al. / Digital Signal Processing 16 (2006) 369–388
383
Fig. 3. Simulated signal zk , and smoothing estimates zˆ k,k+lag , with lag = 1 and for the values p = 0.2 and S = 0.15.
In order to compare the performance of the estimators for different degrees of correlation between the signal and the observation noise, we computed the filtering and smoothing error variances for different values of S, considering p = 0.2. The error variances Σk,k , Σk,k+1 , and Σk,k+5 , for S = 0.15 and are given in Fig. 5, while Fig. 6 shows the error variances for the single-stage smoother, Σk,k+1 , and the fixed-interval smoother Σk,100 , k < 100, and Σ100,100 , for S = 0.15 and 0.3 and for S = 0.15 and 0.2, respectively. These two figures suggest that error variances are smaller as the correlation S increases and, consequently, the performance of the estimators is better; these results are to be expected, since the correlation between the signal and the observations increases with S.
6. Conclusions In this paper, using covariance information, a general expression for the smoother, which leads to fixed-point, fixed-interval, and fixed-lag smoothing recursive algorithms, is derived from randomly delayed measurements. The random delay is modelled by a sequence of independent Bernoulli random variables whose parameters, which represent the probabilities of delay, are known. If a Bernoulli variable has value zero, the corresponding
384
S. Nakamori et al. / Digital Signal Processing 16 (2006) 369–388
Fig. 4. Simulated signal zk , single-stage smoothing estimates zˆ k,k+1 and fixed-interval smoothing estimates zˆ k,100 , for the values p = 0.2 and S = 0.15.
measurement is up-to-date; otherwise, if it takes value one, the measurement is delayed by one sampling period. We consider that the signal is correlated with the observation noise and it is assumed that the autocovariance functions of the signal and the crosscovariance function between the signal and the noise are expressed in a semi-degenerate kernel form. The proposed recursive algorithms are derived by applying the innovation technique and do not require knowledge of the state-space model of the signal, but only the covariance function of the signal and the observation white noise, the crosscovariance function between the signal and noise, the observed values and the probability of delay in the observations. The numerical simulation example considers the same discrete-time system with correlated disturbances and delayed observations as in Nakamori et al. [14], and shows that the smoothers proposed in this paper are computationally feasible.
Acknowledgment This work is partially supported by the ‘Ministerio de Ciencia y Tecnología’ under contract BFM2002-00932.
S. Nakamori et al. / Digital Signal Processing 16 (2006) 369–388
385
Fig. 5. Filter error variances Σk,k , and smoother error variances Σk,k+lag , with lag = 1 and for the value p = 0.2, when S = 0.15 and 0.3.
Appendix A. Proof of Theorem 1 Taking into account expression (2) for j = k, establishing the filter requires us to calT ], for 1 i k. In order to culate the coefficients Sk,i = E[zk νiT ] = E[zk yiT ] − E[zk yˆi,i−1 obtain these, we must first determine a formula for the one-stage predictor of the observation, which, from the innovation approach, is written as yˆi,i−1 =
i−1 E yi νjT Πj−1 νj ,
i 2;
yˆ1,0 = 0.
(A.1)
j =1
From the model hypotheses, it is deduced that E yi νjT = (1 − pi )Si,j + pi Si−1,j , j i − 2, T T = (1 − pi )Si,i−1 + pi Si−1,i−1 + pi (1 − pi−1 ) Di−1 Ci−1 + Ri−1 , E yi νi−1 and substituting these expressions in (A.1), we conclude that the one-stage predictor of the observation satisfies the following relation for i 2, yˆi,i−1 = (1 − pi )ˆzi,i−1 + pi zˆ i−1,i−1 −1 T + pi (1 − pi−1 ) Di−1 Ci−1 + Ri−1 Πi−1 νi−1 .
(A.2)
386
S. Nakamori et al. / Digital Signal Processing 16 (2006) 369–388
Fig. 6. Single-stage smoother error variances, Σk,k+1 , fixed-interval smoother Σk,100 , k < 100, and Σ100,100 , for the value p = 0.2, when S = 0.15 and 0.2.
Hence, T T Sk,i = E zk yiT − (1 − pi )E zk zˆ i,i−1 − pi E zk zˆ i−1,i−1 −1 T − pi (1 − pi−1 )Sk,i−1 Πi−1 Ci−1 Di−1 + Ri−1 , 2 i k. Using now (2) for zˆ i,i−1 and zˆ i−1,i−1 , and taking into account the hypotheses on the model, we have T Sk,i = (Ak | Ck )G(B − i |Di )
i−1 T T Sk,j Πj−1 (1 − pi )Si,j + pi Si−1,j j =1
− Sk,i−1 HiT , 2 i T . Sk,1 = (Ak | Ck )G(B 1 |D1 )
k;
This expression guarantees that Sk,i = (Ak | Ck )Ji ,
1 i k,
where J is a matrix function satisfying
(A.3)
S. Nakamori et al. / Digital Signal Processing 16 (2006) 369–388
T Ji = G(B − i |Di )
i−1 T T Jj Πj−1 (1 − pi )Si,j + pi Si−1,j − Ji−1 HiT ,
387
2 i k;
j =1
T . J1 = G(B 1 |D1 )
(A.4)
Therefore, if we denote Ok =
k
Ji Πi−1 νi ,
k 1;
O0 = 0,
(A.5)
i=1
it is clear, from (A.3) and (A.5), that the signal filter is given by (4). The recursive relation (5) for the vector Ok is immediate from (A.5). By an exactly analogous procedure to that used to obtain (4) the following expression for the linear one-stage predictor of the signal can be deduced zˆ k,k−1 = (Ak | Ck )Ok−1 ,
k1
(A.6)
and using the definition of νk and (A.2), expression (6) for the innovation is obtained. Expression (7) for Jk is established by putting i = k in (A.4), taking into account (A.3) and noting that rk = E[Ok OkT ] satisfies rk =
k
Jj Πj−1 JjT ,
r0 = 0.
(A.7)
j =1
From (A.7), the recursive relation (8) for the matrices rk is obvious. Finally, let us calculate the covariance matrix of the innovation process, Πk = E[νk νkT ], which, as a consequence of the Orthogonal Projection Lemma, can be expressed as T Πk = E yk ykT − E yˆk,k−1 yˆk,k−1 . Formula (9) for Πk is then easily deduced, taking into account the hypotheses on the model for E[yk ykT ], the relation (A.2) for yˆk,k−1 together with (4) and (A.6), the recursive expression (5) for the vectors Ok−1 and using the fact that Ok−2 is orthogonal to νk−1 . 2
References [1] S. Nakamori, R. Caballero-Águila, A. Hermoso-Carazo, J. Linares-Pérez, Linear recursive discrete-time estimators using covariance information under uncertain observations, Signal Process. 83 (2003) 1553– 1559. [2] S. Nakamori, R. Caballero-Águila, A. Hermoso-Carazo, J. Linares-Pérez, Linear estimation from uncertain observations with white plus coloured noises using covariance information, Digital Signal Process. 13 (2003) 552–568. [3] S. Nakamori, R. Caballero-Águila, A. Hermoso-Carazo, J. Linares-Pérez, Fixed-point smoothing with nonindependent uncertainty using covariance information, Int. J. Syst. Sci. 7 (10) (2003) 439–452. [4] S. Nakamori, R. Caballero-Águila, A. Hermoso-Carazo, J. Linares-Pérez, Fixed-interval smoothing algorithm based on covariances with correlation in the uncertainty, Digital Signal Process., in press. [5] J. Nilsson, B. Bernhardsson, B. Wittenmark, Stochastic analysis and control of real-time systems with random time delays, Automatica 34 (1998) 57–64. [6] I.V. Kolmanovsky, T.L. Maizenberg, Optimal control of continuous-time linear systems with a time-varying, random delay, Syst. Control Lett. 44 (2001) 119–126.
388
S. Nakamori et al. / Digital Signal Processing 16 (2006) 369–388
[7] J.S. Evans, V. Krishnamurthy, Hidden Markov model state estimation with randomly delayed observations, IEEE Trans. Signal Process. 47 (1999) 2157–2166. [8] Y.I. Yaz, E.E. Yaz, M.J. Mohseni, LMI-based state estimation for some discrete-time stochastic models, in: Proceedings of the 1998 IEEE International Conference on Control Applications, 1998, pp. 456–460. [9] C.L. Su, C.N. Lu, Interconnected network state estimation using randomly delayed measurements, IEEE Trans. Power Syst. 16 (2001) 870–878. [10] A.S. Matveev, A.V. Savkin, The problem of state estimation via asynchronous communication channels with irregular transmission times, IEEE Trans. Automat. Control 48 (2003) 670–676. [11] W.A. Gardner, A series solution to smoothing, filtering, and prediction problems involving correlated signal and noise, IEEE Trans. Inform. Theory 21 (1975) 698–699. [12] Y. Atse, K. Nakayama, Z. Ma, Performance of single- and multi-reference NLMS noise canceler based on correlation between signal and noise, IEICE Trans. Fundament. Electron. Commun. Comp. Sci. E78-A 11 (1995) 1576–1588. [13] A.G. Bhatt, R.L. Karandikar, Robustness of the nonlinear filter: the correlated case, Stochast. Process. Appl. 97 (2002) 41–58. [14] S. Nakamori, R. Caballero-Águila, A. Hermoso-Carazo, J. Linares-Pérez, Estimation algorithms from delayed measurements with correlation between signal and noise using covariance information, IEICE Trans. Fundament. Electron. Commun. Comput. Sci. E87-A (5) (2004) 1219–1225. [15] T. Kailath, Lectures on Linear Least-Squares Estimation, Springer-Verlag, New York, 1976.
Seiichi Nakamori was born in Kagoshima prefecture in 1951. He received the B.E. degree in electronic engineering from Kagoshima University in 1974, and the Dr. Eng. degree in applied mathematics and physics from Kyoto University in 1982. He was with the Faculty of Engineering, Ohita University, from 1974 to 1987, where he was promoted to a lecturer, Interdisciplinary Chair (Applied Mathematics) in 1985. He joined Kagoshima University in 1987, as an Associate Professor in the Department of Technology (Electric Technology), Faculty of Education, where he was promoted to a Professor of the graduate school of education in 1994. He is mainly interested in the stochastic signal estimation problem. He is a member of the Society of Instrument and Control Engineers, the Institute of Systems, Control and Information Engineers, and the Japanese Society of Technology Education. Aurora Hermoso-Carazo received the M.Sc. degree in Mathematics (Statistics) and her Ph.D. in Likelihood Test in Log Normal Processes, both from the University of Granada (Spain), in 1979 and 1984, respectively. In 1986, she became Associate Professor at the Department of Statistics and Operations Research, University of Granada (Spain). Her current research interests include stochastic systems, filtering and prediction. Josefa Linares-Pérez received the M.Sc. degree in Mathematics (Statistics) and her Ph.D. in Stochastic Differential Equations, both from the University of Granada (Spain), in 1980 and 1982, respectively. She is currently Professor at the Department of Statistics and Operations Research, University of Granada (Spain). Her research interest involves the areas of stochastic calculus and estimation in stochastic systems.
A general smoothing equation for signal estimation using randomly delayed observations in the correlated signal-noise case S. Nakamori a,∗ , A. Hermoso-Carazo b , J. Linares-Pérez b a Department of Technology, Faculty of Education, Kagoshima University, 1-20-6, Kohrimoto,
Kagoshima 890-0065, Japan b Departamento de Estadística e Investigación Operativa, Universidad de Granada, Campus Fuentenueva, s/n,
18071 Granada, Spain Available online 16 June 2005
Abstract This paper treats the least-squares linear smoothing problem for signal estimation using measurements contaminated by additive white noise correlated with the signal, with stochastic delays. We derive a general smoothing equation which is applied to obtain specific smoothing algorithms, which are referred in the signal estimation literature as fixed-point, fixed-interval, and fixed-lag smoothing. Using an innovation approach, the general smoothing equation is derived without requiring the whole knowledge of the state-space model generating the signal, but only covariance information of the signal and the observation noise, as well as the delay probabilities. © 2005 Elsevier Inc. All rights reserved. Keywords: Fixed-point smoothing; Fixed-interval smoothing; Fixed-lag smoothing; Covariance information; Delayed observations
* Corresponding author.
E-mail addresses: [email protected] (S. Nakamori), [email protected] (A. Hermoso-Carazo), [email protected] (J. Linares-Pérez). 1051-2004/$ – see front matter © 2005 Elsevier Inc. All rights reserved. doi:10.1016/j.dsp.2005.04.013
370
S. Nakamori et al. / Digital Signal Processing 16 (2006) 369–388
1. Introduction For a long time, least-squares estimation problems of random signals from observed data have been the subject of intensive research. There are several solution methods, which depend on the process assumptions and the models that represent the possible relationships between the unknown signal and the observable variables. If the observations of the signal are without stopping in time, there are many algorithms capable of obtaining the signal estimates (see, for example, [1–4] in an uncertain measurement context). However, many real-world problems exist in which errors due to transmission delays may be present and, hence, the measurements do not arrive in time but are randomly delayed. Such situations arise frequently in engineering applications and in certain economic studies (for example, the delay between submitting a patent application and actually obtaining the patent). Systems with random delays are considered in communication theory (Nilson et al. [5], discuss modelling and analysis of systems subject to random time-delays in the communication network and they present a method for different control schemes) and also in some control applications (Kolmanovsky and Maizenberg [6] consider a finite horizon optimal control problem for a class of linear systems with a randomly varying time-delay, affecting only the state variable, which is modelled by a Markov process with a finite number of states). In the cases in which measurements are randomly delayed, it is necessary to modify of the standard algorithms for the signal estimation problem in order to incorporate the effects of the random delay. Using different approaches, the state estimation problem for system models with randomly varying time-delays has been investigated and modifications of conventional signal estimation algorithms have been proposed to accommodate the randomly delayed observations. Let us cite, among other papers, the following: Evans and Krishnamurthy [7] consider state estimation for a discrete-time hidden Markov model in which the observations are delayed by a random time and the delay process is modelled as a finite state Markov chain. In Yaz et al. [8], state estimation in linear discrete-time state-space models with stochastic parameters is treated using linear matrix inequalities, and the results are applied to the problem of state estimation with a random sensor delay. Su and Lu [9] designed an extended Kalman filtering algorithm which provides optimal estimates of interconnected network states for systems in which some or all measurements are delayed. More recently, Matveev and Savkin [10] proposed a recursive minimum variance state estimator in a linear discrete-time partially-observed system perturbed by white noises when the observations are transmitted via communication channels with random transmission times and when various measurement signals may incur independent delays. In this paper, we consider the situation in which the state-space model of the signal is not available and the signal observations may be delayed by one sampling time; hence, the observed data may not contain the actual signal (randomly delayed observations). The delay is modelled by a sequence of independent Bernoulli random variables whose values, one or zero, indicate whether the measurements arrive on time or are delayed by one sampling time, respectively. Moreover, since the assumption that the signal process and the observation noise are independent can be restrictive in many important engineering and stochastic control problems [11–13], we consider the situation in which randomly delayed observations are perturbed by a white noise which is correlated with the signal process.
S. Nakamori et al. / Digital Signal Processing 16 (2006) 369–388
371
In Nakamori et al. [14], a recursive algorithm for the least-squares linear one-stage prediction and filtering problems was obtained under the above hypotheses concerning the model; in this paper we are interested in complementing the latter study by considering linear interpolation or smoothing problems, and derive recursive algorithms satisfied by smoothed estimates. A general smoothing equation is deduced from an innovation approach, which simplifies obtention since the innovation process is a white noise; from this equation we obtain the fixed-point, fixed-interval, and fixed-lag smoothing algorithms. The proposed algorithms do not require the whole knowledge of the state-space model generating the signal, but only covariance information of the signal and the observation noise, as well as the delay probabilities. Specifically, it is assumed that both, the autocovariance functions of the signal and the crosscovariance function between the signal and the observation noise are expressed in a semi-degenerate kernel form. The precision of the least-squares estimators is measured by the estimation error covariance matrices, which provide a global measurement of the performance of the estimators. For this reason, the formulas to obtain such matrices (for the fixed-point, fixed-interval, and fixed-lag smoothers) are also included in the present study.
2. Problem formulation We consider observations of a n × 1 signal zk described by the equation y˜k = zk + vk ,
k 0,
where {vk ; k 0} is the observation noise. It is assumed that the real observation at time k can be delayed by one sampling period, with a certain probability pk , or updated, with a probability 1 − pk . Therefore, if {γk ; k 1} denotes a sequence of independent Bernoulli random variables with P [γk = 1] = pk , the received measurement at time k is given by yk = (1 − γk )y˜k + γk y˜k−1 ,
k 1.
(1)
Remark. If γk = 1, then yk = y˜k−1 and the measurement is delayed by one sampling period. Otherwise, if γk = 0, then yk = y˜k , which means that the measurement is up-todate. Therefore, the value pk represents the probability of a delay in the measurement yk . In applications of communication networks, the noise {γk ; k > 1} usually represents the random delay from sensor to controller and the assumption of one-step sensor delay is based on the reasonable supposition that the induced data latency from the sensor to the controller is restricted so as not to exceed the sampling period [14]. To treat the least-squares (LS) estimation problem of the signal from the observations given by Eq. (1), we assume the following hypotheses on the signal and noise processes: (i) The signal process {zk ; k 0} has zero mean and its autocovariance function is expressed in a semi-degenerate kernel form as T Ak BsT , 0 s k, Kz (k, s) = E zk zs = Bk ATs , 0 k s,
372
S. Nakamori et al. / Digital Signal Processing 16 (2006) 369–388
where Ak and Bk are known n × M matrix functions. (ii) The noise process {vk ; k 0} is a zero-mean white sequence with known autocovariance function 0, if s = k, E vk vsT = Rk , if s = k. (iii) The signal process {zk ; k 0} and the noise process {vk ; k 0} are correlated, being T Kzv (k, s) = E zk vsT = Ck Ds , 0 s k, 0, 0 k < s, where Ck and Dk are known n × N matrix functions. (iv) The multiplicative noise {γk ; k 1} is a sequence of independent Bernoulli random variables with known probabilities, P [γk = 1] = pk > 0, ∀k 1. (v) The process {γk ; k 1} is independent of {zk ; k 0} and {vk ; k 0}. Remark. The special form of correlation in hypothesis (iii) allows us to consider the estimation problem of signals from randomly delayed observations when the evolution of the signal is influenced by the observation noise up to the considered time. This fact is shown in the numerical example (Section 5), where the disturbances of the signal and the observation are correlated at consecutive time instants and, consequently, the signal at time k is correlated with the observation noise up to the time k. Our aim is to approach the LS linear estimation problem of the signal zk , based on the observations {y1 , . . . , yL } with L > k; that is, the time k at which it is desired to estimate the signal is less than the time L of the last measurement. More specifically, we are interested in deriving of recursive smoothing algorithms which can be applied to various practical situations; specifically, we seek to obtain the smoothers of the signal at a fixed point, k, when the number of available measurements, L, is increasing (fixed-point smoother), the smoothers of the signal at different instants, k, based on a fixed number, L, of measurements (fixed-interval smoother) and the smoothers of the signal at each time, k, based on the observations up to the time k + λ, where λ is a fixed number (fixed-lag smoother). In order to simplify obtaining the algorithms, we use an innovation approach; the innovation at time i represents the new information which provides the observation at time i and is given by νi = yi − yˆi,i−1 , where yˆi,i−1 denotes the LS linear estimator of the observation yi based on the observations {y1 , . . . , yi−1 }. It is known (Kailath [15]) that the LS linear estimator of zk based on the observations {y1 , . . . , yj }, which is denoted by zˆ k,j , is equal to the LS linear estimator based on the innovations {ν1 , . . . , νj }. Hence, the fact that the innovations constitute a white process leads us to the following expression for the estimator of the signal (Nakamori et al. [4]), zˆ k,j =
j Sk,i Πi−1 νi , i=1
where
Sk,i = E zk νiT
and Πi = E νi νiT .
(2)
S. Nakamori et al. / Digital Signal Processing 16 (2006) 369–388
373
Taking into account expression (2), it is possible to obtain the smoothers, zˆ k,L , L k + 1, from the filter, by adding the remaining terms of the sum in (2), whose coefficients, Sk,i , i k + 1, usually have a common form (Nakamori et al. [4]); however, for the model at hand, the correlation between the signal and the noise specified in hypothesis (iii), together with the delay in the observations, leads to different expressions for the coefficients Sk,k+1 and Sk,i , i k + 2. Hence, in this case, it is more appropriate to start from the single-stage smoother, zˆ k,k+1 , and express the smoothers, zˆ k,L for L > k + 1 as zˆ k,L = zˆ k,k+1 +
L
Sk,i Πi−1 νi ,
L > k + 1.
(3)
i=k+2
Therefore, we first derive an algorithm to obtain zˆ k,k+1 and, afterwards, the different smoothing algorithms (fixed-point, fixed-interval, and fixed-lag) are derived from a general formula for the second term of the right member of (3).
3. Single-stage smoother To obtain the single-stage smoother at time k, we start from (2) for j = k + 1 and, by separating the last term of the sum, it is immediately obtained that that estimator verifies −1 νk+1 . zˆ k,k+1 = zˆ k,k + Sk,k+1 Πk+1
Then, since the filter, zˆ k,k , is required as an input for each k, the single-stage smoother is performed accompanied by the former; hence, we first present the recursive filtering algorithm established in Nakamori et al. [14]. Since some expressions of the proof of the filtering algorithm are used in the deduction of the smoothers, we present a brief proof of this algorithm, in Appendix A. Theorem 1. Consider the observation equation (1) verifying hypotheses (i)–(v); then the signal filter, zˆ k,k , is given by 1 zˆ k,k = (Ak | Ck )Ok ,
k 1,
(4)
where the vector Ok is recursively calculated from Ok = Ok−1 + Jk Πk−1 νk ,
k 1;
O0 = 0,
(5)
where νk is the innovation, which satisfies νk = yk − G(Ak |Ck ) Ok−1 − Hk νk−1 ,
k 1;
ν0 = 0,
(6)
and T T Jk = G(B − rk−1 G(A − Jk−1 HkT , k |Dk ) k |Ck )
k 1;
J0 = 0,
(7)
where the matrices GYk (Yk = (Ak | Ck ), (Bk | Dk )) and Hk are given by −1 T GYk = (1 − pk )Yk + pk Yk−1 , Hk = pk (1 − pk−1 ) Dk−1 Ck−1 + Rk−1 Πk−1 . 1 (U | V ) denotes a partitioned matrix into sub-matrices U and V .
374
S. Nakamori et al. / Digital Signal Processing 16 (2006) 369–388
The matrix rk = E[Ok OkT ] is recursively calculated from rk = rk−1 + Jk Πk−1 JkT ,
k 1;
r0 = 0,
(8)
and Πk , the covariance matrix of the innovation νk , is given by Π = (1 − pk ) Ak BkT + Ck DkT + Dk CkT + Rk T T T + Ck−1 Dk−1 + Dk−1 Ck−1 + Rk−1 + pk Ak−1 Bk−1 T T T − G(Ak |Ck ) rk−1 G(A +Jk−1 HkT − Hk Jk−1 G(A + Πk−1 HkT , k |Ck ) k |Ck )
k 1;
Π0 = 0.
(9)
Proof. The proof of Theorem 1 is given in Appendix A.
2
In the following theorem, we present an algorithm for computing zˆ k,k+1 , the singlestage smoothed estimator of the signal zk , for any given k 1. Theorem 2. Consider the observation equation (1), satisfying hypotheses (i)–(v). Then, the single-stage smoother of the signal zk can be obtained from −1 νk+1 , zˆ k,k+1 = zˆ k,k + Sk,k+1 Πk+1
k 1,
(10)
where 2 T Sk,k+1 = (Bk | 0n×N ) − (Ak | Ck )rk G(A k+1 |Ck+1 ) T − (Ak | Ck )Jk Hk+1 + pk+1 Ck DkT .
(11)
The filter zˆ k,k , the matrices G(Ak |Ck ) , Hk , rk and Jk , the innovations νk and their covariances Πk are obtained from the formulas given in Theorem 1. Proof. As indicated above, Eq. (10) for the single-stage smoother is immediately obtained from (2) for j = k + 1. T ], we use formula (6) for the In order to calculate expression (11) for Sk,k+1 = E[zk νk+1 innovation νk+1 and we obtain T T T − E zk νkT Hk+1 . − E zk OkT G(A Sk,k+1 = E zk yk+1 k+1 |Ck+1 ) T ] and expression (A.5) for O , it is Then, using the hypotheses on the model for E[zk yk+1 k deduced that T Sk,k+1 = (Bk | 0n×N )G(A + pk+1 Ck DkT k+1 |Ck+1 )
−
k T T Sk,j Πj−1 JjT G(A − Sk,k Hk+1 . k+1 |Ck+1 ) j =1
Now, using (A.3) for Sk,j , 1 j k and afterwards (A.7), the formula (11) is proven. 2 0 a×b denotes the a × b zero matrix.
2
S. Nakamori et al. / Digital Signal Processing 16 (2006) 369–388
375
Error covariance matrices. The precision of the filter and the single-stage smoother is measured by the covariance matrices of the corresponding estimation errors, T Σk,j = E zk − zˆ k,j zk − zˆ k,j , j = k, k + 1. Since the error zk − zˆ k,j is orthogonal to the estimator zˆ k,j , it is immediately verified that T Σk,j = Kz (k, k) − E zˆ k,j zˆ k,j . Then, hypothesis (i) on the model and expressions (4) and (10) for the filter and the singlestage smoother, respectively, lead to Σk,k = Ak BkT − (Ak | Ck )rk (Ak | Ck )T , −1 T Σk,k+1 = Σk,k − Sk,k+1 Πk+1 Sk,k+1 ,
k 1,
(12)
k 1.
(13)
4. General expression for the linear smoother Once the algorithm to obtain the single-stage smoother, zˆ k,k+1 , has been deduced, our purpose is to derive a general expression for the remaining smoothers, zˆ k,L , L > k + 1. For
−1 this purpose, as expressed in (3), we must determine the sum L i=k+2 Sk,i Πi νi or, more T specifically, the coefficients Sk,i = [zk νi ], for i k + 2 which, from expression (6) for νi , are given by T T T T G(Ai |Ci ) − E zk νi−1 Hi . Sk,i = E zk yiT − E zk Oi−1 Taking into account the hypotheses on the model and (A.5), it is obtained that T − Sk,i = (Bk | 0n×N )G(A i |Ci )
i−1 T Sk,j Πj−1 JjT G(A − Sk,i−1 HiT , i |Ci )
i k + 2.
j =1
Then, using (A.3) and considering (A.7) for rk , these coefficients are expressed as follows: −1 T T Jk+1 G(Ai |Ci ) Sk,i = (Bk | 0n×N ) − (Ak | Ck )rk − Sk,k+1 Πk+1 −
i−1 j =k+2
T Sk,j Πj−1 JjT G(A − Sk,i−1 HiT , i |Ci )
i > k + 2;
−1 T T T Jk+1 G(Ak+2 |Ck+2 ) − Sk,k+1 Hk+2 . Sk,k+2 = (Bk | 0n×N ) − (Ak | Ck )rk − Sk,k+1 Πk+1 This latter expression for Sk,i allows us to deduce that Sk,i = Mk k,i ,
i k + 2,
(14)
where
−1 T Jk+1 − Sk,k+1 Mk = (Bk | 0n×N ) − (Ak | Ck )rk − Sk,k+1 Πk+1
and k,i verifying
(15)
376
S. Nakamori et al. / Digital Signal Processing 16 (2006) 369–388
k,i = (G(Ai |Ci ) | 0n×n )T −
i−1 j =k+2
T k,j Πj−1 JjT G(A − k,i−1 HiT , i |Ci )
i > k + 2;
T k,k+2 = G(Ak+2 |Ck+2 ) | Hk+2 .
(16)
Hence, if we define qk,L =
L
k,i Πi−1 νi ,
L > k + 1;
qk,k+1 = 0,
(17)
i=k+2
the following general expression for the smoother is obtained from (3), using (14), zˆ k,L = zˆ k,k+1 + Mk qk,L ,
L k + 1.
(18)
General expression for the smoother error covariance matrices. Taking into account that the error covariance matrices of the smoothers can be expressed as Σk,L = Kz (k, k) − T ], the general expression (18) for the estimators leads to the following general E[ˆzk,L zˆ k,L expression for such matrices Σk,L = Σk,k+1 − Mk Qk,L MTk ,
L k + 1,
(19)
where Σk,k+1 , the single-stage smoother error covariance matrices, are given in (13) and T ]. Qk,L = E[qk,L qk,L Next, different smoothing algorithms are obtained from the general expression (18); specifically, the smoothing estimators of the signal are recursively obtained in three different ways: • Fixed-point smoothing; a fixed integer k 1 is considered and the estimators zˆ k,L with L = k + 1, k + 2, . . . , are obtained (recursions for increasing L are proposed in Section 4.1). • Fixed-interval smoothing; the estimators zˆ k,L with L a fixed positive integer and k = L − 1, L − 2, . . . , 1, are derived (recursions in k are established in Section 4.2). • Fixed-lag smoothing; here, the obtained estimators are zˆ k,k+λ , where λ is a positive integer and k = 1, 2, . . . (recursions for increasing k are obtained in Section 4.3). Recursive relationships for the error covariance matrices of the smoothers in the three above cases are obtained in the respective sections mentioned above. 4.1. Fixed-point smoothing algorithm In this section we wish to determine a recursive algorithm to implement the computation of zˆ k,L for a fixed k 1, and for any L k + 1; it is not necessary to know, a priori, the time point of the final observation, since the algorithm proceeds forward recursively in time and can be stopped at will. So, taking into account the general expression (18), to solve this problem we need to obtain a forward recursive formula in L for qk,L . From definition (17) it is immediately deduced that qk,L = qk,L−1 + k,L ΠL−1 νL ,
L > k + 1;
qk,k+1 = 0.
(20)
S. Nakamori et al. / Digital Signal Processing 16 (2006) 369–388
377
Now, if we denote δk,i =
i
k,j Πj−1 JjT ,
i > k + 1;
δk,k+1 = 0,
(21)
j =k+2
expression (16) is rewritten as 3 T k,L = (IM+N | 0(M+N )×n )T − δk,L−1 G(A − k,L−1 HLT , L |CL ) k,k+2 = (G(Ak+2 |Ck+2 ) | Hk+2 ) , T
L > k + 2; (22)
and, from (21), it is clear that the function δ verifies δk,L = δk,L−1 + k,L ΠL−1 JLT ,
L > k + 1;
δk,k+1 = 0.
(23)
So, Eqs. (20), (22) and (23) provide the desired forward recursive formula in L for qk,L . Fixed-point smoother error covariance matrices. The error covariance matrices of the smoother, Σk,L , are given in (19). Now, we need to obtain a forward recursive expression T ]. in L for Qk,L = E[qk,L qk,L Since qk,L−1 and νL are uncorrelated, the following expression is immediately obtained from (20), Qk,L = Qk,L−1 + k,L ΠL−1 Tk,L ,
L > k + 1;
Qk,k+1 = 0.
(24)
Remark. From (2), we have the following alternate form of the fixed-point smoothing algorithm zˆ k,L = zˆ k,L−1 + Sk,L ΠL−1 νL ,
L > k + 1,
(25)
whose initial condition is the single-stage smoother, zˆ k,k+1 , given in Theorem 2. The matrices Sk,L can be calculated from (14); that is, Sk,L = Mk k,L ,
L > k + 1,
where Mk is given in (15) and k,L can be recursively calculated from expressions (22) and (23). Finally, from (25), we also have the following alternate form for obtaining the fixedpoint smoothing error covariance matrices, T Σk,L = Σk,L−1 − Sk,L ΠL−1 Sk,L ,
L > k + 1,
(26)
with initial condition Σk,k+1 given in (13). 4.2. Fixed-interval smoothing algorithm In this section we assume that the measurement data over the fixed interval [1, L] are available and that, for each time point k within the interval, we wish to obtain the LS linear estimator of the signal zk based on all the observations {y1 , . . . , yL }; thus, since 3 I denotes the a × a identity matrix. a
378
S. Nakamori et al. / Digital Signal Processing 16 (2006) 369–388
the estimators of zL and zL−1 are the filter, zˆ L,L , and the single-stage smoother, zˆ L−1,L , respectively, our purpose is to obtain the remaining estimators, zˆ k,L , for k < L − 1. For this purpose, we start from the general expression (18) and the fixed-interval smoothing algorithm, whose starting point is zˆ L−1,L , is obtained by deducing a backward recursive formula in k (k < L − 1) for qk,L . From (16) and (17), it is immediate that −1 qk,L = (G(Ak+2 |Ck+2 ) | Hk+2 )T Πk+2 νk+2 +
L
k,i Πi−1 νi ,
k < L − 2;
i=k+3
qL−2,L = (G(AL |CL ) | HL )T ΠL−1 νL . (1)T
Now, we express k,i = (k,i (a)
(27)
(2)T
| k,i )T and calculate, using (16), the differences
(a)
k,i − k+1,i (a = 1, 2), for i k + 3; by comparing the resultant expressions with (a)
those obtained from (16) for k+1,i , we deduce the following relations for i k + 3, (1) (1) (2) T T k,i = IM+N − G(A Π −1 J T k+1,i − G(A , k+2 |Ck+2 ) k+2 k+2 k+2 |Ck+2 ) k+1,i −1 T T T k,i = −Hk+2 Πk+2 Jk+2 k+1,i − Hk+2 k+1,i . (2)
(1)
(2)
Hence, we have k,i = Φk+2 k+1,i , where Φk+2 =
i k + 3,
(28)
T Π −1 J T IM+N − G(A k+2 |Ck+2 ) k+2 k+2
T −G(A k+2 |Ck+2 )
T Π −1 J T −Hk+2 k+2 k+2
T −Hk+2
,
and by substituting (28) in (27), we obtain the desired backward recursive formula for the vectors qk,L , −1 qk,L = Φk+2 qk+1,L + (G(Ak+2 |Ck+2 ) | Hk+2 )T Πk+2 νk+2 ,
qL−2,L = (G(AL |CL ) | HL )
T
ΠL−1 νL .
k < L − 2; (29)
Fixed-interval smoother error covariance matrices. The formula for the error covariance matrices of the fixed-interval smoother, Σk,L , is obtained from (19), using the following T ], which is obtained from (29) taking into account recursive relation for Qk,L = E[qk,L qk,L that qk+1,L and νk+2 are uncorrelated, T Qk,L = Φk+2 Qk+1,L Φk+2 + (G(Ak+2 |Ck+2 ) | Hk+2 )T −1 × Πk+2 (G(Ak+2 |Ck+2 ) | Hk+2 ),
QL−2,L = (G(AL |CL ) | HL )
T
k < L − 2;
ΠL−1 (G(AL |CL ) | HL ).
(30)
4.3. Fixed-lag smoothing algorithm We wish now to obtain an algorithm for the smoothers zˆ k,k+λ , where k = 1, 2, . . . , and λ 1 is a positive integer which represents a fixed-lag of λ units of time.
S. Nakamori et al. / Digital Signal Processing 16 (2006) 369–388
379
Note that, in the case λ = 1, the request smoothers, zˆ k,k+1 , can be obtained from the filtering algorithm as is proposed in Theorem 2. On the other hand, if we consider λ = 2, in view of general expression (18), using definition (17) for qk,k+2 and (16) for k,k+2 , we have −1 νk+2 , zˆ k,k+2 = zˆ k,k+1 + Mk (G(Ak+2 |Ck+2 ) | Hk+2 )T Πk+2
k 1.
Hence, we must only consider the cases λ 3. Like the fixed-point and fixed-interval smoothing algorithms, the fixed-lag one is again obtained from the general expression (18); hence, the smoothers zˆ k,k+λ are calculated from the following relation: zˆ k,k+λ = zˆ k,k+1 + Mk qk,k+λ ,
k 1,
(31)
using a recursive expression in k for qk,k+λ . For k 2, we can use definition (17) to calculate and compare the vectors qk,k+λ and qk−1,k−1+λ ; then, by considering the expressions (16) for k−1,k+1 and (28) to deduce that k−1,i = Φk+1 k,i , i k + 2, we obtain the following relation between the above two vectors, −1 −1 νk+1 qk−1,k−1+λ − (G(Ak+1 |Ck+1 ) | Hk+1 )T Πk+1 qk,k+λ = Φk+1 −1 + k−1,k+λ Πk+λ νk+λ , k 2. Thus, since −1 Φk+1 =
IM+N
T −G(A H −T k+1 |Ck+1 ) k+1
−1 T −Πk+1 Jk+1
−1 T −T T −[In − Πk+1 Jk+1 G(A ]Hk+1 k+1 |Ck+1 )
and, hence, −1 (G(Ak+1 |Ck+1 ) | Hk+1 )T = (0n×(M+N ) | −In )T , Φk+1
the above equation for qk,k+λ can be rewritten as −1 −1 qk−1,k−1+λ + (0n×(M+N ) | In )T Πk+1 νk+1 qk,k+λ = Φk+1 −1 + k,k+λ Πk+λ νk+λ ,
k 2.
(32)
In order to initiate this recursive relation, the vector q1,1+λ must be introduced into Eq. (32) at time k = 2; this vector can be obtained from expression (20), starting with q1,2 = 0 and processing the observations to obtain q1,3 , q1,4 , . . . , q1,1+λ . In view of (32), the fixed-lag smoothing algorithm also requires a recursive expression in k for k,k+λ , which is obtained using (21), (22), and (28), −1 −1 T δk−1,k−1+λ − (0n×(M+N ) | In )T Πk+1 Jk+1 k,k+λ = (IM+N | 0(M+N )×n )T − Φk+1 −1 T T × G(A − Φk+1 k−1,k−1+λ Hk+λ , k+λ |Ck+λ )
k 2.
(33)
The start point for (33), 1,1+λ , is obtained, starting with 1,3 = (G(A3 |C3 ) | H3 )T , from expressions (22) and (23), to obtain 1,4 , 1,5 , . . . , 1,1+λ , and δ1,1+λ is obtained from (23), starting with δ1,2 = 0. Finally, taking into account (33), it is also necessary to determine a recursive formula in k for δk,k+λ , k 1. From (17) and (A.5), since the innovation is a white process, it is clear
380
S. Nakamori et al. / Digital Signal Processing 16 (2006) 369–388
T ]. Then, by using (32) that δk,k+λ , as defined in (21), verifies that δk,k+λ = E[qk,k+λ Ok+λ for qk,k+λ , (5) for Ok+λ and taking into account that νk+λ and qk−1,k−1+λ are uncorrelated, we deduce the following relation: −1 −1 T δk−1,k−1+λ + (0n×(M+N ) | In )T Πk+1 Jk+1 δk,k+λ = Φk+1 −1 T + k,k+λ Πk+λ Jk+λ ,
k 2,
(34)
where δ1,1+λ is obtained from (23), starting with δ1,2 = 0. Equations (32)–(34) provide the desired recursive expression in k for qk,k+λ . Fixed-lag smoother error covariance matrices. From (19), the formula for the error covariance matrices of the fixed-lag smoother, zˆ k,k+λ , is Σk,k+λ = Σk,k+1 − Mk Qk,k+λ MTk ,
k 1,
(35)
T with Qk,k+λ = E[qk,k+λ qk,k+λ ]. It is clear that, for λ = 1, Qk,k+1 = 0 and, for λ = 2, −1 Qk,k+2 = (G(Ak+2 |Ck+2 ) | Hk+2 )T Πk+2 (G(Ak+2 |Ck+2 ) | Hk+2 ).
The recursive formula to obtain Qk,k+λ for λ 3 is derived from (32), taking into acT ]= T count that E[qk−1,k−1+λ νk+1 k−1,k+1 = (G(Ak+1 |Ck+1 ) | Hk+1 ) and E[qk−1,k−1+λ × T νk+λ ] = 0; so, −1 −T −1 Qk,k+λ = Φk+1 Qk−1,k−1+λ Φk+1 − (0n×(M+N ) | In )T Πk+1 (0n×(M+N ) | In ) −1 Tk,k+λ , + k,k+λ Πk+λ
k 2,
(36)
and the initial condition, Q1,1+λ , is obtained from expression (24), starting with Q1,2 = 0.
5. Numerical example In this section, we present a numerical simulation example to illustrate the application of the different smoothing algorithms proposed in the paper. For this purpose, we have performed a program in MatLab which simulates, at each iteration, the signal and the observation values and provides the smoothing estimates as well as the smoother error covariance matrices. We have considered the same example as in Nakamori et al. [14]; specifically, the signal {zk ; k 0} is a scalar process generated by the following first-order autoregressive model, zk+1 = 0.95zk + wk ,
(37)
where {wk ; k 0} is a zero-mean white Gaussian noise with Var[wk ] = 0.1, for all k. The covariance function of the signal is given by Kz (k, s) = 1.025641 × 0.95k−s ,
0 s k,
which can be expressed in a semi-degenerate kernel form, with Ak = 1.025641 × 0.95k ,
Bk = 0.95−k .
S. Nakamori et al. / Digital Signal Processing 16 (2006) 369–388
381
Fig. 1. Filter error variances, Σk,k , and smoother error variances, Σk,k+lag , with lag = 1 and for the values p = 0.2 and 0.6, when S = 0.15.
The observation model, affected by a random delay, is given by y˜k = zk + vk ,
k 0,
yk = (1 − γk )y˜k + γk y˜k−1 ,
k 1,
where the process {vk ; k 0} is a zero-mean white Gaussian noise with Var[vk ] = 9, for all k; we also assume that E[wk−1 vs ] = 0, for s = k, and E[wk−1 vk ] = S. Hence, the covariance function of the signal and the observation noise is expressed as Kzv (k, s) = S × 0.95k−s ,
0 s k,
and, consequently, Ck = 0.95k and Dk = S × 0.95−k . Finally, the noise {γk ; k 1} is a sequence of independent Bernoulli random variables with P [γk = 1] = p, for all k; that is, we assume that the probability of a delay in the measurements is constant at any time. We have performed one hundred iterations of the proposed algorithms, for different values of the delay probability p (specifically, p = 0.2 and 0.6) and different values of S, which lead to different degrees of correlation between the signal and the observation noise; the smoothing estimates (fixed-point, fixed-interval, and fixed-lag), as well as the corresponding error variances, have been calculated for these values.
382
S. Nakamori et al. / Digital Signal Processing 16 (2006) 369–388
Fig. 2. Single-stage smoother error variances, Σk,k+1 , fixed-interval smoother error variances, Σk,100 , k < 100, and Σ100,100 , for the values p = 0.2 and 0.6, when S = 0.15.
First, for the models corresponding to S = 0.15 and the above-mentioned values of the probability p, the smoother error variances were computed and compared with those of the filter. Figure 1 shows the filter error variances, Σk,k , and the smoother error variances, Σk,k+lag , with lag = 1 and 5. This figure shows that, for both values of probability p, the error variances corresponding to the smoothers are less than the filter ones and, also, that as probability p decreases, so do the error variances and, consequently, the performance of the estimators is better. From this figure, we also deduce that the accuracy of the smoother at each time k is better as the number of available observations increases. For the same models, the fixed-interval smoother error variances, Σk,100 , and the singlestage smoother ones, Σk,k+1 , for k < 100, are shown in Fig. 2 together with the error variances of the filter, Σ100,100 . Results similar to those of Fig. 1 are inferred. For the values p = 0.2 and S = 0.15, Fig. 3 shows a simulated signal together with the smoothing estimates zˆ k,k+lag , with lag = 1 and while Fig. 4 shows a simulated signal together with the single-stage smoothing estimates and the fixed-interval smoothing estimates, zˆ k,100 . According to the results given in Figs. 1 and 2, respectively, these figures show that the smoothing estimates follow the signal evolution better as the number of the observations used for estimation increases.
S. Nakamori et al. / Digital Signal Processing 16 (2006) 369–388
383
Fig. 3. Simulated signal zk , and smoothing estimates zˆ k,k+lag , with lag = 1 and for the values p = 0.2 and S = 0.15.
In order to compare the performance of the estimators for different degrees of correlation between the signal and the observation noise, we computed the filtering and smoothing error variances for different values of S, considering p = 0.2. The error variances Σk,k , Σk,k+1 , and Σk,k+5 , for S = 0.15 and are given in Fig. 5, while Fig. 6 shows the error variances for the single-stage smoother, Σk,k+1 , and the fixed-interval smoother Σk,100 , k < 100, and Σ100,100 , for S = 0.15 and 0.3 and for S = 0.15 and 0.2, respectively. These two figures suggest that error variances are smaller as the correlation S increases and, consequently, the performance of the estimators is better; these results are to be expected, since the correlation between the signal and the observations increases with S.
6. Conclusions In this paper, using covariance information, a general expression for the smoother, which leads to fixed-point, fixed-interval, and fixed-lag smoothing recursive algorithms, is derived from randomly delayed measurements. The random delay is modelled by a sequence of independent Bernoulli random variables whose parameters, which represent the probabilities of delay, are known. If a Bernoulli variable has value zero, the corresponding
384
S. Nakamori et al. / Digital Signal Processing 16 (2006) 369–388
Fig. 4. Simulated signal zk , single-stage smoothing estimates zˆ k,k+1 and fixed-interval smoothing estimates zˆ k,100 , for the values p = 0.2 and S = 0.15.
measurement is up-to-date; otherwise, if it takes value one, the measurement is delayed by one sampling period. We consider that the signal is correlated with the observation noise and it is assumed that the autocovariance functions of the signal and the crosscovariance function between the signal and the noise are expressed in a semi-degenerate kernel form. The proposed recursive algorithms are derived by applying the innovation technique and do not require knowledge of the state-space model of the signal, but only the covariance function of the signal and the observation white noise, the crosscovariance function between the signal and noise, the observed values and the probability of delay in the observations. The numerical simulation example considers the same discrete-time system with correlated disturbances and delayed observations as in Nakamori et al. [14], and shows that the smoothers proposed in this paper are computationally feasible.
Acknowledgment This work is partially supported by the ‘Ministerio de Ciencia y Tecnología’ under contract BFM2002-00932.
S. Nakamori et al. / Digital Signal Processing 16 (2006) 369–388
385
Fig. 5. Filter error variances Σk,k , and smoother error variances Σk,k+lag , with lag = 1 and for the value p = 0.2, when S = 0.15 and 0.3.
Appendix A. Proof of Theorem 1 Taking into account expression (2) for j = k, establishing the filter requires us to calT ], for 1 i k. In order to culate the coefficients Sk,i = E[zk νiT ] = E[zk yiT ] − E[zk yˆi,i−1 obtain these, we must first determine a formula for the one-stage predictor of the observation, which, from the innovation approach, is written as yˆi,i−1 =
i−1 E yi νjT Πj−1 νj ,
i 2;
yˆ1,0 = 0.
(A.1)
j =1
From the model hypotheses, it is deduced that E yi νjT = (1 − pi )Si,j + pi Si−1,j , j i − 2, T T = (1 − pi )Si,i−1 + pi Si−1,i−1 + pi (1 − pi−1 ) Di−1 Ci−1 + Ri−1 , E yi νi−1 and substituting these expressions in (A.1), we conclude that the one-stage predictor of the observation satisfies the following relation for i 2, yˆi,i−1 = (1 − pi )ˆzi,i−1 + pi zˆ i−1,i−1 −1 T + pi (1 − pi−1 ) Di−1 Ci−1 + Ri−1 Πi−1 νi−1 .
(A.2)
386
S. Nakamori et al. / Digital Signal Processing 16 (2006) 369–388
Fig. 6. Single-stage smoother error variances, Σk,k+1 , fixed-interval smoother Σk,100 , k < 100, and Σ100,100 , for the value p = 0.2, when S = 0.15 and 0.2.
Hence, T T Sk,i = E zk yiT − (1 − pi )E zk zˆ i,i−1 − pi E zk zˆ i−1,i−1 −1 T − pi (1 − pi−1 )Sk,i−1 Πi−1 Ci−1 Di−1 + Ri−1 , 2 i k. Using now (2) for zˆ i,i−1 and zˆ i−1,i−1 , and taking into account the hypotheses on the model, we have T Sk,i = (Ak | Ck )G(B − i |Di )
i−1 T T Sk,j Πj−1 (1 − pi )Si,j + pi Si−1,j j =1
− Sk,i−1 HiT , 2 i T . Sk,1 = (Ak | Ck )G(B 1 |D1 )
k;
This expression guarantees that Sk,i = (Ak | Ck )Ji ,
1 i k,
where J is a matrix function satisfying
(A.3)
S. Nakamori et al. / Digital Signal Processing 16 (2006) 369–388
T Ji = G(B − i |Di )
i−1 T T Jj Πj−1 (1 − pi )Si,j + pi Si−1,j − Ji−1 HiT ,
387
2 i k;
j =1
T . J1 = G(B 1 |D1 )
(A.4)
Therefore, if we denote Ok =
k
Ji Πi−1 νi ,
k 1;
O0 = 0,
(A.5)
i=1
it is clear, from (A.3) and (A.5), that the signal filter is given by (4). The recursive relation (5) for the vector Ok is immediate from (A.5). By an exactly analogous procedure to that used to obtain (4) the following expression for the linear one-stage predictor of the signal can be deduced zˆ k,k−1 = (Ak | Ck )Ok−1 ,
k1
(A.6)
and using the definition of νk and (A.2), expression (6) for the innovation is obtained. Expression (7) for Jk is established by putting i = k in (A.4), taking into account (A.3) and noting that rk = E[Ok OkT ] satisfies rk =
k
Jj Πj−1 JjT ,
r0 = 0.
(A.7)
j =1
From (A.7), the recursive relation (8) for the matrices rk is obvious. Finally, let us calculate the covariance matrix of the innovation process, Πk = E[νk νkT ], which, as a consequence of the Orthogonal Projection Lemma, can be expressed as T Πk = E yk ykT − E yˆk,k−1 yˆk,k−1 . Formula (9) for Πk is then easily deduced, taking into account the hypotheses on the model for E[yk ykT ], the relation (A.2) for yˆk,k−1 together with (4) and (A.6), the recursive expression (5) for the vectors Ok−1 and using the fact that Ok−2 is orthogonal to νk−1 . 2
References [1] S. Nakamori, R. Caballero-Águila, A. Hermoso-Carazo, J. Linares-Pérez, Linear recursive discrete-time estimators using covariance information under uncertain observations, Signal Process. 83 (2003) 1553– 1559. [2] S. Nakamori, R. Caballero-Águila, A. Hermoso-Carazo, J. Linares-Pérez, Linear estimation from uncertain observations with white plus coloured noises using covariance information, Digital Signal Process. 13 (2003) 552–568. [3] S. Nakamori, R. Caballero-Águila, A. Hermoso-Carazo, J. Linares-Pérez, Fixed-point smoothing with nonindependent uncertainty using covariance information, Int. J. Syst. Sci. 7 (10) (2003) 439–452. [4] S. Nakamori, R. Caballero-Águila, A. Hermoso-Carazo, J. Linares-Pérez, Fixed-interval smoothing algorithm based on covariances with correlation in the uncertainty, Digital Signal Process., in press. [5] J. Nilsson, B. Bernhardsson, B. Wittenmark, Stochastic analysis and control of real-time systems with random time delays, Automatica 34 (1998) 57–64. [6] I.V. Kolmanovsky, T.L. Maizenberg, Optimal control of continuous-time linear systems with a time-varying, random delay, Syst. Control Lett. 44 (2001) 119–126.
388
S. Nakamori et al. / Digital Signal Processing 16 (2006) 369–388
[7] J.S. Evans, V. Krishnamurthy, Hidden Markov model state estimation with randomly delayed observations, IEEE Trans. Signal Process. 47 (1999) 2157–2166. [8] Y.I. Yaz, E.E. Yaz, M.J. Mohseni, LMI-based state estimation for some discrete-time stochastic models, in: Proceedings of the 1998 IEEE International Conference on Control Applications, 1998, pp. 456–460. [9] C.L. Su, C.N. Lu, Interconnected network state estimation using randomly delayed measurements, IEEE Trans. Power Syst. 16 (2001) 870–878. [10] A.S. Matveev, A.V. Savkin, The problem of state estimation via asynchronous communication channels with irregular transmission times, IEEE Trans. Automat. Control 48 (2003) 670–676. [11] W.A. Gardner, A series solution to smoothing, filtering, and prediction problems involving correlated signal and noise, IEEE Trans. Inform. Theory 21 (1975) 698–699. [12] Y. Atse, K. Nakayama, Z. Ma, Performance of single- and multi-reference NLMS noise canceler based on correlation between signal and noise, IEICE Trans. Fundament. Electron. Commun. Comp. Sci. E78-A 11 (1995) 1576–1588. [13] A.G. Bhatt, R.L. Karandikar, Robustness of the nonlinear filter: the correlated case, Stochast. Process. Appl. 97 (2002) 41–58. [14] S. Nakamori, R. Caballero-Águila, A. Hermoso-Carazo, J. Linares-Pérez, Estimation algorithms from delayed measurements with correlation between signal and noise using covariance information, IEICE Trans. Fundament. Electron. Commun. Comput. Sci. E87-A (5) (2004) 1219–1225. [15] T. Kailath, Lectures on Linear Least-Squares Estimation, Springer-Verlag, New York, 1976.
Seiichi Nakamori was born in Kagoshima prefecture in 1951. He received the B.E. degree in electronic engineering from Kagoshima University in 1974, and the Dr. Eng. degree in applied mathematics and physics from Kyoto University in 1982. He was with the Faculty of Engineering, Ohita University, from 1974 to 1987, where he was promoted to a lecturer, Interdisciplinary Chair (Applied Mathematics) in 1985. He joined Kagoshima University in 1987, as an Associate Professor in the Department of Technology (Electric Technology), Faculty of Education, where he was promoted to a Professor of the graduate school of education in 1994. He is mainly interested in the stochastic signal estimation problem. He is a member of the Society of Instrument and Control Engineers, the Institute of Systems, Control and Information Engineers, and the Japanese Society of Technology Education. Aurora Hermoso-Carazo received the M.Sc. degree in Mathematics (Statistics) and her Ph.D. in Likelihood Test in Log Normal Processes, both from the University of Granada (Spain), in 1979 and 1984, respectively. In 1986, she became Associate Professor at the Department of Statistics and Operations Research, University of Granada (Spain). Her current research interests include stochastic systems, filtering and prediction. Josefa Linares-Pérez received the M.Sc. degree in Mathematics (Statistics) and her Ph.D. in Stochastic Differential Equations, both from the University of Granada (Spain), in 1980 and 1982, respectively. She is currently Professor at the Department of Statistics and Operations Research, University of Granada (Spain). Her research interest involves the areas of stochastic calculus and estimation in stochastic systems.
