Circuits Syst Signal Process DOI 10.1007/s00034-015-9987-6

A Novel ESPRIT-Based Algorithm for DOA Estimation with Distributed Subarray Antenna Yan Ma · Baixiao Chen · Minglei Yang · Yi Wang

Received: 12 July 2014 / Revised: 22 January 2015 / Accepted: 23 January 2015 © Springer Science+Business Media New York 2015

Abstract Distributed subarray antennas (DSAs), which are sparse arrays consisting of two or more far separated subarrays, have many advantages over uniform linear array, especially for direction of arrival (DOA) estimation. However, they are subject to manifold ambiguity, which has significant influence on DOA estimation. In order to solve the manifold ambiguity of uncorrelated sources for DSA, a novel method based on estimation of signal parameters via rotational invariance technique by utilizing optimal subarray partition is proposed in this paper. In the proposed method, an optimized reference estimation is obtained by the rotational invariance between the new subarrays with optimal partition of DSA. The high accuracy and unambiguous DOA estimations are then disambiguated easily according to the optimized reference estimation. In this way, the performance of disambiguated DOA estimation can be enhanced in cases of the low signal-to-noise ratio and the large spacing between the subarrays. Moreover, the ambiguity threshold effect of DOA estimation for DSA is analyzed by means of the maximum a posteriori estimator. Computer simulation results show good performance of the proposed method and the effectiveness of the ambiguity threshold for DSA.

Y. Ma (B) · B. Chen · M. Yang · Y. Wang National Laboratory of Radar Signal Processing, Xidian University, Xi’an, Shaanxi 710071, China e-mail: [email protected] B. Chen e-mail: [email protected] M. Yang e-mail: [email protected] Y. Wang e-mail: [email protected]

Circuits Syst Signal Process

Keywords Distributed subarray antenna (DSA) · Estimation of signal parameters via rotational invariance technique (ESPRIT) · Optimal subarray partition · Ambiguity threshold · Maximum a posteriori (MAP) estimator

1 Introduction The estimation of direction of arrival (DOA) is an important research field in array signal processing, which has wide applications in radar, wireless communication, sonar etc. The array aperture is closely related to the performance of DOA estimation [5]. Larger effective array aperture usually produces more accurate DOA [6]. The effective array aperture can be extended by increasing the number of sensors, but the cost will be high in terms of both the system hardware and the computational complexity. In order to increase the effective array aperture with low number of sensors, the distributed subarray antenna (DSA), which is a sparse array consisting of two or more far separated subarrays, has been proposed in [21]. Lincoln laboratory also proposed coherent distributed apertures for next generation ballistic missile defence (BMD) radar systems. However, the beam pattern of DSA has grating lobes or high sidelobes. It will cause ambiguities in DOA estimation [9,35]. This leads to a performance degradation of DOA estimation [25]. In order to avoid or solve the ambiguities, many methods have been proposed. They can be classified into three kinds. The first kind is based on the so-called virtual array element construction approaches [8,10,15,17,34]. The basic idea is to reduce the sensors spacing of original antenna arrays by interpolating virtual sensors. For instance, a semi-parametric algorithm named minimum weighted norm (MWN) was adopted to perform interpolation between the subarrays [34]. MWN can also be used to extrapolate outside the subarrays to further enhance angular resolution [8]. In [17], a modified array interpolation approach for coherent source localization using DSA is proposed. A major difficulty of these methods is that the computation will remarkably increase with the increase in spacing between the subarrays of DSA. In the second class of methods, the original array is changed by sliding the positions of the sensors or by adding additional medium substrates [3,13,19,20,28,29]. A multiple signal classification (MUSIC)-based algorithm was proposed for resolving ambiguities via sliding elements position in [28]. The authors in [13] proposed a method of setting planar substrates at the front ends of the elements in a sparse array. It breaks down the linearity of the steering vectors of elements, and the ambiguity can be resolved effectively and feasibly. A criterion to choose the interelement spacings for linear prediction (LP)-based algorithms and a method to resolve possible ambiguities in estimating spatial frequencies were proposed by Tufts et al. [29]. The main drawbacks of these methods are the additional hardware cost and the changed array configuration. The third type focuses on the multi-scale algorithms [1,7,14,16,18,32,33,36–39]. The main procedure is based on a combination of two parts, one is the ambiguous high accuracy estimation and the other is a reference which is used to disambiguate the ambiguous estimation. The authors in [14] used the power corresponding to the true

Circuits Syst Signal Process

and spurious DOAs as the reference to eliminate the spurious DOAs. In [32], a simple estimation of signal parameters via rotational invariance technique (ESPRIT)-based DOA estimator for subarray-based sparse arrays is proposed, which is a combination of the ESPRIT algorithm with MUSIC or conventional beamformer. Wong et al. [36] proposed a dual-size ESPRIT (DS-ESPRIT) algorithm for resolving ambiguities. A unitary dual-resolution ESPRIT (U-DR-ESPRIT) algorithm was proposed in [37], and the pairing between the reference and the ambiguous estimation was automatic in this algorithm. These methods will not change the array configuration, and the computation will not increase with the increase in spacing between the subarrays of DSA. Unfortunately, in the disambiguation procedure of these methods, the accuracy of reference may not meet the disambiguation requirements when the signalto-noise ratio (SNR) is low and the spacing between the subarrays of DSA is large. It leads to the DOA estimation performance degradation of these methods in the above case. In this paper, based on the multi-scale algorithms, a novel ESPRIT-based method by utilizing optimal subarray partition is proposed for DSA DOA estimation. The main idea of the present ESPRIT-based method for DSA DOA estimation is as follows. A low accuracy but unambiguous DOA estimation (coarse estimation) is obtained by means of the ESPRIT algorithm. In a similar way, a high accuracy but cyclically ambiguous DOA estimation (fine estimation) can be obtained. Then, the coarse estimation serves as the reference to disambiguate the fine estimation. By the disambiguation procedure, high accuracy and unambiguous DOA estimations can be finally obtained. However, the lower accuracy of coarse estimation may cause a worse result in the disambiguation procedure when SNR is low and the spacing between the subarrays of DSA is large. For this reason, the bias of DOA estimation will greatly increase in the above case. In this paper, we use “an intermediate estimation” instead of the above coarse estimation as the reference estimation in the disambiguation procedure. The intermediate estimation can be creatively obtained by the rotational invariance between the two new subarrays with optimal partition of DSA according to the coarse estimation. And the intermediate estimation has higher accuracy than coarse estimation. In this way, the performance of DOA estimation can be enhanced in the case of lower SNR and larger subarrays spacing. Furthermore, this paper establishes the probability model of the reference estimation in the disambiguated procedure according to the maximum a posteriori (MAP) estimator [2]. Meanwhile, it combines the probability model and the theoretical mean square error (MSE) of the ESPRIT algorithm [23,24] to analyze the ambiguity threshold effect of DOA estimation of the proposed method. This threshold information is very valuable to the DSA system designer. The rest of this paper is organized as follows. The array signal model of DSA is first introduced in Sect. 2. The proposed DOA estimation method for the proposed array is presented in Sect. 3. The performance of the proposed method is analyzed, and the method of optimal partition for the subarrays within DSA is discussed in Sect. 4. Section 5 gives the complexity analysis of the proposed method. Several simulations are shown to demonstrate the effectiveness of the proposed method in Sect. 6. Conclusions are drawn in the final section.

Circuits Syst Signal Process

Fig. 1 Array configuration of the DSA

2 Array Signal Model Consider a linear DSA composed of two far separated identical subarrays as shown in Fig. 1. Both of the two subarrays are uniform linear array composed of M sensors, with an interelement spacing d ≤ λ/2, where λ denotes the signal wavelength. The distance between the phase centers of the two subarrays, named baseline, is D  λ/2. Assume that V arriving uncorrelated narrowband far-field sources impinge on this DSA with unknown DOAs (θv , v = 1, 2, . . . , V ), then (βv , v = 1, 2, . . . , V ) refer to direction cosines of these DOAs. The array manifold of DSA (a(βv ), v = 1, 2, . . . , V ) can be expressed as: T  a(βv ) = 1, e j2π(D/λ)βv T  ⊗ 1, e j2π(d/λ)βv , · · · , e j2π [(M−1)d/λ]βv , v = 1, 2, . . . , V

(1)

where ⊗ denotes Kronecker product, and the superscript T denotes the transpose. Then, the received data of the DSA can be expressed as: x(t) = As(t) + n(t), t = 1, 2, . . . , L

(2)

where x(t) = [x1 (t), x2 (t), . . . , x M (t), x M+1 (t) . . . , x2M (t)]T ∈ C 2M×1 , xm (t) is the output of the sensor, A = [a(β1 ), a(β2 ), . . . , a(βV )] ∈ C 2M×V is the array manifold matrix, the vector s(t) = [s1 (t), . . . , sV (t)]T ∈ C V ×1 represents the incident signals, and n(t) is an additive complex Gaussian white noise. The covariance matrix of x(t) is derived as   R = E x(t)xH (t)

(3)

where the superscript H denotes the Hermitian transpose. Then, the eigendecomposition of R results in R = UUH = U S  S U S H + U N  N U N H = U S  S U S H + σ 2 I

(4)

Circuits Syst Signal Process

where U = [u1 , u2 , . . . , u2M ] is an unitary matrix of 2M columns which are eigenvectors of the covariance matrix R . The first V eigenvectors of U constitute the so-called signal subspace U S = [u1 , u2 , . . . , uV ], and the last 2M − V eigenvectors constitute the noise subspace U N = [uV +1 , uV +2 , . . . , u2M ].  is a 2M × 2M diagonal matrix with elements ψ1 > ψ2 > · · · > ψV > ψV +1 = · · · = ψ2M = σ 2 , where ψv is the eigenvalue corresponding to the vth source.  S = diag(ψ1 , ψ2 , . . . , ψV ) and  N are diagonal matrices containing the signal subspace and noise subspace eigenvalues of R, respectively. Assuming that R is row full rank, then span(A) = span(U S ) , which means the relation between A and U S has the form U S = AT

(5)

where T denotes an unknown V × V nonsingular matrix. In practical situation, there are only a finite number of data samples. Then, we add ˆ = 1 L x(t)xH (t) a hat on the above notation to denote the approximation, i.e., R t=1 L ˆ S denotes the approximation of refers to the estimation of covariance matrix R, and U signal subspace U S . The others hat notations have the similar definition. 3 The Proposed Algorithm The beam pattern of DSA has grating lobes or high sidelobes because the baseline is larger than a half wavelength. It gives rise to DOA estimation ambiguities. In order to obtain high accuracy and unambiguous DOA estimation, the unambiguous coarse and ambiguous fine direction cosine estimation are first obtained in Sect. 3.1. Then, the intermediate estimation is obtained in Sect. 3.2. In Sect. 3.3, the ambiguous fine direction cosine estimation is disambiguated to obtain the unambiguous high accuracy DOA estimation. 3.1 Fine and Coarse DOA Estimation Based on ESPRIT The first and last M − 1 sensors of every subarray in Fig. 2a are used to constitute subarray C1 and C2, respectively. The invariance relation between C1 and C2 has the form (6) JC2 a(βv ) = ej2π(d/λ)βv JC1 a(βv )     where JC1 = I(2) ⊗ I(M−1) 0((M−1)×1) and JC2 = I(2) ⊗ 0((M−1)×1) I(M−1) are selection matrices. For all V sources, the rotation invariance relation in (6) can then be expressed as (7) JC2 A = JC1 AC C C

where C = diag z 1 , z 2 , . . . , z CV , z vC = ej2π dβv /λ . The invariance equation for coarse estimation can be obtained by using the relation U S = AT, as follows. JC2 U S = JC1 U S C

(8)

Circuits Syst Signal Process

(a)

(b) Fig. 2 Choices of two identical subarrays for DOA estimation by using ESPRIT. a The coarse estimation. b The fine estimation

where C = TC −1 C TC , TC denotes an unknown V × V nonsingular matrix. A similar procedure can be adapted to obtain the invariance equation for ambiguous fine estimation. The first and last subarray in Fig. 2b are used to constitute subarray F1 and F2, respectively. The invariance relation between F1 and F2 has the form JF2 a(βv ) = ej2π(D/λ)βv JF1 a(βv )

(9)

    where JF1 = I(M) 0(M×M) and JF2 = 0(M×M) I(M) are selection matrices. For all V sources, the rotation invariance relation in (9) can then be expressed as JF2 A = JF1 AF

(10)

