DOA Estimation for Local Scattered CDMA Signals by Particle Swarm Optimization

This paper deals with the direction-of-arrival (DOA) estimation of local scattered code-division multiple access (CDMA) signals based on a particle swarm optimization (PSO) search. For conventional spectral searching estimators with local scattering, the searching complexity and estimating accuracy strictly depend on the number of search grids used during the search. In order to obtain high-resolution and accurate DOA estimation, a smaller grid size is needed. This is time consuming and it is unclear how to determine the required number of search grids. In this paper, a modified PSO is presented to reduce the required search grids for the conventional spectral searching estimator with the effects of local scattering. Finally, several computer simulations are provided for illustration and comparison.


Introduction
Adaptive array techniques have been developed for enhancing the performance of code-division multiple access (CDMA) systems. In the CDMA system, each user employs a unique pseudo-noise (PN) codeword to identify their data stream. Each user's transmission interferes with all the other users and causes multiple access interference (MAI). The use of antenna arrays as a tool for improving coverage, reducing interference, and increasing capacity in CDMA systems has attracted significant interest [1]. Multiple propagation is common in cellular systems, which may result from local scatters in the vicinity of the sources. In macrocellular environments, high base station antennas will receive these locally scattered signals from the mobile terminal, which are coherent and confined to a narrow OPEN ACCESS angular region [2]. In addition, the observation period is assumed to be short in comparison to the coherence time of the channel so that the channel may be modeled as time-invariant. A mobile terminal with local scattering can be well approximated by a single point source, as seen from the base station. The model can express the mobile and coherent interference from the surrounding environment under non line-of-sight conditions as another arithmetical model. In cellular radio systems, antenna array processing is considered to be a useful method to solve problems such as multipath and co-channel interference and increase capacity [1]. There have been some studies appraising the impact of a local scattering channel model on CDMA system, which propose uniform linear array (ULA) geometry [3]. They use a first-order Taylor expansion to linearize the nonlinear conventional array manifold produced by scatters and develop the generalized array manifold (GAM) model which can obtain better nominal DOA estimation for mobiles. Here, the multiple signal classification (MUSIC) searching function derived by the steering vector of the GAM model is termed as GMUSIC. However, in all these studies, estimation of direction of arrival (DOA) of the desired signals is required. Array outputs aligned with code-matched filter can make the DOA estimation of multiple sources equivalent to that for single source localization problem in a noisy environment. With the advantage of code-matched filter inherent in the CDMA system, it has been proved that the multiple signal classification estimator [4] can obtain an unbiased DOA estimation with low mean-square-error (MSE) [5]. It also contributes to solving the limitation that the number of array elements must be more than the number of impinging sources, but for a conventional spectral searching DOA estimator such as MUSIC, its searching complexity and estimating accuracy strictly depend on the number of search grids used during the search. This is time consuming and it is unclear how to determine the required number of search grids. Thus, in this paper, we will focus on the GAM model of ULA to employ the DOA estimation in multipath reflection scenario from scatters.
Particle swarm optimization (PSO) is a population-based stochastic optimization paradigm, in which each agent, named particle, of the population, named swarm, is thought of as a collision-proof bird and used to represent a potential solution [6]. Like a genetic algorithm (GA), PSO starts by initializing a population of random solutions and searches for optima by updating generations. But, PSO does not use any evolution operators. In PSO, particles fly through the problem space by following their own experience and the best experience attained by the swarm as a whole. In contrast to analytical or general heuristic methods, PSO is computationally efficient and has great capability of escaping local optima [7,8]. In addition, a key characteristic of PSO is that the algorithm itself is highly robust yet remarkably simple to implement, while processing similar capabilities as other evolutionary algorithms such as GA [9]. A maximum likelihood (ML) criterion with hard-constraint PSO based solution applied to DOA estimation is proposed in [10], that explores the potentially superior performance at less computational costs. The same PSO algorithm based on the ML methodology is also derived in [11], where the cost function is an extension of the ML criteria that were originally developed for angle estimation with some of the sensors that are perfectly calibrated, while others are uncalibrated. Recently an adaptive PSO algorithm has been successfully applied in ML optimization solutions [12]. In another report a combined fuzzy adaptive PSO and differential evolution are also used to solve economic dispatch problems [13]. In this paper, a modified adaptive inertia weight PSO (APSO) is proposed to rationally balance the global exploration and local exploitation abilities of the particle to the GMUSIC criteria function [3] for CDMA signals. The resulting estimator is called APSO-GMUSIC, where a simple and effective measure, individual search ability, is defined to indicate whether each particle lacks global exploration ability in each dimension. By employing APSO to treat our optimization problem of the desired signal direction angle . The DOA, obtained by the APSO-GMUSIC estimator, converges to the optimal or near optimal solution. Simulation results show that the proposed APSO-GMUSIC estimator is very suitable to treat the DOA estimation for CDMA signals in a local scattered scenario.

