Direction-of-Arrival Estimation for Arbitrary Array : Combining Spatial Annihilating and Manifold Separation

In this paper, we address the problemof directionof-arrival (DOA) estimation with the arbitrary array. The manifold separation technique (MST) is employed to transform the arbitrary array into a virtual array with Vandermonde manifold on which the spatial annihilating filter reconstruction method can be applied. When building the optimization problem for annihilating filter reconstruction, we propose the general solution modeling which can reduce the truncation error inMST to a negligible level. Finally, the spatial annihilating filter is reconstructed under the structural total least square (STLS) framework with the multiple measurement vectors structural total least norm (MMV-STLN) approach and the DOAs are estimated from the filter coefficients. Numerical simulations have verified the new proposed method adapts well to the low signal-to-noise ratio (SNR), limited snapshots and closely-spaced sources scenarios and can handle the coherent signals.


Introduction
DOA estimation is regarded as one of the key techniques in various areas, such as wireless communication, radar and sonar [1].For the array geometry, the uniform linear array (ULA) is preferred due to the Vandermonde structure of its array manifold.Owing to this specific structure, it is allowed to adopt some fast DOA estimation algorithms , such as root-MUSIC [2] and ESPRIT [3] which belong to the subspace type method, root-WSF (also called MODE) [4], [5] and IQML [6] which are the fast implementations of the maximum likelihood (ML) method.However, we may face other array geometries without Vandermonde structure for certain applications and even the Vandermonde structure of ULA can be destroyed by the array imperfections [7].
Several methods were proposed to map the physical planar array to the virtual array with Vandermonde structure.The first is the beamspace transform [8] which is based on the phase mode excitation principle [9].However, there exists residual error in this model due to the beam aliasing and this residual error can only be reduced by increasing the number of sensors.Besides, this method is only applicable to the uniform circular array (UCA).The second is the array interpolation method [10] which maps the physical array to a virtual array.Unfortunately, for this method, the mapping error will lead to a biased DOA estimation [11].Furthermore, when the field of interest (FOV) is large, sector-bysector processing is needed, which is known to be sensitive to out-of-sector sources [12].
Another approach for array mapping is the wavefield modeling principle which is reported in [13][14][15].This approach has been termed the MST in [16].It takes advantage of the spatial periodicity of the sensor response and approximates each component of the steering vector by its truncated Fourier series.Then the real manifold can be approximated by the product of the array sampling matrix composed of the Fourier series coefficients and a Vandermonde coefficient vector composed of Fourier bases.In [14], a MST based root-MUSIC is proposed and the Fourier series truncation error therein can be reduced to a negligible level by increasing the order of Fourier basis.In [16], the MST based root-MUSIC was further studied to cope with array imperfections.In [17], the MST was extended to convert the MUSIC null spectrum of the arbitrary array to its Fourier domain and proposed a new search-free DOA estimation algorithm.However, the above methods all belong to the subspace type methods.As we know, the subspace type methods perform poorly in cases like low SNR and limited snapshots and they have to be combined with the spatial smoothing [18] to handle the coherent sources, which reduces the effective array aperture.The ML method does not have these shortcomings.For the arbitrary array, with MST, the ML method also has efficient implementations.In [19], it is suggested MST can be combined with root-WSF.Nevertheless, that method requires the pseudoinverse of the array sampling matrix therein should be a left inverse, which means the array sampling matrix cannot be a fat matrix (more columns than rows) and will limit the order of the Fourier basis, leading to a non-negligible model truncation error.Spatial annihilating filter reconstruction is another DOA estimation method originally designed for ULA and it is equivalent to the deterministic maximum likelihood (DML) estimator of DOA [20], which ensures that it performs well both for coherent sources and under low SNR and limited snapshots scenarios [21].In this paper, we extend the spatial annihilating filter reconstruction method to the arbitrary array.The arbitrary array is transformed into a virtual array with Vandermonde manifold by the MST.While building the optimization problem on the virtual array for the annihilating filter reconstruction, the general solution modeling is proposed to reduce the truncation error in MST to a negligible level.Then the optimization problem is solved under the STLS framework with the MMV-STLN approach and the DOAs are estimated from the reconstructed filter coefficients.Compared with the previous MST based subspace type methods, the new method adapts well to the low SNR,limited snapshots and closely-spaced sources scenarios and can handle the coherent signals.Compared with the previous MST based ML methods, the model truncation error in the new method can be neglected.
Notations used in the paper are introduced as follows.(•) T , (•) H and (•) + are denoted as the transpose, conjugate transpose and pseudo-inverse operator, respectively.• 2 and • F denote the 2 norm and the Frobenius norm, respectively.( * ) and ⊗ are the convolution, and Kronecker product operator respectively.vec(•) is the operator that builds a column vector by stacking the column vectors of a matrix below one another and I N is the N × N identity matrix.a j denotes the jth element of a.A i, j denotes the element at ith row and jth column of A.

