# DOA Estimation for Sources with Large Power Differences

Authors: Qingyuan Fang; Mengzhe Jin; Weidong Liu; Yong Han
Journal: International Journal of Antennas and Propagation
 (2021)

---

## Abstract

Sources with large power differences are very common, especially in complex electromagnetic environments. Classical DOA estimation methods suffer from performance degradation in terms of resolution when dealing with sources that have large power differences. In this paper, we propose an improved DOA algorithm to increase the resolution performance in resolving such sources. The proposed method takes advantage of diagonal loading and demonstrates that the invariant property of noise subspace still holds after diagonal loading is performed. We also find that the Cramer–Rao bound of the weak source can be affected by the power of the strong source, and this has not been noted before. The Cramer–Rao bound of the weak source deteriorates as the power of the strong source increases. Numerical results indicate that the improved algorithm increases the probability of resolution while maintaining the estimation accuracy and computational complexity.


---

## Body

## 1. Introduction

Direction-of-arrival (DOA) estimation has been extensively investigated due to its wide applications in radar, sonar, wireless communication, and navigation [1–4]; e.g., DOA estimation is of significance for source location in MIMO radar [5], and the issue of DOA estimation for noncircular sources is valuable in practical communications [6]. Subspace-based methods are preferred among numerous DOA estimation methods due to their high superresolution and ease of implementation [7–9]. For passive radar, the echo from the target may be several dB lower than the direct signal or intentional interference. The estimation performance deteriorates dramatically if the weak target is masked by the strong target [10]. However, the classical MUSIC, Root-MUSIC, and Capon algorithms suffer from severe performance degradation when resolving closely spaced sources with large power differences [11, 12], especially in the threshold region [13].There are two main types of methods for DOA estimation with sources that have power differences. One type involves estimating the weak sources after mitigating the strong interference. Generally, this type of method requires prior knowledge of the strong interference. With prior DOA knowledge of the strong interference, the interference jamming method (IJM) is proposed in [14], where strong interferences are eliminated by transforming the M × N-dimensional array manifold matrix into M × (N − j) dimensions, in which M is the number of array elements, N is the number of incident sources, and j is the number of strong interferences. However, the IJM method suffers from the loss of array aperture because the steering vector belonging to the strong interferences is removed from the array manifold matrix. Another DOA estimation algorithm based on the elimination of interference is proposed in [15], and it shows better performance than the jamming jam method (JJM). Its performance is limited when the power difference between two sources is not large. Provided with accurate prior DOA knowledge, the constrained MUSIC approach [16], the weighted prior-MUSIC approach [17], and the improved L1-SVD method [18] achieve accurate DOA estimation for weak sources, whereas the performance degrades if the prior DOA knowledge of strong interference is inaccurate. Given the number of strong interference, the DOA estimation of weak sources in the presence of strong jamming has been studied based on the extended noise subspace [19]. Prior knowledge of strong interference is needed; otherwise, the extended noise subspace is unavailable. Additionally, a modified spectrum using eigenvector projections has proven valid in dealing with sources exhibiting large power differences based on the number of known sources [20]. However, prior knowledge regarding the input sources is difficult to acquire, especially in noncooperative situations, and limits the applications of methods that need prior knowledge.The other type of DOA estimation method for sources with power differences does not need prior knowledge. However, this type mostly employs numerous iterations at the cost of increased computational complexity. An oblique projection filter can be used to constrain sources from undesired directions without any prior knowledge of strong interference [21]. Multiple singular value decomposition is needed in this method, and this increases the complexity of the algorithm. A hybrid method without any prior knowledge utilizes compressive sensing (CS) to estimate the DOAs of strong sources and then removes the effects of strong sources via orthogonal complements [22]. However, it suffers from resolution loss if the DOAs of strong sources are not accurately estimated. The number of sources and the DOAs of weak sources can be jointly obtained from the Eigenbeam mCapon method [23]. Repeated iterations are employed at the cost of considerable computation complexity, and this is also true for the RELAX algorithm [24]. Recently, sparse theory and sparse arrays have attracted much attention and have provided a new method for DOA estimation [6]. In [25], the authors analyze the loss of weak targets in the sparse asymptotic minimum variance (SAMV) method and propose the robust sparse asymptotic minimum variance (RSAMV) algorithm. A robust orthogonal projection method based on an optimized sparse array is proposed without prior knowledge in which the DOAs of strong and weak targets are precisely estimated [26], but many iterations are still needed.Both iterations and prior knowledge are not required for the DOA estimation algorithm, which utilizes the invariant property of noise subspace (IPNS) [27, 28]. In this paper, we propose an improved DOA estimation algorithm to increase the resolution performance based on the IPNS. The improved algorithm takes advantage of diagonal loading, and we demonstrate that the IPNS still holds after diagonal loading. The Cramer-Rao bound (CRB) of the weak source is studied in the presence of a strong source to examine the influence caused by the power of the strong source. The numerical results indicate that our approach is highly effective in solving two closely spaced sources with large power differences. The rest of this paper is organized as follows. In Section 2, the signal model is introduced. In Section 3, an improved algorithm for DOA estimation for sources with large power differences is proposed. We demonstrate that the IPNS is still valid after diagonal loading and analyze the complexity of the improved algorithm. In Section 4, the CRBs of weak sources under different powers of strong sources are studied. Numerical simulations and discussions are presented in Section 5. Finally, we conclude our work in Section 6.


## 2. Signal Model

ConsiderP independent narrowband sources impinging on a uniform linear antenna array from the far-field, as shown in Figure 1. The distance between any two antenna elements is d = λ0/2, where λ0 is the carrier wavelength. The number of receiver antenna elements is M (M > P). Suppose that the ith source impinges on the uniform linear antenna array from direction θi. The DOAs of all sources are defined as θ=θ1,θ2,…,θP.Figure 1
Uniform linear antenna array for signal reception.The output of the uniform linear antenna array at timet can be expressed as follows:(1)xt=Ast+nt,where xt=x0t,x1t,…,xMtT is the received signal vector and ⋅T denotes the matrix transpose operation. st=s1t,s2t,…,sPtT is the vector of incident signals. nt=n1t,n2t,…,nMtT is the additive and independent Gaussian noise vector whose covariance and mean are σn2I (I denotes the identity matrix) and zeros, respectively. A=aθ1,aθ2,…,aθP is the array steering matrix with aθi=1,e−jwi,e−j2wi,…,e−jM−1wiT and wi=k0dsinθi, where k0=2π/λ0. The covariance matrix of xt is given by the following:(2)R=ExtxtH=ARsAH+σn2I,where E⋅ and ⋅H denote the expectation and the conjugate transpose operators, respectively. Rs=EststH is the covariance matrix of the input signal. Actually, the covariance matrix R is unavailable in real applications and is typically replaced by a finite sample covariance matrix:(3)R^=1N∑t=1NxtxHt,where N is the number of snapshots.


