Direct Position Determination of Non-Circular Sources for Multiple Arrays via Weighted Euler ESPRIT Data Fusion Method

: In recent years, direct position determination (DPD) with multiple arrays for non-circular (NC) signals is a hot topic to research. Conventional DPD techniques with spectral peak search methods have high computational complexity and are sensitive to the locations of the observation stations. Besides, there will be loss when the signal propagates in the air, which leads to different received signal-to-noise ratios (SNRs) for each observation station. To attack the problems mentioned above, this paper derives direct position determination of non-circular sources for multiple arrays via weighted Euler estimating signal parameters viarotational invariance techniques (ESPRIT) data fusion (NC-Euler-WESPRIT) method. Firstly, elliptic covariance information of NC signals and Euler transformation are used to extend the received signal. Secondly, ESPRIT is applied to avoid the high-dimensional spectral function search problem of each observation station. Then, we combine the information of all observation stations to construct a spectral function without complex multiplication to reduce the computational complexity. Finally, the data of each observation station is weighted to compensate for the projection error. The consequence of simulation indicates that the proposed NC-Euler-WESPRIT algorithm not only improves the estimation performance, but also greatly reduces the computational complexity compared with subspace data fusion (SDF) technology and NC-ESPRIT algorithm.


Introduction
As an important research area about signal processing, passive location has attracted extensive attention recently. It has applications in many fields including sonar, communication and astronomy [1][2][3][4]. With the development of radio technology, the electromagnetic environment becomes more and more complex which makes it hard for the conventional single antenna to extract accurate and useful information [5]. Besides, single antenna can only get the direction information of the targets without distance and other information. Instead, multiple arrays have the advantages of flexible beam control [6], higher spatial resolution [7] and higher signal gain [8]. Original two-step location technology needs to estimate the intermediate parameters from the original data, such as time of arrival (TOA) [9], time difference of arrival (TDOA) [10], angle of arrival (AOA) [11]. Then, the locations of targets can be estimated through the spatial geometric relationship [12]. Instead of estimating the intermediate parameters, DPD can directly obtain the estimation of target locations from the original data [13]. Therefore, DPD avoids the error caused by modeling the intermediate parameters and can acquire better positioning performance than original two-step location [14].
According to the elliptic covariance characteristics of the source, sources could be divided into two categories: circular sources and non-circular sources [15]. At present, a (1) The proposed NC-Euler-WESPRIT DPD algorithm takes full advantage of elliptic covariance information of NC signals to expand the virtual array aperture. Compared with original two-step localization technique and SDF algorithm, NC-Euler-WESPRIT DPD increases the available degrees of freedom (DOF). (2) Euler transformation is applied to decrease the complexity by converting complex number calculation into real number operation and ESPRIT is applied to avoid the high-dimensional spectral function search problem of each observation station. Then we combine the information of all observation stations to construct a spectral function without complex multiplication to further reduce the computational complexity. (3) In practice, there will be loss when the signal propagates in the air and SNRs of received signals at each observation station are often different. Therefore, a specific weight is set to compensate for the projection error and get better positioning performance. (4) Complexity analysis, Cramer Rao lower bound (CRLB) and simulation consequence are given to check the effectiveness and superiority of the NC-Euler-WESPRIT DPD algorithm.
The following structure of the article is: A data model is introduced in Section 2. A NC-Euler-WESPRIT DPD algorithm is derived in Section 3. Performance analysis, simulation experiments and conclusions are put in Sections 4 and 5, respectively. It should be noted that vectors and matrices are lower-case bold and upper-case bold, respectively; diag{·} represents diagonal matrix; || · || represents the l 2 norm; | · | indicates absolute value; (·) * , (·) T , (·) H denote conjugate, transposition and conjugate transpose, respectively; Re(·) and Im(·) denote real part and imaginary part, respectively; I N denotes N dimensional unit array; ⊗ denotes Kronecker product; ∂(a)/∂(b) stands for the derivative of a with respect to b; (·) + and (·) −1 denote general Moore-Penrose inverse and the inverse of the matrix.

