Joint Processing of DOA Estimation and Signal Separation for Planar Array Using Fast-PARAFAC Decomposition

A joint processing of direction of arrival (DOA) and signal separation for planar array is proposed in this paper. Through sensor array processing theory, the output data of a planar array can be reconstructed as a parallel factor (PARAFAC) model, which can be decomposed with the trilinear alternating least square (TALS) algorithm. Aiming at the problem of slow speed on convergence for the standard PARAFAC method, we introduce the propagator method (PM) to accelerate the convergence of the TALS method and propose a novel method to jointly separate signals and estimate the corresponding DOAs. Given the initial angle estimates with PM, the number of iterations of TALS can be reduced considerably. The experiments indicate that our method can carry out signal separation and DOA estimation for typical modulated signals well and remain the same performance as the standard PARAFAC method with lower computational complexity, which verifies that our algorithm is effective.


Introduction
Signal separation and direction of arrival (DOA) estimation are significant themes in signal processing and have been investigated in various engineering fields including wireless communication, navigation, radar, and sonar [1][2][3][4][5]. As fundamental issues for signal processing, they have sparked considerable attention of researchers for decades. These two problems involve multiple signals and sensors which receive a mixture of signals [6]. The goal of DOA estimation is to find the source signal location, while the signal separation is aimed at extracting desired source signals. Through the years, many classical methods have been developed to solve these problems. For DOA estimation, subspace-based methods like MUSIC and ESPRIT have been widely adopted [7,8]. The conventional nonparametric Fourier-based methods have also been further developed [9], and the emerging sparse reconstruction-based methods like orthogonal matching pursuit (OMP) and sparse Bayesian inference (SBI) are introduced into DOA estimation [10,11]. For signal separation, the researches focus primarily on blind source separation (BSS) methods, where independent component analysis (ICA) and Joint Approximative Diagonalization of Eigen matrix (JADE) method are the most famous among these methods [12,13] and have been widely applied in the separation of speech and medical signals.
Compared with conventional DOA methods, BSS methods do not require much waveform prior information and are capable of identifying the transmission parameters based on the mixture signals, which has aroused an amount of attention of researchers. Many researchers study to apply blind separation algorithms into array signal model and have made lots of works. A combined complex blind source separation DOA estimation and signal recovery method was proposed for uniform linear array (ULA) in [14], which obtains better performance by exploiting BSS to estimate the array manifold. In [15], a blind DOA estimation method based on the JADE algorithm was proposed, which introduces fourth-order cumulant and has great performance in multipath environment. In [16], the chaotic adaptive firework algorithm was applied for solving the problem of radar emitter mixed signal. In [17], a new EM-based method for broadband DOA estimation and BSS was proposed, which reduces the complexity of traditional methods. In [18], a method based on eigenvalue decomposition for DOA estimation and blind separation of narrow-band independent signals was presented. The above studies were discussed for ULA geometry and did not involve more complex array structures like planar arrays which are more practical in actual applications.
In recent years, tensor technique has taken off in the field of data analysis and signal processing [19], in which trilinear decomposition or the parallel factor (PARAFAC) technique has been extensively investigated in radar and wireless communication fields especially [20][21][22][23][24]. PARAFAC is a common model for low-rank decomposition of a tensor, whose computation can be completed by alternating least squares (ALS). Many models in array signal processing can be represented as trilinear models, which enables us to utilize the trilinear alternating least square (TALS) algorithm to achieve parameter estimation and signal separation. The authors in [25,26] studied multiparameter estimation in bistatic multiple-input multiple-output (MIMO) radar and proposed a joint direction of departure (DOD) and direction of arrival (DOA) estimation using the PARAFAC model. In [27], the authors proposed a novel 2D-DOA estimation for trilinear decomposition-based monostatic cross MIMO radar. In [28], a joint DOA and carrier frequency estimation of narrow-band sources was proposed using the unitary PARAFAC method. These methods based on TALS can separate source signals and obtain automatically paired parameters without spectral peak search, but have relatively high computational complexity. It can be seen that PAR-AFAC has a great potential in DOA estimation and signal separation, but its shortcoming is the standard PARAFAC method has slow speed on convergence.
Motivated by the works mentioned above, in this paper, under the basic framework of PARAFAC, we propose a method of joint two-dimensional DOA estimation and signal separation for planar arrays. We first model the output of a planar array as the PARAFAC model and then utilize the propagator method (PM) to initialize the updated matrices in the TALS method, which effectively simplifies the complexity of the algorithm. Next, perform the TALS algorithm until convergence. Finally, acquire the 2D-DOA estimates and separated signals from the direction matrices and the source matrix which is estimated by TALS. The proposed method can achieve signal separation and DOA estimation for typical modulated signals well and remains the same performance as the standard PARAFAC method but with lower complexity. The experimental results verify the effectiveness of our algorithm.
We briefly summarize our main contributions as follows: (1) We model the output of the uniform rectangular array and reconstruct it into the PARAFAC model (2) We propose the fast-PARAFAC decomposition method for joint 2D-DOA estimation and signal separation, which utilizes PM to initialize the updated matrices in TALS and accelerate convergence (3) The proposed method has better performance of DOA estimation than 2D-PM and 2D-ESPRIT and can accurately separate the source signals with lower complexity compared with the standard PARAFAC approach (4) The proposed method can obtain separated signals and corresponding DOA estimates without an additional pairing procedure The outline of this paper is given as follows. We discuss the data model for uniform planar array and introduce the PARAFAC model briefly in Section 2. In Section 3, the proposed algorithm is described in detail. In Section 4, the complexity analysis and advantages of the proposed method are provided. The results of numerical simulations are given in Section 5, and conclusions are drawn in Section 6.
1.1. Notation. Lower-case and upper-case boldface letters denote vectors and matrices. ℂ denotes the sets of complex numbers. The superscripts ð⋅Þ T , ð⋅Þ * , and ð⋅Þ H represent the transpose, complex conjugate, and conjugate transpose of a vector or matrix, respectively. diag ð⋅Þ denotes a diagonal matrix that consists of the elements of the matrix. ⊗ denotes the Kronecker product. D m ð⋅Þ denotes a diagonal matrix whose diagonal elements are defined with the m-th row of the matrix. angleð⋅Þ denotes phase angle operator. k⋅k 2 and k⋅k F denote the ℓ 2 and Frobenius norms. ð⋅Þ −1 and ð⋅Þ + stand for the inverse and pseudo-inverse of a matrix.

