Particle Filtering for Joint Symbol and Code Delay Estimation in DS Spread Spectrum Systems in Multipath Environment

We develop a new receiver for joint symbol, channel characteristics, and code delay estimation for DS spread spectrum systems under conditions of multipath fading. The approach is based on particle ﬁltering techniques and combines sequential importance sampling, a selection scheme, and a variance reduction technique. Several algorithms involving both deterministic and randomized schemes are considered and an extensive simulation study is carried out in order to demonstrate the performance of the proposed methods.


INTRODUCTION
Direct sequence (DS) spread spectrum systems are robust to many channel impairments, allow multiuser CDMA and low-detectability signal transmission, and, therefore, are widely used in different areas of digital communications. Unlike many other communication systems, however, spread spectrum receivers require additional code synchronization, which can be a rather challenging task under conditions of multipath fading, when severe amplitude and phase variations take place.
The problem of joint symbol, delay, and multipath estimation has been addressed in the literature before (see e.g., [1,2]), and proved to be a difficult one due to its inherited nonlinearity. The previously proposed approaches were mainly based on the use of the extended Kalman filter (EKF). However, many of them concentrated on the channel parameters and delay estimation only; moreover, in a number of cases, when EKF methods were applied, the estimated parameters were divergent [1]. Joint signal detection and channel estimation was performed using deterministic maximum likelihood (DML) methods [3,4]. However, since the unknown parameters of interest were assumed deterministic in this case, a serious drawback of DML-type approaches was the phenomenon of error propagation. Later, a stochastic maximum likelihood (ML) approach for the estimation of channel parameters was adopted with consequent symbol detection using Viterbi algorithms [5]. The space-alternating generalized expectation maximization (SAGE) scheme for maximum a posteriori (MAP) estimation was presented in [6].
In this paper, we propose to estimate the channel parameters, code delays, and symbols jointly using particle filtering techniques-a set of powerful and versatile simulation-based methods recently appeared in the literature (see [7] for a survey). The idea is to approximate the posterior distribution of interest by swarms of N (N 1) weighted points in the sample space, called particles, which evolve randomly in time in correlation with each other and either give birth to offspring particles or die according to their ability to represent the different zones of interest of the state space.
The methods have already been successfully applied to problems arising in digital communications, in particular, demodulation in fading channels [7,8,9] and detection in synchronous CDMA [10]. In all this work, the unknown fading channel characteristics were integrated out and only the symbols needed to be imputed. The algorithm, thus, made use of the structure of the model, and the unknown state involved discrete parameters only. Later investigation [10,11], however, revealed some concerns regarding the efficiency of the standard randomized particle filtering techniques in this context. It has been shown that, for a fixed computational complexity, more efficient deterministic schemes could be designed leading to an improved performance of the receiver. We attempt here to study these results further, and compare various randomized and nonrandomized approaches. Iltis [12] has recently developed a particle filtering method to address a problem closely related to ours. However, in his approach, the unknown symbol sequence is obtained through a standard algorithm, and only channel parameters and code delays are estimated using particle filtering. The problem we are dealing with is more complex, since it involves both discrete (symbols) and continuous-valued (delays) unknowns. The deterministic particle method, unfortunately, is not applicable directly in this case. However, in view of the recent results, we propose to combine it with sequential importance sampling for the mixed, discrete and continuous-valued parameter case, followed by an appropriate selection procedure. The resulting algorithm explores the state space in a more systematic way at little or no extra cost in comparison with the standard particle filtering, employing a suboptimal importance distribution. We develop and test this approach against other deterministic and stochastic schemes, and demonstrate its performance by means of an extensive simulation study.
The remainder of the paper is organized as follows. The model specifications and estimation objectives are stated in Section 2. In Section 3, a particle filtering method is developed for joint symbol/channel coefficients/code delay estimation. This section also introduces and reviews several alternative deterministic and stochastic schemes, with simulation results and comparisons presented in Section 4. Some conclusions are drawn at the end of the paper.

Transmitted waveform
We denote, for any generic sequence κ t , κ i: j (κ i , κ i+1 , . . . , κ j ) T , and let d n be the nth information symbol, and s tx (τ) the corresponding analog bandpass spread spectrum signal waveform transmitted in the symbol interval of duration T: where s n (·) performs the mapping from the digital sequence to waveforms and corresponds to the modulation technique employed, f 0 denotes the carrier frequency, and u(τ) is a wideband pseudonoise (PN) waveform defined by Here, c 1:K is a spreading code sequence consisting of K chips (with values {±1}) per symbol, η(τ − kT c ) is a rectangular pulse of unit height and duration T c , transmitted at (k − 1)T c < τ ≤ kT c , and T c is the chip interval satisfying the relation T c = T/K.

