# Efficient Two-Dimensional Direction Finding via Auxiliary-Variable Manifold Separation Technique for Arbitrary Array Structure

Authors: Guang Hua; Jiu-Dong Wu; Xi-Cheng Zhu; Hou-Xing Zhou; Wei Hong
Journal: Mathematical Problems in Engineering
 (2015)

---

## Abstract

A polynomial rooting direction of arrival (DOA) algorithm for multiple plane waves incident on an arbitrary array structure that combines the multipolynomial resultants and matrix computations is proposed in this paper. Firstly, a new auxiliary-variable manifold separation technique (AV-MST) is used to model the steering vector of arbitrary array structure as the product of a sampling matrix (dependent only on the array structure) and two Vandermonde-structured wavefield coefficient vectors (dependent on the wavefield). Then the propagator operator is calculated and used to form a system of bivariate polynomial equations. Finally, the automatically paired azimuth and elevation estimates are derived by polynomial rooting. The presented algorithm employs the concept of auxiliary-variable manifold separation technique which requires no sector by sector array interpolation and thus does not suffer from any mapping errors. In addition, the new algorithm does not need any eigenvalue decomposition of the covariance matrix and exhausted search over the two-dimensional parameter space. Moreover, the algorithm gives automatically paired estimates, thus avoiding the complex pairing procedure. Therefore, the proposed algorithm shows low computational complexity and high robustness performance. Simulation results are shown to validate the effectiveness of the proposed method.


---

## Body

## 1. Introduction

Direction finding of the propagating plane waves has received much attention in a variety of applications, including radar, sonar, seismology, and mobile communication. Uniform linear array is often exploited in direction of arrival estimation as it has a steering vector matrix with a Vandermonde structure. Some fast direction finding estimators designed for uniform linear array, root-MUSIC, and root-WSF [1–4], for example, can then be applied to give estimates of source bearings relative to the array axis. In [5] a variant of ESPRIT called unitary ESPRIT is proposed for uniform rectangular arrays.A planar array is often employed if the estimates of both source azimuth and elevation are required. The techniques of beamspace transform [6] and array interpolation [7–10] have been developed to map the steering vector of any planar array onto a steering vector of a virtual uniform linear array. The two preprocessing techniques that are mentioned above may introduce mapping errors [11] in the form of bias [3, 12] and excess variance [13]. Hence, the estimates are not statistically optimal.Manifold separation technique (MST) [14] exploits the wavefield modeling formalism [15, 16] to model the steering vector of any array configurations as the product of a sampling matrix and a Vandermonde-structured angular vector. By using MST, conventional root-MUSIC originally designed for ULA can now be applied to arbitrary array configurations. Compared to mapping and interpolation techniques, MST does not require any sector by sector interpolation or transformation, and mapping errors can thus be avoided. Besides, MST can provide smaller fitting errors over the whole coverage area. Fourier-domain root-MUSIC [17] can also extend the root-MUSIC to arbitrary array configurations with improved performance to computational complexity tradeoffs at the cost of increasing the degree of the polynomial. All of the two mentioned algorithms can only provide one-dimensional estimate of the source bearing. As for the two-dimensional (azimuth and elevation) estimation problem, spherical harmonics [15, 18] have to be introduced.Here, a new computationally efficient direction finding algorithm for arbitrary array structure was developed. Firstly, a variant of MST called auxiliary-variable manifold separation technique (AV-MST) is proposed. AV-MST employs a combination of product-to-sum formula and Jacobi-Anger expansion [19] to model the steering vector of arbitrary array structure as the product of a sampling matrix and two Vandermonde-structured wavefield coefficient vectors. This contribution replaces complicated spherical harmonics computations by the Bessel functions of the first kind and allows fast direction finding estimators to be applied in element space. Then propagator operator [20] is introduced and a system of bivariate polynomial equations is formed; this procedure requires no eigenvalue decomposition of the covariance matrix; thus the computational complexity is reduced. Finally, by using resultants [21, 22], the problem of finding the roots of a polynomial system reduces to a matrix computation of generalized eigenvalue problem which involves no iterative procedure and exhausted 2D searching, a closed-form automatically paired two-dimensional angle estimates are eventually derived.The rest of this paper is organized as follows. In Section2, key assumptions are made and the signal model is defined. In Section 3, the novel AV-MST was presented. In Section 4, the proposed root propagator method (root-PM) algorithm for 2D direction finding was derived. In Section 5, simulation results are given to verify the efficiency of the proposed algorithm. In Section 6, conclusions were made in the paper.Throughout the paper, scalar quantities are denoted by regular lowercase letters and lowercase bold type faces are for vectors. Regular uppercase letters are used for matrices. SuperscriptsT, H, *, and † in this paper denote the transpose, conjugate transpose, complex conjugate, and pseudoinverse, respectively. ⊗ symbolizes the Kronecker-product operator. 0m,n and Im denote the m×n zero matrix and m×m identity matrix, respectively.


