DOA and Polarization Estimation Using an Electromagnetic Vector Sensor Uniform Circular Array Based on the ESPRIT Algorithm

In array signal processing systems, the direction of arrival (DOA) and polarization of signals based on uniform linear or rectangular sensor arrays are generally obtained by rotational invariance techniques (ESPRIT). However, since the ESPRIT algorithm relies on the rotational invariant structure of the received data, it cannot be applied to electromagnetic vector sensor arrays (EVSAs) featuring uniform circular patterns. To overcome this limitation, a fourth-order cumulant-based ESPRIT algorithm is proposed in this paper, for joint estimation of DOA and polarization based on a uniform circular EVSA. The proposed algorithm utilizes the fourth-order cumulant to obtain a virtual extended array of a uniform circular EVSA, from which the pairs of rotation invariant sub-arrays are obtained. The ESPRIT algorithm and parameter pair matching are then utilized to estimate the DOA and polarization of the incident signals. The closed-form parameter estimation algorithm can effectively reduce the computational complexity of the joint estimation, which has been demonstrated by numerical simulations.


Introduction
In the past few decades, the direction of arrival (DOA) estimation of incident signals has been demonstrated to play a significant role in array signal processing [1][2][3][4]. Electromagnetic vector sensor arrays (EVSAs) can receive the incident electromagnetic waves in the form of vectors, which include both the polarization domain and the spatial domain information [5,6]. EVSAs have some inherent advantages [7][8][9][10][11] over traditional arrays. Firstly, EVSAs have superior system performances, which the signals can be distinguished based on their polarization characteristics. Secondly, more accurate models of direction finding systems can be established within EVSAs. Finally, EVSAs have stronger ability of anti-fuzzy than scalar sensor arrays, etc.
As a potential solution, EVSA-based techniques have been widely used in many fields, such as radar [12,13], communication [14], sonar [15,16], etc. Various DOA and polarization joint estimation algorithms have been developed recently, including MUSIC [17,18], ESPRIT [19][20][21][22][23], pencil-MUSIC [24,25], Root-MUSIC [26], etc. An enhanced MUSIC algorithm, proposed by Hua, et al. [24], can effectively improve the estimation accuracy. It is quite applicable to arrays with arbitrary geometries. This algorithm, however, involves the spectrum function construction in four dimensions, leading to fairly high computational complexity. ESPRIT algorithm was initially applied to a uniform linear vector sensor array composed of crossed dipoles for multiple-signal joint DOA and polarization parameters estimation in [20][21][22]. Since the ESPRIT-based algorithm is

Array Signal Model
As shown in Figure 1, K non-Gaussian, narrowband, far-field, incoherent plane wave signals imping on a uniform circular array with M array elements, which are uniformly distributed along a circular path with the radius of r. The element located on the positive x axis is denoted as "1"; the remaining elements are uniformly arranged clockwise on the path, successively denoted as "2" to "M". The phase center of each element is always located on the xoy-plane. Each element consists of two dipoles that are spatially co-located and orthogonal to each other, leading to the array composed of N = 2M dipoles.

