A new low-complexity angular spread estimator in the presence of line-of-sight with angular distribution selection

This article treats the problem of angular spread (AS) estimation at a base station of a macro-cellular system when a line-of-sight (LOS) is potentially present. The new low-complexity AS estimator first estimates the LOS component with a moment-based K -factor estimator. Then, it uses a look-up table (LUT) approach to estimate the mean angle of arrival (AoA) and AS. Provided that the antenna geometry allows it, the new algorithm can also benefit from a new procedure that selects the angular distribution of the received signal from a set of possible candidates. For this purpose, a nonlinear antenna configuration is required. When the angular distribution is known, any antenna structure could be used a priori; hence, we opt in this case for the simple uniform linear array (ULA). We also compare the new estimator with other low-complexity estimators, first with Spread Root-MUSIC, after we extend its applicability to nonlinear antenna array structures, then, with a recently proposed two-stage algorithm. The new AS estimator is shown, via simulations, to exhibit lower estimation error for the mean AoA and AS estimation.


I. Introduction
Smart antennas will play a major role in future wireless communications. There exist several smart antenna techniques such as beamforming, antenna diversity, and spatial multiplexing. Future smart antennas will most likely switch from one technique to another according to the channel parameters [1]. One of the most important parameters is the multipath angular spread (AS). For instance, the beamforming technique is to be considered when the AS is relatively small, while antenna diversity is more appropriate in other cases. Moreover, mean angle of arrival (AoA) and AS estimates are required to locate the mobile station [2].
In the last two decades, several algorithms have been developed for the direction of arrival and AS estimation. Based on the concept of generalization of the signal and noise subspaces, 3 multiple signal classification (MUSIC) is the most known mean AoA estimator. For AS estimation, many derivatives have been proposed. DSPE [3] and DISPARE [4] are two generalizations of the MUSIC algorithm for distributed sources. They involve maximizing cost functions that depend on the noise eigenvectors. The mentioned estimators are computationally heavy because of the required multi-dimensional systems resolution. A low-complexity subspace-based method, Spread Root-MUSIC, is presented in [5] where a rank-two model is fitted at each source, using the standard point source direction of arrival algorithm Root-MUSIC. This rank-two model depends indirectly on the parameters that can be estimated using a simple look-up table (LUT) procedure. In [6], a generalized Weighted Subspace Fitting algorithm is proposed. The latter, in contrast to DSPE and DISPARE, gives consistent estimates for a general class of full-rank data models. In [7], a subspace-based algorithm has been formulated that is applicable to the case of incoherently distributed multiple sources. In this algorithm, the total least squares (TLS) estimation of signal parameters via rotational invariance techniques (TLS-ESPRIT) approach is employed to estimate the source mean AoA. Then, the AS is estimated using the LS covariance matrix fitting. However, the performance of this algorithm shows unsatisfactory results under some practical conditions [8]. In [9], a maximum likelihood (ML) algorithm has been proposed for the localization of Gaussian distributed sources. The likelihood function is jointly maximized for all parameters of the Gaussian model. It requires the resolution of a four-dimensional (4D) nonlinear optimization problem. In [9] and [10], LS algorithms are considered to reduce the dimension of the system. The simplified ML algorithm belongs to the covariance matching estimation techniques (COMET) [11]. In [12], a low-complexity algorithm based on the concept of contrast of eigenvalues (COE) has been developed to estimate AS and mean AoA. The authors establish a bijective relationship between the COE of the covariance matrix: the signal-to-noise ratio (SNR) value and the value of the AS. Hence, for each SNR, a LUT is built. The mean AoA is derived using the estimated AS and the number of dominant eigenvalues of the source covariance matrix.
Many of these estimators make assumptions on the shape of the signal distribution, assume narrow spatial spreads, and eigen-decompose the full-rank covariance matrix into a pseudo-signal subspace and a pseudonoise subspace. Most often they result into a multidimensional optimization problem, implying high computational loads.
To overcome this limitation, a low-complexity estimator [5] has been developed. Spread Root-MUSIC consists in a 2D search using the Root-MUSIC algorithm. Another mean AoA and AS estimator based on the same approach as Spread Root-MUSIC was developed in [13]. Indeed, thanks to Taylor series expansions, the estimation of AoA and AS is transformed into a localization of two closely, equi-powered and uncorrelated rays. However, like other estimators, Spread Root-MUSIC considers scenarios without line-of-sight (LOS). A new low-complexity estimator, based on a LUT approach was therefore developed [14]. First, it estimates the LOS component of the Rician correlation coefficient and deduces the Non-LOS (NLOS) component. Then, it extracts the desired parameters from LUTs computed off-line. The new estimator, like most estimators, assumes the a priori knowledge of the angular distribution of the received signal. In this article, we enable this method to select the angular distribution type from a set of possible candidates. For this purpose, a nonlinear array structure is required. We also compare the new technique to other low-complexity AS estimators. The first one is derived by extending the Spread Root-MUSIC algorithm [5] to the considered antenna configuration. The second one is the two-stage approach developed in [13].
The article is organized as follows. In Section 2, we def ne the used notations and describe the data model. In Section 3, we describe the new method for selecting the angular distribution type. Section 4 details the two low-complexity AS estimation methods that will be used to benchmark our newly proposed approach, that is the Spread Root-MUSIC algorithm [5], modified to handle a nonlinear array structure, and the two-stage approach presented in [13]. In Section 5, simulation 5 results are presented and discussed.

II. Notations and data model
In this article, non-bold letters denote scalars. Lowercase bold letters represent vectors. Uppercase bold letters represent matrices. The row-column notation is used for the subscripts of matrix elements. For example, R is a matrix and R ik is the element of that matrix on the ith row and the kth column. The sign ∧ . denotes an estimate. Superscripts between parenthesis are used to differentiate estimates at different stages of the estimation process.
In this article, we consider the single input-multiple output (SIMO) model for the uplink (mobile to base station) transmission. The mobile has a single isotropic antenna surrounded by scatterers. We also assume that the base station is located high enough and far from the mobile to ensure 2D AoAs and to avoid local scattering shadowing. As one example, these conditions are observed in the current GSM and 3G networks where the base station is usually placed on the building roofs. As in [14,5,15,7], we suppose that the base station antenna-elements are isotropic and that the same mean AoA and AS are seen at all antenna-elements of the base station.
We consider the estimation of the AS and mean AoA from estimates over time of the time-varying channel coefficients associated with a single time-differentiable path at the multiple elements of an antenna array. Our model can therefore be associated with a narrowband channel, or with a given time-differentiable path of a wideband channel. Of course, in a wideband channel scenario, the potential presence of a LOS would only be considered for the first time-differentiable path, and knowledge of a zero K-factor could be assumed for the rest of the paths. We consider the following expression for the Rician channel coefficient [16]: where a i is associated to the channel coefficients of the diffuse component (Rayleigh channel) for antennaelement i, Ω is the power of the received signal, K is the Rician factor, F d and g d , are respectively, the Doppler frequency and Doppler angle. l is the wavelength and d 0i is the distance between the antenna reference and the antenna-element i, and θ 0i is the AoA of the LOS, as shown in Figure 1a. Indeed, in our model, we consider uniform clusters, so that the mean AoA corresponds to the AoA of the LOS. Let x i be where x i (n) =x i (nT s ), and T s is the sampling interval. In this study, we consider an arbitrary array geometry. That is why the array model described for instance in [3] and [11] is not adopted here a . Instead, we use the correlation coefficient of the Rician channel coefficients received at the antenna branch (i, k) given by where (.) H is the transconjugate operator. Hence, the Rician correlation matrix associated with the coefficients,RT ik , would be where M is a square matrix defined by m ik = d oi λ sin(θ 0i ) − d 0k λ sin(θ 0k ). The expression for the correlation coefficient of the diffuse component (Rayleigh channel) is [17] where • θ ik is the mean AoA; • σ θ ik is the AS or the standard deviation of the angular distribution; • f c is the carrier frequency; • c is the speed of light; • d ik is the distance between the antenna-element i and the antenna-element k; and  • the function f (θ , θ ik , σ θ ik ) is the power density function with respect to the azimuth AoA θ.
In this article, we consider only the Gaussian and Laplacian angular distributions, the most popular ones in the literature. However, our approach is still valid with other angular distributions.
If we consider the diffuse component and we assume a small AS value (s θ < s small ), then the correlation coefficient R i, k would be [14,18] • Gaussian distribution: • Laplacian distribution: In this study, we are interested in estimating the mean AoA and the AS. In other terms, we determine the mean and the standard deviation of the angular ditribution of the received signal. The proposed algorithm is valid for non linear antenna arrays. Hence, each antenna branch represents different mean AoA and AS estimation values. That is why the parameters in question are function of the indexes i and k which refer to the associated antenna pair (i, k), as shown in Figure 2. As noticed, the two pairs (i, k) and (k, l) represent different mean AoA and AS values, (θ ik , σ θ ik ) and (θ kl , σ θ kl ). Each couple is estimated using the correlation coefficients, R i, k and R k, l, respectively. This model formulation with global parameters can be advantageous in a parameter estimation framework, when evaluating the Cramér Rao bound (CRB), for instance. In the following, we develop a new mean AoA and an AS estimator based on the correlation coefficient defined in (3).

