# An Efficient DOA Estimation Algorithm Based on Diagonal-Symmetric Loading

Authors: Yangyang Xie; Biao Wang; Feng Chen
Journal: Mathematical Problems in Engineering
 (2021)

---

## Abstract

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. Then, 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.


---

## Body

## 1. Introduction

The direction of arrival (DOA) estimation has become an important research subject in the array signal processing [1–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]. These 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]. Therefore, 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). This 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. This 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)
The 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 describedThe paper is organized as follows. In Section2, 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.


## 2. Data Model and Problem Formulation

Notations: this paper usesE[·] to represent the expected operation, [·]T and [·]H represent the transpose and conjugate transpose, respectively, and rank (·) represents the rank operation.


### 2.1. Data Model

The article assumes that the independent narrowband signals (S1, S2, …, Sk) are incident on the M-element uniform linear array (ULA) from θ1,θ2,⋯,θK, which has been shown in Figure 1. The 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.Figure 1
Schematic diagram of source signals at different angles using ULA.By defining a matrix, the received signal can be simply expressed as(1)Xt=AΘSt+Nt,where Xt=x1t,x2t,…,xKtT, AΘ=aθ1,aθ2,…,aθKT, St=s1t,s2t,…,sMtT, and Nt=n1t,n2t,…,nMtT 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(2)aθi=1,e−j2πdsinθi/ϑ,…,e−j2πM−1dsinθi/ϑT,1≤i≤K,where ϑ represents the wavelength of the signal. Therefore, the covariance matrix of the received signal can be expressed as(3)R=EXtXHt=AΘEStSHtAHΘ+ENtNHt=AΘRSAHΘ+σn2I,where RS=EStSHt represents the signal covariance matrix, σn2 is the power of the white Gaussian noise, and I represents the identity matrix.


### 2.2. Problem Formulation

In equation (3), R is a positive definite matrix. In equation (4), λi and ei are the eigenvalue and eigenvector of R. Moreover, λi is specified in the decreasing order:(4)Rei=λiei,i=1,2,…,M,λ1≥λ2≥⋯≥λK≥λK+1=⋯=λM=σn2.The range space ofA (space spanned by columns of A) is called the signal subspace, and it can be verified that(5)RangeA=Spane1,…,eM.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. Thus, they are orthogonal to the last (M−K) eigenvectors ejj=1,2,…,M, and the above subspace spanned by eK+1,…,eM 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. Therefore, in order to reduce the impact of snapshot number on DOA estimation, a covariance matrix is constructed:(6)R¯=R+ασn2I+βσ2aθaHθ,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(7)R¯eI∗=λi∗ei∗,i=1,2,…,M,λ1∗≥λ2∗≥⋯≥λK∗≥λK+1∗=⋯=λM∗,where ei∗ denotes the corresponding eigenvector of λi∗. Since the order of magnitude of β is greater than that of α, then(8)λi∗≥λi.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:(9)λi∗−ασn2=λi=σn2,i=K+1,…,M.Whenθ is not the actual source direction, it does not have the above property. The property stated does not depend on the value of scalars α,β explicitly. Then, the source directions in our proposed algorithm can be estimated as follows:(1)
Compute the eigenvaluesλii=1,…,M of correlation matrix R of equation (3), and they should satisfy equation (4)(2)
Selectα,β by RMSE(3)
ComputeR¯ given in equation (6) and its eigenvalues λi∗i=1,2,…,M(4)
Obtain the direction of sources that satisfy equation (9)Theoretically, the above covariance matrixR 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. Thus, 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^. The main steps of the algorithm are shown as follow.Step 1:
computeR^ by equation (3)Step 2:
compute the eigenvaluesλ^ii=1,…,M of R^ that are in descending orderStep 3:
substituteR by R^ in equation (6), and select αandβ. Then, compute R¯^ for all possible values of theta. That is, R¯^=R^+ασn2I+βσ2aθaHθ.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(10)Fθ=1∑k=K+1Mλ^k∗−λ^k.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 Xt, the maximums of (10) are obtained when the thetas are set to actual direction of sources. This 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. The 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 matrixR¯ defined in (6) is not diagonal. If some sources are fully correlated, that is, rankR¯<M. Then, 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 rankR¯=M. In such a case, the MUSIC method can still be applied; however, its performance highly worsens [26, 27]. Therefore, 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.


### 2.3. Improved K-Means Method for the Source Number Estimation

Source number estimation is equivalent toR¯ 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. The optimization problem then has the form(11)minimizefPk=∑i=1k∑xi∈cidλi,mi,subject toPk=c1,c2,c3,…,ck,where mi=1/nii=1,…,k, ni is the number of data items in the cluster Ci, and dλi,mi represents the distance from xi to mi. The common calculation process of K-means is shown in Figure 2.Figure 2
Common calculation process of K-means.In this paper, two improved points, (a) and (b) in Figure2, are used by combining with the particularity of this scheme (Table 1).(a)
The basic idea of the improved k-means algorithm is to generate a database containing data objects. The number of clusters is set 2, first; two fixed objects are selected as the initial two cluster centers. The 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)
The 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. Therefore, the inverse distance is used to improve the classification accuracy. As shown in Table1, 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%.Table 1
Classification results under the different value ofK and M with 10000 Monte Carlo test.
K-means(K, M)(2, 6)(2, 7)(2, 8)(2, 9)(3, 6)(3, 7)(3, 8)(3, 9)(3, 10)(3, 11)Common method0%0%70.80%100.00%0%0%0%0.02%98.98%100.00%Proposed method100.00%100.00%100.00%100.00%0%9.15%99.97%100.00%100.00%100.00%


### 2.4.α 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:(12)R=σ02aθ0aHθ0+∑i=1M−1σi2aθiaHθi+σn2I,where σ02aθ0aHθ0 is the expected signal component and ∑i=1M−1σi2aθiaHθi 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:(13)R¯=σ02aθ0aHθ0+∑i=1M−1σi2aθiaHθi+σn2I+ασn2I+βσ2aθaHθ.When the symmetric loading part is divided intoθ=θ0 and θ≠θ0, the above formula can be expressed as follows:(14)R¯=σ02+βσ02aθ0aHθ0+∑i=1M−1σi2aθiaHθi+σn2+ασn2I,θ=θ0,R¯=σ02aθ0aHθ0+∑i=1M−1σi2aθiaHθi+σn2+ασn2I+βσI2aθaHθ,θ≠θ0,θ=θI,where σI2 is the symmetrical loading disturbance power.The signal-to-interference-to-noise ratio (SINR) in the covariance matrixR¯ after diagonal-symmetric loading (DSL) is(15)SINR=σ02+βσ02∑i=1M−1σi2+σn2+ασn2,θ=θ0,SINR=σ02∑i=1M−1σi2+σn2+ασn2+βσI2,θ≠θ0,θ=θI,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. The relationship between α and β is further discussed in Section 3 of this article.


## 2.1. Data Model

The article assumes that the independent narrowband signals (S1, S2, …, Sk) are incident on the M-element uniform linear array (ULA) from θ1,θ2,⋯,θK, which has been shown in Figure 1. The 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.Figure 1
Schematic diagram of source signals at different angles using ULA.By defining a matrix, the received signal can be simply expressed as(1)Xt=AΘSt+Nt,where Xt=x1t,x2t,…,xKtT, AΘ=aθ1,aθ2,…,aθKT, St=s1t,s2t,…,sMtT, and Nt=n1t,n2t,…,nMtT 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(2)aθi=1,e−j2πdsinθi/ϑ,…,e−j2πM−1dsinθi/ϑT,1≤i≤K,where ϑ represents the wavelength of the signal. Therefore, the covariance matrix of the received signal can be expressed as(3)R=EXtXHt=AΘEStSHtAHΘ+ENtNHt=AΘRSAHΘ+σn2I,where RS=EStSHt represents the signal covariance matrix, σn2 is the power of the white Gaussian noise, and I represents the identity matrix.