## 3. Improved Algorithm for DOA Estimation for Sources with Large Power Differences

### 3.1. Improved Algorithm Based on the Invariant Property of Noise Subspace

After introducing a virtual source, the noise subspace remains the same as before only when the DOA of the virtual source overlaps with the actual sources. This property, called the invariant property of noise subspace, is demonstrated in [27] and can be used to estimate the DOAs of actual sources. Diagonal loading is a widespread approach to improve robustness against mismatch errors [29, 30]. We find that diagonal loading is also effective in improving the resolution of the IPNS method. Thus, we propose an improved algorithm based on the IPNS. Diagonal loading effectively reduces the power differences between the strong and weak sources, thereby increasing the resolution of the improved algorithm. Next, we certify that the IPNS still holds after diagonal loading so we can use it for DOA estimation.The eigenvalue decomposition (EVD) of the covariance matrixR is given by the following:(4)R=UΣUH=USUNΣSΣNUSUNH,where U=USUN and Σ=diagΣS,ΣN; diag⋅ denotes the diagonal matrix; and US, UN, ΣS, and ΣN are the eigenvectors and eigenvalues of the signal subspace and noise subspace, respectively. The space spanned by the impinging sources and the noise are orthogonal; that is, spanUS⊥spanUN, where span⋅ denotes the spanned space. ΣS is a diagonal matrix ΣS=diagλ1,λ2,…,λP, where(5)λ1≥λ2≥⋯≥λP,are the eigenvalues of the signal subspace. ΣN is also a diagonal matrix ΣN=diagλP+1,λP+2,…,λM, where(6)λP+1=λP+2=⋯=λM=σn2,are the eigenvalues of the noise subspace. Suppose that a virtual source impinges on the receiving antenna array from the direction θ′. The virtual source does not physically exist but is mathematically needed. Then, the new covariance matrix including the virtual source can be defined as follows:(7)R′≜R+kaθ′aθ′H,where k > 0 is the power of the virtual source. The EVD of the new covariance matrix R′ can be given by the following:(8)R′=U′Σ′U′H=US′UN′ΣS′ΣN′US′UN′H,where U′=US′UN′ and Σ′=diagΣS′,ΣN′, US′ and UN′ are the eigenvectors of the signal subspace and the noise subspace, respectively, and ΣS′=diagλ1′,λ2′,…,λP′ and ΣN′=diagλP+1′,,λP+2′,…,λM′ are the signal and noise eigenvalues corresponding to the new covariance matrix R′, respectively. Suppose that the direction of the virtual source θ′ overlaps with the ith actual source, that is, θ′=θi,i∈1,2,3,…,P. Define a selective unit vector ei as follows:(9)ei≜0,0,…,1,…,0T,where ei is a unit vector consisting of P elements and all the elements are 0 except the ith element, which is 1. Then, the new covariance matrix R′ can be written as follows:(10)R′=ARsAH+kaθ′aθ′H+σn2I=ARsAH+kAeieiTAH+σn2I=ARS′AH+σn2I,where RS′m,n=Rsm,n+k when m=n=i; otherwise, RS′m,n=Rsm,n. It is evident in equation (10) that the virtual source does not affect the noise power but only increases the power of the ith source when the direction θ′ overlaps with that of the ith actual source. Thus, we can obtain that the noise subspace eigenvalues ΣN′ remain consistent with ΣN when θ′ overlaps with the ith actual source, that is,(11)λP+1=λP+2=⋯=λM=λP+1′=λP+2′=⋯=λM′=σn2.When the direction of the virtual source is different from those of any of the actual sources, that is,θ′∉θ1,θ2,θ3,…,θP, then R′ is as follows:(12)R′=ARsAH+kaθ′aθ′H+σn2I=Aaθ′Rs00kAHaθ′H+σn2I.In this case, the noise subspace eigenvaluesΣN′ are different from ΣN. The eigenvalues of R′ in descending order are as follows:(13)λ1′≥λ2′≥⋯≥λP′≥λP+1′≥λP+2′=⋯=λM′=σn2.The eigenvalues of the noise subspace in equation (13) are ΣN′=diagλP+2′,λP+3′,…,λM′, which are different from ΣN=diagλP+1,λP+2,…,λM.Diagonal loading has been demonstrated to be an effective approach for improving robustness in beamforming. We employ diagonal loading to improve the resolution of the IPNS method. After employing diagonal loading, the covariance matrix can be expressed as follows:(14)R¯=R+σ2I,where σ2 is the coefficient of the diagonal loading. After the EVD of R¯, we can obtain the following:(15)R¯=UΣ¯UH=USUNΣ¯SΣ¯NUSUNH=UΣUH+σ2I=UΣ+σ2IUH,where Σ¯=diagΣ¯S,Σ¯N. It is apparent in equation (15) that the diagonal matrix Σ¯N of the noise subspace can be written as follows:(16)Σ¯N=ΣN+σ2I=diagλP+1+σ2,…,λM+σ2.After introducing the virtual source, we define the covariance matrixR¯′ based on R¯ as follows:(17)R¯′≜R¯+kaθ′aθ′H=R+σ2I+kaθaθ′H.Similarly, it can also be certified that the virtual source does not affect the noise subspace ofR¯′ compared to that of R¯ if the direction θ′ overlaps with the actual sources. Assume that the virtual source impinges on the receiving antenna array from the direction θ′. Then, R¯′ can be rewritten as follows:(18)R¯′=R+σ2I+kaθ′aθ′H=ARsAH+kAeieiTAH+σn2I+σ2I=ARs′AH+σn2+σ2I=R′+σ2I,where R′=ARS′AH+σn2I. The EVD of R¯′ is given by the following:(19)R¯′=U′¯Σ′¯U′¯H=US′UN′Σ¯S′Σ¯N′US′UN′H=US′UN′ΣS′+σ2IΣN′+σ2IUS′UN′H,where(20)Σ¯N′=diagλ¯P+1′,λ¯P+2′,…,λ¯M′=diagλP+1′+σ2,λP+2′+σ2,…,λM′+σ2.When the directionθ′ overlaps with the actual sources, that is, θ′=θi,i∈1,2,3,…,P, we can obtain from equation (16) that(21)Σ¯N′=diagλP+1+σ2,λP+2+σ2,…,λM+σ2=Σ¯N.In contrast,(22)Σ¯N′=diagλ¯P+2′,λ¯P+3′,…,λ¯M′≠Σ¯N,as the direction of the virtual source is different from those of any of the actual sources θ′∉θ1,θ2,θ3,…,θP.Now, we have demonstrated that the invariant property of the noise subspace still holds after employing diagonal loading. The noise subspace eigenvalues ofR¯′ remain consistent with those of R¯ when θ′=θi,i∈1,2,3,…,P. Hence, we can still use the IPNS for DOA estimation after diagonal loading.According to the above discussion, we can use the IPNS after diagonal loading to estimate the DOAs of the incident sources. The spatial zero spectrumfθ′ can be expressed as follows:(23)fθ′=∑j=P+1Mλ¯j′−λ¯j.The DOAs of the impinging sources are thoseθ^ for which 1/fθ′ reaches their maxima, which can be expressed as follows:(24)θ^=argmaxθs≤θ′≤θe1fθ′,where θs and θe denote the start angle and end angle of the spectral search, respectively.To implement the improved algorithm, the detailed steps are listed in Algorithm1. Generally, the covariance matrix R is not available in practical applications, and the finite sample covariance matrix R^ is used instead of R.Algorithm 1: Pseudocode of the proposed algorithm.
(i)
Step 1. Obtain the output of the uniform linear antenna array xt from equation (1) and compute the sample covariance matrix R^ as in equation (3).(ii)
Step 2. Perform the EVD of R^ to obtain the noise subspace eigenvalue matrix Σ^Ν=diagλ^P+1,λ^P+2,…,λ^M.(iii)
Step 3. Apply diagonal loading on R^ to obtain R¯^=R^+σ2I through equation (14). σ2 is equal to λ^1, which is the largest eigenvalue of the signal subspace. The noise subspace eigenvalues of R¯^ are Σ¯^Ν=diagλ¯^P+1,λ¯^P+2,…,λ¯^M.(iv)
Step 4. Suppose that the virtual source impinges on the antenna array from direction θ′; then, obtain the matrix R¯′^=R¯^+kaθ′aθ′H as in equation (18), where k=trR¯^/M, tr⋅ denotes the trace of R¯^ and M is the number of array elements.(v)
Step 5. Perform the EVD of R¯′^ to obtain the eigenvalues in descending order as λ¯^1′≥⋯≥λ¯^P′≥λ¯^P+1′≥⋯≥λ¯^M′.(vi)
Step 6. For each θ′, compute the spatial zero spectrum fθ′=∑j=P+1Mλ¯^j′−λ¯^j through equation (23). The DOAs of the sources are P sets of θ′ that enable 1/fθ′ to reach P maximum values, as in equation (24).