Data Model
Consider a uniform rectangular array (URA) containing N × M sensors as depicted in Figure 1, where N and M are the numbers of elements along the x-axis and y-axis. The interelement spacings along both the x-axis and y-axis of the array are taken as half the wavelength of the waves, d x = d y = λ/2.
Assume that K uncorrelated far-field signals individually impinge on the array from fðθ k , ϕ k Þjk = 1, 2, ⋯, Kg, where θ k and ϕ k are the corresponding elevation and azimuth angles of the k-th signal (K < N × M, θ k ∈ ð0, 90°Þ, and ϕ k ∈ ð0, 180°Þ). The output of the rectangular array can be represented as follows [29]: whereX ∈ ℂ NM×L is the output data with noise; L denotes the number of snapshots; S = ½s 1 , s 2 , ⋯, s K T ∈ ℂ K×L is the signal matrix of L snapshots; N ∈ ℂ NM×L is the additive white Gaussian noise matrix. The array manifold matrix A ∈ ℂ NM×K consists of the steering vectors and is given by [29] where u k = sin θ k cos ϕ k and υ k = sin θ k sin ϕ k ; a x ðu k Þ and a y ðυ k Þ are the steering vectors of the array, which can be represented as [30] 2 Wireless Communications and Mobile Computing More compactly, (2) can be also written as where A x = ½a x ðu 1 Þ, a x ðu 2 Þ, ⋯, a x ðu K Þ and A y = ½a y ðυ 1 Þ, a y ðυ 2 Þ, ⋯, a y ðυ K Þ. Before further describing the data model, we need to introduce the definition of the tensor outer product and the parallel factor (PARAFAC) model [22].
Definition 2 (PARAFAC [19]). The PARAFAC model is also known as the trilinear decomposition model. A canonical PARAFAC decomposition of a three-order tensor X ∈ ℂ M×N×L can be expressed as where a f , b f , and c f stand for the f -th columns of matrices A ∈ ℂ M×F , B ∈ ℂ N×F , and C ∈ ℂ L×F . For a 3-way tensor X, define its sliced matrices X m ∈ ℂ N×L , X n ∈ ℂ L×M , and X l ∈ ℂ M×N with the element X mðn,lÞ = X nðl,mÞ = X lðm,nÞ = X m,n,l . Then, X m = BD m ðAÞC T can be viewed as "slicing" the 3-D array in a series of "slabs" (2-D arrays) and similarly for others [21].
Based on the PARAFAC model, the noiseless received data for URA can be written as [22] Due to the symmetry of the PARAFAC model, the other two slice matrices can be obtained.
Define X, Y, and Z as the results of the concatenation of matrices X m , Y n , and Z l , respectively, and then, the noisefree received signal matrices X, Y, and Z can be represented as follows: Note that in this paper, we assume that there is no mutual coupling across the sensors. In fact, mutual coupling will degrade the performance of the algorithms. The recent researches in the presence of mutual coupling can be found in [31].

