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

Abstract

To provide rapid direction finding (DF) for unmanned aerial vehicle (UAV) emitters 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 obtain the final direction estimation through rotational invariance between subarrays. RPI-RI causes a significant complexity reduction at the cost of a substantial performance loss. In order to further reduce the complexity and provide good directional measurement results, the rapid power iterative polynomial rooting (RPI-PR) method is proposed, which utilizes the noise subspace combined with the polynomial solution method to obtain 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 RPI-PR method achieves more than two orders of magnitude lower complexity than conventional methods and achieves performance close to the Cramér–Rao Lower Bound (CRLB). Moreover, it is verified that the initial vector and the relative error have a significant impact on the performance with respect to the computational complexity.

1. Introduction

Unmanned aerial vehicles (UAVs) play an important role in wireless communication systems, providing stable and effective wireless connection in areas without extensive communication infrastructure coverage [1,2,3,4]. Due to their economical and flexible arrangement, UAVs are widely used in emergency rescue, traffic control, etc. [5,6]. For example, UAVs can provide emergency communication connection to ground rescue equipment when mountain fires occur. 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, their high mobility leads to rapid changes in channel state information (CSI) and ground equipment requires fast and accurate estimation of CSI to achieve high quality communication [7].

Direction of arrival (DOA) is key information for channel estimation [8,9], which combined with multiple-input multiple-output (MIMO) technology [10,11,12] can achieve secure and energy-efficient UAV information transmission by providing a highly accurate desired signal direction for directional modulation, beamforming, alignment and tracking [13,14]. In [15], the authors proposed a novel 3D framework for UAV localization; the key is to measure the angle of arrival information.

However, due to the fact that the number of antennas tends to be massive in MIMO systems [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 Cramér–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 MSE analysis, while the theoretical MSE is close to the CRLB. Therefore, hybrid analog and digital (HAD) beamforming structures using a parametric method to estimate DOA have emerged, which can achieve a good balance between beamforming computation, circuit cost and complexity, using a parametric method to estimate DOA.

Large antenna arrays using HAD architectures can provide large apertures at low cost and with 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 and switch-based hybrid architectures were 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 fully digital (FD) structures and closed expressions for CRB to evaluate the performance limits of D-ML estimators for different HAD structures.

The authors of [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 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 1-bit ADC was considered. It demonstrated that the MUSIC method could be directly applied in cases without additional preprocessing, while system performance degradation was analyzed. The DOA estimation performance of the low-resolution 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, it 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-fast-tempo 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 as a result 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 subspaces 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. The complexity of massive MIMO based on the covariance matrix decomposition method was extremely high, for example, when the number of antennas is closed to 10,000 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 propose a rapid convergent power iteration (RPI) structure to achieve high performance with low complexity. The main contributions in this paper are summarized as follows:

• 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; in particular, the RPI-PR can dramatically reduce the complexity while maintaining performance close to CRLB.
• To reduce the number of unnecessary iterations and obtain 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, results with fast convergence and fewer iterations can be achieved.

The remainder of this paper is organized as follows. Section 2 provides an in-depth discussion of related work. Section 3 describes the structure of the rapid power iterative estimator for massive/ultra-massive MIMO receiver. In Section 4, two low-complexity estimators are proposed and their performance is also analyzed. The simulation results are presented in Section 5. Finally, we draw conclusions in Section 6.

Notations: Throughout the paper, $\mathbf{x}$ and $\mathbf{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. ${\mathbf{I}}_N$ denotes the $N\times N$ identity matrix. Furthermore, $\mathbb{E}[·]$ represents the expectation operator and $\mathbf{x}\sim \mathcal{CN}(\mathbf{m},\mathbf{R})$ denotes a circularly symmetric complex Gaussian stochastic vector with mean vector $\mathbf{m}$ and covariance matrix $\mathbf{R}$. $\widehat{x}$ represents the estimated value of x.

Definitions of terms: DOA (direction of arrival) is the process of finding the angle or direction from which a signal is arriving at a receiver, with applications in various fields such as radar, sonar, wireless communication and audio signal processing.

HAD (hybrid analog and digital) is a special array structure. Analog components are used to process continuous signals, while digital components are used to process discrete signals. In the receiver array, it can be simply divided into fully connected, sub-connected and switch-based.

D-ML (dynamic maximum likelihood) is a generalized maximum likelihood algorithm for the three HAD structures as well as for the fully digital structure.

Root-MUSIC-HDAPA is a low-complexity DOA estimator. It firstly applies the conventional Root-MUSIC algorithm in the digital domain, followed by digital phase alignment. Finally, analog phase alignment is performed in the analog domain to eliminate pseudo-solutions.

PAST (projection approximation subspace tracking) is a robust and low-complexity signal subspace tracking method that solves the problem by making an appropriate projection approximation to the recursive least squares technique.

2. Related Work

There has been extensive related research conducted on adopting MIMO technology and UAVs to enhance wireless sensor network connectivity and improve user data acquisition. In [35], the authors utilized UAVs as relay networks and employed cooperative MIMO techniques to facilitate information transmission in disconnected regions between sparse WSNs, achieving low latency and high stability of communication. In [36], the authors propose a novel data acquisition scheme by equipping the receiver node with MIMO technology and making it well planned for mobility. The experimental results demonstrate that this approach can significantly enhance system throughput and energy efficiency. To determine the number of radiation nodes, the authors propose a UAV-based comprehensive DOA preprocessing system in [37]. This system consists of two signal detectors and a verifier that can make precise inferences based on the MIMO receiver array, which provides the foundation for the following DOA estimation.

3. System Model

To rapidly achieve the direction finding of the UAV, Figure 1 sketches a low-complexity massive MIMO RPI receiver structure for UAV direction estimation. In this structure, the PI method is considered for application in a uniformly spaced linear array (ULA) with N antennas and an appropriate N-dimensional 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 such as polynomial rooting [38] and rotational invariance [39].

In the proposed framework, it is assumed that the narrow band signals $s\left(t\right)e^{j2\pi f_{c}t}$ transmitted via the UAV emitter will arrive at the array, where $s\left(t\right)$ is the baseband signal and $f_c$ is the carrier frequency. Then the antennas will capture the signal with different time delays dependent on the DOAs. Therefore, the propagation delay of the nth antenna element is expressed as

$$\begin{array}{c}\hfill {\tau }_n={\tau }_0-(d_n/c)sin{\theta }_0,n=1,2,...,N,\end{array}$$

(1)

where ${\tau }_0$ is the propagation delay from the UAV to a reference point on the receive array, ${\theta }_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 ${\tau }_0=0$. Thus, ${\tau }_n=-(d_n/c)sin{\theta }_0.$ The receive signal vector at the array can be expressed as

$$\begin{array}{c}\hfill \mathbf{y}\left(t\right)=\mathbf{a}\left({\theta }_0\right)s\left(t\right)+\mathbf{v}\left(t\right),\end{array}$$

(2)

where $\mathbf{v}\left(t\right)\sim \mathcal{CN}(0,{\sigma }_v^2{\mathbf{I}}_N)$ is the additive white Gaussian noise (AWGN) vector and $\mathbf{a}\left({\theta }_0\right)$ only composed by the phase difference of all antennas is called array manifold, defined as

$$\begin{array}{c}\hfill \mathbf{a}\left({\theta }_0\right)={[e^{j2\pi d_{1}sin{\theta }_0/\lambda },e^{j2\pi d_{2}sin{\theta }_0/\lambda },\cdots ,e^{j2\pi d_{N}sin{\theta }_0/\lambda }]}^T,\end{array}$$

(3)

In practice, although the ideal covariance matrix cannot be obtained directly, its estimated value is given by

$${\mathbf{R}}_y=\frac{1}{K}\sum _{n=1}^K\mathbf{y}\left(n\right){\mathbf{y}}^H\left(n\right),$$

(4)

where K is the number of sampling points.

Make the eigenvalue decomposition of (4)

$$\begin{array}{c}\hfill {\mathbf{R}}_y=\mathbf{U}\mathsf{\Sigma }{\mathbf{U}}^H=\left[{\mathbf{U}}_S{\mathbf{U}}_N\right]\mathsf{\Sigma }{\left[{\mathbf{U}}_S{\mathbf{U}}_N\right]}^H.\end{array}$$

(5)

where $\mathsf{\Sigma }$ is the diagonal matrix consisting of all non-zero eigenvalues of ${\mathbf{R}}_y$ and the matrices ${\mathbf{U}}_S$ and ${\mathbf{U}}_N$ stand for the signal and noise subspaces, respectively.

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 the CRB can be given by

$$\begin{array}{c}\hfill CRB=\frac{{\lambda }^2}{8{\pi }^{2}KSNR{cos}^2\theta \overline{d^2}}\end{array}$$

(6)

where

$$\begin{array}{c}\hfill \overline{d^2}=\sum _{m=1}^N{d}_m^2\end{array}$$

(7)

and SNR is the signal-to-noise ratio of the signal received at each antenna.

4. Proposed Two Rapid Power Iterative Estimators for Massive/Ultra-Massive MIMO Receiver

The subspace-based methods are widely used for direction finding of UAV emitters. However, their eigenvalue decomposition brings horrible complexity when the antenna tends to large-scale. Therefore, two low-complexity estimators are proposed in this section based on the RPI structure, the selection of the initial vector and the computational complexity.

4.1. Proposed RPI-RI Estimator

It is assumed that there exist two subarrays both with $N-1$ antennas and overlapping with each other [39]. 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 $\varphi$.

The following assumes that the received data of the first subarray is ${\mathbf{x}}_1$ and the received data of the second subarray is ${\mathbf{x}}_2$, and by combining two subarrays, the $2N-2\times 1$ received data can be expressed as

$$\overline{\mathbf{x}}\left(t\right)=\left[\begin{array}{c}{\mathbf{x}}_1\left(t\right)\\ {\mathbf{x}}_2\left(t\right)\end{array}\right]=\left[\begin{array}{c}{\mathbf{a}}_1\left({\theta }_0\right)\\ {\mathbf{a}}_1\left({\theta }_0\right)\mathsf{\Phi }\end{array}\right]s\left(t\right)+\tilde{\mathbf{v}}={\overline{\mathbf{a}}}_1\left({\theta }_0\right)+\tilde{\mathbf{v}},$$

(8)

where $\tilde{\mathbf{v}}\left(t\right)\sim \mathcal{CN}(0,{\sigma }_{\tilde{v}}^2{\mathbf{I}}_N)$ is the additive white Gaussian noise (AWGN) vector and ${\mathbf{a}}_1\left({\theta }_0\right)$ is an $N-1$ dimensional array manifold vector, which can be defined as

$$\begin{array}{c}\hfill {\mathbf{a}}_1\left({\theta }_0\right)={\left[1,e^{j\frac{2\pi }{\lambda }dsin{\theta }_0},\cdots ,e^{j\frac{2\pi }{\lambda }(N-2)dsin{\theta }_0}\right]}^T,\end{array}$$

(9)

and $\mathsf{\Phi }=e^{j\varphi }$ is a rotation-invariant value, which contains the direction of arrival information of the incoming wave signal and $\varphi =2\pi dsin\theta /\lambda$.

The covariance matrix $\mathbf{R}$ can be given by

$$\mathbf{R}=\frac{1}{K}\sum _{n=1}^K\overline{\mathbf{x}}\left(n\right){\overline{\mathbf{x}}}^H\left(n\right),$$

(10)

Let us define $M=2N-2$ and assume that the covariance matrix $\mathbf{R}$ has M eigenvalues ${\lambda }_1,{\lambda }_2,...,{\lambda }_M$ with an associated collection of linearly independent eigenvectors $\{{\mathbf{u}}_{\mathbf{1}},{\mathbf{u}}_{\mathbf{2}},...,{\mathbf{u}}_{\mathbf{M}}\}$. Moreover, it is assumed that $\mathbf{R}$ has precisely one eigenvalue ${\lambda }_1$, which is the largest in magnitude, i.e.,

$$\begin{array}{c}\hfill |{\lambda }_1|>|{\lambda }_2|>|{\lambda }_3|>···>|{\lambda }_M|\ge 0,\end{array}$$

(11)

There exists a random vector ${\mathbf{v}}^0\in {\mathbb{R}}^n$ that satisfies

$$\begin{array}{c}\hfill {\mathbf{v}}^0=\sum _{k=1}^M{\alpha }_k{\mathbf{u}}_{\mathbf{k}},\end{array}$$

(12)

where ${\alpha }_1,{\alpha }_2,...,{\alpha }_M$$({\alpha }_1\ne 0)$ are scalars.

Taking the vector ${\mathbf{v}}^0$ as the initial vector of the RPI-RI method and multiplying both sides of this equation by $\mathbf{R}$, ${\mathbf{R}}^2$, ..., ${\mathbf{R}}^n,...$ gives

$$\left\{\begin{array}{c}{\mathbf{v}}^1=\mathbf{R}{\mathbf{v}}^0=\mathbf{R}\sum _{k=1}^M{\alpha }_k{\mathbf{u}}_{\mathbf{k}}=\sum _{k=1}^M{\alpha }_k{\mathbf{R}}_k{\mathbf{u}}_{\mathbf{k}}=\sum _{k=1}^M{\alpha }_k{\lambda }_k{\mathbf{u}}_{\mathbf{k}},\hfill \\ {\mathbf{v}}^2=\mathbf{R}{\mathbf{v}}^1={\mathbf{R}}^2{\mathbf{v}}^0={\mathbf{R}}^2\sum _{k=1}^M{\alpha }_k{\mathbf{u}}_{\mathbf{k}}=\sum _{k=1}^M{\alpha }_k{\lambda }_k^2{\mathbf{u}}_{\mathbf{k}},\hfill \\ ···\hfill \\ {\mathbf{v}}^n=\mathbf{R}{\mathbf{v}}^{n-1}={\mathbf{R}}^n{\mathbf{v}}^0={\mathbf{R}}^n\sum _{k=1}^M{\alpha }_k{\mathbf{u}}_{\mathbf{k}}=\sum _{k=1}^M{\alpha }_k{\lambda }_k^n{\mathbf{u}}_{\mathbf{k}}.\hfill \end{array}\right.$$

(13)

${\mathbf{v}}^n$ can be rewritten as

$$\begin{array}{cc}\hfill {\mathbf{v}}^n& ={\lambda }_1^n({\alpha }_1{\mathbf{u}}_{\mathbf{1}}+\sum _{k=2}^M{\alpha }_k{\left(\frac{{\lambda }_k}{{\lambda }_1}\right)}^n{\mathbf{u}}_{\mathbf{k}})\hfill \\ \hfill & ={\lambda }_1^n({\alpha }_1{\mathbf{u}}_{\mathbf{1}}+{\epsilon }_n),\hfill \end{array}$$

(14)

where

$$\begin{array}{c}\hfill {\epsilon }_n=\sum _{k=2}^M{\alpha }_k{\left(\frac{{\lambda }_k}{{\lambda }_1}\right)}^n{\mathbf{u}}_{\mathbf{k}},\end{array}$$

(15)

Since $|{\lambda }_1|>|{\lambda }_k|$ for all $k=2,3,...,M$, the above function can be further expressed as

$$\begin{array}{c}\hfill \underset{n\to \infty }{lim}{\epsilon }_n=\underset{n\to \infty }{lim}\sum _{k=2}^M{\alpha }_k{\left(\frac{{\lambda }_k}{{\lambda }_1}\right)}^n=0;\end{array}$$

(16)

therefore,

$$\begin{array}{c}\hfill \underset{n\to \infty }{lim}{\mathbf{v}}^n=\underset{n\to \infty }{lim}{\lambda }_1^n{\alpha }_1{\mathbf{u}}_{\mathbf{1}},\end{array}$$

(17)

The eigenvector corresponding to the main eigenvalue ${\lambda }_1$ is

$$\begin{array}{c}\hfill {\mathbf{u}}_{{\lambda }_1}=\underset{n\to \infty }{lim}\frac{{\mathbf{v}}^n}{{\lambda }_1^n}=\underset{n\to \infty }{lim}\frac{{\lambda }_1^n({\alpha }_1{\mathbf{u}}_{\mathbf{1}}+{\epsilon }_n)}{{\lambda }_1^n}={\alpha }_1{\mathbf{u}}_{\mathbf{1}},\end{array}$$

(18)

where ${\lambda }_1$ is expressed by the limit of the ratio of the ith component of vector ${\mathbf{v}}^{n+1}$ to the ith component of vector ${\mathbf{v}}^n$, which has the following form

$$\begin{array}{cc}\hfill & \underset{n\to \infty }{lim}\frac{{\left({\mathbf{v}}^{n+1}\right)}_i}{{\left({\mathbf{v}}^n\right)}_i}=\underset{n\to \infty }{lim}\frac{{\lambda }_1^{n+1}{({\alpha }_1{\mathbf{u}}_{\mathbf{1}}+{\epsilon }_{n+1})}_i}{{\lambda }_1^n{({\alpha }_1{\mathbf{u}}_{\mathbf{1}}+{\epsilon }_n)}_i}\hfill \\ \hfill & =\underset{n\to \infty }{lim}\frac{{\lambda }_1[{\alpha }_1{\left({\mathbf{u}}_{\mathbf{1}}\right)}_i+{\left({\epsilon }_{n+1}\right)}_i]}{{\alpha }_1{\left({\mathbf{u}}_{\mathbf{1}}\right)}_i+{\left({\epsilon }_n\right)}_i}={\lambda }_1,\hfill \end{array}$$

(19)

The estimated value of ${\lambda }_1$ is

$$\begin{array}{c}\hfill {\widehat{\lambda }}_1=\frac{{\left({\mathbf{v}}^{n+1}\right)}_i}{{\left({\mathbf{v}}^n\right)}_i}.\end{array}$$

(20)

In cases for which the power method generates a good approximation of a dominant eigenvector ${\mathbf{v}}_{{\lambda }_1}$, the Rayleigh quotient [40] provides a correspondingly good approximation of the dominant eigenvalue ${\lambda }_1$

$$\begin{array}{c}\hfill {\lambda }_1=R(\mathbf{R},{\mathbf{v}}_{{\lambda }_1})=\frac{{{\mathbf{v}}_{{\lambda }_1}}^T\mathbf{R}{\mathbf{v}}_{{\lambda }_1}}{{{\mathbf{v}}_{{\lambda }_1}}^T{\mathbf{v}}_{{\lambda }_1}}\end{array}$$

(21)

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 (13) and (18), we perform the PI method on the covariance matrix $\mathbf{R}$ and obtain the signal subspace ${\mathbf{u}}_{{\lambda }_1}$. Since the shift invariance of the array implies that ${\mathbf{u}}_{{\lambda }_1}$ can be decomposed as

$$\begin{array}{c}\hfill {\mathbf{u}}_{{\lambda }_1}=\left[\begin{array}{c}{\mathbf{u}}_{{\lambda }_{11}}\hfill \\ {\mathbf{u}}_{{\lambda }_{12}}\hfill \end{array}\right]\end{array}$$

(22)

where the two parts ${\mathbf{u}}_{{\lambda }_{11}}$ and ${\mathbf{u}}_{{\lambda }_{12}}$ correspond to the signal subspaces of the subarrays ${\mathbf{X}}_1\left(t\right)$ and ${\mathbf{X}}_2\left(t\right)$. In accordance with the ESPRIT algorithm [39], a matrix similar to $\mathsf{\Phi }$ can be expressed as

$$\begin{array}{c}\hfill \mathrm{\Psi }={\left({\mathbf{u}}_{{\lambda }_{11}}\right)}^{†}{\mathbf{u}}_{{\lambda }_{12}}={\left({\mathbf{u}}_{{\lambda }_{11}}^H{\mathbf{u}}_{{\lambda }_{11}}\right)}^{-1}{\mathbf{u}}_{{\lambda }_{11}}^H{\mathbf{u}}_{{\lambda }_{12}},\end{array}$$

(23)

Therefore, the corresponding eigenvalues of $\mathsf{\Phi }$, i.e., the diagonal elements, can be given by performing eigenvalue decomposition on $\mathrm{\Psi }$ and the final DOA estimation can be calculated by

$$\begin{array}{c}\hfill \widehat{\theta }=arcsin\left(\frac{\varphi \lambda }{2\pi d}\right),\end{array}$$

(24)

where $\varphi$ is the eigenvalues of $\mathsf{\Phi }$. The specific algorithm steps of the proposed RPI-RI estimator are described in Algorithm 1 as follows

Algorithm 1 Proposed RPI-RI estimator

1: Input subarray 1 and subarray 2 to receive ${\mathbf{X}}_1$ and ${\mathbf{X}}_2$, forming the receive data model $\overline{\mathbf{X}}$ based on (8);

2: Calculate the covariance matrix $\mathbf{R}$ and use it as the object matrix of the power iteration to obtain the signal subspace corresponding to the two subarrays and further find the matrix $\mathrm{\Psi }$ similar to $\mathsf{\Phi }$;

3: The corresponding eigenvalues are given by EVD on $\mathsf{\Phi }$ and the final DOA estimation $\widehat{\theta }$ is calculated based on (24).

4.2. 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, the RPI-PR estimator is proposed in this subsection. From (18), the signal subspace can be given by power iterating over ${\mathbf{R}}_y$ in (4). Furthermore, we construct the noise subspace

$$\begin{array}{c}\hfill {\mathbf{U}}_V=\mathbf{I}-{\mathbf{u}}_{{\lambda }_1}{\left({\mathbf{u}}_{{\lambda }_1}^H{\mathbf{u}}_{{\lambda }_1}\right)}^{-1}{\mathbf{u}}_{{\lambda }_1}^H,\end{array}$$

(25)

where $\mathbf{I}$ is $N\times N$ unit matrix. The spatial spectral function can be defined as

$$\begin{array}{c}\hfill S\left(\theta \right)=\frac{1}{\parallel {\mathbf{a}}^H\left(\theta \right){\mathbf{U}}_V{\parallel }^2}\end{array}$$

(26)

The spectral peak corresponds to the desired DOA estimation. Furthermore, let us define $z=e^{2\pi /\lambda dsin\theta }$; the polynomial equation can be expressed as

$$\begin{array}{cc}\hfill f\left(\theta \right)& =S^{-1}\left(\theta \right)={\mathbf{a}}^H\left(\theta \right){\mathbf{U}}_V{\mathbf{U}}_V^H\mathbf{a}\left(\theta \right)\hfill \\ \hfill & ={\mathbf{a}}^T\left(\frac{1}{z}\right){\mathbf{U}}_V{\mathbf{U}}_V^H\mathbf{a}\left(z\right)\triangleq f\left(z\right)\hfill \end{array}$$

(27)

The above polynomial Equation (27) has $2N-2$ roots, i.e, $z_i,i=1,\cdots ,2N-2$, which implies the existence of multiple emitter directions as follows

$$\begin{array}{c}\hfill {\widehat{\mathsf{\Theta }}}_{RM}=\left\{{\widehat{\theta }}_i,i\in \{1,2,\cdots ,2N-2\}\right\},\end{array}$$

(28)

where

$$\begin{array}{c}\hfill {\widehat{\theta }}_i=arcsin\left(\frac{\lambda arg{z}_i}{2\pi d}\right)\end{array}$$

(29)

The angle corresponding to the root inside the unit circle and closest to it is chosen as the final DOA estimation $\widehat{\theta }$.

The basic steps of the proposed RPI-PR method can be summarized as in Algorithm 2.

Algorithm 2 Proposed RPI-PR method

1: Input array receive signal y(t);

2: Based on (4), calculate the covariance matrix ${\mathbf{R}}_y$;

3: Based on (25), estimate the noise subspace ${\mathbf{U}}_V$ and use it to construct the spatial spectrum $S\left(\theta \right)$

4: Multiple roots are given by the 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 $\widehat{\theta }$.

4.3. Selecting Initial Vector and Relative Error

The initial vector ${\mathbf{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 obtain the desired results. Since orthogonality may cause the iterative process failures to converge, it is necessary to ensure that the initial vector is not orthogonal to the incident wave in each iteration.

From (12) and (14), ${\mathbf{v}}^0$ is formed as

$$\begin{array}{c}\hfill {\left({\mathbf{v}}^0\right)}^H=\frac{\sum _{j=1}^N{\mathbf{R}}_{.j}}{\mathbf{S}\left(\mathbf{R}\right)}\end{array}$$

(30)

where $\mathbf{S}\left(\mathbf{R}\right)$ is the sum of all elements in matrix $\mathbf{R}$,

$$\begin{array}{c}\hfill \mathbf{S}\left(\mathbf{R}\right)=\sum _{i=1}^N\sum _{j=1}^N{\mathbf{R}}_{ij}\end{array}$$

(31)

and

$$\begin{array}{c}\hfill \sum _{j=1}^N{\mathbf{R}}_{.j}={\mathbf{1}}_{N\times 1}^H\mathbf{R}\end{array}$$

(32)

Thus, we can calculate

$$\begin{array}{cc}\hfill \mathbf{R}{\mathbf{v}}^0& =\mathbf{R}\frac{{\left(\sum _{j=1}^N{\mathbf{R}}_{.j}\right)}^H}{\mathbf{S}\left(\mathbf{R}\right)}=\mathbf{R}{\left({\mathbf{1}}_{N\times 1}^H\mathbf{R}\right)}^H/\mathbf{S}\left(\mathbf{R}\right)\hfill \\ \hfill & =\mathbf{R}{\mathbf{R}}^H{\mathbf{1}}_{N\times 1}={\mathbf{1}}_{N\times 1}\hfill \end{array}$$

(33)

Based on the above discussion, the element distribution of initial vector ${\mathbf{v}}^0$ in (30) is consistent with the distribution of matrix $\mathbf{R}$, while the distribution is uniform, which can keep them far away from orthogonality and speed up the convergence of the iteration.

Therefore, let us define the array manifold as

$$\begin{array}{cc}\hfill \mathbf{a}\left({\theta }_0\right)& ={[1,e^{-j2\pi \frac{d}{\lambda }sin{\theta }_0},\cdots ,e^{-j2\pi (N-1)\frac{d}{\lambda }sin{\theta }_0}]}^T\hfill \\ \hfill & ={[1,e^{-j\varphi },e^{-j2\varphi },\cdots ,e^{-j(N-1)\varphi }]}^T,\hfill \end{array}$$

(34)

where

$$\begin{array}{c}\hfill \varphi =2\pi \frac{d}{\lambda }sin{\theta }_0\end{array}$$

(35)

and the initial vector ${\mathbf{v}}^0$ is assumed to be

$$\begin{array}{c}\hfill {\mathbf{v}}^0={[b_0,b_1,\cdots ,b_{N-1}]}^T,\end{array}$$

(36)

The orthogonality equation can be expressed as

$$\begin{array}{c}\hfill \mathbf{a}{\left({\theta }_0\right)}^T{\mathbf{v}}^0=b_0+b_1{e}^{-j\varphi }+...+b_{N-1}{e}^{-j(N-1)\varphi }\end{array}$$

(37)

In order to make the incident wave direction not orthogonal to the initial vector, $\mathbf{a}{\left({\theta }_0\right)}^T{\mathbf{v}}^0\ne 0$ is constant. In the following, eight special initial vectors are discussed and their convergence performance is analyzed.

I. The initial vector is assumed to be

$$\begin{array}{c}\hfill {\mathbf{v}}^0={[1,1,1,\cdots ,1]}^T.\end{array}$$

(38)

where all elements are 1. Based on (37), the corresponding orthogonality equation is given by

$$\begin{array}{cc}\hfill \mathbf{a}{\left({\theta }_0\right)}^T{\mathbf{v}}^0& =1+e^{-j\varphi }+e^{-j2\varphi }+...+e^{-j(N-1)\varphi }\hfill \\ \hfill & =\frac{1-{\left(e^{-j\varphi }\right)}^N}{1-e^{-j\varphi }}\hfill \end{array}$$

(39)

where $\varphi =2\pi \frac{d}{\lambda }sin{\theta }_0$ and it can be discussed in three cases.

(1) When $1-e^{j\varphi }\ne 0$ and $1-{\left(e^{-j\varphi }\right)}^N=0$, $\varphi$ can be calculated as

$$\begin{array}{c}\hfill \varphi =2\pi \frac{d}{\lambda }sin{\theta }_0=\frac{2\pi Z}{N}\end{array}$$

(40)

where Z is an integer. It can be seen from (40) that $\mathbf{a}{\left({\theta }_0\right)}^T{\mathbf{v}}^0=0$ holds when ${\theta }_0=arcsin\left(\frac{\lambda Z}{Nd}\right)$ and $\frac{Z}{N}$ is not an integer.

(2) When $1-e^{j\varphi }=0$, $\varphi$ is given by

$$\begin{array}{c}\hfill \varphi =2\pi \frac{d}{\lambda }sin{\theta }_0=2\pi Z\end{array}$$

(41)

Substituting $\varphi =2\pi Z$ into (39), $\mathbf{a}{\left({\theta }_0\right)}^T{\mathbf{v}}^0\ne 0$ is verified to hold.

(3) When $1-e^{j\varphi }\ne 0$ and $1-{\left(e^{-j\varphi }\right)}^N\ne 0$, it is clear that $\mathbf{a}{\left({\theta }_0\right)}^T{\mathbf{v}}^0\ne 0$ holds. Through the above discussion, initial vector I is not guaranteed non-orthogonal always true.

II. The initial vector is assumed to be

$$\begin{array}{c}\hfill {\mathbf{v}}^0={[1,-1,1,-1,\cdots ]}^T.\end{array}$$

(42)

where the odd element is 1 and the even elements is $-1$. Depending on the value of N, two cases can be discussed.

(1) N is even

$$\begin{array}{cc}\hfill \mathbf{a}{\left({\theta }_0\right)}^T{\mathbf{v}}^0& =1-e^{-j\varphi }+e^{-j2\varphi }+...-e^{-j(N-1)\varphi }\hfill \\ \hfill & =(1+e^{-j2\varphi }+e^{-j4\varphi }+...+e^{-j(N-2)\varphi })\hfill \\ \hfill & -(e^{-j\varphi }+e^{-j3\varphi }+...+e^{-j(N-1)\varphi })\hfill \\ \hfill & =\frac{1-{\left(e^{-j2\varphi }\right)}^{\frac{N}{2}}}{1-e^{-j2\varphi }}-\frac{e^{-j\varphi }(1-{\left(e^{-j2\varphi }\right)}^{\frac{N}{2}})}{1-e^{-j2\varphi }}\hfill \\ \hfill & =\frac{(1-e^{-j\varphi })(1-{\left(e^{-j2\varphi }\right)}^{\frac{N}{2}})}{1-e^{-j2\varphi }}\hfill \end{array}$$

(43)

From the above discussion, it is clear that $\mathbf{a}{\left({\theta }_0\right)}^T{\mathbf{v}}^0=0$ holds when $\varphi$ satisfies the following form

$$\left\{\begin{array}{c}i.\varphi =2\pi N,Nisinteger\hfill \\ ii.\varphi =\frac{2\pi Z}{N},\frac{Z}{N}isnotinteger.\hfill \end{array}\right.$$

(44)

If $\varphi$ selects these cases, it may generate incorrect DOA estimation.

(2) N is odd

$$\begin{array}{cc}\hfill \mathbf{a}{\left({\theta }_0\right)}^T{\mathbf{v}}^0& =1-e^{-j\varphi }+e^{-j2\varphi }+...+e^{-j(N-1)\varphi }\hfill \\ \hfill & =(1+e^{-j2\varphi }+e^{-j4\varphi }+...+e^{-j(N-1)\varphi })\hfill \\ \hfill & -(e^{-j\varphi }+e^{-j3\varphi }+...+e^{-j(N-2)\varphi })\hfill \\ \hfill & =\frac{1-{\left(e^{-j2\varphi }\right)}^{\frac{N+1}{2}}}{1-e^{-j2\varphi }}-\frac{e^{-j\varphi }(1-{\left(e^{-j2\varphi }\right)}^{\frac{N-1}{2}})}{1-e^{-j2\varphi }}\hfill \\ \hfill & =\frac{(1-e^{-j\varphi })(1+e^{-jN\varphi })}{1-e^{-j2\varphi }}\hfill \end{array}$$

(45)

It can be seen that $\mathbf{a}{\left({\theta }_0\right)}^T{\mathbf{v}}^0=0$ holds when $\varphi$ satisfies the following conditions

$$\begin{array}{c}\hfill \varphi =\frac{(2n+1)\pi }{N},\end{array}$$

(46)

where $\frac{2n+1}{N}$ is not an integer and orthogonality only can be avoided by selecting other $\varphi$.

III. The initial vector is assumed to be

$$\begin{array}{c}\hfill {\mathbf{v}}^0={[1,1,1,1,\cdots ,-1,-1,\cdots ,-1]}^T.\end{array}$$

(47)

Assuming N is even. The first half of the elements of vector ${\mathbf{v}}^0$ are 1 and the rest are $-1$.

$$\begin{array}{cc}\hfill \mathbf{a}{\left({\theta }_0\right)}^T{\mathbf{v}}^0& =1+e^{-j\varphi }+...-e^{-j(N-1)\varphi }\hfill \\ \hfill & =(1+e^{-j\varphi }+e^{-j2\varphi }+...+e^{-j(\frac{N}{2}-1)\varphi })\hfill \\ \hfill & -(e^{-j\frac{N}{2}\varphi }+...+e^{-j(N-1)\varphi })\hfill \\ \hfill & =\frac{1-{\left(e^{-j\varphi }\right)}^{\frac{N}{2}}}{1-e^{-j\varphi }}-\frac{e^{-j\frac{N}{2}\varphi }(1-{\left(e^{-j\varphi }\right)}^{\frac{N}{2}})}{1-e^{-j\varphi }}\hfill \\ \hfill & =\frac{(1-e^{-j\frac{N}{2}\varphi })(1-{\left(e^{-j\varphi }\right)}^{\frac{N}{2}})}{1-e^{-j\varphi }}\hfill \end{array}$$

(48)

Therefore, $\mathbf{a}{\left({\theta }_0\right)}^T{\mathbf{v}}^0=0$ holds when $\varphi$ satisfies the following cases

$$\left\{\begin{array}{c}i.\varphi =2\pi N,Nisinteger\hfill \\ ii.\varphi =\frac{4\pi Z}{N},\frac{2Z}{N}isnotinteger.\hfill \end{array}\right.$$

(49)

The above value of $\varphi$ will lead to orthogonality.

IV. Let us define $v_0$ as the normalized DFT vector $\mathbf{f}$ with its elements given by [41]

$$\begin{array}{cc}\hfill {\mathbf{v}}^0=\mathbf{f}& ={[e^{-j\frac{2\pi }{N}/\sqrt{N}},e^{-j\frac{2\pi }{N}2/\sqrt{N}},\cdots ,e^{-j\frac{2\pi }{N}N/\sqrt{N}}]}^T\hfill \\ \hfill & ={[e^{-jA},e^{-j2A},\cdots ,e^{-jNA}]}^T.\hfill \end{array}$$

(50)

where

$$\begin{array}{c}\hfill A=\frac{2\pi }{N}/\sqrt{N}\end{array}$$

(51)

Therefore, (37) can be calculated as

$$\begin{array}{cc}\hfill \mathbf{a}{\left({\theta }_0\right)}^T{\mathbf{v}}^0& =e^{-jA}+e^{-j2A}{e}^{-j\varphi }+...+e^{-jNA}{e}^{-j(N-1)\varphi }\hfill \\ \hfill & =e^{-jA}+e^{-j(2A+\varphi )}+...+e^{-j(NA+(N-1\left)\varphi \right)}\hfill \\ \hfill & =\frac{e^{-jA}(1-{\left(e^{-jA}{e}^{-j\varphi }\right)}^N)}{1-\left(e^{-jA}{e}^{-j\varphi }\right)}\hfill \end{array}$$

(52)

Solving the above equation yields

$$\begin{array}{cc}\hfill & \varphi =\frac{2\pi Z}{N}-A=2\pi \frac{d}{\lambda }sin{\theta }_0\hfill \\ \hfill & sin{\theta }_0=\frac{\lambda (Z\sqrt{N}-1)}{N\sqrt{N}d}.\hfill \end{array}$$

(53)

From the above discussion, when $sin{\theta }_0=\frac{\lambda (Z\sqrt{N}-1)}{N\sqrt{N}d}$, where $\frac{Z}{N}$ is not an integer, we have $\mathbf{a}{\left({\theta }_0\right)}^T{\mathbf{v}}^0=0$; the orthogonality may lead to estimation failure.

V. The initial vector is assumed to be

$$\begin{array}{c}\hfill {\mathbf{v}}^0={[\cdots ,1,0,0,0,\cdots ,1,\cdots ]}^T,m,n=1,\cdots ,N.\end{array}$$

(54)

where the m-th element is 1, the n-th is 1 and the rest of the elements are 0 ($m<n$).

$$\begin{array}{c}\hfill \mathbf{a}{\left({\theta }_0\right)}^T{\mathbf{v}}^0=e^{-jm\varphi }+e^{-jn\varphi }=e^{-jm\varphi }(1+e^{-j(n-m)\varphi })\end{array}$$

(55)

When $\varphi =\frac{(2Z+1)\pi }{n-m}$, $\mathbf{a}{\left({\theta }_0\right)}^T{\mathbf{v}}^0=0$.

VI. The initial vector is assumed to be

$$\begin{array}{c}\hfill {\mathbf{v}}^0={[\cdots ,1,0,\cdots ,1,\cdots ,1,\cdots ]}^T,m,n,k=1,\cdots ,N.\end{array}$$

(56)

where $m,n,k$-th is 1 and the rest of the elements are 0 ($m<n<k$).

$$\begin{array}{cc}\hfill \mathbf{a}{\left({\theta }_0\right)}^T{\mathbf{v}}^0& =e^{-jm\varphi }+e^{-jn\varphi }+e^{-jk\varphi }\hfill \\ \hfill & =e^{-jm\varphi }(1+e^{-j(n-m)\varphi }(1+e^{-j(k-n)\varphi }))\hfill \end{array}$$

(57)

when $\varphi =\frac{Z\pi }{n-m}$ or $\frac{(2Z+1)\pi }{k-n}$, $\mathbf{a}{\left({\theta }_0\right)}^T{\mathbf{v}}^0=0$.

VII. The initial vector is assumed to be

$$\begin{array}{cc}\hfill & {\mathbf{v}}^0={[\cdots ,1,0,\cdots ,1,\cdots ,1,\cdots ,1,\cdots ]}^T,\hfill \\ \hfill & m,n,k,l=1,\cdots ,N.\hfill \end{array}$$

(58)

where $m,n,k,l$-th is 1 and the rest of the elements are 0 ($m<n<k<l$).

$$\begin{array}{cc}\hfill \mathbf{a}{\left({\theta }_0\right)}^T{\mathbf{v}}^0& =e^{-jm\varphi }+e^{-jn\varphi }+e^{-jk\varphi }+e^{-jl\varphi }\hfill \\ \hfill & =e^{-jm\varphi }.\hfill \\ \hfill & (1+e^{-j(n-m)\varphi }(1+e^{-j(k-n)\varphi }(1+e^{-j(l-k)\varphi })))\hfill \end{array}$$

(59)

when $\varphi =\frac{Z\pi }{n-m}$ or $\frac{Z\pi }{k-n}$ or $\frac{(2Z+1)\pi }{l-k}$, $\mathbf{a}{\left({\theta }_0\right)}^T{\mathbf{v}}^0=0$.

So we can obtain the general rule, when $b_1,b_2,b_3,\dots ,b_n$ is 1 and the rest of the elements are 0 ($b_1<b_2<\dots <b_n$), when $\varphi =\frac{Z\pi }{b_2-b_1}$ or $\frac{Z\pi }{b_3-b_2}$ or …$\frac{Z\pi }{b_{n-1}-b_{n-2}}$ or $\frac{(2Z+1)\pi }{b_n-b_{n-1}}$, $\mathbf{a}{\left({\theta }_0\right)}^T{\mathbf{v}}^0=0$.

VIII. The initial vector is assumed to be

$$\begin{array}{c}\hfill {\mathbf{v}}^0={\mathbf{e}}_k={[\cdots ,0,0,\cdots ,1,\cdots ]}^T,k=0,\cdots ,N-1.\end{array}$$

(60)

where the k-th element of ${\mathbf{e}}_k$ is 1 and the rest of the elements are 0. $\mathbf{a}{\left({\theta }_0\right)}^T{\mathbf{v}}^0\ne 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 $\epsilon$ also influences the speed of convergence in the convergence process of power iteration, which can be expressed as

$$\begin{array}{c}\hfill \epsilon =\frac{|{\epsilon }_n|}{|{\alpha }_1{\mathbf{u}}_{\mathbf{1}}|}=\frac{|{\displaystyle \sum _{k=2}^N}{\alpha }_k{\left(\frac{{\lambda }_k}{{\lambda }_1}\right)}^n{\mathbf{u}}_{\mathbf{k}}|}{|{\alpha }_1{\mathbf{u}}_{\mathbf{1}}|}\end{array}$$

(61)

and

$$\begin{array}{c}\hfill |\sum _{k=2}^N{\alpha }_k{\left(\frac{{\lambda }_k}{{\lambda }_1}\right)}^n|\le |\sum _{k=2}^N{\alpha }_2{\left(\frac{{\lambda }_2}{{\lambda }_1}\right)}^n|\end{array}$$

(62)

Therefore, (61) can be further expressed as

$$\begin{array}{cc}\hfill & \epsilon \le (N-1){(\frac{|{\lambda }_2|}{|{\lambda }_1|})}^n\frac{{\alpha }_2}{{\alpha }_1},\hfill \\ \hfill & {log}_2\epsilon \le {log}_2(N-1)+n{log}_2(\frac{|{\lambda }_2|}{|{\lambda }_1|})-{log}_2\frac{{\alpha }_1}{{\alpha }_2}\hfill \end{array}$$

(63)

From ${\alpha }_1\ge {\alpha }_2$,

$$\begin{array}{c}\hfill {log}_2\epsilon \le {log}_2(N-1)+n{log}_2(\frac{|{\lambda }_2|}{|{\lambda }_1|})\end{array}$$

(64)

Since $|{\lambda }_1|>|{\lambda }_2|$,

$$\begin{array}{c}\hfill n\le \frac{{log}_2\epsilon -{log}_2(N-1)}{{log}_2(\frac{|{\lambda }_2|}{|{\lambda }_1|})}\end{array}$$

(65)

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 ratio $\frac{|{\lambda }_2|}{|{\lambda }_1|}$, since it dictates the rate of convergence.

The whole procedure is summarized in Algorithm 3.

Algorithm 3 PI method on subspace

Input: matrix $\mathbf{R}$, tolerence $\epsilon$

Output: ${\lambda }_1$, ${\mathbf{v}}^n$ and n

Initialization: choose an initial vector ${\mathbf{v}}^0$, and n = 1.

repeat

1. Calculate ${\mathbf{v}}^n$ = $\mathbf{R}{\mathbf{v}}^{n-1}$;

2. Calculate ${\mathbf{v}}^n$ = ${\mathbf{v}}^n/max\left({\mathbf{v}}^n\right)$;

3. Calculate ${\lambda }^n$ = ${\left[{\mathbf{v}}^n\right]}^T\mathbf{R}{\mathbf{v}}^n$;

4. Update ${\lambda }_1$ = ${\lambda }^n$;

5. n = n + 1.

until

$\parallel \mathbf{max}\left({\mathbf{v}}^{\mathbf{n}}\right)-\mathbf{max}\left({\mathbf{v}}^{\mathbf{n}-\mathbf{1}}\right)\parallel \le$$\epsilon$.

Return ${\lambda }_1$, ${\mathbf{v}}^n$ and n

n is the number of iterations, ${\lambda }_1$ is the dominant eigenvalue and ${\mathbf{v}}^n$ is the dominant eigenvector corresponding to ${\lambda }_1$.

Notice that in each iteration we compute a single matrix–vector multiplication $\left(O\left(N^2\right)\right)$. We never perform matrix–matrix multiplication, which requires a greater number of operations $\left(O\left(N^3\right)\right)$. 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.

4.4. Complexity Analysis

We analyze the computational complexities of the proposed two estimators with traditional ESPRIT and Root-MUSIC algorithms as a complexity benchmark. Thus, the complexity of ESPRIT is as follows

$$\begin{array}{c}\hfill C_{ESPRIT}=O\{{(2N-2)}^3+2N-3\}\end{array}$$

(66)

FLOPs. The complexity of Root-MUSIC is

$$\begin{array}{c}\hfill C_{Root-MUSIC}=O\{N^3+8N^2+NK(2N+3)-11N+4\}\end{array}$$

(67)

FLOPs. The complexity of RPI-RI is

$$\begin{array}{c}\hfill C_{RPI-RI}=O\{\beta {(2N-2)}^2+2N-3\}\end{array}$$

(68)

FLOPs. The complexity of RPI-PR is

$$\begin{array}{c}\hfill C_{RPI-PR}=O\{(\beta +8)N^2+NK(2N+3)-11N+4\}\end{array}$$

(69)

FLOPs, where $\beta$ is the iteration number of the power iteration. Assuming N is far larger than $K,\beta$, 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.

5. Simulation Results and Discussions

In this section, a simulation is presented 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 Monte Carlo simulations L is 1000. Without loss of generality, we assume that UAV emitter located in ${\theta }_0={50}^{\circ }$, the number of snapshots K is 1000, the antenna distance $d=\lambda /2$ and the number of antenna elements N is 64 in the massive MIMO system.

Figure 2 plots the root mean square error (RMSE) curves using four different methods versus SNR with CRLB as the performance benchmark. Observing Figure 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.

Figure 3 and Figure 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\in [16,272]$ and $K\in [100,3600]$. From Figure 3 and Figure 4, it can be seen that the DOA estimators on the power iteration method achieve a performance close to the conventional algorithms for any number of antennas and snapshot scenarios; in particular, RPI-PR can reach the FD CRLB.

Figure 5 plots curves of complexity analysis versus the number of antennas N with $K=50$, SNR = 0 dB. 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.

To explore the influence of the selection of the initial vector on the number of iterations and the convergence speed, Figure 6 shows the relationship between the optimal eigenvector value ${\mathbf{v}}^{\mathbf{n}}$ and the number of iterations n, given three different initial vectors ${\mathbf{v}}^{\mathbf{0}}$. When the initial vector ${\mathbf{v}}^{\mathbf{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 (30). In addition, when the ${\mathbf{v}}^{\mathbf{0}}$ obeys a random vector distribution, the required number of iterations n is increased, which requires an average of eight iterations to reach the convergence of the $max\left(v^n\right)$.

In order to verify the convergence performance of the initial vector ${\mathbf{v}}^0={\mathbf{e}}_k={[\cdots ,0,0,\cdots ,1,\cdots ]}^T$ in (60), five incident wave directions ${\theta }_0=\{{30}^{\circ },{60}^{\circ },{90}^{\circ }$, ${120}^{\circ },{150}^{\circ }\}$ are assumed in Figure 7 and the number of iterations of RPI methods is given by numerical simulation when SNR = 0. As shown in Figure 7, the proposed RPI methods require about five iterations to converge. The $|{\mathbf{v}}_n|$ in the figure indicates the absolute value of the vector at the nth iteration.

To analyze the effect of SNR on the convergence speed of the proposed RPI methods, Figure 8 plots the iteration error ${\alpha }_n$ of the RPI methods versus the number of iterations, where the iteration error ${\alpha }_n$ represents the absolute value of the difference between the result of the nth iteration $C\left({\theta }_n\right)$ and the result of the $(n-1)$th iteration $C\left({\theta }_{n-1}\right)$, i.e., ${\alpha }_n=|C\left({\theta }_n\right)-C\left({\theta }_{n-1}\right)|$. As shown in Figure 8, it is assumed that the number of antennas $N=1024$ and the incident wave direction ${\theta }_0={50}^o$. The number of iterations is 7, 5 and 4 for SNRs of −30 dB, 0 dB and 30 dB, respectively. the number of iterations of the proposed algorithm decreases sequentially as the SNR increases.

6. Conclusions

To find the direction of UAV emitters 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 significantly reduce the computational complexity compared to the conventional algorithms and achieve fast convergence of the power iteration process. In particular, the RPI-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. However, there are still some issues that deserve to be explored in depth, for example, how to find an initial iteration vector that is closer to the maximum eigenvector to make the iteration speed up further. In addition, it is a challenging problem to obtain the signal subspace via the power iteration method in the scenario of multiple radiation sources and these problems will be the research directions we focus on in the future.