## 2. Signal Model

Consider thatK noncoherent narrowband planer wave source signals, parameterized by {θ1,ϕ1},{θ2,ϕ2},…,{θK,ϕK}, impinge upon an array of N (N≥K) sensors. Source elevation angles θ∈[0,π/2], and azimuth angles ϕ∈[0,2π]. It was assumed that L snapshots are observed by the array. The N×L array output matrix X can be written as (1)X=As+Nwith (2)A=[a(θ1,ϕ1),a(θ2,ϕ2),a(θ3,ϕ3),…,a(θK,ϕK)],where (3)a(θn,ϕn)=ej2π/λr1sinθncos⁡⁡ϕn-γ1,…,ej2π/λrNsinθncos⁡⁡ϕn-γNT,s=s1,s2,…,sLT,where {θn,ϕn} is the bearing of nth source and A is the N×L matrix of the steering vector a(θn,ϕn). s is L-dimensional vector complex amplitude of the sources and N is N×L noise matrix, and the noise is modeled as a stationary, second-order, ergodic, zero-mean, spatially and temporally white, and circular Gaussian process.


## 3. Auxiliary-Variable Manifold Separation Technique

In this section, a novel variant of MST called auxiliary-variable manifold separation technique (AV-MST) is proposed. By using AV-MST, the steering vector can be modeled as the product of a sampling matrix and two Vandermonde-structured wavefield coefficient vectors. The proposed array modeling technology is based on the wavefield modeling method can be expressed in formula [16]. Here an array formed with omnidirectional sensors was considered. This assumption does not restrict the generality of the discussion. In (2), the kth element of the steering vector corresponds to the kth sensor for arrays with arbitrary structure. Using the product-to-sum formula, it can be written as(4)aθ,ϕk=ej(2π/λ)rksinθcos⁡(ϕ-γk)=ej(π/λ)rk(sin(θ+ϕ-γk)+sin(θ-ϕ+γk))=ej(π/λ)rksin(θ+ϕ-γk)ej(π/λ)rksin(θ-ϕ+γk);here put (5)θ*=θ-ϕ,ϕ*=θ+ϕ;substituting (5) into (4), which gives (6)aθ,ϕk=ej(π/λ)rksin(ϕ*-γk)ej(π/λ)rksin(θ*+γk)by using Jacobi-Anger expansion [19], (6) can be written as the following(7)aθ,ϕk=ej(π/λ)rksin(ϕ*-γk)ej(π/λ)rksin(θ*+γk)=∑m=-∞∞e-jmγkJm(πrkλ)ejmϕ*·∑n=-∞∞ejnγkJn(πrkλ)ejnθ*=∑m=-∞∞Gm(rk,γk)ejmϕ*∑n=-∞∞Gn(rk,γk)ejnθ*,where Gm(rk,γk)=e-jmγkJm(πrk/λ) and Gn(rk,γk)=ejnγkJn(πrk/λ). Note that Gm(rk,γk) and Gn(rk,γk) depend on the array configurations only and they are independent of the wavefield. Jm(·) and Jn(·) represent the Bessel functions of the first kind of orders m and n, respectively. For m=[-M,M] and n=[-M,M], M is the maximum mode number and it is proposed that M≈2κr, where r is the radius of the smallest circle centered at the origin of the array, enclosing all the physical components. Then we can define (8)Gk=G1k⊗G2k,V(ϕ*,θ*)=v(ϕ*)⊗v(θ*),where (9)G1k=G2k=G-Mrk,γk,G1-Mrk,γk,…,mmcGM-1(rk,γk),GM(rk,γk),v(ϕ*)=[e-jMϕ*,e-j(M-1)ϕ*,…,ej(M-1)ϕ*,e-jMϕ*],v(θ*)=[e-jMθ*,e-j(M-1)θ*,…,ej(M-1)θ*,e-jMθ*];substituting (8) into (7), the kth element of the steering vector can be modelled by(10)aϕ,θk=GkVϕ*,θ*+εM;then the steering vector can be written as(11)aϕ,θ=GV(ϕ*,θ*)+ε*(M,M)=G1G2⋮GNV(ϕ*,θ*)+ε*(M,M).The quantities ε(M) and ε*(M,M) represent the modeling error due to truncation. Ideally, the resulting truncation error can be made arbitrary small as the maximum number mode M increases. Note that all the modeling procedure that is listed above needs to be done only once and can be done offline; then the direction finding algorithm can be applied to give the 2D source bearing estimates.


## 4. Element-Space 2D Root-PM

### 4.1. Construction of a System of Two Bivariate Polynomial Equations

On the assumption that the incoming signals are noncoherent, the steering matrixA is of full rank which is common for all the subspace-based techniques. So there exists a unique linear propagator operator P and a matrix Q can then be formed whose columns are orthogonal to the steering vector [20], and the DOAs are the solutions of (12)aϕ*,θ*HQQHa(ϕ*,θ*)=0,where Q∈CN×(N-K); from (11) and (12), one can equivalently construct the following DOA estimator: (13)Vϕ*,θ*HGHQQHGV(ϕ*,θ*)=0.In order to find the unique solutions about {ϕ,θ}, a system of two bivariate polynomial equations can be constructed by decomposing the spanned subspace, then the matrix Q has been divided into two submatrices Q1 and Q2 that satisfy(14)Q1∪Q2=Q,Q1≠Q2.Using the two submatrixes Q1 and Q2, a system of bivariate polynomial equations related to the DOA of the incoming signals is given by (15)Vϕ*,θ*HG1HQQHG1V(ϕ*,θ*)=0,Vϕ*,θ*HG2HQQHG2V(ϕ*,θ*)=0.Define two complex variables z and w such that (16)z=ejϕ*,w=ejθ*.Viewing (15) as polynomials in z with coefficients that are polynomials in w [21], it can be written as (17)fzTE1g(w)=0,fzTE2g(w)=0.Here f(z)=[z4M,…,1]T∈C(4M+1)×1 and g(w)=[w4M,…,1]T∈C(4M+1)×1. E1 and E2 which contain the coefficients of the bivariate polynomial in (15) can be obtained as (18)E1(x1,x2)=∑diag{∑diag{G1HQQHG1,x1},x2},E2(x1,x2)=∑diag{∑diag{G2HQQHG2,x1},x2}.Here, ∑diag{E,x} symbolizes the sum of the block element of the xth block diagonal of matrix E.


### 4.2. Azimuth-Elevation Angles Estimation

In order to find the solutions to the systems of bivariate polynomial equations in (17), [22] provides a computationally efficient technique by using a combination of multipolynomial resultants and generalized eigenvalue decomposition. The multipolynomial resultant R(w) can be obtained using the Sylvester matrix [21] and it is a polynomial in w. R(w) linearizes (15) and reduces it to a linear system of the form [22] (19)R(w)t(z)=0,0,…,0T,where t(z)=[z8M-1,z8M-2,…,1]∈C8M×1. The solutions of (17) correspond to the roots of the determinant of multipolynomial resultant R(w), and R(w) in matrix polynomial form [22] can be expressed as (20)Rw=R0+R1w+R2w2+⋯+R4Mw4M.Given the matrix polynomial form of R(w), the roots of the polynomial corresponding to its determinant are the generalized eigenvalues [23] of the matrixes C1 and C2 [22]; that is, (21)C2v=wC1v,where C1 and C2 are obtained by using the coefficient matrix of R(w). If wi is a generalized eigenvalue of matrixes C1 and C2, the corresponding root zi can be derived as (22)zi=vi(1)vi(2),where vi is the ith element of the generalized eigenvector corresponding to the eigenvalue wi. Note also that the proposed algorithm provides automatically paired roots which are obtained from (21) and (22); it is assumed that they are [w,z]. The K auxiliary roots relating to the virtual DOAs can be selected roughly by choosing the paired ones with ∠z=ϕ*=θ+ϕ>0 first; further the paired roots are selected with magnitude that are simultaneously closest to both unit circles. If [wi,zi] is the ith of the paired K roots, then the true DOAs can be obtained according to (5) by (23)θi=ϕi*+θi*2=∠zi+∠wi2,ϕi=ϕi*-θi*2=∠zi-∠wi2.


### 4.3. Implementation of the Proposed Method

Assuming that the number of signals is known or correctly estimated, the proposed element-space 2D root-PM algorithm for arbitrary structure arrays can be implemented via the following steps.(1)
Using AV-MST to model the steering vector as the product of an array sampling matrixG dependent on the array structure and wavefield coefficient vector V(ϕ*,θ*).(2)
Form the received data covariance matrix to compute the propagator operator and a matrixQ [20] whose columns are orthogonal to the steering vector a(ϕ*,θ*) can then be formed. Using (11)–(14), a system of two bivariate polynomial equations can then be constructed in the form of (15).(3)
By using Sylvester matrix [21], multipolynomial resultants R(w) can be computed in the form of matrix polynomial to linearize (15) to the form of  (19).(4)
Then the polynomial rooting problem can be reduced to the generalized eigenvalue problem.(5)
After the paired roots of[w,z] are obtained, by using (23), the true DOAs can be estimated by the phase angles of these paired roots.


## 4.1. Construction of a System of Two Bivariate Polynomial Equations

On the assumption that the incoming signals are noncoherent, the steering matrixA is of full rank which is common for all the subspace-based techniques. So there exists a unique linear propagator operator P and a matrix Q can then be formed whose columns are orthogonal to the steering vector [20], and the DOAs are the solutions of (12)aϕ*,θ*HQQHa(ϕ*,θ*)=0,where Q∈CN×(N-K); from (11) and (12), one can equivalently construct the following DOA estimator: (13)Vϕ*,θ*HGHQQHGV(ϕ*,θ*)=0.In order to find the unique solutions about {ϕ,θ}, a system of two bivariate polynomial equations can be constructed by decomposing the spanned subspace, then the matrix Q has been divided into two submatrices Q1 and Q2 that satisfy(14)Q1∪Q2=Q,Q1≠Q2.Using the two submatrixes Q1 and Q2, a system of bivariate polynomial equations related to the DOA of the incoming signals is given by (15)Vϕ*,θ*HG1HQQHG1V(ϕ*,θ*)=0,Vϕ*,θ*HG2HQQHG2V(ϕ*,θ*)=0.Define two complex variables z and w such that (16)z=ejϕ*,w=ejθ*.Viewing (15) as polynomials in z with coefficients that are polynomials in w [21], it can be written as (17)fzTE1g(w)=0,fzTE2g(w)=0.Here f(z)=[z4M,…,1]T∈C(4M+1)×1 and g(w)=[w4M,…,1]T∈C(4M+1)×1. E1 and E2 which contain the coefficients of the bivariate polynomial in (15) can be obtained as (18)E1(x1,x2)=∑diag{∑diag{G1HQQHG1,x1},x2},E2(x1,x2)=∑diag{∑diag{G2HQQHG2,x1},x2}.Here, ∑diag{E,x} symbolizes the sum of the block element of the xth block diagonal of matrix E.


## 4.2. Azimuth-Elevation Angles Estimation

In order to find the solutions to the systems of bivariate polynomial equations in (17), [22] provides a computationally efficient technique by using a combination of multipolynomial resultants and generalized eigenvalue decomposition. The multipolynomial resultant R(w) can be obtained using the Sylvester matrix [21] and it is a polynomial in w. R(w) linearizes (15) and reduces it to a linear system of the form [22] (19)R(w)t(z)=0,0,…,0T,where t(z)=[z8M-1,z8M-2,…,1]∈C8M×1. The solutions of (17) correspond to the roots of the determinant of multipolynomial resultant R(w), and R(w) in matrix polynomial form [22] can be expressed as (20)Rw=R0+R1w+R2w2+⋯+R4Mw4M.Given the matrix polynomial form of R(w), the roots of the polynomial corresponding to its determinant are the generalized eigenvalues [23] of the matrixes C1 and C2 [22]; that is, (21)C2v=wC1v,where C1 and C2 are obtained by using the coefficient matrix of R(w). If wi is a generalized eigenvalue of matrixes C1 and C2, the corresponding root zi can be derived as (22)zi=vi(1)vi(2),where vi is the ith element of the generalized eigenvector corresponding to the eigenvalue wi. Note also that the proposed algorithm provides automatically paired roots which are obtained from (21) and (22); it is assumed that they are [w,z]. The K auxiliary roots relating to the virtual DOAs can be selected roughly by choosing the paired ones with ∠z=ϕ*=θ+ϕ>0 first; further the paired roots are selected with magnitude that are simultaneously closest to both unit circles. If [wi,zi] is the ith of the paired K roots, then the true DOAs can be obtained according to (5) by (23)θi=ϕi*+θi*2=∠zi+∠wi2,ϕi=ϕi*-θi*2=∠zi-∠wi2.


## 4.3. Implementation of the Proposed Method

Assuming that the number of signals is known or correctly estimated, the proposed element-space 2D root-PM algorithm for arbitrary structure arrays can be implemented via the following steps.(1)
Using AV-MST to model the steering vector as the product of an array sampling matrixG dependent on the array structure and wavefield coefficient vector V(ϕ*,θ*).(2)
Form the received data covariance matrix to compute the propagator operator and a matrixQ [20] whose columns are orthogonal to the steering vector a(ϕ*,θ*) can then be formed. Using (11)–(14), a system of two bivariate polynomial equations can then be constructed in the form of (15).(3)
By using Sylvester matrix [21], multipolynomial resultants R(w) can be computed in the form of matrix polynomial to linearize (15) to the form of  (19).(4)
Then the polynomial rooting problem can be reduced to the generalized eigenvalue problem.(5)
After the paired roots of[w,z] are obtained, by using (23), the true DOAs can be estimated by the phase angles of these paired roots.


## 5. Simulations

This section will verify the effectiveness of the proposed auxiliary-variable manifold separation technique and 2D root-PM algorithm for two-dimensional direction finding. In particular, a uniform circular array of radiusr=0.5λ consisting of N=19 elements is employed, with the maximum mode number being M=6, and the reference point is at the center of the circle. The number of snapshots per trial was L=512 and 1000 independent trials in total were simulated.Figures1 and 2 show the corresponding histograms of the joint elevation and azimuth DOA estimation for two (K=2) uncorrelated equipowered signals with arrival angles θ1,ϕ1=(72.7°,90°) and (θ2,ϕ2)=(50.4°,78°). Figure 1 was simulated at signal-to-noise ratio (SNR) = 5 dB and Figure 2 was simulated at SNR = 20 dB. In Figure 1, it was observed that the histograms spread around the true DOAs, while Figure 2 gives very sharp histograms and two dominant peaks can be seen compared to curves that in Figure 1. It can be concluded that the proposed algorithm gives more accurate DOA estimates when the SNR increases.Figure 1
Histograms of DOA estimates for two signals with parametersθ1=72.7°, ϕ1=90.0°, θ2=50.4°, and ϕ2=78.0°: elevation (a); azimuth (b), SNR = 5 dB.
(a)
(b)Figure 2
Histograms of DOA estimates for two signals with parametersθ1=72.7°, ϕ1=90.0°, θ2=50.4°, and ϕ2=78.0°: elevation (a); azimuth (b), SNR = 20 dB.
(a)
(b)Figures3 and 4 show the root mean squared error (RMSE) of the joint elevation and azimuth DOA estimates versus the SNR in dB for the two sources. It was observed from Figures 3 and 4 that the RMSE of the joint elevation and azimuth DOA estimates decreases as the SNR increasing. The two figures demonstrated that the proposed algorithm provides asymptotically similar performance to the Cramer Rao Bound.Figure 3
RMSEs of DOA estimates versus SNRs for the first signal with parametersθ1=72.7°, ϕ1=90.0°: elevation (a); azimuth (b).
(a)
(b)Figure 4
RMSEs of DOA estimates versus SNRs for the second signal with parametersθ2=50.4°, ϕ2=78.0°: elevation (a); azimuth (b).
(a)
(b)


## 6. Conclusions

In this paper, the problem of azimuth and elevation direction finding for multiple plane waves employing an arbitrary array structure was considered. The presented algorithm requires no sector by sector array interpolation and thus does not suffer from any mapping errors. In addition, the new algorithm does not need any eigenvalue decomposition of the covariance matrix and exhausted search over the two-dimensional parameter space. Moreover, the algorithm gives automatically paired estimates, thus avoiding the complex pairing procedure. Above all, the proposed algorithm shows low computational complexity and performs with high robustness. The simulation results show the accuracy of the proposed algorithm.



---

*Source: 740513-2015-03-03.xml*