Two Rapid Power Iterative DOA Estimators for UAV Emitter Using Massive/Ultra-massive Receive Array

To provide rapid direction finding (DF) for unmanned aerial vehicle (UAV) emitter in future wireless networks, a low-complexity direction of arrival (DOA) estimation architecture for massive multiple input multiple output (MIMO) receiver arrays is constructed. In this paper, we propose two strategies to address the extremely high complexity caused by eigenvalue decomposition of the received signal covariance matrix. Firstly, a rapid power-iterative rotational invariance (RPI-RI) method is proposed, which adopts the signal subspace generated by power iteration to gets the final direction estimation through rotational invariance between subarrays. RPI-RI makes a significant complexity reduction at the cost of a substantial performance loss. In order to further reduce the complexity and provide a good directional measurement result, a rapid power-iterative Polynomial rooting (RPI-PR) method is proposed, which utilizes the noise subspace combined with polynomial solution method to get the optimal direction estimation. In addition, the influence of initial vector selection on convergence in the power iteration is analyzed, especially when the initial vector is orthogonal to the incident wave. Simulation results show that the two proposed methods outperform the conventional DOA estimation methods in terms of computational complexity. In particular, the RPIPR method achieves more than two orders of magnitude lower complexity than conventional methods and achieves performance close to CRLB. Moreover, it is verified that the initial vector and the relative error have a significant impact on the performance of the computational complexity.


