0 Introduction

With regard to the direction⁃of⁃arrival (DOA) estimation using sensor arrays, much work has been developed for decades, among which the subspace⁃based algorithms, e.g., multiple signal classification (MUSIC)^[1], estimation of signal parameters via rotational invariance techniques (ESPRIT)^[2] and pro⁃pagator method (PM)^[3,4] are popular due to the super resolution capability they can provide. MUSIC^[1] exploits the orthogonality between the array manifold and the noise subspace, and can be applied in the arbitrary array geometry. In MUSIC, the DOA estimation is carried out in form of a spectral search. Based on MUSIC, root⁃MUSIC was proposed in Ref.[5] with the benefits of lower complexity and improved accuracy since it estimates DOAs through a polynomial rooting technique instead of the spectral search. The ESPRIT algorithm exploits a different subspace, called as the signal subspace, and can reduce the complexity significantly as it exploits the rotational invariance property for closed⁃form solutions by limiting the sensor array to be shift invariant. Based on ESPRIT, some variants^[6,7,8,9] have been proposed. For example, the generalized ESPRIT (GESPRIT)^[6] extends the original ESPRIT to the arbitrary array case, and the unitary ESPRIT^[7] and beamspace ESPRIT^[9] reduce the computational complexity further. However, these algorithms require the eigen⁃decomposition of the array covariance matrix to compute the subspace. The complexity of calculating the eigen⁃decomposition for array covariance matrix with an M ⁃element array and L snapshots is in the order of O(M3+M2L) . Consequently, when the number of sensors is large, the computational complexity of these algorithms becomes rather extensive.

PM^[3,4] can dispense with the high⁃complexity eigen⁃decomposition since it exploits the propagator, which is a linear operator that can be easily estimated from the received data, to compute the subspace. As shown in Ref.[3], the computational complexity for estimating the propagator was O(MLK) , where K is the number of incident sources. The low⁃complexity characteristic of PM makes it potential for real⁃time and large⁃scale signal processing, e.g., the upcoming massive multi⁃input⁃multi⁃output (MIMO) system^{[10,11,12,13,14]}. Till now, many variants^{[15,16,17,18,19,20,21,22]} of the developed PM can generally be categorized into two types. The first type^[15,16,17] is based on the OPM algorithm^[3], which utilizes the propagator to obtain the noise subspace and estimates DOAs through a similar spectrum search like MUSIC^[1]. The other one^{[18,19,20,21,22]} is based on the method called as Tayem’s PM^[4], which limits the array to be shift invariant, utilizes the propagator to compute the signal subspace and obtains DOA estimates in a same manner as ESPRIT^[2].

For enhancement of DOA estimation accuracy, the utilization of the inherent characteristics of the incoming signals for the estimation accuracy enhancement has aroused considerable attention. Among much work of this type^{[23,24,25,26,27]}, researches^[25,26,27] based on non⁃circularity became popular since the non⁃circular signals have been extensively used in actual communication systems, e.g., binary phase shift keying (BPSK) and amplitude modulation (AM) signals^{[25,26,27,28]}. The notion of circularity stems from the geometrical interpretation of complex random variables. Non⁃circular signal refers to the signal whose elliptic covariance matrix is not equal to zero. The extension of conventional algorithms for non⁃circular source has drawn considerable attentions. For example, Refs.[25,26,27] extended MUSIC^[1], root⁃MUSIC^[5] and ESPRIT^[2], respectively. In Ref.[29], a sparse representation⁃based algorithm was proposed, while its main weakness is the high complexity. In addition to various DOA estimation algorithms, the corresponding performance analysis problem has also been addressed. In Refs.[30,31], the deterministic and stochastic Cramer⁃Rao lower bounds (CRLB) were proposed for DOA estimation of non⁃circular sources, respectively.

This paper proposes a PM⁃based algorithm for DOA estimation of non⁃circular signals using an arbitrary array. We firstly reconstruct the received signal by considering the signal and its conjugate simultaneously. Then the proposed algorithm computes the signal subspace using the propagator and estimates the DOAs through searching a 1⁃D spectrum. Moreover, a computationally⁃efficient search⁃free method is also designed using the polynomial rooting technique. The main contributions of the proposed algorithm can be summarized as follows:

(1) The proposed algorithm develops a new variant of PM unlike the two conventional variants^[3]. It exploits the propagator to compute the same signal subspace as Tayem’s PM^[4]. However, the algorithm achieves DOA estimates by introducing the rank reduction criterion newly. The rank reduction has been considered in several work^{[6,32,33,34,35]} for expanding the application scope of conventional algorithms. For example, Ref.[32] extended the conventional MUSIC to the partly calibrated array case, and Ref.[33] extended OPM^[3] to DOA estimation in the presence of the mutual coupling. With the rank reduction criterion, the proposed algorithm removes the limitation in Tayem’s PM^[4] and can be applied to arbitrary array geometry.

(2) When compared to conventional methods, e.g., ESPRIT^[2], GESPRIT^[6], OPM^[3], and Tayem’s PM^[4], the proposed algorithm can achieve much better estimation accuracy and can handle more sources than the number of sensors due to the utilization of the non⁃circular information in the elliptic covariance.

(3) The proposed algorithm can dispense with the high⁃complexity eigen⁃decomposition procedure of the correlation sample matrix.

(4) The proposed algorithm can work well without knowledge of the non⁃circular phase.

Notations: Lower⁃case (upper⁃case) boldface symbols denote vectors (matrix); the transpose, complex conjugate, Hermitian transpose, inverse of the matrix A are denoted by AT , A* , AH , A-1 , respectively; ☉ is the Schur⁃Hadamard product, respectively; j=-1 represents the imaginary unit; det{A} is the determinant of A ; ⋅ returns the maximum integer that is not bigger than the inside argument; diag{⋅} returns a square matrix with the input arguments on the main diagonal; E{⋅} denotes the statistical expectation operator; [a]m and [A]m,n are the mth and the (m,n)th elements of a and A , respectively; Im denotes the m×m identity matrix; O(a) means that the complexity of the arithmetic is linear in a∈R+ . min{⋅} is to take the minimum.

1 Data Model

1.1 Definition of non⁃circular signal

As depicted in Refs.[25,26,27,28,29,30,31,32,33,34,35,36], circularity is an important property of random variables, whose concept stems from the geometrical interpretation of complex random variables. To be more specific, for a complex random signal s , if both its mean E{s} and elliptic covariance E{ss} equal zero, s is circular. Otherwise, if E{s}=0 and E{ss}≠0 , s is a non⁃circular signal. The non⁃circular rate is defined as ρ=|E{ss}|/E{ss*} , where E{ss*} is the covariance. When ρ=1 , the signal, e.g., the AM and BPSK signals, is characterized as the maximal non⁃circularity rated signal. Meanwhile, signal with 0<ρ<1 is called as the common non⁃circularity rated signal. The classical DOA estimators^[37,38], e.g., MUSIC^[1], ESPRIT^[2] and PMs^[3,4], do not consider the non⁃circular property. When they are applied in the non⁃circular signal scenario, the non⁃circular information in the elliptic covariance E{ss} will not be fully utilized.

Throughout this work, we consider the maximal non⁃circularity rated signal only as some popular work^[25,26,27]. In order to facilitate the analysis, we will re⁃write the non⁃circular signals st according to the non⁃circularity^[25,26,27] as follows

where rt=[s1(t),⋯,sK(t)]T∈RK×1 , Ω=diag{e-jψ1,e-jψ2,⋯,e-jψK} with ψk being the non⁃circular phase of the k th signal.

1.2 Model formulation

As shown in Fig.1, we consider an arbitrary linear array consisting of M omnidirectional elements. The coordinates of sensors are denoted by dmm=1M and the first sensor is selected as the reference one, i.e., d1=0 . Suppose K far⁃field, narrow⁃band, non⁃circular source signals are impinging on the array from distinct DOAs θkk=1K . The received signals can be written in a vector form as^[6]

Fig.1 Array geometry

where t is the time index, θ=[θ1,⋯,θK] , A(θ) denotes the unknown M×K array manifold matrix, st denotes the signals transmitted by non⁃circular sources, nt is the additive white Gaussian noise with zero mean and variance σn2 . Moreover, the array manifold A(θ) is given by A(θ)=[a(θ1),…,a(θK)] , and a(θk) is the steering vector with the structure, shown as

2 Generalized Propagator Algorithm

2.1 Data construction

Under the assumption of non⁃circular sources, we can extend the observed array output by combing the original array output and its conjugate as

where Ψ=[ψ1,⋯,ψK] , B(θ,Ψ)=A(θ)ΨA*(θ)Ψ*∈C2M×K is defined as the extended array manifold matrix with the kth column shown as

From Eqs.(5), (6), it could be readily checked that additional (M-1) degrees of freedom (DOFs) can be provided for DOA estimation by considering the non⁃circular property, which will facilitate the estimation accuracy. However, due to the coupling of the non⁃circular phase, direct applying some classical algorithms, e.g., MUSIC^[1] and OPM^[3], will require a 2⁃D search over the DOAs and non⁃circular phases, which creates a heavy computational burden.

2.2 Construction of signal subspace by propagator

As aforementioned, the proposed algorithm exploits the propagator to compute the signal subspace, instead of the eigen⁃decomposition procedure in MUSIC^[1] and ESPRIT^[2]. It is necessary to make a common assumption like much previous work^[3,4] that the extended array manifold B is of full rank and the K rows of B are linearly independent. Then the other (2M-K) rows of B can be expressed as a linear combination of the former K rows. This paper assumes the first K rows of B are linearly independent. Then we partition the extended array manifold matrix B as

where B1 is a nonsingular K×K matrix containing the first K rows of B and B2 is an (2M-K)×K matrix containing the last (2M-K) rows.

Define an M×K matrix Q=[IK,P]H , where P is the so⁃called propagator and IK is a K×K identity matrix. According to Eq.(7), it is easy to obtain that^[1]

Since B1 is nonsingular, it holds that Q=BB1-1 . It follows that Q spans the same column space as B , i.e., span{Q}=span{B} , which means we have computed the signal subspace with the aid of the propagator P .

2.3 DOA estimation via rank reduction criterion

Partition the entire array into two different subarrays with the same number of sensors. With regard to the choice of the sensors in each subarray, it could be arbitrary and the number of sensors in each subarray could take value from 2 to (M-1) . In this paper, we select (M-1) sensors for the two different subarrays in order to fully utilize the received information and ensure the best accuracy. Without loss of generality, we make the first subarray contain the sensors with coordinates {d1,⋯,dM-1} and the other subarray contain the sensors with coordinates {d2,⋯,dM} . The two selection matrices for two subarrays are respectively defined as

where IM-1 is an (M-1)×(M-1) identity matrix and 0 denotes a zero vector of dimension (M-1)×1 . After data constructing in Ref.[4], the extended array manifold matrices of these two subarrays, namely U and V , are given by

where L1≜blkdiagT1,T1∈R(2M-2)×2M , L2≜blkdiagT2,T2∈R(2M-2)×2M . To be more specific, the kth column of U and V , denoted by uk and vk , are given by

According to Eqs.(11), (12), we could find a shift relationship between them, which is given by

where Γk=diag{e-j2πΔ1sinθk/λ,⋯,e-j2πΔM-1sinθk/λ,ej2πΔ1sinθk/λ,⋯,ej2πΔM-1sinθk/λ} and Δi=di+1-di,i=1,⋯,M-1 . By introducing a new (2M-2)×(2M-2) diagonal matrix Θθ as

we can form a new matrix W(θ)=Ev-ΘθEu.

It can be found that

