An Efficient Algorithm for Direction Finding against Unknown Mutual Coupling

In this paper, an algorithm of direction finding is proposed in the presence of unknown mutual coupling. The preliminary direction of arrival (DOA) is estimated using the whole array for high resolution. Further refinement can then be conducted by estimating the angularly dependent coefficients (ADCs) with the subspace theory. The mutual coupling coefficients are finally determined by solving the least squares problem with all of the ADCs utilized without discarding any. Simulation results show that the proposed method can achieve better performance at a low signal-to-noise ratio (SNR) with a small-sized array and is more robust, compared with the similar processes employing the initial DOA estimation and further improvement iteratively.


Introduction
Direction finding of multiple sources has received tremendous attention in the field of radar, sonar, mobile communication, and so on. High resolution algorithms, such as Multiple Signal Classification (MUSIC) [1] and Estimation of Signal Parameters via Rotational Invariance Techniques (ESPRIT) [2], can distinguish closely-spaced sources by the use of subspace theory and, therefore, have been widely used in the past few decades [3]. In spite of the potential advantages of the eigenstructure methods, their direct application to real systems is difficult, due to the critically required precise knowledge of the array manifold. In other words, the performance of the super-resolution techniques is strongly dependent on the array manifold accuracy. In practice, the array manifold is inevitably affected by mutual coupling, position perturbation and array gain or phase uncertainties. This results in significant distortion of the amplitude and phase of the signals received from the array. The direct employment of the eigenstructure-based methods leads to a serious degradation of direction finding [4][5][6][7]. Array calibration is an effective way to alleviate the deviation of direction of arrival (DOA) estimation.
Several calibration algorithms have been discussed in the last few decades. The earliest investigation was made by Schmidt [8] and Weiss [9]. They measured and stored the array manifold directly by interpolation. However, a large amount of memory is required, and the size and cost of the system could obviously be increased. In order to overcome the drawbacks of the above scheme, a kind of calibration method was proposed using a set of calibration sources in known locations. The maximum-likelihood approach [10] can be used to estimate the calibration matrix for array error compensation. Similarly, the algorithm in [11] was developed utilizing an iterative least mean squares approach. Although these techniques can obtain high accuracy and a large scope of calibration, it is impractical to set a collection of calibration sources and get their DOAs exactly as prior knowledge.
An approach to mitigate the influence of array errors is to calibrate the array by the use of the received signals. Such methods for estimating the DOA, unknown coupling, gain and phase simultaneously are called self-calibration. Friedlander and Weiss proposed to use an iterative process to acquire the parameters and DOA [12,13]. Svantesson formulized it as an optimization problem and solved the problem iteratively to estimate the mutual coupling coefficients for coupling compensation in the linear array of dipoles [14]. However, the result will converge to the local optimum if the initial values deviate far from the real ones. Alternatively, another kind of method was proposed conducting multidimensional search based on subspace fitting. The maximization of a posteriori (MAP) estimator proposed in [15] and [16] is one of these algorithms. Although it does not have the problem of convergence, non-linear multidimensional optimization in these methods is computationally consuming, and the convergence rate is relatively slow.
In recent years, several algorithms have been developed based on the characterization of the mutual coupling matrix (MCM). In [17], Ye obtained the initial estimation of DOA by setting the instrumental sensors on each side of the array. As a result, a one-dimensional search of the spatial spectrum can be performed directly using the original array data. Moreover, the result can be refined iteratively by estimating the mutual coupling coefficients for compensation. In [18], the signals received from the middle subarray were directly exploited for the traditional MUSIC, whereas the MCM is assumed to be a complex symmetric Toeplitz matrix. In order to increase the performance of DOA estimation, another method [19] proposed recently takes advantage of the special structure of MCM to parameterize the steering vector. It achieves the estimation of DOA using the whole array and improves the result by mutual coupling compensation.
In this paper, we assume that the MCM is a band-symmetric Toeplitz matrix with finite non-zero elements. By using the parameterized steering vector in [19], the preliminary DOA estimates can be obtained with the whole array. Further compensation of mutual coupling is then made by estimating the ADC first with the orthogonality of the subspace and the mutual coupling coefficients following with the full information of the ADCs. Simulation results show that the proposed algorithm has satisfactory performance compared with the methods in [17] and [19], especially when the SNR is low and small number of sensors exist.
The remainder of the paper is organized as follows. Section 2 is devoted to the problem formulation. Section 3 addresses the direction finding with initial estimation of DOA and coupling compensation by using the subspace theory and all of the parameters estimated. Section 4 gives the concluding remarks.

