SIS-BAM Algorithm for Angular Parameter Estimation of 2-D Incoherently Distributed Sources

For two-dimensional (2-D) incoherently distributed sources, this paper presents an effective angular parameter estimation method based on shift invariant structure (SIS) of the beamspace array manifold (BAM), named as SIS-BAM algorithm. In the proposed method, a shift invariance structure of the observed vectors is established utilizing a generalized array manifold of a uniform linear orthogonal array. Then, based on Fourier basis vectors and the SIS, a beamspace transformation matrix can be obtained. It projects received signals into the corresponding beamspace, so as to carry out dimension reduction of observed signals in beamspace domain. Finally, according to the SIS of beamspace observed vectors, the closed form solutions of the nominal azimuth and elevation are derived. Compared with the previous works, the presented SIS-BAM method provides better estimation performance, for example: 1) the computational complexity is reduced due to dealing with low-dimension beamspace signals and avoiding spectral search; 2) it can not only improve the angular parameter estimation accuracy but also have excellent robustness to the change of signal-to-noise ratio (SNR) and snapshot number. The theoretical analysis and simulation results confirm the effectiveness of the proposed method.


Introduction
In millimeter wave (mmWave) massive multiple-input multiple-output (MIMO) systems, accurate estimation of direction of arrival (DOA) can not only improve the energy efficiency, but also enhance the performance of key techniques, such as spatial modulation and non-orthogonal multiple access (NOMA) [1,2]. Conventionally, the noise subspace based 1 3 DOA estimation algorithms [3,4], such as multiple signal classification (MUSIC) and estimation of signal parameter via rotational invariance technique (ESPRIT), have found widespread usage and have become very popular due to their ease and simplicity. Most of DOA estimation algorithms are based on the point source model, while diffusion and multipath propagation cannot be neglected in numerous practical applications. Therefore, the distributed source is more suitable than the point one [5]. In general, the distributed sources can be divided into coherently distributed (CD) sources and incoherently distributed (ID) sources [6]. This paper will focus on the the angular parameter estimation of ID sources.
In recent years, several high-resolution angular parameter estimation methods for ID sources have been proposed. The dispersed signal parameter estimator is one of the efficient subspace-based approaches, which exploits the orthogonality between the quasi-signal subspace and the quasi-noise subspace [7] to obtain the nominal DOA and the angular spread estiamtion. Unfortunately, it requires multidimensional spectrum search to estimate the angular parameters, which may result in high computational costs. Therefore, several computationally simpler algorithms are proposed, which estimate the angular parameters via the closed form expression. For example, a low-complexity algorithm based on a novel propagator is given for the nominal DOA estimation, which avoids spectrum search and the eigendecomposition [8]. A generalized shift invariance property and an extended crosscorrelation matrix are used to derive the closed form expressions of the nominal DOAs in [9] and [10], respectively. Besides, via beamspace transformation, the matrix dimension can be reduced, which is also able to reduce the computational cost [11]. However, for large-scale MIMO systems, high computational complexity is undoubtedly a burden.
In large-scale MIMO systems, two-dimensional (2-D) DOAs are necessary for beamforming, wireless location and interference coordination [12]. Some of the existing parameter estimation approaches for one-dimensional (1-D) ID sources can be extended to twodimensional (2-D) sources [13,14]. A low-rank metric is introduced to estimate the 2-D spatial spectrum of ID sources and estimate the angle parameter [15]. In spite of good performance, the computational costs of the above works are relatively high. To overcome this obstacle, several low-complexity parameter estimation solutions for 2-D sources via the rotational invariance technique are proposed [16][17][18][19]. A new approach based on parallel uniform linear arrays is proposed in [16], the presented algorithm firstly estimates the nominal elevation by the modified TLS-ESPRIT method, and then the nominal azimuth is estimated by one-dimensional searching. However, the angular spreads are not taken into account in the aforementioned methods. An ESPRIT-based approach is presented for the estimation of the nominal DOAs and their angular spreads, which enjoys low computational cost since it avoids spectral search [17]. In consideration of the signal non-circular property, a computational efficient method for DOA estimation is given [18], which can enhance the DOA estimation performance such as higher estimation accuracy, without spectrum searching and pairing, etc. A modified conjugate ESPRIT algorithm is presented, which considers the mixed circular and noncircular signals [19]. However, the computational complexity of the aforementioned algorithms may be high due to high dimensional received signals.
Motivated by this, a low-complexity 2-D DOA estimation method based on the beamspace transformation is presented in this paper. The main contributions of the paper are summarized as follows: (i) Via the first order Taylor expansion, a new generalized array manifold (GAM) model is given, in which the nominal DOA is decoupled from the angular spread. And it lays a foundation for the shift invariant structure (SIS). (ii) By the beamspace transformation matrix derived from the SIS, the observed signals are projected into the beamspace to constructed the beamspace data, which has lower dimension than the observed signals. Compared with related references, the dimension of the received signal is reduced, which greatly decreases the computational complexity. (iii) The closed form solutions of the nominal azimuth and the nominal elevation are obtained through the SIS of the beamspace data.
The rest of this paper is organized as follows. Sect. 2 briefly introduces the data model of ID sources. The proposed method is formulated in Sect. 3. In Sect. 4, simulation results are presented to verify the performance of the proposed method. Finally, the conclusion is provided in Sect. 5.