I. INTRODUCTION
Unmanned aerial vehicle (UAV) plays an important role in wireless communication systems, providing stable and effective wireless connection in areas without extensive communication infrastructure coverage [1]- [4].Due to its economical and flexible arrangement, UAV is widely used in emergency rescue, traffic control, etc [5], [6].For example, UAV can provide emergency communication connection to ground rescue equipment when the mountain fire occurs.Compared with traditional ground-based communications, low-altitude UAVs have shorter line-of-sight (LoS) paths, which can effectively avoid interference and achieve better communication performance.However, its high mobility leads to rapid changes in channel state information (CSI), and ground equipment needs fast and accurate estimation of CSI to achieve high quality communication [7].
Direction of arrival (DOA) is the key information for channel estimation [8], [9], which combined with multiple input multiple output (MIMO) technology [10]- [12] can achieve secure and energy-efficient UAV information transmission by providing highly accurate desired signal direction for directional modulation, beamforming, alignment and tracking [13], [14].In [15], the authors proposed a novel 3-D framework for UAV localization and the key of it is to measure the angle of arrival information.
However, due to the fact that the number of antennas tends to be massive in MIMO system [16], [17], the computational complexity and circuit cost are too high for commercial applications.A DOA-aided channel estimation method was proposed in [18] for a hybrid millimeter-wave MIMO system based on a uniform planar array at the base stations, and the theoretical bounds of the mean squared errors (MSEs) and the Cramer-Rao lower bounds (CRLBs) of the joint DOA and channel gain estimation are derived.The simulation results show that the performances of the proposed methods are close to the theoretical MSEs' analysis, while the theoretical MSE is close to the CRLB.Therefore, hybrid analog and digital (HAD) beamforming structures using parametric method to estimate DOA have emerged, which can achieve a good balance between beamforming computation, circuit cost, and complexity, using parametric method to estimate DOA.
Large antenna arrays using HAD architectures can provide large apertures with low cost and hardware complexity, resulting in enhanced DOA estimation and reduced power consumption.The DOA estimation and power consumption tradeoff problem for large antenna arrays with HAD structures was presented in [19], where the fully-connected, sub-connected (SC), and switches-based (SE) hybrid architectures was formulated into a unified expression, with the compression matrix in a time-varying form.Based on this model, the authors derived dynamic maximum likelihood (D-ML) estimators for HAD and conventional all-digital (FD) structures, and closed expressions for CRB to evaluate the performance limits of D-ML estimators for different HAD structures.
Authors in [20] investigated the DOA estimation using HAD structure in the receiver part and proposed two phase alignment (PA) methods: HAD-PA and hybrid digital and analog PA (HDA-PA).Meanwhile, for this hybrid structure, a fast Root-MUSIC-HDAPA method was proposed to achieve an approximate analytical solution and reduce the computational complexity.In [21], a new design of analog phase shifts was proposed to tackle the phase ambiguities.This enables the cross-correlations to be deterministically calibrated and constructively combined for the noise-tolerant estimation of the propagation phase offset between adjacent subarrays.It is obvious from the simulation that the estimation accuracy of the method can be significantly improved by several orders of magnitude and asymptotically approaches the MSELB.For the DF ambiguity problem caused by HAD MIMO, a fast ambiguity phase elimination method was proposed in [22], which uses only two data blocks to achieve DOA estimation.In [23], the DOA estimation problem in the case of 1bit ADC was considered.It demonstrated that the MUSIC method could be directly applied in the case without additional preprocessing, while the system performance degradation was analyzed.Then, the DOA estimation performance of the lowresolution ADC structure was investigated in [24].
A generalized sparse Bayesian learning (Gr-SBL) method was considered in [25] to solve the DOA estimation problem from one-bit quantized measurements in both single and multi snapshot scenarios.By formulating the one-bit DOA estimation in single-fast-tempo is transformed into a generalized linear model inference problem and solved by applying the recently proposed Gr-SBL method.Then, Gr-SBL is extended to multi-fast-tempo scenarios by decoupling the multi-fasttempo single-bit DOA estimation problem into a series of single-fast-tempo subproblems.Simulation results demonstrate the effectiveness of the Gr-SBL method.A DOA estimation method that is suitable for Non-circular signals with a single snapshot was proposed in [26].By utilizing the waveform characteristics of the NC signal, the proposed algorithm can enlarge the virtual array aperture that is twice the length of the physical array, and resultantly enhance DOA estimation accuracy.Finally, the numerical simulation results are provided to demonstrate the effectiveness and superiority of the proposed method.
In order to avoid the high-complexity operation of eigenvalue decomposition (EVD) in DOA estimation, deep learning network (DNN) has been applied to DOA estimation in recent years.A DNN-based DOA and channel estimation schemes was proposed in [27], which achieved better performance.In [28] the authors introduced a low-complexity DNN-based method to a hybrid massive MIMO system with uniform circular array at the base station, which had similar or even better performance compared to the traditional ML method with lower complexity.An ESPRIT-based HAD method was proposed in [29], which considered a machine learning framework to improve the estimation accuracy.
Aiming at rapidly estimating and tracking the main sub-spaces and major components of a vector sequence in [30], a projection approximation and subspace tracking (PAST) method was proposed.Furthermore, the proposed PAST method in [30] was improved in [31].It proved that the improved PAST method was better in both subspace estimation and computational complexity.In [32], an improved power iteration (PI) method for modal analysis was proposed.
The simulation results showed that the method significantly reduced the number of unnecessary iterations with a faster computational speed.An iterative method was also proposed in [33], and good results were obtained.
Inspired by the idea of radar target detection, we considered a new SVD-based passive target detection model in [34], which achieved better detection performance.While the complexity of massive MIMO based on covariance matrix decomposition method was extremely high.For example, when the number of antennas are closed to 10000, and the computational complexity was tera(T) FLOPs.Therefore, how to significantly reduce the computational complexity of direction finding for UAV emitters is an extremely challenging problem, which is the key to its future applications.Therefore, in this paper, we have proposed a rapid convergent power iteration (RPI) structure to achieve high performance with low complexity.The main contributions in this paper are summarized as follows: 1) To significantly reduce the computational complexity of UAV direction finding, two PI-based DOA estimation methods are respectively proposed, which are called RPI-RI and RPI-PR.Here, the sampling covariance of the received signal vector is first computed.An initial vector subjected to power iteration is determined, which replaces the traditional EVD.Then the final DOA estimation is given by the corresponding signal and noise subspaces.The simulation results show that RPI-RI and RPI-PR can achieve better DOA accuracy and lower computational complexity than conventional algorithms, especially the RPI-PR can dramatically reduce the complexity while maintaining performance close to CRLB. 2) To reduce the number of unnecessary iterations and get a faster computation speed, the optional initial vectors are selected which can converge to the desired results.
In each iteration, it must be ensured that the initial vector is not orthogonal to the incident wave, and it is better to keep them away from orthogonality.Through computational analysis, we selected different initial vector values that satisfy the conditions, and analyzed the performance of iterative convergence.Moreover, the computation result often has a great relationship with the relative error.When a good initial vector and relative error are determined, the results with fast convergence and less iterations can be achieved.
The remainder of this paper is organized as follows.Section II describes the structure of the rapid power-iterative estimator for massive/ultra-massive MIMO receiver.In Section III, two low-complexity estimators are proposed and their performance are also analyzed.The simulation results are presented in Section IV.Finally, we draw conclusions in Section V.
Notations: Throughout the paper, x and X in bold typeface are used to represent vectors and matrices, respectively, while scalars are presented in normal typeface, such as x.Signs (•) H and | • | represent conjugate transpose and modulus, respectively.I N denotes the N × N identity matrix.Furthermore, E[•] represents the expectation operator, and x ∼ CN (m, R) denotes a circularly symmetric complex Gaussian stochastic vector with mean vector m and covariance matrix R.
x represents the estimated value of x.To rapidly achieve the direction finding of the UAV, Fig. 1 sketches a low-complexity massive MIMO RPI receiver structure for UAV direction estimation.In this structure, the PI method is considered for apply in a uniformly-spaced linear array (ULA) with N -antennas, and an appropriate Ndimensional initial array manifold is constructed as the input vector for it, which can effectively reduce the number of iterations.Then, the receive covariance matrix of all antennas is exploited as the object matrix to generate an eigenvector V corresponding to the largest eigenvalue.Through constructing signal or noise subspaces, the optimal DOA estimation can be derived in different methods like polynomial rooting [35] and rotational invariance [36].