Problem Formulation
We consider a uniform linear array (ULA) consisting of M sensors with K narrowband far-field sources s 1 (t), s 2 (t), . . . , s K (t) received. The sources come from directions θ 1 , θ 2 , . . . , θ K with respect to the normal line of the array. Assume that the distance between adjacent sensors is d and the wavelength of the carrier is λ. With the interactions between sensors, the mutual coupling effect cannot be ignored. The M × 1 output vector of the array is then given by: . . . , s K (t)] T and n(t) = [n 1 (t), n 2 (t), . . . , n M (t)] T denote the received signal vector, source signal vector and noise vector, respectively. The notation [·] T denotes the transposition. A = [a(θ 1 ), a(θ 2 ), . . . , a(θ K )] is the array manifold matrix, in which a(θ k ) = [1, β(θ k ), β 2 (θ k ), . . . , β M −1 (θ k )] T is the steering vector with β(θ k ) = exp{−j2πd sin θ k /λ}. C ∈ C M ×M is the MCM, which is generally considered to be independent of the angle [20]. It indicates the interactions between arbitrary two sensors of the array. Equation (1) is obtained under the assumption that the additive sensor noises n i (t), i = 1, . . . , M are independent and identically distributed (i.i.d.) white Gaussian with the common variance σ 2 . s i (t) and n i (t) are zero-mean wide-sense stationary random processes. Several studies of the coupling model [12,21] have shown that the coupling between a pair of sensors is nearly the same. Therefore, the MCM is a banded symmetric Toeplitz matrix for the ULA. Furthermore, based on the fact that the mutual coupling between two sensors is inversely proportional to their distance, the coefficient will be zero if two sensors are several wavelength apart. Let c ij = c ji = c |i−j| denote the mutual coupling coefficient between the i-th and the j-th element of the ULA, and assume that there are P distinct non-zero elements in the MCM with the self coupling c 0 normalized as one. Then, the M × M matrix C can be expressed as: where Toeplitz{·} denotes a symmetric Toeplitz matrix constructed by the P × 1 vector c = [1, c 1 , . . . , c P −1 ] T . Define: as the equivalent steering vector of the direction θ. From Equation (1), the covariance matrix of the received signals is: where (·) H denotes the Hermitian transpose operation, is the equivalent array manifold matrix. As demonstrated in [2], the above equation can be rewritten as: using the eigenvalue decomposition (EVD) of R x , where Λ = diag{Λ s , Λ n } is a matrix with M eigenvalues at the diagonal in descending order and zero elsewhere. E s ∈ C M ×K and E n ∈ C M ×(M −K) are the eigenvectors corresponding to the K largest eigenvalues and the M − K smallest eigenvalues, respectively. According to the well-known subspace-based algorithm [1], the signal subspace spanned by the column of A m = CA is orthogonal to the noise subspace spanned by that of E n . Therefore, we have: If the MCM is known, the directions of sources θ 1 , θ 2 , . . . , θ K can be estimated based on the spectrum function P (θ) = 1/ E H n a m (θ) 2 . However, the vector c is usually unknown, under which circumstance the traditional DOA estimation approach cannot be used, and a new method should be investigated for direction finding.

Direction Finding and Mutual Coupling Compensation
In this section, we first reformulate the equivalent steering vector by parameterizing the MCM. Then, we solve the DOA estimation problem in the presence of unknown mutual coupling by applying the whole array data for better performance [19,22]. For further refinement, a new mutual coupling compensation algorithm is proposed by utilizing the special structure of the reformulated MCM with no information lost.

DOA Estimation Using the Whole Array
Let r k represents the k-th element of the M × 1 equivalent steering vector a m (θ). By combining Equations (2) and (3), a m (θ) can be expressed as: where: For notational clarity, we set 0 in the above equations. In the case of g 0 = 0, we can extract g 0 out of r k and describe a m (θ) as [19]: where: and: where µ k and α k are the ADCs. They are functions of c k and β(θ). It has been shown from Equations (8)-(10) that the angularly independent MCM is transformed into an angularly-dependent expression g 0 Γ(θ). Since Γ(θ) is a diagonal matrix with 2(P − 1) unknown variables and a(θ) is a vector, it is feasible to exchange their elements and to reduce the number of ones in Γ(θ) as follows [19]: is an M × (2P − 1) matrix with an (M − 2P + 2) × 1 vector and two (P − 1) × (P − 1) diagonal matrices locatedat the diagonal, respectively. Considering the subspace principle given in Equation (6) and a m (θ) in Equation (11), we have the following equation on condition that g 0 = 0: As mentioned in [19] and [22], when θ is consistent with any one of the K incoming angles, Q(θ) is rank deficient. Therefore, the following spectrum function with the determinant of Q(θ) as the sensitive factor provides an effective means of DOA estimation: It is worth noting that Q(θ) does not contain any information of c. Therefore, the spectrum function given above can be employed even if the mutual coupling is unknown. Compared with the algorithm in [17], this spectrum function takes advantage of the whole array. People do not need to extract the middle array for mutual coupling eliminating. As a result, no information is lost in the spectrum estimation. This method is available only if: in which circumstance, the

Mutual Coupling Compensation
Once we get the DOA estimates, further refinement can be conducted based on the angularly dependent expression of a m (θ) as shown in Equation (11). We assumeθ is the DOA estimated from Equation (14) and that no blind angles exists (i.e., g 0 = 0). Based on the fact that the noise subspace is orthogonal to the column subspace of a m (θ), we have: Notice that the P -th element of v(θ) is one, and the others are ADCs to be estimated. Denote , where q j , j = 1, . . . , 2P − 1 is the j-th column of Q(θ). By replacing the columns of Q(θ) and the associated elements of v(θ), Equation (16) can then be reformed as: where Q (θ) = [q P , q 1 , . . . , q P −1 , q P +1 , . . . , q 2P −1 ] and v (θ) = [1, µ 1 , . . . , µ P −1 , α 1 , . . . , α P −1 ] T . Moreover, extracting q P from Q (θ), Equation (17) can be expressed as: Consequently, we can get the estimates of µ k and α k by solving the above equation as: v (2 : P ) = −Q (:, 2 : 2P − 1) q P where (·) represents the pseudo inverse operation. Now, we estimate the mutual coupling coefficients by utilizing the estimated ADCsμ k ,α k . Based on the observation of µ k , α k in Equations (9) and (10), it is not difficult to find that they are the linear functions of c k , k = 1, . . . , P − 1. Let B(θ) be the coefficient matrix between c = [c 1 , . . . , c P −1 ] T and [μ 1 , . . . ,μ P −1 ,α 1 , . . . ,α P −1 ] T . Then, we have: where: can be obtained easily ifθ is estimated and: From Equations (20) and (22), we can see that g 0 must be determined before the estimation of c k . Although containing the unknown mutual coefficient, g 0 can be easily obtained by employing the particular composition ofμ k ,α k . Notice thatμ 1 andα P −1 have the complementary elements of g 0 . We can therefore get g 0 by:ĝ whereμ 1 andα P −1 have been acquired from Equation (19). Combining Equations (23) and (22), the mutual coupling coefficient vector c can be determined by solving Equation (20) as: For better performance, all of the DOAs estimated in the above subsection can be used to form the extended coefficient matrixB.
and v gi as the matrix and vector evaluated at the i-th estimated DOAθ i . Then, the extension of Equation (20) will be:B Solving Equation (25) by the least squares, we can get a more precise estimation of c as: The above approach provides us a means of mutual coupling compensation. That is to say, once the vector c is determined, the matrix C can be formed by locating its element on the corresponding sub-diagonal. Therefore, DOA estimation can be further obtained by searching the peak of: The performance can be further improved by repeating the above procedure. The proposed algorithm for DOA estimation and mutual coupling compensation can be summarized as follows.
(1) Get L snapshots of the received signal x(t) at t = t 1 , . . . , t L , and form the following matrix as: (2) Generate the covariance matrix using the above data matrix by: (3) Conduct the EVD ofR x , and get the noise subspaceÊ n . (4) Scan the direction from −90 • to 90 • with 1 • as the step size. Calculate the special spectrum using Equation (14), and obtain the initial estimation of DOAθ 1 , . . . ,θ K .
(7) Enhance the DOA estimation with the estimatedĈ and Equation (27). (8) Repeat Step (5) to Step (7) to get a more precise estimation of directions.