Model Formulation
For ease of illustration, Figure 1 shows the direct positioning scene with multiple arrays. Consider that there are Q unknown far-field incoherent narrowband sources impacting G known observation stations. Each observation station is outfitted with a uniform linear array (ULA) of M elements horizontally placed along x-axis and the interval between arrays is d. The number of snapshots is K. In this two-dimensional plane, we assume that the unknown target locations are p q = x q , y q T (q = 1, 2, · · · , Q) and known positions of observation stations are u g = x g , y g T (g = 1, 2, · · · , G). In practice, there will be loss when the signal propagates in the air [31,32]. In other words, when a target impinges on different observation stations, SNR of received signals at each observation station often varies greatly, so we define the path propagation loss coefficient as [33] where P q and P g,q stand for the power of the qth emitter signal and the power of the signal from the qth target received at the gth observation location, respectively. Supposing the emitter signal sent to the gth observation station at the kth sampling time from the qth target is s g,q (k). For simplicity, we only consider the signal whose non-circular rate is 1, i.e., Equation (2) is satisfied [15] where ϕ denotes non-circular phase and (·) T denotes transposition. Then, the signal received by the gth observation station at the kth sampling time can be stated as [33] r g (k) = Q ∑ q=1 α g,q a θ g,q s g,q (k) + n g (k), where θ g,q represents DOA of the qth target received by the gth observation station, a θ g,q = 1, e −j 2πd λ cos θ g,q , · · · , e −j 2πd λ denotes the wavelength and n g (k) denotes noise vector, which obeys Gaussian distribution. According to the previous far-field assumption, the observation station could be approximated to a point. Hence, we rewrite Equation (4).
Then, Equation (3) can be rewritten as where A g (p), Λ g , s g (k) represent the direction matrix, propagation loss matrix and source vector of the gth observation station, respectively.

NC-ESPRIT-DPD
From [25], the non-circular signal can be written in the following form where Φ represents the NC phase matrix where ϕ q represents NC phase of the qth source and s (R) g,q (k) stands for the amplitude of the signal. Thus, Equation (6) can be rewritten according to Equation (11): The received signal can be extended by making full use of the elliptic covariance information of NC sources [34]: where B g can be considered as an extended direction matrix and where row exchange matrix J is Define γ g,q = 2πd ||u g −p q || , then direction matrix B g can be expressed as where Next, the covariance matrix of the gth observation station is calculated by Then, we perform eigenvalue decomposition (EVD) of R g [27] Here, we use λ g,m (m = 1, · · · , Q, · · · , 2M) to denote eigenvalues of R g , which is listed in descending order. Their eigenvectors are denoted with e g,m . Hence, signal subspace and noise subspace in Equation (24) can be stated as U s g = e g,1 , · · · , e g,Q and U n g = e g,Q+1 , · · · , e g,2M , respectively. Σ g = diag λ g,1 , λ g,2 , · · · , λ g,2M .
On the basis of the orthogonality of signal subspace U s g and noise subspace U n g , we can get span B g = span U (s) g , i.e., there is a reversible matrix T that makes Define row exchange matrix where O ∈ R (M−1)×M is a zero matrix. Then, we can obtain where (·) + denotes general Moore-Penrose inverse and Then, γ g,q (q = 1, 2, · · · Q) could be easily acquired according to Equation (32). In view of the above definition, we can know that the phase of γ g,q contains the azimuth information of targets. Then, by combining the azimuth information of G observation stations, we can obtain the position information of targets. The joint spectral function of NC-ESPRIT-DPD algorithm can be constructed as where the qth diagonal element of diagonal matrix Γ g is µ g,q and | · | indicates absolute value.

NC-Euler-ESPRIT-DPD
On the basis of the method of converting complex number calculation into real number operation [29], Equation (15) can be rewritten: where H g can be considered as a direction matrix (·) * , Re(·) and Im(·) indicate conjugate, real part and imaginary part, respectively. As ||u g −p q || is defined above, direction matrix H g can be expressed as where Then, we can obtain the signal subspace U (s) g according to Equations (23), (24) and (34). As signal subspace and noise subspace are orthogonal, we can get span Define reversible matrix where O ∈ R (M−1)×M is a zero matrix. Then, we can obtain where the qth column element of Similarly, where the qth column element of According to Equations (45) and (47), and the sum difference product formula, we can obtain where The following relationship can be obtained from Equations (40) and (49): where (·) + denotes general Moore-Penrose inverse. Then, γ g,q (q = 1, 2, · · · Q) can be acquired by performing EVD on Equation (51). According to the above definition, the phase of γ g,q contains the azimuth information of targets. Then, by combining the azimuth information of G observation stations, we can obtain the position information of targets. The joint spectral function of NC-Euler-ESPRIT-DPD algorithm can be constructed as where ν g,q represents the qth diagonal element of diagonal matrix Ξ g .

