An optimal deconvolution smoother for systems with random parametric uncertainty and its application to semi-blind deconvolution

Signal Processing 92 (2012) 2497–2508

Contents lists available at SciVerse ScienceDirect

Signal Processing journal homepage: www.elsevier.com/locate/sigpro

An optimal deconvolution smoother for systems with random parametric uncertainty and its application to semi-blind deconvolution Chengpu Yu a, Nan Xiao a, Cishen Zhang b,n, Lihua Xie a a b

EXQUISITUS, Centre for E-City, School of Electrical and Electronic Engineering, Nangyang Technological University, Singapore 639798, Singapore Faculty of Engineering and Industrial Sciences, Swinburne University of Technology, Hawthorn, VIC 3122, Australia

a r t i c l e i n f o

abstract

Article history: Received 9 September 2011 Received in revised form 2 February 2012 Accepted 23 March 2012 Available online 30 March 2012

This paper develops a Kalman smoother to estimate white noise input for systems with random parametric uncertainties in the observation equation. The derived input estimator is optimal in terms of the mean square error (MSE) criterion. Convergence analysis for the derived Kalman smoother is provided, which shows that stability of the Kalman filter cannot guarantee that of the designed fixed-point Kalman smoother. Furthermore, the designed smoothing estimator is applied to the semi-blind deconvolution problem, and an optimal solution is obtained. Numerical examples are given to demonstrate the performance of the proposed method in comparison with two typical deconvolution methods. & 2012 Elsevier B.V. All rights reserved.

Keywords: Kalman filtering Kalman smoothing Semi-blind deconvolution

1. Introduction System input estimation (deconvolution or equalization) is a very important problem in a wide range of application areas, such as medical imaging, telecommunications, seismology, ultrasound and radar signal processing, etc. It has been intensively investigated for a long time and many classic deconvolution methods have been successfully applied, such as least-mean-square method, regularized least-mean-square method, maximum-likelihood method [1], Bayesian inference method [2,3], and so on. The Kalman filter, as a recursive least mean squares method, is a well-known linear unbiased estimator for dynamic systems. It has been successfully applied to estimate inputs (white noise sequences) to moving average (MA), autoregressive (AR), and autoregressive moving average (ARMA) systems [4].

n

Corresponding author. E-mail address: [email protected] (C. Zhang).

0165-1684/$ - see front matter & 2012 Elsevier B.V. All rights reserved. http://dx.doi.org/10.1016/j.sigpro.2012.03.013

White noise system input usually occurs in reflective seismic, ultrasound imaging, and telecommunication systems. Mendel [4] presented optimal white noise input estimators with applications to oil seismic exploration based on Kalman filtering [5]. Deng [6] presented a white noise estimator for dynamic systems with correlated system and measurement noises. Sun et al. [7] presented a white noise deconvolution estimator for time-varying systems by utilizing an information fusion technique. However, a limitation of these works is that system uncertainties were not considered in problem formulations and solutions. There have been many works on estimation of dynamic systems with random parameter matrices [8–10]. Most of these works have temporally independent random parameters and the solutions are suboptimal. Although the recursive estimation algorithm in [11] is optimal, its random parameter matrices are temporally independent and the estimation problem can be converted to the standard Kalman filter. It is noted that, in many practical applications such as multi-target tracking [8] and semi-blind deconvolution, the random parameters

2498

C. Yu et al. / Signal Processing 92 (2012) 2497–2508

are correlated and it is a challenging problem to derive an optimal state or input estimation. In this paper, an optimal smoothing estimator is derived for dynamic systems subject to white noise input and random parametric uncertainties. The system uncertainties are time invariant and thus temporally correlated. Furthermore, convergence analysis is carried out for the smoothing algorithm. It is shown that the convergence condition for Kalman filter cannot guarantee the convergence of the corresponding fixed-point Kalman smoothing due to the presence of temporally correlated system matrices and the system input term in the measurement equation. The derived system input estimator is applied to the semi-blind deconvolution problem which is concerned with two system input–output pairs of the same convolution model: one is known (which can also be considered as the training data) and the other system input needs to be estimated from its corresponding system output. To solve this problem, we adopt the cross relation method [12]. When the measurement noise is present in the convolution model, the state-space representation of the system cross relation contains color noises in the measurement equation and time invariant (temporally correlated) multiplicative noises in system matrices. Such a deconvolution problem is hard to deal with using existing results. Our derived system input estimator can resolve this semi-blind deconvolution problem to present an optimal solution in terms of the mean square error (MSE) performance. The rest of this paper is organized as follows. Section 2 formulates the estimation problem for a dynamic system with temporally correlated random parameters. Section 3 presents an optimal system input estimator. Section 4 presents convergence analysis for the estimator. Section 5 applies the derived system input estimator to the semiblind deconvolution problem. In Section 6, numerical examples of the semi-blind deconvolution problem are given. Finally, conclusions are made in Section 7. 2. Problem formulation Consider the following linear discrete time-varying system with random parametric uncertainties in measurement matrices

yðtÞ ¼

CðtÞ þ

i¼1

! C i ðtÞzi xðtÞ þ DðtÞ þ

Remark 2. The dynamical system described in (1) has the system input term and temporally correlated random parameters in the measurement equation which makes it different from the dynamic systems considered in existing input estimation methods [4,6] and state estimation methods [8–11]. The problem under investigation is stated as: Given the system measurements fyð1Þ,yð2Þ, . . . ,yðLÞg, design an optimal input estimator of the input signal u(t) for t ¼ 1; 2, . . . ,L such that the MSE is minimized. 3. Optimal deconvolution smoother The recursive estimation algorithm derived in this section includes two steps: (1) Compute forward estimation via Kalman filtering and (2) estimate the system input via fixedpoint smoothing. The following lemmas, which are independent of system formulations, will be used for designing the smoothing estimator for the system input in (1). Lemma 1 (Mendel [4]). Denote the set of measurement samples up to time t by YðtÞ ¼ ðyð1Þ,yð2Þ, . . . ,yðtÞÞ0 , where the superscript 0 represents the transpose, the best linear unbiased estimation for the state xðtÞ, t r t can be realized by ~ ^ t9tÞ ¼ E½xðtÞ9YðtÞ ¼ E½xðtÞ9Yðt1Þþ E½xðtÞ9yðt9t1ÞE½xð tÞ xð

ð2Þ ~ ^ where yðt9t1Þ ¼ yðtÞE½yðtÞ9Yðt1Þ ¼ yðtÞyðt9t1Þ is the ^ one-step prediction error of y(t), and yðt9t1Þ the one-step prediction estimate of y(t). Lemma 2 (Mendel [4]). Given Y(t), the best posteriori linear ^ unbiased estimation xðt9tÞ for the state x(t), and its estimation error covariance Px ðt9tÞ can be written as ^ ^ ~ xðt9tÞ ¼ xðt9t1Þ þ J x ðtÞS1 ðtÞyðt9t1Þ P x ðt9tÞ ¼ P x ðt9t1ÞJ x ðtÞS1 ðtÞðJ x ðtÞ0

ð3Þ

~ ~ where Jx ðtÞ ¼ E½xðt9t1Þ y~ 0 ðt9t1Þ, SðtÞ ¼ E½yðt9t1Þ y~ ðt9t1Þ, ~ ~ P x ðt9tÞ ¼ E½xðt9tÞ x~ 0 ðt9tÞ, P x ðt9t1Þ ¼ E½xðt9t1Þ x~ 0 ðt9t1Þ, and ~ ^ xðt9t1Þ ¼ xðtÞxðt9t1Þ. 0

xðt þ 1Þ ¼ AðtÞxðtÞ þ BðtÞuðtÞ m X

Remark 1. It is easy to see that x(t) is independent of u(t). In a similar way, u(t) is independent of fyðtÞ, t otg, Zi and zi are independent of x(t). When the white noise processes u(t) and w(t) are finite-time correlated, it can be coped with using the method in [13].

m X

! Di ðtÞZi uðtÞþ wðtÞ

i¼1

ð1Þ where xðtÞ 2 Rn is the state, yðtÞ 2 Rny is the output, zi , Zi ,i ¼ 1, . . . ,m are random variables, u(t) is the input, w(t) is the measurement noise and AðtÞ,BðtÞ,CðtÞ,C i ðtÞ,DðtÞ,Di ðtÞ are known matrices with proper dimensions. It is noteworthy that the random variables Zi and zi are time invariant, so are temporally correlated random parameters. Assumption 1. The system input u is a white noise process with zero mean and unit variance, and xð0Þ,uðtÞ,fzi gm i ¼ 1, fZi gm ,wðtÞ are mutually independent. Eð z z Þ ¼ d ,Eð Z Z ij i j i jÞ ¼ i¼1 dij ,Eðwðt þ iÞwðt þ jÞÞ ¼ dij , where E denotes the mathematical expectation and dij is the Dirac Delta function.

