# DOA and Polarization Parameters Estimation by Exploiting Canonical Polyadic Decomposition of Tensors

Authors: Long Liu; Ling Wang; Jian Xie; Zhaolin Zhang
Journal: Wireless Communications and Mobile Computing
 (2019)

---

## Abstract

A new algorithm to estimate the direction of arrival (DOA) and polarization parameters of signals impinging on an array with electromagnetic (EM) vector-sensors is presented by exploiting the canonical polyadic decomposition (CPD) of tensors. In addition to spatial and temporal diversities, further information from the polarization domain is considered and used in this paper. Estimation errors of these parameters are evaluated by the Cramér-Rao lower bound (CRB) benchmark, in the presence of additive white Gaussian noise (AWGN). The superiority of the proposed algorithm is shown by comparing with the derivative algorithms of MUSIC and ESPRIT. In the proposed algorithm, the parameters can be estimated by virtue of the diversities of the spatial and polarization belonging to the factor matrices, rather than the conventional subspace which is the foundation of MUSIC and ESPRIT. Additionally, the classical CPD algorithm based on Alternating Least Squares (ALS) is introduced to verify the efficacy of the proposed CPD algorithm. Results demonstrate that when the number of snapshots is greater than 50, the proposed algorithm requires a smaller number of snapshots to achieve a high level of performance, compared against the subspace-based algorithms and the ALS-based algorithm. Furthermore, in the matter of the array with a small number of sensors, the discovered advantage concerning the Root Mean Square Error (RMSE) in estimating the DOA and the polarization state of the signal is noteworthy.


---

## Body

## 1. Introduction

Signal parameters estimation is of great significance in many applications such as satellite navigation, wireless communication, radar, and sonar. An array of multiple sensors placed in different spatial locations is used to estimate the signal parameters, which mainly include spatial and polarization parameters of the signal. The spatial parameter refers to the direction of arrival (DOA) of the signal, while the polarization auxiliary angle and the polarization phase difference of the signal can be regarded as the polarization parameters. Among signal parameters estimation algorithms based on matrix operations, some depend on exhaustive search, e.g., multiple signal classification (MUSIC) [1], whereas others do not, e.g., estimation of the signal parameters via rotational invariance techniques (ESPRIT) [2] and root-MUSIC [3]. One of the common features of the matrix-based algorithms is that the matrix can only reflect the two-dimensional diversity of the desired signal. Some methods exploiting the signal-subspace embodied by MUSIC, ESPRIT and their derivative algorithms have been introduced in [4–12], in which the tensor-based data model has been combined in [10–12]. Based on the Higher-Order Singular Value Decomposition (HOSVD) of the measurement tensor, two algorithms of parameters estimation combined with ESPRIT have been developed in [10]. Similarly, in [11, 12], several tensor MUSIC methods based on HOSVD have been derived. The existing research results show that the performance of the subspace-based tensor method is improved compared to the subspace-based matrix method, and the computational complexity is higher in the tensor method than in the matrix method but of the same order. In [13], a greedy algorithm called randomized multiple candidate iterative hard thresholding for DOA estimation has been proposed. A set of potential candidates using the iterative hard thresholding algorithm are first obtained, and then the best candidate can be selected based on the a priori knowledge of the distribution of the signal and noise matrices.The canonical polyadic decomposition (CPD) of the third-order tensor is a minimal decomposition into a sum of rank-1 tensors [14]. The uniqueness of CPD, which is often called essential uniqueness in engineering papers, has been extensively studied in the field of algebraic geometry. The property of uniqueness makes the CPD a basic tool for signal separation, and it is widely used in telecommunication, array processing, machine learning, etc. [15–20]. In [15], a novel multiuser receiver for joint symbol and channel estimation by capitalizing on a tensor modeling of the end-to-end system has been proposed. The multiple invariance sensor array processing has been linked to the uniqueness of CPD in [20]. It shows the uniqueness of single and multiple invariance ESPRIT, which stems from a deterministic decomposition of the third-order tensor signal model. Moreover, it provides the DOA of the source with less restrictions of signal stationarity than the aforementioned statistical approaches. In [21], a method of DOA estimation for seismic plane waves has been considered from a deterministic perspective using CPD. In addition to temporal and spatial information, the different propagation speed of waves is taken into account by using the multidimensional feature of the tensor. In [22], a CPD-based approach for distinguishing the signals with the same DOA and copolarized state has been proposed. However, the computational complexity of the algorithm is very high due to the complexity of the problem model.The main purpose of this paper is to explore deterministic decomposition of the third-order tensor signal model by exploiting the uniqueness of CPD and then estimate the signal parameters. The polarization diversity of the signal is considered in this paper, in addition to the temporal and spatial sampling employed by the most current parameters estimation methods. Note that the proposed algorithm focus on not only the DOA estimation but also the estimation of polarization parameters delivered by the signal. The strength of our approach is its excellent performance compared against the traditional subspace-based algorithms such as MUSIC, HOSVD-MUSIC, and HOSVD-ESPRIT given shorter snapshots and the array with a small number of sensors.The rest of this paper is organized as follows. Section2 reviews the CPD prerequisites. The tensor signal model is derived in Section 3, and the introduction to traditional subspace-based algorithms is also provided. The CPD-based algorithm is proposed in Section 4, followed by the numerical simulations in Section 5. Finally, the last part presents our concluding remarks. The algebra notations involved in the paper are shown in Table 1.Table 1
Algebra notations.
a Scalara a Vectora A MatrixA A TensorA A∗ Conjugate ofA AT Transpose ofA AH Conjugate transpose ofA A⊚B Outer product betweenA and B A⊗B Kronecker product betweenA and B A⊙B Khatri-Rao product betweenA and B A∘B Hadamard product betweenA and B rA Rank of the matrixA R Real number fieldR C Complex number fieldC


## 2. Prerequisites for CPD

