Geometric Algebra-Based ESPRIT Algorithm for DOA Estimation

Direction-of-arrival (DOA) estimation plays an important role in array signal processing, and the Estimating Signal Parameter via Rotational Invariance Techniques (ESPRIT) algorithm is one of the typical super resolution algorithms for direction finding in an electromagnetic vector-sensor (EMVS) array; however, existing ESPRIT algorithms treat the output of the EMVS array either as a “long vector”, which will inevitably lead to loss of the orthogonality of the signal components, or a quaternion matrix, which may result in some missing information. In this paper, we propose a novel ESPRIT algorithm based on Geometric Algebra (GA-ESPRIT) to estimate 2D-DOA with double parallel uniform linear arrays. The algorithm combines GA with the principle of ESPRIT, which models the multi-dimensional signals in a holistic way, and then the direction angles can be calculated by different GA matrix operations to keep the correlations among multiple components of the EMVS. Experimental results demonstrate that the proposed GA-ESPRIT algorithm is robust to model errors and achieves less time complexity and smaller memory requirements.


Introduction
Direction-of-arrival (DOA) estimation of electromagnetic (EM) signals has attracted wide attention in many communication fields, such as radar [1,2], mobile networks [3] and sonar [4]. It is clear that DOA estimation is the basic and essential part in an array signal processing system. For example, a corresponding transmitting or receiving beamformer can be designed to extract signals in the direction of interest and suppress uninteresting interference signals. The electromagnetic vector sensor (EMVS) can catch polarization-related information compared to a conventional scalar sensor, which can further improve the target resolution, anti-interference ability and detection stability for DOA estimation [5][6][7]; therefore, the research for EMVS array direction finding has become a hotspot.
With the appearance of the Long-Vector MODEL (LV-MODEL) [5] (built for EMVS), multiple researchers have proposed various DOA estimators. The existing estimators can be summarized into three categories: (1) research on DOA estimators transplanting from scalar sensor; (2) research based on special array arrangement; (3) research based on advanced mathematical tools.
In terms of transplantation, the classic subspace-based super-resolution algorithm [8] (Multiple Signal Classification-MUSIC) was transplanted to the EMVS [9][10][11] array, but the algorithms often suffer high computational complexity because of the four-dimensional parameter search for two direction angles and two additional polarization angles; therefore, Weiss [12] used the polynomial root to reduce the computational complexity to a certain extent. In addition, another subspace-based super-resolution algorithm [13,14] (Estimation of Signal Parameters via Rotational Invariance Techniques-ESPRIT) was also transplanted 1.
We incorporate the multi-dimensional consistency of GA into ESPRIT, and propose a Geometric Algebra-based ESPRIT algorithm (GA-ESPRIT) for 2D-DOA estimation.

2.
We use the new calculation rules of the high-dimensional algebra system to preserve the correlation among multiple components of EMVS.

3.
Experimental results demonstrate that the proposed GA-ESPRIT algorithm can achieve more accurate, stable and lighter DOA estimation.
The rest of this paper is organized as follows. Section 2 introduces the basics of GA and the EMVS model for narrow-band signals based on GA. Section 3 describes the proposed GA-ESPRIT in detail. Experimental results and analysis are provided in Section 4, followed by concluding remarks in Section 5.

Fundamental of Geometric Algebra
The concept of GA [30] was proposed by David Hestenes in the 1960s, who combined Clifford Algebra with a physical geometric structure. After decades of research, GA has shown its absolute superiority in electromagnetism [31], cosmology [32], multi-channel image [33][34][35] and other physical sciences. The crucial product operation in GA theory is the geometric product [30]. For vectors a and b, the geometric product is denoted by where {·} and {∧} denote the inner product and the outer product, respectively.

Multi-Vector
Let G n = C n,0 , which is the real GA of the quadratic pair (V, Q) where V = R n and Q is the quadratic form of signature (n, 0). There is an orthogonal basis {e 1 , e 2 , . . . , e n } in R n , which generates 2 n basis elements of G n via the geometric product as shown in (2): , . . . , {e 1 e 2 · · · e n } k=n (2) for i, j = 1, 2, . . . , n.
The multi-vector A of G n is defined as where The reverse of multi-vector A is defined as where A k for k = 1, 2, 3, . . . , 7 are all m × n real number matrices. The transpose with reversion of A is denoted by A H where A T i for k = 1, 2, 3, . . . , 7 denotes the transpose.

