sensors Article

A Fast Algorithm for 2D DOA Estimation Using an Omnidirectional Sensor Array Weike Nie 1 , Kaijie Xu 2 , Dazheng Feng 3 , Chase Qishi Wu 1,4 , Aiqin Hou 1, * and Xiaoyan Yin 1 1 2 3 4

*

School of Information Science and Technology, Northwest University, Xi’an 710127, China; [email protected] (W.N.); [email protected] (C.Q.W.); [email protected] (X.Y.) School of Electro-Mechanical Engineering, Xidian University, Xi’an 710071, China; [email protected] National Laboratory of Radar Signal Processing, Xidian University, Xi’an 710071, China; [email protected] Department of Computer Science, New Jersey Institute of Technology, Newark, NJ 07102, USA Correspondence: [email protected]; Tel.: +86-29-8830-8119

Academic Editor: Vittorio M.N. Passaro Received: 15 December 2016; Accepted: 28 February 2017; Published: 4 March 2017

Abstract: The traditional 2D MUSIC algorithm fixes the azimuth or the elevation, and searches for the other without considering the directions of sources. A spectrum peak diffusion effect phenomenon is observed and may be utilized to detect the approximate directions of sources. Accordingly, a fast 2D MUSIC algorithm, which performs azimuth and elevation simultaneous searches (henceforth referred to as AESS) based on only three rounds of search is proposed. Firstly, AESS searches along a circle to detect the approximate source directions. Then, a subsequent search is launched along several straight lines based on these approximate directions. Finally, the 2D Direction of Arrival (DOA) of each source is derived by searching on several small concentric circles. Unlike the 2D MUSIC algorithm, AESS does not fix any azimuth and elevation parameters. Instead, the adjacent point of each search possesses different azimuth and elevation, i.e., azimuth and elevation are simultaneously searched to ensure that the search path is minimized, and hence the total spectral search over the angular field of view is avoided. Simulation results demonstrate the performance characters of the proposed AESS over some existing algorithms. Keywords: sensor array; fast algorithm; spectrum peak diffusion effect; Direction of Arrival estimation; 2D MUSIC

1. Introduction Direction of Arrival (DOA) is a critical parameter in many sensor-based systems using radar, sonar and wireless communications. As a well-known high-resolution DOA estimation algorithm, the MUltiple SIgnal Classification (MUSIC) algorithm [1] has attracted considerable attention in the areas concerned with estimating parameters from observations of array output [2]. However, this algorithm suffers from prohibitively high computational cost for its search-based strategy, especially for two dimensional (2D) DOA estimation. The work by Huang [3] indicated that it is generally difficult to control a fixed stepsize in the adaptation process due to the inherent conflict between fast convergence and estimated accuracy. A large stepsize may simplify calculation at the expense of a worse Root Mean Square Error (RMSE); meanwhile, a small stepsize may require more calculation but with an improved RMSE. Variable stepsize algorithms [4] often include a coarse-grained search and a fine-grained search with a large and small stepsize respectively, which can effectively balance the convergence rate and error. By combining MUSIC results of two decomposed linear arrays, Zhou [5] accurately estimated the DOA with

Sensors 2017, 17, 515; doi:10.3390/s17030515

www.mdpi.com/journal/sensors

Sensors 2017, 17, 515

2 of 14

phase ambiguity removed. In its search procedures of two arrays, an adaptive scheme including coarse-grained and fine-grained search is used to reduce the complexity. 1D processing approaches can also be applied for solving the 2D DOA problem to reduce computational complexity. Wang [6] successfully used 1D MUSIC with a novel tree structure to estimate the azimuth and elevation dependently. By constructing a novel relation between azimuth and elevation, Karthikeyan [7] proposed an elegant 1D search technique to compute the 2D DOA estimation using the MUSIC algorithm in the case of a single source. A robust 2D DOA estimation based on two 1D angles estimations using L-shaped array is proposed by Kikuchi [8], and its pair matching of azimuth and elevation is accurate in low SNR and multi-source scenarios. Search-free [9] methods are an important approach for reducing complexity. Han [10] reconstructed a special antenna array based on the Toeplitz matrix whose rank is only related to the DOA of the signals, hence the subspace can be extracted without being affected by the coherence between the signals. With the estimated subspace, an improved ESPRIT method is provided to implement the DOA estimation without spectrum peak search for coherent signal DOA estimation. Chen [11] extended the algorithm in [10] to a 2D situation and presented another search-free algorithm, where they obtained a block Hankel matrix from correlations and found that the rank of the block Hankel matrix only depends on the number of impinging waves. Two matrix pencils including DOA parameters are derived from the signal subspace of the Hankel matrix, and a pairing free method [12] is used to achieve the 2D parameters estimation. As a well-known search-free fast algorithm, Root-MUSIC [13,14] algorithm converted spectrum search into finding roots of a complex polynomial formed from noise subspace. Generally, Root-MUSIC is more accurate than the normal MUSIC in a larger stepsize such as 1◦ , but for a smaller stepsize such as 0.1◦ , the MUSIC algorithm yields a higher accuracy than Root-MUSIC. In the presence of completely polarized electromagnetic wave, Costa [15] proposed a joint high-resolution 2D DOA and polarization estimator [16,17] based on polynomial rooting techniques and Fast Fourier decomposition of the steering vector. Several antenna arrays including small handheld terminals and channel sounding applications show a significant improvement in performances. Goossens and Rogier [18] proposed an efficient 2D DOA estimation algorithm based on the combination of UCA-RARE and Root-MUSIC. They estimated the azimuth angles by UCA-RARE which decoupled the estimation of azimuths from that of elevations. For each azimuth, they designed a new search-free rooting algorithm by expanding the array manifold into a double Fourier series. It is proved that even in the presence of severe mutual coupling [19,20], this hybrid approach can yield very robust 2D DOA estimation performances. Converting complex-valued computations to real-valued ones can often accelerate DOA estimations. Ren [21] preprocessed the block Hankel matrix in [11] through an averaging like method for unitary transformation, which transforms the eigenvalue and singular value decomposition into real-valued computation. Subsequently, the DOA estimations are solved by 1-D unitary ESPRIT [22]. The 2D unitary ESPRIT algorithm in [23] involves a specific unitary transformation to formulate the eigenvalue decomposition in terms of a real-valued computation. Motivated by these low-complexity algorithms, which have been successfully in applying sensor array processing [24], this paper proposed a fast 2D MUSIC algorithm, which involves three rounds of searching: the first two rounds are 1D searches along a circle and straight lines, respectively, and the third one is a 2D search on several concentric circles. The contribution of this work lies in the observation of spectrum peak diffusion effect and its utilization for finding the approximate direction of sources, while most of the other works in the literature concentrate on how to suppress this effect. The second contribution is the design of a search strategy to minimize the search path and avoid the total spectral search over the angular field of view. Simulation results demonstrate that with the same search stepsize, AESS has a similar estimated RMSE compared with the 2D MUSIC algorithm, while the computing complexity and the actual running time are remarkably reduced. Meanwhile, with the same computation cost, the estimated

Sensors 2017, 17, 515

3 of 14

RMSE of AESS is significantly improved. The 2D DOA estimation success rate is also analyzed versus the computation complexity. To further evaluate the estimation performance, we present the scatter plots of estimated 2D angle, and compare the running time of the proposed AESS with the 2D MUSIC algorithm and PIE algorithm. 2. Problem Formulation Consider a planar array composed of M (m = 1,2, . . . ,M) sensors on which P (p = 1,2, . . . ,P) narrow band non-coherent far field signals sp (t) are impinging with different DOAs. We denote the coordinates of the m-th sensor as (xm ,ym ,0), and the DOA related to the p-th signal is (θ p ,ϕp ), where θ p and ϕp are the azimuth and elevation, respectively. The uniform rectangular array (URA) coordinates system is shown in Figure 1. In the case of the m-th sensor receiving the p-th signal, there is a wave path difference between sensors. As illustrated by the coordinate system in Figure 1, if we choose the element at the origin of the coordinate system as the reference point, this path difference is calculated as: ∆mp = xm cos θ p sin ϕ p + ym sin θ p sin ϕ p

(1)

The complex factor corresponding to phase difference is calculated as: amp = exp (− j2π∆mp /λ)

(2)

where λ is the wavelength. At the t-th snapshot (t = 1,2, . . . ,T), the signal arriving at the m-th sensor is: P

xm (t) =

∑ amp s p (t) + nm (t)

(3)

p =1

where nm (t) is the White Gaussian noise. The array output at the t-th snapshot is modeled as: x(t) = [ x1 (t), x2 (t), · · · , x M (t)] T = As(t ) + n(t )

(4)

where T represents the transpose, A = [a1 ,a2 , . . . ,ap ] is the array direction matrix, ap = [exp(−j2π∆1p /λ), . . . ,exp(−j2π∆Mp /λ)]T is the steering vector, s(t) = [s1 (t),s2 (t), . . . ,sp (t)]T is the vector of source waveforms, and the noise vector n(t) = [n1 (t),n2 (t), . . . ,nm (t)]T is assumed to be zero-mean and independent of the observed signal. The correlation matrix of the received data is defined as: Rx = E[x(t)x H (t)] = ARs A H + σ2 I

(5)

where σ2 and I represent the noise power and unit matrix, respectively, and H represents the complex ˆ n. conjugate transpose. By the eigen-decomposition of Rx, we obtain the estimation of noise subspace U Utilizing the orthogonality between signal and noise subspace [25], the 2D MUSIC algorithm estimates the source DOA by a spectrum search as follows: 1 P(θ, ϕ) = a H (θ, ϕ)U ˆ nU ˆ nH a(θ, ϕ)

(6)

where a(θ,ϕ) is the steering vector corresponding to the azimuth angle θ and elevation angle ϕ. When the noise subspace Un is sufficiently accurate and (θ,ϕ) is equal to the actual DOA, the spatial spectrum function P(θ,ϕ) shows a peak. 2D MUSIC has to perform a 2D search with a certain stepsize over the angular field of view. A large stepsize may simplify the calculations but also lead to a worse RMSE. Furthermore, MUSIC possesses an extremely sharp spectrum peak, and a larger stepsize is more likely to miss the peak, especially in the case of a high signal-to-noise ratio (SNR). A small stepsize may improve RMSE, but may also increase computing overhead.

Sensors 2017, 17, 515

Sensors 2017, 17, 515

4 of 14

Sensors 2017, 17, 515

4 of 14

4 of 14

Z

s p (θ p , ϕ p )

s p ( p ,  p ) ϕp

p

xm

• •

Xxm

o• • • • oθ p    p

 

Y • • • y• m    ym

• • m ( xm , ym , 0) m

( x , y , 0)

m m Figure The coordinate coordinate system. Figure 1. 1. The system.

3. Spectrum Peak Diffusion Effect 3. Spectrum Peak Diffusion Effect

Figure 1. The coordinate system.

From Equation (6), we can see that if signal or noise subspace is sufficiently accurate and the 3. Spectrum Peak Diffusion From Equation (6), we canEffect see that if signal or noise subspace is sufficiently accurate and the estimated (𝜃𝜃�, 𝜑𝜑�) equals to the actual (θ,ϕ), the spatial spectrum will show a peak. Usually, when ˆ ˆ ϕˆ ) From (6), weactual can see thatthe if signal or spectrum noise subspace is sufficiently accurate and the(θ, estimated (θ, ϕˆ )Equation equals to the (θ,ϕ), spatial will show a peak. Usually, when (𝜃𝜃�, 𝜑𝜑�) deviates from the (θ,ϕ), there will still be a spatial spectrum peak although its peak becomes estimated ( , (θ,ϕ), ) equals to the actual (θ,aspatial ), the spatial spectrum will show aits peak. Usually, when deviates from the there still be spectrum peak although peak becomes smaller smaller with the increase of will deviation from the actual DOA. We call this the spectrum peak diffusion ( , increase ) deviates from the (θ,from ), there will still be a We spatial spectrum peak although its peak becomes witheffect. the of deviation the actual DOA. call this the spectrum peak diffusion effect. It is this observation that motivated us to change the traditional search strategy in 2D MUSIC, It is smaller with the increase of deviation from the actual DOA. We call this the spectrum peak diffusion this observation motivated us to change traditional searchatstrategy theDOA diffusion the diffusion that spectrum peak provides us twothe pieces of information least: oneinis2D thatMUSIC, the actual effect. It is this observation that motivated us to change the traditional search strategy in 2D MUSIC, is near,peak and the other is us thetwo approximate of the sources. spectrum provides pieces ofdirection information at least: one is that the actual DOA is near, the diffusion spectrum peak provides us two pieces of information at least: one is that the actual DOA Figureis2the shows the spatial spectrum diffusion angle with 100 Monte Carlo trials. Here and the other approximate directionversus of thethe sources. is near, and the other is the approximate direction of the sources. aFigure uniform rectangular array with three rows and four columns omnidirectional sensors is used andHere a 2 shows the the spatial spectrum versus with100 100Monte Monte Carlo trials. Figure 2 shows spatial spectrum versusthe thediffusion diffusion angle angle with Carlo trials. Here a narrow band non-coherent far field signal comes from (45°,45°). The SNR, snapshot and stepsize uniform rectangular arrayarray withwith three rows and four columns sensors used a uniform rectangular three rows and four columnsomnidirectional omnidirectional sensors is is used andand a are set to be 0 dB, 400 and 1°, respectively. The simulation in◦ Figure 2 manifests that the spatial ◦ narrow band non-coherent far field signal comes from (45(45°,45°). ,45 ). The and stepsize stepsize are a narrow band non-coherent far field signal comes from TheSNR, SNR,snapshot snapshot and spectrum peak is high◦ enough to be detected even the deviation is far from the actual (θ,ϕ). to 400 be 0and dB, 1400 and 1°, respectively. The simulation in Figure 2 manifests that the spatial set toare be 0setdB, , respectively. The simulation in Figure 2 manifests that the spatial spectrum spectrum peak is high enough to be detected even the deviation is far from the actual (θ,  ). peak is high enough to be detected even the deviation is far from the actual (θ,ϕ).

Figure 2. Spatial spectrum versus the diffusion angle. Figure 2. Spatial spectrum versus the diffusion angle. 4. Proposed Fast Algorithm Figure 2. Spatial spectrum versus the diffusion angle.

SupposeFast thereAlgorithm are P narrow band non-coherent far field signals sp(t) (p= 1,2,…,P), the azimuth 4. Proposed

4. Proposed Fast angle Algorithm and elevation of the p-th source are set to be (θp,ϕp), where θp ∈ (0°,Θ) and ϕp ∈ (0°,Φ). Figure 3

there are Pofnarrow band non-coherent far DOA field signals showsSuppose the search scheme the proposed AESS in the 2D plane. sp(t) (p= 1,2,…,P), the azimuth Suppose thereangle are Pofnarrow field signals sp (t) (pand = 1,2, . . ,P), Figure the azimuth and elevation the p-th band sourcenon-coherent are set to be (θfar p,p), where θp  (0°,Θ) p . (0°,Φ). 3 ◦ ,Θ) and ϕ ∈ (0◦ ,Φ). Figure 3 shows theangle searchofscheme of source the proposed 2D DOA plane. and elevation the p-th are setAESS to bein(θthe ,ϕ ), where θ ∈ (0 p p p p

shows the search scheme of the proposed AESS in the 2D DOA plane.

2017, Sensors Sensors 2017, 17, 51517, 515

5 of 14 5 of 14

ϕ Φ



(θ p − Θ 2) tan δ p + Φ 2 ϕ= p

∗ sp

•

∗ s2

(θ p , ϕ p )



o

δp

(Θ 2 , Φ 2)

Ω

∗ s1 ∗ sP 0 Θ

θ

Figure 3. Proposed search scheme. *s1 and *s2 represent signals 1 and 2. Figure 3. Proposed search scheme. * s1 and * s2 represent signals 1 and 2.

Our proposal is a is three rounds search-based In the thefirst firstsearch, search, select point Our proposal a three rounds search-basedscheme. scheme. In wewe select the the point o(Θ/2,Φ/2) as the center a circle.Let LetΩΩbe be the the radius, radius, which not be be tootoo large or too small. o(Θ/2,Φ/2) as the center of of a circle. whichshould should not large or too small. An appropriate should be selected Ω order in order modestlysitting sittingthe thecircle circle on on the the area An appropriate valuevalue should be selected for for Ω in to to modestly area of of 2D plane. Itbe should noted theofunit Ω islength not length degree, parametricequation equation of of this plane.2D It should notedbethat thethat unit Ω isofnot but but degree, so so thethe parametric this circle is: circle is: ( θ= δ+ δ +Θ/2 Ω cos Θ2 θΩ= cos ◦ 0◦ , 360 °] ] (7) (7)  Ω sin δ + Φ/2 , ,δδ∈∈(0(°,360 ϕ= ϕ = Ω sin δ + Φ 2 wherewhere δ is the angle. TheThe first search isisimplemented on this thiscircle, circle,since since is selected before δ iscentral the central angle. first search implemented on ΩΩ is selected before theround first round of search, it is actually a 1D search. The search variable is neither azimuth nor the first of search, it is actually a 1D search. The search variable is neither azimuth nor elevation, elevation, the rotation δ which be Figure seen from Figure 3. The corresponding 1D spectral but the rotationbut angle δ whichangle can be seen can from 3. The corresponding 1D spectral search is search written as: is written as: 1 1 (8) PNO−1 = PNO-1 = a HH (δ)U ˆ nU ˆHnH a(δ) (8) ˆ ˆ a (δ )U n U n a(δ ) After searching on the circle with certain stepsize from δ = 0◦ to δ = 360◦ , and due to the existence After searching on the circle with certain stepsize from = 0° get to δ P = 360°, and due to thewhose existence of the spectrum peak diffusion effect mentioned above, weδcan spectrum peaks search of the spectrum peak diffusion effect mentioned above, we can get P spectrum peaks whose search ˆ ˆ variables are δp (p =̂ 1,2, . . . ,P). Corresponding to each δ̂p , we can get a straight line whose parametric variables are 𝛿𝛿𝑝𝑝 (p= 1,2,…,P). Corresponding to each 𝛿𝛿𝑝𝑝 , we can get a straight line whose parametric equation is: equation is:

_

ϕ p = (θ p − Θ/2) tan δ p + Φ/2, ( p = 1, 2, · · · , P) (θ p − Θ 2) tan δ p + Φ 2 , ( p = 1, 2, , P) ϕ= p And the length of the p-th line (p = 1,2, . . . ,P) is: And the length of the p-th line (p= 1,2,…,P) is:  _    [2sin( ( θ pθ)]p )], ,  Φ/[2Φsin  L p = Lp =  _  Θ [2 cos(θ )],   Θ/[2 cos( θ p )]p ,

_  tan > tan δ δp p Φ >Θ Φ/Θ  tan δ_p < Φ Θ tan δ p < Φ/Θ

(9)

(9)

(10) (10)

Similarly, the unit of Lp is degrees. The second round of searching is implemented along the P straight is also which isround expressed as: Similarly,lines; the itunit of Lap1D is spectral degrees.search The second of searching is implemented along the P

straight lines; it is also a 1D spectral search expressed as: 1 which is , l p ∈ ( 0, L p  , ( p = 1, 2,  , P) PNO-2 = H  (11) H ˆ U ˆ a(l ) a (l p )U  1 n n p , l p ∈ 0, L p , ( p = 1, 2, · · · , P) PNO−2 = (11) a Hwe ˆ nget ˆ nHone (l pcan )U U a(l pspectrum ) Applying Equation (11), peak on each straight line. For each spectrum peak, the corresponding search variable in Equation (11) is marked as 𝑙𝑙̃𝑝𝑝 . Now we can derive the

Applying Equation (11), we can get one spectrum peak on each straight line. For each spectrum peak, the corresponding search variable in Equation (11) is marked as e l p . Now we can derive the e ep ) which corresponds to the spectrum peak on the p-th straight line, azimuth and elevation (θ p , ϕ ep is expressed as: the parametric equation of θep and ϕ

Sensors 2017, 17, 515

6 of 14

 _  θe = e l p cos δ p + Θ/2 p , ( p = 1, 2, · · · , P) _  ϕ e =e l sin δ + Φ/2 p

p

(12)

p

ep ) (p = 1,2, . . . ,P) is adopted, we can perform the third round of searching which Once the (θep , ϕ ep ) (p = 1,2, . . . ,P), is a 2D search on a set of concentric circles, the centers of the circles are (θep , ϕ ◦ ◦ respectively. This time the search variable are rotation angle τ ∈ (0 ,360 ] and the neighborhood radius is ω. The spectrum peak function is given as:  1 , τp ∈ (0◦ , 360◦ ], ω p ∈ 0, Wp PNO−3 = H H a (τp , ω p )U ˆ nU ˆ n a(τp , ω p )

(13)

Performing the above 2D search, we can obtain the search variable corresponding to the spectrum

peaks which can be noted as τˆp , ωˆ p (p = 1,2, . . . ,P), from τˆp , ωˆ p we can obtain the estimation of the actual 2D DOA as: ( θˆ p = ωˆ p cos τˆp + θep , ( p = 1, 2, · · · , P) (14) ep ϕˆ p = ωˆ p sin τˆp + ϕ 5. Computational Complexity Analysis We carefully compute the flops required by the classical 2D MUSIC and the proposed AESS methods. As mentioned above, M, P and T represent the total sensor number, source number and snapshot, respectively. Assume K represents the stepsize. Firstly, we present several basic operations and their computation consume in flops. Computing the multiplication of a M row and P column complex matrix by a P row N column complex matrix takes 6MN(2P − 1) flops. Computing the mathematical expectation of a M row M column matrix takes 6M2 + 1 flops. Computing the eigen-decomposition of a M row M column matrix takes 4M5 − 5M4 + M3 flops. Computing the absolute of a complex scalar takes 3 flops. Computing the reciprocal of a scalar takes 1 flops. For 2D MUSIC, assume the stepsize is K. Computing Rxx = E[xxH ] takes 8M2 T − 2M2 + 1 ˆ nU ˆ nH takes flops. Computing Rxx = UΣUH takes 4M5 − 5M4 + M3 flops. Computing G N = U 1 8M3 − 8M2 P − 2M2 flops. Computing P(θ, ϕ) = H takes (6M2 + 16M + 52)(Θ/K)(Φ/K) |a (θ,ϕ)GN a(θ,ϕ)| flops, where Θ and Φ are the maximum search scope of azimuth and elevation which are depicted in Figure 3. The total computation consumption in flops is: C2D MUSIC

= (6M2 + 16M + 52)(ΘΦ/K2 ) +8M2 T + 4M5 − 5M4 + 9M3 − 8M2 P − 4M2 + 1

(15)

For the proposed AESS algorithm, the first three steps also involve computing ˆ nU ˆ nH , so the corresponding computation consumptions Rxx = E[xxH ], Rxx = UΣUH and G N = U 1 are identical to 2D MUSIC. The following step of computing PNO−1 = in | a H ( δ )G N a( δ ) | 1 2 Equation (8) takes (6M + 16M + 52)(360/K) flops. Computing PNO−2 = |aH (l p )GN a(l p )|
P
2 in Equation (11) takes 6M + 16M + 52 ∑ p=1 L p /K flops, where the meaning of Lp can 1 be seen from Equation (10).Computing PNO−3 = in Equation (13) takes |aH (τp ,ω p )GN a(τp ,ω p )|

6M2 + 16M + 52 ∑ Pp=1 (360/K ) Wp /K flops, where the meaning of Wp can be seen from Equation (13). The total computation consumption in flops is:

" C AESS =

(6M2

P

P

p =1

p =1

+ 16M + 52) (360/K ) + ∑ ( L p /K ) + ∑ (360Wp +8M2 T + 4M5 − 5M4 + 9M3 − 8M2 P − 4M2 + 1

# /K2 )

(16)

Sensors 2017, 17, 515

7 of 14

From Equations (15) and (16), we can see that the sensor number M and stepsize K are the Sensors important 2017, 17, 515 factors in determining the computational complexity. Snapshot T7 and of 14 the two most number of sources P are not the principal factors that increase the calculation. Figure 4 shows From Equations (15) and (16), we can see that the sensor number M and stepsize K are the two the computational complexity varies with the sensor number and the stepsize. We can see that the most important factors in determining the computational complexity. Snapshot T and the number of computational complexity of AESS and 2D MUSIC algorithms both increase with the decrease of sources P are not the principal factors that increase the calculation. Figure 4 shows the computational stepsize. The increase computational complexity of stepsize. 2D MUSIC is much faster than that complexity varies of with the sensor number and the We algorithm can see that the computational of thecomplexity proposed AESS. AESS thus has a remarkable computational advantage over the conventional of AESS and 2D MUSIC algorithms both increase with the decrease of stepsize. The 2D MUSIC method, especiallycomplexity in the caseofof 2D small stepsize. Filik [21] proposed fastthat automatically increase of computational MUSIC algorithm is much faster a than of the AESS. AESS thus has a remarkable computational advantage overthe therotational conventional 2D by pairedproposed 2D DOA estimation algorithm termed PIE, where they exploited matrix MUSICreal method, in the case small stepsize. Filikeigenvalues [21] proposedofa fast pairedhave employing and especially virtual arrays, andofdeduced that the the automatically rotational matrix 2D DOA estimationat algorithm termed PIE,and where theywhich exploited the rotational matrix by the angle information both magnitude phase allows the estimation ofemploying azimuth and real and virtual arrays, and deduced that the eigenvalues of the rotational matrix have the angle elevation angles via closed-form expressions. Here we analyze the computational flops of the information at both magnitude and phase which allows the estimation of azimuth and elevation PIE algorithm with the symbols of [21]. Computing y4 = (B12 + B13 )y1 takes 2M2 + 8M2 P flops, angles via closed-form expressions. Here we analyze the computational flops of the PIE
− algorithm   ˆ = Sˆ H Sˆ 1 1 Sˆ H Sˆ 4 takes computing Ry = E y5 , y H takes 128M5 + 80M4 + 8M3 flops, computing Φ 2 2 with the symbols of [21]. Computing y4 = (B12 + B13)y1 takes 2M + 8M P flops, computing 5 5 1 1𝐑𝐑 𝐲𝐲5 =

2 flops, 5 + 80M 4 + 8M 3 flops, computing 18PM2𝐸𝐸[𝐲𝐲 −56PM 6P3 +128M (8/3)P computing ϕi = 𝚽𝚽 arctan νi𝐒𝐒�)1𝐻𝐻/arccos νi |/2)2 −] 6PM takes+ 6P 9P 3+flops, � = (𝐒𝐒�[1𝐻𝐻arg , 𝐲𝐲5𝐻𝐻 ] +takes 𝐒𝐒�1 )(−1 𝐒𝐒�4 takes(|18PM r n o 8/3P2 flops, computing 𝜑𝜑𝑖𝑖 = arctan[arg(𝜈𝜈 2 𝑖𝑖 )/arccos(|𝜈𝜈𝑖𝑖 |/2) ] takes 29P flops, computing 𝜃𝜃𝑖𝑖 = computing θi = arcsin [λ arg(νi )/2πd] + [(λ arccos(νi /2)/2πd] takes 27P flops. The total

arcsin�{[arg(𝜈𝜈𝑖𝑖 )/2𝜋𝜋𝜋𝜋]2 + [(𝜆𝜆 arccos(𝜈𝜈𝑖𝑖 /2)/2𝜋𝜋𝜋𝜋]2 }

computation consumption consumption in flops is:in flops is:

takes

27P

flops.

The

total

computation

8 28 2 3 2 3 5 M 5 + 804M 4 + 8M 26 P + = 128 + 22))M + 6+P6P +3 + P +P36 + P 36P CPIE =CPIE 128M + 80M + 8M3 ++((26P M2−−6 PM 6PM 3 3

(17) (17)

Figure 4. Computational complexity the sensor sensornumber number and stepsize. Figure 4. Computational complexityversus versus the and stepsize.

6. Simulation Results 6. Simulation Results In this section, simulations are presented to illustrate the performance of the proposed fast 2D

In this section, simulations are presented to illustrate the performance of the proposed fast 2D MUSIC algorithm. A uniform rectangular array (URA) with omnidirectional sensors is used, along MUSIC algorithm. A uniform rectangular array (URA) with omnidirectional sensors is used, along with three narrow band non-coherent far field signals coming from (20°,40°), (75°,50°) and (35°,60°). with three narrow band non-coherent far field signals coming from (20◦ ,40◦ ), (75◦ ,50◦ ) and (35◦ ,60◦ ). 6.1. RMSE and Estimation Success Rate versus Computational Complexity

6.1. RMSE and Estimation Success Rate versus Computational Complexity The root mean square errors (RMSE) is defined as:

The root mean square errors (RMSE) is defined as: 2 2 1 N 1 P ˆ (degree)  + ϕˆ p (n) − ϕ p  v RMSE= ( ) θ n − θ (18) ∑ ∑ p p     u N P o P pn1 = n 1= u1 N    1 2 2 RMSE = t ∑ ∑ θˆ p (n) − θ p + ϕˆ p (n) − ϕ p (degree) (18) N nindependent where N = 100 represents the trials, P = 3 represents the number of sources, 𝜃𝜃�𝑝𝑝 (𝑛𝑛) is the =1 P p =1

{

}

estimation of θp for the n-th Monte Carlo trial. The results have been averaged over 100 Monte Carlo whereruns. N = 100 represents the independent trials, P = 3 represents the number of sources, θˆ p (n) is It is obvious that the RMSE depends on the accuracy of the estimated subspace. Moreover, the the estimation of θ p for the n-th Monte Carlo trial. The results have been averaged over 100 Monte estimated subspace is usually derived from the eigen-decomposition calculation. As this kind of Carlo runs. calculation is critically affected by SNR and snapshots, the performance of RMSE is always closely

Sensors 2017, 17, 515

8 of 14

It is obvious that the RMSE depends on the accuracy of the estimated subspace. Moreover, the estimated subspace is usually derived from the eigen-decomposition calculation. As this kind of calculation is critically affected by SNR and snapshots, the performance of RMSE is always closely related to the SNR and snapshots. In this paper, both the classic 2D MUSIC and the proposed AESS are Sensors 2017, 17, 515 8 of 14 search-based methods, procedure, Sensors 2017, 17, 515 they use the same subspace and their main difference is the search 8 of 14 so SNRrelated and snapshots have the same influence on RMSE. For this reason, our simulations here to the SNR and snapshots. In this paper, both the classic 2D MUSIC and the proposed AESS analyze related to thesensor SNR and snapshots. In this paper, both the classic 2D MUSIC andthe the proposed AESSarray of the RMSE versus number equal stepsizes. We respectively are search-based methods, theyinuse the same subspace and their main use differencerectangular is the search are search-based methods, they use the same subspace and their main difference is the search procedure, so SNR and snapshots have the same influence on RMSE. For this reason, our simulations three row and three column with a total of 3 × 3 = 9 omnidirectional sensors, and the following procedure, so SNR and snapshots have the same influence on RMSE. For this reason, our simulations here analyze the RMSE versus sensor number in equal stepsizes. We respectively use the rectangular are 3 ×here 4 =analyze 12, 3 × 5 = 15, 3 × 6 = 18, 3 × 7 = 21 and 3 × 8 = 24 omnidirectional sensors. the RMSE versus sensor number in equal stepsizes. We respectively use the rectangular From array of three row and three column with a total of 3 × 3 = 9 omnidirectional sensors, and the following Equations (15) and (16),and wethree know thatwith the asensor M and snapshot K take anfollowing important role array of three row column total of number 3 × 3 = 9 omnidirectional sensors, and the are 3 × 4 = 12, 3 × 5 = 15, 3 × 6 = 18, 3 × 7 = 21 and 3 × 8 = 24 omnidirectional sensors. From Equations are 3 × 4 = 12, 3 × 5 = 15, 3 × 6 =of 18,the 3 × 2D 7 = 21 and 3 ×and 8 = 24 omnidirectional sensors. From Equations in the computational MUSIC substituting (15) and (16), wecomplexity know that the sensor number M and proposed snapshot K AESS, take anso important role inthe thevalue of (15) and (16), know that(15) the and sensor number Maand snapshot K takewe an can important role in the sensor computational number into we Equations (16) under certain stepsize, get the corresponding complexity of the 2D MUSIC and proposed AESS, so substituting the value of sensor computational complexity of the 2D MUSIC and proposed AESS, so substituting the value of sensor computational complexity, we(16) canunder derive the RMSE versus complexity as number into Equations hence (15) and a certain stepsize, we the can computation get the corresponding number into Equations (15) and (16) under a certain stepsize, we can get the corresponding computational complexity, hence we can derive the RMSE versus the computation complexity as shown in Figure 5. In Figure 5, SNR is set to be 0 dB, snapshot is set to be 400, and the stepsize is computational complexity, hence we can derive the RMSE versus the computation complexity as shown in Figure 5. In Figure 5, SNR is set to be 0 dB, snapshot is set to be 400, and the stepsize is 0.1°. ◦ 0.1 . From Figure 5, we can see that if the stepsizes are equal, the classic 2D MUSIC and our proposed shown in Figure 5. In Figure 5, SNR is set to be 0 dB, snapshot is set to be 400, and the stepsize is 0.1°. From Figure 5, we can see that if the stepsizes are equal, the classic 2D MUSIC and our proposed AESS possess similar RMSE, while computational complexity of the proposed will be greatly From Figure 5, we can see thatthe if the stepsizes are equal, the classic 2D MUSIC andAESS our proposed AESS possess similar RMSE, while the computational complexity of the proposed AESS will be AESS possess similar RMSE, whilecomplexity the computational complexity of the proposed AESSRMSE will be decreased. With equal computational provided, AESS will derive a better greatly decreased. With equal computational complexity provided, AESS will derive a better RMSEthan the greatly decreased. In With equal computational complexity provided, AESSestimation will derive asuccess better RMSE 2D MUSIC algorithm. this experiment, we also validate the DOA rate versus than the 2D MUSIC algorithm. In this experiment, we also validate the DOA estimation success rate than the 2D MUSIC algorithm. In this experiment, we also validate the DOA estimation success rate versus computational SNR isat fixed at 0and dB and snapshot 400, stepsize stepsize isissetset to to be be 0.1◦ . computational complexity,complexity, the SNR the is fixed 0 dB snapshot atat400, versus computational complexity, the SNR is fixed at 0 dB and snapshot at 400, stepsize is set to be � (𝑛𝑛), 0.1°. In each independent run, we decide that the estimated �𝜃𝜃 (𝑛𝑛)� are successfully resolved 𝜑𝜑 � In each 0.1°. independent run, we decide that the estimated [θˆ p (n�𝜃𝜃�)𝑝𝑝, (𝑛𝑛), ϕˆ p (𝜑𝜑�n𝑝𝑝 )] are successfully resolved if: In each independent run, we decide that the estimated (𝑛𝑛)� are successfully resolved if: if:

𝑝𝑝

𝑝𝑝

o θˆ p (n)ˆ− θ p (n) + ϕˆ p (n) − ϕ p (n) ≤ θˆp (n) − θ p (n) + ϕˆ p (n) − ϕ p (n) ≤ 0.5oo 0.5 θ p (n) − θ p (n) + ϕˆ p (n) − ϕ p (n) ≤ 0.5

Figure 5. RMSE versus computational complexity.

Figure 5. RMSE complexity. Figure 5. RMSEversus versus computational computational complexity.

Figure 6. Success rate versus computational complexity. Figure 6. Success rate versus computational complexity.

Figure 6. Success rate versus computational complexity.

(19) (19)

(19)

Sensors Sensors2017, 2017,17, 17,515 515

9 9ofof1414

Figure 6 demonstrates that the 2D MUSIC algorithm need many more calculation flops to obtain Figuresuccess 6 demonstrates that the algorithm needanother many more flops obtain the same rate compared to 2D theMUSIC proposed AESS. From pointcalculation of view, we canto see AESS the same success rate compared to the proposed AESS. From another point of view, we can see AESS possesses much higher success rate in the same calculation flops provided. The numerical results possesses much rate flops in the provided, same calculation flopscover provided. The numerical show show that withhigher 39 to success 250 million AESS can the success rate fromresults 36% to 94% that with 39 to 250 million flops provided, AESS can cover the success rate from 36% to 94% while 2D while 2D MUSIC needs 554 to 3192 million flops to achieve the same success rate. MUSIC needs 554 to 3192 million flops to achieve the same success rate. 6.2. RMSE versus SNR and Snapshot 6.2. RMSE versus SNR and Snapshot Figures 7 and 8 describe the RMSE versus SNR and snapshot, respectively. Here a rectangular Figures 7 and 8 describe the RMSE versus SNR and snapshot, respectively. Here a rectangular array of four rows and four columns with a total of sixteen omnidirectional sensors are used. In Figure array of four rows and four columns with a total of sixteen omnidirectional sensors are used. In Figure 7 7 the snapshot is set to be 400, the SNR varies from −5 dB to 19 dB with 3 dB interval. the snapshot is set to be 400, the SNR varies from −5 dB to 19 dB with 3 dB interval.

(a)

(b)

(c) Figure 7. RMSE versus SNR. Snapshot is set to be 400, stepsize is set to be 0.1 . (a) The first source Figure 7. RMSE versus SNR. Snapshot is set to be 400, stepsize is set to be 0.1◦ . (a) The first source come from (20 , 40) ; (b) The second source come from (75, 50) ; (c) The third source come from come from (20◦ , 40◦ ) ; (b) The second source come from (75◦ , 50◦ ) ; (c) The third source come from . (35  , 60  ) (35◦ , 60◦ ).

In Figure 8 the SNR is set to be 0 dB, the snapshot varies from 300 to 1000 with 100 interval. The In Figure 8 the SNR is set to be 0 dB, the snapshot varies from 300 to 1000 with 100 interval. stepsize of 2D MUSIC and AESS in this experiment are all set to be 0.1°. The results have been The stepsize of 2D MUSIC and AESS in this experiment are all set to be 0.1◦ . The results have been averaged over 500 Monte Carlo runs. Since our proposed AESS and classic 2D MUSIC both are averaged over 500 Monte Carlo runs. Since our proposed AESS and classic 2D MUSIC both are search-based methods, their RMSEs mainly depend on the stepsize, the two figures show that they search-based methods, their RMSEs mainly depend on the stepsize, the two figures show that they possess almost the same RMSE versus SNR and snapshot. It should be noted that in the higher SNR possess almost the same RMSE versus SNR and snapshot. It should be noted that in the higher SNR scenario, the PIE algorithm derives a significant RMSE improvement. scenario, the PIE algorithm derives a significant RMSE improvement.

Sensors Sensors2017, 2017,17, 17,515 515

1010ofof1414

(a)

(b)

(c) Figure8.8.RMSE RMSEversus versussnapshot. snapshot.SNR SNR isisset set to tobe be00dB, dB,stepsize stepsize is is set set to to be be 0.1 0.1 (a)The Thefirst firstsource source ◦ . .(a) Figure come from ; (b) The second source come from ; (c) The third source come from (20  , 40  ) (75  , 50  ) come from (20◦ , 40◦ ) ; (b) The second source come from (75◦ , 50◦ ) ; (c) The third source come from (35 60◦)).. ◦ , ,60 (35

6.3.Comparison ComparisonScatter ScatterPlots PlotsofofEstimated Estimated2D 2DAngles AnglesininDifferent DifferentSNR SNR 6.3. RMSEcan cancoarsely coarselyreflect reflectthe theestimated estimatedaccuracy. accuracy.For Forfurther furtheranalyzing analyzingthe theperformances performancesofofthe the RMSE estimator, we show the scatter plots of the DOA in Figure 9. Here the SNR of three sources are set estimator, we show the scatter plots of the DOA in Figure 9. Here the SNR of three sources are set toto be−−5 dB, 00 dB dB and and 10 10 dB dB respectively. respectively.A Arectangular rectangulararray arrayof offour fourrow rowand andfour fourcolumn columnwith withtotally totally be 5 dB, sixteen omnidirectional sensors are used. The stepsize of 2D MUSIC and AESS in this experiment are sixteen omnidirectional sensors are used. The stepsize of 2D MUSIC and AESS in this experiment are allset settotobe be0.1 0.1°. Thesnapshots snapshotsofofthe thethree threealgorithms algorithms are 400. The cross symbols Figure ◦ . The all are setset toto bebe 400. The cross symbols in in Figure 9 9 represent the actual position of sources. represent the actual position of sources. Figure9a–d 9a–ddemonstrate demonstratethat thataahigher higherSNR SNRmakes makeseach eachofofthe thealgorithms algorithmsderive deriveaasmaller smallerDOA DOA Figure estimationbias biasand andsmaller smallervariance. variance.On Onthe theother otherhand, hand,lower lowerSNR SNRmakes makesitithard hardtotodetect detectthe theDOAs. DOAs. estimation We can see that AESS has a higher uniformity than the 2D MUSIC and PIE algorithms. That is to We can see that AESS has a higher uniformity than the 2D MUSIC and PIE algorithms. That is to say,say, its its scatter plots concentrate on a smaller area than the other two algorithms. Unfortunately, the scatter plots concentrate on a smaller area than the other two algorithms. Unfortunately, the estimation is the biased from the actual DOAs. The the reason is that ofthe search of our AESSonis isestimation biased from actual DOAs. The reason is that last search ourlast AESS is implemented implemented on several small concentric circles, when the actual source is situated the several small concentric circles, when the actual source is situated outside the biggest outside concentric biggest concentric circle, AESS also take the highest spectrum peak inside thecircle biggest concentric circle circle, AESS also take the highest spectrum peak inside the biggest concentric as the actual DOA. as the actual DOA. Obviously, this spectrum peak is the highest within the small area of concentric Obviously, this spectrum peak is the highest within the small area of concentric circles, but not the circles,spectrum but not the highest spectrum peak field in theoftotal angular field of view. This problem is inevitable highest peak in the total angular view. This problem is inevitable for almost all local for almost all local search methods. A slightly enlargement of the radius of the concentric circles can search methods. A slightly enlargement of the radius of the concentric circles can significantly decrease significantly decrease the estimation bias of AESS from the actual DOAs, which we can see from the estimation bias of AESS from the actual DOAs, which we can see from Figure 9d. The simulation Figure 9d.of The simulation conditions Figure are same as forradius Figureis9cmodified except the search conditions Figure 9d are same as forof Figure 9c9d except the search from 5◦ toradius 7.5◦ . is modified from 5° to 7.5°.

Sensors 2017, 17, 515

11 of 14

Sensors 2017, 17, 515

11 of 14

(a)

(b)

(c)

(d)

Figure 9. Scatter plots of estimated 2D angles when three sources with different SNR level are Figure 9. Scatter plots of estimated 2D angles when three sources with different SNR level are impinging impinging on the URA. (a) Classic 2D MUSIC algorithm; (b) PIE algorithm; (c) The proposed AESS on the URA. (a) Classic 2D MUSIC algorithm; (b) PIE algorithm; (c) The proposed AESS algorithm; algorithm; (d) The proposed AESS algorithm with a slightly enlargement of the radius of the third (d) The proposed AESS algorithm with a slightly enlargement of the radius of the third round searching. round searching.

6.4. 6.4. CPU CPU Time Time Comparison Comparison In fairness In fairness to to both both algorithms, algorithms, we we compare compare the the search search time time of of the the proposed proposed AESS AESS and and the the classic classic 2D MUSIC algorithm under same stepsize conditions. The results are obtained using a PC with 2D MUSIC algorithm under same stepsize conditions. The results are obtained using a PC with an an Intel(R) Intel(R) Core(TM) Core(TM) i5-3470 i5-3470 3.2 3.2 GHz GHz CPU CPU and and 88 GB GB RAM RAM by by running running the the Matlab Matlab codes codes in in the the same same ◦ ◦ ◦ environment. environment. Table Table 11 shows shows that that in in the the different differentstepsizes stepsizessuch suchas as0.1 0.1°,, 0.25 0.25° and and 0.5 0.5°,, the the search search time time of of the the proposed proposed AESS AESS is is much much smaller smaller than than that that of of the the 2D 2D MUSIC MUSIC algorithm. algorithm. It It also also can can be be seen seen that that this this advantage advantage is is more more obvious obvious with with the the decrease decrease of of stepsize. stepsize. As As PIE PIE is is aa search-free search-free method, method, it it needs needs much much less less CPU CPU time time for for the the DOA DOA estimation. estimation. Table 1. Comparison of CPU time (s). Ss in the table means stepsize. Table 1. Comparison of CPU time (s). Ss in the table means stepsize. Configuration Ss 0.1° Ss 0.1◦ Configuration Row Row × Column × Column AESS AESS 3 × 33 × 3 2.4032 2.4032 3 × 43 × 4 2.4063 2.4063 3 × 53 × 5 2.4176 2.4176 3 × 63 × 6 2.4221 2.4221 2.4399 3 × 73 × 7 2.4399 2.5659 3 × 83 × 8 2.5659

Ss0.1° 0.1◦ Ss MUSIC MUSIC 41.9570 41.9570 41.9736 41.9736 43.4240 43.4240 43.9129 43.9129 44.8107 44.8107 45.4299 45.4299

Ss 0.25° 0.25◦ Ss AESS AESS 0.4205 0.4205 0.4221 0.4221 0.4228 0.4228 0.4248 0.4248 0.4398 0.4398 0.4553 0.4553

◦ Ss Ss 0.25 0.25° MUSIC MUSIC 6.6706 6.6706 6.7262 6.7262 6.8137 6.8137 6.9126 6.9126 7.0829 7.0829 7.7317

◦ SsSs 0.50.5° AESS AESS 0.1228 0.1228 0.1236 0.1236 0.1238 0.1238 0.1241 0.1241 0.1351 0.1351 0.1525 0.1525

◦ Ss 0.5 Ss 0.5°

MUSIC PIE PIE MUSIC 1.72630.00430.0043 1.7263 1.73940.00500.0050 1.7394 1.77280.00520.0052 1.7728 1.79550.00780.0078 1.7955 1.8056 1.80560.00890.0089 1.8943 1.89430.00950.0095

6.5. RMSE and CRB versus Source Separation δ To evaluate the resolution, we examine the capability of the proposed AESS and the related algorithms to estimate the DOAs of closely spaced signals. We fixed the first DOA source as (10°,40°),

Sensors 2017, 17, 515

12 of 14

6.5. RMSE and CRB versus Source Separation δ Sensors 2017, 17, 515

12 of 14

To evaluate the resolution, we examine the capability of the proposed AESS and the related algorithms to estimate the DOAs of closely spaced signals. We fixed the first DOA source as (10◦ ,40◦ ), and the second source is situated at (θ2,2) = (10° + δ,40°  δ), where δ = 0°, 1°, 2°, 3°, 5°, 7°, 11°, 15°, and the second source is situated at (θ 2 ,ϕ2 ) = (10◦ + δ,40◦ − δ), where δ = 0◦ , 1◦ , 2◦ , 3◦ , 5◦ , 7◦ , 11◦ , 20°, 25°, respectively. From Figure 10 we can see that when δ = 0°, there is only one source and the 15◦ , 20◦ , 25◦ , respectively. From Figure 10 we can see that when δ = 0◦ , there is only one source and RMSE is lower. When δ = 0°, 1°, 2°, 3°, 5°, 7°, the two source influence each other due to their close the RMSE is lower. When δ = 0◦ , 1◦ , 2◦ , 3◦ , 5◦ , 7◦ , the two source influence each other due to their spacing, and the RMSE in this section are worse. When δ = 11°,◦ 15°,◦ 20°,◦ 25°, the RMSE of all the ◦ close spacing, and17,the of all the Sensors 2017, 515 RMSE in this section are worse. When δ = 11 , 15 , 20 , 25 , the RMSE 12 of 14 algorithms are improved for the ease of distinguishing two far spaced signals. algorithms are improved for the ease of distinguishing two far spaced signals. and the second source is situated at (θ2,ϕ2) = (10° + δ,40° − δ), where δ = 0°, 1°, 2°, 3°, 5°, 7°, 11°, 15°, 20°, 25°, respectively. From Figure 10 we can see that when δ = 0°, there is only one source and the RMSE is lower. When δ = 0°, 1°, 2°, 3°, 5°, 7°, the two source influence each other due to their close spacing, and the RMSE in this section are worse. When δ = 11°, 15°, 20°, 25°, the RMSE of all the algorithms are improved for the ease of distinguishing two far spaced signals.

(a)

(b)

Figure10. 10.RMSE RMSEand andCRB CRBversus versussource sourceseparation. separation.(a) (a)RMSE RMSEofofthe thefirst firstsource; source;(b) (b)RMSE RMSEofofthe the Figure second source. (a) (b) second source. Figure 10. RMSE and CRB versus source separation. (a) RMSE of the first source; (b) RMSE of the

6.6. Performance in the Presence of Mutual Coupling second 6.6. Performance insource. the Presence of Mutual Coupling As the DOA estimation performance is often affected by mutual coupling [18], we simulated the As 6.6. thePerformance DOA estimation performance is often affected by mutual coupling [18], we simulated in the Presence of Mutual Coupling performance of our proposed and the related algorithms in the presence of mutual coupling. A the performance of our proposed and the algorithms in coupling the presenceweofsimulated mutual coupling. the DOA estimation isrelated often affected mutual rectangularAsarray of four rowsperformance and four columns with by a total of sixteen [18], omnidirectionalthe sensors is A rectangular array of four rows and four columns with a total of sixteen omnidirectional is performance of our proposed and the related algorithms in the presence of mutual coupling.sensors A used. The stepsizes of 2D MUSIC and AESS in this experiment are all set to be 0.1°. The snapshots of ◦ rectangular of four rowsand and AESS four columns a total ofare sixteen omnidirectional sensors is used. The stepsizesarray of 2D MUSIC in this with experiment all set to be 0.1 . The snapshots of the three algorithms areofall set to beand 400. Theinmutual coupling coefficient in our is same The stepsizes 2D AESS this experiment all set to be The experiment snapshotsisofsame the threeused. algorithms are all setMUSIC to be 400. The mutual coupling are coefficient in 0.1°. our experiment as as the literature [19], i.e., c x = c y = 0.3527 + 0.4854i, cxy = 0.0927  0.2853i. Scatter plots of the estimated the three algorithms are all set to be 400. The mutual coupling coefficient in our experiment is same the literature [19], i.e., cx = cy = 0.3527 + 0.4854i, cxy = 0.0927 − 0.2853i. Scatter plots of the estimated as and the literature [19], i.e., cx = cy =in 0.3527 + 0.4854i, cxy = 0.0927 − 0.2853i. Scatter plots of the estimated azimuth elevation are plotted Figure 11, where the results have been averaged over 500 Monte azimuth and elevation are plotted in Figure 11, where the results have been averaged over 500 Monte azimuth and elevation are plotted in Figure 11, where the results have been averaged over 500 Monte method Carlo runs. A challenge for further research is to incorporate the mutual coupling processing Carlo runs. challenge for further research the mutual coupling processing CarloAruns. A challenge for further researchisisto to incorporate incorporate the mutual coupling processing methodmethod to provide a more realistic scheme like the excellent algorithms [18–20] for applying sensor arrays. to provide a morea realistic scheme likelike thethe excellent [18–20] applying sensor to provide more realistic scheme excellent algorithms algorithms [18–20] for for applying sensor arrays.arrays.

(a)

(a)

(b) Figure 11. Cont.

Figure 11. Cont.

Figure 11. Cont.

(b)

Sensors 2017, 17, 515

13 of 14

Sensors 2017, 17, 515

13 of 14

(c) Figure 11. Scatter plots of estimated 2D angles in the presence of mutual coupling. (a) Classic 2D

Figure 11. Scatter plots(b) of PIE estimated 2D in theAESS presence of mutual coupling. (a) Classic 2D MUSIC algorithm; algorithm; (c)angles The proposed algorithm. MUSIC algorithm; (b) PIE algorithm; (c) The proposed AESS algorithm. 7. Conclusions

7. Conclusions We have proposed a new search-based algorithm, called AESS, for 2D DOA estimation. By exploiting the spectrum peak diffusion effect, our method is capable of searching the sources along

Wetheir have proposed a new search-based algorithm, called AESS, for 2D DOA estimation. directions. Unlike the well-known 2D MISIC, which fixes one of the 2D angles and searches the By exploiting theadjacent spectrum peak diffusion our method is capable of searching thetosources other, the search points of AESSeffect, have different azimuth and elevation parameters along their directions. Unlike the well-known 2D MISIC, whichreduce fixes one of the 2D angles and searches minimize the search path, which enables AESS to drastically the computation and storage costs. thatof AESS achieves RMSE andand success rate butparameters with a the other, theSimulation adjacent results searchshow points AESS have comparable different azimuth elevation to significant computational advantage over 2D MUSIC. Furthermore, the search process is much faster minimize the search path, which enables AESS to drastically reduce the computation and storage than that of 2D MUSIC in the scenarios of different stepsizes. We measured the RMSE versus SNR costs. Simulation results show that AESS achieves comparable RMSE and success rate but with a and snapshot, presented the CRB of each source for performance comparison with the theoretical significant computational advantage over 2Destimation MUSIC. performance, Furthermore, searchthe process much optimum. To further evaluate the 2D DOA wethe presented scatter is plots of faster than that of 2D MUSIC ininthe scenarios different stepsizes. We measured the variance RMSE versus SNR and estimated 2D angles different SNR,of which show that AESS has smaller estimated although it ispresented a biased estimation. simulation demonstrated thatcomparison this bias can with be decreased by slightly snapshot, the CRB ofOur each source for performance the theoretical optimum. enlarging the search ofestimation the third round. To evaluate the wethe examined capability To further evaluate the 2Dradius DOA performance, weresolution, presented scatterthe plots of estimated of the proposed AESS and the related algorithms to estimate the DOAs of closely spaced signals. In 2D angles in different SNR, which show that AESS has smaller estimated variance although it is a the end, we tested AESS and the related algorithms in the presence of mutual coupling. A challenge biased estimation. Our simulation demonstrated that this bias can be decreased by slightly enlarging for further research is to incorporate the mutual coupling processing method to provide a more the search radius of the evaluateinthe resolution, weuse. examined the capability of the realistic scheme suchthird as theround. excellentTo algorithms [18–20] for practical proposed AESS and the related algorithms to estimate the DOAs of closely spaced signals. In the end, Acknowledgments: This work was supported by the National Natural Science Foundation of China under Grant we tested andNo. the relatedand algorithms in Scholarship the presence of mutual coupling. A challenge for further No.AESS 61373177, 61472320, Fund of China Council No. 201606975002. researchAuthor is to incorporate the mutual coupling processing method to provide a more realistic scheme Contributions: The idea was proposed by Dazheng Feng; Kaijie Xu and Weike Nie simulated the such as algorithm; the excellent algorithms in [18–20] for practical use. Aiqin Hou analyzed the data and Chase Qishi Wu designed the experiments and polish the English; Weike Nie and Aiqin Hou wrote the paper. Thanks are given to Xiaoyan Yin for the suggested corrections.

Acknowledgments: This work was supported by the National Natural Science Foundation of China under Grant Conflicts of Interest: The authors declare no conflict of interest. No. 61373177, No. 61472320, No. 61501372 and Fund of China Scholarship Council No. 201606975002. References Author Contributions: The idea was proposed by Dazheng Feng; Kaijie Xu and Weike Nie simulated the algorithm; Aiqin Hou analyzed the data and Chase Qishi Wu designed the experiments and polish the English; 1. Schmidt, R.O. Multiple emitter location and signal parameter estimation. IEEE Trans. Antennas Propag. 1986, Weike Nie and Aiqin Hou wrote the paper. Thanks are given to Xiaoyan Yin for the suggested corrections. 34, 276–280.

Conflicts2.of Interest: The authors declare noalgorithm conflict DOA of interest. Basikolo, T.; Arai, H. APRD-MUSIC estimation for reactance based uniform circular array. IEEE Trans. Antennas Propag. 2016, 64, 4415–4422. 3. Huang, W.; Yang, X.; Liu, D.; Chen, S. Diffusion LMS with component-wise variable step-size over sensor References networks. IET Signal Process. 2016, 10, 37–45.

1. 2. 3. 4.

Schmidt, Multiple emitter location and signal parameter estimation. Trans.estimation. AntennasSignal Propag. 1986, 4. Ni,R.O. J.; Yang, J. Variable step-size diffusion least mean fourth algorithm for IEEE distributed 34, 276–280. [CrossRef] Process. 2015, 122, 93–97. Basikolo, T.; Arai, H. APRD-MUSIC algorithm DOA estimation for reactance based uniform circular array. IEEE Trans. Antennas Propag. 2016, 64, 4415–4422. [CrossRef] Huang, W.; Yang, X.; Liu, D.; Chen, S. Diffusion LMS with component-wise variable step-size over sensor networks. IET Signal Process. 2016, 10, 37–45. [CrossRef] Ni, J.; Yang, J. Variable step-size diffusion least mean fourth algorithm for distributed estimation. Signal Process. 2015, 122, 93–97. [CrossRef]

Sensors 2017, 17, 515

5.

6.

7.

8. 9. 10. 11. 12.

13.

14. 15. 16. 17.

18.

19. 20. 21. 22. 23. 24. 25.

14 of 14

Zhou, C.; Shi, Z.; Gu, Y.; Shen, X. DECOM: DOA estimation with combined MUSIC for coprime array. In Proceedings of the IEEE International Conference on Wireless Communications and Signal Processing, Hangzhou, China, 24–26 October 2013; pp. 1–5. Wang, Y.Y.; Lee, L.C.; Yang, S.J.; Chen, J.T. A tree structure one-dimensional based algorithm for estimating the two-dimensional direction of arrivals and its performance analysis. IEEE Trans. Antennas Propag. 2008, 56, 178–188. [CrossRef] Karthikeyan, B.R.; Kadambi, G.R.; Vershinin, Y.A. A formulation of 1-D search technique for 2-D DOA estimation using orthogonally polarized components of linear array. IEEE Antennas Wirel. Propag. Lett. 2015, 14, 1117–1120. [CrossRef] Kikuchi, S.; Tsuji, H.; Sano, A. Pair-matching method for estimating 2-D angle of arrival with a cross-correlation matrix. IEEE Antennas Wirel. Propag. Lett. 2006, 5, 35–40. [CrossRef] Gershman, A.B.; Rübsamen, M.; Pesavento, M. One- and two dimensional direction-of-arrival estimation: An overview of search-free techniques. Signal Process. 2010, 90, 1338–1349. [CrossRef] Han, F.M.; Zhang, X.D. An ESPRIT-like algorithm for coherent DOA estimation. IEEE Antennas Wirel. Propag. Lett. 2005, 4, 443–446. [CrossRef] Chen, F.J.; Kwong, S.; Kok, C.W. ESPRIT-Like Two-Dimensional DOA Estimation for Coherent Signals. IEEE Trans. Aerosp. Electron. Syst. 2010, 46, 1477–1484. [CrossRef] Chen, F.J.; Fung, C.C.; Kok, C.W.; Kwong, S. Estimation of 2-Dimensional Frequencies Using Modified Matrix Pencil Method. In Proceedings of the IEEE 6th Workshop on Signal Processing Advances in Wireless Communications, New York, NY, USA, 5–8 June 2005; pp. 718–724. Barabell, A.J. Improving the resolution performance of eigenstructure-based direction-finding algorithms. In Proceedings of the International Conference on Acoustics, Speech and Signal Processing, Boston, MA, USA, 14–16 April 1983; pp. 336–339. Qian, C.; Huang, L.; So, H.C. Improved unitary root-MUSIC for DOA estimation based on pseudo-noise resampling. IEEE Signal Process. Lett. 2014, 21, 140–144. [CrossRef] Costa, M.; Richter, A.; Koivunen, V. DoA and polarization estimation for arbitrary array configurations. IEEE Trans. Signal Process. 2012, 60, 2330–2343. [CrossRef] Weiss, A.; Friedlander, B. Direction finding for diversely polarized signal using polynomial rooting. IEEE Trans. Signal Process. 1993, 41, 1893–1905. [CrossRef] Wong, K.; Li, L.; Zoltowski, M. Root-MUSIC-based direction finding and polarization estimation using diversely polarized possibly collocated antennas. IEEE Antennas Wirel. Propag. Lett. 2004, 3, 129–132. [CrossRef] Goossens, R.; Rogier, H. A hybrid UCA-RARE/Root-MUSIC approach for 2-D direction of arrival estimation in uniform circular array in the presence of mutual coupling. IEEE Trans. Antennas Propag. 2007, 55, 841–849. [CrossRef] Ye, Z.F.; Liu, C. 2-D DOA estimation in the presence of mutual coupling. IEEE Trans. Antennas Propag. 2008, 56, 3150–3158. [CrossRef] Filik, T.; Tuncer, T.E. A fast and automatically paired 2-D direction-of-arrival estimation with and without estimating the mutual coupling coefficients. Radio Sci. 2010, 45, 428–451. [CrossRef] Ren, S.W.; Ma, X.C.; Yan, S.F.; Hao, C.P. 2-D unitary ESPRIT-Like direction-of-arrival (DOA) estimation for coherent signals with a uniform rectangular array. Sensors 2013, 13, 4272–4288. [CrossRef] [PubMed] Haardt, M.; Nossek, J.A. Unitary ESPRIT: How to obtain increased estimation accuracy with a reduced computational burden. IEEE Trans. Signal Process. 1995, 43, 1232–1242. [CrossRef] Zoltowski, M.D.; Haardt, M.; Mathews, C.P. Closed-form 2-D angle estimation with rectangular arrays in element space or beamspace via unitary ESPRIT. IEEE Trans. Signal Process. 1996, 44, 316–328. [CrossRef] Perez-Neira, A.I.; Lagunas, M.; Kirlin, R.L. Cross-coupled DOA trackers. IEEE Trans. Signal Process. 1997, 45, 2560–2565. [CrossRef] Zhou, Y.; Feng, D.Z.; Liu, J.Q. A novel algorithm for two-dimension frequency estimation. Signal Process. 2007, 87, 1–12. [CrossRef] © 2017 by the authors. Licensee MDPI, Basel, Switzerland. This article is an open access article distributed under the terms and conditions of the Creative Commons Attribution (CC BY) license (http://creativecommons.org/licenses/by/4.0/).