Channel model
The signal is passed through a noisy multipath fading channel which causes random amplitude and phase variations on the signal. The channel can be represented by a time-varying tapped-delayed line with taps spaced T s seconds apart, where T s is the Nyquist sampling rate for the transmitted waveform; T s = T c /2 due to the PN bandwidth being approximately 1/T c (see sampling theorem [13]). The equivalent discretetime impulse response of the channel is given by where t is a discrete time index, t = 1, 2, . . . . By L we understand the maximum number of paths (nonzero coefficients of h t ) of the channel [5], f (l) t are the complex-valued timevarying multipath coefficients arranged into the vector f t , and δ t,l denotes the Kronecker delta, which is 1 if t = l, and 0 otherwise.
We assume here that the channel coefficients f t and code delay θ t propagate according to the first-order autoregressive (AR) model: which corresponds to a Rayleigh uncorrelated scattering channel model; here A diag(α (0) , . . . , α (L−1) ), B diag(σ (0) f , . . . , σ (L−1) f ), with σ (l) f being the standard deviation, and α (l) accounting for the Doppler spread (see [2,14] for details and discussion on the use of higher-order AR). In this paper, matrices A, B, and parameters γ and σ θ are assumed known. Directions on the choice of these parameters are given in [2,14].

Received signal
The complex output of the channel sampled at the Nyquist rate, (in which case, t = 2K(n − 1) + 1, . . . , 2Kn samples correspond to the nth symbol transmitted, that is, d n ↔ y 2K(n−1)+1:2Kn ) can, thus, be expressed as where and σ 2 is the noise variance. 1 The noise sequences ϑ t , t , and v (l) t , l = 0, . . . , L − 1 are assumed mutually independent and independent of the initial states The received waveform s rx (τ) is obtained after ideal lowpass filtering of rectangular pulses and is given by [2] s rx (τ)

Estimation objectives
The symbols d n , which are assumed i.i.d., the channel characteristics f t , and the code delay θ t are unknown for n, t > 0. Our aim is to obtain sequentially in time an estimate of the joint posterior probability density of these parameters p(d 1:n , f 0:2Kn , θ 0:2Kn |y 1:2Kn ), and some of its characteristics, such as the MAP estimates of the symbols d 1:n = arg max d1:n p d 1:n |y 1:2Kn , ( 9 ) and the minimum mean square error (MMSE) estimates of the channel characteristics E(f 0:2Kn |y 1:2Kn ) and the delays E(θ 0:2Kn |y 1:2Kn ). This problem, unfortunately, does not admit any analytical solution and, thus, approximate methods must be employed. One of the methods that has proved to be useful in practice is particle filtering, and, in the next section, we propose a receiver based on the use of these techniques.

PARTICLE FILTERING RECEIVER
Particle filtering receivers have already been designed in [7,8,9], although for a much simpler case including symbols estimation only. The problem considered here is more complicated since an additional continuous parameter is involved, and, in this section, the particle filtering algorithm for the joint estimation of all unknown parameters is detailed. We begin our treatment with incorporating a variance reduction technique, namely, Rao-Blackwellisation, and then proceed with the derivation of the particle filtering equations for the estimation of the required posterior distribution. The alternative deterministic and stochastic approaches are considered at the end of the section.

Rao-Blackwellisation
In this paper, we follow a Bayesian approach and, given the measurements y 1:2Kn , base our inference on the joint posterior distribution p(d 1:n , f 0:2Kn , θ 0:2Kn |y 1:2Kn ). A straightforward application of particle filtering would, thus, focus on the estimation of this joint probability distribution, and, consequently, obtain the estimates of d 1:n , f 0:2Kn , and θ 0:2Kn sequentially in time. It is beneficial, however, to improve the standard approach by making most of the structure of the model and applying the variance reduction techniques.
Indeed, similar to [7,8,9,15], the problem of estimating p(d 1:n , f 0:2Kn , θ 0:2Kn |y 1:2Kn ) can be reduced to a one of sampling from a lower-dimensional posterior p(d 1:n , θ 0:2Kn |y 1:2Kn ). If the approximation of p(d 1:n , θ 0:2Kn |y 1:2Kn ) could be obtained, say, via particle filtering: one could compute the probability density p(f 0:2Kn |y 1:2Kn , d 1:n , θ 0:2Kn ) using the Kalman filter associated with (4) and (6). As a result, the posterior p(f 0:2Kn |y 1:2Kn ) could be approximated by a random mixture of Gaussians leading to lower variance of the estimates and, therefore, increased algorithm efficiency [15]. Strictly speaking, we are interested in estimating the information symbols only with the tracking of the channel being naturally incorporated into the proposed algorithm. However, following this approach, the MMSE (conditional mean) estimates of fading coefficients can, of course, be obtained if necessary as follows: with E[f 2K(n−1)+1:2Kn |y 1:2Kn , d (i) 1:n , θ (i) 0:2Kn ] being computed by the Kalman filter, with 2K steps required for each symbol transmitted.