Let the(i,j,k)-th entry D(i,j,k) of the third-order tensor D∈CI×J×K be symbolized by dijk. Herein, we define rD=1; i.e., there exist three nonzero vectors a∈CI, b∈CJ, and c∈CK such that D=a⊚b⊚c, which means that dijk=aibjck stands for all possible values of indices.A third-order tensorF∈CI×J×K with general-rank can be expressed as a sum of rank-one items:(1)F=∑r=1RDr⇔F=∑r=1Rar⊚br⊚cr,where ar∈CI, br∈CJ, and cr∈CK. The decomposition shown in (1) is called the polyadic decomposition (PD). Particularly, if the number of rank-one terms R is minimal, then R is the rank of F (i.e., rF=R), and (1) is defined as the CPD of F.Additionally, (1) can also be written as F=[A,B,C]R, where the matrices A=[a1,…,aR]∈CI×R, B=[b1,…,bR]∈CJ×R, and C=[c1,…,cR]∈CK×R are defined as the first, second, and third factor matrix of F, respectively. Obviously, the entries of F can be written as fijk=∑r=1Rairbjrckr. In addition, here is a potential scenario: the CPD of a general-rank tensor F is not unique. Nevertheless, we still call that the CPD of F is unique when it is only subject to the trivial indeterminacies, which include that the rank-one terms D(r) shown in (1) can be arbitrarily permuted and the vectors belonging to any rank-one term can be arbitrarily scaled. Similarly, the factor matrices for any two CPDs F=[A,B,C]R and F=[A¯,B¯,C¯]R coincide up to column permutation and scaling; for instance, A¯ can be obtained by scaling and permuting the columns of A (see [9] for more details). Herein, for reshaping the tensor F into a matrix F, we introduce the following definition.Definition 1.
The matrix unfolding that yieldsN′∈{1,…,N} and the permutation T{1,…,N}={p1,…,pN}, pn∈{1,…,N} for the tensor A∈CM1×M2×…×MN is defined as AT,N′∈C(Mp1Mp2…MpN′)×(MpN′+1…MpN). The (l1,l2)-th entry of A is given by(2)AT,N′l1,l2=am1m2…mN,where 1⩽mn⩽Mn for n∈{1,2,…,N}, and l1=mpN′+∑n=2N′(mpn-1-1)∏j=nN′Mpj, l2=mpN+∑n=N′+2Nmpn-1-1∏j=nNMpj. Assuming that N′=2 and the permutation T1{1,2,3}={1,2,3}, the corresponding matrix unfolding of (1) can be expressed as(3)FT1,2=F1⋮FI=Bdiaga1CT⋮BdiagaICT=A⊙BCT,where F(T1,2)∈CIJ×K, Fi=[fijk]j,k=1J,K symbolizes the ith horizontal slice of F, ai denotes the ith row of A, and diag(·) is the diagonalization operator. What is more, in the case of N′=2, we can obtain the other matrix unfoldings:(4)FT2,2=B⊙CAT,FT3,2=C⊙ABT,etc,where T2{1,2,3}={2,3,1} and T3{1,2,3}={3,1,2}. This paper mainly uses the algebraic algorithms described in [14, 23] to explore the application of CPD in signal parameters estimation. The key of the algebraic algorithms presented in [14, 21] is the following theorem, which has been derived in [24].Theorem 2.
Define a third-order tensorF=[A,B,C]R, and assume that the factor matrices B and C have full column rank and any two columns selected from the factor matrix A are linearly independent:(5)rB=rC=R,kA⩾2,where kA represents the largest number such that every subset with kA columns of the factor matrix A is linearly independent. Then rF=R; the CPD of F can be computed by algebraic algorithms uniquely.One algebraic algorithm elaborated in [23] is adopted in this paper. Herein, we consider the particular case of F∈C2×R×R that conforms to Theorem 2. By (3), the matrix unfolding of F can be expressed as(6)FT1,2=F1F2=Bdiaga1CTBdiaga2CT,which obtains the equations that(7)F1F2-1=Bdiaga1diaga2-1B-1,F1TF2-T=Cdiaga1diaga2-1︸DC-1.SincekA⩾2, the diagonal elements of the diagonal matrix D are distinct. Therefore, the factor matrices B and C can be uniquely recovered by the eigenvalue decomposition (EVD) of (7), and then the matrix A can be easily identified from (6). However, the factor matrices are not necessarily square matrices in reality. Hence, we need to consider the generalized inverse operation of the matrices to obtain the CPD of the tensor.Herein, in order to avoid ambiguity in signal modeling, some of the main assumptions are stated as follows.A1. Completely polarized wave: the horizontal and vertical electric components located in the wavefront of the transverse electromagnetic (TEM) wave are completely correlated.A2. Far-field assumption: we assume that the distance between the source and the receiving array is much larger than the array aperture. Thus, the source of the signal can be considered point-like. Moreover, the wavefront of the TEM wave can be approximated as a plane at the sensor level.A3. Narrow-band assumption: the reciprocal of the signal bandwidth is much greater than the maximum propagation time of the signal passing through the array aperture. Typically, the signal received by an array is a varying complex envelope, which is modulated at a carrier with high frequency (HF). In order to facilitate the digital signal processing, the intermediate frequency (IF) or baseband signal is usually wanted. Hence, the analog down-conversion for the received HF signal is a more common practice.


## 3. Signal Modeling

This section mainly focuses on the three-dimensional feature, i.e., the vector-sensor enjoys in spatial, polarization, and time diversities and then establishes the third-order tensor-based signal model. In addition, conventional matrix-based solutions are introduced briefly. Firstly, the diversities involved in this paper are shown as follows.Time Diversity. The baseband analog signal obtained by down-conversion is a time-dependent function. Moreover, the discrete-time signal can be obtained by sampling the baseband analog signal in time domain. It is assumed that R signals with the complex amplitudes {sr(k),r=1,…,R} impinge on the vector-sensor array. Thus, the time diversity can be formulized as(8)f1k=∑r=1Rsrk,where k=1,…,K, represents the snapshot at the instant tk.Spatial Diversity. The spatial sampling of the array reflects the spatial diversity. Consider a vector-sensor array composed of L EM vector-sensors located at bl∈R3 for l=1,…,L. Combined with the time diversity, the signal received by the array can be expressed as(9)f2k,l=∑r=1Re-i2πblTϵr/λsrk,where the symbol λ denotes the wave length of the signal impinged on the vector-sensor array. Moreover, ϵr=-sin⁡ϕrcos⁡θrsin⁡ϕrsin⁡θrcos⁡ϕrT denotes the incident direction of the rth signal source associated with the azimuth angle θr∈[0,2π) and the elevation angle ϕr∈[0,π].Polarization Diversity. The EM vector-sensor is usually composed of multiple magnetic and electric short dipoles in a common spatial point. Hence, the EM vector-sensor is capable of measuring multiple electromagnetic components with the number of not less than one. For example, a three-dipole sensor can measure the three components of the electric field in three perpendicular directions. Furthermore, the complete EM vector-sensor has the ability to obtain all electromagnetic components (i.e., three electric components and three magnetic components). In particular, if the vector-sensor array is composed of a monopole, the orientations of the monopoles cannot be the same; otherwise, the array will degrade into a scalar-array. We formulate the polarization diversity combined with the time diversity,(10)f3k,j=∑r=1RajpΨrsrk,where j=1,…,J, J is the number of component measured by an EM vector-sensor, J=3 for three-dipoles, and J=6 for complete EM vector-sensor. Moreover, ajp(Ψr) represents the jth component of the polarization steering vector ap(Ψr), and the symbol Ψr=(θr,ϕr,γr,ηr) denotes the spatial-polarization parameter of the rth signal source with respect to the polarization auxiliary angle γr∈[0,π/2] and the polarization phase difference ηr∈[-π,π).


