Adaptive blind multiuser detection over flat fast fading channels using particle filtering

In this paper, we propose a method for blind multiuser detection (MUD) in synchronous systems over flat and fast Rayleigh fading channels. We adopt an autoregressive-moving-average (ARMA) process to model the temporal correlation of the channels. Based on the ARMA process, we propose a novel time-observation state space model (TOSSM) that describes the dynamics of the addressed multiuser system. The TOSSM allows an MUD with natural blending of low complexity particle filtering (PF) and mixture Kalman filtering (for channel estimation). We further propose to use a more efficient PF algorithm known as the stochastic M-algorithm (SMA), which, although having lower complexity than the generic PF implementation, maintains comparable performance.


INTRODUCTION
When multiuser detection (MUD) was introduced in the eighties, it has received a great deal of attention due to its ability to reduce multiple access interference (MAI) and potential for increasing the capacity of CDMA systems. Since then, numerous detectors have been proposed in the literature for both synchronous and asynchronous transmission and some popular ones include the decorrelating detector, the minimum mean square error (MMSE) detector, the multistage detector, and the decision feedback detector [1].
In practice, distortion in signal strength due to timevarying fading channels must be attended while performing MUD. Even though noncoherent detection methods as proposed in [2] are often appealing owing to their simplicity since no inference on fading channels is needed, coherent detection has been proved to deliver better performance [3]. With coherent detection, estimation of channels can be obtained with or without pilot signals. Between them, significant amount of research has been devoted to schemes without using pilot signals, or blind MUD methods. Blind MUD methods are bandwidth more efficient and the approaches proposed, to name a few, include the recursive least square (RLS) [4,5], subspace-based [6], expectationmaximization [7], genetic algorithm [8] and Kalman filtering [9,10,11,12,13,14]. However, most of the approaches cited above assume slow or quasi-static fading channels.
In this paper, we focus on blind MUD for fast flat Rayleigh fading channels and in synchronous systems. In particular, we assume to know a priori the second-order statistics of the underlying channel, based on which a mathematical tractable approximation using autoregressivemoving-average (ARMA) model is adopted. The approximation enables a dynamic state-space modeling (DSSM) of the problem, which lends itself naturally to a Kalmanfiltering-related detection solution. The use of Kalman filtering for blind MUD on similar modeling has been seen in [11,12,14], where the decision-directed approach was used to estimate the channel variable necessary for the Kalman filtering. One inherent drawback with the decision-directed approach is the error propagation, which greatly limits the performance of such implementation.
Recently, the combined (mixture) Kalman filtering and sequential importance sampling (particle filtering) algorithms have been applied to blind detection of convolutional codes [15], space-time trellis codes [16], and blind MUD [17] over fading channels. The mixture Kalman filtering (MKF) approach is shown to greatly reduce the error propagation of the decision-directed implementations and thus exhibits considerable performance improvement. However, in the proposed use of the MKF to blind MUD in [17], particle filtering (PF) was mainly intended for channel tracking and the embedded MUD at a symbol interval was achieved by an optimum detector, which has exponential complexity with the number of users. Consequently, the proposed MKF algorithm becomes prohibitively complex even for systems with moderate number of users.
In this paper, unlike all existing Kalman filtering detectors, a completely different viewpoint to multiuser systems is taken and we propose a novel time-observation state-space model (TOSSM). Even though the TOSSM is equivalent to the common DSSM, it allows the PF-based multiuser detection to be naturally blended with the mixture Kalman filtering for channel estimation. The new mixture Kalman filtering algorithm samples one user at a time and therefore permits efficient implementation. We further propose to use a more efficient PF algorithm known as the stochastic Malgorithm (SMA), which has shown to attain additional complexity reduction over the generic PF implementation and yet maintain comparable performance.
The rest of the paper is organized as follows. In Section 2, the problem of blind MUD is formulated. In Section 3, a novel TOSSM is described and in Section 4, the optimum solution is discussed. Particle filtering and SMA solutions are proposed in Sections 6 and 7, respectively. The simulation results are presented in Section 8. Section 9 contains some concluding remarks.

