Graph representation for joint channel estimation and ... - IEEE Xplore
Graph Representation for Joint Channel Estimation and Symbol Detection Sau-Hsuan Wu, Urbashi Mitra and C.-C. Jay Kuo Department of Electrical Engineering University of Southern California Los Angeles, California 90089-2564 USA E-mail: {sauhsuaw, ubli}@usc.edu, [email protected] Abstract— A unified structure is proposed for joint channel tracking and symbol detection in multipath fading channels. Based on the expectation maximization (EM) algorithm, a group of recursive stochastic filters is derived for blind channel estimation, using the maximum likelihood and minimum mean squared error criteria. In conjunction with the BCJR algorithm, it is shown that the recursive procedure for joint channel estimation and maximum a posteriori symbol detection can be represented as a message-passing scheme on a factor graph. The graphical model not only provides a generalized implementation architecture for existing trainingbased schemes, it also greatly reduces the complexity of existing algorithms for blind channel tracking.
I. Introduction Joint channel estimation and symbol detection for timevarying channels has been an interesting, yet challenging topic for wireless communications. The challenge lies in the difficulty with channel estimation in the absence of training data and also with data detection in the absence of channel information; yet this is precisely what joint data and channel estimation must accomplish. A common approach for this task is the per survival processing (PSP) technique proposed in [1] where the Kalman filter and the Viterbi algorithm is employed for joint minimum mean squared error (MMSE) channel estimation and maximum likelihood sequence detection (MLSD). A weak point of PSP is that the number of channel estimators increases with the number of survival paths which is closely related to the length of channel memory. The PSP approach is later shown in [2, 3] to be a special realization of the more general EM algorithm, where transmission data are considered as discrete and deterministic unknown system parameters, while channel vectors are designated as systems’ hidden random states. Despite the significant efforts expended to reduce the complexity of PSP, its implementation cost is still high. Nevertheless, the recursive implementation of the EM algorithm in [3] does provide a common platform for the estimation of time-varying channel parameters. Independent of the research in PSP, an earlier work done by Kaleh et al [4] had investigated joint channel estimation and system detection for static channels, using the EM algorithm [5]. Instead of posing data sequences as deterministic unknown parameters, they considered data sequence as systems’ hidden random states. Having the static channel parameter as a deterministic unknown, they came up with a stochastic maximum likelihood (ML) estimator that involved the BCJR algorithm [6] for the evaluation of the a posteriori probability (APP) of transmitted data. This makes the algorithm a joint ML channel estimator and maximum a posteriori (MAP) symbol detector. A significant contribution of the approach in [4] is that only one filter is required for evaluating the channel estimate even without the knowledge of transmitted data. Recognizing the convenience and the advantage of this apThis research has been funded in part by the Integrated Media Systems Center, a National Science Foundation Engineering Research Center, Cooperative Agreement No. EEC-9529152, and in part by the National Science Foundation under Grant NSF/ANI-0087761.
IEEE Communications Society Globecom 2004
proach, in this work, we re-investigate joint estimation and detection problem for a variety of channel characteristics, based on the recursive EM algorithm [3]. By endowing the unknown channel vector with different statistical characteristics, (e.g. random versus deterministic and stationary versus dynamic), we will show that different types of blind channel estimators can evolve from the recursive EM framework without knowing the data sequence in advance. More interestingly, different estimators can be represented as different types of messages passing on Forney-style factor graphs [7] as shown in Fig. 3. To motivate the research results conducted in this work, we will first briefly review the EM algorithm, and give the sequential implementation of the stochastic ML estimator for a static channel vector. Then, the sequential stochastic estimator will be generalized to handle blind channel tracking for time-varying channels. II. Preliminaries We consider joint channel estimation and symbol detection for a Gaussian linear system of the form y(m) = G(ξm )hm + n(m), (1) where y is the received signal vector of dimension N ×1, hm is the channel vector of dimension L × 1 and G is a matrix of N × L depicting the system structure. The noise vector n is zero-mean, complex, Gaussian distributed, denoted by n ∼ N (0, C), with the covariance matrix C of dimension N × N . The symbol vector, ξm = [b(m), . . . , b(m − Lp )], describes the transmitted symbols embedded in the system matrix at symbol interval m. Lp is the length of the channel memory. The number of resolvable transmission paths is denoted by L. The channel vector can be expressed as hm = [h0 (m), h1 (m), · · · , hL−1 (m)]T where each entry relates to the channel coefficient on a certain path at time m. This model is quite general. Many wireless systems, including direct-sequence code-division multiple access (DS-CDMA) and space-time multiple input multiple output (MIMO) systems, can be expressed in this form. Our objective is the joint inference of ξm and hm , ∀m, given the structure of the matrix G(ξm ) and the observations of the system up to time m, denoted by Ym = [y(m), · · · , y(0)]. This problem can be characterized by a set {Ym , Ψm , Θm }, where the system’s hidden states up to time m are defined as Ψm = [ψm , · · · , ψ0 ], and the set of time-varying system parameters up to time m is Θm = [θm , · · · , θ0 ]. The system matrix G(ξm ) is fully characterized by the system state, ξm {b(m), qm }, where qm [b(m − 1), . . . , b(m − Lp )] is referred to as the symbol state. At time instant m, the system takes in b(m) and transits from state ξm−1 to state ξm . For simplicity, y(i) and b(i) are also denoted by yi and bi , respectively. To reduce the complexity for joint estimation and detection, the EM algorithm is often employed to approach the goals iteratively. To this end, we first define a Kullback-Leibler (K-L) measure for Θm at iteration as fol −1 and the AP P (ξi ; Θ −1 ) lows, given the estimate Θ m|m m|m
2329
0-7803-8794-5/04/$20.00 © 2004 IEEE
−1 ), i = 0, . . . , m: p(ξi |Ym ; Θ m|m −1 ) = Qm (Θm |Θ m|m
m
Eξi {− log{π N det(C)}
i=0
−1 }, (2) −[yi − G(ξi )hi ] C−1 [yi − G(ξi )hi ]|Ym ; Θ m|m H
where Eξi {·|Ym ; Θm } stands for taking expectation with i|m respect to (w.r.t) ξi using the AP P (ξi |Ym ; Θm ), and Θ denotes the parameter estimate at time i using information collected up to time m. The EM algorithm is summarized in two steps as follows, −1 ); E-step: Compute Qm (Θm , Θ m|m −1 ). M-step: Θ = arg maxΘ Qm (Θm , Θ m
m|m
m|m
A. The Recursive EM Algorithm The above algorithm provides an iterative framework for performing optimization. The complexity of the algorithm depends on the specific problem to which the algorithm is applied, as well as the underlying statistical assumptions. Even for a static channel where hm = h, ∀m, it is, in general, difficult to obtain the analytical solutions for h and C, simultaneously, except for some special cases, e.g. C = σ 2 I. Despite the mathematical difficulty, we can approximate the joint optimization for h and C, iteratively using the extended EM framework [2]. −1 ˆ −1 , C ˆ −1 }. Let Qm (h, C −1 |θˆ−1 ) be Define θˆ|m {h |m |m |m |m −1 . the K-L measure with C in (2) being replaced by C Maximizing the K-L measure w.r.t h gives ˆ = arg max Qm (h, C −1 |θˆ−1 ). h |m |m |m h
(3)
|m
|m
Similarly, maximizing this K-L measure w.r.t C gives ˆ , C|θˆ−1 ). ˆ = arg max Qm (h (4) C |m |m |m C
ˆ−1 Based on Theorem 1 of [2], we have Qm (θˆ|m |θ|m ) ≥ −1 −1 −1 −1 ˆ , C ˆ ) ≥ Qm (θˆ |θˆ ), which guarantees | θ Qm (h |m |m |m |m |m that the log likelihood (LLK) of θˆ is non-decreasing, de|m
−1 noted by LLK(θˆ|m ) ≥ LLK(θˆ|m ). In the sequel, we will follow this approach for the estimation of hm and C. Our goal is to develop a sequential algorithm for static channels, h, as well as a recursive algorithm for dynamic channels, hm . These can be done with the recursive EM algorithm [3], based on the second-order approximation to the M − step, 2 −1 −1 ) ∂ Qi (Θi |Θ m|m Θ · i|m = Θi−1|m − Θi =Θ i−1|m ∂ 2 Θi −1 ) ∂Qi (Θi |Θ m|m , (5) Θi =Θ i−1|m ∂Θ∗i 2
2
(·|·) ˆ0 0 L where ∂ ∂Q2iΘ(·|·) ∂∂ΘQ∗ i∂Θ T , Θm|m = {θm|m , Θm−1|m−1 }, i i i and L is the total number of iterations at each time step. 0 L For time varying parameters, θˆm|m = f (Θ m−1|m−1 ) is a predicted value given by the dynamic evolution function of θm , [c.f. (21)]. Based on the recursive EM algorithm, we shall discuss the estimators for both static and dynamic channels.
IEEE Communications Society Globecom 2004
|m−1 )|Ym ; θˆ−1 } = EΨm {log p(Ym , Ψm ; h, C |m |m−1 |θˆ−1 ) + SLLm (h, C |m−1 |Ym ; θˆ−1 ),(6) = Qm−1 (h, C |m |m where the synthetic log likelihood function (SLL) of {h, C} at time i is defined as −1 SLLi (h, C|Ym ; θˆ|m ) − log{π N det(C)} −1 −Eξi {[yi − G(ξi )h]H C−1 [yi − G(ξi )h]|Ym ; θˆ|m } −1 H −1 H Si|m yi − h]H Ci|m [Ci|m Si|m yi − h] + Di , (7) = −[Ci|m where Di is a constant term, and Ci|m ≡ Ci|m (C) and Si|m ≡ Si|m (C), which are defined as follows
Ci|m (C)
|X |
−1 [GH (xn )C−1 G(xn )]P (ξi = xn |Ym ; θˆ|m ), (8)
n=1
|m
ˆ , we further define another K-L Upon the acquisition of h |m ˆ , C|θˆ−1 ) by substituting h in (2) by h ˆ . measure Qm (h |m
B. Sequential Stochastic ML Estimator for Static Channels For static channels, the system parameter θ = {h, C} and the hidden state Ψm = {qm , · · · , q0 }. The estimate of ˆ |m , C |m }. θ using observations Ym is denoted by θˆ|m = {h ˆ ˆ We want to sequentially approach θ|m by having θi|m expressed in terms of θˆi−1|m , i = 0, · · · , m, and θˆm|m θˆ|m . ˆ |m , according to the iterative proceWe first evaluate h dure (3) and (4). The K-L measure w.r.t h is given by |m−1 |θˆ−1 ) Qm (h, C |m
Si|m (C)
|X |
−1 C−1 G(xn )P (ξi = xn |Ym ; θˆ|m ). (9)
n=1
As shown in (5), the sequential update of h requires the first and the second derivatives of (6). Define P−1 i|m = |m−1 |θˆ−1 ) ∂ 2 Qi (h,C
|m − With the Matrix Inversion ˆ i−1|m . h ∂2h Lemma, one can show that Pi|m = Pi−1|m − Ki|m Pi−1|m , (10) −1 −1 , (11) Ki|m = Pi−1|m (C i|m + Pi−1|m ) |m−1 ). The sequential stochastic ML i|m Ci|m (C where C (SS-ML) channel estimator can be shown equal to ˆ i|m = h ˆ i−1|m + Ki|m (h i|m − h ˆ i−1|m ), h (12) where, for simplicity, the iteration index is suppressed, i|m C −1 S H and h i|m i|m ym and Si|m Si|m (C|m−1 ). i|m maximizes the SLL, It is obvious from (7) that h −1 ˆ i|m is referred to SLLi (h, C|m−1 |Ym ; θ ). Therefore, h
|m
as the synthetic ML estimate of h at time i, using Ym . |m , the approach of (5) is rather complicated. For C ˆ m|m , C|θˆ−1 ) w.r.t. C Thus, we directly maximize Qm (h |m at the end of each iteration. The result is approximated by M L, |m (1 − 1 )C |m−1 + 1 C (13) C m m m M L is obtained by maximizing where C m ˆ m|m , C|Ym ; θˆ−1 ), and can be expressed as SLLm (h |m ML ˆ m|m ][ym − G ˆ m|m ]H m|m h m|m h [ym − G C m
ˆ m|m h ˆH G m|m h H −G (14) m|m m|m + Rm|m , −1 ˆ where Gi|m Eξi {G(ξi )|Ym ; θ|m } and ˆ m|m h ˆ H GH (ξm )|Ym ; θˆ−1 . m|m Eξ G(ξ ) h R m m m|m |m (15)
2330
0-7803-8794-5/04/$20.00 © 2004 IEEE
III. Graphical Approach for Iterative Stochastic Estimators We have given so far the sequential estimator for a static channel, based on the frequentist approach where channel parameters are considered as deterministic unknowns. The iterative procedure can also be obtained by the Forneystyle factor graph (FFG) [7]. The graph is comprised of the BCJR algorithm for the evaluation of AP P (ξi ), as well as a SS-ML for the inference of θ. The iteration itself is essentially a set of message-passing schemes on the graph. In the sequel, we will first use the FFG approach to develop a sequential Bayesian channel estimator. This estimator will be generalized to a stochastic Kalman-like filter for tracking time-varying random channels without training symbols. A. Sequential Stochastic MMSE Estimator for Static Channels For the Bayesian approach, the channel parameter h is considered as a random vector. Instead of estimating h itself as in the frequentist approach, the goal is to estimate the mean of h conditioned on Ym , denoted by −1 ]. The corresponding conditional coµ|m E[h|Ym ; θˆ|m
variance matrix is denoted by Σ|m E{[h − µ|m ][h − −1 }. Thus, the parameter set is redefined as µ|m ]H |Ym ; θˆ|m θ = {µ|m , C}. In addition to AP P (ξi ), the APP of h is also necessary, therefore Ψm = {h, q0 , · · · , qm } {h, qm }. The K-L measure in this case is given by |m−1 |θˆ−1 ) = EΨ {log p(Ym , q , h)|Ym ; θˆ−1 } Qm (µ|m , C |m
=
m
m
m
|m
−λi (µ|m , µ i|m ) − [µ|m − µ−1 ]H Σ−1 −1 [µ|m − µ−1 ] + E,
i=0
where E is a constant, and i|m [µ|m − µ λi (µ|m , µ i|m ) [µ|m − µ i|m ]H C i|m ],
(16)
−1 S H C i|m i|m yi ,
with µ i|m and µ−1 and Σ−1 being the initial mean and the covariance matrix of h, respectively. Notice that λi (µ|m , µ i|m ) is a Gaussian quadratic form. Thus, µ i|m can be associated with a Gaussian PDF: i|m ). Similarly, we also have P ( µi|m ; µ|m ) N (µ|m , C P (µ−1 ; µ|m ) N (µ|m , Σ−1 ). The EM is equivalent to µ |m = arg max P (µ−1 ; µ|m ) µ|m
m
P ( µi|m ; µ|m ).
(17)
i=0
−1|m = Σ−1 , the optimization is Let µ −1|m = µ−1 and Σ equivalent to finding the sequential expressions of µ i|m ∼ N (µ|m , Σi|m ), i = 0, · · · , m, that satisfy P ( µi|m ; µ|m ) P ( µ−1|m ; µ|m )
i
P ( µj|m ; µ|m )
j=0
µi|m ; µ|m ), = P ( µi−1|m ; µ|m )P (
where K i|m
= µ i−1|m + K i|m ( µi|m − µ i−1|m ), = Σi−1|m − K i|m Σi−1|m , −1 i−1|m (C Σ i|m
IEEE Communications Society Globecom 2004
i−1|m )−1 . +Σ
µˆ i−1|m
(19) (20)
=
µˆ i|m
Fig. 1. The Forney-style Factor Graph representation for the sequential stochastic MMSE channel estimator of h at time i.
The sequential estimator of C is still equal to (14) in this ˆ m|m replaced by µ case, by having h m|m . The FFG representation for the sequential stochastic MMSE (SS-MMSE) channel estimator given in Fig. 1 has local variables, µ i−1|m , µ i|m and µ i|m , and the corresponding local functions, P ( µi−1|m ), P ( µi|m ) and P ( µi|m ), respectively. The PDF P ( µi−1|m ) can be considered as the input LLK function of µ|m at time i − 1, P ( µi|m ) as the input LLK function of µ|m due to the new observation yi , and P ( µi|m ) as the output LLK function at time i. The equality sign denotes that the input and output LLK functions are all for the same channel vector µ|m . B. Recursive stochastic MMSE estimator for Dynamic Channels We use an AR filter to model dynamic channels . Define hm [hTm , · · · , hTm−p+1 ]T and a random vector vm ∼ N (0L×L , IL×L ). The channel dynamics is modeled by hm = Fhm−1 + Bvm F1 F2 . . . Fp B 0 ... 0 I 0 hm−1 + .. .. .. ... ... vm , (21) . . . 0 ... I 0 0 where B and Fi , i = 1, · · · , p, are of dimension L × L, and p is the order of the channel model. To setup a K-L measure of the EM algorithm for this case, the hidden state is given by ψi = {qi , hi }, i = 0, · · · , m, and Ψm = {ψ0 , · · · , ψm }. Let −1 ], µ −1 ] and µi|m E[hi |Ym ; Θ E[hi |Ym ; Θ m|m m|m i|m Σi,j|m , i, j, = 0, · · · , m, be the cross-covariance matrix between hi and hj . For convenience of expression, we also define i [hTi , · · · , hT0 ]T and Ui|m [µTi|m , · · · , µT0|m ]T , with the conditional covariance matrix of i denoted by Ξi|m = [Σi−j,i−l|m ]j+1,l+1 , 0 ≤ j, l ≤ i. The parameter set becomes Θm|m {C, Um|m }, and its estimate is denoted |m , Um|m }. We will focus on Um|m only. The m|m {C by Θ |m is similar to the SS-MMSE estimator. derivation for C The K-L measure in this case is given by −1 |Θ −1 ) Qm (Um|m , C |m m|m −1 } = EΨm {log p(Ym , qm , m )|Ym ; Θ m|m
(18)
and µ |m = µ m|m . This problem can solved with the FFG as shown in Fig. 1 for Gaussian linear systems, (see [8, 9] i|m are given by for details). The expressions of µ i|m and Σ µ i|m Σi|m
µ~ i|m
=
m
{−λi (µi|m , µ i|m ) − ηi (µi|m , µ ¯i|m ) + Ei } + F
i=1
−[µ0|m − µ−1 ]H Σ−1 0|m ),(22) −1 [µ0|m − µ−1 ] − λ0 (µ0|m , µ where Ei and F are constant terms, and i|m [ λi ( µi|m , µi|m ) [ µi|m − µi|m ]H C µi|m − µi|m ] (23) ηi (µi|m , µ ¯i|m ) [µi|m − µ ¯i|m ]H (BBH )−1 [µi|m − µ ¯i|m ], (24)
2331
0-7803-8794-5/04/$20.00 © 2004 IEEE
u~ i|m
,i (
U
+i
FL×pL 0
U
i|m
i −1|m
IiL×iL
+
U
◊ i |m
where Vi−1|m = E{[Ui−1|m − Ui−1|m ][µi−1|m − µ i−1|m ]H }. ¯ Then, it is straightforward to show that U = U and
Ji
=
i|m
U
i|m
JiBL×L
i
i|m
vi
Fig. 2. The Forney-style Factor Graph representation for the RSMMSE estimator for dynamic channel parameters hi , i = 0, · · · , m. The system, Gi , changes with time and Ji [IL×L , 0L×iL ]T .
ˇ [F1 , · · · , Fp ]. ˇ and F µ ¯i|m Fµ i−1|m Notice that, except for the constant terms Ei and F , all other terms are in a Gaussian quadratic form. Thus, we consider the mean vector µi|m as a random hyper-parameter i,i|m ). In this sense, ηm (µi|m , µ ∼ N ( µi|m , Σ ¯i|m ) can be associated with the PDF, P (µi|m |¯ µi|m ) N (¯ µi|m , BBH ). Similarly, −λi ( µi|m , µi|m ) is also associated with a PDF of −1 ). Thus, maximizing the form P ( µi|m |µi|m ) = N (µi|m , C i the K-L measure (22) w.r.t Um|m is equivalent to Um|m = arg max P (µ0|m |µ−1 )P ( µ0|m |µ0|m ) Um|m
m
P ( µi|m |µi|m )P (µi|m |¯ µi|m ),
i|m
= Ξ ¯ i|m + Ji BBH JH , where Ji [IL×L , 0L×iL ]H , Ξ i i|m 0 ≤ i ≤ m. The recursive stochastic MMSE (RS-MMSE) channel estimator is given by Ui|m = U + Ki|m ( µi|m − JH U ). (30)
(25)
i=1
where P (µ0|m |µ−1 ) N (µ−1 , Σ−1 ). This boils down to finding the recursive expression of Um|m that satisfies m|m ) P (Um|m ; Um|m ) N (Um|m , Ξ
i|m Ξ
, − Ki|m JH Ξ = Ξ i i|m i|m −1 Ji (C Ξ i|m i|m
−1 JH i Ξi|m Ji )
i|m
(31) JH i Ui|m
+ and = where Ki|m = ˇ 0|m , and the initial F µi−1|m . The initial mean, U0|m = µ 0|m = Σ 0,0|m , can be obtained by (19) covariance matrix, Ξ −1|m Σ−1 . and (20), with U−1|m µ−1 and Ξ Notice that, for every time index, i, the estimates from time 0 up to time i get updated, and the matrix dimen i|m also increases by L for each time step. Howsion of Ξ −1 + JH Ξ Ji = ever, the dimension of the matrix, C i i|m i|m ˇ H +BBH , that requires inversion is still L×L. ˇΣ F F i−1,i−1|m
The ML estimator of C still has the form of (13) by ˆ m|m with µ m|m . replacing h
IV. Graphical Representations It is clear that implementing stochastic filters requires i|m which are functions the synthetic variables µ i|m and C −1 of AP P (ξi |Ym ; Θm|m ). However, the evaluation of the −1 ) involves the calculation of the LLK AP P (ξi |Ym ; Θ m|m
−1 H −1 −1 −1 µi|m ] (C|m ) [yi − G(ξi ) µi|m ], in function, [yi − G(ξi ) the BCJR algorithm for each possible state of ξi . This makes the algorithm iterative in nature, due to the lack m of exact information of ξi . The iteration is guaranteed to P (U0|m ; U0|m ) P ( µi|m |µi|m )P (µi|m |¯ µi|m ) converge to a local optimum [5]. Therefore, the entire joint i=1 estimation and detection scheme is comprised of two parts: = P (Um−1|m ; Um−1|m )P ( µm|m |µm|m )P (µm|m |¯ µm−1|m ). the BCJR algorithm for AP P (ξi |Ym ) and a stochastic filter for the inference of Θm . The type of the filter depends 0|m ) . This probabilistic on the criterion that is used for the K-L measure of the where U0|m µ0|m ∼ N (U0|m , Ξ model implies that the time evolution characteristic of Ui|m EM algorithm. The factor graph representation of BCJR is shown in the upper part of Fig. 3. This is a standard can be modeled by a linear dynamic system: graphical representation and can be found in many artiˇ (26) cles, e.g. [9]. The FFG representation of the RS-MMSE µi|m = Fµi−1|m + Bνi µ i|m = µi|m + ωi , (27) filter is shown in the lower part of Fig. 3. The content of the system block denoted by Gi in the graph is shown in −1 ). where νi ∼ N (0, IL×L ) and ωi ∼ N (0, C Fig. 2 for the RS-MMSE filter. The differences between i Now the optimization problem in (22) is equivalent to filters lie in the content of the system block, Gi , as well as seeking the recursive expression of Ui|m . This is a standard the message being passed among them. We will introduce Kalman filtering problem. We give a more general solution the message-passing scheme for the RS-MMSE filter in the to this problem using the FFG as shown in Fig. 2. The next paragraph. The message-passing schemes for other filters are similar. system, denoted by Gi , is characterized by the matrix Assume we already have the channel estimates up to time ˇ L Fi FpL×L 0(i−p)L×L . (28) m − 1, denoted by Θ m−1|m−1 , and L is the total number IiL×iL of iterations. We want to advance to time m based on Notice that the dimension of Fi grows with time, which the knowledge accumulated up to m − 1. First, we demakes the system Gi a time-varying linear system. ˆ0 ˆ0 0 0 L fine Θ m|m = {θm|m , Θm−1|m } {θm|m , Θm−1|m−1 }, where ¯ , Ξ T ¯ i|m ). UsµTi|m , Ui−1|m ]T ∼ N (U Define U¯i|m [¯ i|m 0 L θˆm|m = f (Θ m−1|m−1 ) is a predicted value obtained by the ¯ ing the update approach provided in [8], we have U i|m = dynamic evolution function of θm , [c.f. (21)], or simply, T ˇµ 0 L [(F )T , Ui−1|m ]T , and θˆm−1|m−1 if the dynamic function is not accurate θˆm|m i−1|m 0 , is fed into the BCJR H H ˇ ˇ ˇ enough. The initial estimate, Θ FVi−1|m FΣi−1,i−1|m F m|m ¯ i|m , (29) Ξ H 0 ˇ |Y for the evaluation of AP P (ξ Vi−1|m F Ξi−1|m i m ; Θm|m ), i = 0, · · · , m.
IEEE Communications Society Globecom 2004
2332
0-7803-8794-5/04/$20.00 © 2004 IEEE
BCJR q0
α1
αi
β0
β i-1
U
−1|m
)
,0
αi+1
βi
βm-1
P (u~i|m | ui|m )
P (u~0|m | u0|m ) P(
qi
αm
U
P(
0|m
U
) P(
i -1|m
)
U
i|m
)
,m
Uˆ
Uˆ
0|m
U
P(
m|m
−1
10
−2
10
) −3
10
BER of joint estimation and detection BER with perfect channel information
Uˆ
i|m
BER v.s. SNR, L=3, fdTs=5× 10−3
0
10
P (u~m|m | um|m ) P(
,i
qm
1 1 βm = M 1 BER (dB)
p(q00 ) 1 p(q ) α0 = 0 M D p(q0 )
m|m
−4
10
0
2
4
6
8
Stochastic Filter
10 E / N (dB) b
Fig. 3. The graphical representation of the recursive joint channel estimation and symbol detection.
12
14
16
18
20
0
Fig. 5. BER vs. Eb /N0 of joint MAP detection and channel estimation using the RS-MMSE estimator. The number of paths L=3. −3
MSE v.s. SNR, L=3, f T =5× 10 b
0
d s
1
10
Channel Gain (E / N = 20dB, f T =10× 10−3) d s
2
1 0
10
True Channel Estimate of RS−MMSE Estimate of RS−ML
0.5
0 0
50
100
150
200
250
300
350
400
450
MSE (dB)
Amplitude
1.5
500
Channel Phase
100
Degree
−1
True Channel Estimate of RS−MMSE Estimate of RS−ML
0
10
−100
−200
−300 −2
10 −400
0
50
100
150
200 250 300 The Number of Symbols
350
400
450
500
2
4
6
8
10 E / N (dB) b
Fig. 4. The simulation results for the channel tracking over a flat fading channel with the normalized Doppler spread fd Ts = 5 × 10−3 .
0 ), the Upon the reception of new AP P (ξi |Ym ; Θ m|m marginalization for the synthetic probability, P ( µi|m |µi|m ) is executed, which gives µ i|m and Ci|m , i = 0, · · · , m, that required for the update of channel estimates shown in (30). i|m ready, as well as the initial setups for With µ i|m , C −1|m , the marginalization for P (Ui|m ; Ui|m ) U−1|m and Ξ can proceed, going from time 0 to time m sequentially, 1 which gives the channel estimate Um|m as well as the co1 for iteration 1. The iterations are repeated variance Ξ m|m
until convergence to a local optimum occurs. V. Computer Simulations The performance of the proposed recursive estimation schemes is demonstrated using the sliding-window RSMMSE and the sliding-window RS-ML estimators with window size equal to one. The total number of iterations is two. The fading channel is simulated using an ARMA model with normalized Doppler spread of fd Ts = 5 × 10−3 . The order of the AR channel model F used in (21) is 3. Thus, there is a mismatch between the model used for channel generation and channel estimation. Fig. 4 demonstrates the tracking performance for a flat fading channel, L=1, using the RS-MMSE and the RSML estimators. It is shown that the RS-MMSE estimator is more robust than the RS-ML estimator due to the noise term introduced in (21) to compensate for modelling errors. Fig. 5 presents the BER of the joint MAP detection and the RS-MMSE channel estimation for a multipath fading channel with L=3. Fig. 6 shows the MSE of channel estimation, normalized with L. IEEE Communications Society Globecom 2004
0
12
14
16
18
20
0
Fig. 6. MSE vs. Eb /N0 of joint MAP detection and channel estimation using the RS-MMSE estimator. The number of paths L=3.
VI. Conclusions A general recursive approach was proposed for joint channel tracking and MAP symbol detection, using the EM algorithm. Employing the concept of message passing, a FFG model was presented to unify the recursive procedure for the joint MAP detection and various channel estimation schemes derived herein. The graphical model not only provides an integrated and simpler structure for blind estimation, it turns out to be a generalized framework for existing training-based estimation schemes as well. References [1] R. Raheli, A. Polydoros, and C. K. Tzou, “Per-survivor processing: a general approach to MLSE in uncertain environments,” IEEE Trans. on Communications, vol. 43, no. 2/3/4, pp. 354– 364, Feb/Mar/Apr 1995. [2] H. Zamiri-Jafarian and S. Pasupathy, “Adaptive MLSDE using the EM algorithm,” IEEE Trans. on Communications, vol. 47, no. 8, pp. 1181–1193, Aug. 1999. [3] H. Zamiri-Jafarian and S. Pasupathy, “EM-based recursive estimation of channel parameters,” IEEE Trans. on Communications, vol. 47, no. 9, pp. 1297–1302, Sept. 1999. [4] G. K. Kaleh and R. Vallet, “Joint parameter estimation and symbol detection for linear or nonlinear unknown channels ,” IEEE Trans. on Communications, vol. 42, no. 7, pp. 2406–2413, July 1994. [5] A. P. Dempster, N. M. Laird, and D. B. Rubin, “Maximumlikelihood from incomplete data via the EM algorithm,” J. Roy. Statist. Soc., vol. 39, pp. 1–17, 1977. [6] L. R. Bahl, J. Cocke, F. Jelinek, and J. Raviv, “Optimal decoding of linear codes for minimizing symbol error rate,” IEEE Trans. on Information Theory, vol. 20, no. 2, pp. 284–287, Mar. 1974. [7] Jr. G. D. Forney, “Codes on graph: normal realizations,” IEEE Trans. on Information Theory, vol. 47, no. 2, pp. 520–548, Feb. 2001. [8] H.-A. Loeliger, “Least squares and Kalman filtering on Forney graphs,” in Codes, Graphs, and Systems, (festschrift in honour of David Forney on the occasion of his 60th birthday), R.E. Blahut and R. Koetter, eds. Kluwer, 2002, pp. 113–135. [9] F. R. Kschischang, B. J. Frey, and H.-A. Loeliger, “Factor graphs and the sum-product algorithm,” IEEE Trans. on Information Theory, vol. 47, no. 2, pp. 498–519, Feb. 2001.
2333
0-7803-8794-5/04/$20.00 © 2004 IEEE
I. Introduction Joint channel estimation and symbol detection for timevarying channels has been an interesting, yet challenging topic for wireless communications. The challenge lies in the difficulty with channel estimation in the absence of training data and also with data detection in the absence of channel information; yet this is precisely what joint data and channel estimation must accomplish. A common approach for this task is the per survival processing (PSP) technique proposed in [1] where the Kalman filter and the Viterbi algorithm is employed for joint minimum mean squared error (MMSE) channel estimation and maximum likelihood sequence detection (MLSD). A weak point of PSP is that the number of channel estimators increases with the number of survival paths which is closely related to the length of channel memory. The PSP approach is later shown in [2, 3] to be a special realization of the more general EM algorithm, where transmission data are considered as discrete and deterministic unknown system parameters, while channel vectors are designated as systems’ hidden random states. Despite the significant efforts expended to reduce the complexity of PSP, its implementation cost is still high. Nevertheless, the recursive implementation of the EM algorithm in [3] does provide a common platform for the estimation of time-varying channel parameters. Independent of the research in PSP, an earlier work done by Kaleh et al [4] had investigated joint channel estimation and system detection for static channels, using the EM algorithm [5]. Instead of posing data sequences as deterministic unknown parameters, they considered data sequence as systems’ hidden random states. Having the static channel parameter as a deterministic unknown, they came up with a stochastic maximum likelihood (ML) estimator that involved the BCJR algorithm [6] for the evaluation of the a posteriori probability (APP) of transmitted data. This makes the algorithm a joint ML channel estimator and maximum a posteriori (MAP) symbol detector. A significant contribution of the approach in [4] is that only one filter is required for evaluating the channel estimate even without the knowledge of transmitted data. Recognizing the convenience and the advantage of this apThis research has been funded in part by the Integrated Media Systems Center, a National Science Foundation Engineering Research Center, Cooperative Agreement No. EEC-9529152, and in part by the National Science Foundation under Grant NSF/ANI-0087761.
IEEE Communications Society Globecom 2004
proach, in this work, we re-investigate joint estimation and detection problem for a variety of channel characteristics, based on the recursive EM algorithm [3]. By endowing the unknown channel vector with different statistical characteristics, (e.g. random versus deterministic and stationary versus dynamic), we will show that different types of blind channel estimators can evolve from the recursive EM framework without knowing the data sequence in advance. More interestingly, different estimators can be represented as different types of messages passing on Forney-style factor graphs [7] as shown in Fig. 3. To motivate the research results conducted in this work, we will first briefly review the EM algorithm, and give the sequential implementation of the stochastic ML estimator for a static channel vector. Then, the sequential stochastic estimator will be generalized to handle blind channel tracking for time-varying channels. II. Preliminaries We consider joint channel estimation and symbol detection for a Gaussian linear system of the form y(m) = G(ξm )hm + n(m), (1) where y is the received signal vector of dimension N ×1, hm is the channel vector of dimension L × 1 and G is a matrix of N × L depicting the system structure. The noise vector n is zero-mean, complex, Gaussian distributed, denoted by n ∼ N (0, C), with the covariance matrix C of dimension N × N . The symbol vector, ξm = [b(m), . . . , b(m − Lp )], describes the transmitted symbols embedded in the system matrix at symbol interval m. Lp is the length of the channel memory. The number of resolvable transmission paths is denoted by L. The channel vector can be expressed as hm = [h0 (m), h1 (m), · · · , hL−1 (m)]T where each entry relates to the channel coefficient on a certain path at time m. This model is quite general. Many wireless systems, including direct-sequence code-division multiple access (DS-CDMA) and space-time multiple input multiple output (MIMO) systems, can be expressed in this form. Our objective is the joint inference of ξm and hm , ∀m, given the structure of the matrix G(ξm ) and the observations of the system up to time m, denoted by Ym = [y(m), · · · , y(0)]. This problem can be characterized by a set {Ym , Ψm , Θm }, where the system’s hidden states up to time m are defined as Ψm = [ψm , · · · , ψ0 ], and the set of time-varying system parameters up to time m is Θm = [θm , · · · , θ0 ]. The system matrix G(ξm ) is fully characterized by the system state, ξm {b(m), qm }, where qm [b(m − 1), . . . , b(m − Lp )] is referred to as the symbol state. At time instant m, the system takes in b(m) and transits from state ξm−1 to state ξm . For simplicity, y(i) and b(i) are also denoted by yi and bi , respectively. To reduce the complexity for joint estimation and detection, the EM algorithm is often employed to approach the goals iteratively. To this end, we first define a Kullback-Leibler (K-L) measure for Θm at iteration as fol −1 and the AP P (ξi ; Θ −1 ) lows, given the estimate Θ m|m m|m
2329
0-7803-8794-5/04/$20.00 © 2004 IEEE
−1 ), i = 0, . . . , m: p(ξi |Ym ; Θ m|m −1 ) = Qm (Θm |Θ m|m
m
Eξi {− log{π N det(C)}
i=0
−1 }, (2) −[yi − G(ξi )hi ] C−1 [yi − G(ξi )hi ]|Ym ; Θ m|m H
where Eξi {·|Ym ; Θm } stands for taking expectation with i|m respect to (w.r.t) ξi using the AP P (ξi |Ym ; Θm ), and Θ denotes the parameter estimate at time i using information collected up to time m. The EM algorithm is summarized in two steps as follows, −1 ); E-step: Compute Qm (Θm , Θ m|m −1 ). M-step: Θ = arg maxΘ Qm (Θm , Θ m
m|m
m|m
A. The Recursive EM Algorithm The above algorithm provides an iterative framework for performing optimization. The complexity of the algorithm depends on the specific problem to which the algorithm is applied, as well as the underlying statistical assumptions. Even for a static channel where hm = h, ∀m, it is, in general, difficult to obtain the analytical solutions for h and C, simultaneously, except for some special cases, e.g. C = σ 2 I. Despite the mathematical difficulty, we can approximate the joint optimization for h and C, iteratively using the extended EM framework [2]. −1 ˆ −1 , C ˆ −1 }. Let Qm (h, C −1 |θˆ−1 ) be Define θˆ|m {h |m |m |m |m −1 . the K-L measure with C in (2) being replaced by C Maximizing the K-L measure w.r.t h gives ˆ = arg max Qm (h, C −1 |θˆ−1 ). h |m |m |m h
(3)
|m
|m
Similarly, maximizing this K-L measure w.r.t C gives ˆ , C|θˆ−1 ). ˆ = arg max Qm (h (4) C |m |m |m C
ˆ−1 Based on Theorem 1 of [2], we have Qm (θˆ|m |θ|m ) ≥ −1 −1 −1 −1 ˆ , C ˆ ) ≥ Qm (θˆ |θˆ ), which guarantees | θ Qm (h |m |m |m |m |m that the log likelihood (LLK) of θˆ is non-decreasing, de|m
−1 noted by LLK(θˆ|m ) ≥ LLK(θˆ|m ). In the sequel, we will follow this approach for the estimation of hm and C. Our goal is to develop a sequential algorithm for static channels, h, as well as a recursive algorithm for dynamic channels, hm . These can be done with the recursive EM algorithm [3], based on the second-order approximation to the M − step, 2 −1 −1 ) ∂ Qi (Θi |Θ m|m Θ · i|m = Θi−1|m − Θi =Θ i−1|m ∂ 2 Θi −1 ) ∂Qi (Θi |Θ m|m , (5) Θi =Θ i−1|m ∂Θ∗i 2
2
(·|·) ˆ0 0 L where ∂ ∂Q2iΘ(·|·) ∂∂ΘQ∗ i∂Θ T , Θm|m = {θm|m , Θm−1|m−1 }, i i i and L is the total number of iterations at each time step. 0 L For time varying parameters, θˆm|m = f (Θ m−1|m−1 ) is a predicted value given by the dynamic evolution function of θm , [c.f. (21)]. Based on the recursive EM algorithm, we shall discuss the estimators for both static and dynamic channels.
IEEE Communications Society Globecom 2004
|m−1 )|Ym ; θˆ−1 } = EΨm {log p(Ym , Ψm ; h, C |m |m−1 |θˆ−1 ) + SLLm (h, C |m−1 |Ym ; θˆ−1 ),(6) = Qm−1 (h, C |m |m where the synthetic log likelihood function (SLL) of {h, C} at time i is defined as −1 SLLi (h, C|Ym ; θˆ|m ) − log{π N det(C)} −1 −Eξi {[yi − G(ξi )h]H C−1 [yi − G(ξi )h]|Ym ; θˆ|m } −1 H −1 H Si|m yi − h]H Ci|m [Ci|m Si|m yi − h] + Di , (7) = −[Ci|m where Di is a constant term, and Ci|m ≡ Ci|m (C) and Si|m ≡ Si|m (C), which are defined as follows
Ci|m (C)
|X |
−1 [GH (xn )C−1 G(xn )]P (ξi = xn |Ym ; θˆ|m ), (8)
n=1
|m
ˆ , we further define another K-L Upon the acquisition of h |m ˆ , C|θˆ−1 ) by substituting h in (2) by h ˆ . measure Qm (h |m
B. Sequential Stochastic ML Estimator for Static Channels For static channels, the system parameter θ = {h, C} and the hidden state Ψm = {qm , · · · , q0 }. The estimate of ˆ |m , C |m }. θ using observations Ym is denoted by θˆ|m = {h ˆ ˆ We want to sequentially approach θ|m by having θi|m expressed in terms of θˆi−1|m , i = 0, · · · , m, and θˆm|m θˆ|m . ˆ |m , according to the iterative proceWe first evaluate h dure (3) and (4). The K-L measure w.r.t h is given by |m−1 |θˆ−1 ) Qm (h, C |m
Si|m (C)
|X |
−1 C−1 G(xn )P (ξi = xn |Ym ; θˆ|m ). (9)
n=1
As shown in (5), the sequential update of h requires the first and the second derivatives of (6). Define P−1 i|m = |m−1 |θˆ−1 ) ∂ 2 Qi (h,C
|m − With the Matrix Inversion ˆ i−1|m . h ∂2h Lemma, one can show that Pi|m = Pi−1|m − Ki|m Pi−1|m , (10) −1 −1 , (11) Ki|m = Pi−1|m (C i|m + Pi−1|m ) |m−1 ). The sequential stochastic ML i|m Ci|m (C where C (SS-ML) channel estimator can be shown equal to ˆ i|m = h ˆ i−1|m + Ki|m (h i|m − h ˆ i−1|m ), h (12) where, for simplicity, the iteration index is suppressed, i|m C −1 S H and h i|m i|m ym and Si|m Si|m (C|m−1 ). i|m maximizes the SLL, It is obvious from (7) that h −1 ˆ i|m is referred to SLLi (h, C|m−1 |Ym ; θ ). Therefore, h
|m
as the synthetic ML estimate of h at time i, using Ym . |m , the approach of (5) is rather complicated. For C ˆ m|m , C|θˆ−1 ) w.r.t. C Thus, we directly maximize Qm (h |m at the end of each iteration. The result is approximated by M L, |m (1 − 1 )C |m−1 + 1 C (13) C m m m M L is obtained by maximizing where C m ˆ m|m , C|Ym ; θˆ−1 ), and can be expressed as SLLm (h |m ML ˆ m|m ][ym − G ˆ m|m ]H m|m h m|m h [ym − G C m
ˆ m|m h ˆH G m|m h H −G (14) m|m m|m + Rm|m , −1 ˆ where Gi|m Eξi {G(ξi )|Ym ; θ|m } and ˆ m|m h ˆ H GH (ξm )|Ym ; θˆ−1 . m|m Eξ G(ξ ) h R m m m|m |m (15)
2330
0-7803-8794-5/04/$20.00 © 2004 IEEE
III. Graphical Approach for Iterative Stochastic Estimators We have given so far the sequential estimator for a static channel, based on the frequentist approach where channel parameters are considered as deterministic unknowns. The iterative procedure can also be obtained by the Forneystyle factor graph (FFG) [7]. The graph is comprised of the BCJR algorithm for the evaluation of AP P (ξi ), as well as a SS-ML for the inference of θ. The iteration itself is essentially a set of message-passing schemes on the graph. In the sequel, we will first use the FFG approach to develop a sequential Bayesian channel estimator. This estimator will be generalized to a stochastic Kalman-like filter for tracking time-varying random channels without training symbols. A. Sequential Stochastic MMSE Estimator for Static Channels For the Bayesian approach, the channel parameter h is considered as a random vector. Instead of estimating h itself as in the frequentist approach, the goal is to estimate the mean of h conditioned on Ym , denoted by −1 ]. The corresponding conditional coµ|m E[h|Ym ; θˆ|m
variance matrix is denoted by Σ|m E{[h − µ|m ][h − −1 }. Thus, the parameter set is redefined as µ|m ]H |Ym ; θˆ|m θ = {µ|m , C}. In addition to AP P (ξi ), the APP of h is also necessary, therefore Ψm = {h, q0 , · · · , qm } {h, qm }. The K-L measure in this case is given by |m−1 |θˆ−1 ) = EΨ {log p(Ym , q , h)|Ym ; θˆ−1 } Qm (µ|m , C |m
=
m
m
m
|m
−λi (µ|m , µ i|m ) − [µ|m − µ−1 ]H Σ−1 −1 [µ|m − µ−1 ] + E,
i=0
where E is a constant, and i|m [µ|m − µ λi (µ|m , µ i|m ) [µ|m − µ i|m ]H C i|m ],
(16)
−1 S H C i|m i|m yi ,
with µ i|m and µ−1 and Σ−1 being the initial mean and the covariance matrix of h, respectively. Notice that λi (µ|m , µ i|m ) is a Gaussian quadratic form. Thus, µ i|m can be associated with a Gaussian PDF: i|m ). Similarly, we also have P ( µi|m ; µ|m ) N (µ|m , C P (µ−1 ; µ|m ) N (µ|m , Σ−1 ). The EM is equivalent to µ |m = arg max P (µ−1 ; µ|m ) µ|m
m
P ( µi|m ; µ|m ).
(17)
i=0
−1|m = Σ−1 , the optimization is Let µ −1|m = µ−1 and Σ equivalent to finding the sequential expressions of µ i|m ∼ N (µ|m , Σi|m ), i = 0, · · · , m, that satisfy P ( µi|m ; µ|m ) P ( µ−1|m ; µ|m )
i
P ( µj|m ; µ|m )
j=0
µi|m ; µ|m ), = P ( µi−1|m ; µ|m )P (
where K i|m
= µ i−1|m + K i|m ( µi|m − µ i−1|m ), = Σi−1|m − K i|m Σi−1|m , −1 i−1|m (C Σ i|m
IEEE Communications Society Globecom 2004
i−1|m )−1 . +Σ
µˆ i−1|m
(19) (20)
=
µˆ i|m
Fig. 1. The Forney-style Factor Graph representation for the sequential stochastic MMSE channel estimator of h at time i.
The sequential estimator of C is still equal to (14) in this ˆ m|m replaced by µ case, by having h m|m . The FFG representation for the sequential stochastic MMSE (SS-MMSE) channel estimator given in Fig. 1 has local variables, µ i−1|m , µ i|m and µ i|m , and the corresponding local functions, P ( µi−1|m ), P ( µi|m ) and P ( µi|m ), respectively. The PDF P ( µi−1|m ) can be considered as the input LLK function of µ|m at time i − 1, P ( µi|m ) as the input LLK function of µ|m due to the new observation yi , and P ( µi|m ) as the output LLK function at time i. The equality sign denotes that the input and output LLK functions are all for the same channel vector µ|m . B. Recursive stochastic MMSE estimator for Dynamic Channels We use an AR filter to model dynamic channels . Define hm [hTm , · · · , hTm−p+1 ]T and a random vector vm ∼ N (0L×L , IL×L ). The channel dynamics is modeled by hm = Fhm−1 + Bvm F1 F2 . . . Fp B 0 ... 0 I 0 hm−1 + .. .. .. ... ... vm , (21) . . . 0 ... I 0 0 where B and Fi , i = 1, · · · , p, are of dimension L × L, and p is the order of the channel model. To setup a K-L measure of the EM algorithm for this case, the hidden state is given by ψi = {qi , hi }, i = 0, · · · , m, and Ψm = {ψ0 , · · · , ψm }. Let −1 ], µ −1 ] and µi|m E[hi |Ym ; Θ E[hi |Ym ; Θ m|m m|m i|m Σi,j|m , i, j, = 0, · · · , m, be the cross-covariance matrix between hi and hj . For convenience of expression, we also define i [hTi , · · · , hT0 ]T and Ui|m [µTi|m , · · · , µT0|m ]T , with the conditional covariance matrix of i denoted by Ξi|m = [Σi−j,i−l|m ]j+1,l+1 , 0 ≤ j, l ≤ i. The parameter set becomes Θm|m {C, Um|m }, and its estimate is denoted |m , Um|m }. We will focus on Um|m only. The m|m {C by Θ |m is similar to the SS-MMSE estimator. derivation for C The K-L measure in this case is given by −1 |Θ −1 ) Qm (Um|m , C |m m|m −1 } = EΨm {log p(Ym , qm , m )|Ym ; Θ m|m
(18)
and µ |m = µ m|m . This problem can solved with the FFG as shown in Fig. 1 for Gaussian linear systems, (see [8, 9] i|m are given by for details). The expressions of µ i|m and Σ µ i|m Σi|m
µ~ i|m
=
m
{−λi (µi|m , µ i|m ) − ηi (µi|m , µ ¯i|m ) + Ei } + F
i=1
−[µ0|m − µ−1 ]H Σ−1 0|m ),(22) −1 [µ0|m − µ−1 ] − λ0 (µ0|m , µ where Ei and F are constant terms, and i|m [ λi ( µi|m , µi|m ) [ µi|m − µi|m ]H C µi|m − µi|m ] (23) ηi (µi|m , µ ¯i|m ) [µi|m − µ ¯i|m ]H (BBH )−1 [µi|m − µ ¯i|m ], (24)
2331
0-7803-8794-5/04/$20.00 © 2004 IEEE
u~ i|m
,i (
U
+i
FL×pL 0
U
i|m
i −1|m
IiL×iL
+
U
◊ i |m
where Vi−1|m = E{[Ui−1|m − Ui−1|m ][µi−1|m − µ i−1|m ]H }. ¯ Then, it is straightforward to show that U = U and
Ji
=
i|m
U
i|m
JiBL×L
i
i|m
vi
Fig. 2. The Forney-style Factor Graph representation for the RSMMSE estimator for dynamic channel parameters hi , i = 0, · · · , m. The system, Gi , changes with time and Ji [IL×L , 0L×iL ]T .
ˇ [F1 , · · · , Fp ]. ˇ and F µ ¯i|m Fµ i−1|m Notice that, except for the constant terms Ei and F , all other terms are in a Gaussian quadratic form. Thus, we consider the mean vector µi|m as a random hyper-parameter i,i|m ). In this sense, ηm (µi|m , µ ∼ N ( µi|m , Σ ¯i|m ) can be associated with the PDF, P (µi|m |¯ µi|m ) N (¯ µi|m , BBH ). Similarly, −λi ( µi|m , µi|m ) is also associated with a PDF of −1 ). Thus, maximizing the form P ( µi|m |µi|m ) = N (µi|m , C i the K-L measure (22) w.r.t Um|m is equivalent to Um|m = arg max P (µ0|m |µ−1 )P ( µ0|m |µ0|m ) Um|m
m
P ( µi|m |µi|m )P (µi|m |¯ µi|m ),
i|m
= Ξ ¯ i|m + Ji BBH JH , where Ji [IL×L , 0L×iL ]H , Ξ i i|m 0 ≤ i ≤ m. The recursive stochastic MMSE (RS-MMSE) channel estimator is given by Ui|m = U + Ki|m ( µi|m − JH U ). (30)
(25)
i=1
where P (µ0|m |µ−1 ) N (µ−1 , Σ−1 ). This boils down to finding the recursive expression of Um|m that satisfies m|m ) P (Um|m ; Um|m ) N (Um|m , Ξ
i|m Ξ
, − Ki|m JH Ξ = Ξ i i|m i|m −1 Ji (C Ξ i|m i|m
−1 JH i Ξi|m Ji )
i|m
(31) JH i Ui|m
+ and = where Ki|m = ˇ 0|m , and the initial F µi−1|m . The initial mean, U0|m = µ 0|m = Σ 0,0|m , can be obtained by (19) covariance matrix, Ξ −1|m Σ−1 . and (20), with U−1|m µ−1 and Ξ Notice that, for every time index, i, the estimates from time 0 up to time i get updated, and the matrix dimen i|m also increases by L for each time step. Howsion of Ξ −1 + JH Ξ Ji = ever, the dimension of the matrix, C i i|m i|m ˇ H +BBH , that requires inversion is still L×L. ˇΣ F F i−1,i−1|m
The ML estimator of C still has the form of (13) by ˆ m|m with µ m|m . replacing h
IV. Graphical Representations It is clear that implementing stochastic filters requires i|m which are functions the synthetic variables µ i|m and C −1 of AP P (ξi |Ym ; Θm|m ). However, the evaluation of the −1 ) involves the calculation of the LLK AP P (ξi |Ym ; Θ m|m
−1 H −1 −1 −1 µi|m ] (C|m ) [yi − G(ξi ) µi|m ], in function, [yi − G(ξi ) the BCJR algorithm for each possible state of ξi . This makes the algorithm iterative in nature, due to the lack m of exact information of ξi . The iteration is guaranteed to P (U0|m ; U0|m ) P ( µi|m |µi|m )P (µi|m |¯ µi|m ) converge to a local optimum [5]. Therefore, the entire joint i=1 estimation and detection scheme is comprised of two parts: = P (Um−1|m ; Um−1|m )P ( µm|m |µm|m )P (µm|m |¯ µm−1|m ). the BCJR algorithm for AP P (ξi |Ym ) and a stochastic filter for the inference of Θm . The type of the filter depends 0|m ) . This probabilistic on the criterion that is used for the K-L measure of the where U0|m µ0|m ∼ N (U0|m , Ξ model implies that the time evolution characteristic of Ui|m EM algorithm. The factor graph representation of BCJR is shown in the upper part of Fig. 3. This is a standard can be modeled by a linear dynamic system: graphical representation and can be found in many artiˇ (26) cles, e.g. [9]. The FFG representation of the RS-MMSE µi|m = Fµi−1|m + Bνi µ i|m = µi|m + ωi , (27) filter is shown in the lower part of Fig. 3. The content of the system block denoted by Gi in the graph is shown in −1 ). where νi ∼ N (0, IL×L ) and ωi ∼ N (0, C Fig. 2 for the RS-MMSE filter. The differences between i Now the optimization problem in (22) is equivalent to filters lie in the content of the system block, Gi , as well as seeking the recursive expression of Ui|m . This is a standard the message being passed among them. We will introduce Kalman filtering problem. We give a more general solution the message-passing scheme for the RS-MMSE filter in the to this problem using the FFG as shown in Fig. 2. The next paragraph. The message-passing schemes for other filters are similar. system, denoted by Gi , is characterized by the matrix Assume we already have the channel estimates up to time ˇ L Fi FpL×L 0(i−p)L×L . (28) m − 1, denoted by Θ m−1|m−1 , and L is the total number IiL×iL of iterations. We want to advance to time m based on Notice that the dimension of Fi grows with time, which the knowledge accumulated up to m − 1. First, we demakes the system Gi a time-varying linear system. ˆ0 ˆ0 0 0 L fine Θ m|m = {θm|m , Θm−1|m } {θm|m , Θm−1|m−1 }, where ¯ , Ξ T ¯ i|m ). UsµTi|m , Ui−1|m ]T ∼ N (U Define U¯i|m [¯ i|m 0 L θˆm|m = f (Θ m−1|m−1 ) is a predicted value obtained by the ¯ ing the update approach provided in [8], we have U i|m = dynamic evolution function of θm , [c.f. (21)], or simply, T ˇµ 0 L [(F )T , Ui−1|m ]T , and θˆm−1|m−1 if the dynamic function is not accurate θˆm|m i−1|m 0 , is fed into the BCJR H H ˇ ˇ ˇ enough. The initial estimate, Θ FVi−1|m FΣi−1,i−1|m F m|m ¯ i|m , (29) Ξ H 0 ˇ |Y for the evaluation of AP P (ξ Vi−1|m F Ξi−1|m i m ; Θm|m ), i = 0, · · · , m.
IEEE Communications Society Globecom 2004
2332
0-7803-8794-5/04/$20.00 © 2004 IEEE
BCJR q0
α1
αi
β0
β i-1
U
−1|m
)
,0
αi+1
βi
βm-1
P (u~i|m | ui|m )
P (u~0|m | u0|m ) P(
qi
αm
U
P(
0|m
U
) P(
i -1|m
)
U
i|m
)
,m
Uˆ
Uˆ
0|m
U
P(
m|m
−1
10
−2
10
) −3
10
BER of joint estimation and detection BER with perfect channel information
Uˆ
i|m
BER v.s. SNR, L=3, fdTs=5× 10−3
0
10
P (u~m|m | um|m ) P(
,i
qm
1 1 βm = M 1 BER (dB)
p(q00 ) 1 p(q ) α0 = 0 M D p(q0 )
m|m
−4
10
0
2
4
6
8
Stochastic Filter
10 E / N (dB) b
Fig. 3. The graphical representation of the recursive joint channel estimation and symbol detection.
12
14
16
18
20
0
Fig. 5. BER vs. Eb /N0 of joint MAP detection and channel estimation using the RS-MMSE estimator. The number of paths L=3. −3
MSE v.s. SNR, L=3, f T =5× 10 b
0
d s
1
10
Channel Gain (E / N = 20dB, f T =10× 10−3) d s
2
1 0
10
True Channel Estimate of RS−MMSE Estimate of RS−ML
0.5
0 0
50
100
150
200
250
300
350
400
450
MSE (dB)
Amplitude
1.5
500
Channel Phase
100
Degree
−1
True Channel Estimate of RS−MMSE Estimate of RS−ML
0
10
−100
−200
−300 −2
10 −400
0
50
100
150
200 250 300 The Number of Symbols
350
400
450
500
2
4
6
8
10 E / N (dB) b
Fig. 4. The simulation results for the channel tracking over a flat fading channel with the normalized Doppler spread fd Ts = 5 × 10−3 .
0 ), the Upon the reception of new AP P (ξi |Ym ; Θ m|m marginalization for the synthetic probability, P ( µi|m |µi|m ) is executed, which gives µ i|m and Ci|m , i = 0, · · · , m, that required for the update of channel estimates shown in (30). i|m ready, as well as the initial setups for With µ i|m , C −1|m , the marginalization for P (Ui|m ; Ui|m ) U−1|m and Ξ can proceed, going from time 0 to time m sequentially, 1 which gives the channel estimate Um|m as well as the co1 for iteration 1. The iterations are repeated variance Ξ m|m
until convergence to a local optimum occurs. V. Computer Simulations The performance of the proposed recursive estimation schemes is demonstrated using the sliding-window RSMMSE and the sliding-window RS-ML estimators with window size equal to one. The total number of iterations is two. The fading channel is simulated using an ARMA model with normalized Doppler spread of fd Ts = 5 × 10−3 . The order of the AR channel model F used in (21) is 3. Thus, there is a mismatch between the model used for channel generation and channel estimation. Fig. 4 demonstrates the tracking performance for a flat fading channel, L=1, using the RS-MMSE and the RSML estimators. It is shown that the RS-MMSE estimator is more robust than the RS-ML estimator due to the noise term introduced in (21) to compensate for modelling errors. Fig. 5 presents the BER of the joint MAP detection and the RS-MMSE channel estimation for a multipath fading channel with L=3. Fig. 6 shows the MSE of channel estimation, normalized with L. IEEE Communications Society Globecom 2004
0
12
14
16
18
20
0
Fig. 6. MSE vs. Eb /N0 of joint MAP detection and channel estimation using the RS-MMSE estimator. The number of paths L=3.
VI. Conclusions A general recursive approach was proposed for joint channel tracking and MAP symbol detection, using the EM algorithm. Employing the concept of message passing, a FFG model was presented to unify the recursive procedure for the joint MAP detection and various channel estimation schemes derived herein. The graphical model not only provides an integrated and simpler structure for blind estimation, it turns out to be a generalized framework for existing training-based estimation schemes as well. References [1] R. Raheli, A. Polydoros, and C. K. Tzou, “Per-survivor processing: a general approach to MLSE in uncertain environments,” IEEE Trans. on Communications, vol. 43, no. 2/3/4, pp. 354– 364, Feb/Mar/Apr 1995. [2] H. Zamiri-Jafarian and S. Pasupathy, “Adaptive MLSDE using the EM algorithm,” IEEE Trans. on Communications, vol. 47, no. 8, pp. 1181–1193, Aug. 1999. [3] H. Zamiri-Jafarian and S. Pasupathy, “EM-based recursive estimation of channel parameters,” IEEE Trans. on Communications, vol. 47, no. 9, pp. 1297–1302, Sept. 1999. [4] G. K. Kaleh and R. Vallet, “Joint parameter estimation and symbol detection for linear or nonlinear unknown channels ,” IEEE Trans. on Communications, vol. 42, no. 7, pp. 2406–2413, July 1994. [5] A. P. Dempster, N. M. Laird, and D. B. Rubin, “Maximumlikelihood from incomplete data via the EM algorithm,” J. Roy. Statist. Soc., vol. 39, pp. 1–17, 1977. [6] L. R. Bahl, J. Cocke, F. Jelinek, and J. Raviv, “Optimal decoding of linear codes for minimizing symbol error rate,” IEEE Trans. on Information Theory, vol. 20, no. 2, pp. 284–287, Mar. 1974. [7] Jr. G. D. Forney, “Codes on graph: normal realizations,” IEEE Trans. on Information Theory, vol. 47, no. 2, pp. 520–548, Feb. 2001. [8] H.-A. Loeliger, “Least squares and Kalman filtering on Forney graphs,” in Codes, Graphs, and Systems, (festschrift in honour of David Forney on the occasion of his 60th birthday), R.E. Blahut and R. Koetter, eds. Kluwer, 2002, pp. 113–135. [9] F. R. Kschischang, B. J. Frey, and H.-A. Loeliger, “Factor graphs and the sum-product algorithm,” IEEE Trans. on Information Theory, vol. 47, no. 2, pp. 498–519, Feb. 2001.
2333
0-7803-8794-5/04/$20.00 © 2004 IEEE
