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

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.


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][2][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][5][6][7][8][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][12][13][14].NC ESPRIT algorithm and NC Unitary ESPRIT algorithm for DOA estimation were presented in [11,12], respectively.The work [13] developed a lowcomplexity 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 Section 2, 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.I K stands for a K × K identity matrix.diag a is the diagonal matrix with diagonal elements given by a. ⊙, ∘, and ⊕ represent the Khatri-Rao product, outer product, and Hadamard product, respectively.r A is the rank of a matrix A.

System Model
Let us consider a uniform linear array consisting of M 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 x n ∈ ℂ M×1 can be modeled as where is the array manifold matrix with size of is the signal vector at snapshot n. v n ∈ ℂ M×1 stands for the additive Gaussian noise vector with zero-mean and covariance σ 2 I M .From (2), we know that A is a Vandermonde matrix with distinct nonzero generators α 1 , … , α K In this paper, the strictly secondorder noncircular signals [18] are considered.Thus, we have It is assumed that the DOAs of the sources are constant in N snapshots.We collect N snapshots corresponding to x n ; the received signal matrix X ∈ ℂ M×N is given by where S = s 1 … s N T ∈ ℂ 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 = S 0 Ψ [18], where S 0 ∈ ℝ N×K is a real-valued matrix, Ψ = diag e −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

Estimation of DOA via Vandermonde Constrained PARAFAC Decomposition
3.1.PARAFAC Model.Since the noncircular signals are considered, we construct the extended data matrix X ∈ ℂ 2M×N as follows [12,14]:  International Journal of Antennas and Propagation where J ∈ ℝ M × M is the exchange matrix with ones on its antidiagonal and zeros elsewhere.V 1 is the extended noise matrix.Λ = diag e j 2π/λ M−1 d sin θ 1 , … , e j 2π/λ M−1 dsin θ K ∈ ℂ K×K , and is a 2 × K matrix.D i Φ 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 ∈ ℂ M × N × P , whose m, n, p − th element is given by where m = 1, … , M, n = 1, … , N, and p = 1, 2 x m,n,p represents the typical element of the tensor .a m,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 as follows: 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.
where a k , s 0k , and ϕ k represent the kth columns of the matrices A, S 0 , 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 S 0 using the following least squares (LS) fitting: where Ŝ0 stands for the estimation of S 0 .In the same way, while A and S 0 are fixed, Φ can be updated by the following LS fitting: where Φ stands for the estimation of Φ.Similarly, we obtain an updated A via LS fitting as where Â stands for the estimation of A. Repeat updating process of the three matrices until convergence, and then the estimation matrices Â, Ŝ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.

DOA Estimation with
Vandermonde Constrained PARAFAC Method.From the previous analysis, we know that the matrix A is a Vandermonde matrix with distinct nonzero generators α 1 , … , α K Therefore, by taking the Vandermonde factor into account, the model 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 Φ = min 2, K If K > 2, we can conclude that the matrix Φ does not have fullcolumn 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: Denote where L 2 is the smoothing factor.The integers L 1 and L 2 must be chosen such that M = L 1 + L 2 − 1 Taking the advantage of the Vandermonde structure of A, we have where Ω = diag e j 2π/λ dsin θ 1 , … , e j 2π/λ dsin θ Therefore, ( 14) can be rewritten as In this way, the matrix representation Y of the tensor , which involves a rank-deficient matrix Φ (while K > 2), is replaced by a matrix representation Y which only involves two full-column rank matrices A L 1 ⊙ S 0 and A L 2 ⊙ Φ The matrix Y can be seen as a matrix representation of a constrained PARAFAC decomposition model [22][23][24] with three factors A L 1 , S 0 and A L 2 ⊙ Φ When the matrix Y is obtained, a more efficient algorithm can be used to fit the Vandermonde constrained PARAFAC model [23,24].
While the matrices A L 1 ⊙ S 0 and A L 2 ⊙ Φ have fullcolumn rank, we perform the compact singular value decomposition (SVD) of Y, which is given by Y ≈ UΣV H , where Σ ∈ ℂ K × K is a diagonal matrix consisting of the largest K nonzero singular values.U ∈ ℂ L 1 N × K and V ∈ ℂ 2L 2 × K are matrices composed of left and right singular vectors associated with the largest K nonzero singular values, respectively.Since the factor matrices A L 1 ⊙ S 0 and A L 2 ⊙ Φ have fullcolumn 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 A L 1 ⊙ S 0 , so we have Therefore, there exists a unique nonsingular matrix Using the inherent Vandermonde structure of A L 1 , we define the following two submatrices: where , and the relationship between A L 1 , A L 1 , and is that Therefore, there exists a diagonal matrix According to (19), ( 20), ( 21) and ( 22), we can obtain that From (23), we have U 2 = U 1 MZM −1 = U 1 Ẑ, where Ẑ = MZM −1 Assuming that the matrix A L 1 ⊙ S 0 has fullcolumn rank, we can conclude that the matrices U 1 and U 2 have full-column rank and Ẑ can be rewritten as Therefore, the generators z k , k = 1, … , K, can be obtained from the eigenvalue decomposition (EVD) of U † 1 U 2 .Then, the DOA of the kth source θ k can be recovered via 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.
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].

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, 4 International Journal of Antennas and Propagation we discuss the identifiabilitity condition for the proposed PARAFAC method.
then the three factor matrices A, S 0 , and Φ are unique up to permutation and scaling of columns, meaning that any other triple ( Â, Ŝ0 , Φ) is related to (A, S 0 , Φ) via where k A , k S 0 , and k Φ are the k − ranks [27,28] of the matrices A, S 0 , and Φ, respectively.Π is a permutation matrix.Δ i , i = 1, … , 3, are diagonal scaling matrices satisfying Δ 1 Δ 2 Δ 3 = I K If the three factor matrices have full rank, then condition (26) becomes min M, K + min N, K + min 2, K ≥ 2K + 2 28 In this paper, we assume N ≥ 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.
When ∈ ℂ M × N × P is a tensor that admits the Vandermonde constrained PARAFAC decomposition, we have the following uniqueness result.
Theorem 2 [23,24].Let ∈ ℂ M×N×2 be a tensor with matrix representation Y = ⊙ S 0 Φ T , where A ∈ ℂ M × K is a Vandermonde matrix with distinct generators, S 0 ∈ ℝ N × K and Φ ∈ ℂ 2 × K .Consider the matrix representation , and the Vandermonde constrained PARAFAC decomposition of unique.Generically, condition ( 29) is satisfied if and only if When we assume N ≥ M and L 1 ≥ 2, the condition (30) becomes 2L 2 ≥ K. Therefore, the maximum number of sources that can be identified by the Vandermonde constrained PARAFAC algorithm is 2L 2 .If 2L 2 > 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, S 0 , and Φ via where Δ j , j = 2, 3, are diagonal scaling matrices satisfying 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 (L 1 , L 2 ) such that the dimensions of the matrices A L 1 ⊙ S 0 and A L 2 ⊙ Φ are close, with the inequality min 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, L 1 is fixed to 3 and L 2 = M − 2 in this paper.Although it may not yield the best estimation result with the parameters being set as L 1 = 3 and L 2 = 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.

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 of Step 1. Construct the tensor model with the matrix representation Y = A ⊙ S 0 Φ T by exploiting the property of NC signals.
Step 2. Choose the pair (L 1 , L 2 ) subject to L 1 + L 2 = M + 1 and construct the matrix representation Algorithm 1: The Vandermonde constrained PARAFAC algorithm for NC signals DOA estimation.5 International Journal of Antennas and Propagation NK 2 + K 3 The TALS-based PARAFAC method is also illustrated for comparison.The computational cost of TALS-based PARAFAC method is O K 3 + 2NMK n 1 [32], where n 1 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.

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: where σ 2 is the noise power, 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].

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 where θk,q is the estimated values of θ k in the qth trial, q = 1, … , 500.K is the number of sources.
The bias and variance of the estimator versus SNR are shown in Figures 5 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.
Table 2 shows the comparison of average run time (s) between Vandermonde constrained PARAFAC method

Algorithm Computational complexity
The proposed method International Journal of Antennas and Propagation 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 TALSbased PARAFAC method.At last, the impact of the number of sources K 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.

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.

Figure 1 :
Figure 1: Illustration of the array geometry.

2
and d kψ = ∂b k /∂ψ k with b k being the kth column of B.

Figure 6 :
Figure 6: Variance of the estimator versus SNR.

Figure 5 :
Figure 5: Bias of the estimator versus SNR.

Figure 3 :
Figure 3: Angle estimation result of close sources for 30 independent trials, SNR = 15 dB.

Figure 4 :
Figure 4: RMSE performance of different methods versus SNR.

Figure 7 :
Figure 7: RMSE performance of the proposed method for different number of sources K (M = 8, N = 300).

8
International Journal of Antennas and Propagation

Table 1 :
Computational complexity of the proposed algorithm and TALS-based PARAFAC.

Table 2 :
Average run time (s) of respective algorithms.