III. New estimator with angular distribution selection
The idea is to find a simple relationship between the mean AoA and AS, and the Rician correlation coefficient. Since the expression of the Rician correlation coefficient R T ik is complex, our approach is to estimate the LOS component first. Then, the diffused component R ik is deduced, and the AS is extracted from LUTs. For each angular distribution type, a LUT is built off-line using the expression (5) for the NLOS component of the correlation coefficient. Indeed, for all possible values of the mean AoA and AS, the correlation coefficient of the diffuse component is computed using a numerical method (5). In our simulations, we varied the mean AoA from 0 to 90 degrees with a step of 0.1 degree. The AS is varied from 0 to 100 degrees with a step of 0.025 for small ASs (s θ <6 degrees) and a step of 0.1 degree for higher ones. One can argue that the building of the LUT using the considered steps requires a lot of time and an accurate resolution of the integral in (5). However, the LUT is computed once for all off-line and would not affect the realworld execution time of the new algorithm. Besides larger steps would affect the accuracy of the new estimator. The LUT expresses the desired parameter b as a function of the magnitude and phase of the diffuse component R ik . As defined in (4), the LOS component depends only on the Rician K-factor and the AoA of the LOS. In this study, we consider uniform clusters. Hence, the AoA of the LOS coincides with the mean AoA. If we assume small AS values and consider the diffuse component of the correlation coefficient (6) associated to the Gaussian distribution, then the relationship in (4) becomes Considering only antenna-element pairs including the antenna-element reference "0″, both terms of the  correlation coefficient R T 0,k admit the same argument: Hence, the mean AoA is estimated by using the phase of the correlation coefficient associated to the antenna pairs (0, k). By analogy, the same expression is obtained for the Laplacian distribution: where ∠ symbolizes the phase operator and the subscript "0″ refers to the antenna-element reference and the distance separating the antenna-element pair (0, k) is such that d 0k ≈ λ 2 . As one can notice, we use only the antenna-elements pair (0, k) to estimate the AoA LOS. Otherwise, the correlation coefficient of the diffuse component, R i, k , and the correlation coefficient of the LOS component would admit different arguments (see (8)). The final mean AoA estimate,θ m , is the mean ofθ 0k over all antenna-elements pairs {(0, k)} spaced by λ 2 . Indeed, the antenna pairs spaced by d 0k ≫ l give high estimation error since the correlation coefficient does not contain enough information, i.e., R T 0,k is close to zero. It is understood that (10) is valid for antenna configurations having at least two antenna-element spaced by λ 2 . In most references, ULAs spaced by λ 2 are considered. Hence, our condition enlarges the set of possible antenna arrays that can be used. One can argue that 10 this solution does not take into account the left-right ambiguity. Indeed for linear arrays (antenna-element pairs in our case), it is not obvious to determine whether the incident signal is coming from the left side or the right one of the array [19,20]. To avoid this ambiguity, we divide the cell into three or more sectors and the mean AoA estimation is achieved in each sector. In the remainder of this article, (10) is used for antenna-element pairs for which the left-right ambiguity does not arise. In other words, we imply that the arrays are constructed in a way that prevents this ambiguity by considering the cell division approach or other methods as in [19]. Indeed, this condition limits the subset of antenna structures that can be used for the mean AoA estimation, but still allows some flexibility in the design of antenna arrays. Without loss of generality, we consider the antenna configurations illustrated in Figure 1. All structures are supposed to be constructed in a way that prevents the left-right ambiguity. For these symmetrical structures, after a simple mathematical manipulation c , it is observed that (10) is also true for correlation coefficientsR T i,k associated with antennaelement pairs (i, k) spaced by d ik ≈ λ 2 , i.e., the antenna pairs (0,1) and (1,2) of all structures presented in Figure 1. The angles must have the same reference which in this case the normal to the antenna structure, and the clockwise sense is the positive one. The AS estimation is not affected by the choice of the angles measurement reference. Indeed, it measures the angular distribution spreading around the mean AoA. One can argue that relation (10) is only valid for small AS values assumption. However, we empirically find that the mean AoA estimate using (10) is still accurate for high AS values.
For the Rician factor, many K-factor estimators have been developed. In [21], the Kolmogorov-Smirnov statistic is used first to test the envelope of the fading signal for Rician statistics and then to estimate the K-factor. In [22], the K-factor estimator is based on statistics of the instantaneous frequency (IF) of the received signal at the mobile station. In [23], ML estimators that only use samples of both the fading envelope and the fading phase are derived. In [24] and [25], a general class of moment-based estimators which uses the signal envelope is proposed. A K-factor estimator that relies on the in-phase and quadrature phase components of the received signal is also introduced in [24].
We choose to consider the closed-form estimator presented in [24], which is easily implemented and quite accurate. This estimator uses the second-order and fourth-order moments (μ 2 and μ 4 ) of the received signal to estimate the K-factor (shown here for the estimate on the ith antenna): The final K-factor estimate is the mean ofK i over all antenna-elements i.
In [26], the expressions for the second-order and fourth-order moments at antenna-element i are: where Ω and N 0 are, respectively, the signal and the noise powers; and k i;a and k i;ω are, respectively, the Rician and noise kurtosis. In our article, we consider an additive white Gaussian noise (AWGN), i.e., k i;ω = 2, [26]. As in [14], we consider an estimated SNR to reduce the noise bias. The expressions of the secondorder and fourth-order moments becomê In the literature, the value of k i;a is computed by using the Rician K-factor which is unknown at this stage [27]. In our procedure, the Rician kurtosisk i;a is obtained by analyzing the term 2 and is computed as follows: Several SNR estimators can be considered, such as in [28,29] or [30]. In this article, we do not consider a specific SNR estimator, to avoid restricting our algorithm to a particular SNR estimator results. Instead, we consider an estimated SNR expressed in dB, (SNR dB ). The latter is characterized by an estimation error modeled as a zero-mean normally distributed random variable with variance . As shown in [28], the studied estimators offer low estimation errors, especially for long observation windows. For the SNR range considered in our simulations, the variance of the estimation error is around 10 -1 . Hence, we choose extreme cases where σ 2 ε = 0 or 1 (i.e., optimistic and pessimistic bounds).
With AWGN bias reduction, the expression for the estimated Rician correlation coefficient (for i ≠ k) iŝ Once the AoA of the LOS and the Rician K-factor are estimated, the estimated NLOS componentR ik is then deduced: When the antenna-elements separation d 0k > λ 2 , we takeθ ok =θ m . Note that all angle measurements must have a common reference. The AS is extracted from the LUT associated to the considered angular distribution type. Using linear interpolation, we determine which AS value corresponds to the magnitude and phase of the estimated correlation coefficientR i,k . In this article, we treat the case when the a priori knowledge of the angular distribution is assumed. In this case, arbitrary arrays can be used including ULAs. We also propose a new method to determine the angular distribution of the received signal when it is unknown. In this case, a nonlinear array is required. In fact, we select the angular distribution type that fits the array geometry from a set of possible candidates. Different mean AoAs and ASs are obtained for the different antenna branches which is not the case for linear structures. Then, the selected angular distribution is the one associated with the minimum of the estimates' standard deviations. The level of accuracy for small AS values is taken into account as well. Indeed, (6) and (7) are computed assuming small AS values. As a result, we must first rank the AS. Then, if the latter is low, we can apply (6) or (7). In fact, the LUT approach shows low accuracy for small ASs. That is why we present here four variants of the new AS estimator depending on the knowledge of the angular distribution and the desired accuracy of the AS estimation.

