# Estimation of DOA for Noncircular Signals via Vandermonde Constrained Parallel Factor Analysis

Authors: Heyun Lin; Chaowei Yuan; Jianhe Du; Zhongwei Hu
Journal: International Journal of Antennas and Propagation
 (2018)

---

## Abstract

We provide a complete study on the direction-of-arrival (DOA) estimation of noncircular (NC) signals for uniform linear array (ULA) via Vandermonde constrained parallel factor (PARAFAC) analysis. By exploiting the noncircular property of the signals, we first construct an extended matrix which contains two times sampling number of the received signal. Then, taking the Vandermonde structure of the array manifold matrix into account, the extended matrix can be turned into a tensor model which admits the Vandermonde constrained PARAFAC decomposition. Based on this tensor model, an efficient linear algebra algorithm is applied to obtain the DOA estimation via utilizing the rotational invariance between two submatrices. Compared with some existing algorithms, the proposed method has a better DOA estimation performance. Meanwhile, the proposed method consistently has a higher estimation accuracy and a much lower computational complexity than the trilinear alternating least square- (TALS-) based PARAFAC algorithm. Finally, numerical examples are conducted to demonstrate the effectiveness of the proposed approach in terms of estimation accuracy and computational complexity.


---

## Body

## 1. Introduction

Direction-of-arrival (DOA) estimation of signals impinging on an array of sensors is an important and fundamental issue in array signal processing due to its wide applications in wireless communications, geophysics, radar, and so on [1–3]. In this context, many DOA estimators have been developed to solve this problem, such as propagator method (PM) [4], maximum likelihood (ML) methods [5, 6], tensor-based method [7], estimation of signal parameters via rotational invariance technique (ESPRIT) algorithm [8], and multiple signal classification (MUSIC) algorithm [9].Unfortunately, all of the above works [4–9] ignore the characteristics of the incident signals. It has been shown that the accuracy of DOA estimation can be enhanced by taking advantage of the noncircular (NC) property of the impinging signals [10]. Examples of the noncircular signals include binary phase-shift keying (BPSK) and amplitude-modulated (AM) signals which are widely applied in wireless telecommunication systems. While the property of NC signals is utilized, the array aperture is virtually doubled in [10], which yields a better DOA estimation performance than that in [9]. Besides taking the noncircularity of signals into account, the maximum number of sources can potentially exceed the number of array elements. Consequently, solid researches on NC signals for DOA estimation have appeared in many literatures [11–14]. NC ESPRIT algorithm and NC Unitary ESPRIT algorithm for DOA estimation were presented in [11, 12], respectively. The work [13] developed a low-complexity noncircular rational invariance propagator method (NC-RI-PM) for angle estimation, which is used for uniform linear array (ULA). However, the performance of the NC-RI-PM became worse rapidly in the case of low signal-to-noise ratio (SNR). The authors in [14] utilized the parallel factor (PARAFAC) analysis [15] to acquire the two-dimensional (2D) angle estimation of NC signals for uniform rectangular array (URA). The main drawback in [14] is that it has large computational load and is sensitive to the iterative initial value.In this paper, we present a Vandermonde constrained PARAFAC method to improve the DOA estimation accuracy by taking advantages of the property of NC signals and the array manifold matrix with Vandermonde structure. In order to construct a PARAFAC data model, we first build an extended matrix by concatenating the received data and its conjugated component. Then, the extended data matrix can be formulated as a PARAFAC model with Vandermonde constrained factor matrix. Finally, a closed-form solution can be utilized to acquire the angle information by taking the Vandermonde structure into account.The contributions of this paper can be summarized as follows: (1) We combine the property of NC signals and the array manifold matrix with Vandermonde structure, based on which a Vandermonde constrained PARAFAC method is proposed to acquire the DOA of NC signals for uniform linear array. (2) Compared with trilinear alternating least square- (TALS-) based PARAFAC method, the proposed method provides a better angle estimation as well as a lower computational complexity. (3) The proposed approach also has a higher estimation accuracy than the ESPRIT method [8], NC ESPRIT method [11], and NC-RI-PM method [13].The rest of this paper is structured as follows. In Section2, we introduce the system model for ULA. Section 3 presents the Vandermonde constrained PARAFAC method for DOA estimation of noncircular signals. The complexity analysis and Cramér-Rao Bound (CRB) of NC signal DOA estimation are given in Section 4. Numerical results are provided in Section 5. Finally, Section 6 concludes the paper.Notations: scalars, vectors, matrices, and tensors are denoted by lowercase letters, lowercase boldface letters, uppercase boldface letters, and calligraphic letters, respectively. The operators⋅∗,⋅T,⋅H,⋅†, and ⋅F denote the conjugate, transpose, conjugate transpose, pseudoinverse, and Frobenius norm, respectively. IK stands for a K×K identity matrix. diaga is the diagonal matrix with diagonal elements given by a. ⊙, ∘, and ⊕ represent the Khatri-Rao product, outer product, and Hadamard product, respectively. rA is the rank of a matrix A.


## 2. System Model

