Application of the Levenberg-Marquardt Scheme to the MUSIC Algorithm for AOA Estimation

Abstract

We note that the cost function of the MUSIC (multiple signal classification) algorithm is quadratic in an array vector and that it can be expressed in a least squares form. Based on this observation, we present a rigorous Levenberg-Marquardt (LM) formulation of the MUSIC algorithm for simultaneous estimation of an azimuth and an elevation. We show a convergence property and compare the performance of the LM-based MUSIC algorithm with that of the standard MUSIC algorithm via Monte-Carlo simulation. We also compare the performance of the MUSIC algorithm with that of the Capon algorithm both for the standard implementation and for the LM-based implementation.

1. Introduction

There has been much study on the determination of the AOA (angle-of-arrival) of incident signals using array antenna structure [1–15]. Array signal processing algorithm can also be extended to radar target recognition [16, 17] and detection [18, 19].

In this paper, we consider a simultaneous estimation of the azimuth and the elevation using the MUSIC algorithm [2]. In standard MUSIC algorithm, we have to perform two-dimensional exhaustive grid search, which can be computationally intensive when the search step is small. The idea proposed in this paper is to apply the LM algorithm to the numerical optimization of the cost function defined from the MUSIC algorithm. From the numerical results, it is shown that a performance of the proposed scheme is better than that of the standard MUSIC algorithm in terms of the computational burden and the estimation accuracy.

In this paper, we give a rigorous formulation for the improvement of the azimuth/elevation AOA estimation for the UCA structure by applying the LM scheme [20, 21] to the cost function of the MUSIC algorithm.

Assume there are 𝑀 incident signals. The cost function of the maximum likelihood (ML) algorithm [1, 3–7] is 𝑀-dimensional for estimation of azimuth angle only and 2𝑀-dimensional for simultaneous estimation of azimuth angle and elevation angle. For large value of 𝑀, the optimization of 𝑀-dimensional or 2𝑀-dimensional cost function by exhaustive grid search is not feasible. That is why Newton-type algorithms including the LM algorithm have been widely applied to the optimization of the cost function of the ML algorithm.

Note that, in the MUSIC algorithm [1, 2], the cost function is one-dimensional for estimation of azimuth angle and two-dimensional for simultaneous estimation of azimuth angle and elevation angle, which is independent of the actual number of the incident signals 𝑀. When we are only concerned with the estimation of the azimuth angle, the optimization of one-dimensional cost function using exhaustive grid search even with small search step is quite feasible. But, for simultaneous estimation of the azimuth and the elevation, the computational complexity of two-dimensional optimization by exhaustive grid search is much larger than that of one-dimensional optimization by exhaustive grid search, which justifies why we apply the LM algorithm to the optimization of the two-dimensional cost function of the MUSIC algorithm.

There has been much study on applying the Newton-type method to the optimization of the cost function of the ML DOA algorithm [3–7]. To the best of our knowledge, there has been no study on applying the LM algorithm to the optimization of the cost function of the MUSIC algorithm for simultaneous estimation of the azimuth angle and the elevation angle.

2. Standard MUSIC Algorithm

Consider the case of a uniform circular array (UCA) with 𝑁 antenna elements. 𝑀 signals are incident on the array, and we are to estimate the azimuth angle and the elevation angle of each incident signal. Note that, for simultaneous estimation of the azimuth and elevation, a uniform linear array (ULA) cannot be adopted since ULA-based algorithm cannot uniquely estimate the azimuth and the elevation simultaneously due to the ambiguity pertinent to the ULA structure.

Let 𝜃 and 𝜙 represent the azimuth and the elevation, respectively. We denote the noise eigenvector matrix of the array covariance matrix 𝐔. Then the MUSIC spectrum can be usually defined as [4] 𝑃=1𝜃,𝜙𝐚^H𝜃,𝜙𝐔𝐔^H𝐚.𝜃,𝜙 (1)

The array vector can be written as 𝐚𝜃,𝜙=exp𝑗𝜓₁exp𝑗𝜓₂⋯exp(𝑗𝜓_𝑁)^𝑇, (2) where 𝜓_𝑛 for the UCA is defined as follows: 𝜓_𝑛𝑟𝜃,𝜙=2𝜋𝜆()cos𝜙cos𝜃−2𝜋𝑛−1𝑁𝑛=1,…,𝑁. (3)

