Asymptotic Performance Analysis of Maximum Likelihood Algorithm for Direction-of-Arrival Estimation: Explicit Expression of Estimation Error and Mean Square Error

This paper proposes a new method to get explicit expressions of some quantities associated with performance analysis of the maximum likelihood DOA algorithm in the presence of an additive Gaussian noise on the antenna elements. The motivation of the paper is to make a quantitative analysis of the ML DOA algorithm in the case of multiple incident signals. We present a simple method to derive a closed-form expression of the MSE of the DOA estimate based on the Taylor series expansion. Based on the Taylor series expansion and approximation, we get explicit expressions of the MSEs of estimates of azimuth angles of all incident signals. The validity of the derived expressions is shown by comparing the analytic results with the simulation results.

In [8], a performance analysis of the ML DOA estimation algorithm for low SNR and small number of snapshots is considered. A threshold effect in the ML DOA algorithm is exploited, and the authors derive approximations to the mean square error and probability of outlier.
If the noise variance at each sensor in the array antenna system is equal, the noise covariance matrix is considered to be multiples of an identity matrix. In [9], nonuniform white noise, whose covariance matrix can be expressed as an arbitrary diagonal covariance, is considered, and the new ML DOA algorithm for nonuniform noise is proposed, and the performance analysis of the proposed algorithm is also presented.
In [10], the authors addressed the DOA estimation using sparse sensor arrays, where the sensor noises can be uncorrelated between different subarrays due to large intersubarray spacings. The authors proposed a new maximum-likelihood estimator, which can be extended to the uncalibrated arrays with sensor gain and phase mismatch.
A new computationally efficient ML DOA algorithm exploiting spatial aliasing is proposed in [11]. Generally, spatial aliasing is undesirable since it degrades the DOA estimation accuracy. In [11], the authors adopted a nested array structure with a doubly spaced aperture. The computational burden of the ML DOA estimation algorithm is reduced by the highly compressed search range and the small number of candidate angles to be searched. The authors also presented Monte Carlo simulation based mean square error (MSE). However, analytic performance analysis of the proposed scheme is not presented in [11].
A new ML DOA algorithm for use with a uniform linear array is proposed [12]. The scheme is superior to the conventional ML algorithm when the true DOAs of two incident signals are very close to each other. The formulation for the new DOA algorithm is based on an asymptotic approximation of the unconditional maximum likelihood (UML) procedure when two closely space signals are incident on the ULA. Taylor approximation is also adopted for derivation of the new algorithm. Empirical MSE is illustrated to validate the proposed scheme. However, the authors do not present analytic performance analysis of the new DOA algorithm.
To overcome the problem of large computational complexity for implementation of the ML DOA estimation algorithm, based on a spatially overcomplete array output formulation, an efficient ML DOA estimator is proposed in [13]. Empirical performance from Monte Carlo simulation is present to illustrate the superiority of the proposed scheme over the other DOA algorithms. Analytic performance analysis is not given.
Although the ML estimator is known to be optimal in DOA estimation, its computational cost can be quite prohibitive, especially for a large number of incident signals. To solve this problem, in [14], three kinds of natural computing algorithms, differential evolution, clonal selection algorithm, and the particle swarm optimization, are applied for implementation of the multivariable nonlinear optimization of the cost function of the ML DOA estimation algorithm. It turns out that all three natural computing algorithms are capable of optimizing the ML DOA cost function, irrespective of the number of incident signals and their nature. In addition, the number of points evaluated by natural computing algorithms is much smaller than that associated with exhaustive grid search-based algorithms, justifying the application of these natural computing algorithms to the optimization of the cost function of the ML DOA estimation algorithm.
In [15], a new implementation of ML DOA estimation, which outperforms the other DOA algorithms for closely spaced incident signals, is proposed. The concept of Monte Carlo importance sampling is applied. The superiority of the proposed scheme comes from its better convergence to a global maximum in comparison with other iterative approaches. Although analytical performance analysis of the proposed scheme is not presented, empirical performance of the propose algorithm and the other DOA algorithms is given. Note that Monte Carlo simulation is employed to get empirical performance in terms of the MSE of the DOA estimate.
A heuristic optimization algorithm, called gravitational search algorithm, is presented to optimize the cost function of the ML DOA estimation algorithm for a uniform circular array [16]. It is empirically shown that the proposed algorithm is superior to the MUSIC algorithm and particle swarm optimization-based ML algorithm. Analytic performance analysis of the proposed scheme is not presented in [16].
To reduce computational burden of optimizing the cost function of ML DOA estimation algorithm, the artificial bee colony (ABC) algorithm is applied to maximize the cost function of the ML DOA estimation algorithm [17]. It is empirically shown that the proposed scheme is superior to other ML-based DOA estimation methods in the view point of efficiency in computation and statistical performance. Analytic performance analysis of the proposed scheme is not presented in [17].
DOA estimation of narrowband sources in unknown nonuniform white noise is considered in [18]. The stepwise concentration of the log-likelihood function with respect to the signal parameters and noise parameters is obtained by alternating minimization of the Kullback-Leibler divergence. Closed-form expressions for the signal parameters and noise parameters are derived, implying that the proposed scheme results in significant reduction in computational complexity in comparison with exhaustive multidimensional search-based ML DOA algorithms.
In [19], a new wideband ML DOA estimation algorithm for an unknown nonuniform sensor noise is proposed to reduce the performance degradation due to nonuniformity of the noise. Two associated implementation schemes are proposed: one is iterative and the other is non-iterative. Simulation results show that the performance of two processing algorithms is consistent with the Cramer-Rao lower bound. Analytic performance analysis, more specifically the Cramer-Rao lower bound, of the proposed algorithm is presented in [19].
In this paper, we are concerned with quantitative study on how much estimation error is induced due to an additive Gaussian noise on array antennas. More specifically, mean-squared error (MSE) of direction-of-arrival estimation in terms of a standard deviation of an additive noise is derived. In this paper, performance analysis of azimuth estimation using uniform linear array (ULA) is presented.
In this paper, the estimate with no superscript denotes the estimate of the original ML algorithm. Note that no approximation is used in getting the estimate with no superscript. In this paper, the estimate with the superscript (u = 1) denotes the estimate by using the first approximation, and that with the superscript (u = 1, v) represents the estimate by using the first approximation and the second approximation.
The difference between the estimate with no superscript and the estimate with the superscript (u = 1) quantifies the error due to the first approximation since the first approximation is applied in getting the estimate with the superscript (u = 1). Note that no approximation is applied in getting the estimate with no superscript. Similarly, the difference between the estimate with the superscript (u = 1) and the estimate with the superscript (u = 1, v) quantifies the error due to the second approximation since the first approximation and the second approximations are applied in getting the estimate with the superscript (u = 1, v). Based on this intuition, by comparing these three estimates, we can easily determine which approximation results in the dominant approximation error. This inspection cannot be obtained from the scheme presented in the previous study [7][8][9]19].
In this paper, Gaussian noise is used to model measurement uncertainty. The effect of Gaussian noise on the accuracy of the azimuth estimate is rigorously derived. Furthermore, an explicit expression of the MSE of the azimuth estimate is also derived. In comparison with the previous studies on the performance analysis of the maximum likelihood algorithm [7][8][9]19], a more explicit representation of the MSE of the azimuth estimate is proposed in this paper.
Many previous studies on the ML DOA estimation algorithm focused on how the performance of the ML DOA estimation algorithm can be improved by proposing new algorithms or by modifying the ML DOA estimation algorithm [9][10][11][12][13][14][15][16][17][18][19]. Note that our contribution in this paper does not lie in how much improvement can be achieved by proposing an improved ML DOA algorithm. Our contribution lies in a reduction in computational cost in getting the MSE of an existing ML DOA algorithm by adopting an analytic approach, rather than the Monte Carlo simulation-based MSE under measurement uncertainty which is assumed to be Gaussian distributed. That is, the scheme described how analytic MSE can be obtained with much less computational complexity than the Monte Carlo simulation-based MSE.
In this paper, the derivation is based on the Taylor series expansion of the sample covariance matrix since the cost function of the ML DOA estimation algorithm can be explicitly written in terms of the sample covariance matrix. The difference between the sample covariance matrix associated with noisy measurement and that associated with noiseless sample covariance matrix is explicitly expressed in terms of additive noises on the antenna arrays. Azimuth estimation error is explicitly expressed in terms of the additive noises. Finally, the MSE of the azimuth estimate is given in terms of the statistics of an additive noise. To the best of our knowledge, no previous study presented these explicit expressions of the azimuth estimation error and the MSE of the azimuth estimate in terms of the statistics of an additive noise.
The proposed scheme can be used in predicting how accurate the estimate of the ML DOA estimation algorithm is without a computationally intensive Monte Carlo simulation. The performance of the ML algorithm depends on various parameters including the number snapshots, the number antenna elements in the array, inter-element spacing between adjacent antenna elements, and the SNR. Therefore, Monte Carlo simulations for different values of the various parameters can be computationally intensive. Therefore, the scheme presented in this paper can be adopted to predict the accuracy of the ML DOA algorithm for different values of various parameters.