II. SYSTEM MODEL
In the proposed framework, it is assumed that the narrow band signals s(t)e j2πfct transmitted by UAV emitter will arrive at the array, where s(t) is the baseband signal and f c is the carrier frequency.Then the antennas will capture the signal with different time delays depended on the DOAs.Therefore, the propagation delay of the nth antenna element is expressed as where τ 0 is the propagation delay from the UAV to a reference point on the receive array, θ 0 is the direction of the UAV relative to the line perpendicular to the array.d n is the distance from the nth array element to the reference point, and c is the speed of light.Without loss of generality, we can assume that τ 0 = 0. Thus, τ n = −(d n /c)sinθ 0 .The receive signal vector at array can expressed as where v(t) ∼ CN (0, σ 2 v I N ) is the additive white Gaussian noise (AWGN) vector, and a(θ 0 ) only composed by the phase difference of all antennas is called array manifold, defined by In practice, although the ideal covariance matrix can not obtained directly, its estimated value is given by where K is the number of sampling points.Make the eigenvalue decomposition of ( 4) The CRLB is the minimum variance of DOA estimation errors.It can provide a useful characterization of the achievable accuracy of the systems.According to [8], for an N -element linear array we have where and SNR is the signal-to-noise ratio of the signal received at each antenna.

III. PROPOSED TWO RAPID POWER-ITERATIVE ESTIMATOR FOR MASSIVE/ULTRA-MASSIVE MIMO RECEIVER
The subspace-based methods are widely used for direction finding of UAV emitters.But their eigenvalue decomposition brings a horrible complexity when the antenna tends to largescale.Therefore, two low-complexity estimators are proposed in this section based on the RPI structure, the selection of initial vector and computational complexity are followed analyzed.

