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.

1. 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 algorithm can only be applied to uniform rectangular array (URA) composed of six component electromagnetic vector sensors, leading to great restriction on the application. Moreover, the six component electromagnetic vector sensors have serious mutual coupling effect [25] and bring great economic burden. In addition, literature [26] presented an improved MUSIC algorithm based on PSA without ambiguity. It is quite applicable to the arrays with arbitrary geometries. However, the spectrum function of this algorithm is four-dimensional (4D) with respect to DOA and polarization parameters, leading to a 4D spectral peak search process. High computational load severely decelerates the parameter estimation, therefore, the reduction of computational complexity of MUSIC-based algorithm is a problem worth studying.

Root-MUSIC algorithm is a modified MUSIC algorithm which replaces the spectral peak searching by polynomial rooting, so as to reduce the computational complexity of DOA estimation [19,20]. Nevertheless, Root-MUSIC algorithm is generally used only for uniform linear array (ULA), since it requires a complicated preprocessing for other array structures. Reducing the number of searching points can also improve computation efficiency. For MUSIC-based algorithms using PSAs, decoupling the DOA estimation from polarization parameter estimation can reduce the dimensionality of the spectral peak searching and drastically reduce the number of iterations. A MUSIC-based solution using one-dimensional spectral peak searching has been proposed in [27] base on the ULA composed of cross-dipoles. On the other hand, it is well-known that processing the data in real domain is much faster than in complex domain. In subspace-based algorithms based on TSA, by utilizing the forward/backward (FB) averaging technique, the complex covariance matrix is converted into a real-valued matrix, thus the eigenvalue decomposition (EVD) can be calculated in real domain [28,29]. However, the FB averaging cannot be applied directly to the PSAs consist of EMVSs with multiple vibrators.

In this paper, a novel real-valued MUSIC algorithm based on arbitrary centrosymmetric PSA is proposed to achieve two-dimensional DOA estimation with low computational complexity and high accuracy. The merits are reflected in the following four aspects. First of all, the proposed algorithm is applicable to an arbitrary centrosymmetric PSA, even three-dimensional array geometry. Then, the modified FB averaging which can be applied to PSA is presented in the proposed algorithm to process the data in real domain, thus achieving a large reduction of computational cost. Furthermore, through decoupling the DOA estimation from the polarization parameter estimation, the 4D MUSIC spectrum function is dimension-reduced to a 2D function depending only on the DOA parameters. Hence, a joint search of both DOA and polarization parameters is avoided, which significantly reduces the number of iterations in searching the spectral peaks. The azimuth and elevation simultaneous searches (AESS) [30] can be applied to the proposed algorithm to perform a more efficient 2D DOA peak search process. Finally, electric co-centered crossed-dipole (ECCD) pairs with lower mutual coupling effect are used in the array and it will cause less compromise on the estimation performance than most other form of vector sensors. The simulation results demonstrate that the proposed algorithm can get accurate parameter estimation with a low computational complexity.

The remainder of this paper is organized as follows. Section 2 gives the array signal model used throughout this paper. Then, we explicitly describe the proposed algorithm in Section 3, and discuss the estimation performance and the complexity in Section 4. The simulation results are shown in Section 5 to verify the effectiveness of the proposed algorithm. Finally, the conclusion of this paper is given in Section 6.

Notations: We use , and to represent transpose, conjugate, and conjugate transpose of matrices, respectively. Denote operators and as kronecker product and Khatri-Rao (KR) product, respectively. and are used to indicate real part and imaginary part of complex variates, respectively. denotes the expectation operator, and denotes the determinant of matrices. The non-bold letters, lower-case bold letters, and upper-case bold letters denote scalars, vectors and matrices, respectively. Especially, denotes an identity matrix, and denotes an anti-identity matrix with ones on its anti-diagonal and zeros elsewhere.

2. Array Signal Model for DOA Estimation

2.1. 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] where and 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; and denote the auxiliary polarization angle and the polarization phase difference, respectively.

Figure 1

The coordinate system for the k-th source.

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 where means wavelength. For simplicity, we use and to respectively represent and in the following parts. As such, the received signal model (RSM) of array can be expressed as a long-vector form where symbolizes the zero-mean additive complex white Gaussian noise vector at the m-th vector sensor.

2.2. Conjugate Centrosymmetric Signal Model

Suppose that there a permute matrix where is an permute matrix, and it can be obtained readily that . Then, a novel signal model is proposed as follow

According to (3) and (4), it is obvious that the inter-vector-sensor spatial phase factors 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 such that the spatial steering vector of the proposed signal model can be expressed as where denotes a matrix of ones, denotes the top or rows of and denotes the bottom part. Denote the proposed signal model as the conjugate centrosymmetric signal model (CCSM). Furthermore, note that the matrix multiplications involving or 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 corresponding to each array is expressed as where denotes a null matrix.

Figure 2

Two kinds of centrosymmetric array. (a) A UCA with 12 ECCDs; (b) A URA of 3 rows and 4 columns with totally 12 ECCDs. The numbers represent the order of vector sensors.

The noise-subspace and spatial steering vectors of CCSM are orthogonal.

The covariance matrices of RSM and CCSM are 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.,

As such, we can obtain which illustrates that the noise-subspace and the steering vectors of CCSM are still orthogonal. ☐

3. Real-valued 2D DOA Estimation Based on Modified FB Averaging

3.1. The Modified FB Averaging

Considering about finite received array data, the covariance matrix of CCSM is expressed as where L is the number of snapshots and is the CCSM observation matrix.

In contrast to [34,35], define the modified FB covariance matrix as where is the modified FB observation matrix in which is a unitary permute matrix denoted as

It is noteworthy that, different from the original FB averaging [34], the modified FB averaging uses the unitary permute matrix 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 that

Then, let be a unitary matrix which is defined as where denotes an null matrix, and is defined similarly as (23). It is obvious that .

Therefore, we can obtain the real-valued covariance matrix as

By substituting (22) into (26), can be further simplified as

The second equation in (27) is derived readily by using (24). (27) indicates that only depends on the covariance of CCSM.

The real-valued covariance matrix is also symmetric due to the following properties

Thus, the EVD of require only real-valued computation [36]. Furthermore, the eigenvalues and eigenvectors of are real-valued as well.

The EVDs of and are where and are eigenvalue diagonal matrices and and are eigenvector matrices. The subscript S and N represent signal- and noise-subspace, respectively.

Substituting (30) into (26), can be expressed as

Compared (31) with (29), we have the following equations

Define as where is defined as

Then, we can simplify by using (26), (33) and (34)

Note that is a real-valued matrix.

Partition CCSM observation matrix as for even M, or for odd M. Then, substituting , , and yields

It is easy to determine that 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 does not involve general matrix multiplication, but is a simple construction based on (the CCSM data matrix) involving only additions. Therefore, it has significant consequences for processing modified FB covariance matrix.

3.2. 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, 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 only depends on polarization parameters. Therefore, we can obtain that where 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 is rank-deficient only when equals the actual DOA, i.e., . 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 in order to estimate the DOAs only through real-valued operations. Substituting (13) and (32) into (39), it can be obtained that where

Obviously, is a real-valued matrix since and is real-valued. From (39) and (44), it is clear that these two forms of are exactly identical. However, calculating 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 1	Obtain the CCSM observation matrix Y=1/L[y_1(t),⋯,y_L(t)].
Step 2	Compute the real-valued matrix Z_FB by (36).
Step 3	Compute the symmetric real-valued covariance matrix R_RV=Z_FBZFBT.
Step 4	Perform the EVD of R_RV by (30) to get V_N.
Step 5	Compute the spectrum function at each searching point via H(θ,ϕ)=D^T(θ,ϕ)V_NVNTD(θ,ϕ), and estimate 2D azimuth angle and elevation angle of K signals by {θ^_k,ϕ^_k}=argmax1/det(H(θ,ϕ)).

4. Discussion

4.1. 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 , elevation angle , auxiliary polarization angle and polarization phase difference 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 is given by

Hence, the sub-matrix of with respect to azimuth is

Similarly, the sub-matrices with respect to elevation and the cross terms of can be expressed as where denotes the trace of a matrix, denotes the i-th column of the identity matrix. is the covariance matrix of sources which is defined as . and are the derivations of with respect to and with .

Then, the CRB matrix with respect to azimuth and elevation is the inverse of , i.e.,

As such, the CRBs of azimuth and elevation of the k-th source are where denotes the (k,k)-th element of .

4.2. 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 from RSM requires about flops; (2) The real-valued EVD (based on symmetric QR algorithm [40]) of requires about flops; (3) Computing the spectrum function in (44) for each searching point requires about flops. Assume that there are 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. (1) Computing complex covariance matrix from RSM requires about flops; (2) The EVD of complex covariance matrix requires about flops; (3) Computing the spectrum function for each searching point requires about 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

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.

Figure 3

Computational complexity versus (a) the number of snapshots (for 12 sensors and 1° search step); (b) the number of sensors (for 300 snapshots and 1° search step), and (c) the search step of parameters (for 12 sensors and 300 snapshots).

5. 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 and the interspace of URA is . Two uncorrelated far-field narrowband signals whose DOAs and polarization parameters are and , 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 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.

5.1. The Simulation Results Distribution Scatter Plots

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.

Figure 4

Scatter plots of 2D DOA estimation for (a) the proposed algorithm and (b) 2D LV-MUSIC algorithm.

5.2. 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.

Figure 5

RMSE and CRB of (a) azimuth and (b) elevation versus SNR (for 300 snapshots).

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.

Figure 6

RMSE and CRB of (a) azimuth and (b) elevation versus number of snapshots (for SNR = 10 dB).

Figure 5 and Figure 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.

5.3. 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.8 GHz CPU and 8 GB RAM by running the Matlab (Ver. R2016b) codes in the same environment. Two signals are considered here. 200 Monte Carlo trails are carried out with 10 dB SNR and 300 snapshots. The azimuth angle and the elevation angle have been searched in the range of 0° to 360° and 0° to 90°, respectively. The average running time with different numbers of sensors and different search steps are shown in 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.

Table 2

Comparison of the average running time (s) with 200 Monte Carlo trails.

Search Step	Number of Sensors	Proposed Algorithm	LV-MUSIC (2D Search)
0.25	12	17.6096	21.4674
0.5	12	4.8517	5.5413
1	12	1.1539	1.4175
0.5	6	3.9773	4.5010
0.5	18	6.2451	8.5407

6. 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.