Let us consider a uniform linear array consisting ofM omnidirectional sensors. The distance between two adjacent sensors is d. As shown in Figure 1, the first sensor is regarded as the reference point for the array, assuming that there are K uncorrelated narrowband far-field signals impinging on this array from distinct directions θ1,…,θK. We also assume that the number of sources K is known a priori or has been estimated by the methods shown in [16, 17]. Under such a scenario, the received signal vector xn∈ℂM×1 can be modeled as
(1)xn=Asn+vn,where
(2)A=1⋯1α1⋯αK⋮⋱⋮α1M−1⋯αKM−1is the array manifold matrix with size of M×K, αk=e−j2π/λdsinθk,k=1,…,K.λ is the wavelength; sn=s1n,…,sKnT∈ℂK×1 is the signal vector at snapshot n. vn∈ℂM×1 stands for the additive Gaussian noise vector with zero-mean and covariance σ2IM. From (2), we know that A is a Vandermonde matrix with distinct nonzero generators α1,…,αK. In this paper, the strictly second-order noncircular signals [18] are considered. Thus, we have Esn=0K×1,EsnsHn≠0K×K, and EsnsTn≠0K×K.Figure 1
Illustration of the array geometry.It is assumed that the DOAs of the sources are constant inN snapshots. We collect N snapshots corresponding to xn; the received signal matrix X∈ℂM×N is given by
(3)X=x1,…,xN=AST+V,where S=s1…sNT∈ℂN×K stands for the source matrix and V∈ℂM×N contains the samples of the additive noise. Due to the NC property, we can decompose the source matrix as S=S0Ψ [18], where S0∈ℝN×K is a real-valued matrix, Ψ=diage−jψ1,…,e−jψK∈ℂK×K is a diagonal matrix, and ψk is referred to as the phase shift of the kth source. Then (3) can be rewritten as
(4)X=AΨS0T+V.


## 3. Estimation of DOA via Vandermonde Constrained PARAFAC Decomposition

### 3.1. PARAFAC Model

Since the noncircular signals are considered, we construct the extended data matrixX¯∈ℂ2M×N as follows [12, 14]:
(5)X¯=XJX∗=AΨS0TJAΨS0T∗+VJV∗=AΨAΛΨ−1S0T+V1=AD1ΦAD2ΦS0T+V1=Φ⊙AS0T+V1,where J∈ℝM×M is the exchange matrix with ones on its antidiagonal and zeros elsewhere. V1 is the extended noise matrix. Λ=diagej2π/λM−1dsinθ1,…,ej2π/λM−1dsinθK∈ℂK×K, and
(6)Φ=e−jψ1,…,e−jψKej2πλM−1dsinθ1+ψ1,…,ej2πλM−1dsinθK+ψKis a 2×K matrix. DiΦ is the diagonal matrix with diagonal elements given by the ith row of Φ. According to Harshman [15], the signal part of (5) can be seen as a matrix representation of a third-order PARAFAC model X∈ℂM×N×P, whose m,n,p−th element is given by
(7)x¯m,n,p=∑k=1Kam,ks0n,kϕp,k,where m=1,…,M,n=1,…,N, and p=1,2.x¯m,n,p represents the typical element of the tensor X. am,k denotes the m,k−th element of A, and similarly for the others. Using the reshaping operations, we can obtain the other two equivalent matrix representations of X as follows:
(8)Y¯=A⊙S0ΦT+V2∈ℂMN×2,Z¯=S0⊙ΦAT+V3∈ℂ2N×M.Till now, we successfully link the expanded data matrix with a tensor model which satisfies the PARAFAC decomposition. According to (5), the PARAFAC decomposition can be accomplished by solving the following unconstrained optimization problem.
(9)minA,S0,Φχ−∑k=1Kak∘s0k∘ϕkF2,where ak, s0k, and ϕk represent the kth columns of the matrices A, S0, and Φ, respectively. This is most often done by means of the TALS method. As its programming simplicity, the TALS algorithm has been successfully employed to estimate the parameters of several tensor models [7, 19, 20]. The basic idea of the TALS algorithm for fitting the proposed PARAFAC model is stated as follows:With the consideration of (5), while A and Φ are fixed, we update S0 using the following least squares (LS) fitting:
(10)S^0=argminS0X¯−Φ⊙AS0TF=Φ⊙A†X¯T,where S^0 stands for the estimation of S0. In the same way, while A and S0 are fixed, Φ can be updated by the following LS fitting:
(11)Φ^=argminΦY¯−A⊙S0ΦTF=A⊙S0†Y¯T,where Φ^ stands for the estimation of Φ. Similarly, we obtain an updated A via LS fitting as
(12)A^=argminAZ¯−S0⊙ΦATF=S0⊙Φ†Z¯T,where A^ stands for the estimation of A. Repeat updating process of the three matrices until convergence, and then the estimation matrices A^,S^0, and Φ^ can be obtained. It is well known that the TALS algorithm is sensitive to the initial value as well as the array manifold, which may cause the TALS algorithm converging slowly. In the next subsection, a more efficient algorithm is applied to obtain the estimation of DOA for noncircular signals.


### 3.2. DOA Estimation with Vandermonde Constrained PARAFAC Method