Simulation Results
In this section, simulations will be conducted to validate the performance of the proposed method. Consider two independent sources from the far-field incident on the ULA from θ 1 = −10 • and θ 2 = 20 • . Sensors are located in the array with equal spacing d = λ/2. The number of effective mutual coupling coefficients is P = 3. Here, we set c = [1, 0.43301 − 0.25i, 0.14142 − 0.14142i], which is used in the second simulation of [17], to guarantee g 0 = 0 at any direction θ.
In the first simulation, we evaluate the performance of DOA estimation in Step (4) without mutual coupling compensation. Assume that the array number is M = 7 and that the snapshot number is 500. The root mean squared error (RMSE) is used to compare the DOA accuracy of different algorithms. It can be calculated, in general, by: The RMSE as a function of SNR is illustrated in Figure 1a. The method used in our proposed algorithm is superior at low SNRs compared with the method in [17], since the whole array is utilized. As the SNR increases, the RMSE of DOA estimation decreases gradually for all of the methods. When the SNR is greater than 10 dB, the accuracy of the two methods with unknown C is almost the same. Figure 1b shows the effect of the array size on RMSE when SNR = −5 dB. Notice that the choice of M should satisfy Equation (15). It can be seen that our method slightly outperforms the method in [17] when M ≤ 10.
Using the initial DOA estimates, we now proceed to get the mutual coupling vector c. In the second simulation, we consider the same scenario. Define the RMSE of c as Ns n=1 1)], whereĉ k,n represents the estimated coupling coefficient in the n-th Monte Carlo experiment. From Figure 2a, we can see that the proposed method can obtain significant improvement of the coupling coefficient estimation when the SNR is lower than 3 dB. Besides, it is robust compared with the method in [17] and [19]. Now, we keep the SNR at −5 dB and vary the snapshots from 10 to 960. Figure 2b shows that the proposed method can achieve higher accuracy compared with the other two algorithms. From the estimation of the coupling matrix in Steps (5) and (6), it is not difficult to find that the full use ofμ k ,α k leads to the superior performance of the mutual coupling estimation.   Proposed method Ye's method [17] Liao's method [19] (b) In the third simulation, we investigate the DOA estimation obtained by mutual coupling compensation. We consider two sources incoming from θ 1 and θ 2 = θ 1 + ∆θ with SNR = 0 dB arriving at the ULA for M = 7 and L = 500. Define: as the decision condition of whether the anglesθ 1 ,θ 2 can be identified or not. Then, we conduct N s = 200 Monte Carlo experiments for each algorithm and record the distinguishable ones. Denote n s as the number of experiments satisfying the above inequality. The curve of Probability = n s /N s with respect to the angle interval is presented in Figure 3a. Compared with the initial estimation of DOA, as shown in Step (4), the proposed method improves the accuracy dramatically by refinement with mutual coupling compensation. Ye's method [17] is inferior to the new method, because only the middle array is used. Figure 3b-d illustrates the RMSE of the DOA estimation with respect to the SNR, snapshots and array size, respectively. As expected, the results are very promising. It is the full use of x(t),μ k and α k that leads to the good performance of the proposed method. In the above simulations, the number of sensors is assumed to be seven. It is worth mentioning that similar results can be obtained for bigger M , as long as it satisfies M ≤ 10 for P = 3. The only difference is that the superiority is less obvious. Probability with known coupling with unknown C proposed method Ye's method [17] Liao's method [19] (a) with known C with unknown C proposed method Ye's method [17] Liao's method [19] (b) with known C with unknown C proposed method Ye's method [17] Liao's methd [19] (c) In the fourth simulation, we will investigate the effect of mutual coupling on the DOA estimation. First, consider a signal impinging on the ULA from −10 • . Mutual coupling between the adjacent two sensors exists, i.e., P = 2. Then, define r si = C(i, :)a(θ), i = 1, ..., M as the response of the i-th element to the received signal. r s1 is presented in Figure 4a with the amplitude and phase of the coupling coefficient varying. The figure illustrates that when the phase stays near zero, the response becomes higher as the amplitude increases. When the phase grows into the range of [π/2, π], the result will be the inverse to that of [0, π/2). The effect of coupling is to some extent similar to the beamforming. Any change of the amplitude and phase could make the response different.  proposed method−π/3 proposed method−5π/6 with known C−π/3 with known C−5π/6 with unknown C−π/3 with unknown C−5π/6 (b) The RMSE as a function of the amplitude of c 1 is presented in Figure 4b with the phase φ fixed at π/3 and 5π/6, respectively. From this, we can see that the results get better with the increase of |c 1 | on the condition that φ = π/3. The situation gets worse for φ = 5π/6, since the response becomes gradually lower. Figure 4c demonstrates the RMSE as a function of the phase with |c 1 | fixed at 0.15 and 0.35, respectively. The error increases slightly as the phase grows from zero to π, which coincides with the response shown in Figure 4a. Figure Figure 4d we can conclude that the proposed algorithm could achieve better performance when the coupling coefficient has a smaller phase and greater amplitude.
In the fifth simulation, we access the performance of DOA estimation when signals impinge on the ULA with different angle intervals. Figure 5a presents the response of every element to DOA with c = [1, 0.43301 − 0.25i, 0.14142 − 0.14142i]. Figure 5b illustrates the influence of the angle interval to the initial estimation in Step (4). It is shown that not only the angle interval, but also the initial DOA will affect the accuracy of estimation seriously. As long as the DOAs locate in the main lobe of the middle element, the gain of the array will not decline rapidly. Figure 5c gives the RMSE of DOA estimation versus SNR with different angle intervals. With the same initial DOA, the performance obviously decreases when the interval becomes bigger.  In summary, the proposed method can achieve more precise estimates of DOA as opposed to Ye's method [17] and Liao's method [19], especially when the SNR is low and the array size is close to the minimum available value of M = K+2P −1. The reservation of the information obtained guarantees the good performance of the proposed algorithm. The experiments of the mutual coupling effect show that the coupling coefficient, initial DOA and the angle interval could have an influence on the performance of DOA estimation at the same time.

Conclusions
This paper addresses the DOA estimation in the presence of unknown mutual coupling. Based on the subspace theory, the initial DOA is first estimated using the whole array without calibration sources and auxiliary sensors, which leads to high accuracy. With the assumption that no blind angles exist in the space, mutual coupling compensation is further conducted by estimating the coupling coefficients indirectly from the angular-dependent coefficients. Finally, with all of the ADCs utilized without discarding any, the mutual coupling coefficients are determined by solving the least squares problem. Simulations show that the proposed method can achieve better performance at low SNR with a small-sized array. The robustness of the method can be verified, as well.

Author Contributions
Weijiang Wang and Yingtao Ding provided the idea of this work. Weijiang Wang and Shiwei Ren conceived of and designed the experiments. Shiwei Ren performed the experiments and provided all of the figures and data for the paper. Haoyu Wang prepared the literature and analyzed the data. Weijiang Wang and Shiwei Ren wrote the paper. Correspondence and requests for the paper should be addressed to Yingtao Ding.

Conflicts of Interest
The authors declare no conflict of interest.