The Proposed Algorithm
We show how to perform signal separation and DOA estimation using our proposed algorithm in this section. The standard PARAFAC suffers from expensive computation cost due to slow convergence. To handle this problem, we introduce the propagator method (PM) to accelerate TALS by providing the initial angle estimates. Then, alternately update the LS estimates of S, A y , and A x until they converge.

Wireless Communications and Mobile Computing
Finally, obtain the separated signals and the corresponding DOA estimates.
Note that in practice, we need to estimate the number of sources from the received signal first. In this study, we assume that the number of sources is known in advance.
3.1. Initialization with Propagator Method. By exploiting the property of rotational invariance, the propagator method can achieve the angle estimation with relatively low complexity [32,33].
First, compute the data covariance matrixR using the received signal data in (1) and partition it as follows [33]: whereĜ ∈ ℂ MN×K is the first column to the K-th column of R andĤ ∈ ℂ MN×ðMN−KÞ stands for the remaining columns. Then, we can estimate the propagator P bŷ Define [33] where Φ y = diag fexp ð−j2πd y υ 1 /λÞ, ⋯, exp ð−j2πd y υ K /λÞg and I K denotes a K-order identity matrix. The estimates b υ k0 of υ k can be obtained by partitioning the matrixP c and eigenvalue decomposition.
After reconstructing the matrixP c , another matrixP cs can be obtained bŷ where Φ x = diag fexp ð−j2πd x u 1 /λÞ, ⋯, exp ð−j2πd x u K /λÞg. The estimatesû k0 of u k can also be obtained by a similar method.

Trilinear
Alternating Least Square. Trilinear alternating least square (TALS) is the most common method for trilinear model decomposition [21,22]. The standard TALS algorithm utilizes random matrices as the initial load matrices, which usually converges slowly. In this part, the initial estimatesû k0 and b υ k0 provided by PM are used to construct the matrices A∧ x ð0Þ and A∧ y ð0Þ as initial matrices.
Recall that we assume noise is additive Gaussian noise, and it is reasonable to employ the least square principle to estimate S, A x , and A y . The estimation of the matrix S can be conducted by minimizing the following quadratic cost function [21]: : The LS fitting for A y is similar to S.
whereỸ n denotes the data matrix Y n with noise, n = 1, 2, ⋯, N; S∧ ðnÞ denotes the estimate of S according to (15). Then, the LS estimate of A y can be represented aŝ whereZ l denotes the data matrix Z l with noise, l = 1, 2, ⋯, L, and A∧ y ðnÞ is the estimate ofÂ y according to (17).
The estimate of A x can be expressed as : ð19Þ According to (15), (17), and (19), we can repeatedly update the estimates of S, A y , and A x until convergence. Because of the utilization of PM, the proposed algorithm fast converges to the final estimates of S, A x , and A y , noted asŜ f , A f x , andÂ f y . At this point, the task of signal separation is complete. The last part is to perform DOA estimation.
It is worth noting that the TALS algorithm outlined above contains only the simplest steps. Some techniques, like line search [34,35], can be coupling with the basic TALS algorithm, which may improve the rate of convergence further. There is no universally accepted most efficient TALS algorithm for all of the problems. We use the basic implementation of TALS for our issues and compare the complexity of our method with line search schemes [34,35] in Section 5.