Maximum Likelihood Algorithm
In this section, the maximum likelihood (ML) for use with the uniform linear array (ULA) algorithm is briefly described.
In the case of ULA, for the incident signal from θ c , the array vector associated with the m-th antenna can be written as where λ is wavelength and ∆ is the distance between two neighboring elements. Using (1), A (θ 1 , θ 2 , · · · , θ d ) is defined as where the number of incident signals is d.
Projection matrix on to the column space of A (θ 1 , θ 2 , · · · , θ d ) can be expressed as [6] The incident signals on the array antenna elements can be written as, for i = 1, · · · , L, The noisy incident signals can be expressed as It is assumed that the entries of the Gaussian vector are independent and identically distributed Gaussian random variables with the same mean and the same variance. Note that the noise are complex-valued and that the real part and the imaginary part of the noise are independent Gaussian random variables with non-zero mean µ. The variance of the real part is denoted by σ 2 2 , which is equal to the variance of the imaginary part, σ 2 2 . Let L denote the number of the snapshots.R is a sample covariance matrix given bŷ where x (t i ) is given by (4). From (5), we getR δR is defined as δR =R −R.
From (11), we obtain where, from (1), B(θ k , θ l ) is defined as Using (12), can be written as adj (·) denotes adjoint of a matrix. k, l-th element of the adjoint of A (θ 1 , · · · , θ d ) H A (θ 1 , · · · , θ d ) can be expressed as where the determinant can be obtained in many ways, one of which is a cofactor expansion. Generally, the number of incident signals is d. From (14) and (15), an explicit expression of the entry at the k-th row and l-th column of can be obtained. However, it is very complicated to express the determinant in (15) in terms of the entries of A (θ 1 , · · · , θ d ) H A (θ 1 , · · · , θ d ) for all k = 1, . . . , d and l = 1, . . . , d. In addition, due to very complicated expression, expressing each entry of A (θ 1 , · · · , θ d ) H A (θ 1 , · · · , θ d ) −1 in terms of the entries of A (θ 1 , · · · , θ d ) H A (θ 1 , · · · , θ d ) may impair the readability of this paper. Therefore, the number of signals in this paper is set to two.
Using (12) in (3), the entries of P A(θ 1 ,θ 2 ) can be written as where det (θ 1 , θ 2 ) is given in Appendix B. The numerator of P kl (θ 1 , θ 2 ) is defined as Q kl (θ 1 , θ 2 ): Let θ (0) 1 and θ (0) 2 denote the true incident angles of two incident signals and letθ 1 andθ 2 denote the estimates of two incident signals. From (9) and (10), f kl (θ 1 , θ2) and g kl (θ 1 , θ2) in Appendix I, the estimates, θ 1 ,θ 2 , satisfy ∂ ∂θ 1 tr Using (26) and (27) in (25), (25) can be written as The solution of (28) and the associated estimates are given by where the superscript (u = 1) indicates that the first order Taylor expansion is used. This approximation is called U approximation. At high SNR, it is true thatR lk is much larger than δR lk : The solution of (32) and the associated estimates are given by

