Propagator Method using PARAFAC Model for Two Dimensional Source Localization

In this paper, we addressed the problem of estimating the two-dimensional (2D) Direction of Arrival (DOA) elevation and azimuth angles for multiple sources. The proposed method employs Propagator Method (PM) in conjunction with parallel factor (PARAFAC) model using a new antenna array configuration. The proposed method overcomes two main drawbacks in the existing 2D DOA schemes: use of high computation eigenvalue decomposition (EVD) or singular value decomposition (SVD), and complex pair-matching methods for elevation and azimuth angles in case of multiple sources. Therefore, significant reduction in computational load and complexity is achieved. Computer simulations demonstrate the effectiveness of the proposed method.


Introduction
Two-dimensional direction-of-arrival (2D DOA) estimation plays an important role in the field of array signal processing especially in applications such as radar, sonar, mobile communication systems, etc. Conventional methods such as MUSIC [1] and ESPRIT [2] and their variants have been widely used owing to their high-resolution DOA estimation performance.However, these methods incur high computational cost as they require either eigenvalue decomposition (EVD) or singular value decomposition (SVD) operations to decompose the received data matrix into a signal subspace and a noise subspace.In the presence of multiple sources, the problem of computation complexity becomes more acute.Another problem that needs to be contended with is the pair-matching of azimuth and elevation angles for each of the sources.
The problem of computational cost can be addressed by considering simpler and less complex techniques [3][4][5] which do not require either EVD or SVD.The propagator method (PM) [5] is one such method which employs linear operations with a least square criterion.The PM has been applied to different array geometries such as uniform linear arrays, L-shaped and parallel-shaped arrays [6][7][8][9].However, the PM applied to parallel-shaped arrays requires pair-matching of azimuth and elevation angles and also suffers from estimation failure problem when elevation angles are between 70° and 90° [7].The L-shaped array also suffers from pair-matching of the two electric angles estimated from the two orthogonal linear sub-arrays which form the L-shape.The work reported in [8] combines the PM method with ESPRIT algorithm and exploits the conjugate symmetry property of the array manifold matrix for DOA estimation without pair-matching problem.One drawback of this method is its relatively high computational cost due to its use of EVD.
The 2-D DOA estimation algorithm reported in [10], estimates the elevation angle based on the polynomial root methods (such as fast root MUSIC or ESPRIT) using a 1-D uniform linear subarray along the z axis.To estimate azimuth angle, it uses the elevation angle estimated earlier and a 2-D uniform linear subarray along the x axis based on a subspace method (such as 2-D MUSIC).This algorithm does not require pair-matching, but its drawbacks are relatively high estimation errors at low SNR and the need for a large number of snapshots.A cumulant-based direction finding algorithm has been proposed in [11] which resolves the pair-matching and aperture loss problems seen in other methods.Another method [12] which resolves the pairmatching problem does so by employing a trilinear decomposition-based [17] blind 2D DOA estimation algorithm for an L-shaped array.One drawback of this method is that it requires higher number of snapshots for estimation accuracy and higher computation cost for constructing the cross-correlation matrices.The work proposed in [13] employs PM-like method for an L-shaped antenna array.It employs cross-correlation of the received signals and uses linear operations to reduce the computation complexity.However, it does not resolve the pair-matching problem and therefore does not work for multiple sources.The pairmatching problem can also be alleviated by considering a parallel factor analysis (PARAFAC) [14] model which has become popular lately in array signal processing.
PARAFAC was used to solve the problem of parameter pairing (or pair-matching) for the L-shaped array in [15] and for the conformal array in [16].
In this paper, we propose a novel 2D DOA estimation method that employs two parallel linear antenna arrays, and uses the PM in conjunction with a PARAFAC model to avoid pair-matching problem.The methods for 2D DOA estimation employing PARAFAC model in [11], [12], and [15] are based on dividing the whole array into several subarrays and then forming the cross-correlation matrices using second-order statistics.Construction of each crosscorrelation matrix requires O(4N 2 L) FLOPs where N represents the number of antennas and L represents number of snapshots.In the proposed method, the propagator method can be applied directly to the data matrix with an order of O(3(N -1)L) or to a second-order statistics of the autocorrelation and cross-correlation matrices with an order of O(3(N -1) 2 L).In addition, the proposed method requires (N -K)  K  3 dimensional three-way array (TWA) required for PARAFAC model whereas the methods in [11], [12], and [15] require N  N  K.As a drawback, this will increase the complexity, processing and convergence time, and memory storage requirements for the methods in [11], [12], and [15].Similarly, while the method in [16] is accurate in DOA estimation and does not require pair-matching, it suffers from higher computational complexity.Compared with existing methods [10][11][12], [19], the proposed method has lower computational complexity and solves the pairmatching problem with multiple sources.Therefore, significant reduction in computational load and complexity is achieved.
The following notations are used throughout this paper: the operations (.) H , (.) † , (.) T , and (.) -1 represent conjugate transpose, pseudo-inverse, transpose, and inverse, respectively.A scalar is represented by u, a constant by U, a vector by u, a matrix by U, an ij−th member of a matrix U by u ij , and a three-way array (TWA) by U. The operations diag(.)and diag -1 (.) represent the conversion of a vector to diagonal matrix and diagonal matrix to a vector, respectively.