### 3.2. Complexity Analysis of the Improved Algorithm

The computational complexity of the algorithm has a direct relationship with its time consumption. LetJ denote the number of spectral points. We must compute aHθUN2 for each spectral point in the MUSIC method, where ⋅2 is the Frobenius norm. The dimensionality of UN is M-P; then, the computational complexity of the MUSIC method in the spectral search step is JM−PM+1 flops [31, 32]. To obtain the noise subspace eigenvector UN, the EVD of the covariance matrix R is needed. If the FSD method in [33] is employed for the EVD of R, the complexity is given by M2P+2 flops. Thus, the total computational complexity of the classical MUSIC method is given by the following:(25)CMUSIC=JM−PM+1+M2P+2flops.However, the EVD ofR′ is required in each spectral point for the IPNS method. Thus, the computational complexity of the IPNS method in the spectral search step is JM2P+2 flops. Before the spectral search step, JM2 flops are needed for computing kaθ′aθ′H to obtain R′. The total computational complexity of the IPNS method is as follows:(26)CIPNS=JM2P+3+M2P+2flops.For the proposed algorithm, diagonal loading requires anotherM flop. The total computational complexity of the improved algorithm is as follows:(27)CIMPROVED=JM2P+3+M2P+2+Mflops.


## 3.1. Improved Algorithm Based on the Invariant Property of Noise Subspace