Particle filtering algorithm
We can now proceed with the estimation of p(d 1:n , θ 0:2Kn |y 1:2Kn ) using particle filtering techniques. The method is based on the following remark. Suppose N particles can be easily simulated according to an arbitrary convenient importance distribution π(d 1:n , θ 0:n |y 1:n ) (such that p(d 1:n , θ 0:n |y 1:n ) > 0 implies π(d 1:n , θ 0:n |y 1:n ) > 0). Then, using the importance sampling identity, an estimate of p(d 1:n , θ 0:n |y 1:n ) is given by the following point mass approximation: where w (i) n are the so-called normalized importance weights, and y n denotes y n = y 2K(n−1)+1:2Kn (16) for n = 1, 2, . . .. The distribution π(d (i) 1:n , θ (i) 0:n |y 1:n ) has to admit π(d (i) 1:n−1 , θ (i) 0:n−1 |y 1:n−1 ) as a marginal distribution so that one could propagate this estimate sequentially in time without subsequently modifying the past simulated trajectories. The weights w (i) n could also be updated online in this case: The sequential importance sampling described above is combined with a selection procedure when the effective sample size N eff falls below some fraction of N, say N thres (see [15] for details). This helps to avoid the degeneracy of the algorithm by discarding particles with low normalized importance weights and multiplying those with high ones.

Selection step
If N eff < N thres , multiply/discard particles n to obtain N unweighted particles Algorithm 1: Particle filtering algorithm.

Implementation issues
The choice of importance distribution and selection scheme is discussed in [16]; depending on these choices, the computational complexity of the algorithm varies.
If K is long, it is useful to resample the particles at intermediate steps between t = 2K(n − 1) + 1 and t = 2Kn. One can also use Markov chain Monte Carlo (MCMC) steps to rejuvenate the particles and in particular d n .
The importance weights in this case can be calculated as where θ n is drawn from the prior Gaussian distribution with mean γθ n−1 and variance σ 2 θ : The importance weight w (i) n in (27) does not actually depend on d (i) n , and the weights evaluation and selection steps can be done prior to the sampling of d (i) n as in Algorithm 3. For each symbol detection, this procedure requires the evaluation of the M 2K-step-ahead Kalman filters, which is quite computationally expensive. Further research should, therefore, concentrate on development of other more efficient suboptimal importance distributions on a case by case basis.

Selection
As far as the selection step is concerned, a stratified sampling scheme [18] is employed in this paper since it has the minimum variance one can achieve in the class of unbiased schemes [19]. The algorithm is based on generating N points equally spaced in the interval [0, 1], with the number of offspring N i for each particle being equal to the number of points lying between the partial sums of weights q i−1 and q i , t . The procedure can be implemented in O(N) operations.

Deterministic particle filter
The use of the suboptimal importance distribution described in Section 3.3.1 increases the efficiency of the algorithm in comparison with the standard approach using the prior.
However, as shown in [11], if one already opts for (27), and all the calculations have to be performed anyway, it might be better to base our approximation of p(d 1:n , θ 0:n |y 1:n ) directly on p N×M d 1:n , θ 0:n |y 1:n where corresponding weights w (i,m) n are equal to and θ (i) n is still drawn from its prior (28). Indeed, all possible "extensions" of the existing state sequences at each step n are considered in this case, and one does not discard unnecessarily any information by selecting randomly one path out of the M available. In the above expression, w (i) n−1 is the weight of the "parent" particle, which has M "offspring" instead of the usual one, resulting in a total number of N × M particles at each stage. This number increases exponentially with time, and, therefore, a selection procedure has to be employed at each step n.
The simplest way to perform such selection is just to choose the N most likely offspring and discard the others (as, e.g., in [20]). The superiority of this approach over other methods in the fully discrete framework is shown in [10,11]. A more complicated procedure involves preserving the particles with high weights and resampling the ones with low weights, thus reducing their total number to N. An algorithm of this type is presented in [21] but other selection schemes can be designed. Contrary to the case involving the discrete parameters only, in this scenario, a resampling scheme with replacement could be employed, since θ (i) n is chosen randomly. Therefore, stratified resampling could be used in order to select N particles from N × M particles available.
Whether we choose to preserve the most likely particles, employ the selection scheme proposed in [21], or stratified resampling, the computational load of the resulting algorithms at each time step n is that of N × M × 2K Kalman filters, and the selection step in the first two cases is imple- M is large, which is the case in many applications, all these methods are too computationally expensive to be used, and one should employ a standard particle filter.

SIMULATION RESULTS
In the following experiments the bit error rate (BER) and the tracking delay error (θ t − θ t ) were evaluated by means of computer simulations. Gray-encoded M-ary differential phase shift keyed (MDPSK) signals were employed, with mapping function In order to assess the performance of the proposed approaches we first applied them to a simpler case of synchronization in flat fading conditions, L = 1, for a system with no spectrum spreading employed, c 1 = 1, K = 1. In the first experiment, 4DPSK signals were considered with the average signal-to-noise ratio (SNR) varying from 5 to 20 dB. The AR coefficients for the channel, (4), were set to α (0) = 0.999, σ (0) f = 0.01, and the delay model parameters in (5) were chosen to be the same, γ = 0.999 and σ θ = 0.01. The BER obtained by the particle filtering receiver employing prior (PFP) and suboptimal (PFS) importance distributions, and the deterministic receiver preserving N most likely particles (DML) and using stratified resampling (DSR) is presented in Figure 1. The marginal maximum a posteriori estimate (MMAP) d n = arg max dn p d n |y 1:2Kn (32) was employed to obtain the symbols. The number of particles used in these algorithms was equal to N = 100, and little or no improvement in BER was gained by increasing this number for deterministic schemes. For the randomized approaches, the number of particles required to achieve the BER of DSR algorithm was equal to N = 1200. In Figure 2, the mean square delay error (MSE) is presented as a function of the number of particles N for SNR = 10 dB: where L d is a length of the symbol sequence, L d = 1000. The results for the different SNRs are given in Figure 3. As one can see, the deterministic particle filter with stratified resampling slightly outperforms the receiver selecting most likely particles, and is more efficient than both standard particle filtering schemes.
The particle filtering approach was also compared with the EKF method [2] in the second experiment. The algorithm was simplified to consider channel and code delay estimation only (the transmitted symbols were assumed known as, e.g., with pilot symbols being used). Otherwise, simulation set-up was the same (N = 100 particles were employed). The results for the SNR = 10 dB presented in Figure 4 demonstrate good performance of the particle filtering method. Please note that deterministic particle filter is in  a sense a variant of the conventional per-survivor processing (PSP) algorithm combined with the randomized particle filtering procedure for the channel and delay estimation. Thus, this procedure has proved more efficient, even in the situation which is more favorable for EKF, that is, which does not involve the uncertainty associated with the unknown transmitted symbols. Simulations demonstrating the superiority of the particle filtering approach over RAKE receiver for related problems are presented in [22]; for symbol detection, the comparison of the particle filtering methods with other well-known approaches could be found in [8,9,23].  Finally, we applied the proposed algorithm to perform joint symbols/channel coefficients/code delay estimation for DS spread spectrum systems with K = 15, L = 4. A binary DPSK modulation scheme was employed with the multipath channel response and AR coefficients chosen as in Channel B in [2]. As shown in Figure 5, the algorithm employing 100 particles exhibits good BER performance. A tracking error trajectory for 100 information symbols (corresponding to 1500 chips and 3000 channel samples) and an average SNR equal to 10 dB is presented in Figure 6. Figure 7 also illustrates the mean square delay error as a function of SNR for L d = 1000.

CONCLUSION
In this paper, we propose the application of particle filtering techniques to a challenging problem of joint symbols, channel coefficients, and code delay estimation for DS spread spectrum systems in multipath fading. The algorithm is designed to make use of the structure of the model, and incorporates a variance reduction technique. The work is based on the recent results on the superiority of the DML approach in a fully discrete environment [10,11]. The method cannot be applied straightforwardly, however, and several procedures combining both deterministic and randomized schemes are considered. The algorithms are tested and compared. Although computer simulations show that all methods are capable of providing good performance, in this particular case involving additional continuous-valued parameters, the deterministic scheme employing stratified resampling turns out to be the most efficient one. The choice of the algorithm might, however, be application-dependent, so further investigation is necessary. The receiver can be extended to address multiuser DS-CDMA transmission, using the techniques proposed in [24], for example, or simplified to consider channel tracking only since it is naturally incorporated in the proposed algorithm. Future research should concentrate on the development of suboptimal importance distributions and selection schemes capable of increasing the algorithm efficiency.