where C=V-Θ(θ)U=[(Γ1-Θθ)u1,⋯,

According to Eq.(16), when θ=θk , the kth column of C , i.e., (Γk-Θθ)uk , will become zero. In this case, if K≤(2M-2) , since B1 is nonsingular, the matrix W(θ) will drop rank. It follows that the determinant of GHW(θ) will equal to zero, where G is an arbitrary (2M-2)×K full⁃rank matrix. The choice of G has been discussed in Refs.[6,8]. In this work, we adopt the choice in Ref.[8] and let G=WH(θ) . Hence, the following spectral function can be utilized for DOA estimation

Remark 1 From Eq.(17), it can be seen that the proposed algorithm can estimate DOAs with no estimation nor prior knowledge of the non⁃circular phases.

Remark 2 The proposed algorithm can estimate up to (2M-2) sources due to the extension of the virtual array aperture, whereas the classical MUSIC^[1,2] and PM^[3,4] can only estimate (M-1) sources.

2.4 Search⁃free implementation of the proposed algorithm using polynomial rooting

The proposed algorithm has advantage of rather low complexity when compared with some conventional non⁃circular estimators, e.g. NC⁃MUSIC^[25], by using the propagator instead of the eigen⁃decomposition procedure. However, when the array configuration satisfies certain conditions, the computational complexity can be further reduced by using the polynomial rooting technique.

Without loss of generality, we assume 0<Δ1≤⋯≤ΔM-1 . Define z≜e-j2πΔ1sinθ/λ and pi=Δi/Δ1,i=1,⋯,M-1 . Then Eq.(15) can be rewritten as

If the array geometry satisfies

DOAs can be estimated by rooting the polynomial

where W(z)=Ev-ΘzEu and F(z)=EvH-ΘTz-1EuH .

Although there may exist more than K roots for polynomial, we can simply take K roots that are closest to the unit circle as our final results, similarly to root⁃MUSIC^[5]. Then the DOA estimates can be obtained from the phases of these roots. The highest order of Eq.(20) is proved to be 2∑n=2M-K-12M-2qn , where qn=pn+12,n=1,⋯,2M-2 , · returns the maximum integer that is not bigger than the inside argument. The proof is shown in Appendix A for detail. The complexity of finding roots of Eq.(20) is O2∑n=2M-K-12M-2qn3 ^[39], which is much smaller than searching the 1⁃D spectrum Eq.(17). The specific complexity analysis will be given in Section 3.

The array geometry for the proposed rooting method should satisfy that all pii=1M-1 are integers, however, it still applies to a much more general class of array geometries than some classical algorithm based on the rotational invariance property, e.g., ESPRIT^[2], NC⁃ESPRIT^[7] and Tayem’s PM^[4].

2.5 Actual implementation

The propagator P in practice could be estimated directly from the received data or from the sample covariance matrix^[4]. Both two estimation methods of P need to collect a series of snapshots of the received signal as

where L denotes the number of snapshots.

The sample covariance matrix Rˆ is given by

According to Ref.[3], by introducing the following partition of the received data and the sample covariance matrix as

where Y1∈CK×L and Y2∈C(2M-K)×L contain the first K rows and the last (2M-K) rows of Y , Rˆ1∈C2M×K and Rˆ2∈C2M×(2M-K) contain the leftmost K and rightmost (2M-K) columns of Rˆ , we could estimate the propagator Pˆ via

Using Eq.(22) and Eq.(25) (or Eq.(26)), the actual version of the proposed algorithm is given by

where Wˆ(θ)=Eˆv-ΘθEˆu , Wˆ(z)=Eˆv-ΘzEˆu , Fˆ(z)=EˆvH-ΘTz-1EˆuH , Eu=L1Qˆ , Ev=L2Qˆ and Qˆ=[IK,Pˆ]H .

The major steps of the proposed estimator are concluded as follows:

(1) Rearrange the received signal as Eq.(4) and select two subarrays.

(2) Estimate the propagator from the received data via Eq.(25) or via the covariance matrix Eq.(22) as shown in Eq.(26). Compute the corresponding signal subspaces as Eq.(13).

(3) Estimate the DOAs through the spectral peak search over Eq.(27) or rooting Eq.(28) when the array geometry satisfies the condition Eq.(19).

3 Complexity Analysis

This paper mainly considers the complex multiplication operation, which costs the most complexity.

In the proposed algorithm, if we choose to estimate the propagator directly from the received data Eq.(25), it will cost complexity of O((2MK+K2)L+K3) , whereas if we estimate from the sample covariance, it will cost complexity of O((4L+4K)M2+2K2M+K3) . In contrast, for computing the same subspace, the eigen⁃decomposition procedure will require complexity of O(8M3+4LM2) , which is typically larger than the complexity of the propagator⁃based works, especially when we directly estimate the propagator from the received data. Let α denote the number of search time for DOA estimation. In the proposed 1⁃D search algorithm, the complexity of finding the DOA estimates is O(α((2M-2)K2 , (2M-2)K+K3)) ，whereas it is O2∑n=2M-K-12M-2qn3 in the proposed rooting algorithm.

Table1concisely lists the complexity of the proposed algorithm, ESPRIT^[2], NC⁃ESPRIT^[7], GESPRIT^[6] and OPM^[3]. The reason we choose GESPRIT and OPM for comparison is that they can also be applied to arbitrary array geometry, whereas the reason for ESPRIT^[2] and NC⁃ESPRIT^[7] is they are classical low⁃complexity DOA estimators. For intuitive illustration, we consider K=3 non⁃circular source signals impinge on a uniform linear array (ULA) of M=12 sensors. In such a case, the highest order of Eq.(28) is 6. The other parameters are set as L=200,α=1800 .

The corresponding complexity of each algorithm is plotted in Fig.2, where the notation “Proposed⁃search (data)” represents the complexity of the proposed 1⁃D search algorithm in which the propagator is directly estimated from the received data as Eq.(25), whereas the notation “Proposed⁃search (covariance)” denotes the complexity of the proposed 1⁃D search algorithm with the propagator estimated from the sample covariance as Eq.(26). The notation rule is similar to the proposed rooting algorithm and OPM^[3]. It can be clearly seen from Fig.2that because of avoiding the spectral search, the proposed rooting algorithm has an obvious advantage of low complexity especially when the propagator is directly from the received data, and its complexity could be even close to ESPRIT and NC⁃ESPRIT, which also have low complex because they can provide closed⁃form solutions. With regard to the complexity of the algorithms that can be applied to arbitrary array geometry, the proposed 1⁃D search algorithm has close computation cost compared with GESPRIT^[6] and the conventional OPM^[3] has the lowest computational efficiency. However, the proposed algorithm enlarges the virtual array aperture by considering the non⁃circular property and exhibits better estimation accuracy, which will be proven in the following simulations in Section 4.

Fig.2 Complexity comparison when M=12,K=3,L=200,α=1800

4 Simulations

This section provides numerical simulations to demonstrate the performance of the proposed algorithm. Root mean square error (RMSE) is adopted as a measurement of numerous estimators, shown as

where P is the number of Monte⁃Carlo simulations, βˆk,p the estimate of a parameter β of the pth trial corresponding to the kth source. The signal⁃to⁃noise ratio (SNR) is defined as SNR=σs2/σn2 , where σs2 is the power of incoming signals. Totally 500 Monte⁃Carlo runs are used for average. Unless otherwise stated, the number of snapshots is 200, the proposed searching algorithm is adopted and the propagator is estimated from the sample covariance matrix as Eq.(26) in all simulations. The search region for DOA estimation is (-90°,90°).

4.1 Arbitrary array case

We firstly examine the performance of the proposed algorithm when applied to an arbitrary array. Herein we consider a 9⁃element array with coordinates being (0,0.45λ,0.93λ,1.43λ,1.92λ,2.42λ,29λ)throughout all simulations in this subsection. Note that this array is sufficiently arbitrary and imposes no shift invariance property, which means ESPRIT^[2], NC⁃ESPRIT^[7] and Tayem’s PM^[4] cannot be applied here.

In the first example, we assume there are K=12 non⁃circular sources from [-60°,-45°,-35°, -20 °,-10°,0°,10°,25°,35°,45°,60°,75°]. The corresponding non⁃circular phases are [0°,5°,10°,15°,20°,25°,30°,35°,40°,45°,50°,55°]. The spectrum of the proposed algorithm is shown in Fig.3, where SNR is set as 15 dB. It could be seen that clear peaks are formed around the theoretical directions. Fig.4shows the 100 estimation results of SNR=15 dB and 25 dB, respectively, for comprehensive presentation. From Fig.4, we can find that the proposed algorithm can obtain efficient DOA estimates of 12 non⁃circular sources in the 7⁃element arbitrary array case. With the increase of SNR, the estimation results become more accurate. From Figs.3, 4, it is confirmed that the maximum number of sources the proposed algorithm can estimate is (2M-2) , which is twice than that of conventional algorithms without utilization of the non⁃circular property^[1,2,3,4].

Fig.3 Spectrum of the proposed 1⁃D search

Fig.4 Estimation results of 100 Monte⁃Carlo trials of the proposed 1⁃D search algorithm

In the second example, we compare the estimation accuracy of the proposed algorithm with GESPRIT^[6] and OPM^[3]. It is herein assumed that there are K=3 non⁃circular sources from [10°,30°,50°] with non⁃circular sources of [10°,25°,40°]. RMSEs versus SNR for each algorithm are illustrated in Fig.5. The CRLB for DOA estimation of non⁃circular sources^[30] is also displayed as a benchmark. It is seen that the proposed algorithm significantly outperforms GESPRIT^[6] and OPM^[3] since the proposed algorithm fully exploits the non⁃circular information in the elliptic covariance.

Fig.5 RMSEs versus SNR of the proposed algorithm, GESPRIT^[6] and OPM^[3]

In the third example, we illustrate the performance of the proposed algorithm, GESPRIT^[6] and OPM^[3] with various number of snapshots. It is assumed that the number of snapshots increases from 100 to 500, whereas the other parameters are the same as that in the former example. Fig.6displays RMSEs of the three algorithms versus various number of snapshots as well as CRLB at SNR=10 dB and 20 dB, respectively. It can be seen that the increase of the number of snapshots improves the estimation accuracy of all the three algorithms since more sample information can be utilized. Moreover, the proposed algorithm has better estimation performance than both GESPRIT^[6] and OPM^[3].

Fig.6 RMSEs versus the number of snapshots of the proposed algorithm, GESPRIT^[6] and OPM^[3]

4.2 ULA case

ULA scenario is considered in the section to demonstrate the performance of the proposed algorithm, especially the performance of the proposed rooting algorithm. We consider a ULA of M=8 sensors on which there are K=3 non⁃circular source from [10°,23°,35°] impinging. The non⁃circular phases are [10°,30°,50°] correspondingly. The inter⁃element spacing of ULA is half of the wavelength.

In the first example, we compare the RMSE performance of the proposed 1⁃D search algorithm and rooting algorithm with ESPRIT^[2], Tayem’s PM^[4], and NC⁃ESPRIT^[7]. RMSEs versus SNR of the proposed 1⁃D search algorithm, the proposed rooting algorithm and the other algorithms are displayed in Fig.7as well as CRLB. It can be found that the proposed algorithm performs slightly better than the proposed 1⁃D search algorithm since it maximizes the same cost (Eq.(27)) as the 1⁃D search algorithm in a search⁃free way, while the accuracy of the 1⁃D search algorithm depends on the size of the search grid. Meanwhile, the proposed algorithms have better estimation performance than Tayem’s PM^[4] and ESPRIT^[2]. When compared with the NC⁃ESPRIT algorithm, the proposed algorithm performs rather worse in low SNRs and has close performance under high SNRs. This phenomenon can be explained with the analysis in Ref.[40] that using propagator to compute the signal subspace is less robust than the eigen⁃decomposition of the sample covariance matrix in Ref.[7]. However, the eigen⁃decomposition procedure is less computationally efficient and the NC⁃ESPRIT^[7] could not be applied to arbitrary array geometry.

Fig. 7 RMSEs versus SNR of the proposed algorithm, ESPRIT^[2], Tayem’s PM^[4] and NC⁃ESPRIT^[27]

In the second example, we vary the number of sensors and compare the estimation performance of several algorithms, similarly as the former example. We make the number of sensors increase from 7 to 11 and the inter⁃element spacing is always half of the wavelength. RMSEs versus the number of sensors for each algorithm at SNR=10 dB and 20 dB are plotted in Fig.8as well as CRLB. It is demonstrated that RMSEs of all the algorithms degrade with the increase of the number of sensors. This is reasonable since the diversity gain can be achieved with the increase of M . In a relatively low SNR condition of 10 dB, Tayem’s PM still performs worse than ESPRIT^[2], whereas the proposed algorithm can perform close to the NC⁃ESPRIT^[7] due to the utilization of the non⁃circular property. In a relatively high SNR condition of 20 dB, all the PM⁃based algorithms have close estimation accuracy to the ESPRIT algorithms^[2,27]. Meanwhile, the proposed rooting algorithm always performs slightly better than the proposed 1⁃D search algorithm when the number of sensors varies.

Fig.8 RMSEs versus the number of snapshots

In the third example, we vary the number of non⁃circular sources and plot RMSEs of the proposed algorithm versus SNR in Fig.9. In the single source case, DOA is set as 10°. In the two⁃source case, DOAs are [10°,30°]. And in the three⁃source case, DOAs are [10°,30°,50°]. From Fig.9, it could be observed that the estimation results of the proposed algorithm are more accurate with the increase of the number of the non⁃circular sources. It is reasonable since the more sources exist, the more unknowns the proposed algorithm will handle.

Fig.9 RMSEs versus SNR of the proposed algorithm

5 Conclusions

In this paper, we propose a propagator⁃based algorithm for DOA estimation of non⁃circular sources using arbitrary array geometry. The proposed algorithm exploits the propagator for computing the signal subspace and estimates DOA through a newly⁃formed 1⁃D spectral search by introducing the rank reduction criterion. To further reduce the complexity, we also propose a polynomial rooting⁃based algorithm which avoids the spectral search. So, it avoids the high⁃complexity eigen⁃decomposition procedure and works well without information of the non⁃circular sources. By considering the non⁃circular property, the proposed algorithm can estimate twice more sources than the conventional subspace⁃based algorithms, e.g., MUSIC^[1], ESPRIT^[2], GESPRIT^[6] and OPM^[3], and have much better estimation accuracy.

Appendix

Prove that the highest polynomial order of the polynomial Eq.(20) is 2∑n=2M-K-12M-2qn , where qn=pn+12,n=1,⋯,2M-2 and [p1,⋯,pM-1] is in ascending order.

The expression of W(z) can be written as

where ηi,k and ϵi,k(i=1,⋯,M-1,k=1,⋯,K) are complex constants.

Perform row exchange operation on W(z) and obtain W¯(z) as

Similarly, F(z) and F¯(z) can be written as

It should be noted that W¯(z) and F¯(z) are obtained by performing an even number of row exchange operations on W(z) and F(z) . Hence, it follows that

We then partition W¯(z) and F¯(z) as

where W¯1(z)∈C(2M-2-K)×K and W¯2(z)∈CK×K contain the top (2M-2-K) and the upper K rows of W¯(z) . Meanwhile, F¯1(z)∈CK×(2M-2-K) and F¯2(z)∈CK×K contain the right (2M-2-K) columns and the left K columns of F¯(z) , respectively.

The following equation holds

Let R{⋅} denote an operator that returns the highest order of the input polynomial. For each entry in F¯1(z)W¯1(z) and F¯2(z)W¯2(z) , it can be observed that

where u,v=1,⋯,K .

Note that {pi}i=1M-1 is in ascending order, hence q2M-2 is not less than q2M-1-K . It then follows that

Hence, when we compute the determinant of F¯(z)W¯(z) , the entries in F¯2(z)W¯2(z) will play a major role and the following equation holds