Array Signal Model
As shown in Figure 1, K non-Gaussian, narrowband, far-field, incoherent plane wave signals imping on a uniform circular array with M array elements, which are uniformly distributed along a circular path with the radius of r . The element located on the positive x axis is denoted as "1"; the remaining elements are uniformly arranged clockwise on the path, successively denoted as "2" to "M". The phase center of each element is always located on the xoy-plane. Each element consists of two dipoles that are spatially co-located and orthogonal to each other, leading to the array composed of 2 N M = dipoles. In Figure 1, η π ∈ , referred to as the polarization auxiliary angle and polarization phase difference, respectively, to completely describe the polarization state of the incident signal.
Here, assuming that the noise is the additive white Gaussian noise and independent of the incident signals. In real direction finding system, in order to facilitate processing, amplitudes of received data are always compensated and the mutual coupling between array elements [47,48] are calibrated before they are provided to direction finding signal processor, so that the received data of each element is regarded as uniform in term of amplitude and effect of the mutual coupling. Therefore, the output of array can be expressed as: 2 ( 1) j2π( sin cos( ))/ j2π( sin cos( ))/ j2π( sin cos( ))/ , ( ) e , ,e , ,e , j sin cos cos cos ( , , ) cos cos sin sin e k In Figure 1, θ ∈ [0, 2π] and ϕ ∈ [0, π/2] are the azimuth angle and elevation angle of the incident signal, respectively. Here we introduce another two angles, γ ∈ [0, π/2] and η ∈ [0, 2π], referred to as the polarization auxiliary angle and polarization phase difference, respectively, to completely describe the polarization state of the incident signal.
Here, assuming that the noise is the additive white Gaussian noise and independent of the incident signals. In real direction finding system, in order to facilitate processing, amplitudes of received data are always compensated and the mutual coupling between array elements [47,48] are calibrated before they are provided to direction finding signal processor, so that the received data of each element is regarded as uniform in term of amplitude and effect of the mutual coupling. Therefore, the output of array can be expressed as: where x(t) is an N × 1 received signal vector, s(t) is a K × 1 incident signal vector, n(t) is an N × 1 Gaussian white noise vector with zero mean and noise power σ 2 N , and A is the N × K array manifold matrix formed by the set of the K array manifold vectors a(θ k , ϕ k , γ k , η k ) which can be expressed as: where a S (θ k, ϕ k ) is the spatial steering vector and a P (θ k, ϕ k , γ k , η k ) is the polarization-spatial domain steering vector. The expressions of the two vectors are shown as follows: a S (θ k, ϕ k ) = e −j2π(r sin ϕ k cos(−θ k ))/λ k , · · · , e −j2π(r sin ϕ k cos( 2π(m−1) M −θ k ))/λ k , · · · , e −j2π(r sin ϕ k cos( 2π(M−1) According to x(t), the observed output of identical polarization directions can then be obtained as: x n(t) and n y (t) = I M ⊗ e T y n(t), so that A x and A y are respectively formed by the set of the K array steering vectors a x (θ k , ϕ k , γ k , η k ) and a y (θ k , ϕ k , γ k , η k ), which can be expressed as: where P x (θ k , ϕ k , γ k , η k ) and P y (θ k , ϕ k , γ k , η k ) are scalars and denote the two components of polarization-spatial domain steering vector with the polarization direction along the x-axis and y-axis, respectively. For any incident signals, P x (θ k , ϕ k , γ k , η k ) and P y (θ k , ϕ k , γ k , η k ) are two different complex numbers. The module value of them can be regarded as gains of two dipoles generated by polarization reception, respectively, and the phase angles of them can be regarded as additional phase differences of two dipoles generated by polarization receiving, respectively. Then, a renewed received signal vector can be obtained as: where A = e x ⊗ A x (t) + e y ⊗ A y (t), n(t) = e x ⊗ n x (t) + e y ⊗ n y (t). A is formed by the set of the K array steering vectors a k (θ k , ϕ k , γ k , η k ), which can be expressed as:

Virtual Array Extension of the Uniform Circular EVSA Consisting of Orthogonal Dipoles
Assuming that the signal is a zero-mean, non-Gaussian, stationary random process, the fourth-order cumulant can be defined as: The results of Equation (10) can be collected in the form of matrixand denoted by the cumulant matrix R 4 as: where R 4 is the fourth-order cumulant matrix. If the incident signals are independent of each other, R 4 can be written as: where C S = diag(γ 4s k ), and the steering vector matrix B(θ, ϕ, γ, η) is expressed in the following form: Based on the equations above, it can be inferred that the fourth-order cumulant matrix of the received data can be considered as the covariance matrix of the data received by a virtual extended array, and the corresponding array steering vector of the kth signal can be expressed as follows: The coefficients of the virtual array steering vector can be regarded as the gains and additional phase differences of the elements in the virtual extended array. By separating the elements of the array steering vector according to different coefficient from vector b(θ k , ϕ k , γ k , η k ), four vectors can be obtained as follows: where P 2 x (θ, ϕ, γ, η), P x (θ, ϕ, γ, η)P y (θ, ϕ, γ, η), P y (θ, ϕ, γ, η)P x (θ, ϕ, γ, η) and P 2 y (θ, ϕ, γ, η) are the coefficients of the virtual array steering vector, respectively.
Separating the rows of the fourth-order cumulant matrix R 4 , which follows the same order of Equations (15)-(18), we obtain: Then, a new fourth-order cumulant matrix R 4 can be expressed as: where, c S (θ k , ϕ k ) = a S (θ k , ϕ k ) ⊗ a S (θ k , ϕ k ) is the spatial steering vector of the virtual array and then the polarization-spatial domain steering vector c P (θ k , ϕ k , γ k , η k ) can be written as: sin 2 θ k − sin θ k cos θ k cos ϕ k − sin θ k cos θ k cos ϕ k cos 2 θ k cos 2 ϕ k − sin θ k cos θ k cos 2 θ k cos ϕ k − sin 2 θ k cos ϕ k sin θ k cos θ k cos 2 ϕ k − sin θ k cos θ k cos 2 θ k cos ϕ k − sin 2 θ k cos ϕ k sin θ k cos θ k cos 2 ϕ k cos 2 θ k sin θ k cos θ k cos ϕ k sin θ k cos θ k cos ϕ k According to the Equations (25) and (26), it can be observed that R 4xy and R 4yx have identical expressions, thus both of them can be regarded as the output of the same dipole. Therefore, each element of the virtual array consists of three spatially co-located dipoles, and the polarization-spatial domain steering vector of the virtual array has three components (P 2 x (θ, ϕ, γ, η), P y (θ, ϕ, γ, η)P x (θ, ϕ, γ, η) and P 2 y (θ, ϕ, γ, η)) along the corresponding polarization direction. The spatial steering vector of the virtual array c S (θ k , ϕ k ) can be rewritten in the following matrix form: When the array element number is even (M = 2n), the element number of the virtual extended array is 2n 2 + 1, where each element consists of three co-located dipoles. The virtual array consists of one element located at the origin of the coordinate system and n concentric uniform circular arrays with the same element number of 2n. The elements of virtual array corresponding to the elements on the main diagonal of matrix D, constitute the largest uniform circular array of radius 2r, and the elements located on the i-th and M − i-th diagonal which are parallel to the main diagonal of matrix D, constitute a uniform circular array of radius 2r|cos(πi/M)|. The elements located on the n + 1-th diagonal of matrix D correspond to the virtual elements located at the origin of the coordinate system. In order to express this situation more clearly, an example of a uniform circular array composed of 6 elements is illustrated in Figure 3, where the matrix D 6 (left) and the corresponding virtual array (right) are shown.   Figure 3. Matrix D6 and the corresponding virtual array.
As before, the corresponding elements of D6 and the elements of the virtual array have the same color. The virtual array is comprised of one element located at the origin of the coordinate system and three concentric uniform circular arrays. Each uniform circular array has six elements.
The radii of the uniform circular arrays are 0 2 R r = , 1 2 cos(30 ) It can be inferred that the cross-dipole uniform circular array can be extended to a set of uniform circular arrays (or a set of uniform circular arrays with one additional element located at the origin of the coordinate system), whose elements consist of three co-located dipoles, by constructing the fourth-order cumulant of the received data. Thus, some pairs of sub-arrays with  n + , where each element consists of three co-located dipoles. The virtual array consists of one element located at the origin of the coordinate system and n concentric uniform circular arrays with the same element number of 2n . The elements of virtual array corresponding to the elements on the main diagonal of matrix D , constitute the largest uniform circular array of radius 2 r , and the elements located on the i -th and M i − -th diagonal which are parallel to the main diagonal of matrix D , constitute a uniform circular array of radius 2 cos( / ) r i M π . The elements located on the 1 n + -th diagonal of matrix D correspond to the virtual elements located at the origin of the coordinate system. In order to express this situation more clearly, an example of a uniform circular array composed of 6 elements is illustrated in Figure 3, where the matrix 6 D (left) and the corresponding virtual array (right) are shown. 6 (1,1)  Figure 3. Matrix D6 and the corresponding virtual array.
As before, the corresponding elements of D6 and the elements of the virtual array have the same color. The virtual array is comprised of one element located at the origin of the coordinate system and three concentric uniform circular arrays. Each uniform circular array has six elements.
The radii of the uniform circular arrays are 0 2 R r = , 1 2 cos(30 ) It can be inferred that the cross-dipole uniform circular array can be extended to a set of uniform circular arrays (or a set of uniform circular arrays with one additional element located at the origin of the coordinate system), whose elements consist of three co-located dipoles, by constructing the fourth-order cumulant of the received data. Thus, some pairs of sub-arrays with As before, the corresponding elements of D 6 and the elements of the virtual array have the same color. The virtual array is comprised of one element located at the origin of the coordinate system and three concentric uniform circular arrays. Each uniform circular array has six elements. The radii of the uniform circular arrays are R 0 = 2r, R 1 = 2r|cos(30 • )| and R 2 = r, respectively.
It can be inferred that the cross-dipole uniform circular array can be extended to a set of uniform circular arrays (or a set of uniform circular arrays with one additional element located at the origin of the coordinate system), whose elements consist of three co-located dipoles, by constructing the fourth-order cumulant of the received data. Thus, some pairs of sub-arrays with rotational-invariance can be obtained, and then the DOA and polarization of the incident signals can be jointly estimated based on the ESPRIT algorithm.

Proposed Algorithm
For the purposes of discussion, a uniform circular array equipped with eight elements is analyzed as an example in this section. However, the concepts presented in this section are also applicable to uniform circular array with arbitrary element number. The corresponding virtual extended array is drawn as a diagram as shown in Figure 4.

Proposed Algorithm
For the purposes of discussion, a uniform circular array equipped with eight elements is analyzed as an example in this section. However, the concepts presented in this section are also applicable to a uniform circular array with arbitrary element number. The corresponding virtual extended array is drawn as a diagram as shown in Figure 4. In Figure 4, the black elements represent the actual array and the white ones represent the virtual arrays derived from the actual array. The element located at the origin of the Cartesian coordinate system is labeled as A1, and the other elements located on four concentric uniform circular paths are labeled as Bi , Ci

Selection of Rotational Invariant Sub-Array Pairs
In order to avoid peak searching and reduce the computational complexity, we develop an ESPRIT algorithm based on the fourth-order cumulant to exploit the property of rotational invariance. It is therefore necessary to search for those pairs of sub-arrays that are provided with rotational invariance. In order to estimate the two-dimensional DOAs, at least two pairs of sub-arrays are required. The idea behind searching for sub-array pairs available for DOA and polarization estimation is proposed as follows.
First, the spatial phase factors between the virtual elements A1and B1, A1 and B3 are defined as the rotation invariant factors, which can be expressed as:  In Figure 4, the black elements represent the actual array and the white ones represent the virtual arrays derived from the actual array. The element located at the origin of the Cartesian coordinate system is labeled as A1, and the other elements located on four concentric uniform circular paths are labeled as Bi, Ci, Di and Ei (i = 1, 2, · · · , 8) , where the array elements B1, C1, D1 and E1 correspond to the elements D(1, 4), D(1, 3), D(1, 2) and D(1, 1) of the matrix D, respectively.

Selection of Rotational Invariant Sub-Array Pairs
In order to avoid peak searching and reduce the computational complexity, we develop an ESPRIT algorithm based on the fourth-order cumulant to exploit the property of rotational invariance. It is therefore necessary to search for those pairs of sub-arrays that are provided with rotational invariance. In order to estimate the two-dimensional DOAs, at least two pairs of sub-arrays are required. The idea behind searching for sub-array pairs available for DOA and polarization estimation is proposed as follows.
First, the spatial phase factors between the virtual elements A1 and B1, A1 and B3 are defined as the rotation invariant factors, which can be expressed as: Based on these two rotation invariant factors, two pairs of sub-arrays with the property of rotational invariance can be obtained as In order to display the pairs of sub-arrays clearly, the four sub-arrays (Sub-array 1, 2, 3 and 4) are drawn in four different colors, as shown in Figure 5. It is obvious that Sub-array 1 and Sub-array 2 have the property of rotational invariance, and u B1 is the rotation invariant factor. It's the same case for Sub-array 3 and Sub-array 4, where u B3 is the rotation invariant factor.
In addition, any two elements of Bi separated by one array element can be selected as the spatial phase factor, and the two corresponding pairs of sub-arrays then can be used to realize two-dimensional DOA estimation.

Two-Dimensional DOA Estimation
The relationships between the elements of the virtual sub-arrays and the elements of matrix D are listed in Table 1. It can be observed that all of the elements in the virtual array are contained within the two pairs of sub-arrays shown above. As each element consists of three virtual dipoles, the algorithm efficiently takes advantages of all non-redundant spatial information of the fourth-order cumulant matrices. The DOA estimation can be accomplished by using one of the three groups of dipoles, which has the same gain and additional phase difference. In general, to improve the estimation accuracy, the algorithm usually uses more dipoles consequently resulting in higher computation complexity. In order to reduce the computational complexity effectively, dipoles of one group are only used in the proposed algorithm while ensuring the estimation accuracy.
In addition, any two elements of Bi separated by one array element can be selected as the spatial phase factor, and the two corresponding pairs of sub-arrays then can be used to realize two-dimensional DOA estimation.

Two-Dimensional DOA Estimation
The relationships between the elements of the virtual sub-arrays and the elements of matrix D are listed in Table 1. Based on the equation D(i, j) = c S ((i − 1)M + j), and the corresponding relationships between the elements of the virtual extended array and the rows of matrix R 4 , the corresponding selection matrices E 1 , E 2 , E 3 and E 4 can be constructed as follows: where the matrix element corresponding to the n-th element of the m-th sub-array appears as the i submn -th row and j submn -th column of matrix D. N denotes the total number of sub-arrays, which is 14 for the eight-sensor array. e (i submn −1)M+j submn denotes a M 2 × 1 vector with its ((i submn − 1)M + j submn )th element as the only non-zero item. Therefore, we can get: where B 1 (θ, ϕ, γ, η), B 2 (θ, ϕ, γ, η), B 3 (θ, ϕ, γ, η) and B 4 (θ, ϕ, γ, η) denote the steering vectors of sub-array 1, sub-array 2, sub-array 3 and sub-array 4, respectively. The relationships between B 1 (θ, ϕ, γ, η), B 2 (θ, ϕ, γ, η), B 3 (θ, ϕ, γ, η) and B 4 (θ, ϕ, γ, η) can be expressed as follows: where . Furthermore, we have: Then, matrices C 1 and C 2 are constructed based on Equations (40) and (41) as: When using singular value decomposition (SVD), C 1 and C 2 can be denoted as: The subspace spanned by array steering vector is included in the signal subspace, hence there exists a nonsingular matrix T expressed as: Consequently, we have: so, the total least square solutions of Equations (48) and (49) are: θ, ϕ, γ, η) and B 3 (θ, ϕ, γ, η) are of full rank, thus: This implies that the diagonal elements of the diagonal matrices Φ i can be estimated by obtaining the K eigenvalues of each matrix Ψ i , where the corresponding eigenvectors are the column vectors of T i .

Polarization Estimation
It can be inferred from Equations (15) and (16) that: The rotation invariant matrix can then be constructed by the above equation as: The matrix C 3 can be constructed using Equations (19) and (20) as: By using SVD, C 3 can be represented as: Based on the relationship between the subspace spanned by array steering vector and the signal subspace, we have: The total least square solution of Equation (57) is: In the same way, Φ 3 can be calculated by eigen-decomposition of Ψ 3 .

Pair Matching
According to the discussion above, it is observed that the three eigen-decompositions of Ψ i are independent of each other. Therefore, pair matching is required when multiple signals impinge on the array. The matching process can be understood as the problem that whether any two parameters are successfully matched. Here we firstly consider the matching between Φ 1 and Φ 2 . Parameters of one signal correspond to the same signal subspace, and the eigen-vector corresponding to different eigenvalues are mutually orthogonal. Therefore, a matrix can be constructed for ranking based on the following principle: where T 1 (i) denotes the i-th column of T 1 . Therefore: The pair matching between Φ 1 and Φ 3 is also similarly conducted. At the same time, K groups of values u B1k , u B3k and (P x /P y ) k are obtained , where k = 1, 2, · · · , K. The values of the azimuth angles, elevation angles, polarization auxiliary angles and polarization phase differences can be estimated as follows: γ k = arctan sin θ k + angle((P x /P y ) k ) cos θ k cos ϕ k cos θ k − angle((P x /P y ) k ) cos ϕ k sin θ k (63)

Steps of the Proposed Algorithm
The steps of the proposed algorithm are executed as follows: Step 1: Construct the fourth-order cumulant based on Equations (19) and (20).
Step 2: Select two pairs of sub-arrays which display rotational invariance based on the theory introduced in Section 4.1.
Step 3: Construct the rotational invariance matrices C 1 , C 2 and C 3 .
Step 5: Perform pair matching among Ψ 1 , Ψ 2 and Ψ 3 and estimate the DOA and polarization information of the incident signals.
The main contributors of computational complexity of LV-MUSIC method are discussed as follows: (1) The structure of the covariance matrix requires C 1 = LN 2 = 4LM 2 complex multiplications; (2) The EVD of the covariance matrix requires C 2 = 4M 2 (K + 2) complex multiplications; (3) Assuming that peak searching with all identical parameters, the construction of the spectrum function requires C 3 = (1 + 360/n) 2 (1 + 90/n) 2 (4M 2 + 2M) complex multiplications, where each searching point of the spectrum function requires C 31 = N 2 + N = 4M 2 + 2M complex multiplications. 1 + 360/n is the number of searching points along the azimuth angle and polarization phase difference, 1 + 90/n is the number of searching points along the elevation angle and polarization auxiliary angle, and n is the searching step with respect to every parameter.
When the array structure, source number and snapshots are definite, it can be inferred from the analysis above that the computational complexity of the proposed algorithm is a fixed number, but that of the LV-MUSIC algorithm changes with searching steps. In this paper, we assume that the source number is 2 and the number of snapshots is 100. The computational complexity of the proposed algorithm and that of the LV-MUSIC algorithm using different searching steps are shown in Figure 6. that of the LV-MUSIC algorithm changes with searching steps. In this paper, we assume that the source number is 2 and the number of snapshots is 100. The computational complexity of the proposed algorithm and that of the LV-MUSIC algorithm using different searching steps are shown in Figure 6. From Figure 6, it can be observed that the computational complexity of the LV-MUSIC algorithm is higher than that of the proposed algorithm when the searching step is less than 60°. The estimation accuracy of the LV-MUSIC algorithm is known to converge to half of the searching steps, and therefore the LV-MUSIC algorithm becomes invalid when the searching step is too large. In other words, though the LV-MUSIC algorithm can be used to estimate DOA and polarization jointly in theory, its computational complexity is not acceptable, mainly because the LV-MUSIC algorithm requires four dimensional peak searching. Fortunately, the proposed algorithm are not limited to these problems with a fixed computational complexity. The estimation accuracy and the performance of the proposed algorithm are discussed in the following section. From Figure 6, it can be observed that the computational complexity of the LV-MUSIC algorithm is higher than that of the proposed algorithm when the searching step is less than 60 • . The estimation accuracy of the LV-MUSIC algorithm is known to converge to half of the searching steps, and therefore the LV-MUSIC algorithm becomes invalid when the searching step is too large. In other words, though the LV-MUSIC algorithm can be used to estimate DOA and polarization jointly in theory, its computational complexity is not acceptable, mainly because the LV-MUSIC algorithm requires four dimensional peak searching. Fortunately, the proposed algorithm are not limited to these problems with a fixed computational complexity. The estimation accuracy and the performance of the proposed algorithm are discussed in the following section.

Simulation Results
In this section, several numerical simulation results are presented to illustrate the performance of the proposed algorithm. A uniform circular array with eight crossed dipoles elements and array radius of 0.087m is used in the simulations. Complex additive Gaussian white noise is added into the system and the number of Monte Carlo trials is 200.
The root mean square error (RMSE) for the estimated parameters is defined as: where M is the number of Monte Carlo trials, χ i is the perfect value of the estimated parameter corresponding to the i-th incident signal, andχ im is the estimated value of the m-th time Monte Carlo trials.  Figure 7 shows that the points which represent the results of the proposed algorithm always cluster around the true value of the estimated parameters. From Figure 7a, it can be seen that the maximum evaluated error of elevation and azimuth angles are both less than 2 • . From Figure 7b, we can see that the maximum evaluated error of polarization auxiliary angles and polarization phase differences are less than 2 • and 3 • , respectively. It can be observed that the proposed algorithm shows excellent performance.  Figure 7 shows that the points which represent the results of the proposed algorithm always cluster around the true value of the estimated parameters. From Figure 7a, it can be seen that the maximum evaluated error of elevation and azimuth angles are both less than 2°. From Figure 7b, we can see that the maximum evaluated error of polarization auxiliary angles and polarization phase differences are less than 2°and 3°, respectively. It can be observed that the proposed algorithm shows excellent performance.

Performance under Different SNR
In order to demonstrate the remarkable performance of the proposed algorithm, we compare the performance of the proposed algorithm with that of the LV-MUSIC algorithm and the Cramer-Rao bound (CRB) with different values of SNR.As shown in Figure 6 above, the computational complexity of LV-MUSIC algorithm is huge when DOA and polarization parameters

Performance under Different SNR
In order to demonstrate the remarkable performance of the proposed algorithm, we compare the performance of the proposed algorithm with that of the LV-MUSIC algorithm and the Cramer-Rao bound (CRB) with different values of SNR.As shown in Figure 6 above, the computational complexity of LV-MUSIC algorithm is huge when DOA and polarization parameters are estimated within acceptable searching steps. Due to limited computing power, when we focus on the DOA parameters estimation performance of LV-MUSIC algorithm, we suppose that the polarization parameters are known or they are estimated beforehand, and when we focus on the polarization parameters estimation performance of LV-MUSIC algorithm, the DOA parameters are set to be known or have been estimated already. Referring to this comparison, we assume that there are two signals with parameters (θ, ϕ, γ, η) of (60 • , 40 • , 30 • , 40 • ) and (50 • , 50 • , 40 • , 60 • ), the number of snapshots is 200, and the SNR changes from −5 dB to 40 dB. As shown in Figures 8 and 9, it can be seen that as the SNR increases, the RMSE of the estimated parameters decreases gradually for the proposed algorithm, the LV-MUSIC algorithm and the corresponding CRB. Meanwhile, the proposed algorithm and the LV-MUSIC algorithm both achieve good DOA and polarization estimation performance, and the results also match with the CRB perfectly. The estimation accuracy of the proposed algorithm is similar to the LV-MUSIC algorithm. When the SNR increases, the RMSE of the proposed algorithm continuously decreases, and yet that of the LV-MUSIC algorithm remains unchanged, which depends on the searching step. Comparison between Figures 8 and 9 show that the RMSE of polarization parameters estimation is higher than that of DOA parameters estimation for the proposed algorithm for the same SNR. The reason is that the estimation of polarization parameters is carried out based on the estimated values of DOA parameters. The secondary error, thus, affects the estimation accuracy of polarization parameters. (a) (b) Figure 9. RMSE of (a) polarization auxiliary angle and (b) polarization phase difference estimation versus SNR (for 200 snapshots).

Performance for Different Numbers of Snapshots
The performance of the proposed algorithm is compared with the LV-MUSIC algorithm and the Cramer-Rao bound (CRB) versus different number of snapshots. Two signals with parameters ( , , , ) θ ϕ γ η of (60 , 40   Comparison between Figures 8 and 9 show that the RMSE of polarization parameters estimation is higher than that of DOA parameters estimation for the proposed algorithm for the same SNR. The reason is that the estimation of polarization parameters is carried out based on the estimated values of DOA parameters. The secondary error, thus, affects the estimation accuracy of polarization parameters. (a) (b) Figure 9. RMSE of (a) polarization auxiliary angle and (b) polarization phase difference estimation versus SNR (for 200 snapshots).

Performance for Different Numbers of Snapshots
The performance of the proposed algorithm is compared with the LV-MUSIC algorithm and the Cramer-Rao bound (CRB) versus different number of snapshots. Two signals with parameters ( , , , ) θ ϕ γ η of (60 , 40  are considered here, the SNR is set as 20 dB, and the number of snapshots changes from 10 to 2000. From Figures 10 and 11, we can see that the Comparison between Figures 8 and 9 show that the RMSE of polarization parameters estimation is higher than that of DOA parameters estimation for the proposed algorithm for the same SNR. The reason is that the estimation of polarization parameters is carried out based on the estimated values of DOA parameters. The secondary error, thus, affects the estimation accuracy of polarization parameters.

Performance for Different Numbers of Snapshots
The performance of the proposed algorithm is compared with the LV-MUSIC algorithm and the Cramer-Rao bound (CRB) versus different number of snapshots. Two signals with parameters (θ, ϕ, γ, η) of (60 • , 40 • , 30 • , 40 • ) and (50 • , 50 • , 40 • , 60 • ) are considered here, the SNR is set as 20 dB, and the number of snapshots changes from 10 to 2000. From Figures 10 and 11, we can see that the RMSE values of these estimated parameters decrease gradually for the proposed algorithm, the LV-MUSIC algorithm and the CRB as the number of snapshots increases. The RMSE of the proposed algorithm is larger than that of the LV-MUSIC algorithm, and the two algorithms both show good DOA and polarization estimation performance and keep agreement with the CRB perfectly. Obviously, reducing computational complexity by using the proposed algorithm is at the expense of slightly worse performance as shown in the figure. (a) (b) Figure 11. RMSE of (a) polarization auxiliary angle and (b) polarization phase difference estimation versus number of snapshots (for SNR = 20 dB).

Running Time
The  (a) (b) Figure 11. RMSE of (a) polarization auxiliary angle and (b) polarization phase difference estimation versus number of snapshots (for SNR = 20 dB).

Running Time
The

Running Time
The running time of the proposed algorithm is compared against that of the LV-MUSIC algorithm with two dimensional search. Two signals with parameters (θ, ϕ, γ, η) of (60 • , 40 • , 30 • , 40 • ) and (50 • , 50 • , 40 • , 60 • ) are considered here, the SNR is set as 20dBand the number of snapshots is 200. For the LV-MUSIC algorithm, the azimuth angle and the elevation angle have been searched in the range of 0 • to 360 • and 0 • to 90 • with an interval of 0.5 • , respectively. 200 Monte Carlo simulations are carried out. Table 2 shows the average running time of these two algorithms. The running time of the proposed algorithm is shorter than that of the LV-MUSIC algorithm with two dimensional searching. In addition, four dimensional search is required when the LV-MUSIC algorithm is used to jointly estimate of DOA and polarization information. It is assumed that the polarization auxiliary angle and the polarization phase difference of the LV-MUSIC algorithm have been achieved by searching in the range of 0 • to 360 • and 0 • to 90 • with an interval of 0.5 • , respectively. The running time of LV-MUSIC algorithm with four dimensional search is (360 • /0.5 • )(90 • /0.5 • ) = 129, 600 times of the case of two dimensional search. Therefore, the running time of the LV-MUSIC algorithm with four dimensional search is unacceptable, and the proposed algorithm can obtain the four parameters quite efficiently.

Conclusions
A fourth-order cumulant-based ESPRIT algorithm is proposed in this paper, which can achieve the joint estimation of the DOA and polarization information based on a uniform circular EVSA. The proposed algorithm overcomes the limitation of the ESPRIT algorithm of failing in uniform circular EVSAs, and simultaneously achieves a significant reduction in terms of the computation cost. The fourth-order cumulant has been used to virtually extend the original array, and then construct a few sub-arrays with identical shapes. By matching two pairs of sub-arrays with the rotation-invariant structure, the estimation of DOA and polarization information can be carried out using the ESPRIT algorithm. As the algorithm does not require the construction of the spectrum function and does not resort to a multidimensional peak search, the estimation results can be achieved quite efficiently and also ensured the acceptable simulation accuracy. The reduction in computational complexity of the proposed algorithm has been illustrated by comparing against the LV-MUSIC algorithm for different searching steps theoretically and numerically. Numerical simulation results validate that the proposed algorithm has higher calculation efficiency than the LV-MUSIC algorithm. Future work may focus on utilizing a hypercomplex framework, such as a tensor, to re-establish the model of a four-order cumulant matrix, aiming at obtaining higher estimation accuracy in the direction finding algorithm.