After introducing a virtual source, the noise subspace remains the same as before only when the DOA of the virtual source overlaps with the actual sources. This property, called the invariant property of noise subspace, is demonstrated in [27] and can be used to estimate the DOAs of actual sources. Diagonal loading is a widespread approach to improve robustness against mismatch errors [29, 30]. We find that diagonal loading is also effective in improving the resolution of the IPNS method. Thus, we propose an improved algorithm based on the IPNS. Diagonal loading effectively reduces the power differences between the strong and weak sources, thereby increasing the resolution of the improved algorithm. Next, we certify that the IPNS still holds after diagonal loading so we can use it for DOA estimation.The eigenvalue decomposition (EVD) of the covariance matrixR is given by the following:(4)R=UΣUH=USUNΣSΣNUSUNH,where U=USUN and Σ=diagΣS,ΣN; diag⋅ denotes the diagonal matrix; and US, UN, ΣS, and ΣN are the eigenvectors and eigenvalues of the signal subspace and noise subspace, respectively. The space spanned by the impinging sources and the noise are orthogonal; that is, spanUS⊥spanUN, where span⋅ denotes the spanned space. ΣS is a diagonal matrix ΣS=diagλ1,λ2,…,λP, where(5)λ1≥λ2≥⋯≥λP,are the eigenvalues of the signal subspace. ΣN is also a diagonal matrix ΣN=diagλP+1,λP+2,…,λM, where(6)λP+1=λP+2=⋯=λM=σn2,are the eigenvalues of the noise subspace. Suppose that a virtual source impinges on the receiving antenna array from the direction θ′. The virtual source does not physically exist but is mathematically needed. Then, the new covariance matrix including the virtual source can be defined as follows:(7)R′≜R+kaθ′aθ′H,where k > 0 is the power of the virtual source. The EVD of the new covariance matrix R′ can be given by the following:(8)R′=U′Σ′U′H=US′UN′ΣS′ΣN′US′UN′H,where U′=US′UN′ and Σ′=diagΣS′,ΣN′, US′ and UN′ are the eigenvectors of the signal subspace and the noise subspace, respectively, and ΣS′=diagλ1′,λ2′,…,λP′ and ΣN′=diagλP+1′,,λP+2′,…,λM′ are the signal and noise eigenvalues corresponding to the new covariance matrix R′, respectively. Suppose that the direction of the virtual source θ′ overlaps with the ith actual source, that is, θ′=θi,i∈1,2,3,…,P. Define a selective unit vector ei as follows:(9)ei≜0,0,…,1,…,0T,where ei is a unit vector consisting of P elements and all the elements are 0 except the ith element, which is 1. Then, the new covariance matrix R′ can be written as follows:(10)R′=ARsAH+kaθ′aθ′H+σn2I=ARsAH+kAeieiTAH+σn2I=ARS′AH+σn2I,where RS′m,n=Rsm,n+k when m=n=i; otherwise, RS′m,n=Rsm,n. It is evident in equation (10) that the virtual source does not affect the noise power but only increases the power of the ith source when the direction θ′ overlaps with that of the ith actual source. Thus, we can obtain that the noise subspace eigenvalues ΣN′ remain consistent with ΣN when θ′ overlaps with the ith actual source, that is,(11)λP+1=λP+2=⋯=λM=λP+1′=λP+2′=⋯=λM′=σn2.When the direction of the virtual source is different from those of any of the actual sources, that is,θ′∉θ1,θ2,θ3,…,θP, then R′ is as follows:(12)R′=ARsAH+kaθ′aθ′H+σn2I=Aaθ′Rs00kAHaθ′H+σn2I.In this case, the noise subspace eigenvaluesΣN′ are different from ΣN. The eigenvalues of R′ in descending order are as follows:(13)λ1′≥λ2′≥⋯≥λP′≥λP+1′≥λP+2′=⋯=λM′=σn2.The eigenvalues of the noise subspace in equation (13) are ΣN′=diagλP+2′,λP+3′,…,λM′, which are different from ΣN=diagλP+1,λP+2,…,λM.Diagonal loading has been demonstrated to be an effective approach for improving robustness in beamforming. We employ diagonal loading to improve the resolution of the IPNS method. After employing diagonal loading, the covariance matrix can be expressed as follows:(14)R¯=R+σ2I,where σ2 is the coefficient of the diagonal loading. After the EVD of R¯, we can obtain the following:(15)R¯=UΣ¯UH=USUNΣ¯SΣ¯NUSUNH=UΣUH+σ2I=UΣ+σ2IUH,where Σ¯=diagΣ¯S,Σ¯N. It is apparent in equation (15) that the diagonal matrix Σ¯N of the noise subspace can be written as follows:(16)Σ¯N=ΣN+σ2I=diagλP+1+σ2,…,λM+σ2.After introducing the virtual source, we define the covariance matrixR¯′ based on R¯ as follows:(17)R¯′≜R¯+kaθ′aθ′H=R+σ2I+kaθaθ′H.Similarly, it can also be certified that the virtual source does not affect the noise subspace ofR¯′ compared to that of R¯ if the direction θ′ overlaps with the actual sources. Assume that the virtual source impinges on the receiving antenna array from the direction θ′. Then, R¯′ can be rewritten as follows:(18)R¯′=R+σ2I+kaθ′aθ′H=ARsAH+kAeieiTAH+σn2I+σ2I=ARs′AH+σn2+σ2I=R′+σ2I,where R′=ARS′AH+σn2I. The EVD of R¯′ is given by the following:(19)R¯′=U′¯Σ′¯U′¯H=US′UN′Σ¯S′Σ¯N′US′UN′H=US′UN′ΣS′+σ2IΣN′+σ2IUS′UN′H,where(20)Σ¯N′=diagλ¯P+1′,λ¯P+2′,…,λ¯M′=diagλP+1′+σ2,λP+2′+σ2,…,λM′+σ2.When the directionθ′ overlaps with the actual sources, that is, θ′=θi,i∈1,2,3,…,P, we can obtain from equation (16) that(21)Σ¯N′=diagλP+1+σ2,λP+2+σ2,…,λM+σ2=Σ¯N.In contrast,(22)Σ¯N′=diagλ¯P+2′,λ¯P+3′,…,λ¯M′≠Σ¯N,as the direction of the virtual source is different from those of any of the actual sources θ′∉θ1,θ2,θ3,…,θP.Now, we have demonstrated that the invariant property of the noise subspace still holds after employing diagonal loading. The noise subspace eigenvalues ofR¯′ remain consistent with those of R¯ when θ′=θi,i∈1,2,3,…,P. Hence, we can still use the IPNS for DOA estimation after diagonal loading.According to the above discussion, we can use the IPNS after diagonal loading to estimate the DOAs of the incident sources. The spatial zero spectrumfθ′ can be expressed as follows:(23)fθ′=∑j=P+1Mλ¯j′−λ¯j.The DOAs of the impinging sources are thoseθ^ for which 1/fθ′ reaches their maxima, which can be expressed as follows:(24)θ^=argmaxθs≤θ′≤θe1fθ′,where θs and θe denote the start angle and end angle of the spectral search, respectively.To implement the improved algorithm, the detailed steps are listed in Algorithm1. Generally, the covariance matrix R is not available in practical applications, and the finite sample covariance matrix R^ is used instead of R.Algorithm 1: Pseudocode of the proposed algorithm.
(i)
Step 1. Obtain the output of the uniform linear antenna array xt from equation (1) and compute the sample covariance matrix R^ as in equation (3).(ii)
Step 2. Perform the EVD of R^ to obtain the noise subspace eigenvalue matrix Σ^Ν=diagλ^P+1,λ^P+2,…,λ^M.(iii)
Step 3. Apply diagonal loading on R^ to obtain R¯^=R^+σ2I through equation (14). σ2 is equal to λ^1, which is the largest eigenvalue of the signal subspace. The noise subspace eigenvalues of R¯^ are Σ¯^Ν=diagλ¯^P+1,λ¯^P+2,…,λ¯^M.(iv)
Step 4. Suppose that the virtual source impinges on the antenna array from direction θ′; then, obtain the matrix R¯′^=R¯^+kaθ′aθ′H as in equation (18), where k=trR¯^/M, tr⋅ denotes the trace of R¯^ and M is the number of array elements.(v)
Step 5. Perform the EVD of R¯′^ to obtain the eigenvalues in descending order as λ¯^1′≥⋯≥λ¯^P′≥λ¯^P+1′≥⋯≥λ¯^M′.(vi)
Step 6. For each θ′, compute the spatial zero spectrum fθ′=∑j=P+1Mλ¯^j′−λ¯^j through equation (23). The DOAs of the sources are P sets of θ′ that enable 1/fθ′ to reach P maximum values, as in equation (24).