A. Known angular distribution type and low AS estimation accuracy for small AS values
Let us first study a simple case. Consider a pair of antenna-elements (0, 1) spaced by d 01 ≈ λ 2 . We first estimate the LOS component, i.e., the K-factor and the mean AoA (using (10)). Owing to the estimated NLOS componentR 0,1 , we obtain the AS estimate,σ (c) 01 from the LUT associated with the considered distribution type, g. The procedure is summarized as follows: whereμ 2;i andμ 4;i are the estimated second-order and fourth-order moments, respectively. When the antenna array is composed of more than two elements, the procedure is applied to each pair. The final AS estimateσ R is the mean of allσ (c) ik divided by (K+1). The division by (K+1) is employed to recover the NLOS AS from the one associated to the Rician channel. When a uniform linear array (ULA) is used, the estimation error can be reduced even more by averaging the correlation coefficients over all antenna pairs spaced by the same distance, before using the LUT, instead of averaging the individual AS estimates over all antenna pairs.

B. Unknown angular distribution type and low AS estimation accuracy for small AS values
In real scenarios, the distribution type might not be predictable. To the best of our knowledge, there is no existing procedure that finds out the angular distribution type of the received signal. We present here a new method to select the distribution type from a set of possible candidates. The idea is to seek the distribution type that best fits the geometry of the array. A nonlinear array structure such as the one illustrated in Figure 1b, where the closest antenna-elements are spaced by ≈ λ 2 or less, has to be considered. The angle value is not static and can be modified to fit the base station construction constraints. As in the previous case, we estimate the LOS component. Then, for each antenna pair (spaced by d ≈ λ 2 ) and each distribution type g, the mean AoAθ ik (γ ) and ASσ (c) ik (γ ) are extracted from LUT g . In this case, the estimated mean AoA is actually the sum of the received signal mean AoA, θ ik , and the angle of nonlinearity ± . Hence, to recover the desired mean AoA, we substract the angle ± (see the algorithm below). The selected distribution type (γ ) is the one associated to the minimum of the sum of the standard deviations The chosen criterion is motivated by the nonlinearity of the array. For instance, using the configuration illustrated in Figure 2, the mean AoA impinging at the pair (i, k), θ ik -, must be close to the one associated to the branch (k, l), θ kl + . Indeed, we assumed the same parameter values at the different array elements. However, by considering the wrong distribution type, the obtained mean AoAs would be different and as a result show high standard deviation. The same reasoning is adopted for the ASs estimates (19). One can argue that the mean of the AS estimates could be used instead of the standard deviation in (19). Actually, the mean of the obtained estimates would not give us any information about the angular distribution of the received signal. For instance, for the array structure illustrated in Figure 2b, we obtain two AS estimates associated to the Gaussian and Laplacian distributions. In this case, we cannot select the right angular distribution. This is why we consider the standard deviation of the AS estimates.
For the mean AoA estimation, we no longer require antenna pairs including the antenna reference "0", but instead each pair (i, k) separated by l/2. Indeed, once the LOS component is determined and the diffuse component is deduced, we use (6) or (7) to estimate the mean AoA. The considered expressions are not restricted to antenna pairs (0, i) but to all antenna-elements (i, k).
Note that the procedure above estimates the mean AoA twice. In (10), the resulted mean AoA is used to compute the LOS component. Mean AoAs are then extracted from a LUT using the diffuse component of the Rician correlation coefficient. These estimates are employed to select the angular distribution by comparing their standard deviations (18). One can argue that the standard deviations of the AS could be used instead of estimating the mean AoA twice. However, when the AS is small, the angular distribution is close to an impulse function for both distribution types. In fact, the mean AoA standard deviations bear more information concerning the distribution type. In this case, the distribution type selection using the criterion (20) is no longer due to its high error rate. To overcome this limitation, we look for weights that express the importance of one parameter compared to the other, i.e., weights that ensure better selection. After running exhaustive simulations, results show that when the AS is small (s < s threshold ), only the standard deviation of the mean AoA estimates (18) must be considered. Even when the AS is high, the two standard deviations, s aoa and s as , should not be considered with the same importance. Indeed, a larger weight should be affected to the information provided by the standard deviation of the mean AoA estimates. Hence, the optimal weights were empirically set equal to whereσ (c) (γ ) is the mean of the estimated AS associated with the g th angular distribution and c is the function defined by s threshold was set empirically to 6°. In fact, this value depends also on the distribution type. In other words, σ γ 1 threshold = σ γ 2 threshold . In this study, we consider the same s threshold for both considered distribution types, to simplify implementation. Still, for more accuracy, one could use different values for each distribution type. The selected angular distribution type is then the one associated with the minimum of the weighted sum, so that, instead of (20), the following is used: The final estimates for the mean AoA and AS is the mean of the obtained estimates associated with the selected angular distribution. The overall procedure is as follows:

C. Known angular distribution type and high AS estimation accuracy for small AS values
With small AS values, closed forms can be deduced from (6) and (7): • Gaussian distribution: • Laplacian distribution: The mean AoA has the same expression as in (27), and The analysis of (28) and (29) shows that, when the correlation coefficient amplitude is close to one or zero, the AS estimation error is higher. Indeed, in this case, the estimation error of the correlation coefficient has an important impact on the AS estimation. The solution to their problem is to consider distant antenna-elements spaced by d ≫ l. Indeed, in this case, the correlation coefficient amplitude is reduced. To avoid correlation coefficients with a magnitude too close to zero, we set empirically (i.e., by running several simulations) a lower limit of 0.05 to decrease the estimation error. In other terms, we exploit only distant antenna-element pairs for which the correlation coefficient magnitude is higher than 0.05.
To illustrate the overall AS estimation process, we consider the array configuration illustrated in Figure 1c. In this section, we consider the a priori knowledge of the angular distribution type. The procedure is then as follows. After estimating the LOS component and deducing the diffuse one, we consider first the closest pair of antenna-elements (Ant.0-Ant.1). From the 2-D LUT, we estimate the ASσ (c) 01 . If the obtained preliminary AS is larger than s small , (28) and (29) are not to be considered, and the procedure is terminated. Otherwise, we use the distant elements (Ant.1-Ant.2) and the closedforms provided by (28) and (29) to estimate the ASσ (d) 12 . The overall AS estimation procedure is as follows: Ifσ (c) (k) < σ small and |R 1,2 | > 0.05 The function g refers to (28) or (29) Elsê