PROBLEM FORMULATION
Consider a synchronous CDMA system with a processing gain C and K users. Let T denote the symbol duration and s k (t) the normalized deterministic signature waveform assigned to the kth user. Then, at the nth symbol interval, the received signal y(t) can be expressed as a summation of K antipodally modulated synchronous signature waveforms plus noise, that is, where b n,k ∈ {−1, +1} is the BPSK modulated bit transmitted by the kth user, a k,n the CSI (fading coefficient) of the kth user, and u(t) the received zero mean additive complex white Gaussian noise with variance σ 2 . The cross-correlation between the signature waveforms of the users is given by the cross-correlation matrix R, where element r k1,k2 represents the cross-correlation between the signature waveform of the k 1 th and the k 2 th user and is defined as The channel for each user is considered as Rayleigh flat fading channel and ARMA processes can be adopted to model its time correlation with satisfaction [11,15,18]. Given an ARMA(r 1 , r 2 ) process, the CSI of the kth user at the nth interval a k,n can be represented as a n,k + φ k,1 a n−1,k · · · φ k,r1 a n−r1,k = ρ k,0 v n,k + · · · + ρ k,r2 v n−r2,k , where v n,k is an i.i.d. random complex Gaussian process that drives the ARMA process, {φ k,1 , . . . , φ k,r1 } and {ρ k,1 , . . . , ρ k,r2 } are the AR and MA coefficients of the model. We assume that we know a priori the second-order statistics of the underlying fading channel, and therefore the coefficients of the ARMA model can be precomputed so that the power spectral density of the ARMA process matches that of the fading channel. For convenience, we assume that r 1 = r 2 = r; otherwise zeros can be padded to the coefficients to make the orders equal. An equivalent form of (1) consists of a set of sufficient statistics represented by the matched filter output, The set of matched filter outputs y n = [y n,1 , . . . , y K,n ] T , where (·) T stands for matrix transpose, can be represented in vector-matrix form as where A n = diag{a n,1 , . . . , a n,K } is the diagonal matrix of the channel state information, b n = [b n,1 , . . . , b n,K ] T is the user date vector, and u n is the complex Gaussian noise vector with independent real and imaginary components and with covariance matrix equal to σ 2 R. Our objective is to perform sequential symbol detection without knowing the CSI a n,k , that is, blind multiuser detection.

TIME-OBSERVATION STATE-SPACE SYSTEM MODELING
A succinct mathematical representation of a time-varying system is the dynamic state-space model (DSSM). The statespace representation of CDMA systems in flat fading channels can be found in the existing literatures [11] and it can be expressed as In (6), h k,n for all k and b n are the unknowns to be estimated. Note that the observation y n is not linear in h k,n for all k and b n , and therefore the Kalman filter cannot provide the optimum solution. In fact, the optimum solution can be obtained by a so-called splitting Kalman filter, where, at time n, 2 n Kalman filters are required. The complexity of the splitting Kalman filter is exponential with both time and users and thus computational prohibited. Instead, particle filtering can be used to obtain good approximations of the optimum solution with reduced complexity. PF algorithms on (6) incorporated with Kalman filtering were proposed in [17]. However, as mentioned in the introduction, due to the structure of (6), particles of b n must be sampled jointly, and the complexity becomes exponential with the number of users. The prohibitive complexity on large user systems implies that this PF algorithm is infeasible for practical applications. To circumvent this difficulty, in the following we introduce a time-observation state-space model (TOSSM) for the system: In developing the TOSSM, we start with the Cholesky factorization of the cross-correlation matrix R as where F is a uniquely defined K × K lower triangular matrix. Now, right multiplying (F T ) −1 with the matched filter output, we obtain or, equivalently,ȳ n = FB n a n +ū n , where B n = diag{b n,1 , . . . , b n,K } is the diagonal user data matrix, and a n = [a n,1 , . . . , a n,K ] is the K × 1 vector of CSI. Since the covariance matrix ofū n becomes E[ū nū T n ] = σ 2 F − T RF −1 = σ 2 I, where I is an identity matrix,ȳ n is called the whitened matched filter (WMF) output. Next, define a tall channel vector of K(r + 1) × 1 as h n = [h T 1,n · · · h T K,n ] T and the channel transition becomes where We can thus express a n by h n in a compact form by a n = Ph n , where . Now by replacing a n in (11) by (13), we havē y n = FB n Ph n +ū n .
If we denote the kth row of F by f T k , the kth WMF outputȳ n can be written asȳ whereū n,k is the kth element ofū n . Now, instead of considering the system evolving only along time, we imagine a system progressing alternately along the path of time and the WMF observationsȳ n,k . The concept is further illustrated in Figure 1. To describe this new system, we must collapse the time index n and the observation index k into one timeobservation index l, where l = (n − 1)K + k. This conversion is reversible or, in other words, we can also calculate k and n from l by k = mod(l, K) and n = (l − k)/K + 1, where mod(k, K) is the k modulo K operation. In the following description of the TOSSM indexed by l, all k and n are assumed to be obtained from the corresponding l. Now, we introduce a K × K auxiliary matrix B l = diag{b n,1 , . . . , b n,k , 0 . . . , 0}. The state-space representation for the new time-observation system indexed by l can be then constructed as  and we call (16) the TOSSM. Note that (16) and (6) describe the same system. There are, however, key differences between the two models. Unlike (6), the state transitions of h l in the TOSSM are time (or index) varying, that is, at different l, different transition is applied. Specifically, when k = 1 or, equivalently, n increases by 1, h l updates according to the ARMA channel model, and otherwise when k = 1 and n remains unchanged from l − 1, h l is assumed to be static. Additionally, in the TOSSM, the number of the unknown user bits changes with l and especially, only one new unknown signal b n,k is included each time when l is incremented by one. Therefore, if we assume perfect detection at l − 1, that is, b n,1 , . . . , and b n,k−1 are known exactly, then there is only one unknown user bit to be detected. Note that in the conventional DSSM (6), K unknown users bits need to be detected altogether as the system evolves to time n. This is the key of the model that leads to efficient particle filtering solutions. We, however, want to stress that the decision on b n,k (except k = K) is not finalized at l. Since the observations from y l+1 up to y l+r with r = K − k all contain information about b n,k , the final decision is reached only at l + r, or in general, when k = K.

OPTIMUM BAYESIAN BLIND DETECTION
In a Bayesian framework, the optimum decision on b N can be obtained by the marginalized posterior mode (MPM) criterion, which is expressed as where p(b N |ȳ 1:NK ) is the posterior distribution that is essential for computing (17) Considering the TOSSM (16), we found the joint distribution in (8) and the mean m l and variance c l can be calculated sequentially through the Kalman filter. This is equivalent to say that p(ȳ l |b n,1:k , b 1:n−1 ) can be calculated from a run of the Kalman filter. Now, revisiting (18), we see that, to calculate p(b N |ȳ 1:NK ), p(ȳ l |b n,1:k , b 1:n−1 ) must be evaluated for 2 NK combinations of b 1:N , or 2 NK Kalman filters are needed, each of which corresponding to one possible combination. As a result, totally 2 NK Kalman filters are required for the MPM solution. The expansion of the numbers of the Kalman filters with l presents a tree structure illustrated in Figure 2.
The MPM solution has thus a complexity exponentially increasing with both time n and the number of users K. This is apparently a formidable task not possible for real applications. We, therefore, must resort to suboptimum solutions with manageable complexity. One choice is particle filtering.

A DECISION-DIRECTED APPROACH TO BLIND MUD
A decision-directed approach to blind MUD was proposed in [11] based on DSSM (6). We describe in the following a corresponding decision-directed approach on the TOSSM (16). Figure 2: The tree structure of the optimum solution. Each path in the tree represents a run of the Kalman filter.
One distinct feature of the decision-directed approach on the TOSSM is that the decision on only one user's bit is made at each l. Specifically, let b n,k−1 and h l−1 represent the decisions on b n,k−1 and h l−1 at l − 1, then the decision-directed approach at l can be summarized in Algorithm 1. Clearly, the above decision-directed algorithm is equivalent to one run of the Kalman filter, and therefore it is a lot simpler than the optimum MPM solution. Nevertheless, the user bit is determined based on the prediction of the channel states and the decisions on previous users' bits, and thus it is not optimum. Compared with the algorithm based on DSSM (6), at time k with k from 1 to K, the above algorithm makes a decision on one user at a time and updates the channel state vector h l whenever a decision is reached. The updated h l will then influence the decision on b n,k+1 . Therefore, in both a good and a bad way, decisions at early stages (smaller k) would have more impact on decisions at later stages (larger k) than those made by the algorithm on DSSM. If detection error exists in early stages, they will be propagated into later stages. It is therefore beneficial to rank the users according to the estimated SNR. The performance of the decisiondirected algorithm is, however, ultimately limited by error propagation.

PARTICLE FILTERING DETECTOR FOR BLIND MUD
Particle filtering belongs to the family of Monte Carlo sampling which aims at using samples to approximate posterior distribution. However, particle filtering distinguishes itself by employing a sequential importance sampling scheme, and in particular, it is designed for nonlinear and non-Gaussian systems described through state-space modeling such as the problem concerned.
In the context of the proposed problem, when y N , or equivalentlyȳ N , is observed at time N, the objective of particle filtering is to draw, say, J weighted random samples where δ(·) is the Dirac delta function, and hence the MPM solution of b by a simple weighted summation is where µ ( j) KN−1 is the weight update factor. Examining (24) and (25), we find that given w ( j) KN−1 and p(b N,1:K−1 , b 1:N−1 | y 1:NK−1 ), the importance function (24) and the weights (25) are known exactly as long as p(ȳ NK |b ( j) 1:N , y 1:NK−1 ) can be derived. In fact, we have indicated in Section 4 that p(ȳ NK |b 1:N ,ȳ 1:NK−1 ) can be calculated through the Kalman filter as l (i) are calculated the same way as shown in the appendix but for a set of b 1:NK given in (26). We can therefore obtain samples and weights using a recursive algorithm. To put the idea in concrete procedure, we assume that at l − 1, we have obtained from a previous recursion the trajectories (samples) {b . Using the recent observationsȳ l , we update the trajectories and weights as in Algorithm 2. This process of recursively obtaining particles is called particle filtering. After each recursion, the mean η ( j) l and covariance vectors Ξ ( j) l are passed on to the next recursion. From (21), we also see that to calculate all the elements NK is required. Therefore the decision on all the elements can only be made after recursion l = KN and the particles for b N,k for k = 1, 2, . . . , K − 1 must be stored.
In the above derivation of particle filtering, the adopted importance function is known as optimum in the sense that minimizes the variance of the weights. The above particle filtering procedure suffers from particle impoverishment, that is, after several recursions, some weights of the samples become negligible and stop contributing to the overall evaluation. To prevent it, we insert a residue resampling step [15] after every fixed recursion. Particularly, during the resampling at recursion l, the particles for b ( j) n,1:k , the mean vectors, and covariance matrices must be treated as a set in the resampling process.

STOCHASTIC M-DETECTOR FOR BLIND MUD
Recently, a every efficient particle filtering algorithm called stochastic M-algorithm (SMA) was proposed in [19] for problems with discrete unknowns. SMA can provide similar performance as generic particle filtering but with much reduced complexity. SMA can be considered as a particle filtering algorithm with the discrete delta functions as importance functions. In addition, each trajectory produces two samples (−1 and 1) for the binary case rather than one sample as in the generic PF. A key feature with SMA is that no two trajectories are identical, which is however rarely true with the generic PF. As a result, the SMA can provide more sample diversities with less trajectories than the generic PF. Nonetheless, notice that the number of trajectories doubles after each sampling and therefore a selection step is required BER. In Figure 4, we provide the BER versus SNR performance for a higher Doppler frequency of Ω d = 0.05. Similar observations can be drawn as for the previous case even though the overall performance of the detectors is worse, which is reasonable considering that the channels are fading faster.
In Figure 5, we provide the BER versus SNR of the first user for the different algorithms on a scenario of Ω d = 0.03. In addition, the users have different power. The difference between the power of the first user and that of the last user is 10 dB and the other users' powers are equally spaced in between. The genie-aided detector is also included as a lower bound. In this case, the PFDs and SMDs with 32 trajectories are approximately of the same order of magnitude as that of the genie-aided detector at SNRs of the first user less than 30 dB. As in the case of equal power, the results obtained by the SMDs with 4 and 32 trajectories are very close, especially after 30 dB, and comparable to that of the PFD. Again, the SMD with 4 trajectories is obviously more favorable since it requires only about 1/35 of complexity of the PFD. In this experiment, the performance of the decision-directed detector is much worse compared to the performance of the PDF and SMDs. For example, the latter achieves about 11 dB gain over the former at 10 −2 BER. In Figure 6, we provide the BER versus SNR performance for a Doppler frequency of Ω d = 0.05. Since the channels considered are fading faster, the performance of the detectors is worse. However, in general, similar observations to the tested detectors can be drawn. It is important to outline that the performance of the decision-directed detector gets worse in this case, for example, the PFD and SMDs achieve about 20 dB gain over the decision-directed detectors at 10 −2 BER.

CONCLUSION
In this paper, we proposed to solve blind MUD over flat fast fading channels. We constructed a novel time-observation state-space model, based on which efficient particle filtering and stochastic M detectors were proposed. Particularly, the detectors based on the SMA demonstrated greater potential than those using generic PF. The former can provide comparable performance as the latter but with much smaller complexity. where C l = f T k B l P. The second distribution p(h l |b n,1:k−1 , b 1:n−1 ,ȳ 1:l−1 ) is the predictive density which can be obtained from the predictive step of the Kalman filter [21,22], that is, where m l = C l µ l and K l = Σ l C H l /c l with c l = C l Σ l C H l + σ 2 . Now the integration in (A.1) is readily derived as p ȳ l |b n,1:k , b 1:n−1 ,ȳ 1:l = N m l , c l .