# Efficient AM Algorithms for Stochastic ML Estimation of DOA

Authors: Haihua Chen; Shibao Li; Jianhang Liu; Yiqing Zhou; Masakiyo Suzuki
Journal: International Journal of Antennas and Propagation
 (2016)

---

## Abstract

The estimation of direction-of-arrival (DOA) of signals is a basic and important problem in sensor array signal processing. To solve this problem, many algorithms have been proposed, among which the Stochastic Maximum Likelihood (SML) is one of the most concerned algorithms because of its high accuracy of DOA. However, the estimation of SML generally involves the multidimensional nonlinear optimization problem. As a result, its computational complexity is rather high. This paper addresses the issue of reducing computational complexity of SML estimation of DOA based on the Alternating Minimization (AM) algorithm. We have the following two contributions. First using transformation of matrix and properties of spatial projection, we propose an efficient AM (EAM) algorithm by dividing the SML criterion into two components. One depends on a single variable parameter while the other does not. Second when the array is a uniform linear array, we get the irreducible form of the EAM criterion (IAM) using polynomial forms. Simulation results show that both EAM and IAM can reduce the computational complexity of SML estimation greatly, while IAM is the best. Another advantage of IAM is that this algorithm can avoid the numerical instability problem which may happen in AM and EAM algorithms when more than one parameter converges to an identical value.


---

## Body

## 1. Introduction

The localization of multiple signal sources by a passive sensor array is of great importance in a wide variety of fields, such as radar, geophysics, radio-astronomy, biomedical engineering, communications, and underwater acoustics. The basic problem in this context is to estimate direction-of-arrival (DOA) of narrow-band signal sources located in the far field of the array.A number of super-resolution techniques have been introduced, such as the Maximum Likelihood (ML) method [1–7], MUSIC [8, 9], ESPRIT [10], Weighted Subspace Fitting (WSF) [11], and the Bayesian method [12].Among these techniques, the ML technique is well known for its high accuracy of DOA. There are two famous ML criterions, that is, Deterministic or Conditional Maximum Likelihood (DML) [1, 3–5, 7] and Stochastic or Unconditional ML (SML) [2, 5–7]. The difference between them lies in their signal models. In particular, the SML shows much higher resolution than other techniques. Furthermore, the SML technique can solve coherent signals without any preprocessing, such as the spatial smoothing [9]. However, the SML technique does not become popular in practice because the estimation of SML generally involves the multidimensional nonlinear optimization problem, and it requires high computational complexity.To solve the multidimensional nonlinear optimization problem, the Alternating Minimization (AM) algorithm is one of the most classic techniques. It is an iterative technique and usually needs one-dimensional global search in the updating process. There are also some other efficient techniques, such as Alternating Projection (AP) [3], Expectation Maximization (EM) [13], Space Alternating Generalized EM (SAGE) [14], genetic algorithm [15], ant colony algorithm [16], and Particle Swarm Optimization (PSO) [17, 18]. Generally, they are all iterative techniques and usually are formulated for explicit criterions (e.g., these algorithms formulated for the criterions of DML and WSF). As for the SML criterion, the AM algorithm is still the most commonly used algorithm although its computational complexity is a little high [6].Therefore, this paper addresses the issue of reducing computational complexity of SML estimation of DOA based on the conventional AM algorithm.Firstly we show a brief description of SML estimation of DOA and the conventional solving method, AM algorithm. Then using transformation of matrix and properties of spatial projection, we propose an efficient AM (EAM) algorithm by dividing the AM criterion of SML into two components. One depends on a single variable parameter while the other does not. The computational complexity of EAM can be greatly reduced compared to AM algorithm. However, numerical instability may happen in calculation of EAM criterion when more than one parameter converges to an identical value. As a result, oscillation may happen in the convergence phase and extra calculation is needed. To solve this problem and reduce computational complexity further, based on the EAM criterion we get the irreducible form of the EAM criterion (IAM) using a uniform linear array. In this way, the EAM criterion can be written into polynomial forms. The common zeros can be easily canceled in numerator and denominator of the EAM criterion. Furthermore, the computational complexity is also reduced. Finally, simulation results are also shown to demonstrate the validity of the proposed EAM and IAM algorithms.The rest of this paper is organized as follows. In Section2 we introduce the problem formulation of DOA. A brief introduction of exact definition of SML estimation for incoherent signals is shown in Section 3. In Section 4, we introduce the conventional AM algorithm and our proposed two efficient algorithms, that is, EAM and IAM algorithms. Simulation results are shown in Section 5 and conclusion is drawn in Section 6.


## 2. Problem Formulation