### 3.1. Tensor-Based Signal Model

Consider the diversities mentioned above, and assume that all the signals are cofrequency and incoherent. The spatial steering vector of the signal with DOA(θ,ϕ) is expressed as(11)asθ,ϕ=e-i2πb1Tϵ/λ,…,e-i2πbLTϵ/λT.The polarization steering vector of the signal with spatial-polarization parameterΨ=(θ,ϕ,γ,η) is given by(12)apΨ=BΞθ,ϕhγ,η,Ξθ,ϕ=-sin⁡θcos⁡ϕcos⁡θcos⁡θcos⁡ϕsin⁡θ0-sin⁡ϕcos⁡ϕcos⁡θsin⁡θcos⁡ϕsin⁡θ-cos⁡θ-sin⁡ϕ0,where hγ,η=[cos⁡γ,sin⁡γeiη]T and B∈RJ×6 represents the polarization selection matrix of the EM vector-sensor. Particularly, B=I6×6 represents the complete EM vector-sensor, and B=[I3×3,03×3] symbolizes the three-dipole sensor, where I denotes the identity matrix. Additionally, the first three rows of the matrix Ξθ,ϕ represent the projections of the electric field component of the incident signal with DOA (θ,ϕ) in the three axial directions of the coordinate system, and the last three rows of the matrix Ξθ,ϕ represent the projections of the magnetic field component in the three axial directions. The spatial-polarization steering vector can be further defined as(13)aΨ=asθ,ϕ⊗apΨ.Hereby, letn(k) denote the spatial white noise at the instant tk, which is additive and Gaussian complex circular. The vector-output of the array at the instant tk can be expressed as(14)yk=∑r=1RaΨrsrk+nk.Considering that the multidimensional nature of the tensor just caters to the multidomain diversities of the signal, the following spatial-polarization steering matrix can be obtained by performing the outer product operations of (12) and (11):(15)A~Ψ=apΨ⊚asθ,ϕ.where A~(Ψ)∈CJ×L. Let s=[s(1),…,s(K)] denote the consecutive K snapshots of the signal; the tensor-output of K snapshots can be modeled as(16)Y=∑r=1RA~Ψr⊚sr+N=∑r=1RapΨr⊚asθr,ϕr⊚sr+N,where Y∈CJ×L×K. The tensor noise is expressed as N∈CJ×L×K, which is yielded by the tensorization with respect to the noise matrix N=[n(1),…,n(k)].


### 3.2. Traditional Subspace-Based Algorithms

Subspace decomposition methods with respect to EVD and SVD are the mainstream of MUSIC, ESPRIT, and their derivative algorithms. Herein, we introduce the matrix-based MUSIC algorithm and the tensor-based MUSIC algorithm based on HOSVD. The covariance matrix forK snapshots corresponding to (14) can be given by(17)R^yy=1K∑k=1KykyHk,where R^yy∈CLJ×LJ. Then the eigenvector corresponding to the small eigenvalues of R^yy can be written in a matrix form U^n, where U^n∈CLJ×(LJ-R). Moreover, according to (11) and (12), we need to construct a search matrix:(18)Dθ,ϕ=asθ,ϕ⊗BΞθ,ϕ.The solution of the matrix-based MUSIC algorithm can be defined as(19)θr,ϕr=arg⁡maxθ,ϕ⁡det-1Dθ,ϕHU^nU^nHDθ,ϕ,where r=1,…,R. The estimation of the polarization parameter is detailed in Appendix A. For the tensor-based MUSIC algorithm based on HOSVD [11], the following signal model needs to be constructed:(20)Yk=∑r=1RA~Ψrsr+N,where Y(k)∈CJ×L. The covariance tensor for K snapshots is expressed as(21)R^YY=1K∑k=1KYk⊚Y∗k,where R^YY∈CJ×L×J×L. Then, R^YY is decomposed using the HOSVD procedure as(22)R^YY=K^×1U^1×2U^2×3U^3×4U^4,where ×n represents n-mode product, K^ is the core tensor, and U^(n) can be obtained by SVD of n-mode matrix unfolding of tensor R^YY. Finally the parameters of the sources are obtained by maximizing the following criterion:(23)θr,ϕrγr,ηr=arg⁡maxΨA~Ψ×1U^01U^01H×2U^02U^02HF-1,where U^0(n) is a truncated form of U^(n) (see [11] for specific truncation methods). In addition, the HOSVD-ESPRIT algorithm depending on translation invariant property of the array is elaborated in Appendix B.


## 3.1. Tensor-Based Signal Model

Consider the diversities mentioned above, and assume that all the signals are cofrequency and incoherent. The spatial steering vector of the signal with DOA(θ,ϕ) is expressed as(11)asθ,ϕ=e-i2πb1Tϵ/λ,…,e-i2πbLTϵ/λT.The polarization steering vector of the signal with spatial-polarization parameterΨ=(θ,ϕ,γ,η) is given by(12)apΨ=BΞθ,ϕhγ,η,Ξθ,ϕ=-sin⁡θcos⁡ϕcos⁡θcos⁡θcos⁡ϕsin⁡θ0-sin⁡ϕcos⁡ϕcos⁡θsin⁡θcos⁡ϕsin⁡θ-cos⁡θ-sin⁡ϕ0,where hγ,η=[cos⁡γ,sin⁡γeiη]T and B∈RJ×6 represents the polarization selection matrix of the EM vector-sensor. Particularly, B=I6×6 represents the complete EM vector-sensor, and B=[I3×3,03×3] symbolizes the three-dipole sensor, where I denotes the identity matrix. Additionally, the first three rows of the matrix Ξθ,ϕ represent the projections of the electric field component of the incident signal with DOA (θ,ϕ) in the three axial directions of the coordinate system, and the last three rows of the matrix Ξθ,ϕ represent the projections of the magnetic field component in the three axial directions. The spatial-polarization steering vector can be further defined as(13)aΨ=asθ,ϕ⊗apΨ.Hereby, letn(k) denote the spatial white noise at the instant tk, which is additive and Gaussian complex circular. The vector-output of the array at the instant tk can be expressed as(14)yk=∑r=1RaΨrsrk+nk.Considering that the multidimensional nature of the tensor just caters to the multidomain diversities of the signal, the following spatial-polarization steering matrix can be obtained by performing the outer product operations of (12) and (11):(15)A~Ψ=apΨ⊚asθ,ϕ.where A~(Ψ)∈CJ×L. Let s=[s(1),…,s(K)] denote the consecutive K snapshots of the signal; the tensor-output of K snapshots can be modeled as(16)Y=∑r=1RA~Ψr⊚sr+N=∑r=1RapΨr⊚asθr,ϕr⊚sr+N,where Y∈CJ×L×K. The tensor noise is expressed as N∈CJ×L×K, which is yielded by the tensorization with respect to the noise matrix N=[n(1),…,n(k)].