Array Data Model
Consider a DOA estimation scenario in a CDMA system with P users. The data bit b p (t) ∈ {−1, 1} of the pth user is spread by a pseudo-noise (PN) signature waveform c p (t) for p = 1, 2, …, P. The signature waveform of the pth user is composed of a spreading sequence of L chips, i.e., , where is the lth spreading chip for the pth user and p Tc (t) is the spreading pulse of duration . Thus, the transmitted signal of the pth user during a bit interval T b can be represented by: where ε p is the signal amplitude. Let the processing gain of the spreading code be L = T b /T c . To simplify analysis, the PN code during a bit interval at a sample rate 1/T c can be expressed as an L × 1 vector of discrete sequence c p = [c p1 , c p2 ,…, c pL ] T where the superscript T is the transpose operator. Employing the code-matched filter and DOA estimation algorithm, we propose a receiving base station of CDMA communication system with a uniform linear array (ULA). The antenna array consists of M identical omidirectional elements and receives line-of-sight signals from P users that arrive at the array from different bearing angles θ 1 , θ 2 , …, θ p with respect to the broadside of array in a cell/sector. Assume that the array element space is d. The direction vector associated with the pth user is given by: where a m (θ p ) = exp[-jπd(m-1)sin(θ p )/β] denotes the response of the mth sensor array to a signal with unit amplitude arriving from the direction angle θ p , where β is the propagation speed of plane wave. Thus, the received baseband signal across the array can be written in vector form: where , τ p is the time delay, and the observation noise vector n(t) represents the spatially and temporally zero mean white Gaussian noise vector.

MUSIC Algorithms for Local Scattering Channels
To show how the spatial signatures depend on the spatial distribution of the multipath propagation, we assume that the time dispersion produced by the multipath propagation is small compared to the reciprocal of the signal's bandwidth, the time delay may be approximated as a phase shift , where f c is the carrier frequency, and τ ph is the time delay associated with the hth scattered signal from the pth source arriving at the array. Let s p (t) be the signal scattered by the pth user. Due to the multipath propagation near the pth user, its contribution to the array is modeled as a superposition of N p scattered rays. Then, the data received at the array is given by M × 1 vector form: where the reflection coefficient β ph is the complex amplitude of the hth scattered signal from the pth user, a ph is the response vector of the hth scattered signal from the pth user, and . is spread angle of the pth user due to local scattering. It is assumed that the power is normalized so that all rays have equal power 1/N p . Then , where χ ph is i.i.d. and uniform over [3].
Then, the received data vector can be expressed as: where U = [u 1 , u 2 , …, u p ] and s(t) = [s 1 (t), s 2 (t), …, s p (t)] T . The pth column of U, which is denoted as u p , represents the spatial signature associated with the signal s p (t), where and . In order to pick out the pth user's signal, a code-matched filter containing the specified PN code vector c p is applied to x(t). For simplicity, the bit sampling index is dropped to make the equations more readable in the following description. Therefore, the pth user's despread and sampled array vector signal y p can be represented as: (6) where s p = ε p b p and s q = ε p b p κ qp . Hence, and denote the average power of the desired pth user and each interferer user, respectively. The notation is used to denote the expectation operator. κ qp is the inner product of two different PN code vectors c q T c p , and the M × 1 vector n p with zero mean and variance π n is the code-matched filter output of the noise vector n.
It was shown in Appendix A of [3] that κ qp approaches a real Gaussian random variable with zero mean and unit variance for the normalized PN codes. Consequently, the signals of the undesired users through the code-matched filter, i.e., MAI, can be viewed as a noise vector n MAIp with zero mean and variance π MAIp and . The correlation matrix of Equation (6) in the observation interval is given by: where the subscript H denotes conjugate transpose, π np = π MAIp + π n , and I M is an identity matrix. For high angular resolution subspace methods, the MUSIC algorithm is a kind of DOA estimation technique based on eigenvalue decomposition (EVD), which is also called the noise subspace-based method. The eigendecomposition of Equation (7) can be expressed as: (8) where are the eigenvalues of R p and e pm denotes the eigenvector associated with for m = 1, 2, …, M. Moreover, e p1 and are orthogonal and span the signal and noise subspace corresponding to R p , respectively.
is the noise eigenvalue matrix. Furthermore, e p1 spans the same signal subspace as that spanned by a(θ p ) . Thus, we have and . The MUSIC estimates the DOA of the pth user from the highest peak of the following spectrum [3]: (9) where the largest maxima of is taken to be the estimate of the DOA of the pth user.
As the angular spread of the hth scattered ray from the pth source arriving at the array is relatively small, a first-order Taylor expansion of in Equation (4) may be approximated as follows [14]: (10) where and is the first derivative of the steering vector.
Equation (10) is a linear model and regards as the generalized array manifold (GAM) model with scattering, which is the superposition of the pth source steering vectors and their derivatives. Then, the received data vector can be expressed as: Therefore, the pth user's despread and sampled array vector signal can be represented as: The correlation matrix of Equation (12) in the observation interval is given by: The eigendecomposition of Equation (13) can be expressed as:  (9) can be expanded into a MUSIC with local scattering. The MUSIC searching function derived by the steering vector of the GAM model in Equation (10) , which is termed as GMUSIC and can be expressed as: (15) where and . The parameter of ρ can be estimated by using the method of [3]. And, the in Equation (15) is the noise subspace, which is the eigendecomposition of the autocorrelation matrix derived by the linearized steering vector in Equation (10). We can estimate DOA by searching the array manifold for direction vector and the largest peak of the function GSp(θ, ρ) denotes the DOA of the pth user.

Problem Formulation
Recently, a much more computationally efficient and more accurate parameter estimating approach via polynomial rooting method has been proposed to improve the spectral searching approaches for reducing the computational load. Due to the fact the scanning vector of the GMUSIC is not Vandermonde structure, GSp(θ, ρ) cannot be implemented by using the polynomial rooting approach. The computational complexity of the classical subspace approach of MUSIC estimator is about 12M 3 complex multiplications (CM) for computing EVD of a M × M dimension matrix. Let the search number be B. Therefore, the total required CM for computing Equation (15) of the GMUSIC estimator is about 12M 3 + BM 2 . The performance of the abovementioned spectral searching approach is governed by the scanning grid size and the number of search grids while implementing the high-resolution DOA estimation. A major problem of GMUSIC type algorithms is the heavy computational load incurred with spatial spectral search when a smaller grid size is adapted. However, it is time consuming and the optimum search grid size is not clear. Smaller grid size can improve estimate accuracy, but the required computational load also becomes relatively larger (e.g., if grid sizes are set to 1°, 0.1° and 0.01°, then the required search numbers B = 181, 1,801, and 18,001, respectively.). With the number of particle population P N and the maximum number of iterations k max , the computational complexity of our APSO approach is about P N × k max . Hence, the computational complexity of APSO is smaller than the GMUSIC. Therefore, in order to reduce computational load, this paper investigates the feasibility of applying the PSO to replace the spectrum searching approach. As a result, the proposed PSO-based searching GMUSIC estimator does not increase the computational complexity, but significantly reduces the searching process requirement compared with the spectrum searching GMUSIC.

The Proposed PSO Based Searching Algorithms
Due to the fact that the performance of the abovementioned spectral searching GMUSIC estimator is governed by the scanning grid size and the number of search grids while implementing the high-resolution DOA estimation, it is time consuming and the search grid is not clear. In this section, we present PSO based searching approaches by maximizing the fitness function of GMUSIC at each iteration.

HPSO-GMUSIC Estimator
In order to reduce the scanning accurate angle problems of high computation cost, we use the PSO to replace scan in the ULA. For time-space signal processing, each particle can be treated as a point in the P-dimensional problem space and a swarm consisting of D random particles and then searches for best position (solution) by updating generations until exceeding the limit of iteration number. The position of the ith particle is represented as θ i = [θ i1 , θ i2 , …, θ iP ], the rate of the position change (velocity) for ith particle is represented as The best previous position of the ith particle, which gives the best fitness value, is recorded as p i = [p i1 , p i2 , …, p iP ]. And, the index of the best particle among all the particles in the population is represented by p g = [p g1 , p g2 , …, p gp ] and called global best location. In every iteration, the local bests and global bests are determined through evaluating the fitness values of the current population of particles. Two factors characterize a particle status on the search space: its position and velocity. The P-dimensional position for the ith particle in the kth iteration can be denoted as Similarly, the velocity for the ith particle in the kth generation can be described as . The velocity and the position of the ith particle at the (k + 1)th iteration for i = 1, 2, …, D and k = 1, 2, …, k max are updated according to the following equations: where denotes element-wise product and the positive acceleration constants c 1 and c 2 are two positive constants named learning factors. r 1 (k) and r 1 (k) are P-dimensional vectors consisting of independent random numbers uniformly distributed between 0 and 1, which are used to stochastically vary the relative pull of p i and p g in order to simulate the unpredictable component of natural swarm behavior. The inertia weight w(k) based on the linear decreasing strategy is considered critically for the convergence behavior of PSO [10], which is selected to decrease during the optimization process. Thus, PSO tends to have more global search ability at the beginning of the run while having more local search ability near the end of the optimization. Given a maximum value w max and a minimum value w min , w(k) is updated as: (18) where k max is the maximum number of iterations, w max and w min are chosen as 0.9 and 0.4 respectively in this work as in [10]. Because the PSO particle represents a series of priorities that range from −90° to 90°, all parameters of the P-dimensional particle positions, either initialized or updated during search, must be confined within [θ min , θ max ] = [−90°, 90°], avoiding infeasible particle positions that can lead to slow PSO searches. The new particle position is calculated using Equation (17).
A particle position vector is converted to a candidate solution vector in the problem space through a suitable mapping. The score of the mapped vector evaluated by a GMUSIC function is regarded as the fitness of the corresponding particle. For time-space signal processing, the DOA GMUSIC estimation problem in this paper is to find the ith particle position or estimate angle θ ip (k) of the pth user to maximize the following fitness function: (19) At the end of iteration, global best location is the final estimated angle defined as . The final global best position is taken as the GMUSIC estimates of users.
As a result, the PSO with linear decreasing inertia weight, particle velocity limitation, and particle position clipping is termed as hard-constraint PSO GMUSIC (HPSO-GMUSIC). Unfortunately, the HPSO-GMUSIC is not so facile to implement for its overhead on computing the search space mapping, particle velocity limitation and particle position clipping. In this paper, we present a more feasible and efficient modified PSO algorithm for the HPSO-GMUSIC estimator.

The Proposed Multiple Adaptive Inertia Weight
The inertia weight is critical for the performance of PSO, which balances global exploration and local exploitation abilities of the swarm. In other words, the inertia weight w(k) is employed to control the impact of the previous history of velocities on the current velocity, thereby influencing the trade-off between global and local exploration abilities of the "flying points". A large inertia weight tends to facilitate searching new area and global exploration. Conversely, a small inertia weight facilitates local exploitation in the current search area. Suitable selection of the inertia weight provides a balance between global and local exploration abilities and thus requires less iteration on average to find the optimum. During the search every particle dynamically changes its position, so every particle locates in a complex environment and faces different situation. Therefore, each particle along every dimension may have different trade off between global and local search abilities. According to [15], it has been shown that the performance of the PSO algorithm with linearly decreasing inertia weight has the ability to quickly converge, but the PSO may lack global search ability at the end of a run and may fail to find the required optima in cases when the problem to be solved is too complicated and complex. But to some extent, this can be overcome by employing a self-adapting strategy for adjusting the inertia weight. In this subsection, inertia weight is dynamically adapted for every particle along every dimension. A measure, individual search ability, which characterizes the faced situation for every particle is defined. Basing on this measure, the particle could decide to whether to increase or decrease the values of inertia weight by means of the transfer function. The fine strategy of dynamically adjusting inertia weight could lead to improvement in performance of PSO.
For velocity updating, the last two parts of Equation (16) can be viewed as the accelerations parts and can be defined as a ip (k) = c 1 ·rand(•)·(p ip (k) -θ ip (k)) and a gp (k) = c 2 ·rand(•)·(p ip (k) -θ ip (k)), where rand(•) is independent random number with uniformly distributed between [0,1]. So we could consider that the particle is moving with velocity of v ip (k) and acceleration of a ip (k) and a gp (k). But, a gp (k) is the dominant term for improving convergence rate. Suppose that the mass of the ith particle in the pth dimension is normalized to 1 kg. According to the principle of mechanics, a gp (k) = f ip (k), where f ip (k) is an outside force, which is put on the particle comes from the "pulling" of p ip (k) and p gp (k). For ULA, the DOA is a one-dimensional searching problem. In order to make the particle fly towards optimal 1 fitness ( ) , 1, 2, , region quickly, v ip (k) should turn to the direction of f ip (k) as soon as possible. Define the error between p ip (k) and θ ip (k) is z ip (k) = p gp (k) -θ ip (k). Let m ip (k) = |z ip (k)|/θ max which is a number between 0 and 1. The velocity update problem of the ith particle on the pth dimension can be divided into two classes: (1) Firstly, v ip (k) and z ip (k) are in the same direction, z ip (k)/v ip (k) ≥ 0. If m ip (k) is relatively large, it means the particle is in the right direction, but the velocity is too small. Therefore, the particle needs to speed up, and the inertia weight w ip (k) needs to be set larger. If m ip (k) is relatively small, it means the particle has come to the location, that is near the optimal region. So the velocity of this particle should slow down and the neighborhood of the state should be searched carefully, and w ip (k) should be set smaller.
(2) On the other hand, consider v ip (k) and z ip (k) on different direction, relatively large, it means that the particle's state is far from the optimal region. So the particle needs to change its velocity as soon as possible, and the inertia weight w ip (k) needs to be set smaller. If m ip (k) is relatively small, it means that it's not urgent for the particle to change its direction on this dimension, and w ip (k) could be set a large value.
Based on the aforementioned analysis, an adaptive inertia weight strategy is proposed and is shown in Table 1. The individual normalized search ability of the ith particle along the pth dimension is defined as m ip (k). It is noted that the value of w ip (k) is a function of m ip (k). To get a balance of global search and local search ability, w ip (k) cannot be too large or too small, thus w ip (k) is limited in the range of [w min , w max ], which is like the process of normalization. Therefore, we used μ-law algorithm to achieve our strategy for every particle along every dimension, normalization of the m ip (k) to .
The μ-law algorithm is a companding scheme used in telephone network [16]. Increasing the value of u, the dynamic range capability of μ-law can be improved and defined by: where m ip (k) and are normalized input and output values respectively, and the positive constant is the compression parameter. The reciprocal slop of this curve is given by the derivative of m ip (k) with respect to . For case 1, the value for small m ip (k) increases at the expense of the value for large m ip (k). To accommodate these conflicting requirements (i.e., a reasonable values for both small and large m ip (k)), a compromise is usually made in choosing the value of parameter μ for μ-law. The requirement of case 2 is opposite. Finally, based on our strategy and Equation (20), we can define multiple adaptive inertia weight w ip (k) as: In Equation (21), w min be added to avoid particles from stopping moving. The curves of w ip (k) with w min = 0 using different μ can be plotted in Figure 1. Then, we also investigate the sensitivity for APSO-GMUSIC with different values of μ However, it accords with our strategy for different µ. Thus, μ-law algorithm with μ = 100 is chosen. Note that for every particle in population, w ip (k) is unique and can be computed individually. Therefore, the single inertia weight w ip (k) can be replaced by a multiple adaptive inertia weight w ip (k). The proposed APSO-GMUSIC seems to be robust to control parameters due to the intrinsic advantages of the algorithm and the separation of the problem-independent PSO kernel from newly introduced problem-specific features in our design for adaptive multiple inertia weight. Finally, the steps for implementing the APSO-GMUSIC are shown in Figure 2 and described in the list that follows.

Computer Simulations
Several computer simulations for illustration and comparison are presented in this section. We use a 6-element ULA with half wavelength for simulation. Consider an asynchronous CDMA system with PN codes of length 31 and BPSK modulation. The spatial signature is generated by  Table 2.
All PSO algorithms start with random initializations and are terminated if the maximum iteration k max is reached or the global best particle position is not updated in 20 successive iterations. Every simulation result is presented after 200 data bits were processed and it is averaged by F = 10 3 independent Monte Carlo runs with independent noise samples for each run. Comparison results with other estimators, including the GMUSIC, HPSO-GMUSIC and APSO-GMUSIC with µ = 100 to DOA estimation error are presented. Figure 3 depicts the convergence in terms of DOA RMSE versus the number of iterations. As a result, the HPSO-GMUSIC requires more iterations to achieve convergence. Note that the proposed APSO-GMUSIC achieves fast convergence with the selected parameters, which means that it needs less iterations to approach the desired global extreme. Figure 4 shows the required number of calculating fitness function (B) versus number of particles. For the number of particles in the population, more particles may increase success in searching for optima due to sampling state space more thoroughly. However, more particles require more evaluation cost. The HPSO-GMUSIC needs more particles to approach the desired global extreme. It is confirmed that the proposed PSO-based searching approaches can reduce the computational complexity of the GMUSIC due to the searching process. As expected, this figure also provides a great improvement over the convergence rate on optimization problems. In fact, additional adaptive multiple inertia weight operation can improve the searching speed and RMSE performance further. Figure 5 presents the RMSE of DOA estimation versus varying angular spreads. We note that the subspace-based techniques show serious degradation when faced with local scatters. Local scattering may be viewed as a form of model error and gives rise to the perturbation of the noise subspace. Again, these figures show that the proposed APSO-GMUSIC method yields significantly superior performance over the other methods in the presence of local scatters. For comparison, the result of GA-GMUSIC estimator is also provided. The same parameters of GA-GMUSIC estimator are used in [17]. Figure 6 shows the RMSE versus different SNR of the desired user for angular spreads 2Δ p = 1°. For the low SNR case, all of methods may produce highly biased estimates. Clearly, with the compatible searching resolution, the APSO-GMUSIC can save the required number of searching grids and improve the RMSE performance, as compared with the other estimators. The GA-GMUSIC has a poor performance than HPSO-GMUSIC and APSO-GMUSIC. It is well known that premature convergence degrades the performance of GA and reduces the search ability [18]. In addition, it has been shown that the performance of the PSO algorithm with linearly decreasing inertia weight has the ability to quickly converge, the PSO may lack global search ability at the end of a run and may fail to find the required optima in cases when the problem to be solved is too complicated and complex [19]. But to some extent, this can be overcome by employing the proposed adaptive multiple strategy for adjusting the inertia weight. Finally, in Figure 7, we compare the RMSE performance against the number of active user P, given SNR = 20 dBW and angular spread 2Δ p = 1°. Basically, the RMSE is increased quite steadily with the increase of P. It can be observed that the APSO-GMUSIC obtain more performance improvement when the number of users P is reasonably increasing. Among them, the proposed APSO method achieves the lowest RMSE.

Conclusions
In this paper, we have presented a PSO based searching DOA estimation method, named APSO, which uses the GMUSIC searching function for CDMA signals. With the code-matched filter, the MAIs after code-decorrelation appear as noises for CDMA signals. The rooting MUSIC method is suboptimal in the presence of the noise and MAI and the GMUSIC cannot be implemented by using the polynomial rooting approach. However, the proposed techniques reduce the required search grids for the conventional spectral searching estimators. In the previous work on DOA estimation using PSO algorithm [10], hard constraints have been taken, during each iteration of PSO algorithm. In this paper, multiple inertia weight has been incorporated in PSO. In conjunction with a modified PSO for angle searching, the proposed approach can reduce the required computational load for the conventional spectral searching MUSIC estimator with local scattering. Moreover, the convergence of the proposed approach is much faster. Computer simulations have demonstrated the effectiveness of the proposed approach. Furthermore, a common drawback with MUSIC-like technique is that it is a suboptimal estimator and tends to suffer from low performance due to low SNR, small sample size, and correlated sources. Thus, how to develop a cheaper way and work well with correlated or even coherent sources will be important in the future work.