Signal Model and Manifold Separation Technique
Consider an M-element arbitrary array exposed to K far-field narrowband sources.
The sensors are located at the x-y plane of the Cartesian coordinate system and the Cartesian coordinate vector of the mth sensor is [r m cos(φ m ), r m sin(φ m )] where [r m , φ m ] is its polar coordinate vector.The sources come from the x-y plane with the DOAs ϕ = [ϕ 1 , ϕ 2 , . . ., ϕ K ] which are within [−π, π].Set the coordinate origin as the phase reference point and the array output vector at snapshot n can be written as where A(ϕ) = [a(ϕ 1 ), a(ϕ 2 ), . . ., a(ϕ K )] is the manifold matrix, a(ϕ) ∈ C M×1 is the steering vector, ] T is the signal vector, and ε[n] ∈ C M×1 is the Gaussian noise vector which is temporally and spatially white with variance σ 2 .
Since the array geometry is arbitrary, there is no Vandermonde structure in the steering vector a(ϕ).However, noticing the spatial periodicity of 2π of each element in a(ϕ), we can express its mth element a m (ϕ) using the Fourier series expansion as and the Fourier series coefficient the g m,q can be written as If the array is imperfect (i.e., steering vector has no analytical expression), g m,q can be approximated by the discrete Fourier transform (DFT) of spatial uniform samples of a m (ϕ) [16] based on the spectrum concentration of a m (ϕ) according to the effecitve aperture distribution function (EADF) concept [22].
If the array is perfect, we can write a m (ϕ) as where β = 2π/λ and λ is the signal wavelength.Substituting (4) into (3) and using the Jacobi-Anger expansion, we can write g m,q as g m,q = 1 2π where J q (•) is the qth order Bessel function of the first kind.
As in (2), we need infinite order of Fourier basis to represent a m (ϕ), which is not practical.So, we truncate the Fourier basis to order of Q − 1 and substitute (5) into (2).Then we have j q J q (βr m )e −jqφ m e jqϕ + ν m (ϕ) (6) where ν m (ϕ) is the model truncation error.Since the Bessel function J q (βr m ) decays super-exponentially as q → ∞ beyond |q| = βr m [13], the model truncation error ν m (ϕ) can be reduced to a negligible level by increasing Q to a large enough value.Substituting ( 6) into (1), we have which is the MST model for the array output.
and d(ϕ) is called the coefficient vector with d q (ϕ) = exp(jqϕ).So, D(ϕ) has the Vandermonde structure and x[n] = D(ϕ)s[n] can be regarded as the output of the virtual array whose manifold is D(ϕ).The approximate relation in (7) stems from the model truncation error and we can substitute the approximately equal sign with the equal sign if a large enough number is assigned to Q.

Spatial Annihilating for MST based DOA Estimation
To handle the coherent sources and achieve better performance under low SNR and limited snapshots scenarios, the ML DOA estimation methods are preferred.For the MST model in (7), root-WSF had been employed to take advantage of the Vandermonde structure for fast implementations [19].However, this method needs M ≥ (2Q − 1) to ensure G + G = I, i.e., G + should be a left inverse.This limits the order of Fourier basis, leading to a non-negligible model truncation error in (7).

Spatial Annihilating for MST Model
Since x[n] in ( 7) can be seen as the output of the virtual array with Vandermonde manifold, it can be annihilated by a spatial annihilating filter and the DOAs can be obtained from the filter coefficients [21].Consider filter whose coefficient vector is denoted by h Then, after applying the filter h to x[n], we get the filter output and write its pth element as where According to (8), ξ(ϕ k ) = 0.So the output of the filter is zero for the input x[n] and we call the filter h the spatial annihilating filter as elements of x[n] are spatial samples .The convolution in ( 9) can be expressed in a matrix form which is where L(•) is a matrix operator and L(x[n]) is a Toeplitz matrix written as So, we can say a K + 1 elementary spatial filter whose coefficients constitute a polynomial with roots being exp(jϕ k ), k = 1, 2, . . ., K can annihilate the virtual array output, leading to (10).And according to the uniqueness of the roots, if we find the filter h satisfies (10), then DOAs are able to be derived from the filter coefficients [23].

General Solution Modeling
But now, the question is we only have access to the arbitrary array output y[n] instead of the virtual array output x[n].From (7), we can extract x[n] as As we know, (12) holds only if G + G = I 2Q−1 .This condition makes G a tall matrix, eg., 2Q − 1 ≤ M, which leads to a non-negligible model truncation error and deteriorates the final DOA estimation performance.
If larger number is assigned to Q , then G + G I and x[n] calculated by (12) will generally not satisfy the annihilating relation in (10).For larger Q which makes 2Q−1 > M, (7) becomes an underdetermined system and we can extract x[n] as which is the general solution form of the underdetermined system.w[n] is the certain unknown vector and selected to make x[n] calculated by ( 13) satisfy (10).Here we have substituted the approximately equal sign with the equal sign since large enough number can be assigned to Q.
where E is the MMV representation of ε[n] and W is the certain unknown matrix.Similarly, the MMV representation of the annihilating relation in (10) can be achieved as L(X)h = 0 (15) where L(•) is also a matrix operator and The next step is to reconstruct h.Substituting ( 14) into (15), we have In the least square (LS) or total least square (TLS) sense, we need to minimize L(G + E) 2 F subjected to (16).However, a better way for more accurate solution is to minimize L(E) 2 F since it only contains the noise component and L(E) can be obtained by based on GG + = I 2Q−1 since we have 2Q − 1 > M. Furthermore, as the second component in ( 16) also contains a certain unknown matrix W, we treat it and the noise component as a whole and denote Then we have L Together with (17) we will find a interesting result that is which is also based on GG + = I 2Q−1 .Therefore, to reconstruct h, the original optimization problem is turned to min Since all these are based on the general solution of underdetermined system in (13), we call the above method the general solution modeling.

STLS and DOA Estimation
As we know, L(GΣ) in (20) has special structure (a vertical stack of Toeplitz matrices), so the new optimization problem is better to be solved under the STLS framework [24].Then the optimization problem can be transformed into min where ω = h 0 and h 0 is the initialization of h which will be introduced later.ω H h = 1 is to ensure h to be unique.
To solve (21), we adopt the MMV-STLN method [21] which is a variation of the classical STLN approach [25][26][27].Due to the fact that L(Σ) is also a vertical stack of Toeplitz matrices, we have where 1) , and For the MMV-STLN method, the iterative calculations are required for solving the unknown ν and h.In the (i + 1)th iteration the first constraint in (24) is linearized around the last step solution which are represented by ν (i) and h (i) .Let ∆ν represent a small change in ν (i) , ∆h represent a small change in h (i) and r(ν (i) , where is the Jacobian of r(ν, h) with respect to ν T , h T T and ν (i) = vec(Σ (i) ).So, in the (i + 1)th iteration, the unknowns become the ∆ν and ∆h and MMV-STLN needs to solve a stand linear equality constrained LS problem which is min ∆ν,∆h where Λ ∈ C (M N +K+1)×((2Q−1)N +K+1) is a block diagonal matrix with the first block being I N ⊗ G and the second block being a (K + 1)-by-(K + 1) zero matrix .After ( 26) is solved, we obtain the estimated ∆ν and ∆h.Then update rule for ν (i+1) and h (i+1) is The iterative calculation is terminated when the iteration converges or the maximum number of iterations I is reached.The iteration is considered to be converged when where ζ is the iteration termination threshold which is a small number.
After the convergence of the iterative claculations, the spatial annihilating filter coefficients are reconstructed as ĥ = [ ĥ1 , ĥ2 , . . ., ĥK+1 ] T .Then, according to (8), the final estimated DOAs φ = [ φ1 , φ2 , . . ., φK ] T can be derived as For the iterative calculations in ( 27), an initialization is needed, which means we need to assign values to ν (0) and h (0) .Due to the non-convexity of (24), multiple local minima exist in the cost function.So, we need to choose an initialization close to the optimal solution.From (8), we can find that h can be estimated by performing the inverse one-sided Ztransform on a polynomial constituted from the pre-estimated DOAs.The pre-estimated DOAs can be coarsely obtained via the simple conventional beamforming (CBF) (or the Capon Beamformer [28]).Assuming φk , k = 1, 2, . . ., K is the preestimated DOAs and K is the number of peaks of the preestimated spatial spectrum, then h (0) can be achieved as The true number of sources K is assumed to be known or estimated by classical source number estimator such as the Akaike information criterion (AIC) method [29] or the minimum description length (MDL) method [30].The closelyspaced sources may make K < K, resulting in K − K zeros in h (0) according to (30).In the simulation part, we will see this has no influence on the resolution of the method.For the initialization of ν, according to (18), we have where φ is the extended coarse DOA estimation obtained by CBF or Capon Beamformer and consists of { φk − θ beam /4, φk , φk + θ beam /4} K k=1 [31].Here, θ beam is the beamwidth.
At last, the steps of the algorithm are concluded in Algorithm 1.

Remarks
According to [21], the spatial annihilating method is equivalent to the DML estimator for DOA estimation.And when spatial annihilating is combined with MST, we adopt the general solution modeling by which the model truncation error in MST can be reduced to a negligible level.So, the proposed DOA estimation can adapt well to the coherent sources and low SNR and limited snapshots scenarios.Even though we have (2Q −1− K)N equations in the constraint of (20), the matrix whose Frobenius norm is to be miniminzed contains (M − K)N rows like the ULA case in [21].So the maximum resolvable sources of the proposed is K = (M − K)N, and then K = M N/(N + 1) [21].So, if the number of snapshots is greater than or equal to M − 1, the maximum number of resolvable sources will be M − 1.
The STLN approach used to solve ( 21) is essentially a Gauss-Newton method because of the involved linearization operation [25].In general, the Gauss-Newton method will converge to the closet local minimum on condition that the local minimum of the cost function in (26) is small [32], [33].This condition holds for many applications.Furthermore, in the simulation part, the proposed method converges under all simulation settings.
For the computational complexity of the proposed method, the most time consuming part is solving the stand linear equality constrained LS problem in (26) involved in each iteration.To solve the stand linear equality constrained LS problem, Lagrangian multipliers can be used and the com- Even though the calculation can be accelerated by the robust and efficient constrained linear least-squares solvers, such as cgglse function in the LAPACK library and lsqlin function in MATLAB, the computational burden will be very heavy for large number of snapshots N. Inspired by [34], when N > K, we can replace the array output Y with its largest K singular vectors weighted by the corresponding singular values denoted by Y s ∈ C M×K to reduce the computational complexity.Here K represents the number of incoherent sources and when sources are all coherent, we assign K = 1.Assume the singular value decomposition of Y to be Y = UΓV H and then Y s = YVF (32) where F = [I K , 0] T .Then the MMV version of (7) becomes

. , s[N]],
and E s = EVF.As D has Vandermonde structure, the anni- is still satisfied like in (15).So, in Algorithm 1, we only need to replace Y with Y s which is calculated from (32).Then the number of snapshots becomes K and the computational complexity for solving the stand linear equality constrained

Simulation Results
In this section, some simulations are performed to exhibit the performance of the proposed method.In all simulations, the results are averaged upon 500 Monte Carlo experiments.The algorithms we choose to use as comparisons are MST-root-WSF [19], MUSIC [35], MST-root-MUSIC [14], FD-root-MUSIC [17] and BF-interpolation [10] which are all capable for the arbitrary array.The Cramer-Rao Lower Bound (CRLB) [36] is also employed as the performance reference.For the proposed method, the iteration termination threshold and the maximum number of iteration are set as ζ = 10 −8 , I = 100, respectively, and the stand linear equality constrained LS problem in Algorithm 1 is solved by using the lsqlin function in MATLAB.For the array geometry, we choose one large non-uniform array and one small non-uniform array.They both have 7 elements and their layouts are shown in Fig. 1 where the locations are normalized by the wavelength.
In the first simulation, we compare the proposed method with the MST-root-WSF to demonstrate the impact of the model truncation error on DOA estimation.To exclude the influence of other factors, lenient signal conditions are considered.We assume two equal power uncorrelated sources from ϕ = [70 • , 120 • ] impinge on the array and N = 100 snapshots are collected under SNR = 20 dB.For MST-root-WSF, the parameter Q is assigned its maximum allowable value which is Q = 4.For the proposed method, we can set Q = 20 owing to adopting the general solution modeling method.The simulation is performed on the large and small array.The results of the large array is shown in Fig. 2. It can be seen the proposed method works properly and MSTroot-WSF method converges to the wrong directions due to the large model truncation error.The results of the small array is shown in Fig. 3.For small array, the truncation error component J Q (βr m ) becomes small, so the results of MST-root-WSF method get better.However, it is still biased and the variance of DOA estimation is larger than the proposed method.
In the second simulation, we test all the DOA estimation algorithms with the uncorrelated and coherent signals.The simulation settings are the same with the previous except for adding the coherent signals with the coherent coefficient being exp(j5/8π) and MST-root-MUSIC and FD-root-MUSIC methods are also assigned Q = 20.For BF-interpolation, the interpolation sector is within [70  In the third simulation, we want to find the appropriate value of Q for the algorithms.For the signal condition, we turn back to uncorrelated signals and other simulation settings are the same as previous.The FD-root-MUSIC, MST-root-MUSIC, the proposed method and the CRLB are included.The simulation is performed on the large array and Q is varied from −4 to 25.The estimation performance is measured by the root-mean-square error (RMSE) of the DOA estimation which is defined by where P is the number of Monte Carlo simulations.The results are plotted in Fig. 4. It can be found with finite Q these algorithms can arrive the CRLB.The MST-root-MUSIC and the proposed methods need nearly the same Q to obtain the same performance since the array sampling matrices G in their models have the same analytical expression.The FDroot-MUSIC needs larger Q since its array sampling matrix needs to be estimated from the MUSIC spectrum.
In the fourth simulation, we test the algorithm performances under different SNRs from −6 dB to 20 dB.And from now on, we consider the harsh signal conditions.The DOAs of the sources are ϕ = [70 • , 83.7 • ] and N = 50 snapshots are collected.For FD-root-MUSIC, MST-root-MUSIC and the proposed method, we can set Q = 20 and MST-root-WSF is not included due to its limited performance.For BF-interpolation, the interpolation sector is within [70 • , 100 • ] with 15 grids and the number of the virtual ULA elements is 6.The simulation is performed on the large array and the results are exhibited in Fig. 5.It can be found that MUSIC has the largest SNR threshold which is 7 dB and both of threshold of FD-root-MUSIC and MST-root-MUSIC are 1 dB.The BF-interpolation cannot reach the CRLB even under high SNR due to its biasedness.The proposed method has the lowest SNR threshold which is −4 dB.
In the fifth simulation, we vary the number of snapshots from 5 to 100.The simulation settings are the same with the previous except for fixing the SNR = 5 dB.Results are exhibited in Fig. 6.The proposed only needs 10 snapshots to resolve the sources and its rmse is very close to the CRLB.Even the angle separation is as small as 5 • , the proposed method is able to resolve them.The angle separation threshold of MST-root-MUSIC is slightly smaller than that of FD-root-MUSIC and BF-interpolation has lower resolution ability than MUSIC when the angle separation is larger than 15 • .

Conclusion
In this paper, the problem of DOA estimation with the arbitrary array is addressed.The MST is employed to transform the arbitrary array into a virtual array with Vandermonde manifold on which the spatial annihilating filter reconstruction method can be applied.When building the optimization problem for annihilating filter reconstruction, the general solution modeling is proposed to reduce the truncation error in MST to a negligible level.Finally, the spatial annihilating filter is reconstructed under the STLS framework with the MMV-STLN approach and the DOAs are estimated from the filter coefficients.Compared with the previous MST based subspace type methods, the new method adapts well to the low SNR, limited snapshots and closely-spaced sources scenarios and can handle the coherent signals.Compared with the previous MST based ML methods, the model truncation error can be neglected.

Fig. 7 .
Fig. 7. RMSE of DOA estimation versus different angle separations.The snapshot threshold of FD-root-MUSIC and MST-root-MUSIC are 25 and 20 respectively.MUSIC has the largest snapshot threshold.For BF-interpolation, the mapping error still dominates the DOA estimation error under large number of snapshots.In the last simulation, we vary the angle separation.The DOAs are ϕ = [70 • , 70 • + ∆ϕ] and ∆ϕ is varied form 4 • to 26 • .The other simulation settings are the same with the previous except for fixing N = 50 and SNR = 5 dB.Results are exhibited in Fig 7.Even the angle separation is as small as 5 • , the proposed method is able to resolve them.The angle separation threshold of MST-root-MUSIC is slightly smaller than that of FD-root-MUSIC and BF-interpolation has lower resolution ability than MUSIC when the angle separation is larger than 15 • .
• , 120 • ] with 25 grids and the number of the virtual ULA elements is 6.The simulation is performed on the small array.The estimation performance is measured by the successful resolution probability and if both| φ1 − ϕ 1 | and | φ2 − ϕ 2 | are smaller than |ϕ 1 − ϕ 2 |/2,where (ˆ) means the estimated value, we claim a successful resolution.The results are shown in Tab. 1.We can find that with uncorrelated signals all algorithms works normally.However, under coherent scenario, only MST-root-WSF and the proposed method achieve 100% successful resolution probability.FD-root-MUSIC gets the same probability as MUSIC because it is based on the MUSIC spectrum.