A. Proposed RPI-RI estimator
It is assume that there exist two subarrays with both N − 1 antennas and overlapped with each other [36].The first subarray consists of the first N − 1 antennas of all antennas, and the second subarray consists of the last N − 1 antennas of all antennas.Since the structures of the two subarrays are identical, the outputs of the two subarrays have only one phase difference φ.
The following assumes that the received data of the first subarray is X 1 , and the received data of the second subarray is X 2 , Combine two subarrays, the 2N − 2 × 2N − 2 received data can expressed as where Φ = e jφ is a rotational invariance matrix, which contains the direction of arrival information of the incoming wave signal, and φ = 2πd sin θ/λ.The covariance matrix R can given by Let us define M = 2N − 2 and assume that the covariance matrix R has M eigenvalues λ 1 , λ 2 , ..., λ M with an associated collection of linearly independent eigenvectors {u 1 , u 2 , ..., u M }.Moreover, we assume that R has precisely one eigenvalue λ 1 , which is the largest in magnitude, i.e., There exists a random vector v 0 ∈ R n that satisfies where α 1 , α 2 , ..., α M (α 1 = 0) are scalars.
Taking the vector v 0 as the initial vector of the RPI-RI method, and multiplying both sides of this equation by R, R 2 , ..., R n , ... gives The v n can be rewritten as where Since therefore, The eigenvector corresponding to the main eigenvalue λ 1 is where λ 1 is expressed by the limit of the ratio of the ith component of vector v n+1 to the ith component of vector v n , we have The estimated value of λ 1 is In cases for which the power method generates a good approximation of a dominant eigenvector v λ1 , the Rayleigh quotient [37] provides a correspondingly good approximation of the dominant eigenvalue λ 1 The key problem solved by the ESPRIT algorithm is the proper use of the translation-invariant property of the linear array, so that the eigenvalues of the rotation-invariant matrix can be found to estimate the signal incidence angle.
Based on ( 12) and ( 17), we perform PI method on the covariance matrix R and get the signal subspace u λ1 .Since the shift invariance of the array implies that u λ1 can be decomposed as where the two parts u λ11 and u λ12 corresponding to the signal subspaces of the subarrays X 1 (t) and X 2 (t).According to the ESPRIT algorithm [36], the similar matrice of Φ can be expressed as Therefore, the corresponding eigenvalues of Φ, i.e. the diagonal elements, can be given by performing eigenvalue decomposition on Ψ.The final DOA estimation is followed calculated by where φ is the eigenvalues of Φ.The specific algorithm steps of the proposed RPI-RI estimator is described in Algorithm 1 as follows Algorithm 1 Proposed RPI-RI estimator 1: Input subarray 1 and subarray 2 to receive X 1 and X 2 , forming the receive data model X based on (8); 2: Calculate the covariance matrix R and use it as the object matrix of the power iteration to get the signal subspace corresponding to the two subarrays, and further find the similar matrix Ψ of Φ; 3: The corresponding eigenvalues are given by EVD on Φ and the final DOA estimation θ is calculated based on (28).

B. Proposed RPI-PR estimator
The RPI-RI estimator can significantly reduce the computational complexity by using the power iteration and the rotational invariance of the signal subspace between subarrays, but it is difficult to achieve the desired performance and the complexity can be further optimized.
In order to further reduce the computational complexity and achieve better performance, RPI-PR estimator is proposed in this subsection.From ( 17), the signal subspace can given by power iterating over R y in (4).Furthermore, we construct the noise subspace where I is N × N unit matrix.The spatial spectral function can defined as which spectral peak corresponds to the desire DOA estimation.Furthermore, let us define z = e 2π/λdsinθ , the polynomial equation can be expressed as The above polynomial equation ( 26) has 2N − 2 roots, i.e z i , i = 1, • • • , 2N − 2, which implies the existence of multiple emitter directions as follows where The angle corresponding to the root inside the unit circle and closest to it is chosen as the final DOA estimation θ.
And the basic steps of the proposed RPI-PR method can be summarized in Algorithm 2 as

