A Fast DOA Estimation Algorithm Based on Polarization MUSIC

A fast DOA estimation algorithm developed from MUSIC, which also benefits from the processing of the signals’ polarization information, is presented. Besides performance enhancement in precision and resolution, the proposed algorithm can be exerted on various forms of polarization sensitive arrays, without specific requirement on the array’s pattern. Depending on the continuity property of the space spectrum, a huge amount of computation incurred in the calculation of 4-D space spectrum is averted. Performance and computation complexity analysis of the proposed algorithm is discussed and the simulation results are presented. Compared with conventional MUSIC, the proposed algorithm has considerable advantage in aspects of precision and resolution, with a low computation complexity proportional to a conventional 2-D MUSIC.


Introduction
The direction of arriving signals (DOA) estimation problem has raised great interest in various fields, including radar, mobile communications, and microphone array systems. In the past decades, numerous algorithms and techniques have been proposed. Some important classes of algorithms include the Maximum Likelihood Estimation (ML) [1,2], the signal subspace based algorithms represented by Multiple Signal Classification (MUSIC) [3], and the computationally efficient Estimation of Signal Parameters via Rotational Invariance Techniques (ESPRIT) [4]. These algorithms help to enhance the performance of the DOA estimation substantially, and a great number of relevant algorithms are derived from them, for instance, Stochastic Maximum Likelihood Estimation (SML) [5], Root-MUSIC [6], Total Least Squares version of ESPRIT (TLS-ESPRIT) [7], etc.. In recent years, sparse signal presentation has been introduced into DOA estimation [8,9,10], making the DOA estimation with small number of snapshots possible.
Of all these algorithms, ML and its derivate algorithms are of the best performance in both accuracy and resolution, whereas the computation complexity is extremely high, preventing their practical applications. ESPRIT and MUSIC algorithms, proposed later, with precision error approaches to the Cramer-Rao bound (CRB) as snapshot number of array outputs accumulates to a certain threshold, and the calculation incurred in these 2 algorithms are acceptable.
Numerous factors should be considered to select a proper algorithm. First, the DOA algorithm should adapt to the form of the array. Generally, the signals' azimuths need to be estimated, and in certain circumstance, their elevation estimation are also demanded. Well-designed planar arrays, such as uniform circular array (UCA), can meet these requirements. However, ESPRIT and Root-MUSIC can only work on linear array, and algorithms derived from them for other array forms [11,12,13] can only work on particular array patterns, while only 1 direction parameter can be estimated. The MUSIC algorithm, which can be transplanted to more general forms of planar arrays and provides both azimuth and elevation estimation, seems to be more suitable. Current calculation capacity of DSP chips makes it practical as well.
For algorithm selection, another crucial factor is the performance, especially in aspects of precision and resolution. The results presented in [14,15] suggested that the accuracy and resolution in MUSIC are both determined by three factors: SNR, the number of array's snapshots, and the array's aperture. In practical circumstances, the SNR is generally uncontrollable, the aperture is often restricted, and the numbers of samples is limited. From this perspective, the accuracy and resolution of MUSIC algorithm is confined, and further improvement in performance can only be sought within other methods.
The polarization phenomenon of electromagnetic wave has been noticed and studied by physicists for more than a century. However, it was only in 1980s that this phenomenon began to attract the researchers' attention and proved to be of great advantage in signal processing. The polarization sensitive array is introduced into DOA estimation by Nehorai and Li [16,17]. In [17,18], Jian Li integrated the polarization into ESPRIT successfully. The pencil-MUSIC algorithm is devised for DOA and polarization estimation in [19,20,21], which extend the ESPRIT-like algorithm to partially polarized sources. [22,23] provide an algorithm based on sparse ULA, consisting of sensors which can measure all six electrical and magnetic components. [24] provide several algorithms to solve the correlation problem with the array of vector sensors. In [25,26], the polynomial rooting algorithm was adopted in polarization sensitive arrays. Covariance tensor [27,28] is introduced to define the array manifold in a new way and reduce the amount of calculation. These papers exploited the advantage of polarization of signals. However, to reduce the calculation complexity, all these algorithms assumed specific patterns of array, and most of them adopted the sensors with the capacity to detect all 6 electrical and magnetic components, which is not very applicable.
The algorithms presented in [29][30][31][32][33] require no assistance of magnetic components. But the algorithms mentioned in [29,30] assumed that 1 polarization parameter and 1 direction parameter are predetermined, which cannot be satisfied generally. In [31,32,33], the amount of calculation is reduced by dividing the process of estimation into two steps: one concerns with the polarization, and the other with the direction. The drawback of this method is that the second step in the algorithm depends on the success of previous step. When the DOAs of multiple signals are very close, the second step inclines to fail, even when the signals are distinct in polarization domain.
To retain the advantage in resolving close signals offered combination of spatial and polarization domains, a new algorithm, which reveals the space-polarization spectrum of 4-D space in a new way, is proposed in this paper. Unlike the algorithms mentioned above, the proposed algorithm can be used on various forms of polarization sensitive arrays. Besides, in case the multiple signals are close in space domain, they can be resolved as they are distinctive in polarization domain. By analyzing the property of the space spectrum, the algorithm is designed to reduce the complexity of the 4-D space spectrum calculation to be proportional to a 2-D spectrum calculation [34]. The proposed algorithm makes the transplantation of MUSIC from arrays of ideal omnidirectional antennas to polarization sensitive arrays more practical. This paper is organized as follows. In Sec. 2, the fundamental principle of MUSIC is introduced, and the transplantation of MUSIC into polarization sensitive array is discussed. Although the latter one is highly inapplicable because of the immense amount of calculation, it is essential in developing the faster new algorithm. In Sec. 3, the definition of zero spectrum and another form of space spectrum are given, and the continuous and differential properties of them are studied. Based on these properties, calculation complexity in the 2-D peaks searching of polarization domain can be reduced dramatically. In Sec. 4, the analysis on the accuracy and resolution performance is presented, and the analysis of algorithm complexity is discussed as well. In Sec. 5, simulation results are provided to demonstrate the continuous property mentioned in Sec. 3, and the accuracy and resolution performance.

MUSIC's Transplantation into
Polarization Sensitive Array

Principle of Conventional MUSIC
Consider an arbitrary array (not a polarization sensitive array) which consists of M elements, and N narrow-band zero-mean signals of the same central frequency, arriving at the array from different directions. It is assumed that these signals are uncorrelated with each other, and M > N. The output of the ith element can be written as: where s k (t) is the kth signal, n i (t) is the Additive White Gaussian Noise (AWGN), and ω 0 is the central frequency of all the N signals. The parameter τ ki is the delayed time of the kth signal at the ith element, whose form is determined by the direction of impinging signals and the pattern of the array. For a ULA, τ ki is (i − 1)d sin ϕ k /c, where ϕ k represents the azimuth of the kth signal, d is the element spacing and c is the propagation velocity of signals; for a UCA with M elements, τ ki becomes r sin θ k cos(ϕ k − 2πi/M)/c, where θ k and ϕ k are the elevation and azimuth of the kth signal respectively, and r represents the radius of the UCA. For general planar arrays, the time elpase τ ki becomes (x i sin θ k cos ϕ k + y i sin θ k sin ϕ k )/c, where x i and y i are used to indicate the location of the ith element. Using vector forms, (1) can be written as: where the output vector is X(t) = [x 1 (t), . . . , x M (t)] T , the signal vector is S(t) = [s 1 (t), . . . , s N (t)] T , and the noise vector is N(t) = [n 1 (t), . . . , n M (t)] T ; A is the matrix of array manifold, and the element on its kth row and ith column is The operator ( * ) T indicates the transposed matrix.
Then the covariance matrix of X(t) is: where λ i represents the eigenvalue of R X , arranged as λ 1 λ 2 · · · λ N+1 = · · · = λ M = σ 2 N . σ 2 N is the power of white noise, and u i denotes the eigenvector corresponding to λ i . I is the identity matrix. From (3) we have: Equation (4) indicates the fact that the column vectors of the matrix A, which are also referred as steering vectors, are within the N-dimension space spanned by the first N eigenvectors, perpendicular to the noise subspace spanned by the other M − N eigenvectors.
In MUSIC algorithm [3], the form space spectrum is defined as: where the operator * denotes the Euclidean norm, and U N = [u N+1 , u N+2 , . . . , u M ], whose column vectors span the noise subspace. Corresponding to the signal subspace, The form of vector a is similar to that of the column vectors of A. Every component of a is a function of the 2 direction parameters : the azimuth ϕ and elevation θ. The space spectrum P can be considered as a function of them as well. By calculating the spectrum in the space domain and searching the values of the parameters corresponding to the N maxima of P, the directions of the arriving signals can be estimated. The amount of calculation incurred in this process depends on the number of parameters determining the signal's direction. For ULA, all signals are assumed to be within a planar, and a 1-D search within the range of azimuth ϕ can estimate their directions; whereas for a UCA or a planar array, a 2-D search on the azimuth ϕ and elevation θ is required, and the amount of calculation increases considerably.

Polarization MUSIC
Consider a UCA of M polarized sensitive elements (as shown in Fig. 1, where M = 8). Assuming N signals arrive at the array, as mentioned in [17,18], the ith column vector of A in (2) becomes:  where ⊗ represents the Kronecker product, γ i and η i are the polarization parameters of ith signal, and a S (θ i , ϕ i ) is the steering vector of the ith signal on behalf of its spatial direction, written as: In a UCA the q k is: The vector a P (θ i , ϕ i , γ i , η i ) which represents the effect of polarization is: Substitute (6) into (2), following the similar program in [29] we can define a new joint space spectrum P J : Comparing (10) with (5), the only difference is that the space spectrum P J becomes a function of 4 parameters: θ, ϕ, γ and η, specifying the elevation angle, azimuth angle, and the 2 polarization parameters respectively. It should be noticed that, even if polarization states of all arriving signals are the same, as the direction parameters θ i and ϕ i are involved in vector a P (θ i , ϕ i , γ i , η i ), it still brings difference into the steering vector a(θ i , ϕ i , γ i , η i ) and helps to enhance the algorithm's resolution performance. However, when the combination of these 4 parameters makes the polarization states no longer provide assistance in resolving these 2 signals.
Theoretically, a peak search can be executed in a 4-D space spanned by these 4 parameters to obtain the estimation results. Benjamin Friedlander has proved that the performance of resolution will be improved considerably [35]. However, the calculation complexity of the 4-D space spectrum causes an unbearable computation complexity and makes the ordinary polarization MUSIC algorithm unpractical.

Fast Polarization MUSIC
In this part, properties of the space spectrum P J (θ, ϕ, γ, η) of the polarization MUSIC is analyzed, and a more effective algorithm is developed to reduce the great amount of calculation in DOA estimation.

Property of the Joint Space Spectrum
As previously stated, the values of the 4 parameters corresponding to the N maximums of the space spectrum in (10) indicate the directions and polarization states of the N arriving signals. And a maximum in 4-D space spanned by the 4 parameters, whose position is (θ 0 , ϕ 0 , γ 0 , η 0 ), remains to be a maximum in the 2-D subspace (θ 0 , ϕ 0 , γ, η), spanned by the two polarization parameters.
As mentioned in [36], another form of the joint space spectrum in (10) can be derived: where: and: where Z(θ, ϕ, γ, η) is named as zero spectrum, as it approaches to 0 when values of the 4 parameters approach to the direction and polarization state of a signal. At the same time, U(θ, ϕ, γ, η) in (13) approaches to a local maximum, and a steep peak appears in the graph of P J (θ, ϕ, γ, η).
Suppose θ = θ 0 and ϕ = ϕ 0 , and substitute (11) into (15), we have: The operator Re( * ) is the real part of a plural variable, whereas Im( * ) is the imaginary part. u i and v i are determined by the form of the array, and the values of θ 0 and ϕ 0 .
As the covariance matrix R X is determined after the output of the array are sampled, U S can be regarded as a constant matrix. Then the coefficients x r , y r , a, b, c, d and φ are all constants determined by θ 0 and ϕ 0 , irrelevant with the value of γ and η. The detailed form of these coefficients can be deduced without much difficulties and omitted here.
As we know, the maximum of a continuous function with two variables has to meet the following expressions: Substitute (16) into (17), Substitute (16) into (19), Notice that γ ∈ [0, π/2], and U(θ 0 , ϕ 0 , γ, η) will be constant when γ = 0, which rarely happens. Based on (21) and (22), when γ = 0, the value of η corresponding to the maximum is: where η m indicates the value of η at the position of maximum. k should be 0 or 1, to ensure that η m ∈ [0, 2π]. The value of η m is not affected by γ.
It should be noticed that, in the derivation of (23) and (25), no specific pattern of arrays are involved. This is the very reason that the algorithm proposed in this paper has no particular requirements in the form of arrays.
Step 3: Take another 2-D peak search in P D (θ, φ). The direction parameters of the peaks are the estimation of the signals' direction, and the corresponding values of Γ(θ 0 , ϕ 0 ) and H(θ 0 , ϕ 0 ) provide the estimation of the signals' polarization states.

The Procedure of the Fast Algorithm
Obviously, the process mentioned in 3.1 cannot reduce the computation complexity. However, in the first 2-D search, based on the fact that only one peak exists in P J (θ 0 , ϕ 0 , γ, η), we provide a straightforward method to reduce the amount of calculation.
According to (23), η m is determined by φ at the location of the maximum of U(θ 0 , ϕ 0 , γ, η), which is independent of γ m . So we can assign a constant to γ, and perform a 1-D search along the parameter η to get η m ; then set η as η m , run another 1-D search to get γ m . P D (θ, ϕ), Γ(θ, ϕ) and H(θ, ϕ) can be obtained respectively.
Referring to (16), the a, b, c, d, and φ are all continuous functions of θ 0 and ϕ 0 . From (23) and (25), it can be derived that γ and η, i.e., the Γ(θ, ϕ) and H(θ, ϕ), are also continuous functions of θ and ϕ under general conditions. In other words, the Γ(θ, ϕ) and H(θ, ϕ) will change smoothly if the parameters change gradually in the calculation process of P D (θ, ϕ), and the location of the single peak in P J (θ 0 , ϕ 0 , γ, η) does not change dramatically when the search step of θ and ϕ are limited to small values.
In (16), It should be noticed that φ is derived from an inverse tangent function of d/c. This may cause a flip of π in the value of φ, and H(θ, ϕ) might be affected. However, (25) indicates that Γ(θ, ϕ) will remain unaffected. This flip does not happen frequently, and even if the flip occurs, as shown in (23), the values of H(θ, ϕ ± ∆ϕ) or H(θ ± ∆θ, ϕ) corresponding to maximum will be near H(θ, ϕ) ± π, which is predictable.
The main steps of fast polarization MUSIC are summarized as follows: Step 1: Calculate the covariance matrix R X of the array's output, take an eigendecomposition to get the orthogonal basis of the noise subspace.
Step 2: Select a set of direction parameters θ and ϕ randomly, or according to the result of previous estimation. Use (10) to complete a peak search on the two polarization parameters. Record the maximum value and the corresponding values of Γ(θ, ϕ) and H(θ, ϕ).
Step 5: Exert a 2-D peak search in the P D (θ, ϕ) to obtain all maximums. The directions of impinging signals can be estimated from direction parameters of these peaks.

Precision
The algorithm presented in 3.2 is an improved version of polarization MUSIC, so some analysis results for conventional MUSIC can be applied here. The analysis for the performance of MUSIC, given by Stoica and Nehorai in [15,37], are adopted here.
The theorem 3.1 in [15] states that the MUSIC estimation errorsω i − ω i approach to jointly Gaussian distribution. Here ω i indicates the parameter to be estimated in the MUSIC, which can be either θ i or ϕ i ; andω i is the corresponding estimation result. The theorem also states that the means ofω i − ω i are all 0. When only one signal arrives at the array, it becomes Under this condition, λ 1 is proportional the power of the signal, so it is clear that the variance is largely determined by the SNR and the number of snapshots. It is reasonable to conclude that the variance of results obtained by MUSIC and the proposed algorithm are of the same order if only one signal arrives. Computer simulations are presented to demonstrate the precision performance when multiple signals arrive.

Resolution
In this subsection, the resolution performance of the fast polarization MUSIC is discussed.
In the paper written by Benjamin Friedlander [35], a standard has been given to justify whether two close signals can be descried, when using a ULA. If there are more than 2 signals, the analysis will become extremely difficult, so we only consider a 2 signals' scenario, as the analysis presented in [35] did.
With some conclusion from [38], we have (refer to Appendix B for more details): .
It can be seen that the value of the right side of expression (33) is determined by a 0 , a 1 and a 2 (for u 1 and u 2 are also determined by them). Then when the values of snapshot number K and SNR satisfies expression (33), the two close signals can be descried by the algorithm theoretically.
For the conventional MUSIC, some minor adjustment in deduction is required (see Appendix B), and the necessary snapshot number becomes Equations (33) and (34) provide a way to compare the numbers of snapshots needed in the polarization MUSIC with that of a conventional MUSIC to resolve two signals.