NC-Euler-WESPRIT-DPD
However, the joint spectral function in Equation (52) ignores the influence of propagation loss matrix Λ g and treats the data of all observation stations equally. Therefore, the joint spectral function is easily disturbed by the observation station with poor performance. For the sake of balancing the projection error and obtain better performance, we should assign a weight to compensate for the error. As we all know, higher SNR will lead to smaller error. Then, according to the principle of power allocation [30], we can reasonably allocate a weight w g according to the SNR of each observation station.
Here, it is necessary to assume that noise power σ 2 n,g is invariant in the process of observation [28]. According to Equation (54), SNR of different observation stations is proportional to α 2 g,q P q = P g,q , which is unknown. Nevertheless, we note that Equation (54) can be rewritten as The eigenvalues of R g can be stated as where σ 2 s,m are eigenvalues of R s and can denote the estimated value of signal power. Then, we estimate noise power at the gth observation location: The corresponding estimated value of the received signal power iŝ Finally, the joint spectral function of NC-Euler-WESPRIT is constructed as Acquire Q maximum points through spectral peak search, which corresponds to the target positions. Moreover, it can be seen that the spectral function in Equation (59) does not include complex number multiplication, so compared with the original subspace data fusion (SDF) algorithm [27], the computational complexity of NC-Euler-WESPRIT-DPD algorithm is greatly reduced. The main steps of NC-Euler-WESPRIT-DPD algorithm are summarized as Algorithm 1.

Algorithm 1 NC-Euler-WESPRIT-DPD
Input: all the data collected from G observation stations Output: target positions 1: Obtain the extended received signal z g (k) according to Equations (14) and (34). 2: Calculate the covariance matrix R g according to Equation (23). 3: According to Equation (24), perform EVD on R g in step 2. 4: Estimate the noise powerσ 2 n,g and signal powerP g according to Equations (57) and (58). 5: Perform EVD on Equation (51) to obtain diagonal matrix Ξ g . 6: Construct the weighted joint spectral function f NC−Euler−WESPRIT (p) according to Equation (59). 7: The target positions can be acquired by spectral peak search.

Derivation of the CRLB
According to [27,35], we derive the CRLB for the above positioning scene.
The corresponding signal vector s 0 (k) and noise vector n 1 (k) are The array manifold H is: Thus, we can represent the receiving signal vector z(k) = Hs 0 (k) + n 1 (k) Assuming that noise vector n 1 (k) follows a complex Gaussian distribution, then P(z(1), · · · , z(k)) = 1 The log-likelihood function of Equation (65) is Defines 0k = Re[s 0 (k)] ands 0k = Im[s 0 (k)]. The parameters to be estimated are given by Fisher information matrix Ω = E χχ T −1 , where Eventually, the CRLB can be obtained by where

Complexity Analysis
Assuming that Q far-field targets impact G observation stations and each observation location is outfitted with a ULA of M elements. The number of snapshots is K. G x , G y and G θ stand for search grids of x, y and angles. The proposed NC-Euler-WESPRIT-DPD algorithm mainly includes four steps: estimate R g , EVD on R g , EVD on Equation (51) and the calculation of the spectral function. The corresponding computational complexity Since the spectral function of the proposed algorithm does not need complex multiplication and its computational complexity is difficult to measure, the average running time of the algorithm is used to measure the computational complexity. It is noted that the acquisition of weightP ĝ σ 2 n,g only needs simple real number operation and its computational complexity is basically negligible [28]. Table 1 shows the comparison of complexity.