Without loss of generality, consider that there arep sensors and q narrow-band sources far from the array, centered around a known frequency, impinging on the sensor array from distinct directions θ1,θ2,…,θq, with respect to a reference point, respectively.Note that the received signals may be coherent because of multipath propagation. In the case where there are signals coherent, the independent signal number is less thanq. The task of DOA estimation is to detect all q directions. Also note that here we assume that the signals are narrow-band. For wideband signals, the CSM algorithms [19] can be used as a preprocessing technique to change them into the narrow-band.Furthermore, the sensor configuration can be arbitrary and we assume that all the sensors are omnidirectional and not coupled [6, 7].Using complex envelope representation, thep-dimensional vector received by the array can be expressed as(1)xt=∑k=1qaθkskt+nt,where sk(t) is the kth signal received at a certain reference point. n(t) is a p-dimensional noise vector. a(θ) is the “steering vector” of the array towards direction θ, which is represented as(2)aθ=a1θe-jω0τ1θ,…,apθe-jω0τpθT,where ai(θ) is the amplitude response of the ith sensor to a wave-front impinging from the direction θ. τi(θ) is the propagation delay between the ith sensor and the reference point. The superscript T denotes the transpose of a matrix.In the matrix notation, (1) can be rewritten as(3)xt=AΘst+nt,AΘ=aθ1aθ2⋯aθq,st=s1ts2t⋯sqtT,Θ=θ1θ2⋯θq.Suppose that the received vectorx(t) is sampled at M time instants t1,t2,…,tM and define the matrix of the sampled data as(4)X=xt1xt2⋯xtM=AΘS+N, where(5)S=st1,st2,…,stM,N=nt1,nt2,…,ntM.The problem of DOA finding is to be stated as follows. Given the sampled dataX, obtain a set of estimated directions(6)Θ^=θ^1θ^2⋯θ^q of θ1,θ2,…,θq.


## 3. SML Estimation

In this section, brief descriptions of the exact SML criterion are shown [6, 7].To solve the problem of ML estimation of DOA, we make the following assumptions.(A1)
The array configuration is known and anyp steering vectors for different q directions are linearly independent; that is, the matrixA(Θ) has full rank.(A2)
p, q, and M satisfy the condition that a unique solution of DOA exists in the noise-free case. As for a uniform linear array, q≤2ηp/2η+1 and M≥η are the necessary conditions of the uniqueness [20], where η is the rank ofS.(A3)
p and q are known.(A4)
The noise samplesn(ti) are statistically independent Gaussian random vectors with zero mean and the covariance matrix σ2Ip, where Ip is a p×p identity matrix.(A5)
s(ti) are independent samples from a complex Gaussian random vector which has zero mean and signal covariance matrix with rank η and is independent of the noise.According to [6], the SML criterion can be derived based on the process of the array covariance matrixR defined as(7)R=EXXH=AΘSRAΘH+σ2I, where SR=E{SSH} is the signal covariance matrix.The SML criterion is shown as follows [6, 7].(8)Θ^SML=arg⁡minΘ⁡LSMLΘ,(9)LSMLΘ=det⁡RSS×1p-qtr⁡RNNp-q,(10)RNN=VNHΘR^VNΘ,(11)RSS=VSHΘR^VSΘ,where(12)R^=1MXXH is the sample covariance matrix of the sampled data. VS(Θ) is a p×q matrix composed of an orthonormal system of the signal subspace spanned by A(Θ). VN(Θ) is a p×(p-q) matrix composed of an orthonormal system of the noise subspace, which is an orthogonal complement of the signal subspace. The q×q matrix, RSS, corresponds to the covariance matrix of the components for x(t) in the signal subspace. RNN is the covariance matrix of the components for x(t) in the noise subspace.From (8) and (9), we can see that the estimation of SML is to find a set of Θ which minimizes LSML in (9). This is a multidimensional nonlinear optimization problem.Note that there are literatures [21, 22] discussing the preciseness of SML criterion derived above. This paper focuses on the problem of reducing the computational complexity of SML estimation. Our efficient algorithms are derived based on (9).


## 4. Efficient AM Algorithms for SML Estimation

The AM algorithm is the most classic estimation algorithm for multidimensional nonlinear optimization problem in DOA estimation. In this section, we will introduce the conventional AM algorithm firstly, and then we will introduce our proposed efficient AM algorithms.


### 4.1. Conventional AM Algorithm

The AM method is a popular iterative technique for solving a nonlinear multivariate minimization problem with a multimodal criterion [6]. It can be applied to the SML criterion of DOA in the following manner.LetLk(Θ^(k)) be a cost function of (9) for which the signal number is assumed to be k instead of q, where Θ^(k)={θ^1,θ^2,…,θ^k}.Initialization Phase. For a certain value of SNR, first assuming a single signal, k=1, find θ^1 minimizing L1(Θ^(1)) by one-dimensional global search with respect to θ^1. Next, assuming two signals, k=2, and fixing θ^1 at the value obtained for a single signal, find θ^2 minimizing L2(Θ^(2)) by one-dimensional global search with respect to θ^2. Continue in this fashion until all the initial values for θ^k, k=1,2,…,q are computed.Convergence Phase. Repeat the following updating process until all parameters are converged. At each updating process, let one parameter, say θ^k, be variable and let all other parameters be held fixed. Find θ^k minimizing the criterion Lq(Θ^(q)) by one-dimensional global search with respect to θ^k. Change the index k of the parameter to be updated into (kmodq)+1.Although a global minimum is not guaranteed in the AM algorithm, global solutions can be obtained in most cases because of one-dimensional global searches performed in each update process.


### 4.2. Efficient AM (EAM) Algorithm for SML Estimation