## 2.2. Problem Formulation

In equation (3), R is a positive definite matrix. In equation (4), λi and ei are the eigenvalue and eigenvector of R. Moreover, λi is specified in the decreasing order:(4)Rei=λiei,i=1,2,…,M,λ1≥λ2≥⋯≥λK≥λK+1=⋯=λM=σn2.The range space ofA (space spanned by columns of A) is called the signal subspace, and it can be verified that(5)RangeA=Spane1,…,eM.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. Thus, they are orthogonal to the last (M−K) eigenvectors ejj=1,2,…,M, and the above subspace spanned by eK+1,…,eM 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. Therefore, in order to reduce the impact of snapshot number on DOA estimation, a covariance matrix is constructed:(6)R¯=R+ασn2I+βσ2aθaHθ,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(7)R¯eI∗=λi∗ei∗,i=1,2,…,M,λ1∗≥λ2∗≥⋯≥λK∗≥λK+1∗=⋯=λM∗,where ei∗ denotes the corresponding eigenvector of λi∗. Since the order of magnitude of β is greater than that of α, then(8)λi∗≥λi.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:(9)λi∗−ασn2=λi=σn2,i=K+1,…,M.Whenθ is not the actual source direction, it does not have the above property. The property stated does not depend on the value of scalars α,β explicitly. Then, the source directions in our proposed algorithm can be estimated as follows:(1)
Compute the eigenvaluesλii=1,…,M of correlation matrix R of equation (3), and they should satisfy equation (4)(2)
Selectα,β by RMSE(3)
ComputeR¯ given in equation (6) and its eigenvalues λi∗i=1,2,…,M(4)
Obtain the direction of sources that satisfy equation (9)Theoretically, the above covariance matrixR 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. Thus, 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^. The main steps of the algorithm are shown as follow.Step 1:
computeR^ by equation (3)Step 2:
compute the eigenvaluesλ^ii=1,…,M of R^ that are in descending orderStep 3:
substituteR by R^ in equation (6), and select αandβ. Then, compute R¯^ for all possible values of theta. That is, R¯^=R^+ασn2I+βσ2aθaHθ.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(10)Fθ=1∑k=K+1Mλ^k∗−λ^k.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 Xt, the maximums of (10) are obtained when the thetas are set to actual direction of sources. This 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. The 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 matrixR¯ defined in (6) is not diagonal. If some sources are fully correlated, that is, rankR¯<M. Then, 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 rankR¯=M. In such a case, the MUSIC method can still be applied; however, its performance highly worsens [26, 27]. Therefore, 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.


## 2.3. Improved K-Means Method for the Source Number Estimation

Source number estimation is equivalent toR¯ 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. The optimization problem then has the form(11)minimizefPk=∑i=1k∑xi∈cidλi,mi,subject toPk=c1,c2,c3,…,ck,where mi=1/nii=1,…,k, ni is the number of data items in the cluster Ci, and dλi,mi represents the distance from xi to mi. The common calculation process of K-means is shown in Figure 2.Figure 2
Common calculation process of K-means.In this paper, two improved points, (a) and (b) in Figure2, are used by combining with the particularity of this scheme (Table 1).(a)
The basic idea of the improved k-means algorithm is to generate a database containing data objects. The number of clusters is set 2, first; two fixed objects are selected as the initial two cluster centers. The 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)
The 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. Therefore, the inverse distance is used to improve the classification accuracy. As shown in Table1, 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%.Table 1
Classification results under the different value ofK and M with 10000 Monte Carlo test.
K-means(K, M)(2, 6)(2, 7)(2, 8)(2, 9)(3, 6)(3, 7)(3, 8)(3, 9)(3, 10)(3, 11)Common method0%0%70.80%100.00%0%0%0%0.02%98.98%100.00%Proposed method100.00%100.00%100.00%100.00%0%9.15%99.97%100.00%100.00%100.00%


## 2.4.α 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:(12)R=σ02aθ0aHθ0+∑i=1M−1σi2aθiaHθi+σn2I,where σ02aθ0aHθ0 is the expected signal component and ∑i=1M−1σi2aθiaHθi 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:(13)R¯=σ02aθ0aHθ0+∑i=1M−1σi2aθiaHθi+σn2I+ασn2I+βσ2aθaHθ.When the symmetric loading part is divided intoθ=θ0 and θ≠θ0, the above formula can be expressed as follows:(14)R¯=σ02+βσ02aθ0aHθ0+∑i=1M−1σi2aθiaHθi+σn2+ασn2I,θ=θ0,R¯=σ02aθ0aHθ0+∑i=1M−1σi2aθiaHθi+σn2+ασn2I+βσI2aθaHθ,θ≠θ0,θ=θI,where σI2 is the symmetrical loading disturbance power.The signal-to-interference-to-noise ratio (SINR) in the covariance matrixR¯ after diagonal-symmetric loading (DSL) is(15)SINR=σ02+βσ02∑i=1M−1σi2+σn2+ασn2,θ=θ0,SINR=σ02∑i=1M−1σi2+σn2+ασn2+βσI2,θ≠θ0,θ=θI,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. The relationship between α and β is further discussed in Section 3 of this article.


## 3. 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. Then, the compared experiment among MVDR, NSSP, and proposed algorithm is conducted. The 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.Figure 3
Comparison of spatial spectrum of algorithm.


### 3.2.α 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. X-axis and Y-axis represent the logarithm of β and α, respectively.Figure 4
RMSE under two different snapshot numbers (n = 100 (a), 300 (b)).
(a)(b)As seen in Figure4, 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. Thus, 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.


### 3.3. 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+θ) ° by keeping other conditions unchanged. The theta can be expressed as an arithmetic sequence, which starts at 5 and increase by 2 each time until 25. Therefore, the theta can be regarded as the angle difference between the two targets. Then, 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.Figure 5
Comparison of RMSE of algorithms under different DOA intervals.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. Then, the comparative analysis experiment of the three algorithms on RMSE is conducted. The results are shown in Figure6. The 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. The 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. The results of the three algorithms tend to the CRB curve, which shows the effectiveness of these methods.Figure 6
Comparison of RMSE of algorithms under different SNR.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. The RMSE is shown in Figure7. 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. The 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. Thus, the proposed algorithm performs better than the other two methods. The results of the three algorithms tend to the CRB curve, which also proves their effectiveness.Figure 7
Comparison of the RMSE of the algorithm under different snapshots.


## 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. Then, the compared experiment among MVDR, NSSP, and proposed algorithm is conducted. The 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.Figure 3
Comparison of spatial spectrum of algorithm.


## 3.2.α 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. X-axis and Y-axis represent the logarithm of β and α, respectively.Figure 4
RMSE under two different snapshot numbers (n = 100 (a), 300 (b)).
(a)(b)As seen in Figure4, 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. Thus, 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.


## 3.3. 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+θ) ° by keeping other conditions unchanged. The theta can be expressed as an arithmetic sequence, which starts at 5 and increase by 2 each time until 25. Therefore, the theta can be regarded as the angle difference between the two targets. Then, 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.Figure 5
Comparison of RMSE of algorithms under different DOA intervals.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. Then, the comparative analysis experiment of the three algorithms on RMSE is conducted. The results are shown in Figure6. The 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. The 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. The results of the three algorithms tend to the CRB curve, which shows the effectiveness of these methods.Figure 6
Comparison of RMSE of algorithms under different SNR.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. The RMSE is shown in Figure7. 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. The 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. Thus, the proposed algorithm performs better than the other two methods. The results of the three algorithms tend to the CRB curve, which also proves their effectiveness.Figure 7
Comparison of the RMSE of the algorithm under different snapshots.


## 4. Conclusion

This 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. The compared experiments proved that the proposed algorithm has higher resolution and better estimation accuracy. The conclusions are summarized as follows:(1)
The 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)
The 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



---

*Source: 4432040-2021-12-07.xml*