Algorithms Complexity
Two-step [12] O GM 3 + (GK Considering Q = 4, G = 6, K = 300, M = 5, SNR = 15 dB. G x , G y and G θ take the same value and change from 200 to 900. Figure 2 simulates the running time of different algorithms varying with the number of search grids. From Figure 2, NC-Euler-WESPRIT has lower complexity than NC-ESPRIT because of the conversion of complex number calculation into real number operation. Although the complexity of NC-Euler-WESPRIT is higher than two-step localization technology [12], its localization performance is better. Specially, the spectral function of the proposed NC-Euler-WESPRIT algorithm does not need complex multiplication, its computational complexity is significantly reduced compared with SDF algorithm [27].

Simulation Environment
The proposed NC-Euler-WESPRIT DPD algorithm is experimented on a system having 16 processors with 64 GB RAM and 2.5 TB hard Disk. The code is running with Windows 10 and MATLAB R2018b version.

Simulation Results
To verify the effectiveness of NC-Euler-WESPRIT DPD algorithm, we applied Monte Carlo experiments to evaluate the performance. Define root mean square error (RMSE) as where Mon represents the number of Monte Carlo experiments and Q stands for the number of targets. x q,mn ,ŷ q,mn denotes the estimated position of the qth target in the mnth experiment, x q , y q represents the real position of the qth target. See Table 2 for other simulation parameters. Moreover, to show the heteroscedasticity of the received SNR at each observation station [36], we define α g = ᾱ g ,ᾱ g , · · · ,ᾱ g T ∈ R Q×1 (77)   Figure 3 and M = 8 in Figure 4. The spectral peak diagrams show that NC-Euler-WESPRIT algorithm can successfully distinguish the targets even when the number of array elements at each observation station is less than that of targets. Besides, the estimation performance of NC-Euler-WESPRIT algorithm is getting better and better with the increase of the number of array elements.  Simulation 2: Discuss the performance of NC-Euler-WESPRIT algorithm with different snapshots. Considering that five far-field targets impact six observation stations. SNR = 20 dB, M = 5. K = 50 in Figure 5 and K = 300 in Figure 6. The other simulation parameters are the same as simulation 1. The contour figures indicate that NC-Euler-WESPRIT algorithm can successfully distinguish the targets with small snapshots. In addition, the sidelobe of the spectral peak becomes smaller with the increase of snapshots.   Table 2. SNRs in Figures 7 and 8 are 15 dB and 30 dB, respectively. It can be seen from scatter diagrams that NC-Euler-WESPRIT method can distinguish targets at low SNR well. Moreover, the estimation performance of NC-Euler-WESPRIT method is getting better and better with the increase of SNR.    Table 2. From the figures, the performance of NC-Euler-WESPRIT method is superior to original two-step positioning method [12], SDF technology [27], NC-ESPRIT and NC-Euler-ESPRIT. Compared with the method before Euler transformation, the estimation performance of NC-Euler-ESPRIT algorithm is closed to it. In addition, the superior performance of NC-Euler-WESPRIT algorithm is more obvious when the heteroscedasticity of the received SNR at each observation station increases.   Table 2. The figures indicate that the performance of NC-Euler-WESPRIT method is superior to original two-step location method [12], SDF technology [27], NC-ESPRIT and NC-Euler-ESPRIT. Compared with the method before Euler transformation, the estimation performance of NC-Euler-ESPRIT algorithm is closed to it. In addition, the superior performance of the proposed NC-Euler-WESPRIT algorithm is more obvious when the heteroscedasticity of the received SNR at each observation station increases.

Advantages of the Proposed Algorithm
The advantages of NC-Euler-WESPRIT algorithm are summarized as follows: (1) The proposed NC-Euler-WESPRIT algorithm has more available degrees of freedom than original two-step positioning method, SDF method and can distinguish more targets. (2) The proposed NC-Euler-WESPRIT algorithm significantly reduces the computational complexity compared with SDF method and NC-ESPRIT algorithm. (3) By weighting the received data of each observation station, the proposed NC-Euler-WESPRIT algorithm has higher positioning precision than original two-step positioning method and SDF method.

Conclusions
This paper derives direct position determination of non-circular sources for multiple arrays via weighted Euler ESPRIT data fusion. The proposed NC-Euler-WESPRIT algorithm takes full advantage of elliptic covariance information of NC signals to expand the virtual array aperture. By comparison with original two-step positioning technology and SDF technology, NC-Euler-WESPRIT algorithm has more available DOF. Secondly, this paper introduces the idea of Euler-ESPRIT and constructs a spectral function without complex multiplication. According to the running time of the algorithms, the proposed algorithm greatly reduces the complexity compared with the original SDF algorithm and NC-ESPRIT algorithm. In addition, the proposed algorithm reduces the performance loss caused by signal propagation loss in space as much as possible by weighting the received data of each observation station. Simulation results verify that the proposed NC-Euler-WESPRIT algorithm can effectively estimate the targets even when the number of array elements of each observation station is less than that of targets and its positioning performance is superior to original two-step positioning method and SDF algorithm. While our proposed method only considers the situation of uniform linear array, sparse arrays are now the focus of attention, which is the future work due to the fact that sparse arrays will bring great challenge on the non-circular signals.

Conflicts of Interest:
The authors declare no conflict of interest.

Abbreviations
The following abbreviations are used in this manuscript: