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

Authors: Zhongyuan Que; Benzhou Jin; Jianfeng Li
Journal: Wireless Communications and Mobile Computing
 (2021)

---

## Abstract

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.


---

## Body

## 1. 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–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–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 PARAFAC 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 procedureThe outline of this paper is given as follows. We discuss the data model for uniform planar array and introduce the PARAFAC model briefly in Section2. 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. Dm⋅ denotes a diagonal matrix whose diagonal elements are defined with the m-th row of the matrix. angle⋅ denotes phase angle operator. ⋅2 and ⋅F denote the ℓ2 and Frobenius norms. ⋅−1 and ⋅+ stand for the inverse and pseudo-inverse of a matrix.


## 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. Dm⋅ denotes a diagonal matrix whose diagonal elements are defined with the m-th row of the matrix. angle⋅ denotes phase angle operator. ⋅2 and ⋅F denote the ℓ2 and Frobenius norms. ⋅−1 and ⋅+ stand for the inverse and pseudo-inverse of a matrix.


## 2. Data Model

Consider a uniform rectangular array (URA) containingN×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, dx=dy=λ/2.Figure 1
Geometry of URA withN×M sensors. The circle represents a signal source at far field, and the solid points represent the sensors.Assume thatK uncorrelated far-field signals individually impinge on the array from θk,ϕkk=1,2,⋯,K, 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]:
(1)X~=AS+N,where X~∈ℂNM×L is the output data with noise; L denotes the number of snapshots; S=s1,s2,⋯,sKT∈ℂ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]
(2)A=ayυ1⊗axu1,ayυ2⊗axu2,⋯,ayυK⊗axuK,where uk=sinθkcosϕk and υk=sinθksinϕk; axuk and ayυk are the steering vectors of the array, which can be represented as [30]
(3)axuk=1,exp−j2πdxukλ,⋯,exp−j2πN−1dxukλT,ayυk=1,exp−j2πdyυkλ,⋯,exp−j2πM−1dyυkλT.More compactly, (2) can be also written as
(4)A=AxD1AyAxD2Ay⋮AxDMAy,where Ax=axu1,axu2,⋯,axuK and Ay=ayυ1,ayυ2,⋯,ayυ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 1 (Outer product [19]).
The outer product of three vectors,a∈ℂM×1, b∈ℂN×1, and c∈ℂL×1, denoted by a∘b∘c, is a M×N×L tensor whose elements are defined by a∘b∘cm,n,l=ambncl.Definition 2 (PARAFAC [19]).
The PARAFAC model is also known as the trilinear decomposition model. A canonical PARAFAC decomposition of a three-order tensorX∈ℂM×N×L can be expressed as
(5)X=∑f=1Faf∘bf∘cf,where af, bf, and cf 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 Xm∈ℂN×L, Xn∈ℂL×M, and Xl∈ℂM×N with the element Xmn,l=Xnl,m=Xlm,n=Xm,n,l. Then, Xm=BDmACT 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]
(6)Xm=AxDmAyS.Due to the symmetry of the PARAFAC model, the other two slice matrices can be obtained.(7)Yn=STDnAxAyT,(8)Zl=AyDlSTAxT.DefineX, Y, and Z as the results of the concatenation of matrices Xm, Yn, and Zl, respectively, and then, the noise-free received signal matrices X, Y, and Z can be represented as follows:
(9)X=X1X2⋮XM=AxD1AyAxD2Ay⋮AxDMAyS,Y=Y1Y2⋮YN=STD1AxSTD2Ax⋮STDNAxAyT,Z=Z1Z2⋮ZL=AyD1STAyD2ST⋮AyDLSTAxT.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].


## 3. 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 ofS, Ay, and Ax until they converge. 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]:
(10)R^=G^,H^,where G^∈ℂMN×K is the first column to the K-th column of R^ and H^∈ℂMN×MN−K stands for the remaining columns.Then, we can estimate the propagatorP by
(11)P^=G∧HG−1G∧HH,Define [33]
(12)P^c=IKP∧H=AxAxΦy⋮AxΦyM−1Ax−1,where Φy=diagexp−j2πdyυ1/λ,⋯,exp−j2πdyυK/λ and IK denotes a K-order identity matrix.The estimatesυ^k0 of υk can be obtained by partitioning the matrix P^c and eigenvalue decomposition.After reconstructing the matrixP^c, another matrix P^cs can be obtained by
(13)P^cs=AyAyΦx⋮AyΦxN−1Ax−1,where Φx=diagexp−j2πdxu1/λ,⋯,exp−j2πdxuK/λ. The estimates u^k0 of uk can also be obtained by a similar method.


### 3.2. 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 u^k0 and υ^k0 provided by PM are used to construct the matrices A∧x0 and A∧y0 as initial matrices.Recall that we assume noise is additive Gaussian noise, and it is reasonable to employ the least square principle to estimateS, Ax, and Ay. The estimation of the matrix S can be conducted by minimizing the following quadratic cost function [21]:
(14)minS∧nX~1X~2⋮X~M−A∧xn−1D1A∧yn−1A∧xn−1D2A∧yn−1⋮A∧xn−1DMA∧yn−1S∧nF,where X~m denotes the data matrix Xm with noise, m=1,2,⋯,M; A∧xn−1 and A∧yn−1 denote the estimates of Ax and Ay obtained from (n−1)-th iteration.Then, the LS estimate ofS can be obtained as [21]
(15)S∧n=A∧xn−1D1A∧yn−1A∧xn−1D2A∧yn−1⋮A∧xn−1DMA∧yn−1+X~1X~2⋮X~M.The LS fitting forAy is similar to S.
(16)minA∧ynY~1Y~2⋮Y~N−S∧TnD1A∧xn−1S∧TnD2A∧xn−1⋮S∧TnDNA∧xn−1A∧TnyF,where Y~n denotes the data matrix Yn with noise, n=1,2,⋯,N; S∧n denotes the estimate of S according to (15).Then, the LS estimate ofAy can be represented as
(17)A^yTn=S∧TnD1A∧xn−1S∧TnD2A∧xn−1⋮S∧TnDNA∧xn−1+Y~1Y~2⋮Y~N.Similarly, the LS fitting forAx is
(18)minA∧xnZ~1Z~2⋮Z~L−A∧ynD1S∧TnA∧ynD2S∧Tn⋮A∧ynDLS∧TnA∧TnxF,where Z~l denotes the data matrix Zl with noise, l=1,2,⋯,L, and A∧yn is the estimate of A^y according to (17).The estimate ofAx can be expressed as
(19)A∧xTn=A∧ynD1S∧TnA∧ynD2S∧Tn⋮A∧ynDLS∧Tn+Z~1Z~2⋮Z~L.According to (15), (17), and (19), we can repeatedly update the estimates of S, Ay, and Ax until convergence. Because of the utilization of PM, the proposed algorithm fast converges to the final estimates of S, Ax, and Ay, noted as S^f, A^fx, and A^fy. 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.


### 3.3. DOA Estimation

First, we normalize the column vectors ofA^fx and A^fy and make the first element of the column to equal one. Then, compute the phase vector rx by
(20)rx=−angleaxk=0,2πdxλ,⋯,2πN−1dxλTuk,=Bxuk,where axk denotes the k-th column vector of A^fx after normalization.According to LS criterion, calculate the estimates ofuk by
(21)u^k=Bx+rx.In a similar way, we can also get the estimatesυ^k of υk by the following expressions:
(22)υ^k=By+ry,where ry is another phase vector defined as
(23)ry=−angleayk=0,2πdyλ,⋯,2πM−1dyλTυk,=Byυk,where ayk denotes the k-th column vector of A^fy after normalization.Finally, the estimates ofθk and ϕk can be calculated by
(24)θ^k=arcsinu^k+jυ^k,(25)ϕ^k=angleu^k+jυ^k,where ⋅ denotes the modulus of the complex number and θ^k and ϕ^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 ofS^f, A^fx, and A^fy during the TALS decomposition, so the final estimates, θ^k and ϕ^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 estimatesu^k0 and υ^k0 of uk, υk.Step 2.
According to the PARAFAC models (6)–(8), reshape the received signal data to acquire the data matrices X~, Y~, and Z~.Step 3.
Construct the direction matricesA^x and A^y with u^k0 and υ^k0 and use them as initial matrices.Step 4.
According to (15), (17), and (19), update the estimates of S, Ay, and Ax alternately from the data matrices X~, Y~, and Z~ until convergence.Step 5.
According to (20)-(25), calculate the DOA estimates of separated signals, θ^k and ϕ^k.


## 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]:
(10)R^=G^,H^,where G^∈ℂMN×K is the first column to the K-th column of R^ and H^∈ℂMN×MN−K stands for the remaining columns.Then, we can estimate the propagatorP by
(11)P^=G∧HG−1G∧HH,Define [33]
(12)P^c=IKP∧H=AxAxΦy⋮AxΦyM−1Ax−1,where Φy=diagexp−j2πdyυ1/λ,⋯,exp−j2πdyυK/λ and IK denotes a K-order identity matrix.The estimatesυ^k0 of υk can be obtained by partitioning the matrix P^c and eigenvalue decomposition.After reconstructing the matrixP^c, another matrix P^cs can be obtained by
(13)P^cs=AyAyΦx⋮AyΦxN−1Ax−1,where Φx=diagexp−j2πdxu1/λ,⋯,exp−j2πdxuK/λ. The estimates u^k0 of uk can also be obtained by a similar method.


## 3.2. 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 u^k0 and υ^k0 provided by PM are used to construct the matrices A∧x0 and A∧y0 as initial matrices.Recall that we assume noise is additive Gaussian noise, and it is reasonable to employ the least square principle to estimateS, Ax, and Ay. The estimation of the matrix S can be conducted by minimizing the following quadratic cost function [21]:
(14)minS∧nX~1X~2⋮X~M−A∧xn−1D1A∧yn−1A∧xn−1D2A∧yn−1⋮A∧xn−1DMA∧yn−1S∧nF,where X~m denotes the data matrix Xm with noise, m=1,2,⋯,M; A∧xn−1 and A∧yn−1 denote the estimates of Ax and Ay obtained from (n−1)-th iteration.Then, the LS estimate ofS can be obtained as [21]
(15)S∧n=A∧xn−1D1A∧yn−1A∧xn−1D2A∧yn−1⋮A∧xn−1DMA∧yn−1+X~1X~2⋮X~M.The LS fitting forAy is similar to S.
(16)minA∧ynY~1Y~2⋮Y~N−S∧TnD1A∧xn−1S∧TnD2A∧xn−1⋮S∧TnDNA∧xn−1A∧TnyF,where Y~n denotes the data matrix Yn with noise, n=1,2,⋯,N; S∧n denotes the estimate of S according to (15).Then, the LS estimate ofAy can be represented as
(17)A^yTn=S∧TnD1A∧xn−1S∧TnD2A∧xn−1⋮S∧TnDNA∧xn−1+Y~1Y~2⋮Y~N.Similarly, the LS fitting forAx is
(18)minA∧xnZ~1Z~2⋮Z~L−A∧ynD1S∧TnA∧ynD2S∧Tn⋮A∧ynDLS∧TnA∧TnxF,where Z~l denotes the data matrix Zl with noise, l=1,2,⋯,L, and A∧yn is the estimate of A^y according to (17).The estimate ofAx can be expressed as
(19)A∧xTn=A∧ynD1S∧TnA∧ynD2S∧Tn⋮A∧ynDLS∧Tn+Z~1Z~2⋮Z~L.According to (15), (17), and (19), we can repeatedly update the estimates of S, Ay, and Ax until convergence. Because of the utilization of PM, the proposed algorithm fast converges to the final estimates of S, Ax, and Ay, noted as S^f, A^fx, and A^fy. 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.


## 3.3. DOA Estimation

First, we normalize the column vectors ofA^fx and A^fy and make the first element of the column to equal one. Then, compute the phase vector rx by
(20)rx=−angleaxk=0,2πdxλ,⋯,2πN−1dxλTuk,=Bxuk,where axk denotes the k-th column vector of A^fx after normalization.According to LS criterion, calculate the estimates ofuk by
(21)u^k=Bx+rx.In a similar way, we can also get the estimatesυ^k of υk by the following expressions:
(22)υ^k=By+ry,where ry is another phase vector defined as
(23)ry=−angleayk=0,2πdyλ,⋯,2πM−1dyλTυk,=Byυk,where ayk denotes the k-th column vector of A^fy after normalization.Finally, the estimates ofθk and ϕk can be calculated by
(24)θ^k=arcsinu^k+jυ^k,(25)ϕ^k=angleu^k+jυ^k,where ⋅ denotes the modulus of the complex number and θ^k and ϕ^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 ofS^f, A^fx, and A^fy during the TALS decomposition, so the final estimates, θ^k and ϕ^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 estimatesu^k0 and υ^k0 of uk, υk.Step 2.
According to the PARAFAC models (6)–(8), reshape the received signal data to acquire the data matrices X~, Y~, and Z~.Step 3.
Construct the direction matricesA^x and A^y with u^k0 and υ^k0 and use them as initial matrices.Step 4.
According to (15), (17), and (19), update the estimates of S, Ay, and Ax alternately from the data matrices X~, Y~, and Z~ until convergence.Step 5.
According to (20)-(25), calculate the DOA estimates of separated signals, θ^k and ϕ^k.


## 4. 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 isO5K3+2K2L+3K2NM−1+3K2MN−1+NM−KKL. 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 O3K3+3NMKL+2K2NM+NL+ML. The number of iterations is affected by many factors such as array size, iteration accuracy, and signal type. Define the number of iterations is T, and the total computational complexity is OT3K3+3NMKL+2K2NM+NL+ML. Therefore, the complexity of the proposed method is O5K3+2K2L+3K2NM−1+3K2MN−1+NM−KKL+T3K3+3NMKL+2K2NM+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


## 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 isO5K3+2K2L+3K2NM−1+3K2MN−1+NM−KKL. 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 O3K3+3NMKL+2K2NM+NL+ML. The number of iterations is affected by many factors such as array size, iteration accuracy, and signal type. Define the number of iterations is T, and the total computational complexity is OT3K3+3NMKL+2K2NM+NL+ML. Therefore, the complexity of the proposed method is O5K3+2K2L+3K2NM−1+3K2MN−1+NM−KKL+T3K3+3NMKL+2K2NM+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


## 5. Simulation Results

In this section, we employ a URA equipped with8×8 sensors to illustrate the improvement of the performance of 2D DOA estimation and signal separation of the proposed algorithm.Suppose there areK=3 typical modulated signals impinging on the array simultaneously, which are single-frequency signal s1t=cos2π×5×106t, linear frequency modulated signal s2t=cosπ×1012t2+2π×2×106t, and amplitude modulated signal s3t=cos2π×3×105tsin2π×5×106t. The DOAs of the signals are θ1,ϕ1=10°,15°θ2,ϕ2=20°,25°, and θ3,ϕ3=30°,35°, and the sampling frequency is 100 MHz. The noiseless source signal waveforms are demonstrated in Figure 2.Figure 2
The source signal waveform. (L=800).To assess the performance of DOA estimation, root mean square error (RMSE) is used,(26)RMSE=1CK∑c=1C∑k=1Kα∧k,c−αk2,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 α^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=10log10XF2/NF2, where X is noiseless 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,(27)ρ=1CK∑c=1C∑k=1Ksks^k,cHsk22s^k,c22,where sk is the k-th source signal and s^k,c is the corresponding estimate in the c-th trial.


### 5.1. Convergence Analysis

In Figure3, 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.Figure 3
Analysis of mean CPU time and mean number of iterations versus SNR.Define the sum of squared residualsSSR=∑n=1N∑m=1M∑l=1LX~n,m,l−∑k=1KA∧xn,kA∧ym,kS∧k,l2 and DSSR=SSRi−SSRi−1/SSRi−1, where SSRi 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.Figure 4
Analysis of convergence rate.As mentioned in Section4, 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.


### 5.2. 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 matrixA in (1) and compute the LS estimate of S by S^=A+X~.As shown in Figure5, 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 5
RMSE performance of different algorithms versus SNR (L=800).
(a)
Performance of elevation angle estimation(b)
Performance of azimuth angle estimationFigure6 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 6
Performance of signal separation of different algorithms versus SNR (L=800).Figure 7
Separated signals by the proposed algorithm (L=800, SNR=5 dB). Red line denotes the separated signal, and black line denotes the estimation error.Figure8 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.Figure 8
RMSE performance with different values ofL versus SNR.
(a)
Performance of elevation angle estimation(b)
Performance of azimuth angle estimationFigure9 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.Figure 9
Performance of signal separation with different values ofL versus SNR.


## 5.1. Convergence Analysis

In Figure3, 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.Figure 3
Analysis of mean CPU time and mean number of iterations versus SNR.Define the sum of squared residualsSSR=∑n=1N∑m=1M∑l=1LX~n,m,l−∑k=1KA∧xn,kA∧ym,kS∧k,l2 and DSSR=SSRi−SSRi−1/SSRi−1, where SSRi 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.Figure 4
Analysis of convergence rate.As mentioned in Section4, 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.


## 5.2. 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 matrixA in (1) and compute the LS estimate of S by S^=A+X~.As shown in Figure5, 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 5
RMSE performance of different algorithms versus SNR (L=800).
(a)
Performance of elevation angle estimation(b)
Performance of azimuth angle estimationFigure6 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 6
Performance of signal separation of different algorithms versus SNR (L=800).Figure 7
Separated signals by the proposed algorithm (L=800, SNR=5 dB). Red line denotes the separated signal, and black line denotes the estimation error.Figure8 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.Figure 8
RMSE performance with different values ofL versus SNR.
(a)
Performance of elevation angle estimation(b)
Performance of azimuth angle estimationFigure9 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.Figure 9
Performance of signal separation with different values ofL versus SNR.


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



---

*Source: 9963653-2021-09-09.xml*