From the previous analysis, we know that the matrixA is a Vandermonde matrix with distinct nonzero generators α1,…,αK. Therefore, by taking the Vandermonde factor into account, the model X can be seen as a tensor which admits the PARAFAC decomposition with Vandermonde constraint. And then, a more efficient algorithm is applied to fitting the PARAFAC model with Vandermonde constraint.According to (6), we know that rΦ=min2,K. If K>2, we can conclude that the matrix Φ does not have full-column rank, which implies that rΦ is less than K. Then, the rank of Y¯ is also less than K. Consequently, the spatial smoothing technique can be used to handle with this problem [21]. The detailed process of spatial smoothing is presented as follows:DenoteY¯l2 as the l2th subarray of Y¯ such that Y¯l2=Y¯Nl2−1+1:NL1+l2−1,:, l2=1,…,L2, where L2 is the smoothing factor. The integers L1 and L2 must be chosen such that M=L1+L2−1. Taking the advantage of the Vandermonde structure of A, we have
(13)Y¯l2=S0Dl2AS0Dl2+1A⋮S0Dl2+L1−1AΦT=S0D1AS0D2A⋮S0DL1AΩl2−1ΦT,where Ω=diagej2π/λdsinθ1,…,ej2π/λdsinθK∈ℂK×K. Then, denote Y˜=Y¯1…Y¯L2∈ℂNL1×2L2, we can obtain
(14)Y˜=Y¯1⋯Y¯L2=S0D1AS0D2A⋮S0DL1AΦTΩΦT⋯ΩL2−1ΦT=AL1⊙S0ΦTΩΦT⋯ΩL2−1ΦT,where AL1∈ℂL1×K denotes the first L1 rows of A. Note that
(15)ΦTΩΦT⋯ΩL2−1ΦTT=ΦD1AΦD2A⋮ΦDL2A=AL2⊙Φ,where AL2∈ℂL2×K denotes the first L2 rows of A.Therefore, (14) can be rewritten as
(16)Y˜=AL1⊙S0AL2⊙ΦT.In this way, the matrix representationY¯ of the tensor X, which involves a rank-deficient matrix Φ (while K>2), is replaced by a matrix representation Y˜ which only involves two full-column rank matrices AL1⊙S0 and AL2⊙Φ. The matrix Y˜ can be seen as a matrix representation of a constrained PARAFAC decomposition model [22–24] with three factors AL1,S0 and AL2⊙Φ. When the matrix Y˜ is obtained, a more efficient algorithm can be used to fit the Vandermonde constrained PARAFAC model [23, 24].While the matricesAL1⊙S0 and AL2⊙Φ have full-column rank, we perform the compact singular value decomposition (SVD) of Y˜, which is given by Y˜≈UΣVH, where Σ∈ℂK×K is a diagonal matrix consisting of the largest K nonzero singular values. U∈ℂL1N×K and V∈ℂ2L2×K are matrices composed of left and right singular vectors associated with the largest K nonzero singular values, respectively. Since the factor matrices AL1⊙S0 and AL2⊙Φ have full-column rank, we can conclude that the subspace spanned by the column vectors of U is the same as the subspace spanned by the column vectors of the matrix AL1⊙S0, so we have
(17)spanAL1⊙S0=spanU.Therefore, there exists a unique nonsingular matrixM∈ℂK×K such that
(18)UM=AL1⊙S0.Using the inherent Vandermonde structure ofAL1, we define the following two submatrices:
(19)U1M=A¯L1⊙S0=S0D1AS0D2A⋮S0DL1−1A,(20)U2M=A¯L1⊙S0=S0D2AS0D3A⋮S0DL1A,where U1=U1:L1−1N,1:K and U2=UN+1:L1N,1:K, and the relationship between AL1, A¯L1, and A¯L1 is that
(21)AL1=A¯L1lastrow=firstrowA¯L1.Therefore, there exists a diagonal matrixZ∈ℂK×K such that
(22)A¯L1⊙S0Z=A¯L1⊙S0,where Z=diagz1,…,zK with generators zk=e−j2π/λdsinθk, k=1,…,K.According to (19), (20), (21) and (22), we can obtain that
(23)U2M=U1MZ.From (23), we have U2=U1MZM−1=U1Z^, where Z^=MZM−1. Assuming that the matrix A¯L1⊙S0 has full-column rank, we can conclude that the matrices U1 and U2 have full-column rank and Z^ can be rewritten as
(24)Z^=U1†U2=MZM−1.Therefore, the generatorszk, k=1,…,K, can be obtained from the eigenvalue decomposition (EVD) of U1†U2. Then, the DOA of the kth source θk can be recovered via
(25)θ^k=arcsin−λanglezk2πdk=1,…,K,where θ^k stands for the estimation of θk and angle(·) is used for obtaining the phase. The proposed Vandermonde constrained PARAFAC algorithm for NC signals DOA estimation is presented in Algorithm 1.Algorithm 1: The Vandermonde constrained PARAFAC algorithm for NC signals DOA estimation.
Step 1. Construct the tensor model X with the matrix representation Y¯=A⊙S0ΦT by exploiting the property of NC signals.
Step 2. Choose the pair (L1,L2) subject to L1+L2=M+1 and construct the matrix representation Y˜=AL1⊙S0AL2⊙ΦT via spatial smoothing.
Step 3. Perform the compact SVD of Y˜:Y˜≈UΣVH.
Step 4. Construct U1=U1:L1−1N,1:K,U2=UN+1:L1N,1:K.
Step 5. Perform the EVD of U1†U2:U1†U2=MZM−1, Z=diagz1,…,zK with generators zk=e−j2π/λdsinθk, k=1,…,K.
Step 6. Compute θ^k=arcsin−λanglezk/2πd, k=1,…,K.Remark 1.
The proposed method in this paper is developed for noncircular signals. Therefore, we assume that the circular signals are not present in this paper. When the mixed noncircular and circular signals are considered, the DOAs can be estimated by the algorithms reported in [25, 26].


### 3.3. Identifiabilitity Analysis

The most remarkable characteristic of the PARAFAC decomposition is its essential uniqueness property, which makes PARAFAC decomposition a key tool for signal separation. In this subsection, we discuss the identifiabilitity condition for the proposed PARAFAC method.Theorem 1 [7, 15, 27, 28].
LetX∈ℂM×N×2 be a tensor with matrix representation Y¯=A⊙S0ΦT, where A∈ℂM×K,S0∈ℝN×K, and Φ∈ℂ2×K, if
(26)kA+kS0+kΦ≥2K+2,then the three factor matrices A, S0, and Φ are unique up to permutation and scaling of columns, meaning that any other triple (A^,S^0,Φ^) is related to (A,S0,Φ) via
(27)A^=AΠΔ1,S^0=S0ΠΔ2,Φ^=ΦΠΔ3,where kA,kS0, and kΦ are the k−ranks [27, 28] of the matrices A, S0, and Φ, respectively. Π is a permutation matrix. Δi, i=1,…,3, are diagonal scaling matrices satisfying Δ1Δ2Δ3=IK.If the three factor matrices have full rank, then condition (26) becomes
(28)minM,K+minN,K+min2,K≥2K+2.In this paper, we assumeN≥K,K≥2, then the identifiable condition (28) becomes M≥K. Therefore, the maximum number of sources that can be handled by the PARAFAC algorithm is no more than M.WhenX∈ℂM×N×P is a tensor that admits the Vandermonde constrained PARAFAC decomposition, we have the following uniqueness result.Theorem 2 [23, 24].
LetX∈ℂM×N×2 be a tensor with matrix representation Y¯=A⊙S0ΦT, where A∈ℂM×K is a Vandermonde matrix with distinct generators, S0∈ℝN×K and Φ∈ℂ2×K. Consider the matrix representation Y˜=AL1⊙S0AL2⊙ΦT with L1+L2=M+1. If
(29)rA¯L1⊙S0=K,rAL2⊙Φ=K,for some L1+L2=M+1, then K=rX, and the Vandermonde constrained PARAFAC decomposition of X is unique. Generically, condition (29) is satisfied if and only if
(30)minL1−1N,2L2≥K.When we assumeN≥M and L1≥2, the condition (30) becomes 2L2≥K. Therefore, the maximum number of sources that can be identified by the Vandermonde constrained PARAFAC algorithm is 2L2. If 2L2>M, the Vandermonde constrained PARAFAC method can identify more sources than the PARAFAC method under the same conditions. Notably, the scaling does not involve the Vandermonde factor matrix in a PARAFAC model with Vandermode constraint. This means that the estimated matrices A˜,S˜0, and Φ˜ are related to A, S0, and Φ via
(31)A˜=AΠ,S˜0=S0ΠΔ˜2,Φ˜=ΦΠΔ˜3,where Δ˜j,j=2,3, are diagonal scaling matrices satisfying Δ˜2Δ˜3=IK.Remark 2.
The proposed method can be seen as a generalized ESPRIT method. In line with the experiments for ESPRIT in [29, 30], we can choose the pair (L1,L2) such that the dimensions of the matrices AL1⊙S0 and AL2⊙Φ are close, with the inequality minL1−1N,2L2≥K satisfied. According to [29, 30], we know that such a direct method will yield a very good estimation result. Therefore, in order to reduce the computational complexity, L1 is fixed to 3 and L2=M−2 in this paper. Although it may not yield the best estimation result with the parameters being set as L1=3 and L2=M−2, a good estimation result can be still achieved. Note that simulations in Section 5 show that the DOA performance of such a direct approach still outperforms than other methods.


## 3.1. PARAFAC Model

Since the noncircular signals are considered, we construct the extended data matrixX¯∈ℂ2M×N as follows [12, 14]:
(5)X¯=XJX∗=AΨS0TJAΨS0T∗+VJV∗=AΨAΛΨ−1S0T+V1=AD1ΦAD2ΦS0T+V1=Φ⊙AS0T+V1,where J∈ℝM×M is the exchange matrix with ones on its antidiagonal and zeros elsewhere. V1 is the extended noise matrix. Λ=diagej2π/λM−1dsinθ1,…,ej2π/λM−1dsinθK∈ℂK×K, and
(6)Φ=e−jψ1,…,e−jψKej2πλM−1dsinθ1+ψ1,…,ej2πλM−1dsinθK+ψKis a 2×K matrix. DiΦ is the diagonal matrix with diagonal elements given by the ith row of Φ. According to Harshman [15], the signal part of (5) can be seen as a matrix representation of a third-order PARAFAC model X∈ℂM×N×P, whose m,n,p−th element is given by
(7)x¯m,n,p=∑k=1Kam,ks0n,kϕp,k,where m=1,…,M,n=1,…,N, and p=1,2.x¯m,n,p represents the typical element of the tensor X. am,k denotes the m,k−th element of A, and similarly for the others. Using the reshaping operations, we can obtain the other two equivalent matrix representations of X as follows:
(8)Y¯=A⊙S0ΦT+V2∈ℂMN×2,Z¯=S0⊙ΦAT+V3∈ℂ2N×M.Till now, we successfully link the expanded data matrix with a tensor model which satisfies the PARAFAC decomposition. According to (5), the PARAFAC decomposition can be accomplished by solving the following unconstrained optimization problem.
(9)minA,S0,Φχ−∑k=1Kak∘s0k∘ϕkF2,where ak, s0k, and ϕk represent the kth columns of the matrices A, S0, and Φ, respectively. This is most often done by means of the TALS method. As its programming simplicity, the TALS algorithm has been successfully employed to estimate the parameters of several tensor models [7, 19, 20]. The basic idea of the TALS algorithm for fitting the proposed PARAFAC model is stated as follows:With the consideration of (5), while A and Φ are fixed, we update S0 using the following least squares (LS) fitting:
(10)S^0=argminS0X¯−Φ⊙AS0TF=Φ⊙A†X¯T,where S^0 stands for the estimation of S0. In the same way, while A and S0 are fixed, Φ can be updated by the following LS fitting:
(11)Φ^=argminΦY¯−A⊙S0ΦTF=A⊙S0†Y¯T,where Φ^ stands for the estimation of Φ. Similarly, we obtain an updated A via LS fitting as
(12)A^=argminAZ¯−S0⊙ΦATF=S0⊙Φ†Z¯T,where A^ stands for the estimation of A. Repeat updating process of the three matrices until convergence, and then the estimation matrices A^,S^0, and Φ^ can be obtained. It is well known that the TALS algorithm is sensitive to the initial value as well as the array manifold, which may cause the TALS algorithm converging slowly. In the next subsection, a more efficient algorithm is applied to obtain the estimation of DOA for noncircular signals.


## 3.2. DOA Estimation with Vandermonde Constrained PARAFAC Method

From the previous analysis, we know that the matrixA is a Vandermonde matrix with distinct nonzero generators α1,…,αK. Therefore, by taking the Vandermonde factor into account, the model X can be seen as a tensor which admits the PARAFAC decomposition with Vandermonde constraint. And then, a more efficient algorithm is applied to fitting the PARAFAC model with Vandermonde constraint.According to (6), we know that rΦ=min2,K. If K>2, we can conclude that the matrix Φ does not have full-column rank, which implies that rΦ is less than K. Then, the rank of Y¯ is also less than K. Consequently, the spatial smoothing technique can be used to handle with this problem [21]. The detailed process of spatial smoothing is presented as follows:DenoteY¯l2 as the l2th subarray of Y¯ such that Y¯l2=Y¯Nl2−1+1:NL1+l2−1,:, l2=1,…,L2, where L2 is the smoothing factor. The integers L1 and L2 must be chosen such that M=L1+L2−1. Taking the advantage of the Vandermonde structure of A, we have
(13)Y¯l2=S0Dl2AS0Dl2+1A⋮S0Dl2+L1−1AΦT=S0D1AS0D2A⋮S0DL1AΩl2−1ΦT,where Ω=diagej2π/λdsinθ1,…,ej2π/λdsinθK∈ℂK×K. Then, denote Y˜=Y¯1…Y¯L2∈ℂNL1×2L2, we can obtain
(14)Y˜=Y¯1⋯Y¯L2=S0D1AS0D2A⋮S0DL1AΦTΩΦT⋯ΩL2−1ΦT=AL1⊙S0ΦTΩΦT⋯ΩL2−1ΦT,where AL1∈ℂL1×K denotes the first L1 rows of A. Note that
(15)ΦTΩΦT⋯ΩL2−1ΦTT=ΦD1AΦD2A⋮ΦDL2A=AL2⊙Φ,where AL2∈ℂL2×K denotes the first L2 rows of A.Therefore, (14) can be rewritten as
(16)Y˜=AL1⊙S0AL2⊙ΦT.In this way, the matrix representationY¯ of the tensor X, which involves a rank-deficient matrix Φ (while K>2), is replaced by a matrix representation Y˜ which only involves two full-column rank matrices AL1⊙S0 and AL2⊙Φ. The matrix Y˜ can be seen as a matrix representation of a constrained PARAFAC decomposition model [22–24] with three factors AL1,S0 and AL2⊙Φ. When the matrix Y˜ is obtained, a more efficient algorithm can be used to fit the Vandermonde constrained PARAFAC model [23, 24].While the matricesAL1⊙S0 and AL2⊙Φ have full-column rank, we perform the compact singular value decomposition (SVD) of Y˜, which is given by Y˜≈UΣVH, where Σ∈ℂK×K is a diagonal matrix consisting of the largest K nonzero singular values. U∈ℂL1N×K and V∈ℂ2L2×K are matrices composed of left and right singular vectors associated with the largest K nonzero singular values, respectively. Since the factor matrices AL1⊙S0 and AL2⊙Φ have full-column rank, we can conclude that the subspace spanned by the column vectors of U is the same as the subspace spanned by the column vectors of the matrix AL1⊙S0, so we have
(17)spanAL1⊙S0=spanU.Therefore, there exists a unique nonsingular matrixM∈ℂK×K such that
(18)UM=AL1⊙S0.Using the inherent Vandermonde structure ofAL1, we define the following two submatrices:
(19)U1M=A¯L1⊙S0=S0D1AS0D2A⋮S0DL1−1A,(20)U2M=A¯L1⊙S0=S0D2AS0D3A⋮S0DL1A,where U1=U1:L1−1N,1:K and U2=UN+1:L1N,1:K, and the relationship between AL1, A¯L1, and A¯L1 is that
(21)AL1=A¯L1lastrow=firstrowA¯L1.Therefore, there exists a diagonal matrixZ∈ℂK×K such that
(22)A¯L1⊙S0Z=A¯L1⊙S0,where Z=diagz1,…,zK with generators zk=e−j2π/λdsinθk, k=1,…,K.According to (19), (20), (21) and (22), we can obtain that
(23)U2M=U1MZ.From (23), we have U2=U1MZM−1=U1Z^, where Z^=MZM−1. Assuming that the matrix A¯L1⊙S0 has full-column rank, we can conclude that the matrices U1 and U2 have full-column rank and Z^ can be rewritten as
(24)Z^=U1†U2=MZM−1.Therefore, the generatorszk, k=1,…,K, can be obtained from the eigenvalue decomposition (EVD) of U1†U2. Then, the DOA of the kth source θk can be recovered via
(25)θ^k=arcsin−λanglezk2πdk=1,…,K,where θ^k stands for the estimation of θk and angle(·) is used for obtaining the phase. The proposed Vandermonde constrained PARAFAC algorithm for NC signals DOA estimation is presented in Algorithm 1.Algorithm 1: The Vandermonde constrained PARAFAC algorithm for NC signals DOA estimation.
Step 1. Construct the tensor model X with the matrix representation Y¯=A⊙S0ΦT by exploiting the property of NC signals.
Step 2. Choose the pair (L1,L2) subject to L1+L2=M+1 and construct the matrix representation Y˜=AL1⊙S0AL2⊙ΦT via spatial smoothing.
Step 3. Perform the compact SVD of Y˜:Y˜≈UΣVH.
Step 4. Construct U1=U1:L1−1N,1:K,U2=UN+1:L1N,1:K.
Step 5. Perform the EVD of U1†U2:U1†U2=MZM−1, Z=diagz1,…,zK with generators zk=e−j2π/λdsinθk, k=1,…,K.
Step 6. Compute θ^k=arcsin−λanglezk/2πd, k=1,…,K.Remark 1.
The proposed method in this paper is developed for noncircular signals. Therefore, we assume that the circular signals are not present in this paper. When the mixed noncircular and circular signals are considered, the DOAs can be estimated by the algorithms reported in [25, 26].