where F = diag z 1F , z 2F , . . . , z FV , z vF = ej2π Dβv /λ . The invariance equation for fine estimation can be derived as (11) JF2 U S = JF1 U S F where F = TF −1 F TF , TF denotes an unknown V × V nonsingular matrix. Equations (8) and (11) can be solved by using least squares (LS) or total least ˆ F can be obtained by ˆ C and  square (TLS) method. When using the LS algorithm,  ˆ S for U S to solve the Eqs. (8) and (11) as follows [32] substituting U ˆ C = [(JC1 U ˆ S )H (JC1 U ˆ S )]−1 (JC1 U ˆ S )H (JC2 U ˆ S ) = (JC1 U ˆ S )+ JC2 U ˆS  ˆ F = [(JF1 U ˆ S )H (JF1 U ˆ S )]−1 (JF1 U ˆ S )H (JF2 U ˆ S ) = (JF1 U ˆ S )+ JF2 U ˆS  where the superscript + denotes Moore–Penrose inverse.

(12) (13)

Circuits Syst Signal Process

Fig. 3 Sensors partition of DSA

The coarse and fine direction cosine estimation are solved by eigendecomposition ˆ F , respectively, ˆ C and  of   arg zˆ vC , v = 1, 2, . . . , V = 2π d/λ arg(ˆz vF ) , v = 1, 2, . . . , V βvF = 2π D/λ

βvC

(14) (15)

 ˆ C denotes diagonal matrix where arg zˆ vC is the phase of zˆ vC in the range [−π, π ] ,   C ˆ , and whose diagonal elements zˆ v , v = 1, 2, . . . , V consist of the eigenvalues of   F C ˆ F denotes diagonal matrices whose diagonal elements zˆ v , v = 1, 2, . . . , V consist ˆ F. of the eigenvalues of  3.2 DOA Intermediate Estimation The sensors of DSA are partitioned as shown in Fig. 3. The first and last m sensors of every subarray are used to constitute subarrays SA1 and SA2, respectively. The choice of m is closely related to the performance of the intermediate estimation. We often choose the best m to constitute subarrays SA1 and SA2 for better performance (The method to get the best choice of m is discussed in Sect. 4.2). The array manifold relation between SA1 and SA2 can be described as: JSA2 a(βv ) = ej2π((M−m)d/λ)βv JSA1 a(βv )

(16)









where JSA1 = I(2) ⊗ I(m) 0(m×(M−m)) and JSA2 = I(2) ⊗ 0(m×(M−m)) I(m) are selection matrices. In the case of V sources, the rotation invariance relation in (17) can be described as (17) JSA2 A = JSA1 ASA

SA SA SA j2π(M−m)dβv /λ . The invariance equawhere SA = diag z 1 , z 2 , . . . , z SA V , zv = e tion for intermediate estimation can be derived as JSA2 U S = JSA1 U S SA

(18)

where SA = TSA −1 SA TSA , TSA denotes an unknown V × V nonsingular matrix. ˆ SA can be obtained by substituting U ˆ S for U S and then using LS method to solve the  Eq. (18) [32]

Circuits Syst Signal Process

ˆ SA =[(JSA1 U ˆ S )H (JSA1 U ˆ S )]−1 (JSA1 U ˆ S )H (JSA2 U ˆ S)  ˆ S )+ JSA2 U ˆS =(JSA1 U

(19)

The intermediate direction cosine estimation is solved by eigendecomposition of ˆ SA :   arg zˆ vSA SA βv = , v = 1, 2, . . . , V (20) 2π(M − m)d/λ  ˆ SA denotes diagonal matrix whose diagonal elements zˆ vSA , v = 1, 2, . . . , V where  ˆ SA . consist of the eigenvalues of 

3.3 Disambiguation The fine and intermediate direction cosine estimation is cyclically ambiguous because both the baseline D and (M − m)d are larger than a half wavelength. In order to obtain 

high accuracy and unambiguous direction cosine estimation βˆvF , v = 1, . . . , V , we  first disambiguate the intermediate direction cosine estimation βvSA , v = 1, . . . , V  by using the coarse direction cosine estimation βvC , v = 1, . . . , V as the reference.  F Then, the high accuracy but cyclically ambiguous estimation βv , v = 1, . . . , V are  using intermediate disambiguated direction cosine estimation

disambiguated by SA ˆ βv , v = 1, . . . , V as the reference. The disambiguated intermediate direction cosine estimation can be expressed as [36] λ , v = 1, 2, · · · , V (21) βˆvSA = βvSA + gvo (M − m)d       where gvo = arg min βvC − βvSA − gv λ [(M − m)d], (−1 − βvSA )(M − m)d λ ≤ gv    gv ≤ (1 − βvSA )(M − m)d λ , gv is an intermediate variable which is used to determine the coefficient gvo , c(c )denotes the largest (smallest) integer less (greater) than c, arg min denotes the argument of the minimum. The fine disambiguated direction cosine estimation can be expressed as [36] βˆvF = βvF + h ov

λ , v = 1, 2, · · · , V D

(22)

        where h ov = arg min βˆvSA − βvF − h v λ D , (−1 − βvF )D λ ≤ h v ≤ (1 − βvF ) h v   D λ , h v is an intermediate variable which is used to determine the coefficient h ov . Then, the high accuracy and unambiguous DOA estimation can be calculated as follows [36]

 (23) θˆv = sin−1 βˆvF , v = 1, 2, . . . , V

Circuits Syst Signal Process

Fig. 4 Flow chart of the proposed method

3.4 Implementation of the Proposed Algorithm As shown in Fig. 4, the proposed approach is summarized below. ˆ according to received data of the Step 1 Compute the sample covariance matrix R ˆ ˆ S. Eq. (2). Compute the eigenvectors of R to form the signal subspace U ˆ F. ˆ C and  Step 2 Calculate the LS solution of the Eqs. (8) and (11) to obtain  C Then, the unambiguous coarse estimation βv , v = 1, 2, . . . , V and ambiguous fine  estimation βvF , v = 1, 2, . . . , V of direction cosine can be obtained by the eigenvalue ˆ C and  ˆ F , respectively. decomposition of  Step 3 The sensors of DSA are optimally partitioned according to the best choice of m. The first and last m sensors of every subarray are used to constitute subarray SA1 and SA2, respectively (The best choice of m is discussed in Sect. 4.2). ˆ SA , and then, the ambiguStep 4 Calculate the LS solution of the Eq. (18) to obtain  ous intermediate estimation of direction cosine is obtained by the eigendecomposition ˆ SA . of  Step 5 Perform the disambiguation  every source, the dis procedure (21) for SA ambiguated intermediate estimation βˆv , v = 1, . . . , V are obtained. Then, use the disambiguated intermediate estimation as the reference to disambiguate the fine  ambiguous estimation βvF , v = 1, 2, . . . , V . According to the Eq. (22), the high accuracy and unambiguous estimation βˆvF , v = 1, 2, . . . , V of direction cosine is obtained. Step 6 Compute the DOA estimation according to the Eq. (23).

4 The Best Choice of m and Performance Analysis In Sect. 4.1, we first give the theoretical MSE of the DOA estimation by using ESPRIT algorithm. Then, we use the MSE to choose the best value of m in Sect. 4.2. In Sect. 4.3, we give the probability model of the reference estimation in the disambiguated procedure. And then, we combine the MSE and the probability model to analyze the DOA estimation performance.

Circuits Syst Signal Process

4.1 MSE of the ESPRIT Algorithm An error Δz vC in the eigenvalue z vC of the matrix C due to errors in the matrix C can be expressed as

Δz vC = z vC − zˆ vC ˆC C = C −  Δz vC

=

(24) (25)

qvC ΔC yvC

(26)

where yvC is the eigenvector of C corresponding to the eigenvalue z vC , i.e., C yvC = z vC yvC , and qvC is the corresponding left eigenvector, i.e., qvC C = z vC qvC . yvC and qvC also satisfy qvC yvC = 1 [26]. An approximate expression for C can be derived as follows [23]

(JC1 U S + JC1 ΔU S )(C + ΔC ) ≈ JC2 U S + JC2 ΔU S +

+

ΔC ≈ (JC1 U S ) (JC2 ΔU S ) − (JC1 U S ) (JC1 ΔU S )C

(27) (28)

ˆ S is the error matrix of the signal subspace U S . Substituting where ΔU S = U S − U C (28) into (26), Δz v can be written as   Δz vC = qvC (JC1 U S )+ (JC2 ΔU S )yvC − (JC1 ΔU S )C yvC   = qvC (JC1 U S )+ (JC2 ΔU S )yvC − z vC (JC1 ΔU S )yvC

 = qvC (JC1 U S )+ JC2 − z vC JC1 ΔU S yvC 

∗  = −z vC qvC (JC1 U S )+ JC1 − z vC JC2 ΔU S yvC

(29)

where (·)∗ denotes the conjugate of the entity inside. Then, the MSE of z vC can be expressed as 

H 

 H 

∗   C C C + C C C =qv (JC1 U S ) JC1 − z v JC2 E ΔU S yv ΔU S yv E Δz v Δz v

∗ H  H C H qv (JC1 U S )+ JC1 − z vC JC2 (30)

Circuits Syst Signal Process

Let Δuv = uv − uˆ v , then the following expressions can be obtained by means of the asymptotic properties of the error Δuv [27] ⎛   ψk ⎜ ⎜ E Δuk ΔuHj ≈ L ⎝

⎞

V  l=1 l = j

ψl (ψl − ψk )

u uH + 2 l l

σ2 (σ 2 − ψk )

⎟ H⎟ U U n n ⎠ δk j , 2

k, j = 1, 2, · · · , V   ψk ψ j E Δuk ΔuTj ≈ − u j ukT (1 − δk j ), k, j = 1, . . . , V L(ψk − ψ j )2

(31) (32)

where δk j is the Kronecker delta. Substituting γvH = qvC (JC1 U S )+ and (31) into (30) gives the following expression: 

H  C C E Δz v Δz v  =

γvH

V   

∗    C 2  yvw  JC1 − z vC JC2

w=1



∗ H   H γv JC1 − z vC JC2 E (uw − uˆ w )(uw − uˆ w ) ⎛ ⎛ V  V  

∗  ⎜   ψl  C 2 ψw ⎜ ⎜ y = γvH JC1 − z vC JC2 ⎜ u uH   vw ⎝ 2 l l L ⎝ (ψ − ψ ) l w w=1 l=1 ⎞⎞ +

σ2 (σ 2 − ψw )

l =w

∗ H ⎟⎟  H ⎟⎟ J U U − z vC JC2 γv n C1 n ⎠⎠ 2

(33)

C is the wth element of yC . Let Δβ C = β C −β , according to the relationship where yvw v v v v between βvC and z vC , the MSE of the coarse direction cosine estimation can be estimated by [24]:

E

 2  = ΔβvC

1 2π d

2 E



Δz vC

   C ∗ 2  C  C T   C H  − Re zv Δz v E Δz v Δz v 2

(34) {} where Re denotes the real-value part of the entity inside. The expression of   T C C E Δz v (Δz v ) can be derived as follows

Circuits Syst Signal Process

⎛  V  V

T 

 ⎜ C C C =γvH ⎜ J E Δz vC Δz vC y y − z J C2 C1 vw vl v ⎝ w=1 l=1 l =w

⎞T ⎞

⎛

(35)

 ⎜ ⎟ ⎟ ∗ C ⎟ ⎟ E (uw − uˆ w )(ul − uˆ l )T ⎜ ⎝JC2 − z v JC1 ⎠ ⎠ γv     ∗ 2 T Substituting (32) into (35), the expression of (z vC ) E Δz vC (Δz vC ) can be written as

      z vC

∗ 2

E Δz vC Δz vC ⎛

T

V V   ∗ 

∗ ⎜ C C C C C z = γvH ⎜ y y J − z z J C2 vw vl v v v C1 ⎝ w=1 l=1 l =w

−

ψw ψl L(ψw − ψl )

u uT 2 l w

 

⎞ z vC

∗

T ⎟

∗ ∗ JC2 − z vC z vC JC1 ⎟ ⎠ γv

⎛ V  V 

∗  ⎜ C C J = γvH ⎜ y y − z vC JC2 C1 vw vl ⎝ w=1 l=1 l =w

−

ψw ψl L(ψw − ψl )

u uT 2 l w



⎞

∗ T ⎟ ∗ JC1 − z vC JC2 ⎟ ⎠ γv

(36)



 2 can be estiThe MSE of the coarse direction cosine estimation E (ΔβvC ) mated by substituting (33) and (36) into (34). Note that the MSE of intermediate ambiguous direction cosine estimations βvSA , v = 1, . . . , V and the fine ambigu ous direction cosine estimations βvF , v = 1, . . . , V can be obtained by using the above-mentioned method. For the length limit, we here do not give the processing results. 4.2 The Best Choice of m  It can be seen from (21) that the performance of the intermediate estimation βvSA is closely related to the choice of m. According to the estimation MSE offered in

Circuits Syst Signal Process

Sect. 4.1, the best choice of m is now discussed. In order to take full advantage of the information of the sensors in DSA, we define M/2 ≤ m ≤ M − 2 [22]. To find the best choice of m, we obtain m¯ (the best choice of m) corresponding to the minimum MSE of βvSA . The expression of MSE of βvSA is similar to (34). According to the derivation of (34), m¯ can be determined by ! ! ! !  ! 2 2 ! ! ! 1 SA ! ! ! min ! E Δβm,v ! = min m m ! ! 2π(M − m)d !     SA ∗ 2  SA  SA T  !  SA H  ! SA ! − Re Δz m,v E Δz m,v Δz m,v E Δz m,v Δz m,v ! ! (37) × ! 2 ! !

SA = β SA − where m = M/2 , M/2 + 1, . . . , M − 2, v = 1, 2, . . . , V . Δβm,v m,v βv , the subscript m, v denotes the variable or vector when m = M/2 , M/2 + 1, . . . , M − 2, v = 1, 2, . . . , V . According to the derivation of (33), the expression of  SA H  SA can be derived by E Δz m,v Δz m,v



H  SA SA E Δz m,v Δz m,v  V   

H  ∗ 

 SA 2 SA SA = γm,v JSA2 ym,vw  JSA1 − z m,v w=1

  ∗ H 

H SA SA E (um,w − uˆ m,w )(um,w − uˆ m,w ) γm,v JSA1 − z m,v JSA2 ⎛ ⎛ V  V  H  ∗  ⎜

 ψl  SA 2 ψw ⎜ SA SA ⎜ y JSA1 − z m,v = γm,v JSA2 ⎜ u uH   m,vw ⎝ ⎝ 2 l l L (ψ − ψ ) l w w=1 l=1 ⎞⎞ +

l =w

∗ H

⎟⎟  H ⎟⎟ SA SA J U U − z J γm,v n SA1 SA2 n ⎠⎠ m,v (σ 2 − ψw )2 σ2

(38)

where Δz vSA = zˆ vSA − z vSA , yvSA is the eigenvector of SA corresponding to the eigenSA is the wth element of ySA , u value z vSA , yvw m,w is uw when m = M/2 , M/2 + v  SA H 1, . . . , M − 2. γm,v = qvSA (JSA1 U S )+ , qvSA is the left eigenvector of SA corresponding to the eigenvalue z vSA , i.e., qvSA SA = z vSA qvSA . According to the derivation  2  SA  SA T  SA ∗ can be derived by E Δz m,v Δz m,v of (36), the expression of z m,v

Circuits Syst Signal Process



SA z m,v

∗ 2 

T  SA SA Δz m,v E Δz m,v ⎛

V  V  ∗ 

⎜ SA H ⎜ SA SA SA JSA1 − z m,v = (γm,v ) ⎝ ym,vw ym,vl JSA2 w=1 l=1 l =w

−

ψw ψl L(ψw − ψl )

⎞⎛ ⎞∗  ∗ T ⎟ ⎜

⎟ T SA ⎟ ⎜γ SA ⎟ J u u − z J l SA1 SA2 w m,v m,v ⎠⎝ ⎠ 2

(39)

The value of m¯ can be obtained by substituting (38) and (39) into (37). However, the direction cosine (βv ) is unknown and  to be estimated in practical applications. In this case, we use the coarse estimation βvC instead of βv in (37). In other words, we use βvC to establish the signal model and then analyze the MSE of βvSA to get m¯ (the best choice of m which is obtained according to βvC ). Because of the MSEs of βvSA correspond to m¯ and m¯ are approximately equal, we take m¯ served as the best choice of m in practical applications. 4.3 Threshold Effect of DOA Estimation The beam pattern of DSA has grating lobes or high sidelobes because the baseline (D) is much larger than a half wavelength. In the case of the lower SNR and the larger spacing between the subarrays, the DOA estimation may be obtained by the “false” peaks instead of the “true” peak. The “false” peaks correspond to grating lobes or high sidelobes, and the “true” peak corresponds to the mainlobe. Then, the bias of DOA estimation will become larger when DOA estimation is obtained by the “false” peaks. This phenomenon is often called ambiguity threshold effect of DOA estimation [4,30]. The bias of DOA estimation increases with the decrease in SNR when the baseline does not change. The bias begins to increase obviously when SNR decreases to a certain point. We call this point as the ambiguity threshold of SNR [4]. The bias of DOA estimation decreases with the increase in baseline when the grating lobes do not affect the accuracy of DOA estimation. As the baseline increases up to a certain length, the bias increases rather than decreasing. We call this length as the ambiguity threshold of the baseline [4]. Therefore, the ambiguity threshold effect is closely related to the performance of DOA estimation. According to Sect. 3.3, the main idea of the disambiguation is that the 

procedure Ref ˆ to disambiguate unambiguous estimation is used as the reference estimation βv 

Amb , i.e., βvC and βˆvSA the high accuracy but cyclically ambiguous estimation βˆv are the reference estimation βˆvRef , βvSA and βvF are βˆvAmb to be disambiguated. This subsection establishes the probability model of the reference estimation according to the MAP estimator [2]. Then, the probability model and the MSE of the ESPRIT algorithm calculated in Sect. 4.1 are combined to analyze the ambiguity threshold effect of DOA estimation of the proposed method.

Circuits Syst Signal Process



Note that the error of the reference estimation ΔβˆvRef = βˆvRef − βv is zero mean Gaussian processes because the input noise n(t) is an additive complex Gaussian 

white noise [22]. The probability density function (PDF) of the reference estimation βˆvRef can be expressed as ⎛

 

 p βˆvRef  βv = "

1 e  2  2π E βˆvRef − βv

⎝−

 2E

(

1  2 βˆvRef −βv

)

⎞

2 Ref ˆ βv −βv ⎠

(40)

We know that the distance between grating lobes in direction

 cosine is uniform Ref [31], and then, the probability of the reference estimation βˆv can be written as ⎛

# Pz =

βv +(2z+1)/(2Dn ) βv +(2z−1)/(2Dn )

"



2π E

⎞

⎝−

1 βˆvRef − βv

e 2 

 2E

(

1  (β−βv )2 ⎠ 2 βˆvRef −βv

)

dβ (41)

  the grating lobes. Pz denotes the probwhere 1 Dn = λ D is the distance 

between Ref ˆ beside the zth grating lobe. P0 denotes the ability of the reference estimation βv 

Ref beside the mainlobe. probability of the reference estimation βˆv According to [12], Eq. (8.51), these probabilities are given by ⎛ ⎛ ⎞ ⎞ ⎜ ⎜ Pz = Q ⎜ ⎜ ⎝

⎜ ⎟ ⎟ ⎜ ⎟ ⎟ (2z + 1) (2z − 1) ⎟− Q⎜ ⎟ "  " ⎜ ⎟ ⎟    2 ⎠ 

2 ⎝ ⎠ Ref Ref ˆ ˆ 2Dn E βv − βv 2Dn E βv − βv #

 2 1 √ e(−u 2) du 2π x is Q function, and it can be expressed as the following approximation [30]  1 −x 2 2 , x >0 Q(x) ≈ √ e 2π x where

Q (x) =

(42)

∞

(43)

(44)

The MSE of direction cosine estimation is equal to mainlobe error when the reference estimation is beside the  mainlobe. The MSE of direction cosine estimation is approximately equal to (z Dn )2 when the reference estimation is beside the zth grating lobe [2].  Z  2   2  2 2 + 2 MSE(βˆv ) = E (βˆv − βv ) ≈ P0 E βˆvAmb − βv z Pz (45) Dn z=1

where Z is the number of grating lobes included in the approximation.

Circuits Syst Signal Process

According to [2], we find that only the mainlobe and the first grating lobe have significant contribution to the total MSE at the ambiguity threshold. From (42), we also note that P0 ≈ 1 in this region. Substituting (42) and (44) into (45) results in the following approximate expression of MSE:  2  MSE(βˆv ) ≈ E βˆvAmb − βv $ ⎛ ⎞  % 2  % 1 ⎝−  ⎠  % 8E βˆvRef − βv 2 & 1 8E (βˆvRef −βv ) Dn2 e + π Dn

(46)

In our proposed method, the disambiguated βvSA is the reference estimation βˆvRef , and βvF is the ambiguity estimation βˆvAmb   . Then,  using the MSE result calculated  2

2

in Sect. 4.1, i.e., E (ΔβvSA ) and E (ΔβvF ) , we can determine the ambiguity threshold of SNR and baseline according to (46). According to the definition of the ambiguity threshold of baseline, MSE(βˆvF ) decreases with increasing baseline when SNR is unchanged. As the baseline increases up to a certain length, MSE(βˆvF ) increases rather than decreasing. For this reason, according to (46), the ambiguity threshold of baseline is given by the following equation: D = min D

   MSE(βˆvF ) D

(47)

According to the definition of the ambiguity threshold of SNR, MSE(βˆvF ) increases with decreasing SNR when the baseline is unchanged. And it will keep this trend with decreasing SNR. For this reason, we cannot use an equation to determine the ambiguity threshold of SNR. We can plot the curve of MSE(βˆvF ) versus SNR according to (46). MSE(βˆvF ) will remarkably increase when SNR decreases to a certain point. We take this point as the approximation of the ambiguity threshold of SNR. And the approximation is marked as SNR. 5 Complexity Analysis The computational complexity of the proposed algorithm is dominated by (1) formation of the covariance matrix, (2) calculation of its eigendecomposition, (3) the procedure of choose the optimal m, and (4) calculation of the βvC , βvSA , and βvF . According to [11], we can compute that process (1) requires O[(2M)2 L] FLOPS and process (2) requires O[(2M)3 ] FLOPS; process (3)   M−2 requires O (M − M/2 − 1)V (2(V − 1)(2M)2 + (2M)2 (2M − V )) + m=M/2  (4m(2M)2 + 4M(2m)2 ) FLOPS; process (4) requires O[2V 2 (2M −2)+2V 2 (2m)+ ¯ 2V 2 (2M) +V 3 ] FLOPS. Therefore, the computational load of proposed algorithm is mainly O (2M)2 L + (2M)3 + (M − M/2 − 1)V (2(V − 1)(2M)2 + (2M)2

Circuits Syst Signal Process Table 1 Comparison of different DOA estimation algorithms Algorithms

The proposed method

MUSIC-ESPRIT in [32] DS-ESPRIT in [36] U-DR-ESPRIT in [37]

Computational complexity  O (2M)2 L + (2M)3 + (M − M/2 − 1)V [2(V − 1)(2M)2  M−2 +(2M)2 (2M − V )] + m=M/2 (4m(2M)2 + 4M(2m)2 )  ¯ + 2V 2 (2M) + V 3 +2V 2 (2M − 2) +2V 2 (2m)  O (2M)2 L + (2M)3 + 2V 2 (2M) + V 3 +V Ns (2M + 1)(2M − V )}   O (2M)2 L + (2M)3 + 2V 2 (2M − 2) + 2V 2 (2M) + V 3   O (2M)2 2L + (2M)3 + 2V 2 (2M − 2) + 2V 2 (2M) + V 3

 M−2 (4m(2M)2 + 4M(2m)2 ) + 2V 2 (2M − 2) + 2V 2 (2m) ¯ +  m=M/2 2 3 2V (2M) + V . Table 1 shows the complexity comparison of the existing algorithms in [32,36,37] and the proposed method, where Ns is the total search times within the search range, m¯ is the best choice of m. It can be seen that the proposed method has larger computational complexity than that of the other algorithms. The proposed method may not satisfy the requirements of real-time process in the case of more numbers of sources and sensors. The main cause is that the process (3) requires more computation. The m¯ can be approximated as a constant, when DOAs do not significantly change. If DOAs do not significantly change for a while, the process (3) can be done only once. During this period, the computation cost of the proposed method will greatly reduce, thus satisfying the requirements of real-time process. (2M − V )) +

6 Computer Simulation Results In this section, several numerical examples are presented to demonstrate the efficiency of the proposed DOA estimation algorithm. In the following simulations, we assume that each subarray of DSA has 16 sensors, i.e., M = 16, the interelement spacing d = λ/2, and the number of snapshots is L = 100. Figures 5 and 6 show the average root mean square error (RMSE) of DOA estimation versus SNR and the baseline, respectively, compared with the algorithms in [32,36], and [37]. The average RMSE of the DOA estimation is defined as $ % K  V % 1  2 (θˆv (i) − θv ) RMSE =& KV i=1 v=1

where θˆv (i) denotes the DOA estimation of ith trial of vth incident signal, θv denotes the true DOA of vth incident signal, and K is the number of total Monte Carlo trials. There are two uncorrelated equi-powered sources impinging on the DSA with θ1 = 15◦

Circuits Syst Signal Process

Fig. 5 RMSE estimation versus SNR with M = 16, D = 50λ, L = 100, and 5,000 Monte Carlo trials

Fig. 6 RMSEs of DOA estimation versus baseline with M = 16 , L = 100, SNR = 0 dB, and 5,000 Monte Carlo trials

and θ2 = 20◦ . The number of Monte Carlo trials is 5,000. In Fig. 5, we set the baseline D = 50λ. In Fig. 6, we set SNR = 0 dB. It is seen from Fig. 5 that the performance of the proposed method is close to that of the other algorithms at SNR ≥ 0 dB. It is also seen that the DOA estimation performance of the proposed method is much better than that of the other algorithms at low SNRs (e.g., from −10 to −3 dB). It is seen from Fig. 6 that the performance of the proposed method is close to that of the other algorithms at D ≤ 56λ. However, the DOA estimation performance of the proposed method is much better than that of the other algorithms when the baseline is large (e.g., from 92λ to 200λ). Figure 7 shows the MSE of the intermediate estimation versus m. In this simulation, we set the baseline D = 50λ, SNR = 0 dB, and perform 5,000 Monte Carlo trials. There are two uncorrelated equi-powered sources impinging on the DSA from DOAs

Circuits Syst Signal Process

Fig. 7 MSEs of the intermediate estimation versus m with the source at θ1 = 15◦ (β1 = 0.2588) with M = 16, L = 100, D = 50λ, SNR = 0 dB, and 5,000 Monte Carlo trials

(a)

(b)

Fig. 8 The SNR ambiguity threshold of DOA estimation with the source at θ1 = 15◦ (β1 = 0.2588) with M = 16, L = 100, D = 50λ, L = 100, and 20,000 Monte Carlo trials. a DOA = 15◦ (β1 = 0.2588). b DOA = 15◦ , 20◦ (β = 0.2588, 0.342)

θ1 = 15◦ and θ2 = 20◦ ( β1 = 0.2588, β2 = 0.342). Because we use the MSE of the direction cosine to analyze the theoretical MSE, we define the MSE of the direction cosine as: MSEv =

K 1  2 (βˆv (i) − βv ) K i=1

Figure 7 plots the MSE with the source at β1 = 0.2588. It is seen that the MSE approximation obtained by the method presented in Sect. 4.1 is very close to the simulation results. It is also seen that m¯ = 10 is the best choice of m obtained by both the simulation result and the method presented in Sect. 4.2. It implies that the proposed method can obtain the best choice of m. Moreover, it shows that the MSE corresponding to the best choice of m is much better than that obtained by the coarse estimation.