To derive the optimal input estimator, we introduce auxiliary variables zi ðtÞ ¼ zi xðtÞ,i ¼ f1; 2, . . . ,mg and present the following result. Theorem 1. For the dynamic system (1) with random parameter uncertainties, the optimal input estimation using Kalman smoothing can be calculated recursively as follows: (1) Initial values ^ xð191Þ ¼ Eðxð1ÞÞ,

P x ð191Þ ¼ Varðxð1ÞÞ

^ uð191Þ ¼ Eðuð1ÞÞ,

Pu ð191Þ ¼ Varðuð1ÞÞ

z^ i ð191Þ ¼ Eðzi ð1ÞÞ,

P zi ð191Þ ¼ Varðzi ð1ÞÞ

ð4Þ

C. Yu et al. / Signal Processing 92 (2012) 2497–2508

2499

^

(2) Prediction stage

^ ^ ~ uðt9LÞ ¼ uðt9L1Þ þ Nðt9LÞS1 ðLÞyðL9L1Þ

^ þ 19tÞ ¼ AðtÞxðt9tÞ ^ ^ xðt þ BðtÞuðt9tÞ

ð7Þ

z^ i ðt þ 19tÞ ¼ AðtÞz^ i ðt9tÞ ^ þ19tÞ ¼ Cðt þ1Þxðt ^ þ 19tÞ þ yðt

m X

By introducing the auxiliary variable sequences T1 and fT 2i gm i ¼ 1, the Kalman smoother gain Nðt9LÞ can be computed by

C i ðt þ 1Þz^ i ðt þ 19tÞ

T 1 ðt þ1Þ ¼ B0 ðtÞD0 ðtÞ½AðtÞJx ðtÞ þ BðtÞJ u ðtÞ0

i¼1

P x ðt þ 19tÞ ¼ AðtÞP x ðt9tÞA0 ðtÞ þBðtÞPu ðt9tÞB0 ðtÞ

T 2i ðt þ1Þ ¼ D0 ðtÞJzi ðtÞA0 ðtÞ

P zi ðt þ 19tÞ ¼ AðtÞP zi ðt9tÞA0 ðtÞ þ BðtÞB0 ðtÞ Nðt9t þ 1Þ ¼ T 1 ðt þ 1ÞC 0 ðt þ1Þ þ

Sðt þ 1Þ ¼ Cðt þ 1ÞP x ðt þ 19tÞC 0 ðt þ1Þ m X C i ðt þ1ÞP zi ðt þ 19tÞC 0i ðt þ 1Þ þ m X

T 2i ðt þ1ÞC i 0 ðt þ1Þ

i¼1

T 1 ðt þ2Þ ¼ T 1 ðt þ 1ÞA0 ðt þ 1Þ

i¼1

þ Dðt þ 1ÞD0 ðt þ1Þ þ

m X

Nðt9t þ1Þ½Aðt þ 1ÞJ x ðt þ 1Þ þ Bðt þ 1ÞJ u ðt þ 1Þ0

Di ðt þ 1ÞD0i ðt þ 1Þ þI

i¼1

ð5Þ where I is an identity matrix of proper size. (3) Update stage

T 2i ðt þ2Þ ¼ T 2i ðt þ 1ÞA0 ðt þ 1ÞNðt9t þ 1Þ½Aðt þ1ÞJ zi ðt þ 1Þ0 Nðt9t þ 2Þ ¼ T 1 ðt þ 2ÞC 0 ðt þ2Þ þ

m X

T 2i ðt þ2ÞC i 0 ðt þ2Þ

i¼1

^ ^ ~ xðt9tÞ ¼ xðt9t1Þ þ J x ðtÞyðt9t1Þ J x ðtÞ ¼ Px ðt9t1ÞC 0 ðtÞS1 ðtÞ

^

~ ^ uðt9tÞ ¼ J u ðtÞyðt9t1Þ,

T 1 ðLÞ ¼ T 1 ðL1ÞA0 ðL1ÞNðt9L1Þ½AðL1ÞJ x ðL1Þ

Ju ðtÞ ¼ D0 ðtÞS1 ðtÞ

þBðL1ÞJ u ðL1Þ0

~ z^ i ðt9tÞ ¼ z^ i ðt9t1Þ þ Jzi ðtÞyðt9t1Þ J zi ðtÞ ¼ P zi ðt9t1ÞC 0i ðtÞS1 ðtÞ

T 2i ðLÞ ¼ T 2i ðL1ÞA0 ðL1ÞNðt9L1Þ½AðL1ÞJ zi ðL1Þ0

P x ðt9tÞ ¼ P x ðt9t1ÞP x ðt9t1ÞC 0 ðtÞS1 ðtÞCðtÞP x ðt9t1Þ Nðt9LÞ ¼ T 1 ðLÞC 0 ðLÞ þ

P u ðt9tÞ ¼ ID0 ðtÞS1 ðtÞDðtÞ

m X

T 2i ðLÞC i 0 ðLÞ

ð8Þ

i¼1

P zi ðt9tÞ ¼ P zi ðt9t1ÞP zi ðt9t1ÞC 0i ðtÞS1 ðtÞC i ðtÞPzi ðt9t1Þ

Proof. See Appendix.

&

ð6Þ (4) Fixed-point smoothing ~ þ19tÞ ^ ^ uðt9t þ1Þ ¼ uðt9tÞ þNðt9t þ 1ÞS1 ðt þ1Þyðt ^ ^ ~ þ 29t þ 1Þ uðt9t þ2Þ ¼ uðt9t þ 1Þ þ Nðt9t þ 2ÞS1 ðt þ 2Þyðt

Fig. 1 shows the flow diagram of the derived smoothing estimator for the system input. In the Kalman filtering process, the prediction stage and the update stage are alternatively carried out at each discrete time t. In addition, the fixed-point Kalman smoothing is performed after the Kalman filtering.

Fig. 1. The algorithm flow diagram of the smoothing estimator for system input.

2500

C. Yu et al. / Signal Processing 92 (2012) 2497–2508

0

4. Convergence analysis

Li2 ðPÞ ¼ ðAþ K i3 C i ÞP i3 ðA þ K i3 C i Þ0 þ K i3 @CP1 C 0 þ

This section is on convergence analysis of the Kalman smoother with time-invariant parameters. From the estimation error covariance evolution equations in (5) and (6), the following Riccati iterative functions can be derived: P x ðt þ 1Þ ¼ f 1 ðPx ðtÞ,fP zj ðtÞgm j ¼ 1Þ

1 X C j P j3 C 0j AðK i3 Þ0 jai

ð12Þ i fP 1 ,fP 3 gm i ¼ 1 g 4 0,

Suppose there exits a matrix set P ¼ which implies that each element in the set P is positive definite, such i that P 1 4L1 ðPÞ, P 3 4Li2 ðPÞ. We can have,

¼ AðPx ðtÞP x ðtÞC 0 S1 ðtÞCPx ðtÞÞA0 þBðID0 S1 ðtÞDÞB0 i P zi ðt þ 1Þ ¼ f 2 ðP x ðtÞ,fPzj ðtÞgm j ¼ 1Þ

¼ AðP zi ðtÞPzi ðtÞC i 0 S1 ðtÞC i Pzi ðtÞÞA0 þ BB0 , i 2 f1; 2, . . . ,mg x

ðtÞ,Pzi ðtÞ

x

ð9Þ

ðt9t1Þ,P zi ðt9t1Þ

are used to represent P where P P 0 z 0 for simplicity, SðtÞ ¼ CP x ðtÞC 0 þ m i ¼ 1 C i P i ðtÞC i þ DD þ R, where Pm i R ¼ i ¼ 1 Di D0i þ I, and f 1 ðÞ,f 2 ðÞ are functions with respect to variables Px ðtÞ,fP zi ðtÞgm j ¼ 1. The following lemma is on optimality of the Kalman filtering. i m Lemma 3. Denote K ¼ fK 1 ,K 2 ,fK i3 gm i ¼ 1 g, P ¼ fP 1 ,fP 3 gi ¼ 1 g, and ! m X f1 ðK,PÞ ¼ ðA þ K 1 CÞP1 ðAþ K 1 CÞ0 þ K 1 C i Pi3 C i 0 þ DD0 þ R K 01