D. Unknown angular distribution type and high AS estimation accuracy for small AS, values
This case is a mix between the two previous cases, when an accurate estimation is needed for small AS and the angular distribution type is unknown. As one can conclude, the array structure has to have two main properties: the nonlinearity for the distribution type selection and the existence of distant antenna-elements for a high estimation accuracy in the case of small AS. The butterfly configuration presented in Figure 1d is then considered as an example. Other structures satisfying the conditions mentioned above could be used. Once the LOS component is estimated and the diffuse one is deduced, as in the previous cases, we consider first the closest antenna-elements (spaced by~λ 2 ). For each angular distribution g and antenna pair (i, k) spaced by about λ 2 , we extract the associated mean AoAθ ik (γ ) and ASσ (c) ik (γ ) from LUT g . Then, we compute the associated standard deviations, s aoa (g) and s as (g). The selected distribution typeγ is the one associated to the minimum of the weighted sum [see (24)]. If the preliminary AS σ (c) (γ ) = mean (σ (c) ik (γ )) is larger than s small , then closed forms of (28) and (29) are not to be considered, and the procedure is terminated. Otherwise, for accurate AS estimation, we consider the distant antenna-elements (d ≫ l) with a correlation coefficient amplitude higher than 0.05. The latter is chosen empirically after several simulations. A correlation coefficient with a lower module would not have enough information to allow the AS estimation. Then, for each considered pair and each distribution type, we estimate the ASσ (d) ik (γ ) using (28) and (29). During simulations, we noticed that the standard deviations of the AS estimates obtained using distant antenna-elements offer lower error probability of distribution type selection. A second angular distribution selection is therefore considered, for which at least two AS estimates (σ (d) ik (γ )) are needed. If the number of correlation coefficients with a module higher than 0.05 is inferior to 2, then we cannot compute the standard deviation of one AS estimate. In this case, the procedure is terminated, and the final AS is the preliminary AS associated with the selected distribution typeγ. Otherwise, we compute the standard deviations of the estimated AS obtained using the distant antenna-elements (σ The selected angular distribution,γ f , is the one associated with the minimum of the sum (24) (using the standard deviations of the AS estimates of distant elements). The final mean AoA estimate,θ m , is then the mean AoA associated with the selected distribution type,γ f . The final AS estimate is the mean of the AS estimates over distant antenna pairs associated withγ f , i. e., the estimated AS iŝ The overall AS estimation procedure is summarized as follows: The function g refers to (28) or (29) σ as (γ ) = std({σ (d)

IV. Other as estimation methods selected for performance comparison
In this article, we compare the new AS estimator to other low-complexity algorithms. As mentioned before, there exist more robust estimators [2] and [10], but our purpose is to evaluate low complexity estimators. Spread Root-MUSIC is therefore the appropriate candidate for performance evaluation and comparisons. We also compare the new AS estimator with the two-stage approach [13] which is based on the Spread Root-MUSIC principle.

A. Extended spread Root-MUSIC to 2-D arbitrary arrays
The principle is to localize two rays symmetrically positioned around the nominal AoA. Then, the AS is estimated by using a LUT symbolized by the function Λ(s ω ) computed by considering a mean AoA θ m = 0.
In [14], Spread Root-MUSIC for a ULA with inter-element spacing d = λ 2 , in the presence of a LOS, is presented as follows: is the estimated covariance matrix, andν i = 2π d λ sin(θ i ) is the spatial frequency, θ m is the mean AoA,σ R is the AS of the Rician fading channel, and Λ K (s ω ) is the function defined by Note that the function Λ K (s ω ) depends on the Rician K-factor. To reduce the potentially large number of LUTs, we propose to consider the estimation and the extraction of the LOS component for Spread Root-MUSIC, as does our new estimator. In other words, we consider only the NLOS component instead of considering the total estimated covariance matrixR c . In this case, one function Λ(s ω ), with K = 0, is considered. The relationship in (37) becomeŝ As presented in [5], Spread Root-MUSIC estimates the AS and mean AoA for ULAs. In this article, we adapt Spread Root-MUSIC to the butterfly configuration to be able to evaluate the performance of the new method. In [31], those authors propose an extension of Root-MUSIC to 2-D arbitrary arrays. Since Root-MUSIC exploits the Vandermonde structure of the steering vector of ULAs, the idea is to rewrite the steering vector a of nonlinear arrays as the product of a Vandermonde structured vector d and a matrix G characterizing the antenna configuration (manifold separation) [31]: The matrix G can be determined using the least square (LS) method as follows: Once the characteristic matrix is built, the MUSICspectrum is then rewritten as a function of the new steering vector: where E n is the matrix containing the eigenvectors associated to the noise subspace. As is noticed, the new noise subspace is no longer defined by the eigenvectors of the covariance matrix associated to the smallest eigenvalues, but by the product of the characteristic matrix and the eigenvectors E n . The estimated AoAs are then the arguments of the complex roots of the obtained pseudospectrum. One drawback of the extended Root-MUSIC algorithm is the heavy computations of the pseudo-spectrum. For instance, in our case, to ensure the required accuracy, the dimension of the characteristic matrix is set to (360 × 151), thereby increasing the algorithm's complexity significantly.
The modified Root-MUSIC still fulfills the properties of a consistent estimator. Hence, we can apply the Spread F algorithm described in [5]. Our new extended spread Root-MUSIC (ESRM) algorithm can be applied as follows: As noticed, the extended method does not differ from the original Spread F algorithm. The major difference is observed in the Root-MUSIC approach. The extended version of Spread Root-MUSIC computes directly the physical anglesθ i , i = 1, 2 instead of the spatial frequenciesv i , i = 1, 2, (34). The mean AoAθ m is then the mean of the estimated AoAs.
In Figure 3, we compare the functions Λ(s ω ) obtained for both ULA and butterfly configurations (see Figure  1d). Since our nonlinear structure presents distant antenna-elements spaced by 3l, we consider a ULA with seven elements spaced by λ 2 . Note that while a ULA configuration can estimate AS values higher than 10°, the functions Λ(s ω ) of the butterfly and (0-1-3) structures show lower limits.

B. Angular distribution selection using ESRM
As mentioned before, the new AS estimator selects the angular distribution according to the standard deviations of the mean AoA and AS estimates. In this article, we adopt the same approach for the ESRM. As in our method, we need two symmetric structures with at least a 3D correlation matrix. Indeed, ESRM will require two dimensions for the signal subspace and another one for the noise subspace. In this case, one can consider the following antenna-elements combinations (see Figure  1d: (Ant.0-1-3) and (Ant. 2-1-4). To build the function Λ(s ω ), Root-MUSIC requires symmetric structures. In our case, the considered array structures give the following results: where R c is the covariance matrix. We then consider the function Λ(s ω ) define by The resulting function keeps the same properties as for symmetric arrays.
Once we build the function Λ(s ω ), we apply the Spread F algorithm [5] to estimate the AS. At this step, for each angular distribution type (g), two AS estimates are obtained. The firs one, s 013 (g), is associated to the structure Ant. (0-1-3) and the second, s 214 (g), is associated to the structure Ant. (2-1-4), see Figure 1d. The selected angular distribution type (γ ) is the one associated with the minimum of the standard deviations of the AS estimates: Simulations show that the angular distribution selection using ESRM is irrelevant. Whether the distribution type is Gaussian or Laplacian, the same mean AoA is observed. Indeed, (43) does not require the knowledge of the angular distribution type. That is why the standard deviation of the mean AoA estimates is not used in (48). Moreover, as shown in Figure 3, ESRM cannot estimate an AS higher than 4°using the structures (0-1-3) and (2-1-4). Using the same AS value in both structures, the angular distribution selection cannot be achieved considering the standard deviation of AS estimates. These cases are marked "x" in Table 1. This is why the a priori knowledge of the angular distribution type is required for ESRM.

C. The two-stage approach
In [13], a new two-stage approach similar to Spread Root-MUSIC was presented. The estimation of the mean AoA and the AS of the scattered source is transformed there into the localization of two closely spaced point sources. The new approach approximates the function Λ(s ω ) by a linear function. Indeed for a ULA with M antenna-elements spaced by d = λ 2 , Λ(s ω ) ≈ s ω . As noticed, the function Λ(s ω ) no longer depends on the angular distribution type. In other terms, the selection of the distribution type is not required for the twostage approach. The algorithm is as follows: where ∠ (.) and R(.) represent operators that extract, respectively, the angle and the real parts. As other estimators, the approach described in [13] considers only LOS-free scenarios. In this article, we consider, as for ESRM and the new estimator, the NLOS component of the covariance matrix. The method exploits the Toeplitz structure of the covariance matrix by averaging the coefficient of the mth subdiagonals ofR c . It was shown in [13] that the covariance coefficients on the first subdiagonals give better mean AoA estimates. Simulations in [13] showed also that antenna-elements spaced by 2.5l offer better AS estimation. For the butterfly configuration, the algorithm described above can be applied by considering antenna pairs. In other terms, the antennaelements spaced by d = λ 2 are utilized to estimate the mean AoA and the distant ones are considered for the AS estimates. In this article, we select the antenna pairs spaced by 3l, the closest distance to the one used in [13]. Indeed, the pairs (ant.1-ant.3) and (ant.1-ant.4) can be modeled by a ULA composed by six antenna-elements. In this case, the algorithm of the two-stage approach is applied with m = 5.

V. Simulation results
We illustrate the performance of the new AS estimator by means of Monte-Carlo simulations. We assume here that the channel coefficients are obtained through an appropriate channel estimation algorithm, and that the resulting time-varying channel coefficient estimates can be adequately modeled by the sum of the true timevarying channel coefficients with an AWGN component. The accuracy of the channel estimation procedure is then controlled by the variance of the AWGN component. In our simulations, for the diffuse component, we used a non-selective frequency (f at) Rayleigh channel. We considered the Rayleigh channel simulator described in [32]. The azimuth AS distribution for the incoming multipath signals are of Gaussian or Laplacian type. The carrier frequency was set to 1.9 GHz, which results in a wavelength l of 15.79 cm. The mobile speed was set to 80 Km/h (22.2 m/s), which results in a Doppler frequency F D of 140.74 Hz. The sampling interval was set to T s = 1 1500 ms. The SNR of the estimated channel coefficients is 15 dB.
First, we study the new estimator performance when the angular distribution is known. In this case, we consider a ULA with five antenna-elements spaced by a half wavelength. We compare the algorithm described in III. C. with the weighted least square (WLS) method and the stochastic CRB developed in [10]. In this article, we consider the unknown parameter vector defined as η = [Sσ 2 n σ θ θ m ], where S and σ 2 n are the transmitted signal and noise powers. In [10], the variance σ 2 θ is considered instead of s θ . Thereby, we recompute the CRB using (13) of [10]. The normalized mean square error (NRMSE) is used to evaluate both estimators: As shown in Figure 4, for mean AoA estimation the new estimator presents lower NRMSE and is the closest one to the CRB. For the AS estimation, the new estimator and the WLS method present close NRMSE (see Figure 5), for different computational complexities. For instance, the Gauss-Newton method used in the WLS technique converges in around 600 iterations, which increases significantly its computational complexity. One can also notice that the SNR estimation error does not affect the mean AoA and AS estimation. Indeed, whether the SNR is assumed known or with a Gaussian estimation error with variance σ 2 ω = 1, both estimators show the same results.
Second, we study the angular distribution selection and its impact on the new estimator. To this purpose, we consider the butterfly configuration with j = 10°, d 01 = d 12 = λ 2 , and d 13 = d 14 = 3l. Since it does not allow the proper selection of the angular distribution type, the a priori knowledge of the angular distribution type is assumed for ESRM. To study the effect of the Kfactor estimation error and the variance of the estimated σ 2 ε , σ 2 ε , on AS estimation, we consider several scenarios. In the first one, we assume the a priori knowledge of the K-factor, i.e., we use the true value of the K-factor to estimate the diffuse component. In the second case, we consider the true value of the SNR. In the last two cases, an estimated SNR with variance σ 2 ε = 1 is used. As shown in Figure 6, the new estimator shows lower NRMSE for the mean AoA estimation. In Figure 7, the AS estimation with the variation of the K-factor value is studied. As noticed, the new estimator presents high NRMSE when the SNR is assumed unknown. As mentioned before, this is due to the important estimation error exhibited by the K-factor estimator (11). Indeed, If we consider the true value of the K-factor, the new estimator achieves more accurate results. The ESRM and the two-stage approach also show lower NRMSE when the true K-factor is assumed, especially when its value is high.
One can argue that the new estimator fails in the case of unknown SNR. Note that while ESRM requires the eigen-decomposition of the covariance matrix and finding the roots of a polynomial, our method uses only LUTs, simple closed forms, and some logical operations. Indeed, the Spread Root-MUSIC shows high complexity around M 3 log(M) + M 2 (N a + T) + N 2 a N floating point operations, whereas the new estimator and the twostage approach admit almost the same complexity of   floating point operations. Moreover, owing to the definition of the function Λ(s ω ), ESRM cannot estimate an AS greater than a certain limit. Therefore, for a large AS, ESRM exhibits high NRMSE. The two-stage approach, similar to Spread Root-MUSIC, with a linear function Λ(s ω ), is then considered. For the AS estimation (see Figures 8,9), the new estimator shows lower NRMSE than ESRM and the two-stage approach, as expected. In fact, for a high AS, the function Λ(s ω ) (see Figure 3) is not quite linear. Hence, the accuracy of AS estimation is affected. In contrast, the new estimator shows lower NRMSE for large AS values since unlike the ESRM or the two-stage approach.
As shown in Figure 10, for different values of the AS, the new estimator achieves better results then the ESRM and the two-stage approach. However, for low SNR values, the new estimator shows higher NRMSE then the ESRM, when an estimate of the SNR is considered. When the SNR is assumed known, the new estimator shows the best results (see Figure 11).

VI. Conclusion
In this article, we described a new low-complexity AS estimator for Rician fading channels. The new estimator first estimates the LOS component of the correlation coefficient. Then, the desired parameters are extracted from LUTs computed off-line. The estimate of the LOS component of the correlation coefficient requires the use of a K-factor estimator. The second-and fourthorder moments K-factor estimator is considered for its simplicity and relatively good accuracy. To reduce the   impact of the K-factor estimation error on AS estimation, the noise-induced biases in both the correlation coefficient and the moments of the received signal are reduced using an estimated SNR. The new estimator also includes a new method to select the angular distribution type of the received signal, which requires the use of a nonlinear array structure. The performance of the new method was compared with Spread Root-MUSIC extended to a nonlinear antenna array configuration and with the two-stage approach presented in [13]. Simulations showed that the new technique gives lower NRMSE.

VII. Competing interests
The authors declare that they have no competing interests.

VIII. End notes
a The conventional array model described in [3] and [11] can be extended to a nonlinear structure by rewriting the associated steering vector. In this case, our problem can be reformulated using a matrix representation. However, this would only complicate the new algorithm by adding a new step for the determination of the steering vector. That is why we consider the correlation coefficient of each antenna branch instead of the array formulation. b For each angular distribution type, there are two LUTs. The first is for the mean AoA, and the second is for the AS. c We rewrite the correlation coefficient for the different antenna pairs as in (8).