G-MODEL
A compact polarized GA model for the vector-sensor array was proposed in [21], named G-MODEL, which models the six-component outputs of a vector sensor holistically using a multi-vector in G 3 . Suppose there are K narrow-band, far-field and uncorrelated sources with wavelength λ impinging on an array, which includes Q vector sensors. Define θ k ∈ [0, 2π), φ k ∈ [0, π), γ k ∈ [0, π/2) and η k ∈ [−π, π) are the azimuth angle, elevation angle, polarization amplitude angle and phase difference angle of the kth source, respectively.
Define u k = cos θ k sin φ k e 1 + sin θ k sin φ k e 2 + cos φ k e 3 as the unit vector (see Figure 1) of the kth source when it impinges on the sensor at the origin. v k1 = − sin θ k e 1 + cos θ k e 2 and v k2 = cos θ k cos φ k e 1 + sin θ k cos φ k e 2 − sin φ k e 3 are unit multi-vectors. The position vector of the qth sensor is r q = r q1 e 1 + r q2 e 2 + r q3 e 3 . The output of the qth vector sensor in the array is denoted by [21] Y (q) where X q (θ k , φ k ) = e e 123 2π λ (cos θ k sin φ k r q1 +sin θ k sin φ k r q2 +cos φ k r q3 ) is the spatial phase factor of the kth source incident on the qth vector sensor.
In next section, the GA-ESPRIT algorithm is deduced based on the G-MODEL.
Define u k = cos θ k sin φ k e 1 + sin θ k sin φ k e 2 + cos φ k e 3 as the uni of the kth source when it impinges on the sensor at the origin. v k1 = and v k2 = cos θ k cos φ k e 1 + sin θ k cos φ k e 2 − sin φ k e 3 are unit multi-v vector of the qth sensor is r q = r q1 e 1 + r q2 e 2 + r q3 e 3 . The output of th the array is denoted by [21] Y (q) where X q (θ k , φ k ) = e e 123 2π λ (cos θ k sin φ k r q1 +sin θ k sin φ k r q2 +cos φ k r q3 ) is the sp the kth source incident on the qth vector sensor.
In next section, the GA-ESPRIT algorithm is deduced based on

Proposed Algorithm
The basic premise of the ESPRIT algorithm is that there are id spacing between subarrays is known and the structure of subarra satisfies the rotational invariance in space [13]. Uniform linear arrays it comes to one-dimensional DOA estimation using conventional ESP with ULAs, double parallel uniform linear arrays (DPULAs) can iden

Proposed Algorithm
The basic premise of the ESPRIT algorithm is that there are identical subarrays, the spacing between subarrays is known and the structure of subarrays is identical, which satisfies the rotational invariance in space [13]. Uniform linear arrays (ULAs) appear when it comes to one-dimensional DOA estimation using conventional ESPRIT [1,13]. Compared with ULAs, double parallel uniform linear arrays (DPULAs) can identify two-dimensional DOA because of the special construction, which consists of two parallel ULAs [36][37][38]; therefore, the algorithm discussed in this paper is based on DPULAs.

Complex Representation Matrix and Related Calculations
In view of the paucity of research on calculations with multi-vector, the Complex Representation Matrix (CRM) [20] is introduced because of the mature matrix theories. Consider a matrix A ∈ G m×n 3 , the CRM is defined by Ψ(A) Let ν = (−e 1 + e 13 )/2 ∈ G 3 , and its reversion isν = (−e 1 − e 13 )/2 ∈ G 3 . Then, which imply ννν = ν,ννν =ν, (νν) 2 = νν, (νν) 2 =νν. It immediately follows that, for every A ∈ G 3 , we have where in (10) and (11) we have I k denotes the k × k identity matrix. It is not difficult to prove that where {+} denotes the pseudo-inverse. Referring to (10) and (14c), the pseudo-inverse of any A ∈ G 3 is Since e 2 123 = −1 and e 123 commutes with all elements in G 3 , one can identify it with the complex imaginary unit j [20], and so we can view Ψ(A) given in (8) as a complex matrix.