n n om o Z 0, limk-1 ðL1 Þk ðWÞ ¼ 0, (a) For all W ¼ W 1 , W i3 i¼1 limk-1 ðLi2 Þk ðWÞ ¼ 0. i m (b) Let V ¼ fV 1 ,fV i3 gm i ¼ 1 g Z 0, PðkÞ ¼ fP 1 ðkÞ,fP 3 ðkÞgi ¼ 1 g Z0, and consider the linear system P1 ðk þ 1Þ ¼ L1 ðPðkÞÞ þV 1 Pi3 ðk þ 1Þ ¼ Li2 ðPðkÞÞ þV i3 ,

þðB þ K 2 DÞðB þ K 2 DÞ0 þ K 2 CP 1 C 0 þ

ð13Þ

Proof. (a) We choose a scalar 0 rr o 1 such that rP 1 4 i

L1 ðPÞ,rP 3 4 Li2 ðPÞ and choose another scalar 0 r l such that W r lP, which means the inequality holds for correi

sponding element pairs: W 1 r lP 1 ,W i3 r lP 3 . Then 0 rL1 ðWÞ r lL1 ðPÞ o lrP 1

i¼1

m X

i ¼ 1; 2, . . . ,m

Then, the sequence PðkÞ ¼ fP1 ðkÞ,fP i3 ðkÞgm i ¼ 1 g is bounded.

i

0 rLi2 ðWÞ r lðLi2 ÞðPÞ o lrP 3

! C i P i3 C i 0 þ R K 02

i¼1

0 rðL1 Þ2 ðWÞ r lðL1 Þ2 ðPÞ i

2 o lrL1 ðP 1 ,fP 3 gm i ¼ 1 Þ ¼ lrL1 ðPÞ o lr P 1

fi2 ðK,PÞ ¼ ðA þK i3 C i ÞP i3 ðA þ K i3 C i Þ0 0

þK i3 @CP 1 C 0 þ

1 X j 0 0 C j P 3 C j þ DD þ RAðK i3 Þ0 þ BB0

i

ð10Þ

0 rðLi2 Þ2 ðWÞ r lrLi2 ðPÞ o lr 2 P 3 ^

jai

Then, f 1 ðPÞ ¼ f1 ðK~ ,PÞ r f1 ðK,PÞ,

0 rðL1 Þk ðWÞ o lr k P 1 8K i

i

i

i

f 2 ðPÞ ¼ f2 ðK~ ,PÞ r f2 ðK,PÞ,

8K

ð11Þ

where K~ is the Kalman gain defined by K~ ¼ fK x ,K u ,fK zi gm i ¼ 1 g, K x ¼ AP 1 C 0 S1 , K u ¼ BD0 S1 , K zi ¼ APi3 C 0i S1 , S ¼ CP 1 C 0 þ Pm i 0 0 i ¼ 1 C i P 3 C i þDD þ R. i

Proof. Functions f1 ðK,PÞ and f2 ðK,PÞ are quadratic and convex with respect to the variable K. Hence, the minimum solution can be determined from the corresponding first-order optimality condition. The first order optimality condition for the function f1 ðK,PÞ is ! m X @f1 ðK,PÞ i 0 0 0 ¼ 2CP ðA þ K CÞ þ 2 C P C þ DD þ R K 01 ¼ 0: 1 1 i 3 i @K 01 i¼1 Obviously, the optimal solution is K n1 ¼ K x ¼ AP 1 C 0 S1 , which is exactly the Kalman gain. The proof for other results is analogous. & Lemma 4. Denote P

0 rðLi2 Þk ðWÞ o lr k P 3 Obviously, the conclusion can be obtained as k-1. (b) We choose a scalar 0 r l1 such that V r l1 P, and choose another scalar 0 r l2 such that Pð0Þ r l2 P. Then the bound for each iteration can be obtained as follows: P 1 ð1Þ ¼ L1 ðPð0ÞÞ þV 1 r l2 L1 ðPÞ þV 1 o l2 rP 1 þ V 1 i

P i3 ð1Þ ¼ Li2 ðPð0ÞÞ þV i3 o l2 rP 3 þ V i3 P 1 ð2Þ ¼ L1 ðPð1ÞÞ þ V 1 ¼ L1 ðP 1 ð1Þ,fPi3 ð1Þgm i ¼ 1Þ þ V 1 i m ¼ L1 ðL1 ðPð0ÞÞ,fLi2 ðPð0ÞÞgm i ¼ 1 Þ þ L1 ðV 1 ,fV 3 gi ¼ 1 Þ

þ V 1 o l2 r 2 P 1 þ l1 rP 1 þ V 1 i

^

¼ fP 1 ,fP i3 gm i ¼ 1 g,

and define linear operators ! m X i 0 0 C i P 3 C i K 01 L1 ðPÞ ¼ ðA þ K 1 CÞP 1 ðA þ K 1 CÞ þ K 1 þ K 2 CP 1 C 0 þ

m X i¼1

!

P 1 ðkÞ o l2 r k P 1 þ l1

k1 X

rj P 1 þ V 1

j¼0

i¼1

C i P i3 C 0i K 02

i

P i3 ð2Þ ¼ Li2 ðPð1ÞÞ þ V i3 o l2 r 2 P 3 þ l1 rP 3 þ V i3

i

P i3 ðkÞ o l2 r k P 3 þ l1

k1 X j¼0

i

r j P 3 þ V i3

C. Yu et al. / Signal Processing 92 (2012) 2497–2508

Since 0 r r o1, we can establish that the sequences P1 ðkÞ and fP i3 ðkÞgm i ¼ 1 are bounded with their upper bounds dependent on their initial values. & Based on the above lemmas, a sufficient condition for the convergence of the Riccati iterative equation in (9) is shown in the following theorem. Theorem 2. Suppose there exist positive-definite matrices i P~ 1 ,fP~ 3 gm i ¼ 1 , such that 0 ~ 1

P~ 1 4 AðP~ 1 P~ 1 C S

0 ~ 1

C P~ 1 ÞA0 þBðID S

1 i i i i P~ 3 4 AðP~ 3 P~ 3 C i 0 S~ C i P~ 3 ÞA0 þBB0 ,

DÞB0

i 2 f1; 2, . . . ,mg

ð14Þ

0 ~i 0 i ¼ 1 C i P 3 C i þDD þ R.

where S~ ¼ C P~ 1 C þ Then, the Riccati iteration in (9) converges to a fixed point: limt-1 P x ðtÞ ¼ Px and limt-1 P zi ðtÞ ¼ P zi . Proof. See Appendix

&.

The following theorem provides numerical computation methods to compute the fixed point of the Riccati equation (9). Theorem 3. Suppose there exists a unique limiting point fP x , fPzi gm i ¼ 1 g for the iterative Riccati equation (9), then the fixed point can be computed by the following two ways: (a) Compute the Riccati iterations in (9) until it arrives at a steady solution since limt-1 P x ðtÞ ¼ P x and limt-1 P zi ðtÞ ¼ P zi . (b) Solve the following semi-definite convex programming problem: ! m X z x max tr Pi þ P z x

0

0

0

P x aAðP x Px C S1 C P x ÞA þ BðID S1 DÞB 0

0

0

P zi aAðPzi P zi C i S1 C i P zi ÞA þ BB , x

z P^ i

Since P^ and problem, then

0

i 2 f1; 2, . . . ,mg

are feasible solutions of the optimization

0

0

0

P zi oAðP zi P zi C i S1 C i Pzi ÞA þBB ,

i 2 f1; 2, . . . ,mg

ð16Þ

P ^z ^x The above inequalities imply that trð m i ¼ 1 P i þ P Þ otr P z x ^x ð m i ¼ 1 P i þP Þ, which contradicts the hypothesis that P z x z and P^ are optimal solutions. Therefore, P^ and P^ should i

be the solution of (9) and this concludes the proof of the theorem. & Based on the stability of Kalman filtering, the steady estimation error covariances Px ,fP zi gm i ¼ 1 , the steady matrix gains J u ,Jx ,fJzi gm i ¼ 1 in (6) and the steady Kalman gains K u ,K x ,fK zi gm can be determined. The relationship i¼1 between steady matrix gains and steady Kalman gains can be described by K x ¼ AJ x ,

K u ¼ BJ u ,

K zi ¼ AJ zi

ð17Þ

The covariance of the input estimation error of the Kalman smoothing process based on steady Kalman filtering can be computed as follows: P u ðt9tÞ ¼ ID0 S1 D P u ðt9t þ 1Þ ¼ P u ðt9tÞNðt9t þ1ÞS1 N0 ðt9t þ1Þ ^

