Efﬁcient Cumulant-Based Methods for Joint Angle and Frequency Estimation Using Spatial-Temporal Smoothing

: Most non-Gaussian signals in wireless communication array systems contain temporal correlation under a high sampling rate, which can offer more accurate direction of arrival (DOA) and frequency estimates and a larger identiﬁability. However, in practice, the estimation performance may severely degrade in coloured noise environments. To tackle this issue, we propose real-valued joint angle and frequency estimation (JAFE) algorithms for non-Gaussian signals using fourth-order cumulants. By exploiting the temporal correlation embedded in signals, a series of augmented cumulant matrices is constructed. For independent signals, the DOA and frequency estimates can be obtained, respectively, by leveraging a dual rotational invariance property. For coherent signals, the dual rotational invariance is constructed to estimate the generalized steering vectors, which associates the coherent signals into different groups. Then, the coherent signals in each group can be resolved by performing the forward-backward spatial smoothing. The proposed schemes not only improve the estimation accuracy, but also resolve many more signals than sensors. Besides, it is computationally efﬁcient since it performs the estimation by the polynomial rooting in the real number ﬁeld. Simulation results demonstrate the superiorities of the proposed estimator to its state-of-the-art counterparts on identiﬁability, estimation accuracy and robustness, especially for coherent signals.


Introduction
Sensor arrays have been used in radar, sonar, wireless communications, seismology, etc. Joint direction of arrival (DOA) and frequency estimation of incident signals have been attracting considerable attention due to their application to various fields. For instance, the two parameters can be applied to locate the mobile subscribers and allocate pilot tones in space division multiple access systems in wireless communications. Besides, an accurate estimation of these parameters is helpful with channel estimation and thus enhances the system performance.
The maximum likelihood estimator (MLE) can achieve the optimal estimates, but is computationally cumbersome as it inevitably relies on a multidimensional search [1]. To reduce computational expenditure, numerous subspace estimators based on the linear algebraic theory have been proposed to extract the information efficiently from a batch of observations. In [2], both parameters were obtained by performing a two-dimensional (2D) spectral search of the 1. The most important difference is that the operation of FOC with respect to JAFE is developed for the first time, and the corresponding existence conditions, and the number of achieved DOFs is analysed accordingly. 2. In [6][7][8][9][10][11], only the temporal and spatial white Gaussian noise is considered, while in this paper, we study the JAFE in a more general case, i.e., the temporal and spatial coloured Gaussian noise, and propose an FOC-based estimator to jointly determine DOA, as well as frequency robust to white or coloured noise. 3. In [6,7,9], the number of coherent signals that can be dealt with is up to 2M 3 , where M is the number of array elements, which is a strict condition. In practice, one base station serves multiple users simultaneously, where the propagation from each user composes multipath, so the total number of coherent signals may exceed that of array elements. The approach developed in this paper is capable of handling such a tough challenge that cannot be overcome by the previous techniques. 4. To reduce the cumbersome computations arising from the processing by FOC, we transform the complex-valued information to the field of real numbers and perform the resultant RARE estimator by solving the roots of a polynomial, avoiding the multidimensional search.
The remainder of this paper is organized as follows. In Section 2, an array model for signals corrupted by unknown coloured noise is introduced. In Section 3, the RARE estimators for JAFE of uncorrelated and coherent signals using fourth-order cumulants are developed. Section 4 provides numerical examples for demonstrating the validity and efficiency of our proposed algorithms. Finally, some concluding remarks are given in Section 5.
Throughout this paper, the following notations will be used: the operators (·) T , (·) * , (·) H , (·) + , cum{·}, det{·}, ⊗, •, and · 2 denote the operation of the transpose, conjugate, conjugate transpose, pseudo-inverse, cumulant, determinant, Kronecker product, Khatri-Rao (KR) product and Euclidean ( 2 ) norm, respectively. The symbol diag{z 1 , · · · , z N } represents a diagonal matrix with diagonal entries z 1 , · · · , z N , and blkdiag{Z 1 , Z 2 } represents a block diagonal matrix with diagonal entries Z 1 , Z 2 . The symbol I K stands for an identity matrix of size K × K. The symbol Z(a : b, c : d) refers to a constructed submatrix by the entries from a to the bth row and c to the dth column of Z, and the symbol Z(a, b) denotes the entry in the ath row and the bth column of Z.

Signal Model
Consider Q narrowband non-Gaussian signals impinging on a uniform linear array (ULA) with M identical omnidirectional sensors. Let the signal from the ith far-field source be denoted as e j2π( f c + f i )t s i (t), where f c is the carrier frequency, f i is a small frequency offset (central frequency) for the ith signal with f i f c and s i (t) is the amplitude of the ith signal, for i = 1, 2, · · · , Q. After demodulation to intermediate frequency (IF), the signal due to the ith source becomes e j2π f i t s i (t). We assume that there are K groups of coherent signals, which come from statistically independent far-field sources {s k (t)} K k=1 . In the kth coherent group, the signals coming from direction θ kp , p = 1, 2, · · · , P k , correspond to the pth multipath propagation of the source s k (t), and the complex fading coefficient is α kp . It is apparent that the total number of coherent signals is Q = ∑ K k=1 P k . Since the K sources are mutually independent, the signals between coherent groups are also independent of each other, while the signals in each group are coherent. It is worth noting that if the independent sources arrive at the sensor array undergoing a line of sight channel (no multipath propagation), the incident Q signals can be regarded as Q coherent groups each of which only includes one signal, i.e., K = Q and P 1 = P 2 = · · · = P k = 1. The noise is spatially coloured and statistically independent of the signals. Let P be the sampling rate that is much higher than the largest IF, then the M × 1 array output vector is given by [7]: where a(θ) = 1, β(θ), · · · , β M−1 (θ) T ∈ C M is the steering vector, β(θ) = e j 2πd λ sin θ with λ and d being the wavelength of the carrier signal and the spacing between adjacent elements, respectively, A = a(θ 1 ), · · · , a(θ Q ) is the array manifold, Γ = blkdiag{α 1 , · · · , α K } with α k = α k1 , · · · , α kP k T containing attenuation information of the kth coherent group, while a special case Γ = I Q if the incident signals are independent of each other, Φ = diag e j2π f 1 P , · · · , e j2π f K P contains the information of the IF frequencies of incident signals, s( n P ) = s 1 ( n P ), · · · , s K ( n P ) T ∈ C K and w( n P ) is the spatially-coloured Gaussian noise. Besides, we assume that the array manifold A is unambiguous, i.e., the steering vectors {a(θ i )} Q i=1 are linearly independent for any set of distinct {θ i } Q i=1 . Collecting N samples of the array output x( n P ) at a rate P and stacking them into a row, one has: where W ∈ C M×N collects N samples of the noise vector w( n P ).

Temporal Smoothing
In this section, we consider a data stacking technique, referred to as temporal smoothing, which adds structure to the data model for the implementation of the JAFE algorithm. An m-factor temporally-smoothed data matrix is constructed by stacking temporally-shifted versions of the original data matrix. This results in the following matrix of size mM × (N − m + 1) [7]: where W m represents the noise term constructed from W in a similar manner as X m is obtained from X. Assume that the signals are narrowband, i.e., s(t) ≈ s(t + 1 P ) ≈ · · · ≈ s(t + m−1 P ). In this case, all the block rows in the right-hand term of (3) are approximately equal, which means that X m has the following factorization: where A m , referred to as the extended array manifold, is given by: and: is a matrix collecting N − m + 1 samples of the K sources.