Array Geometry and Signal Model
The proposed parallel-shaped array geometry is shown in Fig. 1.The distance between adjacent antenna elements is d, where d = λ/2 with λ being the wavelength of the incident waveform.The arrays are divided into three subarrays: x a , y b , and z c .Each linear subarray consists of N antenna elements.Consider K narrowband non-coherent sources in the far-field region of the antenna array.The received signal vectors at t th snapshot for x a , y b , and z c subarrays from K sources can be represented as: where 2 cos exp j for the k th source.The diagonal matrices  1 and  2 in (2) and (3) contain information about the elevation and azimuth angles, which are defined as: where  for the k th source.The dimensions of  1 and  2 are K  K.

Proposed DOA Estimation Method
In this section, we propose a new DOA estimation method for non-coherent sources by using PARAFAC [14] model for pair-matching.The auto-correlation matrix R xx of received signal vector x a (t) in ( 1) is obtained as: Assuming uncorrelated noise between antenna elements on subarray x a , the above auto-correlation matrix reduces to: where R s = E[s(t)s H (t)] is the auto-correlation matrix of the signal vector s(t) and I is an identity matrix.The crosscorrelation matrix R yx of received signal vectors at subarrays y b and x a , with the same assumption, we get: Similarly, the cross-correlation matrix R zx or received signal vectors at subarrays z c and x a is: where the parameters (,) in  2 and  in  1 are omitted for simplicity.The above correlation matrices are concatenated to form a new matrix F as follows: where the dimensions of F are 3N  N. The matrix F is partitioned into two parts by employing propagator method (PM) [6], [9].This is given as ] T , where the dimensions of F 1 and F 2 are K  N and (3N -K)  N, respectively.The least squares solution for propagator matrix P ̂, similar to [9], is given by: where the dimensions of P ̂ are K  (3N -K).The propagator matrix P ̂ is partitioned as:        P P P P P P (13) where the dimensions of P ̂1, P ̂2, P ̂3, P ̂4, and P ̂5 are (N -K)  K, K  K, (N -K)  K, K  K, and (N -K)  K, respectively.From ( 8), (9), and (10), and following the same procedure in [9], the partitions P ̂1, P ̂3, and P ̂5 are defined as: where A 1 and A 2 are the partitions of (again, the parameter  is omitted in A for simplicity).The dimensions of A 1 and A 2 are K  K and (N -K)  K, respectively.The (N -K)  K  3 dimensional TWA required for PARAFAC model can be obtained by using P ̂1, P ̂3, and P ̂5 as follows: :,:,1 :,:, 2 :,:,3 where Q denotes the noise matrix.Comparing the above TWA with the PARAFAC model in [14], we can write: and the operation  denotes the Kronecker product.The alternating least squares (ALS) method [17] is applied to get the estimates of D, B, and C.
The COMFAC algorithm [18] is used to fit the PARAFAC model to the TWA G.For fast implementation of alternative least squares method to solve PARAFAC model, COMFAC algorithm is employed which speeds up the least squares fitting and reduces the complexity by utilizing a compressed version of the three-way data into smaller matrix dimensions.The algorithm outputs the identification matrices D ̂, B ̂, and Ĉ.The matrix Ĉ contains information about the elevation and azimuth angles.The diagonal matrices 1 and 2 are estimated from the identification matrix Ĉ as The pair-matching is taken care of using the PARAFAC model.We derive the Cramer-Rao Bound (CRB) for the proposed method in a similar manner as shown in [20].The received signals in (1), (2), and (3) are concatenated as shown in the following: CRB can be expressed by using the above signal model as: , m i is the i th column of M, ⊚ represents Hadamard product, and Table 1 shows the complexity in terms of number of FLOPs for major processing steps of the proposed method, universal 2D DOA estimation [10], and joint elevation and azimuth angle estimation [19].It can be inferred from the table that the proposed method has much lower computational complexity compared to other methods, where N, L, and K denote the number of antenna elements, snapshots of the received signals, and number of sources, respectively.
To illustrate the superiority of the proposed method in terms of complexity compared with other methods in [12], [11], [10], and [19], the case of N = 15, K = 2, and L = 100 is considered and actual FLOPs calculated.Table 2 shows the number of FLOPs required for this case.It is clear from the table that the proposed method requires the least number of FLOPs.

Summary of the Proposed Algorithm
As described earlier, the proposed algorithm considers the extended parallel-shaped array as shown in Fig. 1.The problem of pair-matching between two or more sources is avoided using the PARAFAC model.Following steps summarize the proposed algorithm to estimate DOAs for non-coherent sources without the problem of pairmatching.
Step 1: Construct the data matrices x a (t), y b (t), and z c (t) as shown in ( 1), (2), and (3), respectively, from L snapshots of the received signals.
Step 2: Determine the autocorrelation matrix R xx as in (8) and the cross-correlation matrices R yx and R zx as in ( 9) and ( 10), respectively.
Step 3: Concatenate the matrices in Step 2 to form a new matrix F as in (11).
Step 6: Apply alternating least squares (ALS) to ( 17) to get the matrix C which contains information about  1 and  2 .
Step 7: Estimate the elevation and azimuth angles using ˆând   k k as in (19), for the each of the K sources.

Simulation Results
This section investigates the performance of the proposed method through simulations.The performance is compared with the following methods: novel L-shaped method in [12], cumulant-based method in [11], joint elevation and azimuth angles estimation method in [19], and universal 2D DOA estimation method in [10].The sources considered are non-coherent.The assumed spacing between the adjacent antenna elements is d = λ/2 in all the simulations.For fair comparison, same numbers of antennas are considered in all algorithms.The plots are obtained by averaging over several runs of the algorithm.
Figure 2 shows standard deviation (STD) against SNR from 0 to 30 dB for the proposed method, novel Lshaped in [12], cumulant-based method in [11], and joint L-shaped method in [19], universal 2D DOA in [10], and CRB.We observe that the proposed method has the lowest STD and it is comparable with CRB at SNR beyond 15 dB.
Figure 3 shows the standard deviation (STD) against increasing number of snapshots for the proposed method, cumulant-based method in [11], and joint L-shaped method in [19].The number of antenna elements and number of sources considered here are N = 15 and K = 2, respectively.The SNR is set at 10 dB, and snapshots are varied from L = 200 to L = 2000.The figure shows that the proposed method performs slightly better than the other methods.However, the proposed method outperforms these methods in terms of computational complexity as discussed earlier.
Figure 4 shows the comparison of the scatter plot for the proposed method with novel L-shaped method in [12] and cumulant-based method in [11].The number of sources, antenna elements, and snapshots considered here are

Conclusions
A new 2D DOA estimation method has been proposed that solves the pair-matching problem for elevation and azimuth angles for multiple sources, with significantly reduced computational load.Furthermore, the proposed method does not require high complexity spectral search methods.The proposed method employs propagator Method (PM) in conjunction with parallel factor (PARAFAC) model.A new antenna array configuration consisting of two uniform linear arrays is also proposed to estimate 2D DOA for multiple sources.Simulation results demonstrate the effectiveness and better performance of the proposed method compared to existing methods.

Fig. 1 .
Fig. 1.Parallel-shaped array geometry.where s(t) is a K  1 dimensional vector representing received signals from K sources and n a (t), n b (t), and n c (t) are N  1 zero mean additive white Gaussian noise (AWGN) vectors with variance σ 2 .The elevation and azimuth angles are represented by  and , respectively.The matrix A() has dimensions N  K and is defined as: respectively, where the index k denotes the k th source.The elevation angle   k and the azimuth angle  k  for the k th source are then estimated as:

Tab 1 .
Complexity analysis.is not diagonal in the presence of coherent sources.