## 3.2. Complexity Analysis of the Improved Algorithm

The computational complexity of the algorithm has a direct relationship with its time consumption. LetJ denote the number of spectral points. We must compute aHθUN2 for each spectral point in the MUSIC method, where ⋅2 is the Frobenius norm. The dimensionality of UN is M-P; then, the computational complexity of the MUSIC method in the spectral search step is JM−PM+1 flops [31, 32]. To obtain the noise subspace eigenvector UN, the EVD of the covariance matrix R is needed. If the FSD method in [33] is employed for the EVD of R, the complexity is given by M2P+2 flops. Thus, the total computational complexity of the classical MUSIC method is given by the following:(25)CMUSIC=JM−PM+1+M2P+2flops.However, the EVD ofR′ is required in each spectral point for the IPNS method. Thus, the computational complexity of the IPNS method in the spectral search step is JM2P+2 flops. Before the spectral search step, JM2 flops are needed for computing kaθ′aθ′H to obtain R′. The total computational complexity of the IPNS method is as follows:(26)CIPNS=JM2P+3+M2P+2flops.For the proposed algorithm, diagonal loading requires anotherM flop. The total computational complexity of the improved algorithm is as follows:(27)CIMPROVED=JM2P+3+M2P+2+Mflops.


## 4. CRB of a Weak Source in the Presence of a Strong Source

The CRB provides the lower bound for the unbiased estimators and provides fundamental lower limits for the variance of the parameter estimator [34–36]. Thus, studies on the CRB of weak sources in the presence of strong sources with different powers can be used as references for unbiased estimators. For multiparameter estimation, the CRB inequality has the following form [37]:(28)varθ^i≥J−1ii,where θ^i is the DOA estimation of the ith source, J is the Fisher information matrix, and ⋅−1 denotes the matrix inverse operation. The ijth element of the Fisher information matrix J is given by the following [34]:(29)Jij=N⋅trR−1⋅∂R∂θiR−1⋅∂R∂θj.For independent incident sources, the input signal covariance matrix is diagonal, that is,Rs=diagρ1,ρ2,…,ρP, where ρi represents the power of the ith source. Hence, the covariance matrix R can be simplified as follows:(30)R=∑i=1PρiaθiaθiH+σn2I.The partial derivative of the array covariance matrixR with respect to θi is as follows:(31)∂R∂θi=ρi∂aθi∂θiaθiH+aθi∂aθiH∂θi,where the partial derivative of the steering vector aθi with respect to θi is(32)∂aθi∂θi=∂∂θi1qθi⋯qθiM−1T,and qθi=e−jk0dsinθi. dqθim−1/dθi is given by the following:(33)dqθim−1dθ=m−1⋅qθim−2⋅qθi⋅−jk0dcosθi=−jm−1k0dcosθi⋅qθim−1,m=1,2,…,M.When the distance between the receiving antenna elements isλ0/2, we can obtain that dqθim−1/dθi=−jm−1πcosθie−jm−1πsinθi,m=1,2,…,M. Substituting equation (33) into equation (32) and then into equation (31), we obtain ∂R/∂θi. Then, we can obtain the ijth element of the Fisher matrix Jij from equation (29). The CRB of θ=θ1,θ2,…,θP can be obtained, that is,(34)θCRB=J−1ii.


## 5. Numerical Simulation

Computer simulations regarding the CRB of the weak source, the spatial spectrum of the proposed method, the resolution probability, the root-mean-square error (RMSE), and the computational complexity are conducted to demonstrate the performance of our improved algorithm. We run 1000 independent Monte Carlo simulations for each experiment to study the probability of resolution and the RMSE of the proposed method, while the other two algorithms, the MUSIC algorithm and IPNS method, are also discussed as comparative simulations. The input signal-to-noise ratio (SNR) is defined asSNR=10log10ρs2/σn2, where ρs2 is the power of the source. In the following simulations, we set source 1 with SNR1 as the weak source and source 2 with SNR2 as the strong source. The interference-to-signal ratio (ISR) is defined as SNR2- SNR1. The ISR variation is the result of increasing SNR2, while SNR1 does not change.


### 5.1. CRB Analysis

As shown in Figure1, suppose two independent sources impinge on the receiving antenna array with M = 9 element numbers, and the power of the noise is σn2=1. The study of the CRB of the weak source calculated by θCRB is conducted in the following simulations. The CRB increase percentage is defined as θCRB−θCRB0/θCRB0×100%, where θCRB0 is the CRB of the source with ISR = 0 dB, which means that the incident two sources have the same power.In the first simulation, we are interested in the change in the CRB of the weak source when the ISR increases, but SNR1 remains the same. The numerical simulation results of the CRB as a function of the ISR under different SNR1 conditions are shown in Figure 2. Figure 2 shows that the CRB of the weak source rises as the ISR increases and eventually reaches a fixed value, and this indicates that the CRB of the weak source deteriorates as the power of the strong interference increases. We can also see that CRB deteriorations exist for each SNR1, meaning that increasing the power of the strong source has a negative effect on the weak source, especially when the SNR1 of the weak source is low.Figure 2
The CRBs of the weak source against the ISR under different SNR1.The increase percentages of the CRB of the weak source against the ISR under different SNR1 are shown in Figure 3. From Figure 3, we can see that the CRB increases percentage for the weak source increases sharply from 0 to 11% when SNR1 = 0 dB, but it only increases from 0 to 5.5% when SNR1 = 8 dB. This demonstrates that the CRB of the weak source is more susceptible to the power of the strong source when its SNR is low. Furthermore, it can be seen from Figure 3 that the CRB increase percentage increases from 11% to 12% as the ISR increases from 10 dB to 20 dB when SNR1 = 0 dB. This means that the influence of the power of the strong source is limited when its power exceeds that of the weak source by too much.Figure 3
The CRB increase percentages for the weak source against the ISR under different SNR1.In the second simulation, we further study the influence caused by the power increase of the strong source on the CRB of the weak source. The CRBs of the weak source as a function of the SNR1 under different ISRs are shown in Figure 4. It can be seen in Figure 4 that the CRB of the weak source increases as the strong source power increases at each SNR1, and this once again demonstrates that the power of the strong source leads to CRB degradation for the weak source. Moreover, this reveals that the estimation performances of weak sources are affected by the power of strong sources.Figure 4
The CRBs of the weak source against SNR1 under different ISRs.The CRB increase percentage for the weak source against SNR1 under different ISRs is shown in Figure 5. It is clearly seen in Figure 5 that the CRB increase percentage changes from 18.5% to 4.4% as SNR1 increases from -10 dB to 10 dB when ISR = 20 dB. The CRB degradation is most apparent when SNR1 is low, and this can also be found when ISR = 2 dB, ISR = 5 dB, and ISR = 10 dB; therefore, the weak source is more susceptible when its SNR is lower and the power of the strong interference is higher.Figure 5
The increase percentage of the CRB of the weak source against the ISR under different SNR1.In the third simulation, the CRB of the weak source against the angular separation between the two sources is studied. We suppose that two independent sources impinge on the antenna array from±θ; then, the angular separation between the two sources is 2θ. When θ varies from 1° to 7.5°, the CRB of the weak source versus the angular separation is shown in Figure 6. Figure 6 indicates that the CRB of source 1 (the weak source) is lowest when the two sources have the same power. The CRB of the weak source rises when the strong interference power increases under each angular separation.Figure 6
The CRBs of the weak source against the angular separation under different ISRs.From the above discussion, we can see that the CRB of the weak source is susceptible and tends to be influenced by the power of the strong source, especially when its SNR is low; this has not been studied before. However, sources with large power differences are very common in complex electromagnetic environments. The DOA estimation for sources with different powers should receive special attention.


### 5.2. Spatial Spectra

In this simulation, we compare the spatial spectra of different algorithms when resolving two and three sources. Suppose that two independent sources impinge on the receiving antenna array (M=9) from 2° and −2°, respectively. Source 1 is the weak source with SNR1 = 2 dB, and source 2 is the strong source with SNR2 = 12 dB. The spatial spectra are shown in Figure 7(a).Figure 7
The spatial spectra of different methods. (a) Two sources from±2° with SNR1 = 2 dB and SNR2 = 12 dB; (b) three sources from −2°, 2° and 6° with SNR1 = 17 dB, SNR2 = 7 and SNR3 = 12.
(a)(b)Suppose that three independent sources impinge on the receiving antenna array (M=9) from −2°, 2° and 6° with SNR1 = 17 dB, SNR2 = 7 dB, and SNR3 = 12 dB, respectively. The spatial spectra of the different algorithms are shown in Figure 7(b). From Figure 7, we can see that only the strong source is readily estimated by the MUSIC algorithm when two sources exist. Three sources cannot even be resolved using the MUSIC algorithm. The spectra of the improved algorithm and IPNS method are similar. As shown in Figure 7(b), three peaks are more distinct for the improved algorithm than those of the IPNS method.


### 5.3. Probability of Resolution

In this section, we focus on the performance of the improved algorithm in terms of resolving two impinging sources. For the first simulation, we suppose that two independent sources (source 1 and source 2) impinge on the receiving antenna array (M=9) from 1.5° and −1.5°, respectively.The probability of resolution as a function of SNR1 is shown in Figure 8. It can be seen in Figure 8 that the improved algorithm has a higher probability of resolution than those of the other two methods. When SNR1 = 1 dB and ISR = 10 dB, the probability of resolution for the improved algorithm exceeds those of the IPNS method and the MUSIC algorithm by 17.8% and 62.0%, respectively. Furthermore, it is 22.4% and 55.1% higher than those of the IPNS method and the MUSIC algorithm when ISR = 0 dB and SNR1 = 5 dB.Figure 8
The probability of resolution against SNR1 under different ISRs.For the second simulation, we suppose that two independent sources impinge on the receiving antenna array (M=9) with SNR1 = 10 dB. The probability of resolution as a function of the angular separation between the input sources under different ISRs is shown in Figure 9, which indicates that the improved algorithm has better resolution performance than the other methods. When the angular separation is 1.4° and ISR = 10 dB, the probability of resolution for the improved algorithm exceeds those of the IPNS method of the MUSIC algorithm by 18.4% and 42.3%, respectively. Moreover, it is 25.4% and 56.7% higher than those of the IPNS method and the MUSIC algorithm when the angular separation is 2.15° and ISR = 0 dB, respectively.Figure 9
The probability of resolution against the angular separation under different ISRs.The probability of resolution for each method is low when SNR1 is low or the angular separation is small, and the proposed method is no exception. Then, the increase in the probability of resolution yielded by the proposed method is not so apparent when compared with the other two methods. Similarly, the increase in the probability of resolution for the proposed method is not so significant when SNR1 is high or the angular separation is large. The probability of resolution for each method is high under such conditions.From the abovementioned discussion, we can see that the improved algorithm has better resolution performance for sources with larger power differences; even for sources with the same power, the proposed algorithm still outperforms the other two methods. From the abovementioned discussion, we can see that diagonal loading is effective in improving the probability of resolution of the DOA estimation method based on the IPNS. Diagonal loading reduces the ratio ofλ^P+1/λ^M, thus increasing the probability of resolution. However, the increase in the probability of resolution generated by diagonal loading is limited. Diagonal loading also decreases the signal-to-noise ratio (SNR).


### 5.4. Root-Mean-Square Error

In this section, we study the RMSE of the proposed method. Suppose that two independent sources impinge on the receiving antenna array (M=9) from 1.5° and −1.5°, respectively. The RMSE of the weak source as a function of SNR1 is shown in Figure 10. As a comparison, the CRB is also provided in Figure 10. We can find in Figures 10(a)–10(c) that the improved algorithm achieves better RMSE performance for the weak source than those of the other methods under each ISR. More specifically, the RMSE of the improved algorithm is similar but slightly better than that of the IPNS method when ISR = 5 dB and ISR = 10 dB, and this is clearest when ISR = 0 dB.Figure 10
The RMSE of source 1 against SNR1 for different algorithms. (a) ISR = 5 dB; (b) ISR = 10 dB; (c) ISR = 0 dB.
(a)(b)(c)The RMSE of the strong source and the CRB as a function of SNR1 are shown in Figure 11. We can see in Figure 11 that the RMSE of the improved algorithm is nearly the same as those of the other two methods when ISR = 5 dB. The RMSE of the improved algorithm is not as good as those of the other two methods when ISR = 10 dB, which means that the MUSIC algorithm has better estimation accuracy for strong sources. However, when the power of the strong source decreases to ISR = 0 dB, the RMSE of the improved algorithm is the lowest. However, the challenge in DOA estimation with respect to sources that have large power differences is the estimation performance of the weak source. Hence, it is objective to conclude that the RMSE of the proposed algorithm is similar to that of the IPNS method, which is even better than that of the MUSIC algorithm.Figure 11
The RMSE of source 2 against SNR1 for different algorithms. (a) ISR = 5 dB; (b) ISR = 10 dB; (c) ISR = 0 dB.
(a)(b)(c)


### 5.5. Computational Complexity

In this section, we study the computational complexity of the proposed method. If a rough grid of 3° and a fine grid of 1° are employed, the spectral points60J1=180°/3° and 180J2=180°/1° should be computed in the spectral search, respectively. We compare the computational complexities of different algorithms and display the results in Table 1. It is evident in Table 1 that the IPNS method has more computational complexity than the MUSIC algorithm; moreover, the computational complexity of the IPNS method grows faster than that of the MUSIC algorithm as M and P increase. However, the improved algorithm that exploits diagonal loading has nearly the same computational complexity as that of the IPNS method.Table 1
Comparison of the computational complexities of different methods.
ItemsM = 5, P = 1, J1 = 60M = 5, P = 1, J2 = 180M = 9, P = 2, J1 = 60M = 9, P = 2, J2 = 180CMUSIC15154395452412924CIPNS6075180752462473224CIMPROVED6080180802463373233CIPNS/CMUSIC4.00994.11265.44305.6657CIPNS/CIMPROVED0.99910.99970.99960.9998


## 5.1. CRB Analysis

As shown in Figure1, suppose two independent sources impinge on the receiving antenna array with M = 9 element numbers, and the power of the noise is σn2=1. The study of the CRB of the weak source calculated by θCRB is conducted in the following simulations. The CRB increase percentage is defined as θCRB−θCRB0/θCRB0×100%, where θCRB0 is the CRB of the source with ISR = 0 dB, which means that the incident two sources have the same power.In the first simulation, we are interested in the change in the CRB of the weak source when the ISR increases, but SNR1 remains the same. The numerical simulation results of the CRB as a function of the ISR under different SNR1 conditions are shown in Figure 2. Figure 2 shows that the CRB of the weak source rises as the ISR increases and eventually reaches a fixed value, and this indicates that the CRB of the weak source deteriorates as the power of the strong interference increases. We can also see that CRB deteriorations exist for each SNR1, meaning that increasing the power of the strong source has a negative effect on the weak source, especially when the SNR1 of the weak source is low.Figure 2
The CRBs of the weak source against the ISR under different SNR1.The increase percentages of the CRB of the weak source against the ISR under different SNR1 are shown in Figure 3. From Figure 3, we can see that the CRB increases percentage for the weak source increases sharply from 0 to 11% when SNR1 = 0 dB, but it only increases from 0 to 5.5% when SNR1 = 8 dB. This demonstrates that the CRB of the weak source is more susceptible to the power of the strong source when its SNR is low. Furthermore, it can be seen from Figure 3 that the CRB increase percentage increases from 11% to 12% as the ISR increases from 10 dB to 20 dB when SNR1 = 0 dB. This means that the influence of the power of the strong source is limited when its power exceeds that of the weak source by too much.Figure 3
The CRB increase percentages for the weak source against the ISR under different SNR1.In the second simulation, we further study the influence caused by the power increase of the strong source on the CRB of the weak source. The CRBs of the weak source as a function of the SNR1 under different ISRs are shown in Figure 4. It can be seen in Figure 4 that the CRB of the weak source increases as the strong source power increases at each SNR1, and this once again demonstrates that the power of the strong source leads to CRB degradation for the weak source. Moreover, this reveals that the estimation performances of weak sources are affected by the power of strong sources.Figure 4
The CRBs of the weak source against SNR1 under different ISRs.The CRB increase percentage for the weak source against SNR1 under different ISRs is shown in Figure 5. It is clearly seen in Figure 5 that the CRB increase percentage changes from 18.5% to 4.4% as SNR1 increases from -10 dB to 10 dB when ISR = 20 dB. The CRB degradation is most apparent when SNR1 is low, and this can also be found when ISR = 2 dB, ISR = 5 dB, and ISR = 10 dB; therefore, the weak source is more susceptible when its SNR is lower and the power of the strong interference is higher.Figure 5
The increase percentage of the CRB of the weak source against the ISR under different SNR1.In the third simulation, the CRB of the weak source against the angular separation between the two sources is studied. We suppose that two independent sources impinge on the antenna array from±θ; then, the angular separation between the two sources is 2θ. When θ varies from 1° to 7.5°, the CRB of the weak source versus the angular separation is shown in Figure 6. Figure 6 indicates that the CRB of source 1 (the weak source) is lowest when the two sources have the same power. The CRB of the weak source rises when the strong interference power increases under each angular separation.Figure 6
The CRBs of the weak source against the angular separation under different ISRs.From the above discussion, we can see that the CRB of the weak source is susceptible and tends to be influenced by the power of the strong source, especially when its SNR is low; this has not been studied before. However, sources with large power differences are very common in complex electromagnetic environments. The DOA estimation for sources with different powers should receive special attention.


## 5.2. Spatial Spectra

In this simulation, we compare the spatial spectra of different algorithms when resolving two and three sources. Suppose that two independent sources impinge on the receiving antenna array (M=9) from 2° and −2°, respectively. Source 1 is the weak source with SNR1 = 2 dB, and source 2 is the strong source with SNR2 = 12 dB. The spatial spectra are shown in Figure 7(a).Figure 7
The spatial spectra of different methods. (a) Two sources from±2° with SNR1 = 2 dB and SNR2 = 12 dB; (b) three sources from −2°, 2° and 6° with SNR1 = 17 dB, SNR2 = 7 and SNR3 = 12.
(a)(b)Suppose that three independent sources impinge on the receiving antenna array (M=9) from −2°, 2° and 6° with SNR1 = 17 dB, SNR2 = 7 dB, and SNR3 = 12 dB, respectively. The spatial spectra of the different algorithms are shown in Figure 7(b). From Figure 7, we can see that only the strong source is readily estimated by the MUSIC algorithm when two sources exist. Three sources cannot even be resolved using the MUSIC algorithm. The spectra of the improved algorithm and IPNS method are similar. As shown in Figure 7(b), three peaks are more distinct for the improved algorithm than those of the IPNS method.