Independent Signals
First, we consider the case of independent signals, where K = Q, P 1 = P 2 = · · · = P k = 1 and Γ = I Q , so Γ in (1)-(5) can be neglected. Considering the received signals are assumed to be non-Gaussian, one can establish the array FOC matrix between the received data blocks as: where y(n) = H q X m (:, n) and z(n) = H 1 X m (:, n). Herein, H q = 0 M×M(q−1) , I M , 0 M×M(m−q) , q = 1, 2, · · · , m, and X m (:, n) is the nth column of X m . The entries of C q in the [(k 1 − 1) M + k 2 ] th row and the [(l 1 − 1) M + l 2 ] th column are defined as: where y i (n) and z i (n) are the ith entries of y(n) and z(n), respectively. Substituting (3) into (7) and utilizing the properties of cumulants [CP1]-[CP5]in [18], one can further get: where F s,i is the entry in the ith entry of F s (: By the principle of aperture extension via FOC [18], some rows of A • A * ∈ C M 2 ×Q 2 are redundant with respect to the rest in the case of a ULA, resulting in the effective array aperture being 2M − 1 elements. Considering the regular distribution of the redundancy in A • A * , one can reduce the problem dimension by performing a linear transformation as: where: and To take full advantage of the information embedded in C q m q=1 , one can design a new cumulant matrix, of size m(2M − 1) × m(2M − 1), as: It can be readily verified that: whereJ ∈ R m(2M−1)×m(2M−1) denotes an exchange matrix that has unity entries on the cross-diagonal and zeros elsewhere. It is apparent thatJC * xJ = C x , i.e., C x is centro-Hermitian, as ideally, C s is a real diagonal matrix. However, to "double" the number of snapshots and decorrelate possible sample correlation in an arbitrary matrix C s due to a limited observation window, the centro-Hermitian property is sometimes forced by means of the so-called forward-backward (FB) averaging, that is, Then, one can transformC x into the field of real numbers for the sake of computational efficiency. We define a new matrix C r as: where Q u , u = m (2M − 1), is a sparse unitary matrix given by: where I v ∈ R v×v and J v ∈ R v×v are the identity matrix and exchange matrix, respectively, and v is an arbitrary natural number. Referring to Theorem 3 in [21], one can identify that C r is a real matrix.
Performing the singular-value decomposition (SVD) of C r , one has: where The columns of U s U(:, 1 : Q) are the singular vectors corresponding to the Q largest eigenvalues, while the columns of U n U (:, Q + 1 : m (2M − 1)) are the singular vectors corresponding to the m (2M − 1) − Q singular values, which are all zeros. The signal subspace is spanned by the columns of U s , whereas the noise subspace is spanned by , one can construct the following function: where We now show how to estimate DOAs based on the determinant or the smallest eigenvalue of M(θ). It can be found that the size of in general, is of full row rank, and M(θ) is of full rank. However, when θ coincides with any one of the Q desired DOAs, i.e., θ = θ i , i = 1, · · · , Q, the expression in (18) becomes zero by the orthogonality between the signal and noise subspaces. Since g ( f ) = 0, (18) can hold true only if the matrix M(θ) is rank deficient or, equivalently, its determinant (as well as its smallest eigenvalue) is equal to zero. Now, it becomes clear that the determinant of M(θ) can be utilized to estimate DOAs as:θ One solution to (19) is to perform a one-dimensional spectral search to determine the DOA estimates corresponding to the minima of det {M(θ)}, which approximates to zero. However, it involves more computational cost than JAFE, which exploits ESPRIT without any grid search. If the value of m is small, which is the common case in practice [7][8][9], one can turn (19) into a polynomial rooting to avoid the search process and resolve DOA estimates directly. Denote z = β(θ), . As a result, a root-RARE polynomial can be constructed as: where r k is the coefficient of the polynomial. Therefore, the roots of (20) {z i } Q i=1 that are closest to the unit circle indicate the desired DOAs of independent signals. The DOA estimates related to the polynomial roots can be determined as: Once the DOA estimates are obtained, one can perform the frequency estimation as: Equation (22) is a quadratic minimization problem. Rather than spectral grid search, we provide a closed-form solution to frequency estimation that is more computationally efficient. Let e 1 = [1, 0, · · · , 0] T ∈ R m , then one has the constraint e T 1 g ( f ) = 1, which can remove the trivial solution g ( f ) = 0, and (22) can be rewritten as: By using the Lagrange multiplier µ, one can construct the following cost function [22]: Setting the partial derivative of L ( f ) to zero, i.e., we find that when To estimate the frequency f i , one can extract the phases from g ( f i ) as: Since the solution of f i in (28) is overdetermined, one can construct a fitting relationship to find out the desired frequency in the least squares sense, that is, with c i being an augmented parameter, and: The least squares solution of (29) isĉ i = ĉ i ,f i T = P + h i = P T P −1 P T h i , and one can directly obtainf i fromĉ i accordingly.

Coherent Signals
As emphasized earlier, the above method is restricted to the case of independent signals, and its performance would deteriorate severely provided that there exists coherence between signals due to the multipath propagation. In this subsection, another FOC-based estimator is proposed to tackle coherent signals by extending a dual rotation-invariant property to blind signal separation. More specifically, the dual rotation-invariant property in terms of frequency and spatial signature is first utilized to estimate the frequencies and blindly separate the coherent groups. Then, the DOAs of coherent signals in each group can be estimated after rank recovery by FBSS. DenotingĀ = AΓ = [ā 1 , · · · ,ā K ], referred to as the generalized array manifold, whereā k = [a(θ k1 ), · · · , a(θ kPk )] α k is the generalized steering vector, we first construct a series of cumulant matrices, distinct from (7), as follows: n e j2π f n q−1 P F s,n , n e j2π f n q−1 P F s,n , where q = 1, 2, · · · , m,C s = diag γ 1 |Ā(1, 1)| 2 , · · · , γ K |Ā(1, K)| 2 , and D = diag Ā (2,1) A(1,1) , · · · ,Ā (2,K) A(1,K) . Then, one can build three cumulant matrices, of size mM × nM, with these blocks as: .
in an appropriate manner, one has: It is readily identified thatÃ 2 =Ã 1 Φ. Then, one can perform the eigen-decomposition of Υ 2 Υ + 1 to obtain the frequencies, as well as the scaled column vectors ofÃ 1 by the following proposition. Proposition 1. The principal eigenvectors and eigenvalues of the generalized space-frequency matrix Υ 2 Υ + 1 are the generalized steering vectors ofÃ 1 scaled by nonzero coefficients and Φ, respectively, i.e.,Ȗ s =Ã 1 Λ andΣ s = Φ, where the columns ofȖ s ∈ C mM×K are the K eigenvectors, corresponding to the K largest eigenvalues of Υ 2 Υ + 1 that are contained in the diagonal matrixΣ s , and Λ is an arbitrary diagonal matrix with nonzero entries.
Referring to [23], one knows that rank {Z k } = P k after performing FBSS, which implies that the rank deficiency due to the coherence has been successfully recovered. As a result, by applying prevailing high-resolution DOA estimation methods, such as MUSIC or ESPRIT, toZ k , one can get the DOA estimates of the coherent signals. Note that JZ * k J =Z k , i.e.,Z k is the centro-Hermitian matrix, and one can perform the following linear transformation to relieve the computational load as: where Q l is defined in (16). Referring to Theorem 3 in [21], T k is a real matrix, and thus, one can resort to unitary ESPRIT [24] to resolve the DOA estimates of the coherent signals. As expatiated above, the frequency, generalized steering vectors and DOAs are estimated by using the rotation-invariant property, so the proposed JAFE estimation method for coherent signals can be referred to as the generalized ESPRIT-like estimator.

Remark 1.
Since the necessary condition for (19) being valid is m ≤ m (2M − 1) − Q, the maximum number of resolvable independent signals by our root-RARE polynomial estimator is 2m(M − 1), which is twice as much as JAFE. This is because the virtual array manifold B resulting from FOC doubles the DOFs, as well as the effective aperture. By the principle introduced in [23], up to 2M 3 coherent signals in the kth group can be resolved, provided that l ≥ P k + 1 and 2r ≥ P k . Therefore, considering at most M − 1 coherent groups can be separated, the identifiability of the proposed generalized ESPRIT-like estimator is upper bounded by (M − 1) 2M 3 .

Remark 2.
The cumulant requires many more computations than the covariance and, hence, the proposed FOC-based estimator works well on the condition that the observation window is sufficiently long and the target is (quasi-) stationary. We attempted to implement the proposed method in a test bed composed of a Virtex-7 series FPGA and a TMS320C6x series DSP, but found that it is infeasible to complete the joint DOA and frequency estimation within time less than the scale of milliseconds in the case of M = 5, m = 2 and N = 2000 due to the cumbersome calculations of cumulants and the eigen-decomposition of a matrix, of size 18 × 18. As a result, for some applications, like high-speed moving targets, the coherent processing interval of such cumulants is quite short, and it is hard and even unlikely to realize the proposed algorithm by the prevailing hardware available on the market within the period.

Simulation Results and Discussion
In this section, a series of numerical experiments under different conditions is carried out to examine the performance of the proposed algorithms with some state-of-the-art methods. Simulations were conducted for a small number of ULA with half-wavelength spacing between adjacent elements. The existing approaches MR-ESPRIT and U-JAFE and the Cramér-Rao lower bound (CRLB) of joint DOA and frequency estimates [7] were chosen for comparison. For simplicity, we assumed that all signals had an identical power. Similar to the settings in [20,25], the signals were modelled as s(t) = F(t)r(t), where F(t) = diag { f 1 (t), · · · , f K (t)}, r(t) = [r 1 (t), · · · , r K (t)] T with f i (t) and r i (t) being zero-mean Gaussian processes with unit-variance and σ 2 s -variance, respectively. The noise was assumed to be a first-order spatial autoregressive process, and the (a, b)th entry of the noise covariance matrix is given by R(a, b) = σ 2 n 0.8 |a−b| e j π(a−b) 16 [26,27]. The sampling rate P was fixed at 30 MHz. The accuracy of the DOA estimate, the statistical performance of the algorithms, was measured from 800 Monte Carlo runs in terms of the root mean squared error (RMSE), which is defined as: is the estimate of θ i for the nth trial and Q is the number of all signals.

JAFE of Independent Signals
In the first scenario, we considered that three independent sources from [−52 • , −27 • , −10 • ] impinge on a five-element array, and the corresponding small frequency offsets (central frequencies) for the signals were 1.2 MHz, 1.6 MHz and 1.8 MHz, respectively, while the temporal smoothing factor m = 3. The resultant RMSE as a function of SNR and the number of snapshots N in the first scenario is illustrated in Figure 1. It can be observed that the proposed FOC-based method outperformed its two counterparts up to 10 dB, while being inferior to them when the SNR became even larger. For the other two schemes, U-JAFE performed slightly better than MR-ESPRIT at low SNRs and the RMSEs of both techniques tended to be the same as SNR increased. In addition, the performance of the proposed estimator improved with the increase of the number of snapshots, whereas the RMSEs of both U-JAFE and MR-ESPRIT stabilized for all snapshot sizes, even though as many as 10,000 snapshots were collected, implying that the estimation performance of these two algorithms was mainly subject to the SNR while immune to the large number of snapshots. Figure 2 depicts the RMSEs of the frequency estimates varying with changes of the SNR and snapshot size. Clearly, the proposed approach performed the best when the SNR was less than 7 dB or the number of snapshots exceeded 1500, whereas both U-JAFE and MR-ESPRIT were superior to our solution at moderate to high SNRs, and U-JAFE still achieved somewhat better accuracy of estimates than MR-ESPRIT under low SNRs, similar to the case of DOA estimation. Generally speaking, the proposed FOC-based method had a certain advantage over its second-order statistics counterparts in low SNR environments, whereas the latter provided better performance at high SNRs, which is consistent with the general results of the existing work. It should be noted that clear margins arose between the estimates, both DOA and frequency, and the CRLBs. This is mainly due to the fact that the cumulant-based estimator is inherently biased, and algorithms using second-order correlation are vulnerable to coloured noise, which will cause large estimation errors. As discussed earlier, MR-ESPRIT, U-JAFE and the proposed FOC-based method were able to identify more independent signals than physical sensors, so it is necessary to investigate the performance of three approaches in such an underdetermined case. Assume that six independent sources from [−52 • , −27 • , −10 • , 10 • , 27 • , 47 • , ] impinge on a five-element array (Q ≥ M); the corresponding small frequency offsets for the signals are 1.2 MHz, 1.6 MHz, 1.8 MHz, 2 MHz, 2.5 MHz, 3 MHz, respectively, and the temporal smoothing factor m = 3. It can bee seen from Figure 3 that the absolute estimation performance of all techniques became worse, which is also illustrated by the CRLB, while the advantages of our technique over the other two became more obvious than the overdetermined case, and the counterparts U-JAFE and MR-ESPRIT only worked well above 19 dB and 20 dB, very large SNRs, respectively. Besides, the RMSE of the proposed FOC estimator approached the CRLB as the snapshot size increased, while both U-JAFE and MR-ESPRIT obtained almost the same large errors and remained constant throughout the numbers of snapshots, similar to Figure 1b.  Figure 4 exhibits the performance of the three methods for the frequency estimation, as well as CRLB as a function of SNR and the total number of snapshots. One can observe that generally, the proposed method was distinctly better than the other two algorithms except at 20 dB and above, and U-JAFE outperformed MR-ESPRIT for all SNRs and snapshot sizes, which agrees with the results in Figure 2.

JAFE of Coherent Signals
In this section, we investigate the JAFE performance of the proposed cumulant-based algorithm first in the underdetermined coherent signals case, i.e., Q ≥ M. It is assumed that two groups of six coherent signals uniformly distributed between −60 • and 60 • are imping on a six-element ULA, and the frequency offsets are 1.2 MHz and 1.6 MHz, respectively, as each group of coherent signals was derived from an independent source. The fading amplitudes of the coherent signals were [1,  ], respectively. The temporal and spatial smoothing factors were m = 2 and r = 2, and the SNR and number of snapshots were fixed at 15 dB and 5000. In such a scenario, neither U-JAFE nor MR-ESPRIT was able to resolve all DOAs as their DOFs were both upper bounded by four, whereas the proposed method can handle 20 DOAs at most. Figure 5a depicts the DOA estimation results of 50 independent trials by the proposed method. The solid lines mark the true DOAs. It can be seen that each DOA was correctly determined, and the frequency offset estimates were even more accurate, as shown in Figure 5b, since the DOFs used for the frequency estimation were larger than that for the DOA estimation.
Next, we examined the estimation performance of the proposed FOC-based estimator in the presence of multipath propagation. Assume that one group of coherent signals from [−52 • , −27 • , −10 • ] arrives at a six-element array; the corresponding frequency offsets for the signals are 1. ], respectively. The temporal and spatial smoothing factors were m = 2 and r = 2, respectively. Figure 6 shows the RMSE of DOA estimates in the overdetermined case. It can be found that as the SNR and the number of snapshots increased, the RMSE of DOA estimation decreased gradually for all of the methods, where the accuracy of our algorithm improved asymptotically, whereas the RMSEs of U-JAFE and MR-ESPRIT were close to each other at low and moderate SNRs and tended to saturate at 1 • and 0.7 • , respectively. Evidently, the proposed FOC-based estimator was significantly superior to the counterparts for all SNRs and snapshot sizes tested, mainly because it exploited the dual rotational invariance of the frequency constructed by the rotation-invariant property in the generalized array manifold. In Figure 7, we plot the RMSEs of frequency estimates of the three methods as a function of SNR and the total number of snapshots under multipath propagation. It can be seen that the proposed solution had a prominent advantage over the other two algorithms for almost all the SNRs and snapshot sizes, although its lead kept decreasing as the SNR rose and turned out to be lost at 19 dB and above.Unlike the illustration in Figure 2, the performance of frequency estimation of U-JAFE was strictly worse than MR-ESPRIT over the whole range of SNR and snapshot size for coherent signal estimation. Additionally, the RMSEs of the frequency estimates by our method approached the CRLBs more tightly than the results in Figures 2 and 4 since all three coherent signals shared an identical frequency, and the number of frequencies to be estimated turned out to be one less than the cases studied above.

Conclusions
The joint angle and frequency estimation of non-Gaussian signals in coloured noise is addressed in this paper, and cumulant-based estimation algorithms are developed for independent and coherent signals, respectively. The proposed methods make efficient use of the dual rotational invariance of the newly-constructed array manifolds to enhance the DOF of a ULA, thus improving the estimation accuracy and resolving more signals than the number of array elements, either for the independent or coherent signal. For independent signals, we construct an augmented cumulant matrix, providing twice the DOFs of the covariance matrix formed by the state-of-the-art techniques, then transform it to the real number field and resolve the DOA estimates by using rooting of a polynomial resulting from the RARE property, and the frequency estimates can be resolved in a least squares manner given the DOA estimates already acquired, avoiding the spatial spectra search and reducing the consumption of calculations accordingly. For coherent signals, the generalized steering vectors (coherent groups) and the frequencies are blindly separated from each other by exploiting a new dual rotational invariance arising from a series of dedicated cumulant matrices. Then, FBSS is performed to restore the inherent rank deficiency of the generalized steering vectors, and finally, the DOA estimates are determined by using the unitary ESPRIT algorithm. Since the coherent signals are divided and conquered, into their associated group, more coherent signals than array elements can be handled. Compared with previous second-order statistics methods, the proposed solutions are efficient in the sense that by extracting the temporal correlation embedded in the high-order statistics, they further raise the array DOFs and effective aperture, as well as improving the estimation accuracy, especially for coherent signals, in addition to the robustness to coloured noise.
Author Contributions: All authors contributed extensively to the study presented in this manuscript. Y.W. designed the main idea, methods and experiments, interpreted the results and wrote the paper. L.W. and B.W.-H.N. supervised the main idea, edited the manuscript and provided many valuable suggestions to this study. X.Y., J.X. and P.Z. carried out the experiments.