Real-Valued 2D MUSIC Algorithm Based on Modified Forward/Backward Averaging Using an Arbitrary Centrosymmetric Polarization Sensitive Array

Two-dimensional multiple signal classification (MUSIC) algorithm based on polarization sensitive array (PSA) has excellent performance. However, it suffers a high computational complexity due to a multitude of complex operations. In this paper, we propose a real-valued two-dimensional MUSIC algorithm based on conjugate centrosymmetric signal model, which is applicable to arbitrary centrosymmetric polarization sensitive array. The modified forward/backward averaging, which can be applied to the PSA, is presented. Hence, the eigen-decomposition analysis process and spectrum function computation are converted into real domain, prominently reducing the computational complexity. Then, the direction-of-arrival (DOA) estimation is decoupled from the polarization parameter estimation so that the four-dimensional spectral peak search process is avoided. The theoretical computational complexity is discussed and the Cramer-Rao bound (CRB) of DOA estimation is derived in this paper. The simulation results indicate that the proposed algorithm achieves superior accuracy in DOA estimation and has low computational complexity.


Introduction
Parameter estimation is an important area in array signal processing, such as adaptive beamforming [1,2] and direction-of-arrival (DOA) estimation [3][4][5][6][7]. Especially, the two-dimensional (2D) DOA estimation based on polarization sensitive array (PSA) shows great potential [8][9][10], which has gained considerable interests. Compared with traditional scalar array (TSA), the PSA composed of electromagnetic vector sensors (EMVSs) has more inherent merits. Besides spatial domain information, polarization domain information of the electromagnetic signals can also be measured by the PSA [11,12], which can bring greater potential capacity for parameter estimation.
During the past decades, many effective algorithms based on TSA have been extended to PSA, such as multiple signal classification (MUSIC) [13,14], estimation of signal parameters via rotational invariance techniques (ESPRIT) [15][16][17], Root-MUSIC [18][19][20], etc. With the development of research, some specific joint DOA and polarization parameter estimation algorithms for PSA have been presented, including vector cross-product based algorithms [21,22] and hyper complex MUSIC algorithms [23,24]. In almost all of the literature considering PSAs, the data model makes use of the concatenated vector of the received signal data with multi-components, resulting in the so-called "long-vector" approach [24]. An enhanced ESPRIT algorithm which improves estimation accuracy memorably and can solve the ambiguity of DOA estimation was proposed in [22]. However, this where φ k ∈ [0, π] and θ k ∈ [0, 2π) denote the elevation angle measured from the vertical z-axis and the azimuth angle measured from the positive x-axis, respectively, see in Figure 1; γ k ∈ [0, π/2] and η k ∈ [−π, π) denote the auxiliary polarization angle and the polarization phase difference, respectively. identity matrix, and denotes an × anti-identity matrix with ones on its anti-diagonal and zeros elsewhere.

Received Signal Model
Assume that K uncorrelated transverse electromagnetic plane waves, which have traveled through a homogeneous isotropic medium, impinge on the electromagnetic vector sensors. Then, the k-th unit-power completely polarized incident source has the electric-field vector and the magnetic-filed vector , which can expressed in Cartesian coordinates as [31][32][33] ( Assume that a centrosymmetric polarization sensitive array consists of M ECCD pairs is arranged in the spatial plane of the coordinate system which includes x-axis, y-axis and the original point. For brevity, we abbreviate this plane as xoy-plane. The two vibrators of each ECCD pair are deployed parallel to x-axis and y-axis, respectively. Since ECCD pairs only receive components , and , of electric-field vector, (1) is simplified as ( Note that , and , depend only on the source spatial angular locations and the incent signal polarization states, respectively.
The spatial steering vector of the k-th narrow-band incent signal equals where Δ is the inter-vector-sensor spatial phase factor relating the k-th source to the m-th vector sensor at location ( , , ) and is defined as Assume that a centrosymmetric polarization sensitive array consists of M ECCD pairs is arranged in the spatial plane of the coordinate system which includes x-axis, y-axis and the original point. For brevity, we abbreviate this plane as xoy-plane. The two vibrators of each ECCD pair are deployed parallel to x-axis and y-axis, respectively. Since ECCD pairs only receive components e x,k and e y,k of electric-field vector, (1) is simplified as ( Note that v p,k and h r,k depend only on the source spatial angular locations and the incent signal polarization states, respectively. The spatial steering vector of the k-th narrow-band incent signal equals where ∆ km is the inter-vector-sensor spatial phase factor relating the k-th source to the m-th vector sensor at location (x m , y m , z m ) and is defined as where λ means wavelength. For simplicity, we use a p,k and a s,k to respectively represent a p (θ k , φ k , γ k , η k ) and a s (θ k , φ k ) in the following parts. As such, the received signal model (RSM) of array can be expressed as a long-vector form where n m (t) symbolizes the 2 × 1 zero-mean additive complex white Gaussian noise vector at the m-th vector sensor.

Conjugate Centrosymmetric Signal Model
Suppose that there a 2M × 2M permute matrix where P e is an M × M permute matrix, and it can be obtained readily that P e P T e = P T e P e = I 2M . Then, a novel signal model is proposed as follow y(t) = P e x(t) = P e (As(t) + N(t)) = A e s(t) + n e (t), A se = P e A s = [ P e a s,1 , · · · , P e a s,K ] = [a se,1 , · · · , a se,K ]. (12) According to (3) and (4), it is obvious that the inter-vector-sensor spatial phase factors a s,k of two electromagnetic vector-sensors which are centrosymmetric around the origin of coordinate system are conjugate. Thus, for a centrosymmetric array, there must be a calculable and unique P e such that the spatial steering vector of the proposed signal model can be expressed as where 1 K denotes a K × 1 matrix of ones, A . Furthermore, note that the matrix multiplications involving P e or J are not counted as floating operations because they merely represent the permutations of the rows or columns of the matrix being multiplied [34]. Without loss of generality, consider a uniform circular array (UCA) with 12 vector sensors and a uniform rectangular array (URA) with 3 rows and 4 columns which are shown in Figure 2. Consequently the permute matrix P e corresponding to each array is expressed as where 0 6 denotes a 6 × 6 null matrix.
where and denote the eigenvector matrix and eigenvalue diagonal matrix, subscript x and y denote the RSM and CCSM, subscript and denote the signal-and noise-subspace, respectively.
According to (10) and (15), the covariance matrix of CCSM can be also expressed as Compared the expressions of in (16) and (17), we have Note that the elements in , and , are just a displacement of the elements in , and , due to the permuted matrix . Based on the MUSIC algorithm, the signal subspace is orthogonal to the noise subspace. The signal subspace is spanned by the steering vectors, which means that the steering vectors of sources are also orthogonal to the noise subspace, i.e., , = , = 1, ⋯ , .
As such, we can obtain which illustrates that the noise-subspace and the steering vectors of CCSM are still orthogonal. □

The Modified FB Averaging
Considering about finite received array data, the covariance matrix of CCSM is expressed as Proof. The covariance matrices of RSM and CCSM are where U and Λ denote the eigenvector matrix and eigenvalue diagonal matrix, subscript x and y denote the RSM and CCSM, subscript S and N denote the signal-and noise-subspace, respectively. According to (10) and (15), the covariance matrix of CCSM can be also expressed as Compared the expressions of R y in (16) and (17), we have Note that the elements in U y,S and U y,N are just a displacement of the elements in U x,S and U x,N due to the permuted matrix P e . Based on the MUSIC algorithm, the signal subspace is orthogonal to the noise subspace. The signal subspace is spanned by the steering vectors, which means that the steering vectors of sources are also orthogonal to the noise subspace, i.e., As such, we can obtain which illustrates that the noise-subspace and the steering vectors of CCSM are still orthogonal.

The Modified FB Averaging
Considering about finite received array data, the covariance matrix of CCSM is expressed aŝ where L is the number of snapshots and Y = 1/ √ L[y 1 (t), · · · , y L (t)] is the CCSM observation matrix. In contrast to [34,35], define the modified FB covariance matrix aŝ It is noteworthy that, different from the original FB averaging [34], the modified FB averaging uses the unitary permute matrix J 2M conforming to (23) instead of an anti-identity matrix. This significant change avoids the confusion of the different components of the electric-field vector of the incident signals, and allows the modified FB averaging to be applied to PSAs consist of ECCD pairs. In addition, it is straightforward to certify thatR Then, let Q be a unitary matrix which is defined as where 0 (M−1)×2 denotes an (M − 1) × 2 null matrix, and J M−1 is defined similarly as (23). It is obvious that Q = J 2M Q * . Therefore, we can obtain the real-valued covariance matrix as By substituting (22) into (26), R RV can be further simplified as The second equation in (27) is derived readily by using (24). (27) indicates that R RV only depends on the covariance of CCSM.
The real-valued covariance matrix R RV is also symmetric due to the following properties Thus, the EVD of R RV require only real-valued computation [36]. Furthermore, the eigenvalues and eigenvectors of R RV are real-valued as well. The EVDs ofR FB and R RV arê where M = diag[µ 1 , µ 2 , · · · , µ M ] and Π = diag[π 1 , π 2 , · · · , π M ] are eigenvalue diagonal matrices and E = [e 1 , e 2 , · · · , e M ] and V = [v 1 , v 2 , · · · , v M ] are eigenvector matrices. The subscript S and N represent signal-and noise-subspace, respectively. Substituting (30) into (26),R FB can be expressed aŝ Compared (31) with (29), we have the following equations Define Z FB as where P is defined as Then, we can simplify R RV by using (26), (33) and (34) Note that Z FB is a real-valued matrix.
for odd M. Then, substituting Q H , P, and Y FB yields It is easy to determine that Z FB is always real-valued. According to the previous conclusion, the matrix multiplication involving permute matrix does not contain any floating operations, but only exchanges the corresponding rows or columns. Thus, from (36), we can deduce that calculating Z FB does not involve general matrix multiplication, but is a simple construction based on Y (the CCSM data matrix) involving only additions. Therefore, it has significant consequences for processing modified FB covariance matrix.

2D DOA Estimation with Real-Valued Spectrum Function
Applying MUSIC to the modified FB covariance matrix, the joint spatial-polarized spectrum function is In contrast with the joint spectrum function of the customary sensor-vector MUSIC algorithm, the signal steering vectors of CCSM and the noise-subspace of the modified FB covariance matrix are used in (37). Note that calculating the noise-subspace of modified FB covariance matrix is almost a real operation, thus reducing the substantial computational amount. However, f (θ, φ, γ, η) is a four-dimensional function with respect to both 2D DOA parameters and 2D polarization parameters. Hence, the process of parameter estimation is actually a four-dimensional iterative peak search process which leads to large computation cost. If the four-dimensional joint spectrum function can be turned into a two-dimensional function only with respect to the 2D DOA parameters, the four-dimensional iterative search will be transformed to a two-dimensional iterative search, namely, the DOA parameters can be estimated without polarization parameters estimation.
Substituting (7), (11) and (12) into (37), the reciprocal of the joint spectrum function equals where only depends on DOA parameters, while h r only depends on polarization parameters. Therefore, we can obtain that where (θ k , φ k , γ k , η k ) are the DOA and polarization states of the k-th source, respectively. According to the MUSIC algorithm, it is known that Furthermore, by using (39) and (41), we have It is obvious that the solution to (40) and (42) is unique, which indicates that H(θ, φ) is rank-deficient only when (θ, φ) equals the actual DOA, i.e., (θ, φ) ∈ {(θ k , φ k )} K k=1 . Thus, the DOA can be estimated by solving the following optimization problem Therefore, the direction parameter estimation is successfully decoupled from the polarization estimation, thus reducing the number of variables that need to be optimized.
Then, exchange the expression of H(θ, φ) in order to estimate the DOAs only through real-valued operations. Substituting (13) and (32) into (39), it can be obtained that where Obviously, H(θ, φ) is a real-valued matrix since D(θ, φ) and V N is real-valued. From (39) and (44), it is clear that these two forms of H(θ, φ) are exactly identical. However, calculating H(θ, φ) according to (44) only involves the operations between real-valued matrices. Therefore, (44) is a much better choice to calculate the spectrum function and estimate the azimuth angle and the elevation angle.
The whole process of the proposed algorithm is summarized in Table 1. Table 1. Implementation steps of the entire approach.

Input The Received Signal Model (RSM) of Array x(t) = (A s A P )s(t) + n(t).
Output 2D DOA Estimation.

Step 2
Compute the real-valued matrix Z FB by (36).
Step 3 Compute the symmetric real-valued covariance matrix R RV = Z FB Z T FB .
Step 4 Perform the EVD of R RV by (30) to get V N .

Cramer-Rao Bound (CRB) Analysis
Cramer-Rao bound (CRB) is an important parameter in discussing the performance of a system [37]. In this subsection, the CRB of 2D DOA estimation is discussed. The array received signal model is shown in Equation (5).
The Fisher information matrix (FIM) with respect to azimuth angle θ = [θ 1 , · · · , θ K ], elevation angle φ = [φ 1 , · · · , φ K ], auxiliary polarization angle γ = [γ 1 , · · · , γ K ] and polarization phase difference η = [η 1 , · · · , η K ] is written as Here, the CRB of 2D DOA estimation is solely considered. Thus, the sub-matrix of FIM with respect to azimuth and elevation is expressed as According to literatures [38,39], the (i,j)-th element of F θθ is given by Hence, the sub-matrix of F DOA with respect to azimuth is Similarly, the sub-matrices with respect to elevation and the cross terms of F DOA can be expressed as where tr{·} denotes the trace of a matrix, e i denotes the i-th column of the identity matrix. R s is the covariance matrix of sources which is defined as are the derivations of A with respect to θ i and φ j with i, j = 1, 2, · · · , K. Then, the CRB matrix C DOA with respect to azimuth and elevation is the inverse of F DOA , i.e., As such, the CRBs of azimuth and elevation of the k-th source are where C DOA (k, k) denotes the (k,k)-th element of C DOA .

Theoretical Computational Complexity Analysis
To illustrate the computational effectiveness of the proposed algorithm, the computational complexity is explicitly discussed in this subsection comparing with LV-MUSIC algorithm [24]. The LV-MUSIC algorithm is an improved MUSIC algorithm making use of long-vector data form. It has the advantages of high resolution and superior accuracy, and can be applied to the PSAs with arbitrary structure composed of vector sensors. In view of the fairness of comparison, when considering of LV-MUSIC algorithm, assume that the polarization parameters are priori or have been estimated accurately, so that the search dimensions of the two algorithms are consistent. We abbreviate the LV-MUSIC algorithm with 2D DOAs peak searching as 2D LV-MUSIC temporarily.
The computational complexity of the proposed algorithm mainly includes the following three parts. (1) Computing the real-valued covariance matrix R RV from RSM requires about (16L − 4)M 2 + 6LM flops; (2) The real-valued EVD (based on symmetric QR algorithm [40]) of R RV requires about 72M 3 flops; (3) Computing the spectrum function in (44) for each searching point requires about 16M 3 + (28 − 8K)M 2 + 9M + 8 flops. Assume that there are N SP searching points in the interested search space. Thus, the computational complexity of the proposed algorithm is Similarly, the computational complexity of 2D LV-MUSIC algorithm mainly consists of following three parts.
(57) Figure 3 shows the complexity of the proposed algorithm and 2D LV-MUSIC algorithm versus the number of snapshots, the number of sensors and the search step of parameters. From Figure 3a-c, we can find that the theoretical computational cost of both the algorithms increases as the number of snapshot and the number of the sensors become larger, while reduces as the search step of the parameters decreases. In addition, it is obvious that the computational complexity of the proposed algorithm is always lower than 2D LV-MUSIC algorithm.
following three parts. (1) Computing complex covariance matrix from RSM requires about 8(3 2 + − 1) 2 flops; (2) The EVD of complex covariance matrix requires about 288 3 flops; (3) Computing the spectrum function for each searching point requires about 64 3 + (56 − 43 ) 2 + 88 + 21 flops. Under the assumption that searching points are used in the interested search space, we can obtain the computational complexity of 2D LV-MUSIC algorithm is C LV−MUSIC = 288 3 + 8(3 2 + − 1) 2 + (64 3 + (56 − 32 ) 2 + 88 + 21). (57) Figure 3 shows the complexity of the proposed algorithm and 2D LV-MUSIC algorithm versus the number of snapshots, the number of sensors and the search step of parameters. From Figure 3a-c, we can find that the theoretical computational cost of both the algorithms increases as the number of snapshot and the number of the sensors become larger, while reduces as the search step of the parameters decreases. In addition, it is obvious that the computational complexity of the proposed algorithm is always lower than 2D LV-MUSIC algorithm.

Simulation Results
In this section, several simulation results are presented to illustrate the performance of the proposed algorithm. A UCA and a URA as shown in Figure 2 are used in the following simulations. The radius of UCA is 2λ and the interspace of URA is λ. Two uncorrelated far-field narrowband signals whose DOAs and polarization parameters (θ, φ, γ, η) are 50 • , 10 • , 40 • , −120 • and 100 • , 20 • , 45 • , 70 • , respectively, are assumed to impinge on the array. The first source is left-circularly elliptically polarized, whereas the second source is right-circularly elliptically polarized. Complex additive Gaussian white noise is added into the system. The root mean square error (RMSE) of the k-th source is defined as where N is the number of independent Monte Carlo trails. ς k is the parameter of the k-th source, and ς k,n is the estimated values of ς k in the n-th trail.  Figure 4 shows the simulation results of the elevation and the azimuth of the proposed algorithm and 2D LV-MUSIC algorithm in 200 trials. The signal-noise ratio (SNR) is 10 dB, and the number of snapshots is 300. The search step is set as 0.1 • . The green asterisk symbols represent the actual DOA values, while the blue points represent the estimated DOA values. The small windows inside each figure are the enlarged version of the dotted areas. In Figure 4, the estimated DOA values always cluster around the actual DOA values, which demonstrate that both two algorithms derive a small estimation bias and variance. Note that in Figure 4, all the simulation results are distributed on the searching points, which is a characteristic of the spectral peak search-based algorithms. In addition, from Figure 4a, it can be seen that the maximum evaluated error of elevation and azimuth angles of proposed algorithm are less than 0.2 • and 0.5 • , respectively; from Figure 4b, we can see that the maximum evaluated error of elevation and azimuth angles of the 2D LV-MUSIC algorithm are less than 0.2 • and 0.9 • , respectively. That is, the scatter plots of proposed algorithm concentrate on a smaller area than 2D LV-MUSIC algorithm, which manifests that the proposed algorithm has a better estimation performance.

The Simulation Results Distribution Scatter Plots
where is the number of independent Monte Carlo trails. is the parameter of the k-th source, and ̂ , is the estimated values of in the n-th trail. Figure 4 shows the simulation results of the elevation and the azimuth of the proposed algorithm and 2D LV-MUSIC algorithm in 200 trials. The signal-noise ratio (SNR) is 10 dB, and the number of snapshots is 300. The search step is set as 0.1°. The green asterisk symbols represent the actual DOA values, while the blue points represent the estimated DOA values. The small windows inside each figure are the enlarged version of the dotted areas. In Figure 4, the estimated DOA values always cluster around the actual DOA values, which demonstrate that both two algorithms derive a small estimation bias and variance. Note that in Figure 4, all the simulation results are distributed on the searching points, which is a characteristic of the spectral peak search-based algorithms. In addition, from Figure 4a, it can be seen that the maximum evaluated error of elevation and azimuth angles of proposed algorithm are less than 0.2° and 0.5°, respectively; from Figure 4b, we can see that the maximum evaluated error of elevation and azimuth angles of the 2D LV-MUSIC algorithm are less than 0.2° and 0.9°, respectively. That is, the scatter plots of proposed algorithm concentrate on a smaller area than 2D LV-MUSIC algorithm, which manifests that the proposed algorithm has a better estimation performance.

The Estimation Performance versus SNR and the Number of Snapshots
In this subsection, the RMSEs of the proposed algorithm and 2D LV-MUSIC algorithm versus SNR and the number of snapshots are compared. The azimuth angle and elevation angle are

The Estimation Performance versus SNR and the Number of Snapshots
In this subsection, the RMSEs of the proposed algorithm and 2D LV-MUSIC algorithm versus SNR and the number of snapshots are compared. The azimuth angle and elevation angle are searched in the range of 0 • to 360 • and 0 • to 90 • with a step of 0.1 • , respectively. 500 Monte Carlo trails are performed to acquire the RMSE.
We first examine the RMSE versus SNR where the SNR of the signals varies from −10 dB to 30 dB with a step of 5 dB and the number of snapshots is set as 300. As shown in Figure 5, it is obvious that for the estimation of the azimuth, the proposed algorithm can achieve a lower RMSE than the 2D LV-MUSIC especially for the UCA. As to the estimation of the elevation, in low SNR region, the proposed algorithm using UCA outperforms the 2D LV-MUSIC. In high SNR region, the performance of the proposed algorithm and the 2D LV-MUSIC is similar for estimating the elevation. However, as analyzed in Section 4.2, the computational complexity of the proposed algorithm is much lower than the 2D LV-MUSIC. when the modified FB observation matrix is obtained by the CCSM observation matrix, the size of the data matrix yields an effective doubling, which achieves an expected improvement in the estimator variances. The improvement makes the proposed algorithm have preferable robustness and higher estimation accuracy of direction finding. Furthermore, a better estimation performance of the proposed algorithm can be obtained when UCA is used.  Then, we will examine the DOA estimation performance of the proposed algorithm and 2D LV-MUSIC algorithm versus number of snapshots by varying the number of snapshots from 50 to 500 with a step of 50. The SNR is set as 10 dB. The simulation results are plotted in Figure 6. We can easily find that for the estimation of both azimuth and elevation, the proposed algorithm can get a better estimation result than the 2D LV-MUSIC.
algorithm is much lower than the 2D LV-MUSIC.
Then, we will examine the DOA estimation performance of the proposed algorithm and 2D LV-MUSIC algorithm versus number of snapshots by varying the number of snapshots from 50 to 500 with a step of 50. The SNR is set as 10 dB. The simulation results are plotted in Figure 6. We can easily find that for the estimation of both azimuth and elevation, the proposed algorithm can get a better estimation result than the 2D LV-MUSIC. Figures 5 and 6 verified the effectiveness of the proposed algorithm. This is mainly because when the modified FB observation matrix is obtained by the CCSM observation matrix, the size of the data matrix yields an effective doubling, which achieves an expected improvement in the estimator variances. The improvement makes the proposed algorithm have preferable robustness and higher estimation accuracy of direction finding. Furthermore, a better estimation performance of the proposed algorithm can be obtained when UCA is used.  Figures 5 and 6 verified the effectiveness of the proposed algorithm. This is mainly because when the modified FB observation matrix is obtained by the CCSM observation matrix, the size of the data matrix yields an effective doubling, which achieves an expected improvement in the estimator variances. The improvement makes the proposed algorithm have preferable robustness and higher estimation accuracy of direction finding. Furthermore, a better estimation performance of the proposed algorithm can be obtained when UCA is used.

Running Time Comparison
The runtime of the proposed algorithm and LV-MUSIC algorithm under the same conditions is compared in this subsection. The simulation results are obtained using a PC with an Inter(R) Pentium(R) G2010 2.  Table 2. It can be seen that the running time of the proposed algorithm is lower than the 2D LV-MUSIC, namely, the proposed algorithm is more efficient than the 2D LV-MUSIC in computational complexity. This is mainly because the proposed algorithm transforms the computation of the spectrum function from the complex space into the real-valued space, so that the computing time is significantly reduced at each searching point. It is worth noting that the number of searching points in the LV-MUSIC algorithm with 4D search is several orders of magnitude larger than that of the algorithm with 2D search. Thus, its running time is unacceptable, and it is not considered as a comparison here.

Conclusions
This paper proposes a real-valued 2D MUSIC algorithm based on CCSM data matrix, which ensures high computational efficiency and can be applied to arbitrary centrosymmetric PSA. Due to the modified FB averaging presented, the complex covariance matrix and spectrum function are converted into real domain. As a result, the computational cost is reduced significantly. By utilizing the rank deficiency theorem, the spectrum function which only depends on the DOA parameters is obtained. As such, the 2D azimuth and elevation can be estimated without estimating polarization parameters. Furthermore, in practical implementation, the search space and search step can be set reasonably under different circumstances in order to minimize the time of DOA estimation with the satisfaction of estimation accuracy. Taking UCA and URA for examples, the DOA estimation performance and computational complexity are discussed by several simulations. The simulation results validate that the proposed algorithm outperforms the 2D LV-MUSIC algorithm and has a lower running time.