## 5.3. Probability of Resolution

In this section, we focus on the performance of the improved algorithm in terms of resolving two impinging sources. For the first simulation, we suppose that two independent sources (source 1 and source 2) impinge on the receiving antenna array (M=9) from 1.5° and −1.5°, respectively.The probability of resolution as a function of SNR1 is shown in Figure 8. It can be seen in Figure 8 that the improved algorithm has a higher probability of resolution than those of the other two methods. When SNR1 = 1 dB and ISR = 10 dB, the probability of resolution for the improved algorithm exceeds those of the IPNS method and the MUSIC algorithm by 17.8% and 62.0%, respectively. Furthermore, it is 22.4% and 55.1% higher than those of the IPNS method and the MUSIC algorithm when ISR = 0 dB and SNR1 = 5 dB.Figure 8
The probability of resolution against SNR1 under different ISRs.For the second simulation, we suppose that two independent sources impinge on the receiving antenna array (M=9) with SNR1 = 10 dB. The probability of resolution as a function of the angular separation between the input sources under different ISRs is shown in Figure 9, which indicates that the improved algorithm has better resolution performance than the other methods. When the angular separation is 1.4° and ISR = 10 dB, the probability of resolution for the improved algorithm exceeds those of the IPNS method of the MUSIC algorithm by 18.4% and 42.3%, respectively. Moreover, it is 25.4% and 56.7% higher than those of the IPNS method and the MUSIC algorithm when the angular separation is 2.15° and ISR = 0 dB, respectively.Figure 9
The probability of resolution against the angular separation under different ISRs.The probability of resolution for each method is low when SNR1 is low or the angular separation is small, and the proposed method is no exception. Then, the increase in the probability of resolution yielded by the proposed method is not so apparent when compared with the other two methods. Similarly, the increase in the probability of resolution for the proposed method is not so significant when SNR1 is high or the angular separation is large. The probability of resolution for each method is high under such conditions.From the abovementioned discussion, we can see that the improved algorithm has better resolution performance for sources with larger power differences; even for sources with the same power, the proposed algorithm still outperforms the other two methods. From the abovementioned discussion, we can see that diagonal loading is effective in improving the probability of resolution of the DOA estimation method based on the IPNS. Diagonal loading reduces the ratio ofλ^P+1/λ^M, thus increasing the probability of resolution. However, the increase in the probability of resolution generated by diagonal loading is limited. Diagonal loading also decreases the signal-to-noise ratio (SNR).


## 5.4. Root-Mean-Square Error

In this section, we study the RMSE of the proposed method. Suppose that two independent sources impinge on the receiving antenna array (M=9) from 1.5° and −1.5°, respectively. The RMSE of the weak source as a function of SNR1 is shown in Figure 10. As a comparison, the CRB is also provided in Figure 10. We can find in Figures 10(a)–10(c) that the improved algorithm achieves better RMSE performance for the weak source than those of the other methods under each ISR. More specifically, the RMSE of the improved algorithm is similar but slightly better than that of the IPNS method when ISR = 5 dB and ISR = 10 dB, and this is clearest when ISR = 0 dB.Figure 10
The RMSE of source 1 against SNR1 for different algorithms. (a) ISR = 5 dB; (b) ISR = 10 dB; (c) ISR = 0 dB.
(a)(b)(c)The RMSE of the strong source and the CRB as a function of SNR1 are shown in Figure 11. We can see in Figure 11 that the RMSE of the improved algorithm is nearly the same as those of the other two methods when ISR = 5 dB. The RMSE of the improved algorithm is not as good as those of the other two methods when ISR = 10 dB, which means that the MUSIC algorithm has better estimation accuracy for strong sources. However, when the power of the strong source decreases to ISR = 0 dB, the RMSE of the improved algorithm is the lowest. However, the challenge in DOA estimation with respect to sources that have large power differences is the estimation performance of the weak source. Hence, it is objective to conclude that the RMSE of the proposed algorithm is similar to that of the IPNS method, which is even better than that of the MUSIC algorithm.Figure 11
The RMSE of source 2 against SNR1 for different algorithms. (a) ISR = 5 dB; (b) ISR = 10 dB; (c) ISR = 0 dB.
(a)(b)(c)


## 5.5. Computational Complexity

In this section, we study the computational complexity of the proposed method. If a rough grid of 3° and a fine grid of 1° are employed, the spectral points60J1=180°/3° and 180J2=180°/1° should be computed in the spectral search, respectively. We compare the computational complexities of different algorithms and display the results in Table 1. It is evident in Table 1 that the IPNS method has more computational complexity than the MUSIC algorithm; moreover, the computational complexity of the IPNS method grows faster than that of the MUSIC algorithm as M and P increase. However, the improved algorithm that exploits diagonal loading has nearly the same computational complexity as that of the IPNS method.Table 1
Comparison of the computational complexities of different methods.
ItemsM = 5, P = 1, J1 = 60M = 5, P = 1, J2 = 180M = 9, P = 2, J1 = 60M = 9, P = 2, J2 = 180CMUSIC15154395452412924CIPNS6075180752462473224CIMPROVED6080180802463373233CIPNS/CMUSIC4.00994.11265.44305.6657CIPNS/CIMPROVED0.99910.99970.99960.9998


## 6. Conclusions

In this paper, we propose an improved algorithm to increase the resolution performance with regard to resolving sources with large power differences and further study the Cramer-Rao bound of the weak source in the presence of the strong source. We demonstrate that the IPNS still holds after diagonal loading. Simulation results show that the proposed algorithm has a better probability of resolution than the IPNS and MUSIC algorithms. Moreover, the proposed method does not result in increased computational complexity and maintains the estimation accuracy. We also find that the Cramer-Rao bound of the weak source suffers from degradation as the power of the strong source increases, especially when the SNR of the weak source is low. Our work provides a new way to estimate sources with different powers, and this is valuable in noncooperative situations, such as passive radar applications.



---

*Source: 8862789-2021-03-12.xml*