Circuits Syst Signal Process

(a)

(b)

Fig. 9 The baseline ambiguity threshold of DOA estimation with the source at θ1 = 15◦ (β1 = 0.2588) with M = 16, L = 100, SNR = 0 dB, and 20,000 Monte Carlo trials. a DOA = 15◦ (β1 = 0.2588). b DOA = 15◦ , 20◦ (β = 0.2588, 0.342)

Figure 8a, b show the SNR ambiguity threshold. In this simulation, we set D = 50λ. Monte Carlo trials are equal to 20,000. A single source with θ1 = 15◦ (β1 = 0.2588) and two uncorrelated equi-powered sources with θ1 = 15◦ and θ2 = 20◦ ( β1 = 0.2588, β2 = 0.342) are used in the simulation result shown in Fig. 8a, b, respectively. Figure 8b plots the MSE with the source at β1 = 0.2588. It is seen from Fig. 8a, b that the MSE approximation obtained by the method presented in Sect. 4.3 is very close to the simulation result. It is also seen from these figures that the ambiguity threshold of SNR obtained by the simulation result is equal to that of the method presented. These results demonstrate that the method of SNR ambiguity threshold analysis presented in Sect. 4.3 is effective. Figure 9a, b show the baseline ambiguity threshold. In this simulation, we set SNR = 0 dB, and perform 20,000 Monte Carlo trials. The sources localization is same to the last simulation ones. Figure 9b plots the MSE with the source at β1 = 0.2588. It is seen from Fig. 9a, b that the MSE approximation obtained by the method presented in Sect. 4.3 is very close to the simulation result. It is also seen from these figures that the ambiguity threshold of the baseline obtained by the simulation result is equal to that of the method presented. These results demonstrate that the method of the baseline ambiguity threshold analysis presented in Sect. 4.3 is effective. 7 Conclusions In this paper, we have studied the DOA estimation problem for the DSA which is a sparse array consisting of two or more far separated subarrays. The DSA is capable of extending array aperture without adding antennas but leads to the problem of ambiguous estimation. In order to avoid or solve the ambiguities, a novel ESPRITbased method by using subarray partition for DSA DOA estimation has been proposed. The proposed method has better DOA estimation performance than the conventional algorithms do in cases of the low SNR and the large spacing between the subarrays. It is in view of the fact that we use the optimized intermediate estimation instead of the coarse estimation as the reference estimation in the disambiguation procedure. The optimized intermediate estimation is obtained by the rotational invariance between the

Circuits Syst Signal Process

new subarrays with optimal partition of DSA. Moreover, we have given a method of calculating the ambiguity threshold of SNR and baseline. This threshold information is very valuable to the DSA system designer. Acknowledgments The authors would like to sincerely thank all the anonymous reviewers, the Associate Editor, and the Editor-in-Chief for their valuable comments and suggestions, which have greatly improved the quality of this paper. This work is supported by National Natural Science Foundation of China( 61101244).