Closed-Form Expression of Mean Square Error ofθ
From (34), analytic mean square errors (MSEs) of δθ     The number of snapshots   The number of snapshots    Similarly, the results with 'Simulation E

For all the results in Figures 1-4, the difference between Simulation
, which implies that U approximation results in much greater error than V approximation. Therefore, to improve DOA estimation performance, second-order Taylor expansion, which corresponds to u = 2, can be used.
Actually, in all the results in Figures 1-4 ' are approximately equal.
Instead, to quantify how computationally efficient the proposed algorithm is, execution time is obtained both for analytically derived MSE and for the Monte Carlo simulation-based MSE.
The number of incident signals is two, where two signals are incident from 20 • and 40 • . The number of antenna elements is 10, and the number of snapshots is 1000. In getting the Monte Carlo simulation-based MSE, since the computational complexity is nearly proportional to the number of repetitions, the number of repetitions varies from 100 to 1000 in increments of 100. Figure 5 shows how computationally efficient the proposed algorithm is. The execution times of the simulation-based MSEs and analytic MSEs are illustrated with respect to the number of repetitions. Note that the execution time for analytically derived MSE is essentially independent of the number of repetitions. It is clearly shown in Figure 5 that execution time for the Monte Carlo simulation-based MSE is much greater than that for the analytically derived MSE even for the number of repetitions of 100. Figure 5 illustrates that getting analytically derived MSE is much less computationally intensive than getting Monte Carlo simulation based MSE, which justifies why the analytically derived MSEs should be employed for performance analysis.