C. Selecting initial vector and relative error
The initial vector v 0 has a direct effect on the speed of convergence and determines the number of iterations.A randomly generated vector is usually selected as the iterative initial vector, but not all the optional initial vectors can converge to get the desired results.Since orthogonality may cause the iterative process fails to converge.Therefore, it is Algorithm 2 Proposed RPI-PR method 1: Input array receive signal y(t); 2: Based on (4), calculate the covariance matrix R y ; 3: Based on (24),to estimate the noise subspace U V and use it to construct the spatial spectrum S(θ) 4: Multiple roots are given by polynomial method for the spatial spectral function, the angle corresponding to the root inside the unit circle and closest to it is chosen as the final DOA estimation θ. necessary to ensure that the initial vector is not orthogonal to the incident wave in each iteration.
From ( 11) and ( 13), the v 0 is formed as where S(R) is the sum of all elements in matrix R, and Thus, we can calculate Based on the above discussion, the element distribution of initial vector v 0 in ( 29) is consistent with the distribution of matrix R, while the distribution is uniform, which can keep them far away from orthogonality and speed up the convergence of iteration.
Therefore, let us define the array manifold as where and the initial vector v 0 is assumed as The orthogonality equation can be expressed as In order to make the incident wave direction not orthogonal to the initial vector, i.e., the a(θ 0 ) T v 0 = 0 is constant.In the following, eight special initial vectors are discussed and their convergence performance is analyzed.I.The initial vector is assumed as where all elements are 1.Based on (36), corresponding orthogonality equation is given by where φ = 2π d λ sinθ 0 and it can be discussed in three cases.(1) When 1 − e jφ = 0 and 1 − (e −jφ ) N = 0, φ can be calculated as where Z is integer.It can be seen from ( 39) that a(θ 0 ) T v 0 = 0 holds when θ 0 = arcsin( λZ N d ) and Z N is not integer.(2) When 1 − e jφ = 0, φ is given by Substituting φ = 2πZ into (38), a(θ 0 ) T v 0 = 0 is verified to hold.
(3) When 1 − e jφ = 0 and 1 − (e −jφ ) N = 0, it is clear that a(θ 0 ) T v 0 = 0 holds.Through the above discussion, initial vector I is not guaranteed non-orthogonal always true.II.The initial vector is assumed as where the odd elements is 1 and the even elements is -1.
Depending on the value of N , two cases can be discussed.
(2) N is odd It can be seen that a(θ 0 ) T v 0 = 0 holds when φ satisfies the following conditions N is not integer and orthogonality only can be avoided by selecting other φ.
III.The initial vector is assumed as Assuming N is even.The first half of the elements of vector v 0 are 1, and the rest are -1.
where m-th element is 1, the n-th is 1 and the rest of the elements are 0 (m < n).
When φ = (2Z+1)π n−m , a(θ 0 ) T v 0 = 0. VI.The initial vector is assumed as where m, n, k-th is 1 and the rest of the elements are 0 (m < n < k).
The initial vector is assumed as where m, n, k, l-th is 1 and the rest of the elements are 0 (m < n < k < l).VIII.The initial vector is assumed as where k-th element of e k is 1, and the rest of the elements are 0. a(θ 0 ) T v 0 = 0 is always true.Among the above eight vector forms, the eighth vector form can be used as the initial vector of the power iteration due to its non-orthogonality, while the previous seven cannot guarantee the convergence of the power iteration.
In addition, the relative error ε also influence the speed of convergence in the convergence process of power iteration, which can be expressed as and Therefore, we have When the iteration approaches the final convergence, the convergence speed will be relatively slow and stable.In practice, the usefulness of the power method depends upon the ration |λ2| |λ1| , since it dictates rate of convergence.The whole procedure is summarized in Algorithm 3. where Algorithm 3 PI method on subspace Input: R and ε Output: λ 1 , v n and n Initialization: choose a initial vector v 0 , and n=1.repeat n is the number of iterations, λ 1 is dominant eigenvalue, and v n is the dominant eigenvector corresponding to the λ 1 .
Notice that in each iteration we compute a single matrixvector multiplication (O(N 2 )).We never perform matrixmatrix multiplication, which requires greater number of operations (O(N 3 )).If the matrix R is sparse (only a small portion of the entries of A are non-zero), matrix-vector multiplication can be performed very efficiently.Therefore, the power method is practical even if N is very large, such as in Google's Page Rank algorithm.

D. Complexity analysis
we make an analysis of computational complexities of the proposed two estimators with traditional ESPRIT and Root-MUSIC algorithms as a complexity benchmark.Thus, the complexity of RPI-RI is as follows FLOPs.The complexity of RPI-PR is FLOPs, where β is iteration number of the power iteration.
Assuming N is far larger than K, β, compared with the conventional FD estimator, the complexity of the proposed two estimators is significantly reduced as the number of antennas tends to large-scale.