## 3.3. Identifiabilitity Analysis

The most remarkable characteristic of the PARAFAC decomposition is its essential uniqueness property, which makes PARAFAC decomposition a key tool for signal separation. In this subsection, we discuss the identifiabilitity condition for the proposed PARAFAC method.Theorem 1 [7, 15, 27, 28].
LetX∈ℂM×N×2 be a tensor with matrix representation Y¯=A⊙S0ΦT, where A∈ℂM×K,S0∈ℝN×K, and Φ∈ℂ2×K, if
(26)kA+kS0+kΦ≥2K+2,then the three factor matrices A, S0, and Φ are unique up to permutation and scaling of columns, meaning that any other triple (A^,S^0,Φ^) is related to (A,S0,Φ) via
(27)A^=AΠΔ1,S^0=S0ΠΔ2,Φ^=ΦΠΔ3,where kA,kS0, and kΦ are the k−ranks [27, 28] of the matrices A, S0, and Φ, respectively. Π is a permutation matrix. Δi, i=1,…,3, are diagonal scaling matrices satisfying Δ1Δ2Δ3=IK.If the three factor matrices have full rank, then condition (26) becomes
(28)minM,K+minN,K+min2,K≥2K+2.In this paper, we assumeN≥K,K≥2, then the identifiable condition (28) becomes M≥K. Therefore, the maximum number of sources that can be handled by the PARAFAC algorithm is no more than M.WhenX∈ℂM×N×P is a tensor that admits the Vandermonde constrained PARAFAC decomposition, we have the following uniqueness result.Theorem 2 [23, 24].
LetX∈ℂM×N×2 be a tensor with matrix representation Y¯=A⊙S0ΦT, where A∈ℂM×K is a Vandermonde matrix with distinct generators, S0∈ℝN×K and Φ∈ℂ2×K. Consider the matrix representation Y˜=AL1⊙S0AL2⊙ΦT with L1+L2=M+1. If
(29)rA¯L1⊙S0=K,rAL2⊙Φ=K,for some L1+L2=M+1, then K=rX, and the Vandermonde constrained PARAFAC decomposition of X is unique. Generically, condition (29) is satisfied if and only if
(30)minL1−1N,2L2≥K.When we assumeN≥M and L1≥2, the condition (30) becomes 2L2≥K. Therefore, the maximum number of sources that can be identified by the Vandermonde constrained PARAFAC algorithm is 2L2. If 2L2>M, the Vandermonde constrained PARAFAC method can identify more sources than the PARAFAC method under the same conditions. Notably, the scaling does not involve the Vandermonde factor matrix in a PARAFAC model with Vandermode constraint. This means that the estimated matrices A˜,S˜0, and Φ˜ are related to A, S0, and Φ via
(31)A˜=AΠ,S˜0=S0ΠΔ˜2,Φ˜=ΦΠΔ˜3,where Δ˜j,j=2,3, are diagonal scaling matrices satisfying Δ˜2Δ˜3=IK.Remark 2.
The proposed method can be seen as a generalized ESPRIT method. In line with the experiments for ESPRIT in [29, 30], we can choose the pair (L1,L2) such that the dimensions of the matrices AL1⊙S0 and AL2⊙Φ are close, with the inequality minL1−1N,2L2≥K satisfied. According to [29, 30], we know that such a direct method will yield a very good estimation result. Therefore, in order to reduce the computational complexity, L1 is fixed to 3 and L2=M−2 in this paper. Although it may not yield the best estimation result with the parameters being set as L1=3 and L2=M−2, a good estimation result can be still achieved. Note that simulations in Section 5 show that the DOA performance of such a direct approach still outperforms than other methods.


## 4. Performance Analysis

### 4.1. Computational Complexity Analysis

In this subsection, we discuss the computational complexity of the proposed Vandermonde constrained PARAFAC method. The main computational tasks of the Vandermonde constrained PARAFAC method are computing the compact SVD ofY˜ which requires O2NL1L2K [31], pseudoinverse operation of matrix U1 which requires O2L1−1NK2, and the EVD of U1†U2 which requires OK3. Therefore, the main computational complexity of the proposed Vandermonde constrained PARAFAC method is O2NL1L2K+2L1−1NK2+K3. The TALS-based PARAFAC method is also illustrated for comparison. The computational cost of TALS-based PARAFAC method is OK3+2NMKn1 [32], where n1 is the number of iterations. We summarize the main computational complexity of the proposed Vandermonde constrained PARAFAC method and the TALS-based PARAFAC method in Table 1.Table 1
Computational complexity of the proposed algorithm and TALS-based PARAFAC.
Algorithm
Computational complexity
The proposed method
O
2
N
L
1
L
2
K
+
2
L
1
−
1
N
K
2
+
K
3
TALS-based PARAFAC
O
K
3
+
2
N
M
K
n
1


### 4.2. Cramér-Rao Bound (CRB)

From [33], we can derive the Cramér-Rao Bound (CRB) of noncircular signal DOA estimation for ULA as follows:
(32)CRB=σ22NReDHΠB⊥D⊕P^T−1,where σ2 is the noise power, B=Φ⊙A=b1,…,bK,ΠB⊥=I−BBHB−1BH,P^=1/N∑n=1Ns¯0ns¯0Hn,s¯0n=s0ns0nT,s0n is the nth row of S0, D=d1θ,…,dKθ,d1ψ,…,dKψ,dkθ=∂bk/∂θk, and dkψ=∂bk/∂ψk with bk being the kth column of B.Remark 3.
Since the proposed Vandermonde constrained PARAFAC solution can be seen as a generalized ESPRIT method, hence it could be better optimized by weighting subspace fitting methods such as [34, 35].


## 4.1. Computational Complexity Analysis

In this subsection, we discuss the computational complexity of the proposed Vandermonde constrained PARAFAC method. The main computational tasks of the Vandermonde constrained PARAFAC method are computing the compact SVD ofY˜ which requires O2NL1L2K [31], pseudoinverse operation of matrix U1 which requires O2L1−1NK2, and the EVD of U1†U2 which requires OK3. Therefore, the main computational complexity of the proposed Vandermonde constrained PARAFAC method is O2NL1L2K+2L1−1NK2+K3. The TALS-based PARAFAC method is also illustrated for comparison. The computational cost of TALS-based PARAFAC method is OK3+2NMKn1 [32], where n1 is the number of iterations. We summarize the main computational complexity of the proposed Vandermonde constrained PARAFAC method and the TALS-based PARAFAC method in Table 1.Table 1
Computational complexity of the proposed algorithm and TALS-based PARAFAC.
Algorithm
Computational complexity
The proposed method
O
2
N
L
1
L
2
K
+
2
L
1
−
1
N
K
2
+
K
3
TALS-based PARAFAC
O
K
3
+
2
N
M
K
n
1


## 4.2. Cramér-Rao Bound (CRB)

From [33], we can derive the Cramér-Rao Bound (CRB) of noncircular signal DOA estimation for ULA as follows:
(32)CRB=σ22NReDHΠB⊥D⊕P^T−1,where σ2 is the noise power, B=Φ⊙A=b1,…,bK,ΠB⊥=I−BBHB−1BH,P^=1/N∑n=1Ns¯0ns¯0Hn,s¯0n=s0ns0nT,s0n is the nth row of S0, D=d1θ,…,dKθ,d1ψ,…,dKψ,dkθ=∂bk/∂θk, and dkψ=∂bk/∂ψk with bk being the kth column of B.Remark 3.
Since the proposed Vandermonde constrained PARAFAC solution can be seen as a generalized ESPRIT method, hence it could be better optimized by weighting subspace fitting methods such as [34, 35].


## 5. Simulation Examples

The simulations are presented to evaluate the DOA estimation performance of the proposed method in this section. The root mean square error (RMSE) is used to evaluate the DOA estimation performance, which is defined as(33)RMSE=1K∑k=1K1500∑q=1500θ^k,q−θk2,where θ^k,q is the estimated values of θk in the qth trial, q=1,…,500. K is the number of sources.Figure2 shows the DOA estimation result of the proposed Vandermonde constrained PARAFAC method. In this simulation, we set M=8, N=300, SNR=5dB. 30 Monte Carlo trials are carried out. Moreover, we assume that there exist K=4 noncircular signals locating at angles θ1,θ2,θ3,θ4=10°,20°,30°,40°, and the noncircular phases are set as ψ1,ψ2,ψ3,ψ4=20°,30°,50°,60°. It can be seen from Figure 2 that our Vandermonde constrained PARAFAC method has a good DOA estimation performance even at a low SNR level. Therefore, the effectiveness of the proposed method is verified.Figure 2
Angle estimation result for 30 independent trials withSNR=5dB.Figure3 depicts the DOA estimation result of the proposed Vandermonde constrained PARAFAC method for close sources. In this simulation, we assume M=8, N=300, SNR=15dB, K=4, θ1,θ2,θ3,θ4=10°,20°,22°,40°, and ψ1,ψ2,ψ3,ψ4=20°,30°,50°,60°. The number of the Monte Carlo trials is set to 30. From Figure 3, it can be observed that the proposed Vandermonde constrained PARAFAC method can separate close sources effectively. Moreover, Figures 2 and 3 present that the DOA estimation performance could be improved with SNR increasing.Figure 3
Angle estimation result of close sources for 30 independent trials,SNR=15dB.Figure4 presents the RMSE of the proposed Vandermonde constrained PARAFAC method versus SNR. We assume that K=4, θ1,θ2,θ3,θ4=5°,15°,45°,60°,ψ1,ψ2,ψ3,ψ4=10°,20°,30°,40°,M=8, and N=200. The TALS-based PARAFAC algorithm, ESPRIT algorithm in [8], NC-ESPRIT algorithm in [11], NC-RI-PM algorithm in [13], and the CRB are also illustrated for comparison in Figure 4. From Figure 4, we can see that the Vandermonde constrained PARAFAC method can achieve a better estimation performance than ESPRIT, NC-ESPRIT, and NC-RI-PM. Notably, by taking the Vandermonde structure into account, our Vandermonde constrained PARAFAC method outperforms the TALS-based PARAFAC approach.Figure 4
RMSE performance of different methods versus SNR.The bias and variance of the estimator versus SNR are shown in Figures5 and 6, respectively. We consider there exist K=4 noncircular signals locating at the angles of θ1,θ2,θ3,θ4=10°,20°,40°,60°, and the noncircular phases are set as ψ1,ψ2,ψ3,ψ4=20°,30°,40°,50°. The number of antennas is set as M=8, the snapshot N is fixed to 200, while SNR is varied from 0 to 30 dB. As expected, we can see from Figures 5 and 6 that the bias and variance of the estimator decrease as SNR increases, respectively.Figure 5
Bias of the estimator versus SNR.Figure 6
Variance of the estimator versus SNR.Table2 shows the comparison of average run time (s) between Vandermonde constrained PARAFAC method and TALS-based PARAFAC method. In this simulation, SNR is fixed at 15 dB, M=8, K=3, θ1,θ2,θ3=10°,20°,40°,ψ1,ψ2,ψ3=20°,30°,60°, while N is varied from 50 to 300. We can see that average run time of each method grows as the number of snapshots N increases. It can be also observed that the Vandermonde constrained PARAFAC method is more computationally efficient than the TALS-based PARAFAC method.Table 2
Average run time (s) of respective algorithms.
N
50
100
150
200
250
300
TALS-based PARAFAC
0.0641
0.0945
0.1120
0.1265
0.1397
0.1915
Vandermonde constrained PARAFAC
0.0068
0.0080
0.0093
0.0112
0.0143
0.0176At last, the impact of the number of sourcesK is investigated. Figure 7 presents the RMSE performance of our Vandermonde constrained PARAFAC method versus SNR with different K, when M=8 and N=300. It is observed from Figure 7 that the estimation accuracy decreases as the source number K increases.Figure 7
RMSE performance of the proposed method for different number of sourcesK (M=8, N=300).


## 6. Conclusions

The problem of DOA estimation of noncircular signals for uniform linear array was considered in this paper. Due to the property of noncircular signals, we constructed an extended matrix by concatenating the received data and its conjugated component. And then, a Vandermonde constrained PARAFAC model was derived from the extended matrix. Finally, estimation of the angle information was obtained via eigenvalue decomposition of two submatrices. Compared with some existing methods, the proposed method consistently provided a better DOA estimation performance. Simulation result illustrates the effectiveness of our method in terms of estimation accuracy and computational complexity. Utilizing the weighting subspace fitting to optimize the proposed method will be carried out in our future work. The case of mixed noncircular and circular signals will be also a topic for future research.



---

*Source: 4612583-2018-01-15.xml*