Conclusions and Summary
In Figure   To the best of our knowledge, no previous study used explicit equations (A) and (B) to derive the azimuth error in equation (C). One of the novelties of this paper is that the derivation is based on the observation that the azimuth error in the subscript in equation (B) can be analytically derived since equation (A) is true for noiseless covariance matrix.
In summary, applying U approximation to equation (B) and using the constraint in equation (A), the estimation error in equation (C) can be obtained. The MSE of the azimuth estimation error in equation (D) is obtained from equation (C) and V approximation. Note that, to get equation (D) from equation (C), the statistics of an additive noise should also be exploited.
Quantitative study on the estimation error for direction-of-arrival estimation in terms of standard deviation of an additive noise has been addressed in this paper.
In this paper, in case of estimating azimuths of multiple incident signals, the closed-form expression of the MSE of the DOA estimate for the ML algorithm has been derived by stepwise approximations. The type of antenna array is assumed as ULA. The closed-form of the MSE has been derived by using the Taylor series approximation which linearizes the nonlinear parts of the array vector and additional approximation based on the assumption that the estimation error is very small at high SNR. The closed-form of the MSE has been verified through numerical results. The closed-form of the MSE has been verified through numerical results. All of the stepwise approximated simulation results and the results obtained from the closed-form of the MSE show good agreement.
Although the formulation and the numerical results for two incident signals are presented in this paper, an extension to multiple incident signals is quite intuitively clear and straightforward.
In this paper, we rigorously derive how the MSE of the ML algorithm for direction-of-arrival can be expressed in terms of various parameters, which include the number of sensor elements, the number of incident signals, the number of snapshots, and the variance of additive noises on the antenna elements. Although, for convenience, the statistics of additive noise is assumed to be non-zero-mean Gaussian distributed, the derivation in the Appendices can be extended to the case where noise can be modeled as any other random variable as long as the moments of the random variable are analytically available.

Conflicts of Interest:
The authors declare no conflict of interest.

(·) H
Hermitian matrix transpose (·) Noisy quantity δ (·) Difference between noisy quantity and the corresponding noiseless quantity A (θ 1 , θ 2 ) Matrix whose columns are array vector for θ 1 and θ 2 θ c Estimate of the c-th incident signal c = 1, 2 Noisy signals on the antenna array at t − tî R Sample covariance matrix of the noiseless signal R lk The entry at the k-th row and the l-th column ofR R Sample covariance matrix of the noisy signal R lk The entry at the k-th row and the l-th column of theR δR Difference betweenR andR δR lk Difference betweenR lk andR lk (·) (u=1) First, order U approximation of (·) based on Taylor expansion associated with w-th repetition out of W repetitions. c = 1, 2

Abbreviations
The following abbreviations are used in this manuscript: DOA Note that the first and the second order derivatives of det (θ 1 , θ 2 ) are given in Appendix B and first and second order derivatives of B (θ 1 , θ 2 ) are given in Appendix C.

Appendix C. First and Second Order Derivatives of
The first order and the second order derivatives of a k (θ 1 ) and a k (θ 2 ) are given in Appendix D.

Appendix G. Derivation of E δR lk δR
where the second moments, third moments, and the fourth moments in (A56) are derived in Appendices H-J, respectively. Finally, when i = i , E (δR) lk (δR) * l k is given by for l = k and k = l and l = k .
for l = k and k = l and l = k .
for l = k and k = k and l = k .
for l = l and l = k and k = k.
for l = k and l = k and l = k .
for l = k and k = l and l = k .
for l = l and k = l and k = k .
for l = k and k = k and k = l .
for l = k and k = l and l = k and k = l.
for l = l and l = k and k = k and k = l.
for l = k and k = k and k = l and l = l.
for k = l and l = l and l = k and k = k.
for k = k and k = l and l = l and l = k.
for l = k and k = l and l = k and k = l .
otherwise. When i = i , E (δR) lk (δR) * l k is given by for l = k and k = l and l = k .
for l = k and k = l and l = k .
for l = l and l = k and k = k .
for l = k and k = l and l = k. otherwise.

Appendix I. Third Order Non-Central Moment of Non-Zero-Mean Complex Gaussian Random Variable with Variance σ 2
We define ten cases depending on how a, b, c, d, and e are related: