# Joint Angle and Frequency Estimation in Linear Arrays Based on Covariance Reconstruction and ESPRIT

Authors: Shihong Chen; Qingchang Tao; Zhongtian Yang; Xudong Wang; Sijia Liu; Wei Xu
Journal: Mathematical Problems in Engineering
 (2021)

---

## Abstract

Joint angle and frequency estimation, one of the key technologies in wireless communication and radar science, has been extensively studied by scholars. For linear arrays, this paper proposes a joint angle and frequency estimation method based on covariance reconstruction and the estimation of signal parameters via rotational invariance techniques (CR-ESPRIT). We first use the received conjugate signal to reconstruct a covariance matrix. Then, we use the least squares-ESPRIT (LS-ESPRIT) algorithm to estimate the desired frequencies. Finally, we estimate the angles according to the reconstructed matrix. The proposed method can estimate signal parameters via automatic pairing and without an additional parameter pairing process under the condition of a uniform or a nonuniform array. Moreover, this method has high estimation accuracy, excellent and stable anti-noise performance, and strong algorithmic robustness. Through a computer simulation analysis, we can confirm the reliability and validity of the proposed parameter estimation method. A comparison with other methods further proves the performance advantages of the developed method. The method in this paper can be easily applied to many signal processing contexts, such as electronic reconnaissance and wireless communication.


---

## Body

## 1. Introduction

The joint angle and frequency estimation of received signals submerged in Gaussian white noise has important applications in wireless communication [1], audio and speech signal processing [2], and other fields [3, 4]. For example, in a wireless communication system, accurate and robust joint angle and frequency estimation can help provide better channel information, thereby improving the link quality and anti-interference ability of the system [1]. Especially in electronic reconnaissance [5–8], we often use the operating frequencies and directions of arrival (DOAs) [9–13] of noncooperative radar radiation source signals to describe the main parameters of radar signal characteristics [14–16]. Therefore, to effectively obtain the parameters of noncooperative radar source signals, it is necessary to study a joint DOA and frequency estimation method for such signals submerged in Gaussian white noise.Regarding the joint DOA and frequency estimation of noisy signals, researchers worldwide have proposed various methods [16–22]. In 1986, Schmidt [17] proposed the multiple signal classification (MUSIC) algorithm for parameter estimation. Although the algorithm has good estimation performance, it has high computational complexity since it needs to search for spectral peaks to obtain the estimated values. Lemma et al. [18] presented a joint angle and frequency estimation method based on the multidimensional estimation of signal parameters via rotational invariance techniques (ESPRIT). Nevertheless, this algorithm has low parameter estimation accuracy under low signal-to-noise ratios (SNRs). To effectively improve the accuracy of estimated DOA and frequency results, in 2010, Wang proposed a joint angle and frequency estimation technique using multiple-delay outputs (MDJAFE) [16] based on the ESPRIT algorithm. However, this method cannot realize automatic parameter pairing when performing the joint estimation of signal parameters. Since the propagator method (PM) shows good performance in parameter estimation, it has attracted the attention of scholars. Sun et al. [19] proposed a joint DOA and frequency estimation based on the improved PM. Although the complexity of the algorithm is low and it can realize the automatic pairing of DOA and frequency parameter estimations, its parameter estimation accuracy is not high. Wang et al. [20] proposed an improved ESPRIT algorithm using the multidelay output of a uniform linear antenna (ULA). Although the algorithm’s complexity is greatly reduced, this method is greatly affected by noise, and the estimation accuracy of this method is still very limited when the SNR is low. Based on the extended orthogonal matching pursuit (EOMP) algorithm, Gao et al. [21] proposed an approach to jointly estimate DOAs and frequency, whereas this method has high computational complexity. Xu et al. [22] proposed a joint DOA and narrowband source carrier frequency estimation method based on parallel factor (PARAFAC) analysis. The computational complexity of this method is relatively high, and the hardware cost is also high.Due to the wide range of possible SNRs, frequency and DOA estimation algorithms have unstable anti-noise performance and limited estimation accuracy. We propose a method for the joint DOA and frequency estimation of signals submerged in Gaussian white noise. The algorithm involves a three-step estimation procedure. First, we preprocess the received signal. Second, we use the least squares-ESPRIT (LS-ESPRIT) algorithm to estimate the frequency parameters of the signal. Finally, according to the unique relationship between the signal angle and its frequency, we estimate the DOAs. Computer simulations and comparisons with other methods prove the excellent performance of the proposed method.The main contributions of our work can be summarized as follows:(1)
We improve upon the estimation process in [20]. Under the condition of a uniform or a nonuniform array, the method proposed in this paper can estimate the required parameters by performing automatic pairing without an additional parameter pairing process. Moreover, this method has good estimation accuracy, stable anti-noise performance, and robustness. Therefore, the method proposed in this paper is more suitable than other approaches for the parameter estimation of noncooperative radar radiation sources in an external field, which usually contains a complex electromagnetic environment.(2)
This paper proposes a joint angle and frequency estimation method based on covariance reconstruction and ESPRIT (CR-ESPRIT). Within the SNR range from -15 dB to 15 dB (step: 2 dB), its performance is better than that of the PM, the covariance reconstruction and propagator method (CR-PM), the ESPRIT method [16], and the improved ESPRIT method [20].The remainder of this paper is structured as follows. The materials and methods are presented in Section2; Section 3 contains the results and a discussion, and Section 4 is the summary of the paper.Notations.•H, •∗, •−1, and •+ denote the conjugate transpose, complex conjugation, inverse, and Moore–Penrose inverse (pseudoinverse) operations, respectively. Matrices and vectors are represented by boldfaced capital letters and lowercase letters, respectively.


## 2. Materials and Methods

### 2.1. Signal Model

Consider an antenna array that consists ofM array elements arranged in a straight line at equal distances, where the distance between each pair of array elements is d [23]. We suppose that there exist K (K < M) far-field source narrowband signals (the center frequency is fk), which are incident on the antenna array. Therefore, we can regard the signals as plane waves when they reach the array. Then, we can express the received signal of the mth antenna as follows [24]:(1)ymt=∑k=1Kexp−j2πm−1dfksinθkcskt+nmt,m=1,2,…,M,where skt is the kth incident far-field source signal, c is the speed of light (m/s), θk and fk are the DOA and frequency of the kth signal, respectively, and nmt is the zero-mean additive white Gaussian noise on the mth antenna. We can express the output signal of the linear array as(2)Y0=y1ny2n,…,yMnT,n=1,2,…,N.We assume that the signal is uniformly sampled by a period that conforms to the Nyquist sampling rate and that the number of snapshots isN. Therefore, we can transform the signal model studied in this paper into a joint DOA and frequency estimation model for multiple source signals, where N sampling points are obtained for each source signal.We assume that the number of signal sourcesK is known; thus, we can rewrite output state vector (2) in the following matrix form:(3)Y0=AS+N0,where S=s1,s2,…,sKT∈CK×N, N0=n1,n2,…,nMT∈CM×N, and(4)A=11…1exp−jα1exp−jα2…exp−jαK…………exp−jM−1α1exp−jM−1α2…exp−jM−1αK.In equation (4), αk=2πdfksinθk/c,k=1,…,K. To realize the joint DOA and frequency estimation model, we take (P-1) delays [25, 26] for the signal received from the antenna arrays shown in Figure 1. In addition, we set 0<P−1τ<1/maxfk. Therefore, we can obtain the delay signal with the delay value τ as(5)ymt−τ=∑k=1Kexp−j2πm−1dfksinθkcskt−τ+nm′t=∑k=1Kexp−j2πm−1dfksinθkcsktexp−j2πfkτ+nm′t.Figure 1
Received signals with multilevel delays.We can transform equation (5) into the following form:(6)Y1=AΦS+N1,where βk=2πfkτ,k=1,2,…,K and Φ=diagexp−jβ1,exp−jβ2,…,exp−jβK.When the delay value ispτ, we can express the delay signal as(7)ymt−pτ=∑k=1Kexp−j2πm−1dfksinθkcskt−pτ+nm′t=∑k=1Kexp−j2πm−1dfksinθkcsktexp−j2πfkτp+nm′t.Then, we can also express equation (7) as(8)Yp=AΦpS+Np,p=0,1,2,…,P−1.After reorganizing the equations, we can obtain the following expression:(9)Y=Y0Y1…YP−1=AAΦ…AΦP−1S+N0N1…NP−1.


### 2.2. The Proposed Method

In this paper, inspired by the improved ESPRIT method [20], we propose a joint angle and frequency estimation method based on CR-ESPRIT. In a real space, the improved ESPRIT method is not suitable for complex electromagnetic environments. Due to the noncooperative characteristics of radiation sources, we generally believe that there is no prior information available regarding the parameters. Moreover, in a complex and harsh electromagnetic environment, the detected radiation source signals are very weak. Therefore, in a situation with a low SNR, the developed method not only needs to distinguish useful signals and noise effectively but also needs to have good estimation performance, noise immunity and robustness. Additionally, it also needs to have the ability to automatically pair the relevant parameters without an additional parameter pairing process under the condition of a uniform or a nonuniform array.We first preprocess the received signal in Section2.2.1. Second, we use the LS-ESPRIT algorithm to estimate the frequency parameters of the received signal in Section 2.2.2. Third, according to the relationship between the DOA and frequency in the signal model, we reconstruct the received signal, and then we estimate the DOAs in Section 2.2.3. In Section 2.2.4, we provide the detailed steps of the proposed method. Finally, we provide the detailed steps of the proposed method under the condition of a nonuniform array in Section 2.2.5.


#### 2.2.1. The Preprocessing Procedure

First, we obtain the covariance matrixRY=YYH of the received signal. To make full use of the conjugate information contained in the received signal, we define the permutation matrix J [27]:(10)J=0…01…0100………10…0PM×PM.Therefore, we can construct the following matrix:(11)RJ=JY∗Y∗HJH.We add the covariance matrixRY and RJ from equation (11), and then we average them. The form of the obtained covariance matrix is shown as follows:(12)R=RY+RJ2.Through analysis, we can obtain that the new total covariance matrixR is a Hermitian matrix (PM × PM) [28]. Therefore, we can apply eigenvalue decomposition, and then we can reconstruct the signal subspace Ess. In a no-noise situation, Ess can be approximately expressed as(13)Ess=AAΦ…AΦP−1F,where F is a full-rank matrix with K×K dimensions.Remark 1.
As mentioned earlier, the new total covariance matrixR is already a Hermitian matrix. According to Hermitian matrix characteristics, we assume that the diagonal matrix of the eigenvalues of R is G. Then, there exists a unitary matrix U, which assures RU = UG. Therefore, we can treat R as the unitary matrix U by using this correlation feature to further reduce the complexity of the proposed method and then propose a much lower complexity method.


#### 2.2.2. Frequency Estimation

We define the following parameters:(14)E1=AAΦ…AΦP−2F,(15)E2=AΦAΦ2…AΦP−1F.Therefore, equations (14) and (15) have the following relationship:(16)E2=AΦAΦ2…AΦP−1F=AAΦ1…AΦP−2FF−1ΦF=E1F−1ΦF.LetΨ=F−1ΦF and Ψ=E1+E2. According to the LS-ESPRIT algorithm, we can estimate Φ by the eigenvalue decomposition of Ψ, and we can also estimate the matrix F−1 by the eigenvector of Φ. In a no-noise situation, we define(17)F^−1=F−1Θ,Φ^=Θ−1ΦΘ,where Θ is a fuzzy column matrix. Since Ψ and Φ have the same eigenvalues, we can obtain the eigenvalues λkk=1,2,…,K from matrix Ψ. As shown in equation (6) βk=2πfkτ,k=1,2,…,K, it is obvious that we can estimate the frequency parameter f^k,k=1,2,…,K:(18)f^k=12πτangleλk.


#### 2.2.3. DOA Estimation

AΦP−1 has the following expression:(19)AΦP−1=exp−jP−1β1exp−jP−1β2…exp−jP−1βKexp−jα1exp−jP−1β1exp−jα2exp−jP−1β2…exp−jαKexp−jP−1βK…………exp−jM−1α1exp−jP−1β1exp−jM−1α2exp−jP−1β2…exp−jM−1αKexp−jP−1βK.According to the above estimationF^−1, we can define the following expression by reconstructing equation (13):(20)EQ=BBΤ…BTM−1Θ,where(21)B=11…1exp−jβ1exp−jβ2…exp−jβK…………exp−jP−1β1exp−jP−1β2…exp−jP−1βK,Τ=diagexp−jα1,exp−jα2,…,exp−jαK∈CK×K,αk=2πdfksinθkc,k=1,2,…,K.According to reconstructed equation (20), we can use the method described below to estimate the DOA.We define the following matrices:(22)EQ1=BBΤ…BΤM−2Θ,EQ2=BTBT2…BTM−1Θ.We can also define a matrixD since D=EQ1+EQ2. Then, according to the definitions of EQ1 and EQ2, we can express D in a no-noise situation:(23)D=Θ−1TΘ.Therefore, we can take the diagonal elements ofD, and then, we can obtain ϖkk=1,2,…,K, where αk=2πdfksinθk/c(k=1,2,…,K, to obtain the estimation of the DOA:(24)θ^k=arcsinc2πf^kdangleϖk,k=1,2,…,K.


#### 2.2.4. The Steps of the Proposed Method

Thus far, we have given the complete process for automatically pairing DOA and frequency estimations in a linear array. The main steps required to implement the method proposed in this paper are as follows:(i)
Step 1: according to permutation matrixJ and equation (12), we reconstruct the covariance matrix R.(ii)
Step 2: we apply eigenvalue decomposition toR, and then we reconstruct the signal subspace Ess. According to equations (14) and (15), we construct matrices E1 and E2, respectively.(iii)
Step 3: we use equationΨ=E1+E2 for eigenvalue decomposition to obtain F−1 and Φ^. Finally, we estimate the frequency parameter f^k according to equation (18).(iv)
Step 4: we can obtain matrixEQ according to the reconstruction of Ess in equation (13). Then, we can also construct matrices EQ1 and EQ2.(v)
Step 5: we calculateD=EQ1+EQ2 to obtain matrix D. Finally, we estimate the DOA parameter θ^k according to equation (24).


#### 2.2.5. The Condition of a Nonuniform Array

In this section, we first present the method proposed in this paper when the distances between the array elements are not equal. Then, we present the main steps for implementing the method in the case of a nonuniform array.We assume that the first elementd1 = 0 and that the distance between the mth element and the first element is dm. Then, we can transform equation (4) into the following form:(25)A1=11…1exp−jd2η1exp−jd2η2…exp−jd2ηK…………exp−jdMη1exp−jdMη2…exp−jdMηK,where ηk=2πfksinθk/ck=1,2,…,K. At the same time, equation (19) undergoes the following transformation:(26)A1ΦP−1=exp−jP−1β1exp−jP−1β2…exp−jP−1βKexp−jd2η1exp−jP−1β1exp−jd2η2exp−jP−1β2…exp−jd2ηKexp−jP−1βK…………exp−jdMη1exp−jP−1β1exp−jdMη2exp−jP−1β2…exp−jdMηKexp−jP−1βK.Similarly, we can reconstruct equation (13) and define the following expression:(27)EQQ=BJ1BJ2…BJMΘ,where Jm=diagexp−jdmη1,exp−jdmη2,…,exp−jdmηK∈CK×K, m=1,2,…,M.We define the following matrix:(28)EQQ1=BJ1BJ2…BJM−1Θ=EQQ11EQQ12…EQQ1M−1,EQQ2=BJ2BJ3…BJMΘ=EQQ22EQQ23…EQQ2M.We also define the matrixQm, where Qm=EQQ1m−1+EQQ2m, m=2,…,M. Therefore, we can take the diagonal elements of Qm and then obtain(29)νm=diagexp−jdm−dm−1η1,exp−jdm−dm−1η2,…,exp−jdm−dm−1ηK∈CK×K,m=2,…,M.We sort the diagonal elements and then define the following matrix:(30)V=11…1exp−jd2−d1η1exp−jd2−d1η2…exp−jd2−d1ηK…………exp−jdm−dm−1η1exp−jdm−dm−1η2…exp−jdm−dm−1ηK.According to equation (30), we can obtain the estimation of the DOA:(31)θ^k=1M−1∑m=2Marcsinc2πf^kdm−dm−1angleVmk,k=1,2,…,K.The main steps for implementing the method in this paper under the condition of a nonuniform array are as follows:(i)
Step 1: according to equations (10), (12), and (25), we reconstruct the new covariance matrix.(ii)
Step 2: we apply eigenvalue decomposition to the new covariance matrix, and then we reconstruct the signal subspaceEss. According to equations (14) and (15), we construct matrices E1 and E2, respectively.(iii)
Step 3: we use equationΨ=E1+E2 to perform eigenvalue decomposition and obtain F−1 and Φ^. Finally, we estimate the frequency parameter f^k according to equation (18).(iv)
Step 4: we can obtain the matrixEQQ according to equation (27). Then, we can also construct matrices EQQ1 and EQQ2.(v)
Step 5: we can obtain matrixV according to equation (30). Finally, we estimate the DOA parameter θ^k according to equation (31).


## 2.1. Signal Model

Consider an antenna array that consists ofM array elements arranged in a straight line at equal distances, where the distance between each pair of array elements is d [23]. We suppose that there exist K (K < M) far-field source narrowband signals (the center frequency is fk), which are incident on the antenna array. Therefore, we can regard the signals as plane waves when they reach the array. Then, we can express the received signal of the mth antenna as follows [24]:(1)ymt=∑k=1Kexp−j2πm−1dfksinθkcskt+nmt,m=1,2,…,M,where skt is the kth incident far-field source signal, c is the speed of light (m/s), θk and fk are the DOA and frequency of the kth signal, respectively, and nmt is the zero-mean additive white Gaussian noise on the mth antenna. We can express the output signal of the linear array as(2)Y0=y1ny2n,…,yMnT,n=1,2,…,N.We assume that the signal is uniformly sampled by a period that conforms to the Nyquist sampling rate and that the number of snapshots isN. Therefore, we can transform the signal model studied in this paper into a joint DOA and frequency estimation model for multiple source signals, where N sampling points are obtained for each source signal.We assume that the number of signal sourcesK is known; thus, we can rewrite output state vector (2) in the following matrix form:(3)Y0=AS+N0,where S=s1,s2,…,sKT∈CK×N, N0=n1,n2,…,nMT∈CM×N, and(4)A=11…1exp−jα1exp−jα2…exp−jαK…………exp−jM−1α1exp−jM−1α2…exp−jM−1αK.In equation (4), αk=2πdfksinθk/c,k=1,…,K. To realize the joint DOA and frequency estimation model, we take (P-1) delays [25, 26] for the signal received from the antenna arrays shown in Figure 1. In addition, we set 0<P−1τ<1/maxfk. Therefore, we can obtain the delay signal with the delay value τ as(5)ymt−τ=∑k=1Kexp−j2πm−1dfksinθkcskt−τ+nm′t=∑k=1Kexp−j2πm−1dfksinθkcsktexp−j2πfkτ+nm′t.Figure 1
Received signals with multilevel delays.We can transform equation (5) into the following form:(6)Y1=AΦS+N1,where βk=2πfkτ,k=1,2,…,K and Φ=diagexp−jβ1,exp−jβ2,…,exp−jβK.When the delay value ispτ, we can express the delay signal as(7)ymt−pτ=∑k=1Kexp−j2πm−1dfksinθkcskt−pτ+nm′t=∑k=1Kexp−j2πm−1dfksinθkcsktexp−j2πfkτp+nm′t.Then, we can also express equation (7) as(8)Yp=AΦpS+Np,p=0,1,2,…,P−1.After reorganizing the equations, we can obtain the following expression:(9)Y=Y0Y1…YP−1=AAΦ…AΦP−1S+N0N1…NP−1.


## 2.2. The Proposed Method

In this paper, inspired by the improved ESPRIT method [20], we propose a joint angle and frequency estimation method based on CR-ESPRIT. In a real space, the improved ESPRIT method is not suitable for complex electromagnetic environments. Due to the noncooperative characteristics of radiation sources, we generally believe that there is no prior information available regarding the parameters. Moreover, in a complex and harsh electromagnetic environment, the detected radiation source signals are very weak. Therefore, in a situation with a low SNR, the developed method not only needs to distinguish useful signals and noise effectively but also needs to have good estimation performance, noise immunity and robustness. Additionally, it also needs to have the ability to automatically pair the relevant parameters without an additional parameter pairing process under the condition of a uniform or a nonuniform array.We first preprocess the received signal in Section2.2.1. Second, we use the LS-ESPRIT algorithm to estimate the frequency parameters of the received signal in Section 2.2.2. Third, according to the relationship between the DOA and frequency in the signal model, we reconstruct the received signal, and then we estimate the DOAs in Section 2.2.3. In Section 2.2.4, we provide the detailed steps of the proposed method. Finally, we provide the detailed steps of the proposed method under the condition of a nonuniform array in Section 2.2.5.


### 2.2.1. The Preprocessing Procedure

First, we obtain the covariance matrixRY=YYH of the received signal. To make full use of the conjugate information contained in the received signal, we define the permutation matrix J [27]:(10)J=0…01…0100………10…0PM×PM.Therefore, we can construct the following matrix:(11)RJ=JY∗Y∗HJH.We add the covariance matrixRY and RJ from equation (11), and then we average them. The form of the obtained covariance matrix is shown as follows:(12)R=RY+RJ2.Through analysis, we can obtain that the new total covariance matrixR is a Hermitian matrix (PM × PM) [28]. Therefore, we can apply eigenvalue decomposition, and then we can reconstruct the signal subspace Ess. In a no-noise situation, Ess can be approximately expressed as(13)Ess=AAΦ…AΦP−1F,where F is a full-rank matrix with K×K dimensions.Remark 1.
As mentioned earlier, the new total covariance matrixR is already a Hermitian matrix. According to Hermitian matrix characteristics, we assume that the diagonal matrix of the eigenvalues of R is G. Then, there exists a unitary matrix U, which assures RU = UG. Therefore, we can treat R as the unitary matrix U by using this correlation feature to further reduce the complexity of the proposed method and then propose a much lower complexity method.


### 2.2.2. Frequency Estimation

We define the following parameters:(14)E1=AAΦ…AΦP−2F,(15)E2=AΦAΦ2…AΦP−1F.Therefore, equations (14) and (15) have the following relationship:(16)E2=AΦAΦ2…AΦP−1F=AAΦ1…AΦP−2FF−1ΦF=E1F−1ΦF.LetΨ=F−1ΦF and Ψ=E1+E2. According to the LS-ESPRIT algorithm, we can estimate Φ by the eigenvalue decomposition of Ψ, and we can also estimate the matrix F−1 by the eigenvector of Φ. In a no-noise situation, we define(17)F^−1=F−1Θ,Φ^=Θ−1ΦΘ,where Θ is a fuzzy column matrix. Since Ψ and Φ have the same eigenvalues, we can obtain the eigenvalues λkk=1,2,…,K from matrix Ψ. As shown in equation (6) βk=2πfkτ,k=1,2,…,K, it is obvious that we can estimate the frequency parameter f^k,k=1,2,…,K:(18)f^k=12πτangleλk.


### 2.2.3. DOA Estimation

AΦP−1 has the following expression:(19)AΦP−1=exp−jP−1β1exp−jP−1β2…exp−jP−1βKexp−jα1exp−jP−1β1exp−jα2exp−jP−1β2…exp−jαKexp−jP−1βK…………exp−jM−1α1exp−jP−1β1exp−jM−1α2exp−jP−1β2…exp−jM−1αKexp−jP−1βK.According to the above estimationF^−1, we can define the following expression by reconstructing equation (13):(20)EQ=BBΤ…BTM−1Θ,where(21)B=11…1exp−jβ1exp−jβ2…exp−jβK…………exp−jP−1β1exp−jP−1β2…exp−jP−1βK,Τ=diagexp−jα1,exp−jα2,…,exp−jαK∈CK×K,αk=2πdfksinθkc,k=1,2,…,K.According to reconstructed equation (20), we can use the method described below to estimate the DOA.We define the following matrices:(22)EQ1=BBΤ…BΤM−2Θ,EQ2=BTBT2…BTM−1Θ.We can also define a matrixD since D=EQ1+EQ2. Then, according to the definitions of EQ1 and EQ2, we can express D in a no-noise situation:(23)D=Θ−1TΘ.Therefore, we can take the diagonal elements ofD, and then, we can obtain ϖkk=1,2,…,K, where αk=2πdfksinθk/c(k=1,2,…,K, to obtain the estimation of the DOA:(24)θ^k=arcsinc2πf^kdangleϖk,k=1,2,…,K.


### 2.2.4. The Steps of the Proposed Method

Thus far, we have given the complete process for automatically pairing DOA and frequency estimations in a linear array. The main steps required to implement the method proposed in this paper are as follows:(i)
Step 1: according to permutation matrixJ and equation (12), we reconstruct the covariance matrix R.(ii)
Step 2: we apply eigenvalue decomposition toR, and then we reconstruct the signal subspace Ess. According to equations (14) and (15), we construct matrices E1 and E2, respectively.(iii)
Step 3: we use equationΨ=E1+E2 for eigenvalue decomposition to obtain F−1 and Φ^. Finally, we estimate the frequency parameter f^k according to equation (18).(iv)
Step 4: we can obtain matrixEQ according to the reconstruction of Ess in equation (13). Then, we can also construct matrices EQ1 and EQ2.(v)
Step 5: we calculateD=EQ1+EQ2 to obtain matrix D. Finally, we estimate the DOA parameter θ^k according to equation (24).


### 2.2.5. The Condition of a Nonuniform Array

In this section, we first present the method proposed in this paper when the distances between the array elements are not equal. Then, we present the main steps for implementing the method in the case of a nonuniform array.We assume that the first elementd1 = 0 and that the distance between the mth element and the first element is dm. Then, we can transform equation (4) into the following form:(25)A1=11…1exp−jd2η1exp−jd2η2…exp−jd2ηK…………exp−jdMη1exp−jdMη2…exp−jdMηK,where ηk=2πfksinθk/ck=1,2,…,K. At the same time, equation (19) undergoes the following transformation:(26)A1ΦP−1=exp−jP−1β1exp−jP−1β2…exp−jP−1βKexp−jd2η1exp−jP−1β1exp−jd2η2exp−jP−1β2…exp−jd2ηKexp−jP−1βK…………exp−jdMη1exp−jP−1β1exp−jdMη2exp−jP−1β2…exp−jdMηKexp−jP−1βK.Similarly, we can reconstruct equation (13) and define the following expression:(27)EQQ=BJ1BJ2…BJMΘ,where Jm=diagexp−jdmη1,exp−jdmη2,…,exp−jdmηK∈CK×K, m=1,2,…,M.We define the following matrix:(28)EQQ1=BJ1BJ2…BJM−1Θ=EQQ11EQQ12…EQQ1M−1,EQQ2=BJ2BJ3…BJMΘ=EQQ22EQQ23…EQQ2M.We also define the matrixQm, where Qm=EQQ1m−1+EQQ2m, m=2,…,M. Therefore, we can take the diagonal elements of Qm and then obtain(29)νm=diagexp−jdm−dm−1η1,exp−jdm−dm−1η2,…,exp−jdm−dm−1ηK∈CK×K,m=2,…,M.We sort the diagonal elements and then define the following matrix:(30)V=11…1exp−jd2−d1η1exp−jd2−d1η2…exp−jd2−d1ηK…………exp−jdm−dm−1η1exp−jdm−dm−1η2…exp−jdm−dm−1ηK.According to equation (30), we can obtain the estimation of the DOA:(31)θ^k=1M−1∑m=2Marcsinc2πf^kdm−dm−1angleVmk,k=1,2,…,K.The main steps for implementing the method in this paper under the condition of a nonuniform array are as follows:(i)
Step 1: according to equations (10), (12), and (25), we reconstruct the new covariance matrix.(ii)
Step 2: we apply eigenvalue decomposition to the new covariance matrix, and then we reconstruct the signal subspaceEss. According to equations (14) and (15), we construct matrices E1 and E2, respectively.(iii)
Step 3: we use equationΨ=E1+E2 to perform eigenvalue decomposition and obtain F−1 and Φ^. Finally, we estimate the frequency parameter f^k according to equation (18).(iv)
Step 4: we can obtain the matrixEQQ according to equation (27). Then, we can also construct matrices EQQ1 and EQQ2.(v)
Step 5: we can obtain matrixV according to equation (30). Finally, we estimate the DOA parameter θ^k according to equation (31).


## 2.2.1. The Preprocessing Procedure

First, we obtain the covariance matrixRY=YYH of the received signal. To make full use of the conjugate information contained in the received signal, we define the permutation matrix J [27]:(10)J=0…01…0100………10…0PM×PM.Therefore, we can construct the following matrix:(11)RJ=JY∗Y∗HJH.We add the covariance matrixRY and RJ from equation (11), and then we average them. The form of the obtained covariance matrix is shown as follows:(12)R=RY+RJ2.Through analysis, we can obtain that the new total covariance matrixR is a Hermitian matrix (PM × PM) [28]. Therefore, we can apply eigenvalue decomposition, and then we can reconstruct the signal subspace Ess. In a no-noise situation, Ess can be approximately expressed as(13)Ess=AAΦ…AΦP−1F,where F is a full-rank matrix with K×K dimensions.Remark 1.
As mentioned earlier, the new total covariance matrixR is already a Hermitian matrix. According to Hermitian matrix characteristics, we assume that the diagonal matrix of the eigenvalues of R is G. Then, there exists a unitary matrix U, which assures RU = UG. Therefore, we can treat R as the unitary matrix U by using this correlation feature to further reduce the complexity of the proposed method and then propose a much lower complexity method.


## 2.2.2. Frequency Estimation

We define the following parameters:(14)E1=AAΦ…AΦP−2F,(15)E2=AΦAΦ2…AΦP−1F.Therefore, equations (14) and (15) have the following relationship:(16)E2=AΦAΦ2…AΦP−1F=AAΦ1…AΦP−2FF−1ΦF=E1F−1ΦF.LetΨ=F−1ΦF and Ψ=E1+E2. According to the LS-ESPRIT algorithm, we can estimate Φ by the eigenvalue decomposition of Ψ, and we can also estimate the matrix F−1 by the eigenvector of Φ. In a no-noise situation, we define(17)F^−1=F−1Θ,Φ^=Θ−1ΦΘ,where Θ is a fuzzy column matrix. Since Ψ and Φ have the same eigenvalues, we can obtain the eigenvalues λkk=1,2,…,K from matrix Ψ. As shown in equation (6) βk=2πfkτ,k=1,2,…,K, it is obvious that we can estimate the frequency parameter f^k,k=1,2,…,K:(18)f^k=12πτangleλk.


## 2.2.3. DOA Estimation

AΦP−1 has the following expression:(19)AΦP−1=exp−jP−1β1exp−jP−1β2…exp−jP−1βKexp−jα1exp−jP−1β1exp−jα2exp−jP−1β2…exp−jαKexp−jP−1βK…………exp−jM−1α1exp−jP−1β1exp−jM−1α2exp−jP−1β2…exp−jM−1αKexp−jP−1βK.According to the above estimationF^−1, we can define the following expression by reconstructing equation (13):(20)EQ=BBΤ…BTM−1Θ,where(21)B=11…1exp−jβ1exp−jβ2…exp−jβK…………exp−jP−1β1exp−jP−1β2…exp−jP−1βK,Τ=diagexp−jα1,exp−jα2,…,exp−jαK∈CK×K,αk=2πdfksinθkc,k=1,2,…,K.According to reconstructed equation (20), we can use the method described below to estimate the DOA.We define the following matrices:(22)EQ1=BBΤ…BΤM−2Θ,EQ2=BTBT2…BTM−1Θ.We can also define a matrixD since D=EQ1+EQ2. Then, according to the definitions of EQ1 and EQ2, we can express D in a no-noise situation:(23)D=Θ−1TΘ.Therefore, we can take the diagonal elements ofD, and then, we can obtain ϖkk=1,2,…,K, where αk=2πdfksinθk/c(k=1,2,…,K, to obtain the estimation of the DOA:(24)θ^k=arcsinc2πf^kdangleϖk,k=1,2,…,K.


## 2.2.4. The Steps of the Proposed Method

Thus far, we have given the complete process for automatically pairing DOA and frequency estimations in a linear array. The main steps required to implement the method proposed in this paper are as follows:(i)
Step 1: according to permutation matrixJ and equation (12), we reconstruct the covariance matrix R.(ii)
Step 2: we apply eigenvalue decomposition toR, and then we reconstruct the signal subspace Ess. According to equations (14) and (15), we construct matrices E1 and E2, respectively.(iii)
Step 3: we use equationΨ=E1+E2 for eigenvalue decomposition to obtain F−1 and Φ^. Finally, we estimate the frequency parameter f^k according to equation (18).(iv)
Step 4: we can obtain matrixEQ according to the reconstruction of Ess in equation (13). Then, we can also construct matrices EQ1 and EQ2.(v)
Step 5: we calculateD=EQ1+EQ2 to obtain matrix D. Finally, we estimate the DOA parameter θ^k according to equation (24).


## 2.2.5. The Condition of a Nonuniform Array

In this section, we first present the method proposed in this paper when the distances between the array elements are not equal. Then, we present the main steps for implementing the method in the case of a nonuniform array.We assume that the first elementd1 = 0 and that the distance between the mth element and the first element is dm. Then, we can transform equation (4) into the following form:(25)A1=11…1exp−jd2η1exp−jd2η2…exp−jd2ηK…………exp−jdMη1exp−jdMη2…exp−jdMηK,where ηk=2πfksinθk/ck=1,2,…,K. At the same time, equation (19) undergoes the following transformation:(26)A1ΦP−1=exp−jP−1β1exp−jP−1β2…exp−jP−1βKexp−jd2η1exp−jP−1β1exp−jd2η2exp−jP−1β2…exp−jd2ηKexp−jP−1βK…………exp−jdMη1exp−jP−1β1exp−jdMη2exp−jP−1β2…exp−jdMηKexp−jP−1βK.Similarly, we can reconstruct equation (13) and define the following expression:(27)EQQ=BJ1BJ2…BJMΘ,where Jm=diagexp−jdmη1,exp−jdmη2,…,exp−jdmηK∈CK×K, m=1,2,…,M.We define the following matrix:(28)EQQ1=BJ1BJ2…BJM−1Θ=EQQ11EQQ12…EQQ1M−1,EQQ2=BJ2BJ3…BJMΘ=EQQ22EQQ23…EQQ2M.We also define the matrixQm, where Qm=EQQ1m−1+EQQ2m, m=2,…,M. Therefore, we can take the diagonal elements of Qm and then obtain(29)νm=diagexp−jdm−dm−1η1,exp−jdm−dm−1η2,…,exp−jdm−dm−1ηK∈CK×K,m=2,…,M.We sort the diagonal elements and then define the following matrix:(30)V=11…1exp−jd2−d1η1exp−jd2−d1η2…exp−jd2−d1ηK…………exp−jdm−dm−1η1exp−jdm−dm−1η2…exp−jdm−dm−1ηK.According to equation (30), we can obtain the estimation of the DOA:(31)θ^k=1M−1∑m=2Marcsinc2πf^kdm−dm−1angleVmk,k=1,2,…,K.The main steps for implementing the method in this paper under the condition of a nonuniform array are as follows:(i)
Step 1: according to equations (10), (12), and (25), we reconstruct the new covariance matrix.(ii)
Step 2: we apply eigenvalue decomposition to the new covariance matrix, and then we reconstruct the signal subspaceEss. According to equations (14) and (15), we construct matrices E1 and E2, respectively.(iii)
Step 3: we use equationΨ=E1+E2 to perform eigenvalue decomposition and obtain F−1 and Φ^. Finally, we estimate the frequency parameter f^k according to equation (18).(iv)
Step 4: we can obtain the matrixEQQ according to equation (27). Then, we can also construct matrices EQQ1 and EQQ2.(v)
Step 5: we can obtain matrixV according to equation (30). Finally, we estimate the DOA parameter θ^k according to equation (31).


## 3. Results and Discussion

### 3.1. Performance Analysis of the Proposed Method

#### 3.1.1. Method Complexity

In this section, we focus on the performance analysis with respect to complexity. Complexity is mainly measured by the number of complex multiplications and the running time required by a given method. For the ESPRIT method in [16], the complexity is OM2P2N+M3P3+2K2MP−1+8K3+2K2M−1. For the improved ESPRIT method in [20], the complexity required to calculate the covariance matrix RY is OM2P2N. The complexity required for eigenvalue decomposition is OM3P3. The complexity of calculating Ψ=E1+E2 is O2K2MP−1+2K3. Then, the eigenvalue decomposition complexity of Ψ=E1+E2 is OK3. When estimating the DOA, the complexity is O2K3+2K2M−1P. Therefore, the complexity of the improved ESPRIT method is OM2P2N+M3P3+2K2MP−1+5K3+2K2M−1P.For the proposed method, the preprocessing complexity isOM2P2N+M3P3. The complexity of frequency estimation is O2K2MP−1+3K3. In addition, the complexity of DOA estimation is O2K3+2K2M−1P. Therefore, the complexity of the proposed method is OM2P2N+M3P3+2K2MP−1+5K3+2K2M−1P+1. It should be noted that in the preprocessing of this paper, we only need to calculate the covariance matrix RY and the eigenvalue decomposition, which means that we do not require additional calculations to construct the matrix RJ.The reason is that, according to equation (11), RJ has the following expression:(32)RJ=JY∗Y∗HJH=0…01…0100………10…0Y0Y1…YP−1∗Y0Y1…YP−1∗H0…01…0100………10…0=YP−1∗…Y1∗Y0∗Y0Y1…YP−10…01…0100………10…0,=YP−1∗Y0YP−1∗Y1…YP−1∗YP−1…………Y1∗Y0Y1∗Y1…Y1∗YP−1Y0∗Y0Y0∗Y1..….Y0∗YP−10…01…0100………10…0=YP−1∗YP−1…YP−1∗Y1YP−1∗Y0…………Y1∗YP−1…Y1∗Y1Y1∗Y0Y0∗YP−1…Y0∗Y1Y0∗Y0,while the covariance matrix RY(33)RY=YYH=Y0Y1…YP−1Y0Y1…YP−1H=Y0Y1…YP−1Y0∗Y1∗…YP−1∗=Y0Y0∗Y0Y1∗…Y0YP−1∗Y1Y0∗Y1Y1∗…Y1YP−1∗………YP−1Y0∗YP−1Y1∗…YP−1YP−1∗.By observing equations (32) and (33), we find that through simple moment transformation, we can transform the matrix RY into the matrix RJ. Therefore, when preprocessing, we do not need additional complex multiplications to reconstruct the matrix RJ.For the PM, the complexity of frequency estimation isOM2P2N+4K3+M2P2K+PMK2+2K2M−K. In addition, the complexity of DOA estimation is OK2M−K+2K3+2K2M−1. Therefore, the complexity of the PM is OM2P2N+6K3+M2P2K+PMK2+3K2M−K+2K2M−1. The CR method is also applicable to the PM. Therefore, the complexity of the CR-PM is OM2P2N+6K3+M2P2K+PMK2+3K2M−K+2K2M−1.Figures2 and 3 present the complexity comparison of these algorithms versus the number of signal sources K and the number of snapshots N with M = 12 and P = 3, respectively. Table 1 compares the running time of these algorithms under the condition of an i7-8550U CPU with K = 3, N = 200, and 2000 Monte Carlo simulations. Figures 2 and 3 show that the complexity of the method proposed in this paper is almost the same as that of the ESPRIT method in [16] and that of the improved ESPRIT method in [20] and is much higher than that of the PM and that of the CR-PM. In addition, the running time of the proposed method does not increase much. Moreover, through subsequent analysis, within the SNR range from −15 dB to −1 dB, the advantages of the proposed method are more obvious. In particular, when SNR = −15 dB, compared with the improved ESPRIT method, the frequency estimation accuracy of the method proposed in this paper is an approximately 25.50% improvement; the DOA estimation accuracy of the method proposed in this paper is an approximately 31.95% improvement. Therefore, we can confirm that by increasing the utilization of the originally received data, we can improve the parameter estimation accuracy and the noise robustness of the proposed method.Figure 2
Comparison of algorithm complexity under different snapshotsN.Figure 3
Comparison of algorithm complexity under different signal sourcesK.Table 1
Running time of these methods.
MethodsRunning time (s)PM2.3291CR-PM2.6634ESPRIT2.8223Improved ESPRIT2.8454CR-ESPRIT3.1177


#### 3.1.2. The Advantages of the Proposed Method

In this section, we summarize the advantages of the proposed method in this paper as follows:(1)
Under the condition of a uniform or a nonuniform array, the method can effectively estimate the DOAs and frequencies of source signals. It can also realize automatic pairing without an additional parameter pairing process since the method has the same fuzzy column matrix for both parameters.(2)
For incoherent signal sources whose angles are close together, this method can perform effective identification and parameter estimation.(3)
Compared with those of the PM, the CR-PM, the ESPRIT method [16], and the improved ESPRIT method [20], the frequency and DOA estimation accuracies of the proposed method are greatly improved, and this method has superior estimation performance. Moreover, the proposed method has better anti-noise performance and stronger robustness.


### 3.2. Numerical Simulation

In the simulation, we assume that the array receives signals emitted byK incoherent far-field sources. We also use the root mean square error (RMSE) metric to evaluate the DOA and frequency estimation performances of the proposed method; we define the RMSEs as(34)RMSEDOA=1K∑k=1K1L∑lLθ^k,l−θk2,RMSEfrequency=1K∑k=1K1L∑lLf^k,l−fk2,where θ^k,l and f^k,l are the estimated values of θk and fk, respectively, in the lth Monte Carlo simulation and L is the number of Monte Carlo simulations. In this paper, we set L = 2000.


#### 3.2.1. Performance Analysis of the Proposed Method in a Uniform Array

In this section, we assume that the array receives signals emitted by three incoherent far-field sources. The DOAs and operating frequencies of the signals are (θ1, f1) = (15°, 1 MHz), (θ2, f2) = (40°, 2.1 MHz), and (θ3, f3) = (50°, 3.1 MHz). SNR = 0 dB, M = 12 is the number of array elements, P = 3 is the number of delay values, d = 50 denotes the distances between the array elements, and N = 400 and K = 3 are the numbers of snapshots and signal sources, respectively. The scatter diagram of the joint frequency and DOA estimation of the proposed method in this paper is shown in Figure 4. Figure 4 shows that the proposed method is efficient in estimating the frequency and DOA results for a uniform array.Figure 4
Parameter estimation for a uniform array.


#### 3.2.2. Performance Analysis under Different Numbers of Array ElementsM

We setd = 50 m, K = 3, P = 3, and N = 400. We also set different numbers of array elements (M = 8, 12, and 16). The SNR range is from −15 dB to 15 dB (step: 2 dB), and the RMSEs of the frequency and DOA estimates of the method proposed in this paper are shown in Figures 5 and 6, respectively.Figure 5
RMSEs offk under different numbers of array elements.Figure 6
RMSEs of the DOA under different array elements.We can see from Figures5 and 6 that the method proposed in this paper can achieve high estimation performance within the SNR range of −15 dB to 15 dB (step: 2 dB) under different numbers of array elements. The estimation performance is stable under the condition of a low SNR. Moreover, we can see that the SNR has a great impact on the estimation accuracies of the frequency and DOA. The higher the SNR is, the higher the parameter estimation accuracies of the method for these two parameters. With the increase in the number of array elements, the DOA and frequency estimation accuracies of the method proposed in this paper improve. Furthermore, the RMSEs of the proposed method are greatly reduced. This is because as the number of array elements increases, the space diversity gain increases [29].


#### 3.2.3. Performance Analysis under Different Numbers of SnapshotsN

We setd = 50 m, K = 3, P = 3, and M = 12. We also set different numbers of snapshots (N = 100, 400, and 800). The SNR range is from −15 dB to 15 dB (step: 2 dB), and the RMSEs of the frequency and DOA estimations of the method proposed in this paper are shown in Figure 7 and 8, respectively. We can see from Figure 7 and 8 that when the SNR is within the range of −15 dB to 15 dB (with steps of 2 dB), the RMSEs of the proposed method demonstrate that with different snapshots, the algorithm can still maintain high estimation performance. Even in a situation with a low SNR, the estimation performance is still stable. As the number of snapshots increases, the estimation accuracy of the method proposed in this paper is enhanced, the performance is more precise, and the RMSEs of the frequency and DOA estimations of the proposed method decrease.Figure 7
RMSEs offk under different numbers of snapshots.Figure 8
RMSEs of the DOA under different numbers of snapshots.


#### 3.2.4. Performance Analysis under Different Delay ValuesP

We setd = 50 m, K = 3, N = 400, and M = 12. We also set different delay values (P = 2, 3, and 4). The range of the SNR is from −15 dB to 15 dB (step: 2 dB), and the RMSEs of the frequency and DOA estimations of the method proposed in this paper are shown in Figures 9 and 10, respectively. In Figures 9 and 10, under different delay values, the proposed method maintains high DOA and frequency estimation performance when the SNR ranges from −15 dB to 15 dB. The estimation performance is stable even in a situation with a low SNR. As the delay value increases, the estimation accuracy of the method proposed in this paper is enhanced, the performance is more precise, and the RMSEs of the DOA and frequency estimations of the method decrease.Figure 9
RMSEs offk under different delay values.Figure 10
RMSEs of the DOA under different delay values.


#### 3.2.5. Performance Analysis under Different Numbers of Signal SourcesK

We setd = 50 m, P = 3, N = 400, and M = 12. We also set different numbers of signal sources (K = 2, 3, and 4). The range of the SNR is from −15 dB to 15 dB (step: 2 dB), and the RMSEs of the frequency and DOA estimations of the method proposed in this paper are shown in Figures 11 and 12, respectively. In Figures 11 and 12, under different numbers of signal sources, the proposed method maintains high DOA and frequency estimation performance when the SNR ranges from −15 dB to 15 dB. The estimation performance is stable even in a situation with a low SNR. As the number of signal sources increases, the estimation accuracy of the method proposed in this paper makes the performance more imprecise, and the RMSEs of the DOA and frequency estimations of the proposed method increase. As the number of signal sources increases, the interference between sources increases, and the frequency and DOA estimation performances deteriorate [30].Figure 11
RMSEs offk under different numbers of signal sources.Figure 12
RMSEs of the DOA under different numbers of signal sources.


#### 3.2.6. Identification Performance Analysis for Signal Sources with Close Angles

We setd = 50 m, K = 2, P = 4, N = 200, and M = 12. In this section, we focus on exploring the recognition and identification abilities of the proposed method when the signal sources are at relatively close angles. The DOAs and operating frequencies of the signals are (θ1, f1) = (15°, 1 MHz) and (θ2, f2) = (17°, 2.1 MHz). We also set SNR = 5 dB. As shown in Figure 13, for signal sources with close angles, the proposed method can also perform effective identification and parameter estimation.Figure 13
The identification ability of the proposed method.


#### 3.2.7. Performance Analysis of the Proposed Method in a Nonuniform Array

In an actual field receiving system, the assumed reception model is different from the true model even after a calibration procedure [31]. Therefore, in this section, we mainly discuss the performance analysis under the condition of a nonuniform array, such as array element position deviation [32] and uneven distance between array elements.In this section, we assume that the array receives signals emitted by three incoherent far-field sources. The DOAs and operating frequencies of the signals are (θ1, f1) = (15°, 1 MHz), (θ2, f2) = (20°, 1.9 MHz), and (θ3, f3) = (30°, 2.8 MHz). SNR = 0 dB, M = 12, P = 5, N = 300, and K = 3. For analyzing array element position deviation, we set d = [0; 50; 101; 149; 200; 251; 300; 352; 401; 449; 501; 548] m. For uneven distance between array elements, we set d = [0; 40; 100; 150; 195; 235; 310; 365; 395; 450; 500; 540] m. The scatter diagram of the joint frequency and DOA estimation of the proposed method in this paper is shown in Figures 14 and 15. Figures 14 and 15 show that the proposed method is efficient in estimating the frequency and DOA results for both nonuniform array conditions.Figure 14
Parameter estimation for array element position deviation.Figure 15
Parameter estimation for uneven distances between array elements.


#### 3.2.8. Analysis of the Performances of Different Methods

In this section, we focus on analyzing the performances of different methods. We assume that the array receives signals emitted by two incoherent far-field sources. The DOAs and operating frequencies of the signals are (θ1, f1) = (15°, 1 MHz) and (θ2, f2) = (40°, 2.1 MHz). We set d = 50 m, K = 2, P = 2, N = 400, and M = 12. The range of the SNR is from −15 dB to 15 dB (step: 2 dB), and the RMSEs of the Cramer–Rao lower bound (CRLB), the PM, the CR-PM, the ESPRIT method [16], the improved ESPRIT method [20], and the method proposed in this paper with respect to the frequency and DOA estimations are shown in Figures 16 and 17, respectively.Figure 16
RMSEs offk under different methods.Figure 17
RMSEs of the DOA under different methods.To quantitatively illustrate, under the condition of a low SNR, compared with the improved ESPRIT method [20], the estimate improvement of the method proposed in this paper, we define the relative improvement ratio as(35)ratiofrequency=1−RMSEfrequencyCR−ESPRITRMSEfrequencyImproved ESPRIT×100%,(36)ratioDOA=1−RMSEDOACR−ESPRITRMSEDOAImproved ESPRIT×100%.According to the definitions of equations (35) and (36), we show the relative improvement ratio in Figures 18 and 19.Figure 18
The relative improvement ratios of the frequency estimation.Figure 19
The relative improvement ratios of the DOA estimation.As shown in Figures16 and 17, when the SNR is within the range of −15 dB to 15 dB (in steps of 2 dB), the estimation accuracy of the proposed method is better than that of the PM, the CR-PM, the ESPRIT method [16], and the improved ESPRIT method [20] in terms of both the DOA and frequency. Among them, the ESPRIT method has extremely poor angle estimation accuracy since it cannot automatically pair parameters.As shown in Figures18 and 19, when SNR = −15 dB to −1 dB, compared with the improved ESPRIT method, the estimation accuracy of the proposed method is greatly improved. In particular, when SNR = −15 dB, compared with the improved ESPRIT method, the frequency estimation accuracy of the method proposed in this paper is an approximately 25.50% improvement; the DOA estimation accuracy of the method proposed in this paper is an approximately 31.95% improvement. However, when SNR = −1 dB to 15 dB, the relative improvement ratio fluctuates around zero. This result illustrates that the estimation accuracy of the proposed method is almost the same as that of the improved ESPRIT method. For fluctuation, we surmise that the reason for this phenomenon may be the result of too few simulations in this paper.In summary, a comprehensive analysis of Figures16–19 shows that the estimation accuracy of the proposed method is improved over that of the PM, the CR-PM, the ESPRIT method, and the improved ESPRIT method. The results further verify that the method proposed in this paper has good anti-noise performance and stability under different SNRs. Therefore, compared to the PM, the CR-PM, the ESPRIT method, and the improved ESPRIT method, the method proposed in this paper is more suitable for use in a complex electromagnetic environment.


## 3.1. Performance Analysis of the Proposed Method

### 3.1.1. Method Complexity

In this section, we focus on the performance analysis with respect to complexity. Complexity is mainly measured by the number of complex multiplications and the running time required by a given method. For the ESPRIT method in [16], the complexity is OM2P2N+M3P3+2K2MP−1+8K3+2K2M−1. For the improved ESPRIT method in [20], the complexity required to calculate the covariance matrix RY is OM2P2N. The complexity required for eigenvalue decomposition is OM3P3. The complexity of calculating Ψ=E1+E2 is O2K2MP−1+2K3. Then, the eigenvalue decomposition complexity of Ψ=E1+E2 is OK3. When estimating the DOA, the complexity is O2K3+2K2M−1P. Therefore, the complexity of the improved ESPRIT method is OM2P2N+M3P3+2K2MP−1+5K3+2K2M−1P.For the proposed method, the preprocessing complexity isOM2P2N+M3P3. The complexity of frequency estimation is O2K2MP−1+3K3. In addition, the complexity of DOA estimation is O2K3+2K2M−1P. Therefore, the complexity of the proposed method is OM2P2N+M3P3+2K2MP−1+5K3+2K2M−1P+1. It should be noted that in the preprocessing of this paper, we only need to calculate the covariance matrix RY and the eigenvalue decomposition, which means that we do not require additional calculations to construct the matrix RJ.The reason is that, according to equation (11), RJ has the following expression:(32)RJ=JY∗Y∗HJH=0…01…0100………10…0Y0Y1…YP−1∗Y0Y1…YP−1∗H0…01…0100………10…0=YP−1∗…Y1∗Y0∗Y0Y1…YP−10…01…0100………10…0,=YP−1∗Y0YP−1∗Y1…YP−1∗YP−1…………Y1∗Y0Y1∗Y1…Y1∗YP−1Y0∗Y0Y0∗Y1..….Y0∗YP−10…01…0100………10…0=YP−1∗YP−1…YP−1∗Y1YP−1∗Y0…………Y1∗YP−1…Y1∗Y1Y1∗Y0Y0∗YP−1…Y0∗Y1Y0∗Y0,while the covariance matrix RY(33)RY=YYH=Y0Y1…YP−1Y0Y1…YP−1H=Y0Y1…YP−1Y0∗Y1∗…YP−1∗=Y0Y0∗Y0Y1∗…Y0YP−1∗Y1Y0∗Y1Y1∗…Y1YP−1∗………YP−1Y0∗YP−1Y1∗…YP−1YP−1∗.By observing equations (32) and (33), we find that through simple moment transformation, we can transform the matrix RY into the matrix RJ. Therefore, when preprocessing, we do not need additional complex multiplications to reconstruct the matrix RJ.For the PM, the complexity of frequency estimation isOM2P2N+4K3+M2P2K+PMK2+2K2M−K. In addition, the complexity of DOA estimation is OK2M−K+2K3+2K2M−1. Therefore, the complexity of the PM is OM2P2N+6K3+M2P2K+PMK2+3K2M−K+2K2M−1. The CR method is also applicable to the PM. Therefore, the complexity of the CR-PM is OM2P2N+6K3+M2P2K+PMK2+3K2M−K+2K2M−1.Figures2 and 3 present the complexity comparison of these algorithms versus the number of signal sources K and the number of snapshots N with M = 12 and P = 3, respectively. Table 1 compares the running time of these algorithms under the condition of an i7-8550U CPU with K = 3, N = 200, and 2000 Monte Carlo simulations. Figures 2 and 3 show that the complexity of the method proposed in this paper is almost the same as that of the ESPRIT method in [16] and that of the improved ESPRIT method in [20] and is much higher than that of the PM and that of the CR-PM. In addition, the running time of the proposed method does not increase much. Moreover, through subsequent analysis, within the SNR range from −15 dB to −1 dB, the advantages of the proposed method are more obvious. In particular, when SNR = −15 dB, compared with the improved ESPRIT method, the frequency estimation accuracy of the method proposed in this paper is an approximately 25.50% improvement; the DOA estimation accuracy of the method proposed in this paper is an approximately 31.95% improvement. Therefore, we can confirm that by increasing the utilization of the originally received data, we can improve the parameter estimation accuracy and the noise robustness of the proposed method.Figure 2
Comparison of algorithm complexity under different snapshotsN.Figure 3
Comparison of algorithm complexity under different signal sourcesK.Table 1
Running time of these methods.
MethodsRunning time (s)PM2.3291CR-PM2.6634ESPRIT2.8223Improved ESPRIT2.8454CR-ESPRIT3.1177


### 3.1.2. The Advantages of the Proposed Method

In this section, we summarize the advantages of the proposed method in this paper as follows:(1)
Under the condition of a uniform or a nonuniform array, the method can effectively estimate the DOAs and frequencies of source signals. It can also realize automatic pairing without an additional parameter pairing process since the method has the same fuzzy column matrix for both parameters.(2)
For incoherent signal sources whose angles are close together, this method can perform effective identification and parameter estimation.(3)
Compared with those of the PM, the CR-PM, the ESPRIT method [16], and the improved ESPRIT method [20], the frequency and DOA estimation accuracies of the proposed method are greatly improved, and this method has superior estimation performance. Moreover, the proposed method has better anti-noise performance and stronger robustness.


## 3.1.1. Method Complexity

In this section, we focus on the performance analysis with respect to complexity. Complexity is mainly measured by the number of complex multiplications and the running time required by a given method. For the ESPRIT method in [16], the complexity is OM2P2N+M3P3+2K2MP−1+8K3+2K2M−1. For the improved ESPRIT method in [20], the complexity required to calculate the covariance matrix RY is OM2P2N. The complexity required for eigenvalue decomposition is OM3P3. The complexity of calculating Ψ=E1+E2 is O2K2MP−1+2K3. Then, the eigenvalue decomposition complexity of Ψ=E1+E2 is OK3. When estimating the DOA, the complexity is O2K3+2K2M−1P. Therefore, the complexity of the improved ESPRIT method is OM2P2N+M3P3+2K2MP−1+5K3+2K2M−1P.For the proposed method, the preprocessing complexity isOM2P2N+M3P3. The complexity of frequency estimation is O2K2MP−1+3K3. In addition, the complexity of DOA estimation is O2K3+2K2M−1P. Therefore, the complexity of the proposed method is OM2P2N+M3P3+2K2MP−1+5K3+2K2M−1P+1. It should be noted that in the preprocessing of this paper, we only need to calculate the covariance matrix RY and the eigenvalue decomposition, which means that we do not require additional calculations to construct the matrix RJ.The reason is that, according to equation (11), RJ has the following expression:(32)RJ=JY∗Y∗HJH=0…01…0100………10…0Y0Y1…YP−1∗Y0Y1…YP−1∗H0…01…0100………10…0=YP−1∗…Y1∗Y0∗Y0Y1…YP−10…01…0100………10…0,=YP−1∗Y0YP−1∗Y1…YP−1∗YP−1…………Y1∗Y0Y1∗Y1…Y1∗YP−1Y0∗Y0Y0∗Y1..….Y0∗YP−10…01…0100………10…0=YP−1∗YP−1…YP−1∗Y1YP−1∗Y0…………Y1∗YP−1…Y1∗Y1Y1∗Y0Y0∗YP−1…Y0∗Y1Y0∗Y0,while the covariance matrix RY(33)RY=YYH=Y0Y1…YP−1Y0Y1…YP−1H=Y0Y1…YP−1Y0∗Y1∗…YP−1∗=Y0Y0∗Y0Y1∗…Y0YP−1∗Y1Y0∗Y1Y1∗…Y1YP−1∗………YP−1Y0∗YP−1Y1∗…YP−1YP−1∗.By observing equations (32) and (33), we find that through simple moment transformation, we can transform the matrix RY into the matrix RJ. Therefore, when preprocessing, we do not need additional complex multiplications to reconstruct the matrix RJ.For the PM, the complexity of frequency estimation isOM2P2N+4K3+M2P2K+PMK2+2K2M−K. In addition, the complexity of DOA estimation is OK2M−K+2K3+2K2M−1. Therefore, the complexity of the PM is OM2P2N+6K3+M2P2K+PMK2+3K2M−K+2K2M−1. The CR method is also applicable to the PM. Therefore, the complexity of the CR-PM is OM2P2N+6K3+M2P2K+PMK2+3K2M−K+2K2M−1.Figures2 and 3 present the complexity comparison of these algorithms versus the number of signal sources K and the number of snapshots N with M = 12 and P = 3, respectively. Table 1 compares the running time of these algorithms under the condition of an i7-8550U CPU with K = 3, N = 200, and 2000 Monte Carlo simulations. Figures 2 and 3 show that the complexity of the method proposed in this paper is almost the same as that of the ESPRIT method in [16] and that of the improved ESPRIT method in [20] and is much higher than that of the PM and that of the CR-PM. In addition, the running time of the proposed method does not increase much. Moreover, through subsequent analysis, within the SNR range from −15 dB to −1 dB, the advantages of the proposed method are more obvious. In particular, when SNR = −15 dB, compared with the improved ESPRIT method, the frequency estimation accuracy of the method proposed in this paper is an approximately 25.50% improvement; the DOA estimation accuracy of the method proposed in this paper is an approximately 31.95% improvement. Therefore, we can confirm that by increasing the utilization of the originally received data, we can improve the parameter estimation accuracy and the noise robustness of the proposed method.Figure 2
Comparison of algorithm complexity under different snapshotsN.Figure 3
Comparison of algorithm complexity under different signal sourcesK.Table 1
Running time of these methods.
MethodsRunning time (s)PM2.3291CR-PM2.6634ESPRIT2.8223Improved ESPRIT2.8454CR-ESPRIT3.1177


## 3.1.2. The Advantages of the Proposed Method

In this section, we summarize the advantages of the proposed method in this paper as follows:(1)
Under the condition of a uniform or a nonuniform array, the method can effectively estimate the DOAs and frequencies of source signals. It can also realize automatic pairing without an additional parameter pairing process since the method has the same fuzzy column matrix for both parameters.(2)
For incoherent signal sources whose angles are close together, this method can perform effective identification and parameter estimation.(3)
Compared with those of the PM, the CR-PM, the ESPRIT method [16], and the improved ESPRIT method [20], the frequency and DOA estimation accuracies of the proposed method are greatly improved, and this method has superior estimation performance. Moreover, the proposed method has better anti-noise performance and stronger robustness.


## 3.2. Numerical Simulation

In the simulation, we assume that the array receives signals emitted byK incoherent far-field sources. We also use the root mean square error (RMSE) metric to evaluate the DOA and frequency estimation performances of the proposed method; we define the RMSEs as(34)RMSEDOA=1K∑k=1K1L∑lLθ^k,l−θk2,RMSEfrequency=1K∑k=1K1L∑lLf^k,l−fk2,where θ^k,l and f^k,l are the estimated values of θk and fk, respectively, in the lth Monte Carlo simulation and L is the number of Monte Carlo simulations. In this paper, we set L = 2000.


### 3.2.1. Performance Analysis of the Proposed Method in a Uniform Array

In this section, we assume that the array receives signals emitted by three incoherent far-field sources. The DOAs and operating frequencies of the signals are (θ1, f1) = (15°, 1 MHz), (θ2, f2) = (40°, 2.1 MHz), and (θ3, f3) = (50°, 3.1 MHz). SNR = 0 dB, M = 12 is the number of array elements, P = 3 is the number of delay values, d = 50 denotes the distances between the array elements, and N = 400 and K = 3 are the numbers of snapshots and signal sources, respectively. The scatter diagram of the joint frequency and DOA estimation of the proposed method in this paper is shown in Figure 4. Figure 4 shows that the proposed method is efficient in estimating the frequency and DOA results for a uniform array.Figure 4
Parameter estimation for a uniform array.


### 3.2.2. Performance Analysis under Different Numbers of Array ElementsM

We setd = 50 m, K = 3, P = 3, and N = 400. We also set different numbers of array elements (M = 8, 12, and 16). The SNR range is from −15 dB to 15 dB (step: 2 dB), and the RMSEs of the frequency and DOA estimates of the method proposed in this paper are shown in Figures 5 and 6, respectively.Figure 5
RMSEs offk under different numbers of array elements.Figure 6
RMSEs of the DOA under different array elements.We can see from Figures5 and 6 that the method proposed in this paper can achieve high estimation performance within the SNR range of −15 dB to 15 dB (step: 2 dB) under different numbers of array elements. The estimation performance is stable under the condition of a low SNR. Moreover, we can see that the SNR has a great impact on the estimation accuracies of the frequency and DOA. The higher the SNR is, the higher the parameter estimation accuracies of the method for these two parameters. With the increase in the number of array elements, the DOA and frequency estimation accuracies of the method proposed in this paper improve. Furthermore, the RMSEs of the proposed method are greatly reduced. This is because as the number of array elements increases, the space diversity gain increases [29].


### 3.2.3. Performance Analysis under Different Numbers of SnapshotsN

We setd = 50 m, K = 3, P = 3, and M = 12. We also set different numbers of snapshots (N = 100, 400, and 800). The SNR range is from −15 dB to 15 dB (step: 2 dB), and the RMSEs of the frequency and DOA estimations of the method proposed in this paper are shown in Figure 7 and 8, respectively. We can see from Figure 7 and 8 that when the SNR is within the range of −15 dB to 15 dB (with steps of 2 dB), the RMSEs of the proposed method demonstrate that with different snapshots, the algorithm can still maintain high estimation performance. Even in a situation with a low SNR, the estimation performance is still stable. As the number of snapshots increases, the estimation accuracy of the method proposed in this paper is enhanced, the performance is more precise, and the RMSEs of the frequency and DOA estimations of the proposed method decrease.Figure 7
RMSEs offk under different numbers of snapshots.Figure 8
RMSEs of the DOA under different numbers of snapshots.


### 3.2.4. Performance Analysis under Different Delay ValuesP

We setd = 50 m, K = 3, N = 400, and M = 12. We also set different delay values (P = 2, 3, and 4). The range of the SNR is from −15 dB to 15 dB (step: 2 dB), and the RMSEs of the frequency and DOA estimations of the method proposed in this paper are shown in Figures 9 and 10, respectively. In Figures 9 and 10, under different delay values, the proposed method maintains high DOA and frequency estimation performance when the SNR ranges from −15 dB to 15 dB. The estimation performance is stable even in a situation with a low SNR. As the delay value increases, the estimation accuracy of the method proposed in this paper is enhanced, the performance is more precise, and the RMSEs of the DOA and frequency estimations of the method decrease.Figure 9
RMSEs offk under different delay values.Figure 10
RMSEs of the DOA under different delay values.


### 3.2.5. Performance Analysis under Different Numbers of Signal SourcesK

We setd = 50 m, P = 3, N = 400, and M = 12. We also set different numbers of signal sources (K = 2, 3, and 4). The range of the SNR is from −15 dB to 15 dB (step: 2 dB), and the RMSEs of the frequency and DOA estimations of the method proposed in this paper are shown in Figures 11 and 12, respectively. In Figures 11 and 12, under different numbers of signal sources, the proposed method maintains high DOA and frequency estimation performance when the SNR ranges from −15 dB to 15 dB. The estimation performance is stable even in a situation with a low SNR. As the number of signal sources increases, the estimation accuracy of the method proposed in this paper makes the performance more imprecise, and the RMSEs of the DOA and frequency estimations of the proposed method increase. As the number of signal sources increases, the interference between sources increases, and the frequency and DOA estimation performances deteriorate [30].Figure 11
RMSEs offk under different numbers of signal sources.Figure 12
RMSEs of the DOA under different numbers of signal sources.


### 3.2.6. Identification Performance Analysis for Signal Sources with Close Angles

We setd = 50 m, K = 2, P = 4, N = 200, and M = 12. In this section, we focus on exploring the recognition and identification abilities of the proposed method when the signal sources are at relatively close angles. The DOAs and operating frequencies of the signals are (θ1, f1) = (15°, 1 MHz) and (θ2, f2) = (17°, 2.1 MHz). We also set SNR = 5 dB. As shown in Figure 13, for signal sources with close angles, the proposed method can also perform effective identification and parameter estimation.Figure 13
The identification ability of the proposed method.


### 3.2.7. Performance Analysis of the Proposed Method in a Nonuniform Array

In an actual field receiving system, the assumed reception model is different from the true model even after a calibration procedure [31]. Therefore, in this section, we mainly discuss the performance analysis under the condition of a nonuniform array, such as array element position deviation [32] and uneven distance between array elements.In this section, we assume that the array receives signals emitted by three incoherent far-field sources. The DOAs and operating frequencies of the signals are (θ1, f1) = (15°, 1 MHz), (θ2, f2) = (20°, 1.9 MHz), and (θ3, f3) = (30°, 2.8 MHz). SNR = 0 dB, M = 12, P = 5, N = 300, and K = 3. For analyzing array element position deviation, we set d = [0; 50; 101; 149; 200; 251; 300; 352; 401; 449; 501; 548] m. For uneven distance between array elements, we set d = [0; 40; 100; 150; 195; 235; 310; 365; 395; 450; 500; 540] m. The scatter diagram of the joint frequency and DOA estimation of the proposed method in this paper is shown in Figures 14 and 15. Figures 14 and 15 show that the proposed method is efficient in estimating the frequency and DOA results for both nonuniform array conditions.Figure 14
Parameter estimation for array element position deviation.Figure 15
Parameter estimation for uneven distances between array elements.


### 3.2.8. Analysis of the Performances of Different Methods

In this section, we focus on analyzing the performances of different methods. We assume that the array receives signals emitted by two incoherent far-field sources. The DOAs and operating frequencies of the signals are (θ1, f1) = (15°, 1 MHz) and (θ2, f2) = (40°, 2.1 MHz). We set d = 50 m, K = 2, P = 2, N = 400, and M = 12. The range of the SNR is from −15 dB to 15 dB (step: 2 dB), and the RMSEs of the Cramer–Rao lower bound (CRLB), the PM, the CR-PM, the ESPRIT method [16], the improved ESPRIT method [20], and the method proposed in this paper with respect to the frequency and DOA estimations are shown in Figures 16 and 17, respectively.Figure 16
RMSEs offk under different methods.Figure 17
RMSEs of the DOA under different methods.To quantitatively illustrate, under the condition of a low SNR, compared with the improved ESPRIT method [20], the estimate improvement of the method proposed in this paper, we define the relative improvement ratio as(35)ratiofrequency=1−RMSEfrequencyCR−ESPRITRMSEfrequencyImproved ESPRIT×100%,(36)ratioDOA=1−RMSEDOACR−ESPRITRMSEDOAImproved ESPRIT×100%.According to the definitions of equations (35) and (36), we show the relative improvement ratio in Figures 18 and 19.Figure 18
The relative improvement ratios of the frequency estimation.Figure 19
The relative improvement ratios of the DOA estimation.As shown in Figures16 and 17, when the SNR is within the range of −15 dB to 15 dB (in steps of 2 dB), the estimation accuracy of the proposed method is better than that of the PM, the CR-PM, the ESPRIT method [16], and the improved ESPRIT method [20] in terms of both the DOA and frequency. Among them, the ESPRIT method has extremely poor angle estimation accuracy since it cannot automatically pair parameters.As shown in Figures18 and 19, when SNR = −15 dB to −1 dB, compared with the improved ESPRIT method, the estimation accuracy of the proposed method is greatly improved. In particular, when SNR = −15 dB, compared with the improved ESPRIT method, the frequency estimation accuracy of the method proposed in this paper is an approximately 25.50% improvement; the DOA estimation accuracy of the method proposed in this paper is an approximately 31.95% improvement. However, when SNR = −1 dB to 15 dB, the relative improvement ratio fluctuates around zero. This result illustrates that the estimation accuracy of the proposed method is almost the same as that of the improved ESPRIT method. For fluctuation, we surmise that the reason for this phenomenon may be the result of too few simulations in this paper.In summary, a comprehensive analysis of Figures16–19 shows that the estimation accuracy of the proposed method is improved over that of the PM, the CR-PM, the ESPRIT method, and the improved ESPRIT method. The results further verify that the method proposed in this paper has good anti-noise performance and stability under different SNRs. Therefore, compared to the PM, the CR-PM, the ESPRIT method, and the improved ESPRIT method, the method proposed in this paper is more suitable for use in a complex electromagnetic environment.


## 3.2.1. Performance Analysis of the Proposed Method in a Uniform Array

In this section, we assume that the array receives signals emitted by three incoherent far-field sources. The DOAs and operating frequencies of the signals are (θ1, f1) = (15°, 1 MHz), (θ2, f2) = (40°, 2.1 MHz), and (θ3, f3) = (50°, 3.1 MHz). SNR = 0 dB, M = 12 is the number of array elements, P = 3 is the number of delay values, d = 50 denotes the distances between the array elements, and N = 400 and K = 3 are the numbers of snapshots and signal sources, respectively. The scatter diagram of the joint frequency and DOA estimation of the proposed method in this paper is shown in Figure 4. Figure 4 shows that the proposed method is efficient in estimating the frequency and DOA results for a uniform array.Figure 4
Parameter estimation for a uniform array.


## 3.2.2. Performance Analysis under Different Numbers of Array ElementsM

We setd = 50 m, K = 3, P = 3, and N = 400. We also set different numbers of array elements (M = 8, 12, and 16). The SNR range is from −15 dB to 15 dB (step: 2 dB), and the RMSEs of the frequency and DOA estimates of the method proposed in this paper are shown in Figures 5 and 6, respectively.Figure 5
RMSEs offk under different numbers of array elements.Figure 6
RMSEs of the DOA under different array elements.We can see from Figures5 and 6 that the method proposed in this paper can achieve high estimation performance within the SNR range of −15 dB to 15 dB (step: 2 dB) under different numbers of array elements. The estimation performance is stable under the condition of a low SNR. Moreover, we can see that the SNR has a great impact on the estimation accuracies of the frequency and DOA. The higher the SNR is, the higher the parameter estimation accuracies of the method for these two parameters. With the increase in the number of array elements, the DOA and frequency estimation accuracies of the method proposed in this paper improve. Furthermore, the RMSEs of the proposed method are greatly reduced. This is because as the number of array elements increases, the space diversity gain increases [29].


## 3.2.3. Performance Analysis under Different Numbers of SnapshotsN

We setd = 50 m, K = 3, P = 3, and M = 12. We also set different numbers of snapshots (N = 100, 400, and 800). The SNR range is from −15 dB to 15 dB (step: 2 dB), and the RMSEs of the frequency and DOA estimations of the method proposed in this paper are shown in Figure 7 and 8, respectively. We can see from Figure 7 and 8 that when the SNR is within the range of −15 dB to 15 dB (with steps of 2 dB), the RMSEs of the proposed method demonstrate that with different snapshots, the algorithm can still maintain high estimation performance. Even in a situation with a low SNR, the estimation performance is still stable. As the number of snapshots increases, the estimation accuracy of the method proposed in this paper is enhanced, the performance is more precise, and the RMSEs of the frequency and DOA estimations of the proposed method decrease.Figure 7
RMSEs offk under different numbers of snapshots.Figure 8
RMSEs of the DOA under different numbers of snapshots.


## 3.2.4. Performance Analysis under Different Delay ValuesP

We setd = 50 m, K = 3, N = 400, and M = 12. We also set different delay values (P = 2, 3, and 4). The range of the SNR is from −15 dB to 15 dB (step: 2 dB), and the RMSEs of the frequency and DOA estimations of the method proposed in this paper are shown in Figures 9 and 10, respectively. In Figures 9 and 10, under different delay values, the proposed method maintains high DOA and frequency estimation performance when the SNR ranges from −15 dB to 15 dB. The estimation performance is stable even in a situation with a low SNR. As the delay value increases, the estimation accuracy of the method proposed in this paper is enhanced, the performance is more precise, and the RMSEs of the DOA and frequency estimations of the method decrease.Figure 9
RMSEs offk under different delay values.Figure 10
RMSEs of the DOA under different delay values.


## 3.2.5. Performance Analysis under Different Numbers of Signal SourcesK

We setd = 50 m, P = 3, N = 400, and M = 12. We also set different numbers of signal sources (K = 2, 3, and 4). The range of the SNR is from −15 dB to 15 dB (step: 2 dB), and the RMSEs of the frequency and DOA estimations of the method proposed in this paper are shown in Figures 11 and 12, respectively. In Figures 11 and 12, under different numbers of signal sources, the proposed method maintains high DOA and frequency estimation performance when the SNR ranges from −15 dB to 15 dB. The estimation performance is stable even in a situation with a low SNR. As the number of signal sources increases, the estimation accuracy of the method proposed in this paper makes the performance more imprecise, and the RMSEs of the DOA and frequency estimations of the proposed method increase. As the number of signal sources increases, the interference between sources increases, and the frequency and DOA estimation performances deteriorate [30].Figure 11
RMSEs offk under different numbers of signal sources.Figure 12
RMSEs of the DOA under different numbers of signal sources.


## 3.2.6. Identification Performance Analysis for Signal Sources with Close Angles

We setd = 50 m, K = 2, P = 4, N = 200, and M = 12. In this section, we focus on exploring the recognition and identification abilities of the proposed method when the signal sources are at relatively close angles. The DOAs and operating frequencies of the signals are (θ1, f1) = (15°, 1 MHz) and (θ2, f2) = (17°, 2.1 MHz). We also set SNR = 5 dB. As shown in Figure 13, for signal sources with close angles, the proposed method can also perform effective identification and parameter estimation.Figure 13
The identification ability of the proposed method.


## 3.2.7. Performance Analysis of the Proposed Method in a Nonuniform Array

In an actual field receiving system, the assumed reception model is different from the true model even after a calibration procedure [31]. Therefore, in this section, we mainly discuss the performance analysis under the condition of a nonuniform array, such as array element position deviation [32] and uneven distance between array elements.In this section, we assume that the array receives signals emitted by three incoherent far-field sources. The DOAs and operating frequencies of the signals are (θ1, f1) = (15°, 1 MHz), (θ2, f2) = (20°, 1.9 MHz), and (θ3, f3) = (30°, 2.8 MHz). SNR = 0 dB, M = 12, P = 5, N = 300, and K = 3. For analyzing array element position deviation, we set d = [0; 50; 101; 149; 200; 251; 300; 352; 401; 449; 501; 548] m. For uneven distance between array elements, we set d = [0; 40; 100; 150; 195; 235; 310; 365; 395; 450; 500; 540] m. The scatter diagram of the joint frequency and DOA estimation of the proposed method in this paper is shown in Figures 14 and 15. Figures 14 and 15 show that the proposed method is efficient in estimating the frequency and DOA results for both nonuniform array conditions.Figure 14
Parameter estimation for array element position deviation.Figure 15
Parameter estimation for uneven distances between array elements.


## 3.2.8. Analysis of the Performances of Different Methods

In this section, we focus on analyzing the performances of different methods. We assume that the array receives signals emitted by two incoherent far-field sources. The DOAs and operating frequencies of the signals are (θ1, f1) = (15°, 1 MHz) and (θ2, f2) = (40°, 2.1 MHz). We set d = 50 m, K = 2, P = 2, N = 400, and M = 12. The range of the SNR is from −15 dB to 15 dB (step: 2 dB), and the RMSEs of the Cramer–Rao lower bound (CRLB), the PM, the CR-PM, the ESPRIT method [16], the improved ESPRIT method [20], and the method proposed in this paper with respect to the frequency and DOA estimations are shown in Figures 16 and 17, respectively.Figure 16
RMSEs offk under different methods.Figure 17
RMSEs of the DOA under different methods.To quantitatively illustrate, under the condition of a low SNR, compared with the improved ESPRIT method [20], the estimate improvement of the method proposed in this paper, we define the relative improvement ratio as(35)ratiofrequency=1−RMSEfrequencyCR−ESPRITRMSEfrequencyImproved ESPRIT×100%,(36)ratioDOA=1−RMSEDOACR−ESPRITRMSEDOAImproved ESPRIT×100%.According to the definitions of equations (35) and (36), we show the relative improvement ratio in Figures 18 and 19.Figure 18
The relative improvement ratios of the frequency estimation.Figure 19
The relative improvement ratios of the DOA estimation.As shown in Figures16 and 17, when the SNR is within the range of −15 dB to 15 dB (in steps of 2 dB), the estimation accuracy of the proposed method is better than that of the PM, the CR-PM, the ESPRIT method [16], and the improved ESPRIT method [20] in terms of both the DOA and frequency. Among them, the ESPRIT method has extremely poor angle estimation accuracy since it cannot automatically pair parameters.As shown in Figures18 and 19, when SNR = −15 dB to −1 dB, compared with the improved ESPRIT method, the estimation accuracy of the proposed method is greatly improved. In particular, when SNR = −15 dB, compared with the improved ESPRIT method, the frequency estimation accuracy of the method proposed in this paper is an approximately 25.50% improvement; the DOA estimation accuracy of the method proposed in this paper is an approximately 31.95% improvement. However, when SNR = −1 dB to 15 dB, the relative improvement ratio fluctuates around zero. This result illustrates that the estimation accuracy of the proposed method is almost the same as that of the improved ESPRIT method. For fluctuation, we surmise that the reason for this phenomenon may be the result of too few simulations in this paper.In summary, a comprehensive analysis of Figures16–19 shows that the estimation accuracy of the proposed method is improved over that of the PM, the CR-PM, the ESPRIT method, and the improved ESPRIT method. The results further verify that the method proposed in this paper has good anti-noise performance and stability under different SNRs. Therefore, compared to the PM, the CR-PM, the ESPRIT method, and the improved ESPRIT method, the method proposed in this paper is more suitable for use in a complex electromagnetic environment.


## 4. Conclusions

For linear arrays, this paper proposes a joint angle and frequency estimation method based on CR-ESPRIT. We first preprocess the received signal by taking full advantage of the conjugate information contained in the originally received data, and we reconstruct a new total covariance matrix. Then, we use the LS-ESPRIT algorithm to estimate the frequency parameter. According to the unique relationship between angles and frequencies, we estimate the DOAs based on the reconstructed received signal. The complexity of the method proposed in this paper is almost the same as that of the ESPRIT and the improved ESPRIT. Numerical simulations and comparisons with the PM, the CR-PM, the ESPRIT method, and the improved ESPRIT method prove the superiority of the proposed method. In a real space environment, under the condition of a uniform or a nonuniform array, this method can realize the automatic pairing of the estimated DOAs and frequencies of radiation source signals without an additional parameter pairing process. Moreover, this method has high accuracy and strong anti-noise performance when conducting parameter estimation.



---

*Source: 5477848-2021-07-09.xml*