In this section, we propose an efficient version of the AM algorithm. We call it the EAM algorithm. Using transformation of matrix and properties of spatial projection, the EAM algorithm divides the SML criterion into two components. One depends on a variable parameter and the other one is independent of the variable parameter. In order to simplify expressions, parameters to be estimated are represented without the accent hat, and the argumentΘ or Θ^ is omitted.In each updating process, letθl be a variable parameter and define(13)Θl=θ1,θ2,…,θl-1,θl+1,…,θq,Al=aθ1,…,aθl-1,aθl+1,…,aθq.V S l is an orthonormal system of the subspace spanned by Al.V N l is an orthonormal system of the orthogonal complement of the subspace spanned by Al.(14)PAl=AlAlHAl-1AlH=VSlVSlH,(15)PAl⊥=I-PAl=VNlVNlH,(16)RSl=VSlHRVSl,(17)RNl=VNlHRVNl,(18)vlθl=PAl⊥aθlPAl⊥aθl.Note that when the value of θl changes, only vl(θl) varies and all others above are fixed. From these definitions, we have(19)span⁡VNl=span⁡vlθl⊕span⁡VN,span⁡VS=span⁡vlθl⊕span⁡VSl,where ⊕ represents the direct sum of subspaces. It follows from (19) that there exist a p-q+1×p-q+1 unitary matrix T1 and a q×q unitary matrix T2 which satisfy(20)VNl=vlθlVNT1,(21)VS=vlθlVSlT2.Substituting (20) into (17), we have(22)RNl=T1HvlHθlRvlθlvlHθlRVNVNHRvlθlVNHRVNT1=T1HvlHθlRvlθlvlHθlRVNVNHRvlθlRNNT1.Taking the trace of both sides in (22), we have(23)tr⁡RNl=tr⁡RNN+vlHθlRvlθl. Substituting the definition of PAl⊥ and vl(θl) into (23), we have(24)tr⁡RNl=tr⁡RNN+aHθlVNlRNlVNlHaθlaHθlVNlVNlHaθl.Substituting (21) into (11), we have(25)RSS=T2HvlHθlRvlθlvlHθlRVSlVSlHRvlθlVSlHRVSlT2=T2HvlHθlRvlθlvlHθlRVSlVSlHRvlθlRSlT2.Taking the determinant of both sides in (25) and using the definition of (16), we have(26)det⁡RSS=det⁡RSl×vlHθlRvlθl-vlHθlRVSlRSl-1VSlHRvlHθl. Substitute the definition of PAl⊥ and vl(θl) into (26), and define(27)W=RNl-VNlHRVSlRSl-1VSlHRVNl, and then we have(28)det⁡RSS=det⁡RSlaHθlVNlWVNlHaθlaHθlVNlVNlHaθl.Define(29)uθl=VNlHaθland then the final form of the proposed algorithm, EAM algorithm, can be derived from (24) and (28) as follows:(30)tr⁡RNNΘ=tr⁡RNl-uHθlRNluθluθl2,(31)det⁡RSSΘ=det⁡RSluHθlWuθluθl2.In (29), (30), and (31), all quantities except for u(θl) are fixed and can be computed before starting the one-dimensional search with respect to θl. Therefore main computations of each step in the one-dimensional search are a product of the matrix VNlH and the vector a(θl) in (29) and evaluation of two Hermitian forms in (30) and (31). Evaluation of u(θl) requires O(p×p-q+1) multiplications. Then evaluation of u(θl)2 requires O(p-q+12) multiplications. Therefore, at each step of one-dimensional search with respect to the variable parameter, θl, O(p×(p-q+1)) multiplications are required. Therefore, computational complexity can be greatly reduced since, in the conventional AM algorithm, we need to calculate the whole SML criterion in each step of one-dimensional search.(  1) Oscillation Problem in EAM Algorithm. Define the subspace(32)UΘ=span⁡aθ1aθ2⋯aθq.The linear combination(33)aθj-aθiθj-θi always belongs to U(Θ). Even when θj→θi, we have dim⁡UΘ∣θi→θj=q, where dim⁡ represents the dimension of the subspace.On the other hand, whenθj=θi, we have dim⁡UΘ∣θi=θj<q. This implies that the criterion has discontinuous points in the parameter space.Next, we consider calculation of the EAM criterionvl(θl) in (18). When the variable parameter θl approaches the value of a fixed parameter θi, vl(θl) vanishes and both the numerator and denominator in (18) become zero. Then vl(θl) becomes indefinite. Therefore the calculation of vl(θl) becomes numerically instable. This can be verified in Section 5.In the case that the DOA can be solved, the numerical instability does not occur, since each parameter in the convergence phase of the EAM criterion comes apart from others. However, at the threshold region, when more than one signal approaches an identical value, the numerical instability becomes significant. In practice, when this case happens, the sequence of DOA obtained in the convergence phase of the EAM criterion shows oscillation that is because the estimated directions can not converge well due to the numerical instability and would oscillate around that identical value.Let us give an example. In the simulation there are two sources located in 0 and 8 degrees (the true DOAs). When SNR = 0 dB and with specific noise samples, the solutions of SML criterion are 3.999 and 4.001 degrees (the solutions of SML). However, in this case, when calculatingθ1=θ2=4 degrees, the oscillation happens. As a result, the estimated DOAs may be 4 and 4 degrees (the estimated DOAs). In this case, because of oscillation, extra computation is required. Furthermore, the estimated DOAs are wrong (because they are not the solutions of SML criterion) although they are very close to the solutions of SML criterion. As a result, the estimation accuracy is also affected. These explanations will be shown in Section 5.To solve these problems, the key point is to cancel the common zeros in the numerator and denominator of the EAM criterion. Next we try to establish the irreducible form of the EAM criterion using a uniform linear array.


### 4.3. Irreducible Form of Efficient AM Criterion (IAM)