Model for DPULAs
Consider a DPULA with 2M + 2 sensors, as shown in Figure 2, in which d and M refer to the spacing between two adjacent sensors and the number of sensors in per subarray, respectively. The array is divided into three subarrays. The 1st to Mth sensors on the x-axis compose the first subarray, the 2nd to (M + 1)th sensors form the second subarray and the (M + 2)th to (2M + 1)th that located on a straight line parallel to the x-axis make up the third subarray. The reason for the division can be found in Figure 3, that is, there are two unknown DOA parameters in the model, which need two rotational invariance relations.  Since the three subarrays have the same structure each output of them has only one phase difference for the subarray one, two and three are defined as Y 1 EH , Y 2 EH and the above array model, the outputs of the three subarrays Since the three subarrays have the same structure and the same number of sensor, each output of them has only one phase difference for the same signal. Signals received by subarray one, two and three are defined as Y 1 EH , Y 2 EH and Y 3 EH , respectively. According to the above array model, the outputs of the three subarrays at time t are as follows where and According to (18), we find that the DOA information is contained in matrix A, F and G. Because F and G are diagonal matrices that only contain direction information of incident signals, the focus is the two matrices, i.e., Clearly, it is easy to figure out the DOA in the light of (19) if we obtain the two ideal matrices F and G. From the rules of subarray division, we can see that the latter (M − 1) sensors of subarray one and the former (M + 1) sensors of subarray two are overlapped. Thus, in order to reduce the computational complexity, subarray one and subarray two can be merged to form a new matrix Y EH , that is, After merging, the (2M + 2)th redundant sensor is added to subarray three to form a new subarray P EH , so that the third subarray has the same dimension as Y EH Let A be the array flow pattern of Y EH , then Y EH and P EH can be written as where Then, B(t) is defined as where Finally, the output of the whole array is denoted by

Algorithm Details
It is assumed that the sources received by the vector-sensor array are random signals which are independent and uncorrelated. In the same way, the measuring noise on six antennas of each sensor is white noise with the same power.

Subspace Separation
Under the above assumption, theoretically, the covariance matrix of the array output is given by where E{·} stands for the mathematical expectation operator, σ 2 is the noise power on each vector antenna, R S = E S(t)S H (t) .
Since the geometric product is non-commutativity, the Eigenvalue Decomposition (ED) is different from the conventional real methods but similar to the quaternion case. In other words, there are two possible eigenvalues, namely the left and the right eigenvalue for G 3 matrix. In the proposed algorithm, the right eigenvalue is selected because the right ED of G 3 matrix can be converted to the right ED of its CRM [20].
The ED of R is denoted by According to the principle of subspace separation, U s is the signal subspace corresponding to K larger eigenvalues, and Σ s is a diagonal matrix composed of K larger eigenvalues. In addition, U n is orthogonal to U s and it is the noise subspace corresponding to the remaining 4(M + 1) − K small eigenvalues. Similarly, Σ n is a diagonal matrix composed of the remaining small eigenvalues.
In the actual processing, the received signal is usually sampled. So, for a certain number of snapshots N, (26) and (27) can be rewritten aŝ Because the space formed by the eigenvectors corresponding to the larger eigenvalues is the same as the space formed by the steering multi-vectors of the incident signals, that is, span{U s } = span{C}, there exists a unique non-singular matrix T, which satisfies The rotational invariance relations exist among three subarrays, but U s is the signal subspace of the whole array; therefore, after obtaining U s , the signal subspace of three subarrays must be separated. By the arrangement of sensor array, we find that the signal subspace of three subarrays can be calculated by where U s1 , U s2 and U s3 are signal subspaces of subarray one, subarray two and subarray three, respectively.

Rotation Invariance
From (30), the pivotal matrices F and G can be found. So, let in the same way, It is discovered that the eigenvalues of Ψ x and Ψ y are diagonal elements of F and G, respectively.
Equations (32) and (33) are equations themselves, and are usually solved by the Least Squares (LS) method [29,36,38,39]; however, LS only takes the error on the left side of the equation into account, it ignores that the coefficient matrix also has an error; therefore, in order to reduce the error caused by solving the equation as much as possible, this paper considers a more accurate method-TLS [13]. Next, the solution of the equation is obtained by taking (32) as an example.
Combining the idea of TLS with the orthogonal property of subspace, we define a new matrix U s12 = [U s1 U s2 ]. In fact, the main aim is to seek a unitary matrix D ∈ G M×2K 3 , which is orthogonal to U s12 . In other words, the space formed by D is orthogonal to the space formed by the column vectors of U s1 or U s2 . So the D can be obtained from the ED where Λ is the diagonal matrix whose diagonal elements are composed by K multi-vectors that only have 0-grade-vector (can regard as non-zero real number) and 3K multi-vectors that equal to 0. E can be written as Let E N = E 12 E 22 , which is composed by eigenvectors whose eigenvalues are 0 and form the noise subspace. Since U s12 is signal subspace, we find that D = E N , i.e., Then, The pseudo-inverse of G 3 matrix E 22 can be found in (15).

Angle Estimation
The azimuth and elevation angle of K signals are included in F and G. In theory, the eigenvectors obtained by ED of these two matrices are both T; however, in the actual calculation process, the two eigenvalue decomposition operations are carried out independently, which can not ensure that the arrangement of eigenvectors in them is reflected well; therefore, the diagonal elements of F and G should be matched.
Suppose that T1 and T2 are eigenvector matrices derived from GA-ED of Ψ x and Ψ y , respectively. Then where {| · |} is the operator that gets magnitude of every multi-vector in a matrix. For the same signal, the eigenvectors in T1 and T2 corresponding to matched f k and g k are related; therefore, the order of diagonal elements in F and G can be adjusted by the coordinate of the largest element in each row (or column) of O to complete matching. After observing (19), f and g are multi-vectors that only have scalar and 3-gradevector, if we replace e 123 with the imaginary unit j of complex number, f and g can be regarded as complex numbers. Finally, we calculate θ k and φ k with f k and g k , that is, where angle(·) is the operator for getting phase angle. In conclusion, the steps of the GA-ESPRIT algorithm are: 1.
The original data received from three subarrays are integrated into the measurement model of the whole array according to (25); 2.
Calculate the covariance matrixR, and then the ED in GA ofR is performed and the signal subspace U s can be obtained by the larger eigenvalues; 3.
According to (30), the signal subspace U s of the whole array is divided into three subspaces U s1 , U s2 and U s3 ; 4. Ψ x and Ψ y can be obtained using TLS in GA, and details can be found in (34)-(37); 5.
The ED of Ψ x and Ψ y is performed to obtain matrices F and G; 6.
The eigenvalues are matched in line with (38) and then taken them into Equation (39) to calculate K pairs direction angles.
Further, the corresponding relationship between the logic flow and steps of GA-ESPRIT is shown in Figure 4.

Complexity Analysis
As discussed in [21,22,27], the estimation of the data covariance matrix is an important factor to illustrate the complexity of ESPRIT algorithm and another one is ED, because they imply many repetitive operations and results, which mean heavy computational burden and memory requirements. Thus, we evaluate the time complexity of the two processes and space complexity in terms of real value memory requirements.
Suppose that an array composed of M vector sensors, and N snapshots are taken. LV-ESPRIT [13] and GA-ESPRIT consider six-component measurements of each vector sensor, whereas Q-ESPRIT [25] only records two-component measurements (electric field on x-axis and y-axis.); therefore, we compare the complexity between LV-ESPRIT and GA-ESPRIT. The output of each vector sensor for each signal consists of six complex numbers in LV-ESPRIT, while GA-ESPRIT only has one multi-vector with vector and bivector parts.
The geometric product of two multi-vectors received by the array output implies 36 real multiplications [21], which is nine times as many real multiplications as two complex numbers. As mentioned in Section 2, the ED of a G 3 matrix is calculated by its CRM; therefore, the time complexity of the two algorithms is shown in Table 1. As for space complexity, the memory requirements of a real number is used to measure [21]. In the following two tables, CM is the Covariance Matrix, R represents real number.
The complexity comparison of these two algorithms can be found in Table 2, where CM and R represent covariance matrix and real number, respectively. Observing the time complexity in Table 2, it is not difficult to find that the computational burdens of CM and ED in GA-ESPRIT are a quarter and 1/27 of these in LV-ESPRIT, respectively. As for space complexity, GA-ESPRIT achieves such a significant reduction, more than 1.5 times compared to LV-ESPRIT, which means that the memory pressure is alleviated, especially for the large data size. The reason for the above comparison results is the natural advantage of GA matrix operations. In detail, because the six-dimensional measurement data in LV-MODEL (stored as 12 real numbers) are mapped into a multi-vector in the G-MODEL (stored as six real numbers), the amount of calculation will be reduced to varying degrees with different matrix operations, which will also bring fewer data storage requirements.
The superior description and calculation ability of GA for multi-dimensional signals make GA-ESPRIT a very notable method for direction finding.

Method
Time Complexity Space Complexity

Simulation Results and Analysis
In this section, we simulate and analyze the proposed GA-ESPRIT based on DPULAs with d = λ/2, discuss its feasibility and performance compared with LV-ESPRIT [14] (in complex number field) and Q-ESPRIT [24] (in quaternion field). The estimation accuracy is evaluated by Root Mean Square Error (RMSE), which is calculated by the average of 200 Monte Carlo simulation experiments.
The RMSE of DOA estimation is defined as where K, ∆θ k and ∆φ k denote the number of incident signals and errors between the result calculated by DOA algorithm and direction angle initially defined in the experiment, respectively. In actual applications, the sensor model errors [27,[41][42][43] cannot be ignored, which main include sensor-position error, gain error and phase error. The sensor-position error, as defined in [27], is the error between the actual position and the ideal position of each vector sensor. In the simulation experiment, the sensor-position error is modeled as additive noise with uniform distribution in a certain range, that is, where k m and k m are the actual position and ideal position of the mth sensor in vectorsensor array, respectively. ε mx , ε my and ε mz are uniformly distributed noise terms. P pe represents the perturbation power of sensor-position error and the larger P pe means the greater deviation of the sensor from its ideal position. Further, referring to [43], the array output with the gain and phase error is denoted by where Π = diag(η 1 , η 2 , . . . , η 2M+2 ), Ξ = diag(exp(e 123 ξ 1 ), . . . , exp(e 123 ξ 2M+2 )), in which η i and ξ i (i = 1, 2, 3, . . . , 2M + 2) are gain error and phase disturbance, respectively. In this paper, we also model them as additive noise. In addition, the six components of all EMVSs are added with noise according to the Signal-to-Noise ratio (SNR) in the following experiments. The SNR is defined as SNR = 10lg(P s /P n ), in which P s and P n are the power of signal and noise on each component, respectively.
In the first experiment, we consider three far-field, narrow-band and uncorrelated signals with parameters Γ = {160  Figure 5a shows the estimation results of three algorithms when ideal Gaussian white noise is added, whereas, the noise in Figure 5b is related. It can be concluded that the three algorithms have very close accuracy of calculating DOA at high levels of SNR from Figure 5a,b, while with the lower SNR, GA-ESPRIT has higher accuracy over the other two and can achieve remove the correlation of noise partially.  Figure 6a shows the performance of the three algorithms when sensor-position error exists with different intensities. Meanwhile, we set SNR to 10 dB and the snapshot number is 200. The sensor-position error of the array sensor is changed by the value of P pe , whose range is 0-0.07. It can be seen in Figure 6b that we fix P pe = 0.02 to observe the estimation of the three algorithms by altering SNR from −10 dB to 20 dB. Figure 6a,b both imply that accuracy of GA-ESPRIT is highest in the presence of the sensor-position error, so the conclusion is that GA-ESPRIT has the strongest robustness against sensor-position errors among the three algorithms. The third experiment is also designed for two cases. Case one is that only gain error exists (see Figure 7a), while for case two, only phase error exists (see Figure 7b). Other conditions are the same as experiment two except that there is no position error. The gain error is constructed by the random numbers, whose mean value is 1 and variance is 0.2, and the phase error is constructed by the random numbers with zero-mean and 0.005 variance. We can learn from Figure 7a,b that, whether there is gain error or phase error, GA-ESPRIT can maintain the estimation accuracy very well, especially in low SNR. In general, it is because Q-ESPRIT only takes part of the array output information into consideration that makes large RMSE. The reason for LV-ESPRIT's poor accuracy in the face of the sensor-model error would be that its "long vector" destroys the orthogonality of the signal components. The improvement of detection robustness of GA-ESPRIT largely results from the fact that it can effectively preserve the orthogonality of the signal components and guarantee the completeness of the information.

Conclusions
In this paper, considering that the GA representation contains physical interpretations and complete information of incident signals, we use the idea of the traditional ESPRIT algorithm to find multiple EM signals in the direction finding method in GA. In particular, the model for DPULAs was built in GA and GA-ESPRIT was successfully derived using new calculation rules to achieve the two-dimensional DOA estimation. Compared with the previous ESPRIT algorithms, due to the robustness to sensor-model error and correlated noise, our proposed approach has potential in many practical situations, such as military radar in difficult environments. According to the experimental results, we have confirmed that the GA-ESPRIT has improved accuracy in two-dimension DOA estimation and can resist environmental interference to some extent. More importantly, the proposed algorithm achieves a reduction of more than 1/3 of the memory requirements while the time complexity is also greatly decreased.
Future works on the GA-ESPRIT will include polarization parameter estimation by optimizing matrix operations in GA and the ability of DOA recognition when facing coherent EM signals. It is expected that the proposed GA-ESPRIT will be an efficient DOA estimator.

Conflicts of Interest:
The authors declare no conflict of interest.