Data Model
Consider a uniform linear orthogonal array (ULOA), which consists of two uniform linear arrays (ULAs) X and Y as shown in Fig. 1. Both X and Y have M sensors with interelement spacing d. Assume that K(K ≤ M) uncorrelated narrowband ID sources impinge on the ULOA. The array output vectors x(t) and y(t) of arrays X and Y can be expressed as [20] where s k (t) is the complex-valued signal transmitted by the kth source. L k denotes the number of rays inside the kth signal. kl (t) stands for the complex-valued ray gain. kl (t) and kl (t) represent the azimuth angle and the elevation angle of the lth ray from the kth signal, respectively. n x (t) and n y (t) stand for the additive Gaussian white noise with zero-mean and equal power 2 n in arrays X and Y, respectively. (1) ID source k Fig. 1 Array geometry of uniform linear orthogonal array a x ( kl (t), kl (t)) and a y ( kl (t), kl (t)) are the array manifold vectors of the arrays X and Y, which can be expressed as a x ( kl (t), kl (t)) = [1, e j 2 and , respectively. is the wavelength of the signal. The superscript (⋅) T represents the transpose operation.
The azimuth angle and the elevation angle kl (t) and kl (t) can be expressed as , respectively, where k and k stand for the nominal azimuth DOA and the nominal elevation DOA for the kth signal, which are the mean values of the azimuth DOA and elevation DOA of all path signals of the distributed source, respectively. ̃ kl (t) and ̃ kl (t) denote the corresponding angular deviations.

Problem Formulation
In this section, the GAM model is firstly established to decouple the nominal DOAs from their angular spreads, so as to estimate the nominal DOAs and the angular spreads separately. Then, via the SIS, the beamspace transformation matrix is derived to perform dimensionality reduction for the observed vectors. Finally, via the SIS of the beamspace received signals, the closed form solutions are given for the estimation of the nominal azimuth DOA and elevation DOA.

GAM Model and SIS
Using the first order Taylor series expansion around k and k , the array manifold vectors a x ( kl (t), kl (t)) and a y ( kl (t), kl (t)) can be approximated as where a � x ( k , k ) and a � y ( k , k ) stand for the partial derivatives of a x ( k , k ) and a y ( k , k ) with respect to k , respectively. a � x ( k , k ) and a � y ( k , k ) represent the partial derivatives of a x ( k , k ) and a y ( k , k ) with respect to k , respectively.
Substituting (2) into (1), the received signals x(t) and y(t) can be rewritten as . It is obvious that h k,1 (t) , h k,2 (t) and h k,3 (t) in (3) are irrelevant to the nominal DOAs, hence (3) can be reformulated into GAM model as It is noteworthy that the generalized array manifold matrix A x ( , ) and A y ( , ) depend only on the nominal azimuth DOA and the nominal elevation DOA, while h(t) is independent of the nominal DOAs and relates to the angular spread. Accordingly, the decoupling of the nominal DOAs and the angular spreads is realized via GAM model.
As shown in Fig. 1, it is obvious that the received signal in array X has similar property to that in array Y , hence taking the received signal x(t) as an example to illustrate the proposed method.
Similar to ESPRIT-based approaches, the ULA X with M sensors need to be divided into two subarrays X 1 and X 2 with equivalent number of sensors in the proposed method, so as to construct the SIS of the observed vectors. The subarray X 1 consists of the first M − 1 sensors and the subarray X 2 is composed of the last M − 1 sensors.
According to the aforementioned partition of the array X, the generalized array manifold matrices A 1 ( , ) and A 2 ( , ) of the subarrays X 1 and X 2 can be respectively expressed as . Note that the array manifold vectors of X 1 and X 2 , namely, a 1 ( k , k ) and a 2 ( k , k ) have the following relationship where ( k , k ) = exp(j 2 dcos k cos k ).
The partial derivative of (6) around k and k are respectively given by where (⋅) � refers to the first-order derivative. Therefore, referring to (5) - (7), the SIS between the generalized array manifold matrixes A 1 ( , ) and A 2 ( , ) can be given by denotes the shift invariant matrix with

Beamspace Transformation Matrix
Note that the implementation of the SIS-based algorithms usually requires eigendecomposition, which may result in high computational cost due to the high dimensional observed signals. Beamspace transformation is one way of reducing observed signal dimension, thereby, an effective beamspace transformation matrix will be given. Assume that there is a M × P(P < M) beamspace transformation matrix W , which can be used for reducing the dimension of the received signal x(t) from M × N to P × N . N is the number of snapshots. The beamspace received signal z(t) can be written as where B = W H A x ( , ) stands for the beamspace array manifold matrix. (⋅) H stands for the conjugate transpose.
Referring to the theorem in [22], it can be known that the beamspace transformation matrix W is a unitary matrix and satisfies the following constrains where F ∈ ℂ P×P is a nonsingular matrix.
Furthermore, if there is a matrix Q ∈ ℂ P×P with Qw H M = 0 p×1 and QF H w H 1 = 0 p×1 , then where w 1 and w M denotes the first and last rows of the beamspace transformation matrix W , respectively. Substituting (9) and (11) into (12), it follows that Therefore, the following equation can be easily obtained It is obvious that when F and Q are determined, x can be derived to estimate the nominal azimuth DOA and the nominal elevation DOA. As shown in [22], when standard Fourier basis vectors are used to form W = [c 1 , ⋯ , c P ] , namely Accordingly, F and Q can be constructed as When the beamspace transformation matrix W is construct based on (16), the dimension of the received signal can be reduced to P × N . Then, the nominal azimuth DOA and the nominal elevation DOA can be estimated by utilizing the beamspace received signal z(t) and the above matrices F and Q.

Estimation of the Nominal Azimuth and Elevation
Let R denote the covariance matrix of the beamspace received signal z(t) . Assume that all the signals are uncorrelated, the covariance matrix R can be written as where E{⋅} denotes the mathematical expectation. E s and E n are referred to as the signal subspace and noise subspace, which are formed by the eigenvectors of R corresponding to the 3K largest and (M − 3K) smallest eigenvalues, respectively. s and n are the diagonal matrices, which are composed of the 3K largest and (M − 3K) smallest eigenvalues of R , respectively.
After constructing the beamspace transformation matrix W , the closed form expressions of the nominal azimuth DOA and the nominal elevation DOA can be derived via (15). However, the beamspace array manifold matrix B( , ) cannot be obtained directly, hence based on the subspace theory, the signal subspace E s spans the same column space of B( , ) , namely, where T is a nonsingular matrix.
Substituting (19) into (15), (15) can be rewritten as Similarly, assume that the eigenvalue transformation matrix of array Y is y . Then performing eigenvalue decomposition on x and y yields where U x and U y represent the eigenvector matrices of x and y , respectively. 1, 2, ⋯ , 3K) being the eigenvalue of x . y = diag{λ y,1 , λ y,2 , ⋯ , λ y,3K } with y,q (q = 1, 2, ⋯ , 3K) being the eigenvalue of y . Let x,3(k−1)+r (k = 1, 2, ⋯ , K, r = 1, 2, 3) and y = diag{ y,1 , y,2 , ⋯ , y,K } with . It is noteworthy that the eigenvalues in x and y are estimated independently based on the received signals of subarrays X and Y. Therefore, they may be mismatched and should be paired using the matching technique in [18].
Based on the above discussion, Table 1 provides the proposed angular parameter estimation method for ID source. Apparently, the closed form solutions of the nominal azimuth DOA and the nominal elevation DOA are derived in Step.5. (1) Construct the beamspace transformation matrix W via (16) and obtain the beamspace received signal z(t) using (10); (2) Calculate the covariance matrix R of the beamspace received signal z(t) and perform the eigendecomposition of R to obtain the signal subspace E s ; (3) Solve (20) utilizing the TLS criterion to get x . Similarly, y can be denoted in the same way; (4) Perform the eigendecomposition on x and y to obtain x and y , which is matched using the matching technique in [18]; (5) Estimate the nominal azimuth DOA and elevation DOA ̂k and ̂k via (22) and (23).

Computational Complexity
In this section, the complexity of the proposed algorithm is analyzed and compared with three prevalent algorithms. In order to facilitate the description, the algorithms in [17], [21] and [24] are called ESPRIT-2D, BS-2D and URA-2D here, respectively. The computational complexity of the proposed algorithm mainly includes: 1) the computation of the beamspace received signal z(t) , of order O (PMN) ; 2) the eigendecomposition of an P × P covariance matrix, of order O(P 3 ) ; 3) the calculation of the eigenvalue transformation matrixes x and y , of order O(6K 3 ) . 4) the pair matching of eigenvalues, of order O (18MK 4 ) . Consequently, the overall complexity of the proposed estimation algorithm is of order O(PMN). Table 2 compares the computational costs of the proposed algorithm with ESPRIT-2D, BS-2D and URA-2D. Evidently, the complexity of the proposed algorithm is much lower than other algorithms in the table.

Results and Discussion
In this section, several simulation results are presented to evaluate the performance of the proposed SIS-BAM algorithm, in which the ULOA is composed of two ULAs and equipped with 36 isotropic sensors, namely, M = 18 . Meanwhile, this paper compares the performance of the proposed algorithm with three related algorithms, namely, ESPRIT-2D, BS-2D and URA-2D. For ESPRIT-2D, BS-2D and URA-2D, consider a uniform rectangular array with M = 6 , namely, the total number of sensors is 36. Unless stated otherwise, the SNR and the number of snapshots are set to 20 dB and 512. Two uncorrelated ID sources with angular parameters ( 1 , 1 , 1 , 1 ) = (10 • , 1 • , 30 • , 1 • ) and ( 2 , 2 , 2 , 2 ) = (50 • , 1 • , 40 • , 1 • ) are considered. The number of rays inside the ID sources is set to 200. The root-mean-square error (RMSE) is adopted to measure the estimation performance, which is defined as where r and r are the RMSEs of the nominal azimuth and the nominal elevation estimation, respectively. R denotes the number of Monte Carlo simulations, which is set to  R = 200 in this paper. k (r) and k (r) stand for the estimation of k and k for the rth Monte Carlo trial (r = 1, 2, ⋯ , R) , and K is the number of signals.
In the first example, Fig. 2 gives the RMSE curves of the nominal azimuth and the nominal elevation estimates for the proposed method, ESPRIT-2D, BS-2D and URA-2D under the hypothesis that the SNR ranges from 0 dB to 30 dB. From Fig. 2, it clearly indicates that the proposed algorithm shows better the nominal azimuth and elevation estimation than BS-2D and URA-2D. This is because only one sensor is not employed in the subarray of the proposed method, while for BS-2D and URA-2D, there are 2M − 1 sensors are unused. Compared with ESPRIT-2D, the proposed algorithm is able to provide better nominal azimuth estimation performance, but inferior nominal elevation estimation performance in the high SNR region (SNR>10 dB) . Besides, the proposed method is insensitive to the variation of the SNR, which means that the proposed method is robust to the change of SNR.  The impact of the snapshot number on the proposed method is investigated in Fig. 3 when the number of snapshots varies from 128 to 1024. It can be seen that for both nominal azimuth and nominal elevation estimation, the performances of the proposed and related algorithms become better gradually as the number of snapshots increases. In addition, the proposed algorithm can achieve the lowest RMSE of the nominal azimuth, and the RMSE of the nominal elevation is lower than BS-2D and URA-2D. When the number of snapshots is greater than 256, the RMSE of the nominal elevation of the proposed algorithm is higher than that of ESPRIT-2D. The reason is that the proposed algorithm uses the beamspace transformation matrix to perform low-dimensional projection of array received data, which will cause part of the data loss, and the loss rises as the number of snapshots increases. Significantly, the proposed algorithm offers a steady performance over the whole range of the number of snapshots , which shows that the proposed algorithm is robust to the change of the snapshot number.
In the last example, the influence of the angular spread is shown in Fig. 4, which plots the RMSEs of the nominal azimuth and the nominal elevation estimates. It can be seen that the estimation of the nominal azimuth and the nominal elevation by our method is more accurate than that by the other three methods as the angular spread increases. Meanwhile, the RMSEs of the aforementioned algorithms ascend when the angular spread increases. This is very reasonable as the approximation of the array manifold is under the assumption of small angular spread, and the degree of signal energy diffusion rises with the increase of angular spread.

Conclusions
This paper presents an effective SIS-BAM method for 2-D ID source angular parameter estimation. In the proposed method, the first-order Taylor expansion is performed on the array manifold vector to decouple the nominal DOA and the angular spread. Then, by the beamspace transformation matrix, the low-dimension beamspace data is obtained, whose SIS can be established utilizing generalized array manifold. Finally, the closed-form The theory analysis shows that the complexity of the proposed algorithm is lower than that of ESPRIT-2D, BS-2D and URA-2D. Simulation results indicate that the proposed algorithm is higher than other related algorithms in terms of accuracy and is robust to the change of SNR and the number of snapshots.