In this section, we derive the irreducible form of the EAM criterion using a uniform linear array. With a uniform linear array, the EAM criterion can be written into polynomial forms. Then we can easily cancel the common zero in both numerator and denominator. Thus numerical instability never happens in IAM criterion. Furthermore, we can find that the IAM algorithm can reduce the computational order of EAM criterion from square to one in each updating process.The array configuration is uniform linear array composed of omnidirectional sensors, of which steering vector is represented as(34)aθ=1e-jϕθ⋯e-jp-1ϕθT,(35)ϕθ=2πΔλsin⁡θ, where λ is the wavelength of signals impinging on the array and Δ is the sensor spacing between two adjacent sensors. As a necessary condition that a unique direction θ is determined by the phase parameter ϕ, Δ≤λ/2 is imposed on the array configuration. In this paper, Δ=λ/2.Using the uniform linear array, we derive the irreducible form of (30) and (31). Define(36)flθ=uHθlRNluθluθl2,f~lθ=uHθlWuθluθl2which are the varying parts in (30) and (31).First we derive the irreducible form of (30). Substituting (15) and (29) into fl(θ), we get(37)flθ=aHθlPAl⊥RPAl⊥aθlaHθlPAl⊥aθl.Using the uniform linear array defined above, the steering vectora(θ) in (34) can be represented as follows:(38)aHθ=ζpz=1z⋯zp-1z=e-jϕθ. Then fl(θ) can be rewritten into the form of a rational function(39)flθ=ζpzPAl⊥RPAl⊥ζp∗zζpzPAl⊥ζp∗z=NzDzz=e-jϕθ, where ζp∗(z) is the paraconjugate of ζp(z) defined as(40)ζp∗z=ζpH1z∗=1z-1⋯z-p-1T and the superscript ∗ is the complex conjugate of a complex number.Letθl be variable and all other parameters are held fixed as well as in last subsection. As we have discussed in the problem of EAM, when θl becomes equal to θi, both the polynomials N(z) and D(z) have double zeros at z=e-jϕθi in the complexz-plane, since it holds ζpe-jϕθiPAl⊥=0. Without canceling these common zeros, fl(θ) is indefinite at θ=θi.The irreducible form of the EAM criterion offl(θ) can be derived by canceling these common zeros. First, we define the polynomial Wl(z) having zeros at z=e-jϕθi, i=1,2…,l-1,l+1,…,q,(41)Wlz=∏i=1,i≠lqz-e-jϕθi=w0+w1z+⋯+wq-1zq-1.Using the coefficients ofWl(z), define the following p×p-q+1 matrix as(42)Wl=w00⋯0w1w0⋯0⋮⋮⋱0wq-1wq-2⋱w00wq-1⋱w1⋮⋮⋱⋮00⋯wq-1. Then we have aHθiWl=0 for i=1,…,l-1,l+1,…,q. Since the column vectors in Wl are all orthogonal to a(θi), the projection matrix PAl⊥ can be written as(43)PAl⊥=WlGl-1WlH,Gl=WlHWl. Here, note that Wl(z) represents a polynomial, while Wl represents a matrix.Using the expressions(44)ζp-q+1zWlz=WlzzWlz⋯zp-qWlz=Wlzζp-q+1z,ζp-q+1z=1z⋯zp-q, the irreducible form of fl(θ) is derived as(45)flθ=Wlzζp-qzNlζp-q∗zWlzWlzζp-qzDlζp-q∗zWlzz=e-jϕθ=ζp-qzNlζp-q∗zζp-qzDlζp-q∗zz=e-jϕθ, where Wl(z) is the common zero factor at z=e-jϕ(θ), i≤q, i≠l, and(46)Nl=Gl-1WlHRWlGl-1=ni,ji,j=0p-q,Dl=Gl-1=di,ji,j=0p-q.As for the irreducible form off~l(θ), it can be derived like this form similarly,(47)f~lθ=ζp-qzN~lζp-q∗zζp-qzDlζp-q∗zz=e-jϕθ, where(48)N~l=Nl-Gl-1WlHRVSlRSl-1VSlHRWlGl-1=n~i,ji,j=0p-q.Define the following polynomials(49)nlz=n0+n1z+⋯+np-qzp-q,dlz=d0+d1z+⋯+dp-qzp-q, where(50)n0=∑i=0p-qni,i,nm=2∑i=mp-qni,i-m,d0=∑i=0p-qdi,i,dm=2∑i=mp-qdi,i-m,m=1,2,…,p-q,and similarly define n~l(z) to calculate the matrix N~l which realizes the same function as nl(z).Then we have the final form offl(θ) and f~l(θ) shown as follows:(51)flθ=Re⁡nlejθRe⁡dlejθ,f~lθ=Re⁡n~lejθRe⁡dlejθ,where Re⁡{} represents the real part of the complex value.Therefore the irreducible form of efficient AM criterion (IAM) is shown like this.(52)tr⁡RNN=tr⁡RNl-Re⁡nlejθRe⁡dlejθ,det⁡RSS=det⁡RSlRe⁡n~lejθRe⁡dlejθ.Furthermore, we can find thatRNl, RSl, Nl, N~l, and Dl can be calculated before each updating process. Therefore, at each updating process we only need to evaluate Re⁡nlejθ, Ren~lejθ, and Re{dl(ejθ)}. Evaluation of them only needs O(p-q) multiplications. Compared to O(p(p-q)) multiplications in each updating process of EAM algorithm, the computational complexity can be reduced further.Here we should note that the IAM algorithm is only applicable to the ULA that is because the IAM algorithm is formulated based on (34) and (35). The steering vector can only be written into this form when the array is uniform linear array (ULA). Based on (34) and (35), the EAM criterion can be written into polynomial form. And then the irreducible form of EAM can be derived by canceling the common zero in both numerator and denominator. Therefore, the IAM is only applicable to the ULA.