DOA Estimation.
First, we normalize the column vectors ofÂ f x andÂ f y and make the first element of the column to equal one. Then, compute the phase vector r x by where a xk denotes the k-th column vector ofÂ f x after normalization. According to LS criterion, calculate the estimates of u k byû In a similar way, we can also get the estimates b υ k of υ k by the following expressions: where r y is another phase vector defined as where a yk denotes the k-th column vector ofÂ f y after normalization. Finally, the estimates of θ k and ϕ k can be calculated by where j⋅j denotes the modulus of the complex number and b θ k and b ϕ k are the estimates of the elevation and azimuth angles of the k-th signal. Note that there are the same permutation effects for the estimation ofŜ f ,Â f x , andÂ f y during the TALS decomposition, so the final estimates, b θ k and b ϕ k , are automatically paired.
3.4. The Procedure of the Proposed Algorithm. We summarize the major steps of our algorithm as follows: Step 1. Exploit the propagator method to calculate the initial estimatesû k0 and b υ k0 of u k , υ k .
Step 3. Construct the direction matricesÂ x andÂ y withû k0 and b υ k0 and use them as initial matrices.
Step 4. According to (15), (17), and (19), update the estimates of S, A y , and A x alternately from the data matrices X,Ỹ, andZ until convergence.

Performance Analysis
4.1. Complexity Analysis. Since complex multiplication requires the most computation time and resources, we use the time of complex multiplication to evaluate the complexity of the algorithm. The algorithm proposed in this paper adopts PM for initial estimation, whose complexity is Oð5K 3 + 2K 2 L + 3K 2 NðM − 1Þ + 3K 2 MðN − 1Þ + ðNM − KÞKLÞ. The complexity of the TALS method is related to the number of iterations and the complexity of a single iteration, and the complexity of each iteration is easily obtained as Oð3K 3 + 3NMKL + 2K 2 ðNM + NL + MLÞÞ. The number of iterations is affected by many factors such as array 5 Wireless Communications and Mobile Computing size, iteration accuracy, and signal type. Define the number of iterations is T, and the total computational complexity is OðT ð3K 3 + 3NMKL + 2K 2 ðNM + NL + MLÞÞÞ. Therefore, the complexity of the proposed method is Oð5K 3 + 2K 2 L + 3K 2 N ðM − 1Þ + 3K 2 MðN − 1Þ + ðNM − KÞKL + Tð3K 3 + 3NMKL + 2K 2 ðNM + NL + MLÞÞ. Due to PM initialization, the iterations of the proposed method are greatly reduced compared with the standard PARAFAC method, which we can see in the next section.

4.2.
Advantages. The advantages of the proposed algorithm are as follows: (1) The proposed method has lower computational cost than the standard PARAFAC method due to introducing PM (2) The proposed method outperforms 2D-ESPRIT and 2D-PM in the aspect of angle estimation performance for planar array (3) The proposed method can obtain separated signals and corresponding DOA estimation without an additional pairing procedure
To assess the performance of DOA estimation, root mean square error (RMSE) is used,  where C is the total number of Monte-Carlo trials, α k is the true value of the elevation or azimuth angle of k-th signal, and b α k,c is the estimate of the angle α k in the c-th trial. For each simulation, we set C = 1000. The signal to noise ratio is defined by SNR = 10 log 10 ðkXk 2 F /kNk 2 F Þ, where X is noise-less received data matrix and N is zero-mean white Gaussian noise matrix. Besides, to qualify the performance of signal separation, the average similar coefficient between the source signal and the separated signal is adopted, where s k is the k-th source signal andŝ k,c is the corresponding estimate in the c-th trial.