i¼1

2 6 4

subject to "

i

the optimization problem, but not the solution of (9), i.e.

i

Remark 3. Theorem 2 provides a sufficient condition for the convergence of the Riccati evolution equation (9). In addition, the given condition is easy to check using the LMI toolbox in MATLAB.

P i ,P Z 0

solution of the semi-definite programming problem is the x z fixed point of (9). Suppose, it is not, i.e. P^ and P^ solves

x 0 0 0 0 P^ oAðP x Px C S1 C P x ÞA þ BðID S1 DÞB

Pm

0

2501

APx A0 P x þ BB0

APx C 0

x 0

CP A

S

DB0

0

AP zi A0 P zi þ BB0

AP zi C i 0

C i P zi A0

S

BD0

3

7 0 5Z0 S

#

ð15Þ Z0,

i ¼ f1; 2, . . . ,mg

ð18Þ P u ðt9t þ lÞ ¼ P u ðt9t þ l1ÞNðt9t þ lÞS1 N0 ðt9t þ lÞ P z 0 0 where S ¼ CP x C 0 þ m i ¼ 1 C i P i C i þ DD þ R, and Nðt9t þ lÞ can be represented in terms of T 1 ðt þ lÞ and T 2i ðt þ lÞ: Nðt9t þ lÞ ¼ T 1 ðt þ lÞC 0 þ

m X

ðT 2i ðt þlÞÞC 0i

i¼1

2

where S ¼ CP x C 0 þ

Pm

i¼1

C i Pzi C i 0 þDD0 þ R.

Proof. Here, we only prove that the fixed point of the Riccati equation in (9) can be computed by solving the semi-definite programming problem. The optimization problem can be obtained using the Schur complement decomposition on the coupled Riccati equations (9) 0

0

0

P x r AðP x Px C S1 C P x ÞA þBðID S1 DÞB 0 0 0 P zi rAðP zi P zi C i S1 C i Pzi ÞA þBB , x

0

i 2 f1; 2, . . . ,mg

Clearly, the solutions of P and Pzi belong to the feasible set of the optimization problem. We now show that the

3 C0 6 0 7 6 C1 7 7 ¼ ½T 1 ðt þ lÞ T 21 ðt þ lÞ    T 2m ðt þ lÞ6 6 ^ 7 ¼ Tðt þ lÞC 4 5 0 Cm

where Tðt þ lÞ ¼ ½T 1 ðt þlÞ T 21 ðt þ lÞ    T 2m ðt þlÞ and C ¼ ½C C 1    C m 0 Furthermore, T 1 ðt þ lÞ and T 2i ðt þ lÞ follow the following evolution equation: Tðt þ lÞ ¼ ½T 1 ðt þlÞ T 21 ðt þ lÞ    T 2m ðt þ lÞ ¼ ½T 1 ðt þ l1Þ T 21 ðt þ l1Þ    T 2m ðt þl1ÞM ¼ Tðt þl1ÞM ¼ Tðt þ 1ÞMl1

ð19Þ

2502

C. Yu et al. / Signal Processing 92 (2012) 2497–2508

where M is a constant matrix which is defined in terms of steady matrix gains J u ,J x ,fJzi gm i ¼ 1: 2 3 ðAJ z1 CÞ0  ðAJ zm CÞ0 ðAðAJ x þ BJ u ÞCÞ0 6 x z z u 0 0 0 7 6 ððAJ þ BJ ÞC 1 Þ ðAAJ 1 C 1 Þ    ðAJm C 1 Þ 7 7 M¼6 6 7 ^ ^ ^ 4 5 x z z u 0 0 0 ðAJ1 C m Þ    ðAAJ m C m Þ ððAJ þ BJ ÞC m Þ 2

ðA þ K x C þ K u CÞ0

6 6 ðK x C 1 þ K u C 1 Þ0 ¼6 6 ^ 4 ðK x C m þ K u C m Þ0

3

ðK z1 CÞ0



ðK zm CÞ0

ðA þ K z1 C 1 Þ0



ðK zm C 1 Þ0



ðA þ K zm C m Þ0

^

^

ðK z1 C m Þ0

7 7 7 7 5

be simplified as a single term A þ K x C þK u C and the steady Riccati equation (9) can be written as P x ¼ ðA þK x CÞP x ðAþ K x CÞ0 þ ðK u CÞP x ðK u CÞ0 þ K x ðDD0 þIÞðK x Þ0 þðB þ K u DÞðB þ K u DÞ0 þ K u ðK u Þ0

ð21Þ

The above equation can be written as a vector form: VecðPx Þ ¼ ½ðA þ K x CÞ  ðAþ K x CÞ þK u C  K u CVecðP x Þ þ    where the terms irrelevant to Px are omitted, VecðÞ is a vectorizing operator, and  denotes the Kronecker product. From the above equation, we can find that the Kalman filtering is stable if

rð½ðA þ K x CÞ  ðAþ K x CÞ þK u C  K u CÞ o 1 Theorem 4. Under the assumption that the Kalman filtering is stable, i.e. Riccati equation in (9) converges to a fixed point, then the corresponding fixed-point smoothing is stable if the spectral radius of the matrix M is less than one, i.e. rðMÞ o1. Proof. The estimation error covariance of the system input can be derived iteratively as follows: P u ðt9t þ lÞ ¼ Pu ðt9t þl1ÞNðt9t þ lÞS1 N0 ðt9t þ lÞ ¼ P u ðt9tÞ

l X

Nðt9t þ jÞS1 N0 ðt9t þ jÞ

j¼1 u

¼ P ðt9tÞ

l X j¼1

Tðt þ jÞCS1 C0 T 0 ðt þ jÞ 0

u

¼ P ðt9tÞTðt þ 1Þ@

l X

1 M

j1

1 0

CS

C ðM

j1 0 A 0

Þ

T ðt þ 1Þ

j¼1

ð20Þ In order to prove that the smoothing process converges as l-1, we need to prove that the second term of the above equation is convergent. Here, an auxiliary Lyapunov equation is constructed: M ¼ MMM 0 þ CS1 C0 When rðMÞ o1, the above equation has a unique solution of M which is given by M¼

1 X

M j1 CS1 C0 ðM j1 Þ0

j¼1

Therefore, when rðMÞ o1, the smoothing process is convergent, and the estimation error covariance for the system input will converge to P u ðt91Þ ¼ P u ðt9tÞTðt þ 1ÞMTðt þ1Þ0 . Thus, the proof is completed. & Corollary 1. Suppose there exist no matrix uncertainties and no system input in the measurement equation, i.e. C i ¼ 0, D ¼ 0,Di ¼ 0, the stability of Kalman filtering can result in the stability of Kalman smoothing [4]. When C i ¼ 0,D ¼ 0,Di ¼ 0, the Riccati equation (9) can be simplified as a standard Riccati equation and the matrix M is simplified as a single term A þ K x C. As is well known, when ½C,A is detectable, rðA þ K x CÞ o 1 holds. Thus, the stability of Kalman filtering can result in the stability of Kalman smoothing in this case. In general, however, we cannot have this conclusion even for the case of C i ¼ 0,Di ¼ 0. When C i ¼ 0,Di ¼ 0, the matrix M can

However, this condition cannot result in rðMÞ ¼ rðA þ K x C þ K u CÞ o 1. Thus, the stability of Kalman filtering cannot guarantee the stability of Kalman smoothing in this case. 5. Application of optimal smoother to semi-blind deconvolution The semi-blind deconvolution problem concerns two input–output pairs in the same convolution model. One input–output pair is known and it requires to estimate the other system input from corresponding system observations. Two typical deconvolution methods are usually applied to solve such kind of semi-blind deconvolution problem: (1) the cross relation method which does not need to estimate the impulse response explicitly [12,14] and (2) the cascade deconvolution method which first estimates the impulse response function using the known input–output pair and then estimates the another system input using its corresponding output. The cross relation method can be carried out without involving statistical properties of the impulse response function, which may lead to great estimation error. The cascade deconvolution method can achieve optimal estimations at each step, but it cannot result in a global optimal solution even though it is optimal at each step. In this section, a smoothing estimator is derived for the semi-blind deconvolution problem to achieve an optimal solution. Two system input–output pairs can be represented in one convolution model as follows: y1 ¼ x1 nh þ n1 y2 ¼ x2 nh þ n2

ð22Þ

where x1 2 Rn ,y1 2 Rðn þ l1Þ is a known system input–output pair, h 2 Rl is an impulse response function, and n1 2 Rðn þ l1Þ is a Gaussian white noise sequence with unit variance. In the second equation, x2 ,y2 ,n2 of proper size is another set of system input, system output and measurement noise. The semi-blind deconvolution problem concerned in this paper can be stated as: Assuming x1 and y1 are known, estimate the system input x2 from the observation y2 . It is different from total blind channel estimation problem [14] that the coprime condition for all channel functions is unnecessary for the semi-blind deconvolution problem. To solve the semi-blind deconvolution problem,

C. Yu et al. / Signal Processing 92 (2012) 2497–2508

2503

a cross relation equation can be constructed as follows:

D ¼ ½y1 ð0Þ11 ,

ðy1 n1 Þnx2 ¼ ðy2 n2 Þnx1

Combining the state-space equations in (25) and (26), an augmented system can be constructed: " # " # " #" #" # xðt þ1Þ xðtÞ uðtÞ A 0 B 0 ¼ þ ~ þ1Þ ~ xðt xðtÞ vðtÞ 0 A0 0 B0 |fflfflfflfflfflfflffl{zfflfflfflfflfflfflffl} |fflfflfflfflfflffl{zfflfflfflfflfflffl} |fflfflfflffl{zfflfflfflffl} |fflfflfflfflfflffl{zfflfflfflfflfflffl} |fflfflfflffl{zfflfflfflffl}

ð23Þ

It can be rewritten as y2 nx1 ¼ ðy1 n1 Þnx2 þx1 nn2 i.e. n1 X

nX þ l2

y2 ðttÞx1 ðtÞ ¼

t¼0

xðt þ 1Þ

x2 ðtaÞðy1 ðaÞn1 ðaÞÞ

þ

n2 ðtbÞx1 ðbÞ

ð24Þ

b¼0

C

where sequences x1 , y1 and y2 are known, n1 and n2 are Gaussian white noises, and x2 is a sequence of Gaussian white random variables to be estimated. It is observed that the term n1 ðaÞ does not depend on the time t, so the parametric uncertainty is fixed with respect to the time t. Furthermore, the second term on the right hand side is a color noise process driven by the noise sequence n2 . The convolution equation in (24) can be formulated into a state-space system xðt þ 1Þ ¼ AðtÞxðtÞ þ BuðtÞ yðtÞ ¼

C

nX þ l2

! C i n1 ðn þl1iÞ xðtÞ

i¼1

þðD þ D1 n1 ð0ÞÞuðtÞ þ wðtÞ

ð25Þ

where w(t) is a color noise generated from the following system: ~ þ B0 vðtÞ ~ þ 1Þ ¼ A0 xðtÞ xðt

0 0 0  2 0 1 0  60 0 1  6 6 A0 ¼ 6 6^ 6 40 0 0  0

0

0



2 3 0 607 6 7 6 7 7 B¼6 , 6^7 6 7 405 1 ðn þ l2Þ1

ð26Þ

yðtÞ ¼

nX þ l2 i¼1

uðtÞ

B

1 " # C xðtÞ ½C i 0 n1 ðn þ l1iÞA ~ |fflffl{zfflffl} xðtÞ Ci |fflfflfflffl{zfflfflfflffl}

0

xðtÞ

1

# B C" B C uðtÞ C þB ½D D  þ½D 0 n ð0Þ 0 1 1 B|fflfflffl{zfflfflffl} |fflfflffl{zfflfflffl} C vðtÞ @ A |fflfflfflffl{zfflfflfflffl} D

ð27Þ

D1

uðtÞ

It is noted that the multiplicative noise terms fn1 ðn þ l1iÞgni ¼þ 0l1 are mutually independent with each other and are invariant with respect to the time t. In addition, the system input u(t) is a zero-mean Gaussian white process with an identity covariance matrix. Corollary 2. The semi-deconvolution problem (22) can be properly resolved by the input estimator described (5)–(7) after it is transformed to the state-space system (27) which is a special case of the dynamic system (1). 6. Simulation results

~ þ D0 vðtÞ wðtÞ ¼ C 0 xðtÞ where uðtÞ ¼ x2 ðtÞ, and 2 0 1 0  60 0 1  6 6 A¼6 6^ 6 40 0 0 

B yðtÞ ¼ @½C C 0  þ |fflfflffl{zfflfflffl}

D0 ¼ ½x1 ð0Þ11

xðtÞ

A

0

a¼0 n1 X

D1 ¼ ½111 ,

Pn1

t ¼ 0 y2 ðttÞx1 ðtÞ,

3 0 07 7 7 7 7 7 15 0 ðn þ l2Þðn þ l2Þ 3 0 07 7 7 7 7 7 15 0 ðn1Þðn1Þ 2 3 0 607 6 7 6 7 7 B0 ¼ 6 6^7 6 7 405 1 ðn1Þ1

vðtÞ ¼ n2 ðtÞ,

We now present simulation results to illustrate the performance of our proposed deconvolution algorithm and compare it with the cross relation method [12,14,15] and a cascade deconvolution method which estimates the impulse response function first followed by estimating the other system input. In the cascade deconvolution method, we consider the impulse response sequence as a sequence of white Gaussian variables with zero mean and unit variance (actually the impulse response sequence is generated by the MATLAB command RANDN in the simulation). The first convolution equation in (22) can be formulated as a matrix–vector form y1 ¼ X1 hþ n1 2

3

x1 ð0Þ

6 x ð1Þ 6 1 6 6 ^ 6 X1 ¼ 6 6 x1 ðn1Þ 6 6 6 6 4

x1 ð0Þ x1 ð1Þ

&

^

&

x1 ðn1Þ

C ¼ ½y1 ðn þl2Þ y1 ðn þ l3Þ . . . y1 ð1Þ1ðn þ l2Þ C i ¼ ½0   ffl0} 1 0    01ðn þ l2Þ |fflfflffl{zfflffl i1

C 0 ¼ ½x1 ðn1Þ x1 ðn2Þ    x1 ð1Þ1ðn1Þ

&

7 7 7 7 7 x1 ð0Þ 7 7 7 x1 ð1Þ 7 7 7 ^ 5 x1 ðn1Þ

ð28Þ

ðn þ l1Þl

where X1 is a Toeplitz matrix of full column rank since x1 is not a zero vector. Then, the maximum a posteriori (MAP) estimate of the impulse response is h^ ¼ ðX01 X1 þ IÞ1 X01 y1

ð29Þ

C. Yu et al. / Signal Processing 92 (2012) 2497–2508

In the same way, x^ 2 can be easily estimated based on the available information of h^ and y2 . In the cross relation method (without using any knowledge of the impulse response h), the cross relation equation in the absence of noise can be written as y1 nx2 ¼ y2 nx1 or

Y1 x2 ¼ Y2 x1

ð30Þ

where Y1 and Y2 are Toeplitz matrices of the same form as X1 . Hence, the vector x^ 2 can be computed as follows: x^ 2 ¼ ðY01 Y1 Þ1 Y01 Y2 x1

ð31Þ

Comparing with the cascade method, the cross relation method is less affected by the dimension of the training sequence y1 , because it does not need to estimate the impulse response explicitly. To have a fair comparison with the above two methods, the impulse response sequence h of fixed length is generated by the command RANDN on MATLAB. For the known system input–output pair, its system input x1 is randomly generated and its system output y1 is produced by taking the convolution of the impulse response sequence and the system input followed by adding standard Gaussian white noise n1 . Another set of system input x2 and system output y2 can be generated in the same way. Here, our task is to estimate the sequence x2 when the sequences x1 , y1 and y2 are available. From the definition of convolution in (24), it can be observed that the current system input contributes to the values of future output samples of the same length as that of the impulse response sequence. Thus, in order to reduce the computational amount in a reasonable way, we apply the fixed-lag smoothing method for the input estimation where the fixed-lag is determined by the length of the impulse response sequence. When the length of h is l, we can find that the sequence ðy1 n1 Þ in (24) is of length n þ l1. Therefore, n þ l1 future output samples are utilized to estimate the current system input. Using the fixed-point smoothing algorithm as shown in Theorem 1, the whole input sequence in (24) can be estimated iteratively and the computational load is not heavy in each iteration. On the other hand, the cascade deconvolution method and the cross relation method compute the whole input sequence in one step using matrix–vector operations, and their computational load will be quite heavy (especially for the inverse operation of a large matrix) when the input sequence is too long. To evaluate performances of different deconvolution methods, the root mean square error (RMSE) criterion is adopted. Here, we define RMSE value at time t as vffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi u N u1 X ðiÞ RMSEðtÞ ¼ t ½uðtÞu^ ðtÞ2 ð32Þ Ni¼1 ðiÞ

where N is the total number of samples at time t, u^ ðtÞ denotes the ith sample of the estimated system input at time t. Next, three sets of simulation results will be shown to illustrate the validity and the reliability of our developed method: (1) the performance in terms of RMSE at

different time steps; (2) the performance as a function of the mean energy of the impulse response at a fixed time step; (3) the performance as a function of the length of the impulse response at a fixed time step. The first simulation result shows the overall performance, while the last two simulation results show the influence of the impulse response on the deconvolution performance. In Fig. 2, RMSE curves of three different methods at different time steps are computed by averaging N ¼600 Monte Carlo trials. In this simulation, the length of the impulse response is set to l¼ 7, the length of x1 is 6, and the standard deviation of the impulse response is set to two. Since the impulse response is generated as a Gaussian white noise, the standard deviation can be used to measure the average energy. From Fig. 2, it can be found that our proposed deconvolution algorithm outperforms the cascade deconvolution method and the noise-free cross relation method. The cascade deconvolution method is optimal in each step but is not globally optimal, so its performance is worse than that of our proposed method. In addition, the worse performance of the cross relation method is obvious since it does not take the noise into account. In Fig. 3, performances (in terms of RMSE) of the different algorithms with respect to different conditions of the impulse response are shown at a fixed time step. The concerned conditions of the impulse response include the length and the mean energy, where the mean energy is measured by the standard deviation. Increasing the mean energy will lead to the increase of the signal noise ratio, and increasing the length of the impulse response will increase the dimension of training samples. From Fig. 3, it can be found that the deconvolution performance is improved when the length or the mean energy of the impulse response increases. On the right hand side of this figure, the performance of the cross relation algorithm is better than that of the cascade method when the mean energy of the impulse response is high enough. It is

0.85

Proposed method Cascade method Cross relation method

0.8 0.75 0.7 RMSE

2504

0.65 0.6 0.55 0.5 0.45 0.4 0

10

20

30

40

50

60

70

80

90

Time Step

Fig. 2. Comparison of performances of different deconvolution methods in terms of RMSE: cascade deconvolution method (dot-dash line), noisefree cross relation method (dash curve), and our proposed method (solid curve). The horizontal coordinate axis represents the sample index and the vertical coordinate axis represents the RMSE value.

C. Yu et al. / Signal Processing 92 (2012) 2497–2508

2505

1.6

0.75 Proposed method Cascade method Cross relation method

0.7

Proposed method Cascade method Cross relation method

1.4

0.65 1.2

0.55

RMSE

RMSE

0.6

0.5

1

0.8

0.45 0.6 0.4 0.4

0.35

0.2 7

8

9

10

11

12

13

14

1

Length of the Impulse Response

2

3

4

5

6

7

8

Standard Deviation of the Impulse Response

Fig. 3. Performances of different deconvolution methods with respect to the length (Left) and the mean energy (Right) of the impulse response: cascade deconvolution method (dot-dash line), noise-free cross relation method (dash curve), and our proposed method (solid curve). For the left figure, the standard deviation of the impulse response is two; for the right figure, the length of the impulse response is seven.

because the performance of the cascade method is constrained by the dimension of the training samples in the model (28); however, the cross relation method estimates the system input x2 without explicitly estimating the impulse response as did the cascade method, hence it performs better than the cascade method when the signal noise ratio is high enough.

and y(t) are uncorrelated. We can use Lemmas 1 and 2 to obtain

7. Conclusion

~ þ19tÞ ¼ xðt þ 1Þxðt ^ þ19tÞ ¼ AðtÞxðt9tÞ ~ ~ xðt þ BðtÞuðt9tÞ

In this paper, an optimal system input estimator in terms of the MSE has been derived for systems with random parametric uncertainties in the observation equation and its stability condition has been given. It has been shown that for systems with multiplicative noise terms and the input term in the measurement equation, the stability of the Kalman filter cannot guarantee that of the designed Kalman smoother. Furthermore, the derived input estimator has been applied to solve the semi-blind deconvolution problem whose state-space system contains temporally correlated multiplicative noise and color additive noise in its measurement equation. The optimal solution has been obtained for the semi-blind deconvolution problem and simulation result has been presented to demonstrate the improved system performance as compared to two typical deconvolution methods.

P x ðt þ 19tÞ ¼ AðtÞP x ðt9tÞA0 ðtÞ þBðtÞPu ðt9tÞB0 ðtÞ

Eðzi uðtÞ9YðtÞÞ ¼ Eðzi uðtÞ9yðtÞÞ ¼ Eðzi uðtÞÞ þ E½zi uðtÞy0 ðtÞ½VarðyðtÞÞ1 yðtÞ ¼ 0 In the prediction stage, we can get ^ þ19tÞ ¼ Eðxðt þ 1Þ9YðtÞÞ ¼ AðtÞxðt9tÞ þ BðtÞuðt9tÞ xðt

z^ i ðt þ 19tÞ ¼ AðtÞz^ i ðt9tÞ þ BðtÞEðzi uðtÞ9YðtÞÞ ¼ AðtÞz^ i ðt9tÞ z~ i ðt þ 19tÞ ¼ AðtÞz~ i ðt9tÞ þ BðtÞzi uðtÞ P zi ðt þ19tÞ ¼ AðtÞP zi ðt9tÞA0 ðtÞ þ BðtÞB0 ðtÞ

ðA:2Þ

where the last equation is based on that z~ i ðt9tÞ and zi uðtÞ are uncorrelated: Efz~ i ðt9tÞu0 ðtÞz0i g ¼ Ef½zi ðtÞz^ i ðt9tÞu0 ðtÞz0i g ¼ 0 In addition, the prediction of system output can be deduced as follows: ^ þ 19tÞ ¼ Efyðt þ 1Þ9YðtÞg yðt ^ þ 19tÞ þ ¼ Cðt þ1Þxðt

m X

C i ðt þ 1Þz^ i ðt þ 19tÞ

i¼1

Appendix A. Proof of Theorem 1 According to Assumption 1, the following statistical properties can be achieved: EðuðtÞÞ ¼ 0, 0

Eðzi Þ ¼ 0,

E½uðtÞx ðtÞ ¼ 0, E½zi xðtÞ ¼ 0,

EðxðtÞÞ ¼ 0,

E½zi uðtÞ ¼ 0, 0

E½yðtÞu0 ðtÞzi  ¼ 0

~ þ 19tÞ ¼ Cðt þ1Þxðt ~ þ 19tÞ þ yðt " þ Dðt þ 1Þ þ

EðyðtÞÞ ¼ 0

m X i¼1

m X

C i ðt þ 1Þz~ i ðt þ19tÞ #

Di ðt þ1ÞZi uðt þ1Þ þwðt þ 1Þ

i¼1

E½zi yðtÞ ¼ 0

~ þ 19tÞy~ 0 ðt þ 19tÞg ¼ Cðt þ 1ÞPx ðt þ19tÞC 0 ðt þ1Þ Sðt þ 1Þ ¼ Efyðt ðA:1Þ

From the above properties, it can be verified that zi uðtÞ

þ

m X i¼1

C i ðt þ1ÞP zi ðt þ 19tÞC 0i ðt þ 1Þ þ Dðt þ 1ÞDðt þ 1Þ0

2506

C. Yu et al. / Signal Processing 92 (2012) 2497–2508

þ

m X

"

(

Di ðt þ1ÞD0i ðt þ 1Þ þ I

ðA:3Þ

~ þ 29t þ 1Þþ ¼ E uðtÞ Cðt þ 2Þxðt

#0 ) C i ðt þ2Þz~ i ðt þ29t þ 1Þ

i¼1

i¼1

(

where the last equation is based on that:

m X

i¼1

#0 ) C i ðt þ2ÞAðt þ 1Þz~ i ðt þ19t þ1Þ "

~ þ 19tÞ ¼ E uðtÞ Cðt þ 2ÞAðt þ1Þxðt

In the measurement update stage, the recursive equations in (6) can be obtained using Lemma 2. The difficulty to derive these equations is to determine the matrix gains. One example is shown as follows:

~ ~ ¼ E xðt9t1Þ CðtÞxðt9t1Þ þ

þ

m X

(

~ þ 19tÞz~ 0i ðt þ 19tÞg ¼ Ef½xðt þ 1Þ Efxðt ^ þ 19tÞx0 ðt þ 1Þzi 0 g ¼ 0: xðt

~ J x ðtÞ ¼ Efxðt9t1Þ y~ 0 ðt9t1Þg " (

"

~ þ 19t þ 1Þ þBðt þ 1Þuðt ~ þ19t þ1ÞÞ ¼ E uðtÞ Cðt þ 2ÞðAðt þ 1Þxðt

~ þ19tÞ and (1) uðt þ1Þ is uncorrelated with both xðt z~ i ðt þ19tÞ; ~ þ 19tÞ is uncorrelated with z~ i ðt þ 19tÞ since (2) xðt

#0 )

þ

m X

C i ðt þ2ÞAðt þ 1Þz~ i ðt þ19tÞ

i¼1

"

 Cðt þ 2ÞAðt þ 1ÞJx ðt þ1Þ þ Cðt þ 2ÞBðt þ 1ÞJu ðt þ 1Þ

þ

C i ðtÞz~ i ðt9t1Þ

m X

0 ~ ~ g ¼ P ðt9t1ÞC 0 ðtÞ ¼ Efxðt9t1Þ½CðtÞ xðt9t1Þ

ðA:4Þ

where the second equation uses the fact that u(t) is orthogonal to x(t) and yðtÞ, t o t; the third equation uses ~ þ19tÞ and z~ i ðt þ 19tÞ are uncorrelated. In the fact that xðt the same way, other matrix gains Ju ðtÞ,fJ zi ðtÞgm i ¼ 1 can be computed. In the smoothing stage, the first-order smoother can be derived by

T 1 ðt þ 2Þ ¼ EfuðtÞx~ 0 ðt þ 29t þ 1Þg ¼ T 1 ðt þ1ÞA0 ðt þ 1ÞNðt9t þ 1ÞðAðt þ 1ÞJx ðt þ1Þ þBðt þ1ÞJ u ðt þ1ÞÞ0 T 2i ðt þ 2Þ ¼ EfuðtÞz~ 0i ðt þ 29t þ 1Þg ¼ T 2i ðt þ1ÞA0 ðt þ 1ÞNðt9t þ 1ÞðAðt þ 1ÞJzi ðt þ 1ÞÞ0 Nðt9t þ2Þ ¼ T 1 ðt þ 2ÞC 0 ðt þ 2Þ þ

^ ~ þ19tÞg uðt9t þ1Þ ¼ EfuðtÞ9Yðt þ 1Þg ¼ EfuðtÞ9YðtÞg þ EfuðtÞ9yðt

~ þ19tÞ þ ¼ E uðtÞ Cðt þ 1Þxðt (

m X

#0 )

#0 ) C i ðt þ1ÞAðtÞz~ i ðt9tÞ

0

0

m X

#0 )

~ C i ðt þ 1ÞAðtÞJ zi ðtÞÞyðt9t1Þ

i¼1 x

u

Appendix B. Proof of Theorem 2

x z m f 1 ðPx ðtÞ,fP zi ðtÞgm i ¼ 1 Þ Z f 1 ðQ ðtÞ,fQ i ðtÞgi ¼ 1 Þ

i¼1

  ¼ E uðtÞ Cðt þ1ÞBðtÞuðtÞðCðt þ 1ÞAðtÞJ x ðtÞ þ Cðt þ1ÞBðtÞJ u ðtÞ þ

ðA:7Þ

This proof mainly follows that of Theorem 1 in [16]. First, based on the monotonicity property of the Lyapunov operator, monotonicity properties for the iterative Riccati equation in (9) are obvious. If Px ðtÞ Z Q x ðtÞ,Pzi ðtÞ Z Q zi ðtÞ,i 2 f1; 2, . . . ,mg, then

~ ~ ¼ E uðtÞ Cðt þ 1ÞAðtÞxðt9tÞ þCðt þ 1ÞBðtÞuðt9tÞ m X

T 2i ðt þ 2ÞC i 0 ðt þ 2Þ

i¼1

C i ðt þ1Þz~ i ðt þ 19tÞ

i¼1

"

m X

By induction, the higher-order smoothers can be derived in the same way.

~ þ 19tÞ ^ ¼ uðt9tÞ þNðt9t þ1ÞS1 ðt þ 1Þyðt Nðt9t þ1Þ ¼ EfuðtÞy~ 0 ðt þ19tÞg " (

#0 )

# ~ þ19tÞ C i ðt þ2ÞAðt þ 1ÞJ zi ðt þ 1Þ yðt

i¼1

i¼1 x

þ

m X

0

If Px ð1Þ Z Px ð0Þ,Pzi ð1Þ Z Pzi ð0Þ,i 2 f1; 2, . . . ,mg, then x z m f 1 ðPx ðt þ1Þ,fP zi ðt þ 1Þgm i ¼ 1 Þ Z f 1 ðP ðtÞ,fP i ðtÞgi ¼ 1 Þ

0

¼ B ðtÞC ðt þ 1Þ½ðAðtÞJ ðtÞ þBðtÞJ ðtÞÞDðtÞ C ðt þ 1Þ m X ðAðtÞJ zi ðtÞDðtÞÞ0 C 0i ðt þ1Þ 

x z m f 2 ðPx ðtÞ,fP zi ðtÞgm i ¼ 1 Þ Z f 2 ðQ ðtÞ,fQ i ðtÞgi ¼ 1 Þ

ðA:5Þ

i¼1

Using the auxiliary sequences T 1 ,T 2i , they can be simplified as T 1 ðt þ 1Þ ¼ EfuðtÞx~ 0 ðt þ 19tÞg ¼ ½BðtÞðAðtÞJx ðtÞ þ BðtÞJ u ðtÞÞDðtÞ0

x z m f 2 ðPx ðt þ1Þ,fP zi ðt þ 1Þgm i ¼ 1 Þ Z f 2 ðP ðtÞ,fP i ðtÞgi ¼ 1 Þ

Second, according to the condition (14) and Lemma 4, there exists a set of linear operators ! m X i 0 0 0 ~ ~ ~ ~ ~ ~ ~ C i P C i K~ L 1 ðPÞ ¼ ðA þ K 1 CÞP 1 ðA þ K 1 CÞ þ K 1 3

1

i¼1

T 2i ðt þ 1Þ ¼ EfuðtÞz~ 0i ðt þ 19tÞg ¼ ½AðtÞJ zi ðtÞDðtÞ0 Nðt9t þ1Þ ¼ T 1 ðt þ 1ÞC 0 ðt þ 1Þ þ

m X

m X

i 0 C i P~ 3 C i 0 K~ 2

i¼1

T 2i ðt þ 1ÞC 0i ðt þ1Þ

i¼1

And, the second-order smoother can be derived by Nðt9t þ2Þ ¼ EfuðtÞy~ 0 ðt þ29t þ1Þg

þ K~ 2 C P~ 1 C 0 þ

!

ðA:6Þ

0 i ~ ¼ ðA þ K~ i C i ÞP~ i ðA þ K~ i C i Þ0 þ K~ i @C P~ 1 C 0 þ L~ 2 ðPÞ 3 3 3 3

1 X j i C j P~ 3 C 0j AðK~ 3 Þ0 jai

ðB:1Þ

C. Yu et al. / Signal Processing 92 (2012) 2497–2508

~ P~ i 4 L ~ where K~ 1 ¼ AP~ 1 C 0 S~ 1 , ~ 1 ðPÞ, ~ i ðPÞ, such that P~ 1 4 L 3 2 1 1 i i K~ 2 ¼ BDS~ , K~ 3 ¼ AP~ 3 C i 0 S~ . Then,

þK 1

m X

2507

z C i ðRzi ðtÞP i ðtÞÞC i 0

~k (1) For all W ¼ fW 1 ,fW i3 gm i ¼ 1 gZ 0, limk-1 L 1 ðWÞ ¼ 0, limk-1 ~ i Þk ðWÞ ¼ 0. ðL

(2) Let PðtÞ ¼ fP system

x

ðtÞ,fP zi ðtÞgm i ¼ 1g

and consider the linear 0

P x ðt þ 1Þ ¼ L~ 1 ðPðtÞÞ þ K~ 1 ðDD0 þ RÞK~ 1 0 þ ðB þ K~ 2 DÞðB þ K~ 2 DÞ0 þ K~ 2 RK~

m X

i

i

ðB:2Þ

initialized at Pð0Þ Z 0. Then, the sequence P(t) is bounded. Then, the following proof proceeds in several steps. First we show that the estimation error covariance sequence converges when the initial value is zero. Second we show that the estimation error covariance sequence converges to the same limiting point when the initial value is quite large. We finally show that the estimation error covariance sequence converges with arbitrary positive semi-definite initial values. First step: Let Q ðtÞ ¼ fQ x ðtÞ,fQ zi ðtÞgm i ¼ 1 g. When the iterative Riccati equation (9) is initialized at Q ð0Þ ¼ 0, a nondecreasing sequence can be generated according to the monotonicity of this iterative Riccati equation: 0 ¼ Q ð0Þ rQ ð1Þ r    r Q ðtÞ Then, we will show that this sequence is bounded. Let Pð0Þ Z 0. Then the sequence generated from the iterative equation (B.2) is bounded. Further, considering the optimality of the Kalman filter and Lemma 3, it can be shown that Q ðtÞ r PðtÞ. Thus, Q(t) is a nondecreasing sequence of matrices bounded above and converges to a fixed point, i.e.

ðB:4Þ where f1 ðÞ is defined in (10), KR is the Kalman gain associated with R(t), and K is the Kalman gain associated with P. In a similar way, we can get i

in (B.1) depend on K~ , and functions Functions i L 1 ðÞ,L 2 ðÞ in (B.4)–(B.5) depend on K . In addition, K is the steady optimal feedback gain of the steady state system described in (B.3). According to the optimality of the Kalman filter, the steady state Kalman filtering can achieve the minimum estimation error covariance:

f1 ðK ,PÞ r f1 ðK~ ,PÞ 0

f1 ðK ,PÞ ¼ L 1 ðPÞ þ K 1 ðDD0 þ RÞK 1 0

þ ðB þ K 2 DÞðB þK 2 DÞ0 þK 2 RK 2 0

f1 ðK~ ,PÞ ¼ L~ 1 ðPÞ þ K~ 1 ðDD0 þ RÞK~ 1 0

þ ðB þ K~ 2 DÞðB þ K~ 2 DÞ0 þ K~ 2 RK~ 2

ðB:6Þ

fi2 ðK ,PÞ r fi2 ðK~ ,PÞ i

i

i

i

i

i

fi2 ðK ,PÞ ¼ L 2 ðPÞ þ K 3 ðDD0 þ RÞðK 3 Þ0 þBB0 fi2 ðK~ ,PÞ ¼ L~ 2 ðPÞ þ K~ 3 ðDD0 þ RÞðK~ 3 Þ0 þBB0

ðB:7Þ

~k According to Lemma 4, for any P Z 0, it has limk-1 L 1 ~ i Þk ðPÞ ¼ 0, and f ðK~ ,PÞ, fi ðK~ ,PÞ are ðPÞ ¼ 0, limk-1 ðL 1

2

2

i

bounded. Thus, f1 ðK ,PÞ, f2 ðK ,PÞ are bounded as well. According to the optimality of the Kalman gain K , we can k

lim Q ðtÞ ¼ P 1

t-1

i

x

inequality in (B.4), we can arrive at limt-1 ðRx ðtÞP Þ ¼ 0 z limt-1 ðRzi ðtÞP i Þ ¼

i lim Q z ðtÞ ¼ P 3 t-1 i

This fixed point is the solution of the following equations: j

P 1 ¼ f 1 ðP 1 ,fP 3 gm j ¼ 1Þ j

i

ðB:5Þ

~ 1 ðÞ, L~ i ðÞ L 2

get limk-1 L 1 ðPÞ ¼ 0 and limk-1 ðL 2 Þk ðPÞ ¼ 0. Consider the

x

i

0

i¼1

z

~ ðPðtÞÞ þ K~ ðDD0 þRÞðK~ Þ0 þ BB0 P zi ðt þ 1Þ ¼ L 3 3 2

!

z

C i ðRzi ðtÞP i ÞC i 0 K 2 ¼ L1 ðRðtÞPÞ

0 r Rzi ðt þ 1ÞP i r L 2 ðRðtÞPÞ

2

i

0

K1

i¼1

þK 2 CP1 C 0 þ 2

!

P 3 ¼ f 2 ðP 1 ,fP 3 gm j ¼ 1Þ

ðB:3Þ

Second step: Next, we will show that the iterative Riccati j equation (9) initialized at Rx ð0Þ Z P 1 , Rzj ð0Þ Z P 3 also conx verges to the same limit. Denoting RðtÞ ¼ fR ðtÞ,fRzj ðtÞgm j ¼ 1g j and P ¼ fP 1 ,fP 3 gm j ¼ 1 g, then Rx ð1Þ ¼ f 1 ðRð0ÞÞ ZP 1 ¼ f 1 ðPÞ i

i

i

Rzi ð1Þ ¼ f 2 ðRð0ÞÞ ZP 3 ¼ f 2 ðPÞ A simple inductive argument establishes that Rx ðtÞ Z P 1 , i Rzi ðtÞ Z P 3 . Observe that x

0 r ðRx ðt þ 1ÞP Þ ¼ f 1 ðRðtÞÞf 1 ðPÞ ¼ f1 ðK R ,RðtÞÞf1 ðK ,PÞ x

r f1 ðK ,RðtÞÞf1 ðK ,PÞ ¼ ðA þ K 1 CÞðRx ðtÞP ÞðA þK 1 CÞ0

0. and Third step: We now establish that the Riccati iteration in (9) converges to P for all initial positive semi-definite conditions. Denote PðtÞ ¼ fP x ðtÞ,fPzj ðtÞgm j ¼ 1 g, Q ð0Þ ¼ 0 and Rð0Þ ¼ Pð0Þ þ P, and consider the monotonicity properties stated beforehand, we have Q ðtÞ r PðtÞ rRðtÞ,8t. Since both Q(t) and R(t) converge to P, we can obtain P ¼ lim Q ðtÞ r lim PðtÞ r lim RðtÞ ¼ P t-1

t-1

t-1

to complete the proof. References [1] R. Lane, Methods for maximum-likelihood deconvolution, Journal of the Optical Society of America A—Optics Image Science and Vision 13 (October) (1996) 1992–1998. [2] A. Doucet, P. Duvaut, Bayesian estimation of state-space models applied to deconvolution of Bernoulli–Gaussian processes, Signal Processing 57 (March) (1997) 147–161. [3] G. Pillonetto, B.M. Bell, Bayes and Empirical Bayes semi-blind deconvolution using eigenfunctions of a prior covariance, Automatica 43 (October) (2007) 1698–1712.

2508

C. Yu et al. / Signal Processing 92 (2012) 2497–2508

[4] J.M. Mendel, Lessons in Estimation Theory for Signal Processing, Communications, and Control, Prentice Hall, 1995. [5] J.M. Mendel, White-noise estimators for seismic data-processing in oil-exploration, IEEE Transactions on Automatic Control 22 (5) (1977) 694–706. [6] Z. Deng, Unifying and universal optimal white noise estimators for time-varying systems, Control Theory and Applications 21 (2003). [7] X. Sun, Y. Gao, Z. Deng, Information fusion white noise deconvolution estimators for time-varying systems, Signal Processing 88 (May) (2008) 1233–1247. [8] Y. Bar-Shalom, X.-R. Li, T. Kirubarajan, Estimation with Applications to Tracking and Navigation, Wiley, New York, 2001. [9] F. Wang, V. Balakrishnan, Robust Kalman filters for linear timevarying systems with stochastic parametric uncertainties, IEEE Transactions on Signal Processing 50 (April) (2002) 803–813. [10] Y. Luo, Y. Zhu, D. Luo, J. Zhou, E. Song, D. Wang, Globally optimal multisensor distributed random parameter matrices Kalman filtering fusion with applications, Sensors 8 (December) (2008) 8086–8103.

[11] W. Dekoning, Optimal estimation of linear discrete-time-systems with stochastic parameters, Automatica 20 (1) (1984) 113–115. [12] P. Sarri, G. Thomas, E. Sekko, P. Neveux, Myopic deconvolution combining Kalman filter and tracking control, in: Proceedings of the 1998 IEEE International Conference on Acoustics, Speech and Signal Processing, 1998, vol. 3, May 1998, pp. 1833–1836. [13] P. Jiang, J. Zhou, Y. Zhu, Globally optimal Kalman filtering with finite-time correlated noises, in: IEEE Conference on Decision and Control, December 2010, pp. 5007–5012. [14] G. Xu, H. Liu, L. Tong, T. Kailath, A least-squares approach to blind channel identification, IEEE Transactions on Signal Processing 43 (December) (1995) 2982–2993. [15] Y. Hua, Fast maximum likelihood for blind identification of multiple FIR channels, IEEE Transactions on Signal Processing 44 (March) (1996) 661–672. [16] B. Sinopoli, L. Schenato, M. Franceschetti, K. Poolla, M. Jordan, S. Sastry, Kalman filtering with intermittent observations, IEEE Transactions on Automatic Control 49 (September) (2004) 1453–1464.