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.
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 Section 2 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 are sensors and narrow-band sources far from the array, centered around a known frequency, impinging on the sensor array from distinct directions , 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 than . The task of DOA estimation is to detect all 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, the -dimensional vector received by the array can be expressed aswhere is the th signal received at a certain reference point. is a -dimensional noise vector. is the “steering vector” of the array towards direction , which is represented aswhere is the amplitude response of the th sensor to a wave-front impinging from the direction . is the propagation delay between the th sensor and the reference point. The superscript denotes the transpose of a matrix.
In the matrix notation, (1) can be rewritten as
Suppose that the received vector is sampled at time instants and define the matrix of the sampled data as where
The problem of DOA finding is to be stated as follows. Given the sampled data X, obtain a set of estimated directions of .
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 any steering vectors for different directions are linearly independent; that is, the matrix A has full rank.(A2), , and satisfy the condition that a unique solution of DOA exists in the noise-free case. As for a uniform linear array, and are the necessary conditions of the uniqueness [20], where is the rank of S.(A3) and are known.(A4)The noise samples are statistically independent Gaussian random vectors with zero mean and the covariance matrix , where is a identity matrix.(A5) 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 matrix R defined as where is the signal covariance matrix.
The SML criterion is shown as follows [6, 7].where is the sample covariance matrix of the sampled data. is a matrix composed of an orthonormal system of the signal subspace spanned by . is a matrix composed of an orthonormal system of the noise subspace, which is an orthogonal complement of the signal subspace. The matrix, , corresponds to the covariance matrix of the components for in the signal subspace. is the covariance matrix of the components for in the noise subspace.
From (8) and (9), we can see that the estimation of SML is to find a set of which minimizes 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.
Let be a cost function of (9) for which the signal number is assumed to be instead of , where .
Initialization Phase. For a certain value of SNR, first assuming a single signal, , find minimizing by one-dimensional global search with respect to . Next, assuming two signals, , and fixing at the value obtained for a single signal, find minimizing by one-dimensional global search with respect to . Continue in this fashion until all the initial values for , are computed.
Convergence Phase. Repeat the following updating process until all parameters are converged. At each updating process, let one parameter, say , be variable and let all other parameters be held fixed. Find minimizing the criterion by one-dimensional global search with respect to . Change the index of the parameter to be updated into .
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 be a variable parameter and define
is an orthonormal system of the subspace spanned by .
is an orthonormal system of the orthogonal complement of the subspace spanned by .Note that when the value of changes, only varies and all others above are fixed. From these definitions, we havewhere represents the direct sum of subspaces. It follows from (19) that there exist a unitary matrix and a unitary matrix which satisfySubstituting (20) into (17), we haveTaking the trace of both sides in (22), we have Substituting the definition of and into (23), we have
Substituting (21) into (11), we haveTaking the determinant of both sides in (25) and using the definition of (16), we have Substitute the definition of and into (26), and define and then we have
Defineand then the final form of the proposed algorithm, EAM algorithm, can be derived from (24) and (28) as follows:
In (29), (30), and (31), all quantities except for are fixed and can be computed before starting the one-dimensional search with respect to . Therefore main computations of each step in the one-dimensional search are a product of the matrix and the vector in (29) and evaluation of two Hermitian forms in (30) and (31). Evaluation of requires multiplications. Then evaluation of requires multiplications. Therefore, at each step of one-dimensional search with respect to the variable parameter, , 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
The linear combination always belongs to . Even when , we have , where represents the dimension of the subspace.
On the other hand, when , we have . This implies that the criterion has discontinuous points in the parameter space.
Next, we consider calculation of the EAM criterion in (18). When the variable parameter approaches the value of a fixed parameter , vanishes and both the numerator and denominator in (18) become zero. Then becomes indefinite. Therefore the calculation of 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 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 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 , is imposed on the array configuration. In this paper, .
Using the uniform linear array, we derive the irreducible form of (30) and (31). Definewhich are the varying parts in (30) and (31).
First we derive the irreducible form of (30). Substituting (15) and (29) into , we get
Using the uniform linear array defined above, the steering vector in (34) can be represented as follows: Then can be rewritten into the form of a rational function where is the paraconjugate of defined as and the superscript is the complex conjugate of a complex number.
Let 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 becomes equal to , both the polynomials and have double zeros at in the complex z-plane, since it holds . Without canceling these common zeros, is indefinite at .
The irreducible form of the EAM criterion of can be derived by canceling these common zeros. First, we define the polynomial having zeros at , ,
Using the coefficients of , define the following matrix as Then we have for . Since the column vectors in are all orthogonal to , the projection matrix can be written as Here, note that represents a polynomial, while represents a matrix.
Using the expressions the irreducible form of is derived as where is the common zero factor at , , , and
As for the irreducible form of , it can be derived like this form similarly, where
Define the following polynomials whereand similarly define to calculate the matrix which realizes the same function as .
Then we have the final form of and shown as follows:where represents the real part of the complex value.
Therefore the irreducible form of efficient AM criterion (IAM) is shown like this.
Furthermore, we can find that , , , , and can be calculated before each updating process. Therefore, at each updating process we only need to evaluate , , and . Evaluation of them only needs multiplications. Compared to 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
The Root-Mean-Square-Error (RMSE) is defined as where is the estimation of at the th trial. The number of trails in our simulation is 100.
In Figure 1, the scenario is , , , , . Two true sources are independently located at 0 and 8 degrees. The estimated bearing is fixed at 10 degrees, while 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 degrees because it is indefinite. As for the IAM criterion, we can find that the value becomes monotonic and it is numerically stable.
Due to the numerical instability, oscillation may happen as Figure 2 shows. The scenario is the same as Figure 1. In Figures 2(a) and 2(b), the estimated two bearings and 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 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.
(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 algorithms
Figure 3 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.
Figures 4 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, , , , . The true DOA is 0 and 8 degrees. In Figure 5 , , , . The true DOA is 0, 8, 16, 24, and 32 degrees. Note that in Figure 5 , which means that there are correlated signals (coherent case).
In Figures 4 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.
Table 1 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 () is from −90° to 90°. As a result, the rough search needs 180/(0.5) times of calculation of . And then the fine search needs times of calculation of . Since the number of iterations is 127, the main computational complexity of AM algorithm is 83,820 times of calculation of . 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 every time. As a result, the computational complexity of EAM and IAM is much smaller than that of AM algorithm.
Furthermore, from Figures 4 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, . 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.
Competing Interests
The authors declare that they have no competing interests.
Acknowledgments
This paper is supported by the following funds: The National Natural Science Funds, China, with no. 61601519; Fundamental Research Funds for the Central University, China, with nos. 15CX02047A and 15CX05025A; Shandong Provincial Natural Science Foundation, China, with no. ZR2014FM017; Science and Technology Development Plan of Qingdao with no. 15-9-1-80-jch.