In the standard MUSIC algorithm, to evaluate (1), two-dimensional exhaustive grid search with search step of [Δ𝜃,Δ𝜙] has to be implemented.

3. LM Implementation of the MUSIC Algorithm

Let 𝐹(𝜃,𝜙) be defined from the MUSIC spectrum in (1) as 𝐹=1𝜃,𝜙𝑃𝜃,𝜙=𝐚^H𝜃,𝜙𝐔𝐔^H𝐚𝜃,𝜙. (4)

It is quite clear that the local maxima of 𝑃(𝜃,𝜙) correspond to the local minima of 𝐹(𝜃,𝜙).

Let a vector valued function 𝐟(𝜃,𝜙) be in the form of 𝐟𝜃,𝜙=𝑓₁𝜃,𝜙𝑓₂𝜃,𝜙⋯𝑓_{2(𝑁−𝑀)}𝜃,𝜙^𝑇. (5)

To apply the LM algorithm to the optimization of the cost function of the MUSIC algorithm in (4), what we have to do first is to define 𝐟(𝜃,𝜙) so that ‖𝐟(𝜃,𝜙)‖²₂ corresponds to the cost function of 𝐹(𝜃,𝜙) in (4): 𝐹𝜃,𝜙=𝐟(𝜃,𝜙)²₂=^{2(𝑁−𝑀)}_𝑚=1𝑓_𝑚(𝜃,𝜙)². (6)

The MUSIC algorithm can be formulated as 𝜃,𝜙=argmin_𝜃,𝜙𝐹𝜃,𝜙=argmin_𝜃,𝜙𝐟𝜃,𝜙²₂=argmin_𝜃,𝜙^{2(𝑁−𝑀)}_𝑚=1𝑓_𝑚𝜃,𝜙². (7)

We will show how to express the cost function in (4) in least squares form. From (4), we have 𝐹𝜃,𝜙=𝐚^H𝜃,𝜙𝐔𝐔^H𝐚𝜃,𝜙=𝐔^H𝐚𝜃,𝜙². (8)

In (8), note that 𝐔^H𝐚(𝜃,𝜙) is a column vector with 𝑁−𝑀 entries. Let the 𝑖th entry of 𝐔^H𝐚(𝜃,𝜙) be denoted by 𝑒_𝑖(𝜃,𝜙), 𝑖=1,2,…,𝑁−𝑀: 𝑒_𝑖𝜃,𝜙≡𝐔^H𝐚𝜃,𝜙_𝑖=^𝑁_𝑛=1Re𝐮_𝑖𝑛−𝑗Im𝐮_𝑖𝑛×cos𝜓_𝑛𝜃,𝜙+𝑗sin𝜓_𝑛𝜃,𝜙. (9)

Note that, since 𝐚(𝜃,𝜙) is complex-valued, 𝑒_𝑖(𝜃,𝜙) is also complex valued. Let 𝐞(𝜃,𝜙) be defined as follows: 𝐞𝜃,𝜙=𝑒₁𝜃,𝜙⋯𝑒_𝑁−𝑀𝜃,𝜙^𝑇=𝐮^H₁𝐚𝜃,𝜙⋯𝐮^H_𝑁−𝑀𝐚𝜃,𝜙^𝑇. (10)

Using (8) and (10), it can be easily shown that the cost function of the MUSIC algorithm can be written as 𝐹𝜃,𝜙=𝐞(𝜃,𝜙)²₂=^𝑁−𝑀_𝑖=1𝑒_𝑖(𝜃,𝜙)². (11)

Let 𝑓_2𝑖−1(𝜃,𝜙) and 𝑓_2𝑖(𝜃,𝜙) in (5) be defined by the real part and the imaginary part of 𝑒_𝑖(𝜃,𝜙), respectively. From (9), we have, for 𝑖=1,2,…,𝑁−𝑀, 𝑓_2𝑖−1𝜃,𝜙≡Re𝑒_𝑖𝜃,𝜙=Re𝐮^H_𝑖𝐚=𝜃,𝜙_𝑁𝑛=1Re𝐮_𝑖𝑛cos𝜓_𝑛𝐮𝜃,𝜙+Im_𝑖𝑛sin𝜓_𝑛,𝑓𝜃,𝜙_2𝑖𝜃,𝜙≡Im𝑒_𝑖𝜃,𝜙=Im𝐮^H_𝑖𝐚=𝜃,𝜙_𝑁𝑛=1Re𝐮_𝑖𝑛sin𝜓_𝑛𝐮𝜃,𝜙−Im_𝑖𝑛cos𝜓_𝑛𝜃,𝜙. (12)

From the definition of 𝑓_2𝑖−1(𝜃,𝜙) and 𝑓_2𝑖(𝜃,𝜙), it is clear that the following is true: 𝑒_𝑖𝜃,𝜙²=𝑓_2𝑖−1𝜃,𝜙²+𝑓_2𝑖𝜃,𝜙²𝑖=1,…,𝑁−𝑀. (13)

From (11) and (13), the cost function can also be written in terms of 𝑓_𝑚(𝜃,𝜙)^′s as follows: 𝐹𝜃,𝜙=^𝑁−𝑀_𝑖=1𝑒_𝑖𝜃,𝜙²=^{2(𝑁−𝑀)}_𝑖=1𝑓_𝑖𝜃,𝜙². (14)

The Jacobian matrix corresponding to the matrix 𝐟(𝜃,𝜙) in (5) is defined as 𝐉𝜃,𝜙=𝜕𝑓₁𝜕𝜃𝜕𝑓₁𝜕𝜙⋮⋮𝜕𝑓_{2(𝑁−𝑀)}𝜕𝜃𝜕𝑓_{2(𝑁−𝑀)}𝜕𝜙. (15)

Using (12) in (15), the explicit expressions for the entries of the matrix 𝐉(𝜃,𝜙) are, for 𝑖=1,2,…,𝑁−𝑀, 𝜕𝑓_2𝑖−1=𝜕𝜃^𝑁_𝑛=1𝜓_𝑛𝜃,𝜙𝐮𝜕𝜃−Re_𝑖𝑛sin𝜓_𝑛𝐮𝜃,𝜙+Im_𝑖𝑛cos𝜓_𝑛𝜃,𝜙,𝜕𝑓_2𝑖−1=𝜕𝜙^𝑁_𝑛=1𝜓_𝑛𝜃,𝜙𝐮𝜕𝜙−Re_𝑖𝑛sin𝜓_𝑛𝐮𝜃,𝜙+Im_𝑖𝑛cos𝜓_𝑛𝜃,𝜙,𝜕𝑓_2𝑖=𝜕𝜃^𝑁_𝑛=1𝜓_𝑛𝜃,𝜙𝐮𝜕𝜃Re_𝑖𝑛cos𝜓_𝑛𝐮𝜃,𝜙+Im_𝑖𝑛sin𝜓_𝑛𝜃,𝜙,𝜕𝑓_2𝑖=𝜕𝜙^𝑁_𝑛=1𝜓_𝑛𝜃,𝜙𝐮𝜕𝜙Re_𝑖𝑛cos𝜓_𝑛𝐮𝜃,𝜙+Im_𝑖𝑛sin𝜓_𝑛𝜃,𝜙. (16)

In Figure 1, we outline the Levenberg-Marquardt algorithm. In the problem at hand, 𝐱 can be written as 𝐱=(𝜃,𝜙), and 𝑘 is the number of the iterations [20], and 𝐹(𝜃,𝜙) is defined in (8) or (14).

Figure 1 

Levenberg-Marquardt algorithm [20].

In the LM algorithm, the updates of the estimates are obtained from 𝐉𝜃,𝜙^𝑇𝐉𝜃,𝜙+𝜇𝐈𝜃_LM𝜙_LM^𝑇=−𝐉𝜃,𝜙^𝑇𝐟𝜃,𝜙. (17)

Due to 𝜇>0, the coefficient matrix is positive definite, which makes 𝜃_LM𝜙_LM^𝑇 be a descent direction. When 𝜇 is very large, we get 𝜃_LM𝜙_LM^𝑇1≃−𝜇𝐠 (18) which is a short step in the steepest descent direction. Therefore, we choose large 𝜇 when the current estimates are assumed to be very far from the solution. This is good if the current iterate is far from the solution. If 𝜇 is very small, then the updates of the LM method are approximately equal to the updates of the Newton method, which is a good step in the final stages of the iteration [20].