## 4.1. Conventional AM Algorithm

The AM method is a popular iterative technique for solving a nonlinear multivariate minimization problem with a multimodal criterion [6]. It can be applied to the SML criterion of DOA in the following manner.LetLk(Θ^(k)) be a cost function of (9) for which the signal number is assumed to be k instead of q, where Θ^(k)={θ^1,θ^2,…,θ^k}.Initialization Phase. For a certain value of SNR, first assuming a single signal, k=1, find θ^1 minimizing L1(Θ^(1)) by one-dimensional global search with respect to θ^1. Next, assuming two signals, k=2, and fixing θ^1 at the value obtained for a single signal, find θ^2 minimizing L2(Θ^(2)) by one-dimensional global search with respect to θ^2. Continue in this fashion until all the initial values for θ^k, k=1,2,…,q are computed.Convergence Phase. Repeat the following updating process until all parameters are converged. At each updating process, let one parameter, say θ^k, be variable and let all other parameters be held fixed. Find θ^k minimizing the criterion Lq(Θ^(q)) by one-dimensional global search with respect to θ^k. Change the index k of the parameter to be updated into (kmodq)+1.Although a global minimum is not guaranteed in the AM algorithm, global solutions can be obtained in most cases because of one-dimensional global searches performed in each update process.


## 4.2. Efficient AM (EAM) Algorithm for SML Estimation

In this section, we propose an efficient version of the AM algorithm. We call it the EAM algorithm. Using transformation of matrix and properties of spatial projection, the EAM algorithm divides the SML criterion into two components. One depends on a variable parameter and the other one is independent of the variable parameter. In order to simplify expressions, parameters to be estimated are represented without the accent hat, and the argumentΘ or Θ^ is omitted.In each updating process, letθl be a variable parameter and define(13)Θl=θ1,θ2,…,θl-1,θl+1,…,θq,Al=aθ1,…,aθl-1,aθl+1,…,aθq.V S l is an orthonormal system of the subspace spanned by Al.V N l is an orthonormal system of the orthogonal complement of the subspace spanned by Al.(14)PAl=AlAlHAl-1AlH=VSlVSlH,(15)PAl⊥=I-PAl=VNlVNlH,(16)RSl=VSlHRVSl,(17)RNl=VNlHRVNl,(18)vlθl=PAl⊥aθlPAl⊥aθl.Note that when the value of θl changes, only vl(θl) varies and all others above are fixed. From these definitions, we have(19)span⁡VNl=span⁡vlθl⊕span⁡VN,span⁡VS=span⁡vlθl⊕span⁡VSl,where ⊕ represents the direct sum of subspaces. It follows from (19) that there exist a p-q+1×p-q+1 unitary matrix T1 and a q×q unitary matrix T2 which satisfy(20)VNl=vlθlVNT1,(21)VS=vlθlVSlT2.Substituting (20) into (17), we have(22)RNl=T1HvlHθlRvlθlvlHθlRVNVNHRvlθlVNHRVNT1=T1HvlHθlRvlθlvlHθlRVNVNHRvlθlRNNT1.Taking the trace of both sides in (22), we have(23)tr⁡RNl=tr⁡RNN+vlHθlRvlθl. Substituting the definition of PAl⊥ and vl(θl) into (23), we have(24)tr⁡RNl=tr⁡RNN+aHθlVNlRNlVNlHaθlaHθlVNlVNlHaθl.Substituting (21) into (11), we have(25)RSS=T2HvlHθlRvlθlvlHθlRVSlVSlHRvlθlVSlHRVSlT2=T2HvlHθlRvlθlvlHθlRVSlVSlHRvlθlRSlT2.Taking the determinant of both sides in (25) and using the definition of (16), we have(26)det⁡RSS=det⁡RSl×vlHθlRvlθl-vlHθlRVSlRSl-1VSlHRvlHθl. Substitute the definition of PAl⊥ and vl(θl) into (26), and define(27)W=RNl-VNlHRVSlRSl-1VSlHRVNl, and then we have(28)det⁡RSS=det⁡RSlaHθlVNlWVNlHaθlaHθlVNlVNlHaθl.Define(29)uθl=VNlHaθland then the final form of the proposed algorithm, EAM algorithm, can be derived from (24) and (28) as follows:(30)tr⁡RNNΘ=tr⁡RNl-uHθlRNluθluθl2,(31)det⁡RSSΘ=det⁡RSluHθlWuθluθl2.In (29), (30), and (31), all quantities except for u(θl) are fixed and can be computed before starting the one-dimensional search with respect to θl. Therefore main computations of each step in the one-dimensional search are a product of the matrix VNlH and the vector a(θl) in (29) and evaluation of two Hermitian forms in (30) and (31). Evaluation of u(θl) requires O(p×p-q+1) multiplications. Then evaluation of u(θl)2 requires O(p-q+12) multiplications. Therefore, at each step of one-dimensional search with respect to the variable parameter, θl, O(p×(p-q+1)) multiplications are required. Therefore, computational complexity can be greatly reduced since, in the conventional AM algorithm, we need to calculate the whole SML criterion in each step of one-dimensional search.(  1) Oscillation Problem in EAM Algorithm. Define the subspace(32)UΘ=span⁡aθ1aθ2⋯aθq.The linear combination(33)aθj-aθiθj-θi always belongs to U(Θ). Even when θj→θi, we have dim⁡UΘ∣θi→θj=q, where dim⁡ represents the dimension of the subspace.On the other hand, whenθj=θi, we have dim⁡UΘ∣θi=θj<q. This implies that the criterion has discontinuous points in the parameter space.Next, we consider calculation of the EAM criterionvl(θl) in (18). When the variable parameter θl approaches the value of a fixed parameter θi, vl(θl) vanishes and both the numerator and denominator in (18) become zero. Then vl(θl) becomes indefinite. Therefore the calculation of vl(θl) becomes numerically instable. This can be verified in Section 5.In the case that the DOA can be solved, the numerical instability does not occur, since each parameter in the convergence phase of the EAM criterion comes apart from others. However, at the threshold region, when more than one signal approaches an identical value, the numerical instability becomes significant. In practice, when this case happens, the sequence of DOA obtained in the convergence phase of the EAM criterion shows oscillation that is because the estimated directions can not converge well due to the numerical instability and would oscillate around that identical value.Let us give an example. In the simulation there are two sources located in 0 and 8 degrees (the true DOAs). When SNR = 0 dB and with specific noise samples, the solutions of SML criterion are 3.999 and 4.001 degrees (the solutions of SML). However, in this case, when calculatingθ1=θ2=4 degrees, the oscillation happens. As a result, the estimated DOAs may be 4 and 4 degrees (the estimated DOAs). In this case, because of oscillation, extra computation is required. Furthermore, the estimated DOAs are wrong (because they are not the solutions of SML criterion) although they are very close to the solutions of SML criterion. As a result, the estimation accuracy is also affected. These explanations will be shown in Section 5.To solve these problems, the key point is to cancel the common zeros in the numerator and denominator of the EAM criterion. Next we try to establish the irreducible form of the EAM criterion using a uniform linear array.