IV. SIMULATION RESULTS AND DISCUSSIONS
In this section, we provide the simulation to analyze the performance of the proposed massive MIMO DOA estimators for UAV and two conventional algorithms are used as a comparison.Furthermore, we consider the effect of SNR, the number of antennas and the number of snapshots on the proposed methods in the digital ADC architecture.In each simulation figure, the number of estimation value L is 2000.Without loss of generality, we assume that UAV emitter located in θ 0 = 50 o , the number of snapshots K is 1000, antennas distance d = λ/2, and the number of antenna elements N is 64 in the massive MIMO system.performance benchmark.Observing Fig. 2, it is clear that the proposed RPI-RI method can achieve a similar performance to ESPRIT but with some performance loss, while the proposed RPI-PR method can achieve the corresponding CRLB.Fig. 3 and Fig. 4 present the RMSE of four different methods versus the number of antennas N and snapshots K. Without loss of generality, we assume that N ∈ [16,272] and K ∈ [100, 3600].From the Fig. 3 and Fig. 4, it can be seen that the DOA estimators on the power iteration method achieves a performance close to the conventional algorithms for any number of antennas and snapshots scenarios, especially RPI-PR can reach the FD CRLB.Fig. 5 plots curves of complexity analysis versus the number of antennas N with K = 50, SNR=0dB.From this figure, it can be seen that the complexity of all methods gradually increases as the total number of antennas increases.However, the computational complexity of our proposed methods is two to three orders of magnitude lower when N = 1024 compared to the conventional methods.In particular, the Proposed RPI-PR method can achieve a performance close to CRLB.

RM SE
To explore the influence of the selection of the initial vector on the number of iterations and the convergence speed, Fig. 6 shows the relationship between the optimal eigenvector value ν n and the number of iterations n, given three different initial vectors ν 0 .When the initial vector ν 0 is infinitely close to the signal subspace, only two iterations are needed to complete the convergence, and the result is similar to that of the initial vector selected in (29).In addition, when the ν 0 obeys a random vector distribution, the required number of iterations n is more, which requires an average of eight iterations to reach the convergence of the max(ν n ).
In order to verify the convergence performance of the initial vector v 0 = e k = [• • • , 0, 0, • • • , 1, • • • ] T in (59), five incident wave directions θ 0 = {30 o , 60 o , 90 o , 120 o , 150 o } are assumed in Fig. 7 and the number of iterations of RPI methods is given by numerical simulation when SNR=0.As shown in Fig. 7, To analyze the effect of SNR on the convergence speed of the proposed RPI methods.Fig. 8 plots the iteration error α n of the RPI methods versus the number of iterations, where the iteration error α n represents the absolute value of the difference between the result of the nth iteration C(θ n ) and the result of the (n − 1)th iteration C(θ n−1 ), i.e.α n = |C(θ n )−C(θ n−1 )|.As shown in Fig. 8, assuming the number of antennas N = 1024 and the incident wave direction θ 0 = 50 o .The number of iterations is 7, 5 and 4 for SNRs of -30dB, 0dB and 30dB, respectively.the number of iterations of the proposed algorithm decreases sequentially as the SNR increases.

V. CONCLUSIONS
To find the direction of UAV emitter rapidly and accurately, two low-complexity DOA estimators based on large-scale MIMO arrays are proposed.By determining good initial vector and relative errors, the proposed two algorithms can achieve fast convergence and lower complexity than conventional algorithms.In particular, the PI-PR method achieves more than two orders of magnitude complexity reduction and maintains performance close to CRLB.Adopting the proposed algorithms makes fast direction finding of UAV based on massive MIMO receiver feasible for future practical applications.

Fig. 2 Fig. 4 .
Fig.2plots the root mean square error (RMSE) curves using four different methods versus SNR with CRLB as the

Fig. 5 .
Fig. 5. Computational complexity over the number of antennas N with SNR=0dB and K = 1000

Fig. 8 .
Fig. 8. Curves of the number of iterations with different SNR