A fourth-order cumulant orthonormal propagator rooting method based on Toeplitz approximation

A novel Toeplitz fourth-order cumulant (FOC) orthonormal propagator rooting method (TFOC ‐ OPRM) of direction-of-arrival (DOA) estimation for uniform linear array (ULA) is proposed in this paper. Specifically, the modified (i.e., reduced-dimension) FOC  (MFOC) matrix is achieved at first via removing the redundant information encompassed in the primary FOC matrix; then, the TFOC matrix which possesses Toeplitz structure can be recovered by utilizing the Toeplitz approximation method. To reduce the computational complexity, an effective method based on the polynomial rooting technology is adopted. Finally, the DOAs of incident signals can be estimated by exploiting orthonormal propagator rooting method. The theoretical analysis coupled with simulation results show that the proposed resultant algorithm can reduce the computational complexity significantly, as well as improve the estimation performance in both spatially white noise environment and spatially color noise environment.


Introduction
Direction-of-arrival estimation based on antenna array is one of the important directions of research hotspot in array signal processing, which has wide prospects of application in military and civil fields such as wireless communications, radar, passive sonar, biomedicine, and seismic exploration [1][2][3][4]. Various high-resolution algorithms, such as multiple signal classification (MUSIC) algorithm [5] and estimating signal parameter via rotational invariance technique (ESPRIT) [6] approaches, have been proposed to achieve direction-of-arrival (DOA) estimation of narrowband far-field signal sources. However, these subspace-based DOA estimation algorithms described above are not only very sensitive to the noise, but also require the noise's characteristics of the sensors in advance. Furthermore, it is restricted that the total number of sources acting on the array must be less than or equal with sensors [7]. When the constrained condition cannot be met in practical environments, the estimation performance of those aforementioned algorithms may run into a stone wall. Fortunately, much more attentions have been paid to this issue, and much more efforts have been made to overcome the above drawbacks. Motivated by the truth that the high-order cumulant-based (HOC) is asymptotically insensitive to Gaussian noise, which can be recognized as a promising technique for direction finding by adopting sensor array [8][9][10]. Besides, another key motivation of using HOC is the ability to resolve more number of sources than or equal to that of the array elements [11]. However, the process of eigenvalue decomposition (EVD) or singular value decomposition (SVD) requires large amount of calculation and time taken, which greatly affects the development of rapid source location. Marcos and co-workers [12,13] firstly proposed so-named propagator method (PM) to obtain the signal and noise subspaces by executing a linear-partition operation, which can decrease the computational complexity effectively. Specifically, the algorithm under the conditions of medium and high signal-to-noise ratio can achieve the same performance as that of the traditional high-resolution algorithms but with higher calculation efficiency. Base on [12,13], numerous modifications of PM have been proposed, such as [14,15], to achieve the low-complexity DOA estimation. In [16], an efficient HOPM algorithm was proposed by making full use of intrinsic multi-dimensional characteristics and affordable computability. AFOCbased and OPMlike (FOC -OPM) algorithm [17] was proposed to gain good location performance. However, the computational complexity of this method is high due to a great number of redundant information still existing in the FOC matrix. To mitigate this shortcoming, the improved FOC algorithm [18] was proposed to reduce the computational complexity. However, the performance of the algorithm cannot be asymptotically optimal due to the estimation error of the FOC matrix. Zhang et al. [19] derived a root -MUSIC method using a co-prime linear array to improve the estimation accuracy with low complexity. In [20][21][22], a similar polynomial root-based method was chosen to realize low-complexity for DOA estimation. In [23][24][25][26][27][28], the fractal structure modeling based was utilized to improve the performance of the fractal antennas. In [29,30], the wavelet analysis-based method was proposed to deal with application in signal processing.
In this paper, a novel TFOC -OPRM algorithm is introduced. The contributions of this paper are twofold: Firstly, the reduced dimension matrix is obtained to reduce the computational complexity by removing a large number of redundant elements from the original FOC matrix while maintaining the effective aperture of the virtual array in unchanged state. Secondly, the Toeplitz structure is recovered by the Toeplitz operation of the reduced dimension FOC matrix, and the DOA estimation of the recovered Toeplitz structure matrix is performed based on the polynomial root method.

Data model
Consider M narrowband far-field sources s i (t), (i = 1, ⋯, M) impinging on a uniform linear array(ULA) with N equispaced omnidirectional sensors, where the distance between adjacent sensors is equal to half the wavelength. Assume that the incoming sources are stationary and mutually independent. The noise is the additive white/color Gaussian one and statistically independent of the sources. Let the first sensor be the reference, and then, the observed data received in time t at the kth sensor can be expressed as where s i (t) is the ith source, n k (t) is the Gaussian noise at the kth sensor, and a k (θ i ) is the response of the kth sensor corresponding to the ith source and can be expressed as where λ is the central wavelength and d is the spacing between two adjacent sensors. Therefore, the matrix form of (1) can be expressed as where Assuming that the source signals are zero-mean stationary random process, the FOC can be defined as where x k m ðm ¼ 1; 2; 3; 4Þ is the stochastic process. Apparently, cumðk 1 ; k 2 ; k Ã 3 ; k Ã 4 Þ has N 4 values with the change of k 1 , k 2 , k 3 , k 4 . For simplicity, Eq. (4) can be written in matrix form, which is denoted by cumulant matrix C 4 , and cumðk 1 ; k 2 ; k Ã 3 ; k Ã 4 Þ appears as the [(k 1 − 1)N + k 2 ]th row and [(k 3 − 1)N + k 4 ]th column of C 4 .
where B and C S represent the extended array manifold and the FOC matrix of incident source signals, respectively. B = A ⊗ A and each column of B is bðθÞ ¼ aðθÞ aðθÞ. It is obvious that bðθÞ is a N 2 × 1 vector, which means that the array aperture of ULA is extended. That is, the number of resolved source signals is no less than that of sensors.

The effective array aperture extended
As proven in [31], an array of N arbitrary identical omnidirectional sensors can be extended to at most of N 2 − N + 1. Especially, the number of virtual elements is 2N − 1 for ULA according to [31]. In order to discuss the effective aperture of ULA, four real elements (N = 4) are considered, and bðθÞ can be expressed in detail as follows: where z = exp(j2π(d/λ) sin θ). Equation (6) shows that there is a lot of redundancy in expanded steering vector bðθÞ. That is, only from 1st to Nth and all kNth (k = 2, ⋯, N) items of the bðθÞ are valid, while others are redundant ones. To eliminate these repetitive elements, a (2N − 1) × (2N − 1) matrix R 4 is defined firstly. Next, the 1st to Nth and all kNth (k = 2, ⋯, N) rows of C 4 are taken out in sequence, and then these rows are stored in the 1st to (2N − 1)th row of the new matrix R 4 . The same operation is performed on the 1st to Nth and all kNth (k = 2, ⋯, N) columns of C 4 to obtain the 1st to (2N − 1)th columns of R 4 . Similar to Eq. (5), R 4 can be expressed as where D denotes the extended array manifold without redundancy, and each column of D has the form of d(θ) = [1, ⋯, z 2N − 2 ] T . Therefore, the reduced-dimension R 4 not only contains all of the information about original matrix C 4 , but also keeps the extended array aperture unchanged.

The TFOC-OPRM method
When the incident targets are considered as statistically independent signal sources, the ideal R 4 has a Toeplitz structure. However, in practical applications, for example, due to finite sampling snapshots and the low SNR, the matrix R 4 obtained at this time does not meet the Toeplitz structure anymore; instead, it becomes a diagonally dominant matrix. The happening of such condition will have a negative impact on the performance of the final DOA estimation. In order to improve the DOA estimation accuracy of the antenna array, the first task is to recover the Toeplitz structure of matrixR 4 , that is, to get Toeplitz matrixR 4T . Then, a R 4T of Toeplitz matrix can be approached to the real reduced dimension by solving the following optimization problem: where S T represents Toeplitz matrices, and the entries of the Toeplitz matrix R 4T can be written as where the element r p(p + h − 1) denotes the pth row and (p + h − 1)th column of R 4 , h ∈ [1, ⋯, 2N − 1]. And then R 4T can be obtained by the following Toeplitization operator: where Toep stands for the Toeplitization operator.
Although conventional algorithms, such as MUSIC and ESPRIT, can be applied to estimateDOAs based on the R 4T , the computational burden is much heavier due to the EVD and SVDinvolved. Therefore, we apply OPM for estimating the DOAs to reduce the complex computations effectively. The presented propagator method is based on the following partition where the dimensions of propagator matrix P is defined as a unique linear operator which satisfies the following condition Define , and combine with equation (11) Equation (13) shows that the R 4T is orthogonal to the columns of Q H , and the propagator matrix P can be obtained by minimizing the cost function ξ(P) where ‖•‖ F indicates the Frobenius norm, and the optimal solution P is given by In order to introduce the orthonormalization, the orthonormalized matrix Q 0 is obtained as follows Therefore, the following spectral function p(θ) can be formed to estimate the DOAs of source signals It can be seen from function (17) that the MDOAs of the incoming signals can be obtained by means of one-dimensional (1 − D) spectrum-peak search over θ. However, to further reduce the computational burden, we can improve function (17) to derive a more efficient search-free modification estimator in computation based on polynomial rooting [32]. In order to further reduce the computational complexity of the algorithm, the method based on polynomial roots is used to improve the spatial spectrum estimation function, so as to obtain more efficient estimators in the calculation, with the specific description of the algorithm given as follows.
Set z = exp(j2π(d/λ) sin θ), Then, the denominator of the estimator (17) can be re-expressed with the following polynomial format In an ideal condition, there should be exactly M numbers of roots, that is, z 1 , z 2 , ⋯, z M distributing over the unit circle, and these M numbers of roots are exactly the roots of the polynomial f(z). However, in practical application, due to the influence of various complex factors in the environment, the M roots of the equation f(z) cannot be strictly distributed on the unit circle. In this case, only M roots close to the unit circle need to be selected, similarly to the Root -MUSIC approach [32][33][34]. After M roots {z 1 , ⋯, z i , ⋯z M } are obtained, the DOA estimation of the incident target signal source can then be completed by the following formula So far as it is concerned, the specific operational steps of the proposed Toeplitz fourth-order cumulant orthogonal propagation method based on polynomial roots under limited sampling snapshots can be summarized as follows: Step 1 Estimate C 4 from the received data by (5).
Step 2 Obtain the dimension reduction matrix R 4 by removing the redundant items from the expanded matrix C 4 .
Step 3 Reconstruct the Toeplitz matrix R 4T by performing Toeplitz approximation on R 4 as formulas (9) and (10).
Step 4 Estimate the linear operator P according to Eqs. (14) and (15), then calculate the standard orthonormalized matrix Q 0 based on Eq. (16).
Step 5 Obtain the polynomial function (19), and further to it, obtain the M roots closest to the unit circle, i.e., the roots of the f(z).
Step 6 Obtain the direction estimation of the incoming wave of the incident target signal source from (20).

Complexity analysis
As for the analysis of computational complexity, the main parts of computation are considered, that is, the construction of the cumulant matrix, the linear operation, the spectral peak search operation, the Toeplitz operation, and the polynomial rooting operation. To further prove the superiority of the TFOC -OPRM algorithm in terms of computational complexity, FOC -OPM and MFOC -OPM are used as the comparative algorithms.
For the FOC -OPM technique, the main operation amount comes from three major parts, that is, to calculate the N 2 × N 2 cumulant matrix, to perform the linear operator of cumulant matrix, and to execute one spectral search. Therefore  For the proposed TFOC -OPRM algorithm, the major computational complexity comes from forming one (2N − 1) × (2N − 1) cumulant matrix, to perform Toeplitz operation, to perform the linear operator of cumulant matrix, and to execute once polynomial rooting operation. Therefore, the computational complexity isO(9(2N − 1) 2 From the above analysis, it can be obviously seen that the computational complexity of TFOC -OPRM algorithm proposed in this paper is significantly lower than that of both FOC -OPM algorithm and MFOC -OPM algorithm. The main reason is that the polynomial root method has been involved to reduce the computational complexity further.

Results and discussion
In this section, the proposed TFOC -OPRM algorithm, as well as FOC -OPM [17] andMFOC -OPM [18] algorithms that are used for the purpose of comparison are simulated in the environment of spatial white noise and spatial color noise respectively to verify the superiority of the proposed algorithm. In the simulation experiment, the ULA composed of three antenna elements(N = 3) is used, in which the interval between adjacent antenna elements is d = λ/2. It is assumed that there are three far-field narrow-band statistically independent target signal sources (M = 3), whose incident angles are {− 45°, 15°, 40°} respectively, with the Gaussian white/color noise being considered. Both the proposed TFOC -OPRM algorithm and these two comparative algorithms take 500 Monte-Carlo simulations each time as their estimated performance value. Two performance indexes, namely, normalized probability of success (NPS) and estimated root-mean-square errors (RMSEs), are defined to evaluate the performance of these three algorithms. And the RMSEs and NPS are respectively expressed as whereθ n ðiÞ refers to the estimated value of the reference value θ n in the ith time Monte Carlo trial. The ϒ suc and Τ total denote the times of successes and Monte Carlo trial, respectively. Furthermore, it should be noted that the defined success of a simulation experiment satisfies maxðjθ n − θ n jÞ < ε, and ε in the formula equals to 0.8 and 1.5 for experiments 2 and 3, respectively.

Experiment 1: the spatial spectrum estimation
In the first experiment, the input SNR and the number of snapshots are set to be 10 dB and 500, respectively. Figure 1 shows the spatial spectrum of the proposed TFOC -OPRM, FOC -OPM, and MFOC -OPM algorithms in both spatially white noise and spatially color noise environments. It can be observed from the curves in Fig. 1 that all of the three algorithms have successfully located the sharp peak corresponding to the incident angle. Note that the angular resolution performance of the three algorithms in spatially white noise situation provides better than that of spatially color noise situation.
Further analysis indicates that no matter whether in spatially white noise situation or in spatially color noise situation, the angular resolution of the proposed TFOC -OPRM algorithm is much higher than that of both MFOC -OPM and FOC -OPM algorithms. The reason is that the proposed TFOC -OPRM algorithm can recover the Toeplitz structure of R 4 , making the Toeplitz matrix R 4T closer to the real situation.

Experiment 2: RMSEs and NPS versus SNR
The main objective of this experiment is to evaluate the performance of the TFOC -OPRM algorithm, FOC -OPM algorithm and MFOC -OPM algorithm in terms of RMSEs and NPS with the change of input SNR. The number of sampling snapshots is L = 2000, the input SNRchanges from 8 to 24 dB, with the step being 2 dB. Figures 2  and 3 plot the RMSEs and NPS of DOA estimation with the proposed TFOC -OPRM algorithm and the comparison algorithms as the input SNR changes, respectively. As illustrated in Fig. 2, the RMSEs of the three algorithms decrease monotonically with the increase of the input SNR. Further analysis shows that in the environment of spatial white noise, with the increase of input SNR, the RMSEs performance curve of TFOC -OPRM algorithm is better than that of FOC -OPM algorithm and that of MFOC -OPM algorithm; in the environment of spatial color noise, the RMSEs performance curve of TFOC -OPRM algorithm is better than that of MFOC -OPM algorithm. That is, no matter whether in spatially white noise situation or in spatially-color noise situation, the RMSEs performance of the proposed TFOC -OPRM algorithm achieves better than that of both MFOC -OPM and FOC -OPM algorithms. In addition, when the input SNR changes between 8 and 14 dB, the proposed TFOC -OPRM algorithm manages to achieve almost the same RMSE performance as the FOC -OPM algorithm. But when the input SNR is higher than 14 dB, the performance of TFOC -OPRM becomes better than that of FOC -OPM. From Fig. 3, it can be concluded that the NPS performance of the proposed TFOC -OPRM algorithm is better than that of the FOC -OPM algorithm and MFOC -OPM algorithm in the case of low input SNR (between 8 and 10 dB) in spatially white noise situation. With the increase of input SNR, theNPS of all of the three algorithms ultimately is 1. It is also worth noting that when the SNR changes between 8 and 16 dB in spatially color noise situation, the performance of the proposed algorithm is better than that of the compared algorithms. In addition to that, the proposed algorithm not only removes a lot of redundant data in the original FOC, but also restores the Toeplitz structure of the reduced dimensional FOC. Moreover, it adopts the method of finding roots of polynomials. Therefore, the proposed TFOC - OPRM algorithm not only reduces the computational complexity, but also improves the accuracy of DOA estimation.

Experiment 3: RMSEs and NPS versus snapshots
The In other words, the estimated matrixR 4 deviates greatly from the ideal matrix R 4 due to the number of sampling snapshots that is relatively small. With the increasing number of sampling snapshots, we can see that the performance curves tend to be stable gradually. At the same time, it can be observed that the TFOC -OPRM algorithm proposed in this paper achieves more satisfactory estimation performance thanMFOC -OPM and FOC -OPM algorithms, either in the condition of spatial-white noise or in the condition of spatial-color noise. And the estimation performance of the three algorithms in spatially-white noise situation provide better than that of spatially color noise situation. Note that the computational complexity of proposed algorithm is significantly lower than that of the FOC -OPM algorithm due to the fact that the redundant information of the original cumulant matrix is removed. Moreover, the Toeplitz approxi- mate method is performed on the reduced-rank R 4 to improve estimation performance.
Meanwhile, compared to MFOC -OPM method, the TFOC -OPRM algorithm has lower computational burden, which exploits polynomial rooting instead of spectral search.

Experiment 4: the calculation complexity versus snapshots
In this simulation experiment, we further verify the advantages of TFOC -OPRM algorithm in terms of computational complexity, also by comparing the algorithms with FOC -OPM and MFOC -OPM. The number of incident target signal sources and the number of array elements of ULA are set as M = 3 and N = 3 respectively, with the interval of angular scanning being defined as Δθ = 0.01. Figure 6 shows the calculation complexity of the proposed TFOC -OPRMalgorithm and the comparison algorithms as the number of sampling snapshots changes (the number of sampling snapshots changes from L = 400 to L = 2000). Viewing from the simulation results in Fig. 6, with the increasing number of sampling snapshots, the computational complexity of the proposed TFOC -OPRM algorithm is significantly far lower than that of the FOC -OPM algorithm and the MFOC -OPM algorithm, and this advantage will be more obvious with the further increase of the number of sampling snapshots. The reason is that the proposed TFOC -OPRM algorithm not only eliminates a large number of redundant data in the original FOC, but also adopts the polynomial root method. This is consistent with the theoretical analysis given in Section 3.3 and testify the high-efficiency of the proposed TFOC -OPRM algorithm.

Conclusions
In this paper, a novel low computational complexity TFOC -OPRM localization algorithm for DOA estimation has been proposed in the presence of spatially white noise and spatially color noise environments. Specifically, we reconstruct a new Toeplitz matrix, which is close to the Toeplitz structure information in ideal condition via the Toeplitz approximate method. By exploiting the polynomial root method, the proposed TFOC -OPRM localization algorithm does not include large amount of EVD or SVD computation load, which are required in conventional DOA estimation algorithms such as the MUSIC or ESPRIT algorithms. All simulation results validate the superiority of the proposed TFOC -OPRM localization algorithm. Moreover, the simulation results indicate that the proposed TFOC -OPRM localization algorithm achieves lower computational complexity and better accuracy than the FOC -OPM algorithm and the MFOC -OPM algorithm both in spatially white noise and spatially color noise situations.