## 4.3. Irreducible Form of Efficient AM Criterion (IAM)

In this section, we derive the irreducible form of the EAM criterion using a uniform linear array. With a uniform linear array, the EAM criterion can be written into polynomial forms. Then we can easily cancel the common zero in both numerator and denominator. Thus numerical instability never happens in IAM criterion. Furthermore, we can find that the IAM algorithm can reduce the computational order of EAM criterion from square to one in each updating process.The array configuration is uniform linear array composed of omnidirectional sensors, of which steering vector is represented as(34)aθ=1e-jϕθ⋯e-jp-1ϕθT,(35)ϕθ=2πΔλsin⁡θ, where λ is the wavelength of signals impinging on the array and Δ is the sensor spacing between two adjacent sensors. As a necessary condition that a unique direction θ is determined by the phase parameter ϕ, Δ≤λ/2 is imposed on the array configuration. In this paper, Δ=λ/2.Using the uniform linear array, we derive the irreducible form of (30) and (31). Define(36)flθ=uHθlRNluθluθl2,f~lθ=uHθlWuθluθl2which are the varying parts in (30) and (31).First we derive the irreducible form of (30). Substituting (15) and (29) into fl(θ), we get(37)flθ=aHθlPAl⊥RPAl⊥aθlaHθlPAl⊥aθl.Using the uniform linear array defined above, the steering vectora(θ) in (34) can be represented as follows:(38)aHθ=ζpz=1z⋯zp-1z=e-jϕθ. Then fl(θ) can be rewritten into the form of a rational function(39)flθ=ζpzPAl⊥RPAl⊥ζp∗zζpzPAl⊥ζp∗z=NzDzz=e-jϕθ, where ζp∗(z) is the paraconjugate of ζp(z) defined as(40)ζp∗z=ζpH1z∗=1z-1⋯z-p-1T and the superscript ∗ is the complex conjugate of a complex number.Letθl be variable and all other parameters are held fixed as well as in last subsection. As we have discussed in the problem of EAM, when θl becomes equal to θi, both the polynomials N(z) and D(z) have double zeros at z=e-jϕθi in the complexz-plane, since it holds ζpe-jϕθiPAl⊥=0. Without canceling these common zeros, fl(θ) is indefinite at θ=θi.The irreducible form of the EAM criterion offl(θ) can be derived by canceling these common zeros. First, we define the polynomial Wl(z) having zeros at z=e-jϕθi, i=1,2…,l-1,l+1,…,q,(41)Wlz=∏i=1,i≠lqz-e-jϕθi=w0+w1z+⋯+wq-1zq-1.Using the coefficients ofWl(z), define the following p×p-q+1 matrix as(42)Wl=w00⋯0w1w0⋯0⋮⋮⋱0wq-1wq-2⋱w00wq-1⋱w1⋮⋮⋱⋮00⋯wq-1. Then we have aHθiWl=0 for i=1,…,l-1,l+1,…,q. Since the column vectors in Wl are all orthogonal to a(θi), the projection matrix PAl⊥ can be written as(43)PAl⊥=WlGl-1WlH,Gl=WlHWl. Here, note that Wl(z) represents a polynomial, while Wl represents a matrix.Using the expressions(44)ζp-q+1zWlz=WlzzWlz⋯zp-qWlz=Wlzζp-q+1z,ζp-q+1z=1z⋯zp-q, the irreducible form of fl(θ) is derived as(45)flθ=Wlzζp-qzNlζp-q∗zWlzWlzζp-qzDlζp-q∗zWlzz=e-jϕθ=ζp-qzNlζp-q∗zζp-qzDlζp-q∗zz=e-jϕθ, where Wl(z) is the common zero factor at z=e-jϕ(θ), i≤q, i≠l, and(46)Nl=Gl-1WlHRWlGl-1=ni,ji,j=0p-q,Dl=Gl-1=di,ji,j=0p-q.As for the irreducible form off~l(θ), it can be derived like this form similarly,(47)f~lθ=ζp-qzN~lζp-q∗zζp-qzDlζp-q∗zz=e-jϕθ, where(48)N~l=Nl-Gl-1WlHRVSlRSl-1VSlHRWlGl-1=n~i,ji,j=0p-q.Define the following polynomials(49)nlz=n0+n1z+⋯+np-qzp-q,dlz=d0+d1z+⋯+dp-qzp-q, where(50)n0=∑i=0p-qni,i,nm=2∑i=mp-qni,i-m,d0=∑i=0p-qdi,i,dm=2∑i=mp-qdi,i-m,m=1,2,…,p-q,and similarly define n~l(z) to calculate the matrix N~l which realizes the same function as nl(z).Then we have the final form offl(θ) and f~l(θ) shown as follows:(51)flθ=Re⁡nlejθRe⁡dlejθ,f~lθ=Re⁡n~lejθRe⁡dlejθ,where Re⁡{} represents the real part of the complex value.Therefore the irreducible form of efficient AM criterion (IAM) is shown like this.(52)tr⁡RNN=tr⁡RNl-Re⁡nlejθRe⁡dlejθ,det⁡RSS=det⁡RSlRe⁡n~lejθRe⁡dlejθ.Furthermore, we can find thatRNl, RSl, Nl, N~l, and Dl can be calculated before each updating process. Therefore, at each updating process we only need to evaluate Re⁡nlejθ, Ren~lejθ, and Re{dl(ejθ)}. Evaluation of them only needs O(p-q) multiplications. Compared to O(p(p-q)) multiplications in each updating process of EAM algorithm, the computational complexity can be reduced further.Here we should note that the IAM algorithm is only applicable to the ULA that is because the IAM algorithm is formulated based on (34) and (35). The steering vector can only be written into this form when the array is uniform linear array (ULA). Based on (34) and (35), the EAM criterion can be written into polynomial form. And then the irreducible form of EAM can be derived by canceling the common zero in both numerator and denominator. Therefore, the IAM is only applicable to the ULA.


## 5. Simulations

In this section, we show some simulation results to demonstrate the validity of the EAM and IAM algorithms. In simulation, the array configuration is a uniform linear array as discussed above.The SNR is defined as(53)SNRk=10log10⁡Eskt2σ2.The Root-Mean-Square-Error (RMSE) is defined as(54)RMSE=1qN∑k=1q∑n=1Nθ^k,n-θk2, where θ^k,n is the estimation of θk at the nth trial. The number of trails in our simulation is 100.In Figure1, the scenario is p=3, q=2, η=2, SNR=10, M=100. Two true sources are independently located at 0 and 8 degrees. The estimated bearing θ2 is fixed at 10 degrees, while θ1 varies from 9.99999 to 10.00001. The dashed line represents the value of SML with EAM criterion, while the solid line represents the IAM criterion. It shows clearly that numerical instability occurs when the EAM criterion is used. In particular, the value changes violently around the point θ2=θ1=10 degrees because it is indefinite. As for the IAM criterion, we can find that the value becomes monotonic and it is numerically stable.Figure 1
Numerical instability in EAM and stability in IAM criterion.Due to the numerical instability, oscillation may happen as Figure2 shows. The scenario is the same as Figure 1. In Figures 2(a) and 2(b), the estimated two bearings θ1 and θ2 converge to an identical value represented by the solid line and the dashed line. The convergence condition is that when the variation of each bearing is less than 10-5 or when the iteration reaches the maximum number, 800. When the EAM criterion is used, the iteration does not stop until it reaches the maximum number because of oscillation. As for the IAM criterion, it converges well for less than 100 iterations. Figure 2(c) shows the oscillation rate of the two criterions in 100 independent trials. We can find that there is no oscillation when the IAM criterion is used. Note that the oscillation rate of EAM decreases when SNR increases as shown in Figure 2(c) that is because the oscillation happens when more than one parameter converges to an identical value. In simulations, two sources are independently located in 0 and 8 degrees. When SNR is relatively low, the case that the estimated two bearings are extremely close (e.g., 3.990 and 4.010) or even equal exists. In this case, oscillation happens. When SNR increases, the estimated two bearings will significantly separate each other. As a result, the oscillation rate will decrease.Figure 2
Oscillation in EAM algorithm while not in IAM.
(a)
The iteration does not stop until it reaches the maximum number, 800, using EAM criterion (b)
The iteration stops at about 85 using IAM criterion (c)
Oscillation rate of EAM and IAM algorithmsFigure3 shows the comparison of RMSE between SML, MUSIC, and ESPRIT. The scenario is the same as Figure 2 except for SNR. In simulation we use the original AM and our proposed EAM and IAM algorithms for SML estimation of DOA. We find that the RMSE of all the three algorithms are completely the same. The reason is that the proposed IAM and EAM algorithms are just transformations for SML criterion and there is no any modification. All these algorithms can find the solution of SML successfully (of course we have to delete the oscillation cases for AM and EAM criterions. Obviously, oscillation affects the estimation accuracy because the bearings do not converge to their optimal value). Also Figure 3 shows that the solution of SML is much better than that of MUSIC and ESPRIT.Figure 3
Comparison of RMSE between SML, MUSIC, and ESPRIT.Figures4 and 5 and Table 1 show the comparison of computational complexity of AM, EAM, and IAM. Figures 4 and 5 show the average amount of operations of each algorithm. In Figure 4, p=3, q=2, η=2, M=100. The true DOA is 0 and 8 degrees. In Figure 5 p=8, q=5, η=3, M=100. The true DOA is 0, 8, 16, 24, and 32 degrees. Note that in Figure 5 η<q, which means that there are correlated signals (coherent case).Table 1
Comparison of computational complexity of AM, EAM, and IAM. The scenario is the same as Figure5 when SNR = 5 dB.
Average iteration times Main computational process AM 127 127 ∗ ( 180 / 0.5 + 0.5 ∗ 2 / 0.001 + 0.001 ∗ 2 / 0.00001 ) = 83,820 times of calculation of LSML (det⁡RSS and tr⁡RNN) EAM 127 127 times of calculation ofdet⁡RSl, tr⁡RNl and W + 83,820 times of calculation of u(θl) IAM 127 127 times of calculation ofdet⁡RSl,tr⁡RNl,Nl,N~l and Dl + 83,820 times of calculation of Re⁡nlejθ, Re⁡n~lejθ, and Re⁡dlejθFigure 4
Average amount of operations for AM, EAM, and IAM algorithms.p=3, q=2, η=2.Figure 5
Average amount of operations for AM, EAM, and IAM algorithms.p=8, q=5, η=3.In Figures4 and 5, “amount of arithmetic operations” represents the summation of all the complex addition, subtraction, multiplication, and division. These two figures show clearly that both EAM and IAM can reduce the computational complexity of SML estimation greatly, while IAM is the best.Table1 shows the detailed comparison of computational complexity including average iteration times and main calculation process of Figure 5 when SNR = 5 dB for 30 independent trials. For all these three algorithms, in the one-dimensional global search of each updating process, we use different step-sizes for searching. At first, we have a relatively large step-size (0.5 degrees) for rough search (note that this step-size should not be too large; otherwise the one-dimensional global minimum may be missed) and then much smaller step-sizes (0.001 and 0.00001 degrees) for fine search. For example, for the AM algorithm in the one-dimensional global search of an updating process, the searching range for the variable parameter (θi) is from −90° to 90°. As a result, the rough search needs 180/(0.5) times of calculation of LSML. And then the fine search needs (0.5∗2/(0.001)+0.001∗2/(0.00001)) times of calculation of LSML. Since the number of iterations is 127, the main computational complexity of AM algorithm is 83,820 times of calculation of LSML. For the EAM and IAM algorithms, the main computational complexity is also shown in Table 1 in the same manner. As we have discussed above, for the EAM and IAM algorithm, before each updating process, we can calculate many components in advance (the detailed components are shown in Table 1). We do not need to calculate the whole cost function LSML every time. As a result, the computational complexity of EAM and IAM is much smaller than that of AM algorithm.Furthermore, from Figures4 and 5, we can find two phenomena. First, the efficiency of EAM and IAM is more obvious than that of AM when the number of parameters increases. Here the number of parameters represents the number of incident signals, that is, q. The task of DOA is to find the directions of all the incident signals. When there are more parameters, the AM algorithm needs more iterations to solve the multidimensional nonlinear optimization problem of SML. For our proposed algorithms (EAM and IAM), in each updating process the computational complexity is greatly reduced. As a result, the efficiency is more obvious when the number of parameters increases (in Figure 5 there are 5 parameters to be estimated, and in Figure 4 there are 2). Second, the computational complexity of all these three methods changes with SNR that is because the computational complexity of all these methods mainly depends on the iteration times for the initial value converges to the estimated value. As we show in Section 4.1, the initial value is determined by the [Initialization Phase]. Obviously, the initial value will change according to different SNRs. Therefore, the total iteration times will change with SNR. As a result, computational complexity also changes with SNR.In a huge number of simulations, we have confirmed the efficiency of our proposed EAM and IAM algorithms, while IAM is the best.


## 6. Conclusions

In this paper, to reduce the computational complexity of SML estimation, based on the AM algorithm, we propose two more efficient algorithms, that is, EAM and IAM algorithms. The EAM algorithm mainly uses transformation of matrix and properties of spatial projection to divide the SML criterion into two components. One is variable, while the other is fixed. Computational complexity can be greatly reduced since the fixed part can be calculated once in advance and only the varying part should be calculated in each one-dimensional updating process. To avoid numerical instability of EAM (note that because of the numerical instability, wrong estimation of DOA may be got) and to reduce computational complexity further, we derive the irreducible form of EAM, that is, IAM algorithm. The main idea of IAM is to use a uniform linear array and rewrite the EAM criterion into polynomial forms. Then the irreducible form can be got by canceling the common zero factor of the polynomial form. Simulation results show that the IAM algorithm can avoid the numerical instable problem of EAM and reduce computational complexity further.



---

*Source: 4926496-2016-11-28.xml*