Complexity
In this subsection, the computation complexity of the proposed algorithm is discussed. In Subsection 3.2, the process of the algorithm are divided into 5 major steps.
Step 1 and step 5 concern with the eigendecomposition and peak search in a 2-D matrix. Comparing with the amount of calculation incurred in step 2 to step 4, the computation in step 1 and step 5 is negligible. So we focus on the computation from step 2 to step 4.
The steering vector a(θ, ϕ, γ, η) can be calculated offline, and the noise spectrum matrix U N is calculated only once. Each value of space spectrum in (10) can be obtained by two matrix multiplications. Besides, the computation complexity for each point of spectrum is the same for both conventional MUSIC and polarization MUSIC, if the numbers of their array antennas are the same. So the overall computation complexity is determined by the numbers of spectrum values to be calculated. Elevation angle θ, azimuth angle ϕ and the polarization parameters γ and η are all continuous, whereas they need to be sampled in the process of space spectrum calculation. Assume that the sample numbers of 4 parameters are N 1 , N 2 , N 3 and N 4 , respectively. To obtain the required space spectrum for a conventional 2-D MUSIC, the number of sampled space spectrum is N 1 N 2 , whereas for the polarization MUSIC mentioned in 2.2, it becomes N 1 N 2 N 3 N 4 . For the fast polarization MUSIC algorithm, the required number of spectrum values is discussed as follows.
In subsection 3.2, it is stated that two 1-D peak searches along η and γ would be sufficient to obtain the maximum of U(θ 0 , ϕ 0 , γ 0 , η 0 ). When the worst condition occurs, we need to go through both the ranges of η and γ to locate the peak, and the numbers of spectrum value required to be calculated in step 2 or step 3 would be no more than N 3 + N 4 .
Other conditions are omitted here, since their deductions are similar.
Based on the discussion above, it can be concluded that to obtain the value of each P D (θ, ϕ), at least 5 points of the spectrum are required. So the overall number of spectrum values calculated from step 2 to step 4 is between 5N 1 N 2 and (N 3 + N 4 )N 1 N 2 . In Section 5, simulation results show that the required number of spectrum is approximately 6N 1 N 2 .

Character of Polarization Spectrum
The fast polarization MUSIC algorithm can be used on various sorts of polarization sensitive arrays. In this section, a UCA of 8 polarization sensitive elements with a radius of 0.667 λ are adopted in the following simulation. The output of the array is contaminated by uncorrelated AWGN. Most of the following simulation results are based on these conditions. The directions of the arriving signals can be expressed in various forms, and the form adopted in this section is illustrated in Fig. 1, with a minor difference with those in (8) and (9). The azimuth and elevation here are written as ϕ and θ , respectively, and the relationship between the two forms can be expressed as: sin θ = sin θ cos ϕ.
In this simulation, 3 signals impinge onto the array. The number of snapshots is 500, and the SNR is 15 dB. P(0, 0, γ, η) and P(2.5, 2.5, γ, η) are plotted in Fig. 2(a) and Fig. 2(b), respectively. Fig. 2 indicates that regardless of the number of arriving signals and their polarization, in P(θ 0 , ϕ 0 , γ, η), one and only one maximum exists. From the contour graphs, the symmetrical character of parameter η can be observed more clearly. Based on this, the maximum can be located after two 1-D searches.

5.2
The Continuity in Γ(θ , ϕ ) and H(θ , ϕ ) Fig. 3 presents Γ(θ , ϕ ) and H(θ , ϕ ) with 3 arriving signals of various parameters, to illustrate the continuity property of them. Fig. 3(a) presents the graph of Γ(θ , ϕ ), whereas Fig. 3(b) present H(θ , ϕ ). The result of 3 signals, whose direction and polarization parameter written in the form of (θ , ϕ , γ, η) are (−3, −4, 12.5, 30), (3, −3, 45.5, 50) and (−3, 4, 70.5, 70), is drawn in Fig. 3(a) and 3(b). It can be observed that in the two figures, the Γ(θ , ϕ ) and H(θ , ϕ ) are continuous in most area, whereas the steep slopes occur near the directions of the impinging signals. Fig. 4, Fig. 5 and Fig. 6 provide stimulation results to compare the performances of the fast polarization MUSIC algorithm and other algorithms. The performance of the fast polarization MUSIC is compared with conventional MUSIC, and the rank reduction (RR) algorithm in [25,26], which can be extended to arbitrary form of array as well. The performances of these algorithms are presented through Monte Carlo simulations.   Fig. 4 draws snapshot numbers required to descry 2 close signals with various angle differences, calculated by (33) and (34). Their azimuth angles are 0 • , one signal's elevation angle is set to 0 • , and the other is modified. The simulation is based on a UCA of 16 sensors, whose radius is 0.667 λ. The polarization parameters are (30,15) and (50, 100). In Fig. 4, the calculated results of the fast polarization MUSIC and conventional 2-D MUSIC are given according to the SNR. It can be seen that less snapshots are required to descry 2 close signals in the proposed algorithm, particularly with relatively low SNR.     Fig. 6 shows that the success rates of the fast polarization MUSIC are generally higher than that of conventional 2-D MUSIC and RR algorithm, particularly with small snapshot numbers.

Precision and Resolution
Especially, the success rate of the proposed algorithm is approximately 80% in the 4-signal scenario, higher than that of RR algorithm and the conventional 2-D MUSIC. Figure 7 gives the space spectrum P(θ , ϕ ) of the conventional 2-D MUSIC, the RR algorithm and the graph of P D (θ , ϕ ) of the fast polarization MUSIC . The parameters of the 4 incoming signals are the same as the 4-signal scenario in Fig. 6, with 500 snapshots. Fig. 7(a) shows that the peaks of the conventional 2-D MUSIC are eclipsed severely and cannot be distinguished. For the RR algorithm, the peaks in Fig. 7(b) can be located, although not very distinct. But in the space spectrum of the fast polarization MUSIC, all 4 peaks are quite clear.

Computation Complexity
In simulation results presented in this article, the ranges of azimuth and elevation are sampled with an interval of 0.5 • , and the ranges of azimuth and elevation angles are both limited from −15 • to 15 • , i.e., θ ∈ [−15, 15] and ϕ ∈ [− 15,15]. Then the 2-D space spanned by azimuth and elevation will be sampled with 61 × 61 = 3721 samples, i.e., if the conventional MUSIC is adopted, 3721 points of the space spectrum are required to be calculated to estimate the arriving signals' directions.
In the fast polarization MUSIC, suppose γ and η are sampled with intervals of 3 degrees and 6 degrees, respectively. Based on the analysis in Subsection 4.3, the lower limit of sample number needed will be 5N 1 N 2 = 18605. When the worst case happens, the number becomes (N 3 + N 4 )N 1 N 2 = 338611, whereas the number of samples needed in the conventional 4-D polarization MUSIC is N 1 N 2 N 3 N 4 = 6921060, approximately 20 times more than the worst case of the proposed algorithm.
Tab. 1 presents the average number of samples of space spectrum calculated in polarization MUSIC. The Monte Carlo experiments number is 100. The impinging signals' parameters and the configuration of the array are the same as Fig. 5 and Fig. 6. The data in Tab. 1 shows that the average sample numbers is approximately 6N 1 N 2 , and is insensitive to the snapshot number and/or the signals' number. So the calculation amount in the fast polarization MUSIC is as 6 times as that of conventional 2-D MUSIC needs. The fast polarization MUSIC has the same performance in precision and resolution as the original 4-D search polarization MUSIC, but the low computation complexity property makes it more attractive in practical application.
The computation complexity of the RR algorithm, is approximately 4 times of the conventional 2-D MUSIC [26], whereas the amounts of processed points of them are same. It should be noticed that the computation complexity of fast polarization MUSIC is a bit higher than that of RR algorithm, but they are of the same order of magnitude, whereas the former has a better performance in precision and resolution.

Conclusion
In this paper, a novel fast polarization MUSIC algorithm, which could reduce the considerable amount of computation, is proposed. Based on the fact that the difference of polarization between two close directions is minor or predictable, the daunting 4-D search is simplified by converting into two 2-D peak searches. The performance of fast polarization MUSIC is evaluated by theoretical analyses and simulations. It delivers improved precision and resolution with low computation complexity. The fact that this proposed algorithm is not confined to any specific pattern of polarization sensitive arrays expands its range of application.