Convergence Analysis.
In Figure 3, we give the result of the mean CPU time and mean number of iterations based on our algorithm and other PARAFAC algorithm, where L = 800. PARAFAC and PM-PARAFAC denote the standard PARAFAC without modification and our algorithm, respectively. Besides, we also compare our method with line search schemes as "LS-PARAFAC" [34] and "ELSCS-PARAFAC" [35]. Owing to the initialization with PM, the mean CPU time and the mean number of iterations required are reduced considerably. The standard PARAFAC requires 153.2 iterations on average, while the proposed algorithm requires 10.8 iterations, which means that the proposed algorithm is ten times faster than the standard PARAFAC. It can also be seen that LS-PARAFAC and ELSCS-PARAFAC are faster than the standard PARAFAC but slower than our algorithm. Define the sum of squared residuals SSR = ∑ N n=1 ∑ M m=1 ∑ L l=1 ½X n,m,l − ∑ K k=1 A∧ x ðn, kÞA∧ y ðm, kÞS∧ðk, lÞ 2 and DSSR = where SSR i is the SSR after the i-th iteration. Figure 4 shows typical curves of the evolution of DSSR, where SNR = 10 and L = 800. We can also observe that the proposed method has faster convergence in Figure 4. As mentioned in Section 4, the complexity of the TALS method is related to many factors, like the scale of the array, the signal-to-noise ratio, and the types of signals. Although a slight change in configures can cause a significant difference in the execution time, the proposed method can always reduce the computational cost compared with the standard TALS method.

Comparison of Performance.
To verify the improvement of the proposed algorithm, the proposed method is compared with 2D-PM, 2D-ESPRIT, and the standard PARAFAC method. Note that the original 2D-PM and 2D-ESPRIT do not have the capability of signal separation, so we use their results of DOA estimation to construct the matrix A in (1) and compute the LS estimate of S byŜ = A +X .
As shown in Figure 5, it is evident that our method has the same angle estimation performance as the standard PARAFAC method, which surpasses 2D-ESPRIT and 2D-PM. Figure 6 shows the performance of signal separation of different algorithms with different SNR. From Figure 6, the signal separation performance of our algorithm, standard PARAFAC method, and 2D-ESPRIT are approximately the same with different SNR, while the 2D-PM algorithm has a slightly weaker signal separation performance with low SNR. Figure 7 is the separated signal diagram by the proposed algorithm. From Figure 7, it can be seen that the error between the source signal and the separated signal is very small, which is consistent with the high average similarity coefficient observed in Figure 6. Figure 8 illustrates the DOA estimation performance as a function of SNR with different numbers of snapshots. It is seen from Figure 8 that when the number of snapshots L increases, angle estimation performance can be improved.   Wireless Communications and Mobile Computing Figure 9 presents the performance of signal separation as a function of SNR with different numbers of snapshots. As seen from Figure 9, the average similar coefficients with different numbers of snapshots are close and increase with increasing of SNR.

Conclusions
In this paper, a joint processing of direction of arrival estimation and signal separation for planar array based on the fast-PARAFAC model is proposed. We model the output of planar array as the PARAFAC model and combine PM with the TALS method for DOA estimation and signal separation. The angle estimates by PM are used for the initialization of the TALS method. Then, the TALS method is used to separate the source signal and accurately estimate DOA. The proposed method not only inherits the advantages of the TALS method in signal separation but also takes advantage of the low complexity of the PM algorithm, which greatly reduces the total iterations of the standard PARAFAC algorithm. The results show that, as compared with the conventional DOA estimation approaches such as the 2D-PM and 2D-ESPRIT algorithm, the proposed method has better performance in signal separation and DOA estimation for planar array.

Data Availability
The data used to support the findings of this study are available from the corresponding author upon request.