An Efficient DOA Estimation Algorithm Based on Diagonal-Symmetric Loading

In order to solve the problem that the subspace-like direction of arrival (DOA) estimation performs poor due to the error of sources number, this paper proposes a new super-resolution DOA estimation algorithm based on the diagonal-symmetric loading (DSL). Specifically, orthogonality principle of the minimum eigenvector of the specific covariance matrix and the source number estimation based on the improved K-means method were adopted to construct the spatial spectrum. +en, by considering the signal-to-interference-to-noise ratio (SINR), the theoretical basis for selecting parameters was given and verified by numerical experiment. To evaluate the effectiveness of the proposed algorithm, this paper compared it with the methods of minimum variance distortionless response (MVDR) and new signal subspace processing (NSSP). Experimental results prove that the proposed DSL has higher resolution and better estimation accuracy than the MVDR and NSSP.


Introduction
e direction of arrival (DOA) estimation has become an important research subject in the array signal processing [1][2][3][4][5] because it plays an important role in radar, sonar, wireless communication, and many other application systems. DOA estimation aims to detect the direction information of sources from the received signals [6,7]. Researchers have studied how to implement DOA estimation by weighted subspace fitting (WSF) [8,9], multiple signal classification (MUSIC) [10,11], or estimated signal parameter rotation invariance (ESPRIT) [12,13]. ese methods mainly estimated the target angle by considering the orthogonality of signal subspace and noise subspace, which made certain progress.
By analyzing the DOA estimation process, it is essential to obtain the target sources number accurately, but it is not easy to get the target sources number directly [14]. erefore, some researchers adopted a number of different estimation techniques, including Akaike information criterion (AIC) [15], minimum description length (MDL) [16,17], second-order statistic of eigenvalues (SORTE) [18], property of the variance of the rotational submatrix (VTRS) [19], and characteristic threshold method (eigenthreshold, ET) [20]. However, these techniques are required under certain conditions. Besides these techniques, Capon et al. [21] also the proposed minimum variance distortionless response (MVDR) method that can obtain better resolution and strong anti-interference ability, but MVDR is very sensitive to steering vector mismatch. Zhang et al. [22] designed a narrowband DOA estimation algorithm that does not need to consider the number of sources, which can achieve the best performance of the MUSIC algorithm. Qian et al. [23] described a DOA algorithm that uses Toeplitz matrix reconstruction and the maximum eigenvalue without the number of sources, but the Toeplitz matrix reconstruction is required. Qian et al. [24] propose a coherent DOA estimation algorithm based on a new principal-singular-vector utilization for modal analysis (PUMA). is method uses linear prediction theory to transform the DOA estimation problem into a univariate polynomial root problem.
Although the above works have obtained some valuable results, there are still some improvement in the resolution and accuracy. is paper makes full use of the orthogonality between the steering vector and the noise subspace [25] and constructs a new spatial spectrum by considering the characteristics of the trace of the received data covariance matrix, which can achieve super-resolution effect and better estimation accuracy when the number of sources is unknown. Compared to the existing works, the main contributions of this paper are described as follows: (1) A new super-resolution direction of arrival estimation algorithm based on the diagonal-symmetric loading is proposed (2) e source number estimation is based on the improved K-means method (3) By considering the signal-to-interference-to-noise ratio (SINR), the theoretical basis for selecting parameters in the proposed algorithm is described e paper is organized as follows. In Section 2, the signal and noise model is presented and the DOA estimation problem is formulated; on this basis, the proposed algorithm is presented, and its parameter asymptotic properties are discussed. In Section 3, the simulation experiment of the proposed algorithm is carried out, and the comparison results are given. In Section 4, the work of this paper is summarized and the main conclusions are given.

Data Model and Problem Formulation
Notations: this paper uses E[·] to represent the expected operation, [·] T and [·] H represent the transpose and conjugate transpose, respectively, and rank (·) represents the rank operation.

Data Model.
e article assumes that the independent narrowband signals (S 1 , S 2 , . . ., S k ) are incident on the M-element uniform linear array (ULA) from θ 1 , θ 2 , · · · , θ K , which has been shown in Figure 1. e numbers of signal K are less than the numbers of array elements M, and the element spacing d of the uniform linear array is half a wavelength. Here, the proposed method in the article is also applicable to multivariate arrays of other shapes.
By defining a matrix, the received signal can be simply expressed as where , n 2 (t), . . . , n M (t)] T represent the output matrix, direction response matrix, incident signal matrix, and noise matrix of the ULA array, respectively. For the ULA, the steering vector can be expressed as where ϑ represents the wavelength of the signal. erefore, the covariance matrix of the received signal can be expressed as where represents the signal covariance matrix, σ 2 n is the power of the white Gaussian noise, and I represents the identity matrix. (3), R is a positive definite matrix. In equation (4), λ i and e i are the eigenvalue and eigenvector of R. Moreover, λ i is specified in the decreasing order:

Problem Formulation. In equation
e range space of A (space spanned by columns of A) is called the signal subspace, and it can be verified that In equation (5), each column of A is a steering vector corresponding to a source direction and is equal to the subspace spanned by the first M eigenvector of R. us, they are orthogonal to the last (M−K) eigenvectors e j (j � 1, 2, . . . , M), and the above subspace spanned by e K+1 , . . . , e M is called the noise subspace.
For the sources at fixed directions, the last (M−K) eigenvalues and their corresponding eigenvectors of the new covariance matrix are the same as before when changing their powers (variances). In other word, the noise subspace is invariant to the power of sources.
However, when the number of snapshots is few, the noise and signal spaces are not completely independent. erefore, in order to reduce the impact of snapshot number on DOA estimation, a covariance matrix is constructed: where R is defined by equation (3), a(θ) is M ×1 array response vector of θ in equation (1), σ 2 is the power of the signal, β is the positive constant scalars, and α is a undetermined constant. If λ * i denotes eigenvalue in the decreasing order of R, which can be expressed as where e * i denotes the corresponding eigenvector of λ * i . Since the order of magnitude of β is greater than that of α, then As shown in equation (9), another important property of R is that the remaining eigenvalues of R and R are the same when θ in equation (6) is set to one of the source directions: When θ is not the actual source direction, it does not have the above property. e property stated does not depend on the value of scalars α, β explicitly. en, the source directions in our proposed algorithm can be estimated as follows: (1) Compute the eigenvalues λ i (i � 1, . . . , M) of correlation matrix R of equation (3), and they should satisfy equation (4) (2) Select α, β by RMSE (3) Compute R given in equation (6) and its eigenvalues λ * i (i � 1, 2, . . . , M) (4) Obtain the direction of sources that satisfy equation (9) eoretically, the above covariance matrix R is continuously available, and the direction of sources can be obtained accurately. However, in practice the covariance matrix is estimated from finite number of snapshots (N), so it is discrete. us, the (M−K) smaller eigenvalues of the covariance matrix in equation (4), equation (5), and equation (9) are different.
Based on the above analysis, an efficient DOA algorithm is proposed. Here, the discrete covariance matrix is calculated by equation (3) and defined as R. e main steps of the algorithm are shown as follow.
Step 1: compute R by equation (3) Step 2: compute the eigenvalues λ i (i � 1, . . . , M) of R that are in descending order Step 3: substitute R by R in equation (6), and select α and β. en, compute R for all possible values of theta. at is, R � R + ασ 2 n I + βσ 2 a(θ)a H (θ).
Step 4: for each value of theta, compute the eigenvalues λ * i (i � 1, . . . , M) of R in descending order and calculate spatial spectrum F(θ) as Step 5: the values of θ are the direction angles of sources and correspond to d largest maximums of F(θ).
By virtue of (8), it can be confirmed that F(θ) is a positive function. Combining (8) and (9) together, the denominator of (10) is zero when R � R. Considering the discrete property of the signal X(t), the maximums of (10) are obtained when the thetas are set to actual direction of sources.
is proves the asymptotic consistency of the proposed method. When the number of samples (N) used for estimating R in equation (6) and (7) approaches infinity, the maximums of F(θ) will reach the exact direction of sources. Compared to MUSIC algorithm, the proposed method not only is applicable to all array configuration but also shows insensitive to the power level differences of closely spaced sources. e followings will explain in detail.
When there is a correlation between the directions of sources, such as multipath propagation, or there exist smart jammers in communication applications; their covariance matrix R defined in (6) is not diagonal. If some sources are fully correlated, that is, rank(R) < M. en, the conventional MUSIC and proposed methods fail to predict the DOA of sources. Usually, with symmetric-array configurations employed, both conventional MUSIC and proposed methods can be applied after certain preprocessing such as forward-backward smoothing. As the sources are partially correlated, the matrix R is still not diagonal, but rank(R) � M. In such a case, the MUSIC method can still be applied; however, its performance highly worsens [26,27]. erefore, the proposed algorithm is more suitable for the case of partially correlated sources, and its performance is much better even though the sources are highly correlated.

Improved K-Means Method for the Source Number
Estimation. Source number estimation is equivalent to R eigenvalue classification problem. Hence, K-means can be utilized to solve this problem. It is assumed that, in the K-means clustering problem, there is a set of M eigenvalue values Λ � λ 1 , λ 2 , . . . , λ n to be clustered. e optimization problem then has the form where m i � (1/n i )(i � 1, . . . , k), n i is the number of data items in the cluster C i , and d(λ i , m i ) represents the distance from x i to m i . e common calculation process of K-means is shown in Figure 2.
In this paper, two improved points, (a) and (b) in Figure 2, are used by combining with the particularity of this scheme (Table 1). (a) e basic idea of the improved k-means algorithm is to generate a database containing data objects. e number of clusters is set 2, first; two fixed objects are selected as the initial two cluster centers. e distances between the remaining samples and each cluster center are calculated, and the sample is classified into the nearest cluster center, the new cluster center is calculated by the average method. (b) e usual spatial clustering algorithms are based on various distances, such as Euclidean distance, Manhattan distance, and Mincus distance. However, the commonly used distance may cause misclassification. erefore, the inverse distance is used to improve the classification accuracy. As shown in Table 1, when the number of sources (K) is less than or equal to 1/3 of the number of array elements (M), the classification result reaches 100%.

α and β Value Analysis.
In the proposed algorithm, α and β are required to discuss. According to the feature subspace decomposition theory, the feature decomposition of R is obtained: where σ 2 0 a(θ 0 )a H (θ 0 ) is the expected signal component and is the interference component. Actually, R is often sampled by a limited number of snapshots, which can reduce the performance of DOA resolution. Compared with the ideal interference plus noise covariance matrix, the sampling covariance matrix has two disadvantages. (1) R is obtained by a finite number of snapshots, and the small eigenvalues corresponding to the noise subspace will be disturbed. It can cause that the sidelobes rise and the DOA resolution performance decreases. (2) R contains the desired signal component, which will reduce the performance of DOA resolution, especially in many nonideal situations [28]. It can be seen that suppressing the disturbance of the small eigenvalues corresponding to the noise subspace can not only reduce the sidelobes but also improve the DOA resolution performance and reduce the expected signal component.
According to (6) and (12), the eigendecomposition of R is obtained: When the symmetric loading part is divided into θ � θ 0 and θ ≠ θ 0 , the above formula can be expressed as follows: where σ 2 I is the symmetrical loading disturbance power. e signal-to-interference-to-noise ratio (SINR) in the covariance matrix R after diagonal-symmetric loading (DSL) is It can be seen from (15): (1) Symmetrical loading is equivalent to adding signal and interference components into the covariance matrix at the same time. For the symmetrical loading, when β becomes larger, the main eigenvalue corresponding to the signal subspace is strengthened, and the proportion of interference components will be increased. When θ � θ 0 , the expected signal component and interference component in the covariance matrix are increased by β times simultaneously, the SINR ratio in equation (15) will become larger, and the peak can be obtained well. When θ ≠ θ 0 , the interference component is moved to the denominator, and the proportion of the denominator will increase, which can also effectively suppress the increase in sidelobes. (2) Diagonal loading is equivalent to enhancing the noise component in the covariance matrix. In order to further analyze the relationship between symmetry and diagonal, two extreme cases are considered. (i) When the value of β is very small, only the diagonal loading part is considered temporarily; when α ⟶ − 1, theoretically, the small eigenvalue disturbance corresponding to the noise subspace can be suppressed, and an effective peak can be formed. In fact, the number of snapshots is limited, and the formed noise covariance matrix is not a unit matrix and is greatly affected by the randomness of noise. As a result, sidelobes tend to appear frequently under different noises. (ii) When the value of β is relatively large, the influence of α in equation (15) becomes weak, the sidelobes gradually disappear, and the peak is likely to be significant. At this time, SINR is close to a stable value. When β is much greater than 1, the original information of the signal will be overwritten by the diagonal loading component. In this paper, the appropriate range is decided by experiments.
It can be seen from the above analysis that a suitable α can increase the resolving power of the main peak, and increasing the value of β reasonably can improve the robustness of the algorithm. e relationship between α and β is further discussed in Section 3 of this article.

Simulation Experiment
3.1. Algorithm Feasibility Analysis. We assume two equal power and uncorrelated far-field narrowband signals are incidents on a uniform linear array, which is composed of 6 sound pressure sensors at the arrival angles of 0°and 10°. d � 0.5λ is the distance between the array elements, the number of snapshots is 300, and the SNR is 0 dB. en, the compared experiment among MVDR, NSSP, and proposed algorithm is conducted. e result is shown in Figure 3. It can be found that when the angle difference between the two targets is small, the MVDR algorithm cannot accurately identify the two targets, but the proposed algorithm and the NSSP algorithm can accurately identify the two targets. It indicates that the proposed algorithm and the NSSP algorithm have higher resolution than the MVDR algorithm. Moreover, at this time, the proposed algorithm has relatively flat sidelobes and no wrong peaks, which proves the correctness of (11) and further demonstrates the correctness and feasibility of the proposed algorithm theory.

α and β Value Estimation.
According to the α and β value analysis in Section 2.4, there is a limit to SINR. To estimate the limit situation, the influence of α and β on the root mean square error (RMSE) is analyzed under two different snapshot numbers, which are shown in Figure 4. Xaxis and Y-axis represent the logarithm of β and α, respectively.
As seen in Figure 4, the results are affected by the number of snapshots evidently. When the number of snapshots is low, the effect of β is obvious. While in the case of the high snapshot number, the accuracy of RMSE is improved, and the influence of randomness is aggravated and α is sensitive. Meanwhile, it can be found that the influence of α is very small, and the RMSE gradually decreases with the increase of β. Lower sample frequency generates more random information, which causes the noise subspace and signal subspace are not completely orthogonal. With the decrease of the sampling frequency, the random information takes the place of noise information. us, the addition of α is approximately regarded as a disturbance to the noise. According to (15), the increase of β advances the value of the SINR, which can decrease the noise effect and improve the robustness of the proposed algorithm. In addition, considering the estimation of α in equation (9), α and β are set as −1 and 10, respectively.

Algorithm Statistical Performance Analysis.
In order to verify the performance of the algorithm, 500 independent Monte Carlo experiments were performed. We set the arrival angle of one target to 0°, and the other angle to (0+ θ)°b y keeping other conditions unchanged. e theta can be expressed as an arithmetic sequence, which starts at 5 and increase by 2 each time until 25. erefore, the theta can be regarded as the angle difference between the two targets. en, this paper set the SNR to 0 dB and conducted 500 Monte Carlo experiments to obtain the RMSE of the algorithm under different angle differences. When the target angle difference is less than 10°, it can be analyzed that the intersection between two signal subspaces becomes large. To decrease this effect, this study added β in equation (15) to improve the identification results. Meanwhile, it can be also seen that the performance of the proposed algorithm is the best among the three algorithms in Figure 5. When the angle increases, the RMSE of the three algorithms are the same. Overall, the performance of the proposed algorithm is better than the MVDR algorithm and can obtain the super-resolution effect without knowing the number of sources.
Next, by keeping other conditions unchanged, we set the signal-to-noise ratio (SNR) starts from −5 Db and increased it by 2 dB each time until 15 dB. en, the comparative analysis experiment of the three algorithms on RMSE is conducted. e results are shown in Figure 6. e RMSE of the proposed algorithm and the NSSP algorithm decrease  when the signal-to-noise ratio increases. However, when the SNR of the MVDR algorithm is less than 7 dB, the RMSE decreases more slowly. When the SNR is greater than 7 dB, the RMSE of the algorithm decreases sharply, which means that when the MVDR algorithm is less than 7 dB, the two targets of 0°and 10°cannot be distinguished. e reason is that the decrease of SNR can make DOA estimation more difficult. To solve this problem, this study increases SINR to improve the algorithm performance by introducing β. It also shows that the proposed algorithm has a super-resolution effect. e results of the three algorithms tend to the CRB curve, which shows the effectiveness of these methods. Furthermore, to analyze the impact of the number of snapshots on the algorithm, by keeping other conditions unchanged, we set the target incident angle to 0°and 15°and fix the numbers of snapshots are 100, which will be increased by 200 each time until 900. Finally, 500 Monte Carlo experiments are conducted. e RMSE is shown in Figure 7. As the number of snapshots increases, the RMSE of the MVDR algorithm also decreases. However, the performance of MVDR is the worst among the three algorithms. When the numbers of snapshots are less than 400, the proposed algorithm performs better than the NSSP method. When the numbers of snapshots are over 400, the RMSE of the proposed algorithm is relatively stable and closer to the NSSP algorithm. e above results show that the orthogonality is worse when the snapshots are not many, which is the same with the conclusion from Figure 5, that is, the proposed algorithm is more suitable to the few snapshots. us, the proposed algorithm performs better than the other two methods. e results of the three algorithms tend to the CRB curve, which also proves their effectiveness.

Conclusion
is paper presented an efficient DOA estimation algorithm based on diagonal-symmetric loading, which can hold the invariance property of noise subspace to the source powers.
e compared experiments proved that the proposed algorithm has higher resolution and better estimation accuracy. e conclusions are summarized as follows: (1) e improved K-means method performs better than the original method, especially when the number of array elements is small (2) When the numbers of sources are unknown, the proposed algorithm is more effective than the MVDR and NSSP algorithms (3) e theoretical analysis of α and β selections is given, which can overcome the drawback of false peaks appearing on the scanning spectrum when the feature vector corresponding to the minimum feature value of the noise is orthogonal to the scan steering vector Data Availability e data are available at https://download.csdn.net/ download/Aireyi/20355374.