References 1. Y.I. Abramovich, N.K. Spencer, A.Y. Gorokhov, Resolving manifold ambiguities in direction-of-arrival estimation for nonuniform linear antenna arrays. IEEE Trans. Signal Process. 47, 2629–2643 (1999) 2. F. Athley, Direction-of-arrival estimation using separated subarrays, in The 34th asilomar conference on signals, systems, and computer (Pacific Grove, CA, 2000), pp. 585–589 3. F. Athley, Optimization of element positions for direction finding with sparse arrays, in Proceedings of the 11th IEEE signal processing workshop on statistical signal processing (Singapore, 2001), pp. 516–519 4. F. Athley, Threshold region performance of maximum likelihood direction of arrival estimators. IEEE Trans. Signal Process. 53, 1359–1373 (2005) 5. J.A. Cadzow, Y.S. Kim, D.C. Shiue, General direction-of-arrival estimation: a signal subspace approach. IEEE Trans. Aerosp. Electron. Syst. 25, 31–47 (1989) 6. B. Champagne, Source detection and DOA estimation from two observations of a finite line aperture, in IEEE international conference acoustics, speech, and signal processing, vol 4, pp. 61–64 (1993) 7. G.H. Chen, B.X. Chen, Eigenstructure-based ambiguity resolution algorithm for distributed subarray antennas VHF radar. Electron. Lett. 48, 788–789 (2012) 8. W. Chen, X. Xu, S. Wen, Super-resolution direction finding with far-separated subarrays using virtual array elements. IET Radar Sonar Navig. 5, 824–834 (2011) 9. S. Coutts, K. Cuomo, F. Robey, Distributed coherent aperture measurements for next generation BMD radar, in IEEE sensor array and multichannel signal processing workshop, vol 1, pp. 390–393 (2005) 10. B.K. Feng, C.J. David, Grating lobe suppression for distributed digital subarrays using virtual filling. IEEE Antennas Wirel. Propag. Lett. 12, 1323–1326 (2013) 11. G.H. Golub, C.F. Van Loan, Matrix Computation, 3rd edn. (The John Hopkins University Press, Baltimore, 1996) 12. S. Haykin, An Introduction to Analog and Digital Communications (Wiley, New York, 1986) 13. Z. He, Z. Zhao, Z. Nie, P. Ma, Q.H. Liu, Resolving manifold ambiguities for sparse array using planar substrates. IEEE Trans. Antennas Propag. 60, 2558–2562 (2012) 14. Z. He, Z. Zhao, Z. Nie, P. Tang, J. Wang, Q.H. Liu, Method of solving ambiguity for sparse array via power estimation based on MUSIC algorithm. Signal Process. 92, 542–546 (2012) 15. P. Hyberg, M. Jansson, B. Ottersten, Array interpolation and DOA MSE reduction. IEEE Trans. Signal Process. 53, 4464–4471 (2005) 16. A.N. Lemma, A.J. van der Veen, E.F. Deprettere, Multiresolution ESPRIT algorithm. IEEE Trans. Signal Process. 47, 1722–1726 (1999) 17. B. Li, B. Xu, Y.S. Yuan, Preestimation-based array interpolation approach to coherent source localization using multiple sparse subarrays. IEEE Signal Process. Lett. 16, 81–84 (2009) 18. B. Liao, S.C. Chan, Direction-of-arrival estimation in subarrays-based linear sparse arrays with gain/phase uncertainties. IEEE Trans. Aerosp. Electron. Syst. 49, 2268–2280 (2013) 19. J. Liu, Z. Zhao, Z. He, Z. Nie, Q.H. Liu, Resolving manifold ambiguities for direction-of-arrival estimation of sparse array using semi-circular substrates. IET Microw. Antennas Propag. 7, 1016– 1020 (2013) 20. G. Motti, J.W. Anthony, Array geometry for ambiguity resolution in direction finding. IEEE Trans. Antennas Propag. 44, 889–895 (1996) 21. J.E. Nilsson, H. Warston, Radar with separated subarray antennas, in Proceedings of IEEE radar conference (Australia, 2003), pp. 194–199

Circuits Syst Signal Process 22. B. Ottersten, M. Viberg, T. Kailath, Performance analysis of the total least squares ESPRIT algorithm. IEEE Trans. Signal Process. 39, 1122–1134 (1991) 23. B. Rao, K.V.S. Hari, Performance analysis of ESPRIT and TAM in determining the direction of arrival of plane waves in noise. Trans. Acoust. Speech Signal Process. 37, 1990–1995 (1989) 24. B. Rao, K.V.S. Hari, Performance analysis of root-music. IEEE Trans. Acoust. Speech Signal Process. 37, 1939–1949 (1989) 25. M. Rbsamen, A.B. Gershman, Subspace-based direction-of-arrival estimation for more sources than sensors using planar arrays, in IEEE sensor array and multichannel signal processing workshop (SAM), pp. 21–24 (2010) 26. R.E. Skelton, D.A. Wagie, Minimal root sensitivity in linear systems. J. Guidance 7, 570–574 (1984) 27. P. Stoica, A. Nehorai, Performance comparison of subspace rotation and MUSIC methods for direction estimation. IEEE Trans. Signal Process. 39, 446–453 (1991) 28. H.L. Sun, H.H. Tao, H.R. Chang, Method of resolving ambiguity for sparse array via modified sparse even array based on MUSIC algorithm, in Second asian-pacific conference synthetic aperture radar, pp. 246–249 (2009) 29. D.W. Tufts, H. Ge, R. Kumaresan, Resolving ambiguities in estimation spatial frequencies in sparse linear array, in IEEE international conference acoustics, speech, and signal processing, vol 2, pp. 345–348 (1994) 30. H.L. Van Trees, Detection, Estimation, and Modulation Theory, Part I: Detection, Estimation and Linear Modulation Theory (Wiley, New York, 2002) 31. H.L. Van Trees, Detection, Estimation, and Modulation Theory, Part IV: Optimum Array Processing (Wiley, New York, 2002) 32. V.I. Vasylyshyn, Direction of arrival estimation using ESPRIT with sparse arrays, in European radar conference (Rome, Italy, 2009) pp. 246–249 33. V.I. Vasylyshyn, Unitary ESPRIT-based DOA estimation using sparse linear dual size spatial invariance array, in European radar conference (Paris, France, 2005), pp. 157–160 34. Q. Wang, R.B. Wu, M.D. Xing, Z. Bao, A new algorithm for sparse aperture interpolation. IEEE Geosci. Remote Sens. Lett. 4, 480–484 (2007) 35. M.J. Wilson, R. McHugh, Sparse-periodic hybrid array beamformer. IET Radar Sonar Navig. 1, 116– 123 (2007) 36. K.T. Wong, M.D. Zoltowski, Direction-finding with sparse rectangular dual-size spatial invariance array. IEEE Trans. Aerosp. Electron. Syst. 34, 1320–1336 (1998) 37. G.M. Zheng, B.X. Chen, Unitary dual-resolution ESPRIT for joint DOD and DOA estimation in bistatic MIMO radar. Multidim. Syst. Signal Process. (2013). doi:10.1007/s11045-013-0244-5 38. G.M. Zheng, B. Wu, Y. Ma, B.X. Chen, Direction of arrival estimation with a sparse uniform array of orthogonally oriented and spatially separated dipole-triads. IET Radar Sonar Navig. 8, 885–894 (2014) 39. M.D. Zoltowski, K.T. Wong, Closed-form eigenstructure-based direction finding using arbitrary but identical subarrays on a sparse uniform Cartesian array grid. IEEE Trans. Signal Process. 48, 2205– 2211 (2000)