𝜏 can be chosen to be small value of 𝜏=10⁻⁶ if the initial estimate, [𝜃₀,𝜙₀], is assumed to be close to the solution. Otherwise, 𝜏 is chosen to be larger value of 𝜏=10⁻³ [20].

𝐋(𝜃_LM,𝜙_LM) in Figure 1 is defined as [20, 21] 𝐋𝜃_LM,𝜙_LM=𝐹𝜃,𝜙+𝜃_LM𝜙_LM𝐉𝜃,𝜙^𝑇𝐟+1𝜃,𝜙2𝜃_LM𝜙_LM𝐉𝜃,𝜙^𝑇𝐉𝜃,𝜙𝜃_LM𝜙_LM^𝑇. (19)

4. Comment on How to Get Initial Estimates

The alternating projection (AP) algorithm [22, 23] is very popular algorithm for getting the initial estimates of the AOA estimation since it is computationally efficient. The shortcoming of the AP scheme is that the initial estimates obtained from the AP scheme are less accurate in comparison with those obtained from the exhaustive search algorithm. Recently, there is a study [15] on the extension of the original AP scheme provided in [22, 23].

5. Numerical Results

In numerical results, the UCA is used. Therefore, 𝐟(𝜃,𝜙) is obtained from (12), and 𝐉(𝜃,𝜙) is given by (16) by substituting (3) for 𝜓_𝑛(𝜃,𝜙).

The parameters for the LM algorithm are defined as follows: both 𝜖₁ and 𝜖₂ are chosen to be 10⁻⁵. 𝜏 is set to be 1.0. The maximum number of the iteration which corresponds to 𝑘_max is 100. The RMSE (root mean square error) is obtained from ten repetitions.

In the first example, the LM algorithm is initialized on a grid of azimuth-elevation to check which initial angles finally achieve the global convergence. Therefore, initial estimates used for the LM algorithm are chosen arbitrarily to see whether each of the initial estimates actually leads to the global minimum or not. The cost function is evaluated at each two-dimensional grid point with azimuth increment of Δ𝜃=4^∘ and the elevation increment of Δ𝜙=4^∘ followed by the LM optimization.

Since the scheme is based on the LM formulation, it is shown that once the initial estimates are within the main lobe, the initial estimates are supposed to converge to the local minimum corresponding to the main lobe which is actually a global minimum.

The initial azimuth angles and the initial elevation angles which result in the global convergence are marked in Figure 2. The criterion for the global convergence is that the RMSE is arbitrarily chosen to be less than two degrees.

(a) 𝑟=0.2𝜆

(b) 𝑟=0.4𝜆

(c) 𝑟=0.6𝜆

(d) 𝑟=0.8𝜆

Figure 2 

Initial azimuth angles and the initial elevation angles resulting in the global convergence for various radii at SNR = 4 dB.

For the fixed number of the antenna elements, as the radius of the antenna array becomes larger, the number of the spurious peaks of 𝑃(𝜃,𝜙) increases. The initial estimate close to the spurious peak of 𝑃(𝜃,𝜙) will converge to the spurious peak, which is not the global minimum of 𝐹(𝜃,𝜙). Therefore, the larger radius of the antenna array results in smaller domain of initial angles leading to global convergence as shown in Figure 2, which reflects that the proposed scheme is very efficient in the case of the smaller radius of the antenna array.

In Figure 3, for two incident signals, we compared the performance of the LM optimization following the MUSIC algorithm of coarse exhaustive grid search with that of the standard MUSIC algorithm of fine exhaustive grid search. The purpose is to illustrate the fact that the proposed scheme still works for more than one incident signal.

(a) [Δ𝜃fine,Δ𝜙fine]=[2∘,2∘], [Δ𝜃coarse,Δ𝜙coarse]=[4∘,4∘]

(b) [Δ𝜃fine,Δ𝜙fine]=[2∘,2∘], [Δ𝜃coarse,Δ𝜙coarse]=[7∘,7∘]

(c) [Δ𝜃fine,Δ𝜙fine]=[4∘,4∘], [Δ𝜃coarse,Δ𝜙coarse]=[7∘,7∘]

Figure 3 

Comparison of the standard MUSIC algorithm with fine exhaustive grid search and the MUSIC algorithm with coarse exhaustive grid search followed by LM optimization (𝑟=0.679𝜆, five antennas, [𝜃(true)1𝜙(true)1]=[−21∘30∘], [𝜃(true)2𝜙(true)2]=[120.5∘72.5∘]).

[Δ𝜃^coarse,Δ𝜙^coarse] and [Δ𝜃^fine,Δ𝜙^fine] denote the search step of the coarse exhaustive grid search and the fine exhaustive grid search, respectively. We arbitrarily choose Δ𝜃^coarse=Δ𝜙^coarse and Δ𝜃^fine=Δ𝜙^fine. We can see that the LM implementation of the MUSIC algorithm is superior to the standard MUSIC algorithm for all three cases in terms of the RMSE and the computational complexity.

In Figure 4, we compare the performances of the MUSIC algorithm and the Capon algorithm [11]. Each algorithm is implemented in two ways: the first way is the standard MUSIC algorithm and the standard Capon algorithm of fine exhaustive grid search. The second way is the LM optimization following the MUSIC and Capon algorithms of coarse exhaustive grid search.

(a) [𝜃(true)1𝜙(true)1]=[−21∘30∘], [𝜃(true)2𝜙(true)2]=[120.0∘72.0∘]

(b) [𝜃(true)1𝜙(true)1]=[85.5∘63.5∘], [𝜃(true)2𝜙(true)2]=[120.5∘72.5∘]

Figure 4 

Comparison of the standard MUSIC algorithm and the standard Capon algorithm with fine exhaustive grid search and the MUSIC algorithm and the Capon algorithm with coarse exhaustive grid search followed by LM optimization (𝑟=0.679𝜆, five antennas, [Δ𝜃fine,Δ𝜙fine]=[4∘,4∘], [Δ𝜃coarse,Δ𝜙coarse]=[7∘,7∘]).

In Figure 4(a), it is shown that the performances of the MUSIC algorithm and the Capon algorithm are nearly equal for the case that the incident signals are separated enough. In Figure 4(b), for the case that two incident signals are close to each other, it is shown that the performance of the MUSIC algorithm is better than that of the Capon algorithm, both for the standard implementation and for the LM-based implementation. The results in Figure 4(b) can be justified due to the fact that the MUSIC algorithm is a superresolution algorithm, while the Capon algorithm is not a superresolution algorithm.

In the standard implementation of the Capon algorithm, when two incident signals are close to each other in Figure 4(b), the Capon algorithm fails to resolves two closely incident signals at low SNRs of 4 dB and 8 dB, which is why the RMSEs of the first signal at SNRs of 4 dB and 8 dB are very large. On the other hand, the corresponding RMSEs for the MUSIC algorithm are much smaller than those for the Capon algorithm, which reflects that the MUSIC algorithm successfully resolve two closely incident signals at low SNRs, while the Capon algorithm cannot resolve. For high SNRs of 12 dB to 40 dB in increment of 4 dB, the RMSEs of the first signal for the Capon algorithm and for the MUSIC algorithm are small, which reflects that both algorithms can resolve two closely incident signals at high SNRs.

6. Conclusion

In this paper, we propose to apply the LM formulation to the optimization of the cost function of the MUSIC algorithm. Using the LM algorithm, we implemented the nonlinear multivariable optimization for azimuth/elevation AOA (angle-of-arrival) estimation based on the MUSIC algorithm.

We have compared the performance of the LM implementation of the MUSIC algorithm with that of the standard MUSIC algorithm in terms of the estimation accuracy and the computational complexity.

We illustrated that the performance of the MUSIC algorithm is better than that of the LM-based implementation as well as the standard implementation when there are two incident signals that are very close to each other.

We consider the array structure of the UCA to be able to uniquely estimate the azimuth and elevation simultaneously, but it is quite straightforward to extend the proposed scheme to an arbitrary array structure by simply modifying the array vector consistently with the specific array structure as long as the adopted array structure is able to uniquely estimate both the azimuth angle and the elevation angle. In addition, the proposed scheme can also be applied to the other AOA estimation algorithm by simply modifying the cost function consistently with the adopted AOA estimation algorithm as long as the corresponding cost function can be reformulated in the least squares form.

Acknowledgment

This research was supported by Basic Science Research Program through the National Research Foundation of Korea (NRF) funded by the Ministry of Education (2013R1A1A2A10012245).