## 3.2. Traditional Subspace-Based Algorithms

Subspace decomposition methods with respect to EVD and SVD are the mainstream of MUSIC, ESPRIT, and their derivative algorithms. Herein, we introduce the matrix-based MUSIC algorithm and the tensor-based MUSIC algorithm based on HOSVD. The covariance matrix forK snapshots corresponding to (14) can be given by(17)R^yy=1K∑k=1KykyHk,where R^yy∈CLJ×LJ. Then the eigenvector corresponding to the small eigenvalues of R^yy can be written in a matrix form U^n, where U^n∈CLJ×(LJ-R). Moreover, according to (11) and (12), we need to construct a search matrix:(18)Dθ,ϕ=asθ,ϕ⊗BΞθ,ϕ.The solution of the matrix-based MUSIC algorithm can be defined as(19)θr,ϕr=arg⁡maxθ,ϕ⁡det-1Dθ,ϕHU^nU^nHDθ,ϕ,where r=1,…,R. The estimation of the polarization parameter is detailed in Appendix A. For the tensor-based MUSIC algorithm based on HOSVD [11], the following signal model needs to be constructed:(20)Yk=∑r=1RA~Ψrsr+N,where Y(k)∈CJ×L. The covariance tensor for K snapshots is expressed as(21)R^YY=1K∑k=1KYk⊚Y∗k,where R^YY∈CJ×L×J×L. Then, R^YY is decomposed using the HOSVD procedure as(22)R^YY=K^×1U^1×2U^2×3U^3×4U^4,where ×n represents n-mode product, K^ is the core tensor, and U^(n) can be obtained by SVD of n-mode matrix unfolding of tensor R^YY. Finally the parameters of the sources are obtained by maximizing the following criterion:(23)θr,ϕrγr,ηr=arg⁡maxΨA~Ψ×1U^01U^01H×2U^02U^02HF-1,where U^0(n) is a truncated form of U^(n) (see [11] for specific truncation methods). In addition, the HOSVD-ESPRIT algorithm depending on translation invariant property of the array is elaborated in Appendix B.


## 4. CPD-Based Algorithm

We hereby execute the CPD of the tensor-based signal model to obtain the factor matrices, which are then used to extract the parameters of the signals. LetY=[Ap,As,S]R, where Ap=[ap(Ψ1),…,ap(ΨR)]∈CJ×R, As=[as(θ1,ϕ1),…,as(θR,ϕR)]∈CL×R, and S=[s1,…,sR]∈CK×R. Assume that at least two factor matrices of Y have full column rank and the k-rank of the rest factor matrix is not less than 2, which are consistent with the Theorem 2. Moreover, there are conditions in which J, L, and K are not less than R, respectively, and this cannot be ignored here to ensure the uniqueness of the factor matrices.


### 4.1. Computation of CPD

As is well known, CPD of a third-order tensor can be achieved by the classical Alternating Least Squares (ALS) algorithm, which is detailed in AppendixC. Differently, we hereby propose an EVD-based algebraic method to achieve CPD of tensor. The tensor-based signal model (16) can be unfolded as(24)YT1,2=Y1⋮YJ=Asdiagap1ST⋮AsdiagapJST,where apj denotes the jth row of the Ap, j=1,…,J. According to (24), the following relationship can be easily derived:(25)Yj1Yj2†=Asdiagapj1diagapj2-1As†,Yj1TYj2T†=Sdiagapj1diagapj2-1S†,where j1≠j2, j1,j2=1,…,J, and (·)† denotes the Moore-Penrose inverse. From (24), we have Yj=Asdiag(apj)ST, and j=1,…,J. Hence, for any j, diag(apj)=As†Yj(ST)†. The features of the CPD are demonstrated as follows:(1)
The columns of the factor matrixAs are the eigenvectors of Yj1Yj2†. The corresponding eigenvalues are the first R larger eigenvalues, which are arranged in a descending order.(2)
The eigenvectors corresponding to theR-large eigenvalues of Yj1T(Yj2T)† constitute the factor matrix S.(3)
The rows of factor matrixAp can be recovered thanks to the relationship diag(apj)=As†Yj(ST)†.Obviously, whenJ>2, there are various combinations (j1,j2) that can be used for EVD. In order to make full use of the statistical properties of (24), we select CJ2 possible combinations for EVD to obtain Asj′, Sj′, where j′=1,…,CJ2. We follow the convention that the (j1,j2)-tuples are ordered lexicographically: the j1′th tuple (j1,j2) is preceding the j2′th tuple (j3,j4) if and only if j1⩽j3, j2⩽j4. According to the trivial indeterminacies aforementioned for CPD uniqueness, we know that Asj1′ and Asj2′, j1′≠j2′, coincide up to column permutation and scaling. In order to achieve the addition of multiple sets of results, we need to match the multiple sets of results. That is to say, the eigenvectors in each factor matrix obtained by the EVD need to be arranged in a descending order according to the corresponding eigenvalues, and all the eigenvectors need to be normalized. After the matching process, A~sj′ and S~j′ can be obtained, and then the following statistical operations are performed:(26)As=1CJ2∑j′=1CL2A~sj′,S=1CJ2∑j′=1CL2S~j′.Next, the factor matrixAp can be recovered according to the results of (26).


### 4.2. Parameters Estimation from the Factor Matrices

From the aforementioned tensor-based signal model, we can see that the spatial-polarization parameters of signals can be reflected by the factor matricesAp and the factor matrix As that only contains the spatial parameters of signals. Hence, once the factor matrices are obtained, we first need to extract the spatial parameters from As.


#### 4.2.1. Estimation of Spatial Parameters

Spatial DOA(θr,ϕr) of the rth signal impinging on the array can be calculated directly by solving the following equations:(27)As1r=e-i2πb1Tϵr/λ⋮AsLr=e-i2πbLTϵr/λObviously, any two equations belonging to (27) can be used to derive the DOA of the rth signal. Therefore, we can get CL2 sets of solutions from (27). Improve the accuracy of the estimated parameters by averaging multiple sets of solutions:(28)θr,ϕr=1CL2∑m=1CL2θr,ϕrm,where (θr,ϕr) represents the solution of the mth subequations {[As]l1r,[As]l2r}. We follow the convention that the (l1,l2)-tuples are ordered lexicographically: the mth tuple (l1,l2) is preceding the nth tuple (l1′,l2′) if and only if l1⩽l1′, l2⩽l2′. Furthermore, we can construct the optimization problem as follows to estimate the DOAs of the signals:(29)θr,ϕr=arg⁡minθ,ϕ⁡ζsr∘asθ,ϕ-as:,r2,where as:,r represents the rth column of the factor matrix As obtained by CPD, ζsr represents the amplitude vector of as:,r, and the operator · denotes the 2-norm. The solution method to (29) and its performance can be used as a further work, and this paper discusses only the performance based on the solution of (28).


#### 4.2.2. Estimation of Polarization Parameters

Once the estimation of DOA is obtained, the polarization parameters(γr,ηr) of the rth signal can be extracted from(30)ap:,r=BΞθr,ϕrcos⁡γr,sin⁡γreiηrT,where ap:,r represents the rth column of the factor matrix Ap obtained by CPD. Hereby, the estimated value of hγr,ηr can be expressed as h^γr,ηr=Ξθr,ϕr-1B-1ap:,r, and the estimated values of the polarization parameters can be obtained by(31)γr=arctan⁡h^γr,ηr2h^γr,ηr1,ηr=arg⁡h^γr,ηr2h^γr,ηr1,where arctan(·) is the arctangent operator and the operator arg(·) denotes phase acquisition.Obviously, the estimation accuracy of the DOA only depends on the factor matrixAs obtained by CPD, and the estimation accuracy of the polarization parameters is affected by both the estimation accuracy of the DOA and the matrix Ap obtained by CPD. It is well known that, in order to satisfy the validity of the parameters estimation, the algorithm must follow the premise of L>R. According to [14], there exists Kruskal bound ensuring the CPD uniqueness R<(J+L+K-2)/2. With the noise present, the accuracy of the CPD is severely impaired, when R is very close to the bound. However, as (J+L+K-2)/2 increases, the dependence of the CPD accuracy on (J+L+K-2)/2 will rapidly weaken. In addition, as can be seen from (26) and (28), the proposed algorithm can obtain multiple sets of parameters estimation results in each CPD and perform statistical operations on these results. However, what is different is that the subspace-based methods requires sufficient statistical information, which is usually obtained by increasing the number of sensors or snapshots, to ensure the accuracy of the parameters estimation, even if the number of the parameters to be estimated is very small. We, therefore, conclude that the proposed algorithm can achieve relatively desirable performance by using a smaller number of snapshots as well as a smaller number of sensors compared with the subspace-based methods. The subsequent simulation results also justify the above theoretical analysis.


### 4.3. Computational Complexity

The parameters estimation algorithm using EVD-based CPD of tensors is summarized in Algorithm1, which is easy to implement by row-by-row access. Note that the computational complexity scales linearly to the number of loops, since the procedure of the algorithm runs serially.Algorithm 1:EVD-based CPD for parameters estimation.
Require:  Y∈CJ×L×K, L, J, K, R. 1: Reshape the matrciesY1,…,YJ∈CL×K by unfolding the tensor Y; 2: Forj′=1 to CJ2 do 3:Calculate Yj1Yj2†, Yj1T(Yj2T)†; 4:Obtain Asj′ by executing the EVD of Yj1Yj2†; 5:Obtain Sj′ by executing the EVD of Yj1T(Yj2T)†; 6:Obtain A~sj′ and S~j′ by matching processing. 7: End for 8:As=1CJ2∑j′=1CL2A~sj′,S=1CJ2∑j′=1CL2S~j′; 9: RecoverAp according to diag(apj)=As†Yj(ST)†; 10: Forr=1 to R do 11:For m=1 to CL2 do 12:Calculate (θr,ϕr)m by solving equations {[As]l1r,[As]l2r}; 13:End for 14:θr,ϕr=1CL2∑m=1CL2(θr,ϕr)m; 15:h^γr,ηr=Ξθr,ϕr-1B-1ap:,r; 16:γr=arctanh^γr,ηr2h^γr,ηr1,ηr=argh^γr,ηr2h^γr,ηr1. 17: End for Ensure:  Ψr=(θr,ϕr,γr,ηr), r=1,…,R.In fact, the EVDs of (25) follow from the generalized EVDs (GEVDs) of the matrix pencils (Yj1,Yj2) and (Yj1T,Yj2T). In order to evaluate the complexity of the proposed algorithm, we first need to know the number of multiplications required by GEVD. An effective GEVD solution is given by the method of QZ iteration [25, pp.413-414], in which pLKR multiplications are required for two L×K matrices with rank-R, where the constant p depends on the selected iterative strategy, e.g., for Hessenberg-Triangular Reduction, p=8. The Moore-Penrose inverse of the M×N matrix or N×M matrix involves 4MN2+8N3 multiplications, where M>N. Additional J×(RLK+R2K) multiplications are also required to recover the factor matrix Ap. Therefore, the computational complexity of the EVD-based CPD does not exceed (pCJ2+8+J)K·[O(L2)], since the premise of L,J,K>R has been established. Furthermore, in order to compare the computational complexity for the proposed algorithm, we analyze the required number of multiplications for subspace-based methods. For the aforementioned matrix-based method, the computational complexity of obtaining subspaces is approximately (L+K)·[O(L2)]. Considering the operation process of HOSVD based on orthogonal iterations [25, pp.454-455], the computational complexity of subspace for tensor model is about (KJ2+8J2+pJ)·[O(L2)]. Obviously, the conclusion that the computational complexity of the CPD algorithm is lower than that of the tensor-based HOSVD methods can be obtained through the above analysis. Additionally, for the extraction process of the parameters, the computational complexity of the search strategies is generally much higher than that of the equation-solving methods, under the premise of ensuring the estimation accuracy. Therefore, the computational complexity of the proposed algorithm is lower than that of the search-based methods.Additionally, according to the procedure of the ALS-based CPD described in AppendixC, the computational complexity of the factor matrix S is not less than (5+L+K)·[O(R3)], since the premise of L,J,K>R has been established. Similarly, the computational complexity of the factor matrices Ap and As is not less than (5+2K)·[O(R3)]. Therefore, the computational complexity of the ALS-based CPD is about (15+L+5K)G·[O(R3)], where G represents the number of iterations. We use the Matlab 2016a to implement the after-mentioned simulations on a computer with Intel Core i5 CPU6200U 2.3GHz and 8GB memory running Windows 10. From the perspective of computing time, the time required for ALS-based CPD is about twice that of EVD-based CPD to obtain the subsequent simulation results.


## 4.1. Computation of CPD

As is well known, CPD of a third-order tensor can be achieved by the classical Alternating Least Squares (ALS) algorithm, which is detailed in AppendixC. Differently, we hereby propose an EVD-based algebraic method to achieve CPD of tensor. The tensor-based signal model (16) can be unfolded as(24)YT1,2=Y1⋮YJ=Asdiagap1ST⋮AsdiagapJST,where apj denotes the jth row of the Ap, j=1,…,J. According to (24), the following relationship can be easily derived:(25)Yj1Yj2†=Asdiagapj1diagapj2-1As†,Yj1TYj2T†=Sdiagapj1diagapj2-1S†,where j1≠j2, j1,j2=1,…,J, and (·)† denotes the Moore-Penrose inverse. From (24), we have Yj=Asdiag(apj)ST, and j=1,…,J. Hence, for any j, diag(apj)=As†Yj(ST)†. The features of the CPD are demonstrated as follows:(1)
The columns of the factor matrixAs are the eigenvectors of Yj1Yj2†. The corresponding eigenvalues are the first R larger eigenvalues, which are arranged in a descending order.(2)
The eigenvectors corresponding to theR-large eigenvalues of Yj1T(Yj2T)† constitute the factor matrix S.(3)
The rows of factor matrixAp can be recovered thanks to the relationship diag(apj)=As†Yj(ST)†.Obviously, whenJ>2, there are various combinations (j1,j2) that can be used for EVD. In order to make full use of the statistical properties of (24), we select CJ2 possible combinations for EVD to obtain Asj′, Sj′, where j′=1,…,CJ2. We follow the convention that the (j1,j2)-tuples are ordered lexicographically: the j1′th tuple (j1,j2) is preceding the j2′th tuple (j3,j4) if and only if j1⩽j3, j2⩽j4. According to the trivial indeterminacies aforementioned for CPD uniqueness, we know that Asj1′ and Asj2′, j1′≠j2′, coincide up to column permutation and scaling. In order to achieve the addition of multiple sets of results, we need to match the multiple sets of results. That is to say, the eigenvectors in each factor matrix obtained by the EVD need to be arranged in a descending order according to the corresponding eigenvalues, and all the eigenvectors need to be normalized. After the matching process, A~sj′ and S~j′ can be obtained, and then the following statistical operations are performed:(26)As=1CJ2∑j′=1CL2A~sj′,S=1CJ2∑j′=1CL2S~j′.Next, the factor matrixAp can be recovered according to the results of (26).


## 4.2. Parameters Estimation from the Factor Matrices

From the aforementioned tensor-based signal model, we can see that the spatial-polarization parameters of signals can be reflected by the factor matricesAp and the factor matrix As that only contains the spatial parameters of signals. Hence, once the factor matrices are obtained, we first need to extract the spatial parameters from As.


### 4.2.1. Estimation of Spatial Parameters

Spatial DOA(θr,ϕr) of the rth signal impinging on the array can be calculated directly by solving the following equations:(27)As1r=e-i2πb1Tϵr/λ⋮AsLr=e-i2πbLTϵr/λObviously, any two equations belonging to (27) can be used to derive the DOA of the rth signal. Therefore, we can get CL2 sets of solutions from (27). Improve the accuracy of the estimated parameters by averaging multiple sets of solutions:(28)θr,ϕr=1CL2∑m=1CL2θr,ϕrm,where (θr,ϕr) represents the solution of the mth subequations {[As]l1r,[As]l2r}. We follow the convention that the (l1,l2)-tuples are ordered lexicographically: the mth tuple (l1,l2) is preceding the nth tuple (l1′,l2′) if and only if l1⩽l1′, l2⩽l2′. Furthermore, we can construct the optimization problem as follows to estimate the DOAs of the signals:(29)θr,ϕr=arg⁡minθ,ϕ⁡ζsr∘asθ,ϕ-as:,r2,where as:,r represents the rth column of the factor matrix As obtained by CPD, ζsr represents the amplitude vector of as:,r, and the operator · denotes the 2-norm. The solution method to (29) and its performance can be used as a further work, and this paper discusses only the performance based on the solution of (28).


### 4.2.2. Estimation of Polarization Parameters

Once the estimation of DOA is obtained, the polarization parameters(γr,ηr) of the rth signal can be extracted from(30)ap:,r=BΞθr,ϕrcos⁡γr,sin⁡γreiηrT,where ap:,r represents the rth column of the factor matrix Ap obtained by CPD. Hereby, the estimated value of hγr,ηr can be expressed as h^γr,ηr=Ξθr,ϕr-1B-1ap:,r, and the estimated values of the polarization parameters can be obtained by(31)γr=arctan⁡h^γr,ηr2h^γr,ηr1,ηr=arg⁡h^γr,ηr2h^γr,ηr1,where arctan(·) is the arctangent operator and the operator arg(·) denotes phase acquisition.Obviously, the estimation accuracy of the DOA only depends on the factor matrixAs obtained by CPD, and the estimation accuracy of the polarization parameters is affected by both the estimation accuracy of the DOA and the matrix Ap obtained by CPD. It is well known that, in order to satisfy the validity of the parameters estimation, the algorithm must follow the premise of L>R. According to [14], there exists Kruskal bound ensuring the CPD uniqueness R<(J+L+K-2)/2. With the noise present, the accuracy of the CPD is severely impaired, when R is very close to the bound. However, as (J+L+K-2)/2 increases, the dependence of the CPD accuracy on (J+L+K-2)/2 will rapidly weaken. In addition, as can be seen from (26) and (28), the proposed algorithm can obtain multiple sets of parameters estimation results in each CPD and perform statistical operations on these results. However, what is different is that the subspace-based methods requires sufficient statistical information, which is usually obtained by increasing the number of sensors or snapshots, to ensure the accuracy of the parameters estimation, even if the number of the parameters to be estimated is very small. We, therefore, conclude that the proposed algorithm can achieve relatively desirable performance by using a smaller number of snapshots as well as a smaller number of sensors compared with the subspace-based methods. The subsequent simulation results also justify the above theoretical analysis.


## 4.2.1. Estimation of Spatial Parameters

Spatial DOA(θr,ϕr) of the rth signal impinging on the array can be calculated directly by solving the following equations:(27)As1r=e-i2πb1Tϵr/λ⋮AsLr=e-i2πbLTϵr/λObviously, any two equations belonging to (27) can be used to derive the DOA of the rth signal. Therefore, we can get CL2 sets of solutions from (27). Improve the accuracy of the estimated parameters by averaging multiple sets of solutions:(28)θr,ϕr=1CL2∑m=1CL2θr,ϕrm,where (θr,ϕr) represents the solution of the mth subequations {[As]l1r,[As]l2r}. We follow the convention that the (l1,l2)-tuples are ordered lexicographically: the mth tuple (l1,l2) is preceding the nth tuple (l1′,l2′) if and only if l1⩽l1′, l2⩽l2′. Furthermore, we can construct the optimization problem as follows to estimate the DOAs of the signals:(29)θr,ϕr=arg⁡minθ,ϕ⁡ζsr∘asθ,ϕ-as:,r2,where as:,r represents the rth column of the factor matrix As obtained by CPD, ζsr represents the amplitude vector of as:,r, and the operator · denotes the 2-norm. The solution method to (29) and its performance can be used as a further work, and this paper discusses only the performance based on the solution of (28).


## 4.2.2. Estimation of Polarization Parameters

Once the estimation of DOA is obtained, the polarization parameters(γr,ηr) of the rth signal can be extracted from(30)ap:,r=BΞθr,ϕrcos⁡γr,sin⁡γreiηrT,where ap:,r represents the rth column of the factor matrix Ap obtained by CPD. Hereby, the estimated value of hγr,ηr can be expressed as h^γr,ηr=Ξθr,ϕr-1B-1ap:,r, and the estimated values of the polarization parameters can be obtained by(31)γr=arctan⁡h^γr,ηr2h^γr,ηr1,ηr=arg⁡h^γr,ηr2h^γr,ηr1,where arctan(·) is the arctangent operator and the operator arg(·) denotes phase acquisition.Obviously, the estimation accuracy of the DOA only depends on the factor matrixAs obtained by CPD, and the estimation accuracy of the polarization parameters is affected by both the estimation accuracy of the DOA and the matrix Ap obtained by CPD. It is well known that, in order to satisfy the validity of the parameters estimation, the algorithm must follow the premise of L>R. According to [14], there exists Kruskal bound ensuring the CPD uniqueness R<(J+L+K-2)/2. With the noise present, the accuracy of the CPD is severely impaired, when R is very close to the bound. However, as (J+L+K-2)/2 increases, the dependence of the CPD accuracy on (J+L+K-2)/2 will rapidly weaken. In addition, as can be seen from (26) and (28), the proposed algorithm can obtain multiple sets of parameters estimation results in each CPD and perform statistical operations on these results. However, what is different is that the subspace-based methods requires sufficient statistical information, which is usually obtained by increasing the number of sensors or snapshots, to ensure the accuracy of the parameters estimation, even if the number of the parameters to be estimated is very small. We, therefore, conclude that the proposed algorithm can achieve relatively desirable performance by using a smaller number of snapshots as well as a smaller number of sensors compared with the subspace-based methods. The subsequent simulation results also justify the above theoretical analysis.


## 4.3. Computational Complexity

The parameters estimation algorithm using EVD-based CPD of tensors is summarized in Algorithm1, which is easy to implement by row-by-row access. Note that the computational complexity scales linearly to the number of loops, since the procedure of the algorithm runs serially.Algorithm 1:EVD-based CPD for parameters estimation.
Require:  Y∈CJ×L×K, L, J, K, R. 1: Reshape the matrciesY1,…,YJ∈CL×K by unfolding the tensor Y; 2: Forj′=1 to CJ2 do 3:Calculate Yj1Yj2†, Yj1T(Yj2T)†; 4:Obtain Asj′ by executing the EVD of Yj1Yj2†; 5:Obtain Sj′ by executing the EVD of Yj1T(Yj2T)†; 6:Obtain A~sj′ and S~j′ by matching processing. 7: End for 8:As=1CJ2∑j′=1CL2A~sj′,S=1CJ2∑j′=1CL2S~j′; 9: RecoverAp according to diag(apj)=As†Yj(ST)†; 10: Forr=1 to R do 11:For m=1 to CL2 do 12:Calculate (θr,ϕr)m by solving equations {[As]l1r,[As]l2r}; 13:End for 14:θr,ϕr=1CL2∑m=1CL2(θr,ϕr)m; 15:h^γr,ηr=Ξθr,ϕr-1B-1ap:,r; 16:γr=arctanh^γr,ηr2h^γr,ηr1,ηr=argh^γr,ηr2h^γr,ηr1. 17: End for Ensure:  Ψr=(θr,ϕr,γr,ηr), r=1,…,R.In fact, the EVDs of (25) follow from the generalized EVDs (GEVDs) of the matrix pencils (Yj1,Yj2) and (Yj1T,Yj2T). In order to evaluate the complexity of the proposed algorithm, we first need to know the number of multiplications required by GEVD. An effective GEVD solution is given by the method of QZ iteration [25, pp.413-414], in which pLKR multiplications are required for two L×K matrices with rank-R, where the constant p depends on the selected iterative strategy, e.g., for Hessenberg-Triangular Reduction, p=8. The Moore-Penrose inverse of the M×N matrix or N×M matrix involves 4MN2+8N3 multiplications, where M>N. Additional J×(RLK+R2K) multiplications are also required to recover the factor matrix Ap. Therefore, the computational complexity of the EVD-based CPD does not exceed (pCJ2+8+J)K·[O(L2)], since the premise of L,J,K>R has been established. Furthermore, in order to compare the computational complexity for the proposed algorithm, we analyze the required number of multiplications for subspace-based methods. For the aforementioned matrix-based method, the computational complexity of obtaining subspaces is approximately (L+K)·[O(L2)]. Considering the operation process of HOSVD based on orthogonal iterations [25, pp.454-455], the computational complexity of subspace for tensor model is about (KJ2+8J2+pJ)·[O(L2)]. Obviously, the conclusion that the computational complexity of the CPD algorithm is lower than that of the tensor-based HOSVD methods can be obtained through the above analysis. Additionally, for the extraction process of the parameters, the computational complexity of the search strategies is generally much higher than that of the equation-solving methods, under the premise of ensuring the estimation accuracy. Therefore, the computational complexity of the proposed algorithm is lower than that of the search-based methods.Additionally, according to the procedure of the ALS-based CPD described in AppendixC, the computational complexity of the factor matrix S is not less than (5+L+K)·[O(R3)], since the premise of L,J,K>R has been established. Similarly, the computational complexity of the factor matrices Ap and As is not less than (5+2K)·[O(R3)]. Therefore, the computational complexity of the ALS-based CPD is about (15+L+5K)G·[O(R3)], where G represents the number of iterations. We use the Matlab 2016a to implement the after-mentioned simulations on a computer with Intel Core i5 CPU6200U 2.3GHz and 8GB memory running Windows 10. From the perspective of computing time, the time required for ALS-based CPD is about twice that of EVD-based CPD to obtain the subsequent simulation results.


## 5. Simulation Results

The superiorities of the proposed algorithm will be demonstrated by the results of numerical simulations in this section. We consider the following scenario with a uniform linear array (ULA) separated by half wavelength. The array consists of complete EM vector-sensors, and two signals with the polarization parameters(γ1,η1)=(10∘,90∘) and (γ2,η2)=(30∘,90∘) impinge on the array from (θ1,ϕ1)=(50∘,90∘) and (θ2,ϕ2)=(120∘,90∘), respectively. Both signals satisfy the aforementioned assumptionsA1 toA3. The modulation modes of the signals are BPSK and QBPSK, respectively, and the symbol transmission rate of both signals is 320B. Additionally, modulation frequency and sampling frequency are 16kHz and 32kHz, respectively. The noise component is assumed to be zero-mean additive white Gaussian noise. Performance of the EVD-based CPD algorithm is compared with three other subspace-based algorithms, MUSIC, HOSVD-MUSIC, HOSVD-ESPRIT, and the ALS-based CPD algorithm. The performances of the algorithms are evaluated by the Root Mean Square Error (RMSE):(32)RMSEμ=1NR∑n=1N∑r=1Rμ^rn-μr2,where μr denotes one of the parameters {θr,ϕr,γr,ηr}, and μ^rn is the estimation of μr in the nth trial. Throughout all simulations, results are averaged by 500 Monte-Carlo trails and compared with the Cramér-Rao lower bound (CRB) benchmark, which is described in Appendix D.Example 3.
Assume that the number of complete EM vector-sensors is 4, and the signal-to-noise ratio (SNR) is 30dB. The proposed EVD-based CPD algorithm is compared with the subspace-based algorithms and the ALS-based CPD algorithm in terms of the DOA RMSE versus the different number of snapshotsK. Figure 1 illustrates a superior performance of the proposed EVD-based CPD algorithm with respect to the other algorithms when the number of snapshots is larger than 50. It is worth noting that, when the number of snapshots is greater than 100, the performance change of the ALS-based CPD algorithm with the increase of snapshots is negligible. In the case where the number of snapshots is extremely small, the main reason for the poor performance of the EVD-based CPD algorithm is that the number of the signals to be estimated is close to the uniqueness condition for CPD.Figure 1
DOA RMSE versusK for L=4 and SNR=30dB.Example 4.
Scenarios with various numbers of sensors (L∈{4,6,8}) are considered in this example. Each of these scenarios is tested versus SNR, and the number of snapshots is fixed to K=200. The following conclusions can be obtained by observing Figure 2.(1)
For all algorithms in each case, RMSEs decrease gradually with the increase of SNR. Note that the EVD-based CPD algorithm performs better than the other algorithms from case to case when SNR is greater than10dB.(2)
In the case of low SNR, i.e., SNR is less than10dB, the problem of the poor estimation accuracy for the proposed EVD-based CPD algorithm is effectively improved with the increase of the number of sensors.(3)
The accuracy of the DOA estimation for MUSIC algorithm changes better obviously, as the number of sensors increases. This is due to the fact that the MUSIC algorithm determines the DOA by executing an exhaustive search [21].Figure 2
DOA RMSE versus SNR forK=200 and various numbers of sensors. (a) Array of L=4, (b) array of L=6, and (c) array of L=8.
(a) (b) (c)Example 5.
Once the estimation of DOA is obtained in Example4, the polarization parameters can be estimated. Figure 3 shows the estimation accuracy of the parameter γ by the EVD-based CPD and the other algorithms in the scenarios in Example 4. The differences in performance between algorithms are similar to the tendencies emerged in Figure 2. Similarly, the performance of the ALS-based CPD is significantly improved as the number of sensors increases. However, the RMSEs of the algorithms for γ are far away from the CRB boundary, especially for L=4 and L=6. This loss can be explained by the 2-step method of the estimation process, which consists of estimating DOA parameters firstly, and then corresponding polarization parameters in a second step. In other words, the accumulation of errors in the 2-step calculation leads to this loss.Figure 3
Polarization parameter’s RMSE versus SNR forK=200 and various numbers of sensors. (a) Array of L=4, (b) array of L=6, and (c) array of L=8.
(a) (b) (c)


## 6. Conclusion

In this paper, we proposed an EVD-based CPD approach for tensors to tackle the problem of signal parameters estimation. The procedure of the EVD-based CPD was described, followed by expatiations and simulations of the EVD-based CPD algorithm applied to the signal parameters estimation. The traditional subspace-based algorithms, MUSIC, HOSVD-MUSIC HOSVD-ESPRIT, and the classical ALS-based CPD algorithm, were introduced to verify the efficacy of the proposed EVD-based CPD algorithm.The strength of the proposed algorithm lies in the estimation of the factor matrices from the tensor model of the signal, which is then used to extract the signal parameters. The signal parameters can be estimated by virtue of the diversities of the spatial and polarization belonging to the factor matrixes. Furthermore, when the number of sensors is small in the array, the proposed algorithm achieves a lower RMSE.



---

*Source: 7389306-2019-02-03.xml*