Cumulant-Based DOA Estimation of Noncircular Signals against Unknown Mutual Coupling

To effectively find the direction of non-circular signals received by a uniform linear array (ULA) in the presence of non-negligible perturbations between array elements, i.e., mutual coupling, in colored noise, a direction of arrival (DOA) estimation approach in the context of high order statistics is proposed in this correspondence. Exploiting the non-circularity hidden behind a certain class of wireless communication signals to build up an augmented cumulant matrix, and carrying out a reformulation of the distorted steering vector to extract the angular information from the unknown mutual coupling, by exploiting the characteristic of mutual coupling, i.e., a limited operating range and an inverse relation of coupling effects to interspace, we develop a MUSIC-like estimator based on the rank-reduction (RARE) technique to directly determine directions of incident signals without mutual coupling compensation. Besides, we provide a solution to the problem of coherency between signals and mutual coupling between sensors co-existing, by selecting a middle sub-array to mitigate the undesirable effects and exploiting the rotation-invariant property to blindly separate the coherent signals into different groups to enhance the degrees of freedom. Compared with the existing robust DOA methods to the unknown mutual coupling under the framework of fourth-order cumulants (FOC), our work takes advantage of the larger virtual array and is able to resolve more signals due to greater degrees of freedom. Additionally, as the effective aperture is virtually extended, the developed estimator can achieve better performance under scenarios with high degree of mutual coupling between two sensors. Simulation results demonstrate the validity and efficiency of the proposed method.


Introduction
Direction of arrival (DOA) estimation, an important research area of sensor array signal processing, has attracted a large amount of attention because of its wide applications to electromagnetic, acoustic, seismic sensing, etc [1][2][3][4][5][6]. However, non-negligible interaction between sensors, i.e., mutual coupling, makes prevailing DOA estimation algorithms invalid since extra unknowns are introduced. Weiss and Friedlander [7] first discuss the structure of the mutual coupling matrix (MCM) in uniform linear and circular arrays and estimated the DOAs and mutual coupling coefficients in an iterative way. This work is followed by Sellone et al. whose approach is more robust to the detrimental effects and does not need a preliminary estimate [8]. Both schemes suffer from high computational complexity, which is inevitable given their reliance on multidimensional search processes. To circumvent the adverse effects of mutual coupling, the technique in [9,10] uses a subset of auxiliary sensors to make the "middle sub-array" mutual coupling free, and prevalent super-resolution algorithms array that is coupling-free, and then perform a DOA estimation in each coherent group via the forward-backward spatial smoothing (FBSS) technique. The estimation performance in terms of identifiability and computational complexity is analyzed in Section 4. Simulation results in Section 6 show the significant improvement of the proposed method. Finally, some concluding remarks are given in Section 7. Throughout this submission, the following notations will be adopted: the operators (·) T , (·) * , (·) H , (·) + , | · |, E[·], det{·}, ⊗, •, and · 2 denote the operation of transpose, conjugate, conjugate transpose, pseudo-inverse, modulus, expectation, 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 , blkdiag{Z 1 , Z 2 } symbolizes a block diagonal matrix with diagonal entries Z 1 and Z 2 , Toeplitz {·} represents a symmetric Toeplitz matrix constructed by a vector in the brace, and cum{z 1 , z 2 , z 3 , z 4 } symbolizes a cumulant calculated from data z 1 , z 2 , z 3 , and z 4 . The symbol I K stands for an identity matrix of size K × K. The symbol Z(a : b, c : d) refers to a constructed sub-matrix by the entries from a to b-th row and c to d-th column of Z, and the symbol Z(a, b) denotes the entry in the a-th row and b-th column of Z.

Strictly Second-Order Non-Circular Signals
For a signal s(t), if its elliptic covariance is E[s 2 (t)] = ρe jφ E[|s(t)| 2 ] = ρe jφ σ 2 s = 0, where φ is the deterministic non-circularity phase, and ρ is non-circularity rate satisfying 0 < ρ ≤ 1, then we say s(t) is non-circular. The special case ρ = 1 is common in wireless communications, and the received baseband signals satisfying such a condition is referred to as strictly second-order non-circular, like amplitude modulation (AM), binary phase shift keying (BPSK), amplitude shift keying (ASK), offset quadrature phase shift keying (OQPSK), and pulse amplitude modulation (PAM), etc. [25][26][27][28]30,38,39]. To aid the understanding of such a particular class of signals, we take the baseband signal of BPSK or ASK modulation as an example. In the framework of DOA estimation, the analog received BPSK or ASK modulated signals are bandpass filtered and down converted to baseband, the in-phase and quadrature components are matched-filtered, sampled, and paired to obtain complex signal s(t) that can be factorized into s(t) = e jϕ s 0 (t) where s 0 (t) is a real-valued symbol and ϕ is an arbitrary phase shift, due to the initial phase of the transmitted carrier, that can be different for each signal but constant with time. Therefore, the baseband signal s(t) at the receiver side is complex-valued. This model has been used by numerous authors who have studied DOA estimation of non-circular signals [25][26][27][28]30,38,39]. Examining the covariance and elliptic covariance of s(t), one gets E [s(t)s * (t)] = E e jϕ s 0 (t)e −jϕ s 0 (t) = E[s 2 0 (t)] = σ 2 s and E [s(t)s(t)] = E e jϕ s 0 (t)e jϕ s 0 (t) = e j2ϕ E[s 2 0 (t)] = e j2ϕ σ 2 s = 0. By the definition of non-circular signal, one can readily find that ρ = 1 and φ = 2ϕ, i.e., the baseband signal of BPSK or ASK modulation is strictly non-circular, and the non-circularity phase is twice as much as the initial phase. For the case of the received baseband signals with unknown modulation, one can resort to two kinds of classification approaches to identify the modulation of the signals. One class is the likelihood-based approach [40][41][42] and the other is the so-called feature-based approach [43][44][45][46][47][48]. If the signals are recognized as non-circular, such as BPSK, ASK, OQPSK, or PAM, the proposed method in the sequel is applicable, otherwise it fails to work. The research on modulation classification is beyond the scope of this manuscript, so we assume that the modulation of the received baseband signals is known a priori or can be recognized by using existing classification approaches if not specified, and the signals to be dealt with are strictly second-order non-circular.

Array Model for Non-Circular Signals
Consider narrowband non-circular signals s 0,i (t) in the far filed impinge on an M-elements uniform linear array (ULA) from N different angles θ i , i = 1, 2, · · · , N, and the mutual coupling between sensors cannot be neglected, then at the time index t the corresponding array observation can be expressed as where a(θ) = 1, β(θ), · · · , β M−1 (θ) T ∈ C M is the steering vector, β(θ) = e j 2πd λ sin θ , λ, and d are the carrier wavelength and the spacing between adjacent sensors, respectively, C denotes the MCM quantifying the degrees of electromagnetic coupling among elements, A = a(θ 1 ), · · · , a(θ N ) is the T ∈ R N , and n(t) is the noise following Gaussian distribution but spatially colored. In addition, it is assumed that A is unambiguous, i.e., the steering vectors {a(θ i )} N i=1 are linearly independent for any set of distinct {θ i } N i=1 . As described in [9][10][11][49][50][51], one can sufficiently model the perturbed ULA by taking the electromagnetic characteristic of mutual coupling, such as a limited operating range and an inverse relation of coupling effects to interspace, into consideration. To be specific, If the mutual coupling acting on two sensors covers P interelement spacing at most such that the resultant MCM has finite non-zero coefficients, then the resultant MCM is banded symmetric Toeplitz, which can be formulated as where 0 < |c 1 |, |c 2 |, · · · , |c P−1 | < c 0 = 1 are the mutual coupling coefficients.

Proposed Non-Circular FOC-Based Estimator
To increase the effective aperture as well as the DOFs, one can design an augmented cumulant matrix, of size 2M 2 × 2M 2 , as

DOA Estimation Without Mutual Coupling Compensation
The augmented matrix C x has the following factorization through the singular value decomposition (SVD) where Σ = diag {λ 1 , · · · , λ 2M 2 } is composed of 2M 2 singular values that can be sorted in descending order as λ 1 ≥ · · · ≥ λ N > λ N+1 = · · · = λ 2M 2 = 0. The matrix U s U(:, 1 : N) collects the singular vectors corresponding to the N largest singular values, while U n U(:, N + 1 : 2M 2 ) contains the rest of the singular vectors corresponding to the 2M 2 − N zero singular values. It is known that U s spans the signal subspace, while its orthogonal space, namely the noise subspace, is spanned by U n , thus one has However, due to the unknown C one cannot directly obtain the DOA estimates of the desired non-circular signals via 1-D search for making Equation (11) hold. For the purpose of decoupling the angular information from the electromagnetic impact between sensors of the ULA, we reparameterize the distorted steering vector as [11] Ca(θ) = T(θ)α (12) where with and It should be noted that τ(θ) in Equation (14) is generally not equal to zero, otherwise "blind angles" cannot be identified [11]. Substituting Equation (12) back to Equation (11), one has generally speaking, then the matrixT H (θ)U n has a full row rank and M(θ) has a full rank accordingly. Nevertheless, for some special cases that θ is exactly the same with any one of the N true angles, i.e., θ = θ i , i = 1, · · · , N, the expression in Equation (21) becomes zero. Since υ = 0, Equation (21) is valid only when M(θ) is singular that it takes as a zero-value determinant. As a result, the DOA estimation is now dependent on the determinant of M(θ). Since M(θ) does not contain any mutual coupling coefficients, det {M(θ)} is insensitive to the nuisances. Therefore, one can estimate the bearings aŝ

Performance Analysis
In this section, the identifiability as well as the computational complexity of the proposed method is discussed, with a comparison of the existing work in the context of FOC, e.g., Liu's [21] and Liao's [24] approaches.

Identifiability of DOA Estimation
As 2(2P − 1) 2 ≤ 2M 2 − N makes M(θ) be of full rank, which is a necessary condition of the proposed algorithm, one can deduce that N ≤ 2M 2 − 2(2P − 1) 2 , which implies that the maximum DOFs of the developed estimator are 2M 2 − 2(2P − 1) 2 . By contrast, existing methods cannot estimate as many sources as our technique. Liao's method is based on ESPRIT and as a result, the upper bound of estimation identifiability depends on the dimension of signal subspace. It is known that the signal subspace produced by Liao's method has a size of M(M − 2P + 1) × N, and there are linearly dependent rows due to the property of the middle sub-array. Therefore, the row dimension of the signal subspace is 2P(M − 2P + 1), in other words, Liao's method can estimate up to 2P(M − 2P + 1) signals. In Liu's method, only the middle sub-array utilized, which means the extended aperture by FOC is 2(M − 2P + 2) − 1 elements. Therefore, Liu's method can handle at most 2(M − 2P + 1) DOAs. It is evident that the proposed method has the largest estimation capacity among the three algorithms, followed by Liao's method and then Liu's method. Consequently, Liu's method suffers from the reduced effective aperture. Besides, if the coupling between sensors is prominent, i.e., P is large, or equivalently M is not sufficiently large, then the size of the middle sub-array will be very small. For instance, if P = M 2 and N = 4, Liu's method is not valid. This is because the effective array aperture [18] is 3 elements in this case, and at most 2 DOAs can be estimated.

Computational Complexity
Next, we compare the major computations of the proposed estimator, Liao's and Liu's methods, which involve in the fourth-order statistical matrices construction, SVD, and spectral search. For Liao's approach, the main calculation cost is to form M − 2P + 2 FOC matrices of the same size M × M and to perform the SVD of the augmented flops in total where δ is the sampling grid spacing. By comparison, the proposed scheme establishes two M 2 × M 2 cumulant matrices, carries out SVDs of the augmented 2M 2 × 2M 2 matrix, and executes a one-dimensional spectral search on determinants of a 2 (2P − 1) 2 matrix. The resulting flops taken in our work are in order of O 18M 4 L + 8M 6 + 1440 δ (2P − 1) 6 . As in general L is relatively large, say several thousands, the proposed method has the highest computational complexity among the three algorithms, followed by Liao's approach and then Liu's technique.

Solution to the Case of Coherent Signals
It should be noted that the aforementioned method is applicable to incident signals being independent, and the coherency between signals may lead to severe performance degradation or even total failure due to the rank deficiency. This motivates us to provide a solution to the knotty problem, i.e., multiple groups of coherent signals received by a tightly coupled ULA in colored noise, and another direction finding algorithm based on FOC is proposed in this section, exploiting rotation-invariant property to blindly separate the coherent signals into different groups, and then resolving the DOAs in each coherent group via rank restoration techniques such as FBSS.
LetÃ =ĀΓ = [b 1 , · · · , b K ], being the so-called generalized array manifold, where b k = [a(θ k1 ), · · · , a(θ kPk )]α k is the generalized steering vector, one can get four cumulant matrices as follows (1, m)e j φm 2 s 0,m (t), b n e j φn 2 s 0,n (t), (1, m)e j φm 2 s 0,m (t), b n e j φn 2 s 0,n (t), wherex 1 (t) andx 2 (t) are the first and second columns ofx(t), respectively. Then, one can take the eigen-decomposition to Ξ 2 Ξ + 1 and Ξ 4 Ξ + 3 respectively to work out the scaled generalized steering vectors. To see how it works, one can perform the following operations: where in the third equality we have used the identityÃ +Ã = I K sinceÃ is of full column rank. Through eigen-decomposition, one has Ξ 2 Ξ + 1 =Ẽ sΣsẼ H s , whereẼ s contains the eigenvectors of Ξ 2 Ξ + 1 whilẽ Σ s is a diagonal matrix comprising eigenvalues. Equivalently, one can deduce that Ξ 2 Ξ + 1Ẽ s =Ẽ sΣs . As a result, it is natural to revealẼ s =ÃΛ, where Λ is an arbitrary diagonal matrix with nonzero entries, andΣ s =D. Following the same principle, one can also obtain the eigenvectors of Ξ 4 Ξ + 3 as the scaled generalized steering vectors. Denotingũ k andȗ k are the the eigenvectors of Ξ 2 Ξ + 1 and Ξ 4 Ξ + 3 , respectively, one can get a more robust estimate of the scaled generalized steering vectors by averaging them, i.e.,b Then one can perform FBSS tob k to recover the rank deficiency where W i = 0 m×(i−1) , I m , 0 m×(q−i) is the selection matrix for the i-th sub-array, with m = M − 2P + 3 − q being is the number of sensors in each sub-array, and J ∈ R m×m is the exchange matrix. For the k-th group of coherent signals, one now can apply ESPRIT to B f b k to determine the DOA estimates.

Remark 1.
The parameter q plays a significant role in rank restoration since rank B f b k after FBSS becomes 1 + 2q, and it also restricts the DOFs and effective aperture for estimation by m = M − 2P + 3 − q. We consider two extreme cases to discuss the choices of q. If q = P i Remark 2. A cumulant requires many more computations than a covariance and, hence, the proposed FOC-based estimator of non-circular signals against mutual coupling works well on condition that the observation window is sufficiently long and the target is (quasi-) stationary. We attempted to implement the proposed method in a testbed composed of a Virtex-7 series FPGA and a TMS320C6x series DSP, but found that it is infeasible to complete the DOA estimation within a time less than the scale of milliseconds in the case of M = 8 and L = 2000 due to the cumbersome calculations of cumulants and the SVD of the augmented cumulant matrix C x , of size 16 × 16. As a result, for some applications, like maneuverable targets, the coherence time of such cumulants is quite short, and it is hard and even unlikely to realize the proposed algorithm by the prevailing hardware available in the market. On the other hand, for some applications that expect not much real-time quality, such as radio, hydrological, or meteorological environment monitoring where the incident signals have cyclostationarity, cumulant-based signal selective algorithms can be implemented for location estimation of far-field signals; the price to be paid is the need for the large number of computations and large data lengths for reliable estimation of the cumulants.

Remark 3.
For the case of wideband signals, one can divide the observations at each channel into some (possibly overlapping) segments, where for each segment, a number of frequency sub-bands are computed by the short-time Fourier transform (STFT). The idea of frequency-domain processing is to decouple the wideband model into a multitude of narrowband models, and then at each frequency subband one can apply the proposed method to obtain M( f , θ) constructed by Equation (21). In the frequency-domain approach, the final step is to combine M( f , θ) at various frequencies to obtain a DOA spectrum fusion, that is,θ = arg min 1 2 is the discrete normalized frequency band of the received signals. The combination in the spectrum fusion equation follows the principle of incoherent signal subspace method (ISSM) in [55]. By contrast, our solution is inapplicable to the focusing-based approaches [56,57], which is another category of DOA estimation algorithms for wideband signals, since initial DOA estimates required therein cannot be obtained by conventional beamforming in the presence of unknown mutual coupling.

Remark 4.
The array perturbations, such as mutual coupling, gain-phase errors, and sensor location uncertainties, can be calibrated by placing signal sources at known positions, which is the so-called active calibration [58][59][60][61][62][63][64][65]. However, the process of measuring the array manifold including various perturbations can be time consuming and expensive, and it is inconvenient or even infeasible to make the calibration source available in some cases. In addition to the deployment of signals of opportunity, keeping the calibration effective is another essential issue. Plenty of factors contribute to the variation of array response over time: gradual changes in the behaviors of the sensor itself, the electronic circuitry, and the analog-to-digital converter (due to thermal effects, aging of components, etc.), changes in the electromagnetic environment (e.g., metal objects beside an antenna array cause a distortion of the beam pattern), and changes in the sensor locations (e.g., an antenna array mounted on the vibrating wing of an aircraft or a hydrophone array towed behind a ship). The active calibration scheme has its inherent shortcomings that are difficult to be overcome and, hence, in this paper we resort to the self-calibration solution that is not reliant on the array manifold measurement as well as the calibration sources whose positions are known in advance.

Remark 5.
The received data can be prewhitened with a square-root of the covariance matrix of the colored noise, and then the covariance matrix of colored noise is reduced to an identity matrix accordingly, i.e., the colored noise becomes "isotropic". It should be noted that this operation relies on a prerequisite that the noise covariance matrix is available in advance. Generally speaking, noise in the internal circuitry of the antenna elements is slowly time-varying, and the prewhitening technique works well if the colored noise is wide-sense stationary over a period of time. However, this prerequisite may not hold in some scenarios. For instance, in wireless communication systems, the cell sites equipped with antenna arrays always receive signals from users, so the noise covariance matrix cannot be obtained separately; the reverberation noise in sonar has a fast time-variance due to the underwater dynamic media, so the noise covariance matrix at the next moment varies dramatically from that at the previous moment. As a result, the prewhitening technique is invalid under these adverse circumstances while fourth-order cumulants can handle such challenges.

Simulation Results and Discussion
In this section, numerous simulations are offered to assess the performance of the proposed non-circular FOC method that exploits the non-circular structure via actual steering vector reparameterization. Specifically, we compare our solution with its circular counterparts, Liu's [21] and Liao's approaches [24] that are based on the FOC. A ULA with half-wavelength spacing between adjacent elements is considered. Similar to the settings in [39], it is assumed that the BPSK modulated incident signals are statistically independent and have identical power. The noise is assumed to be spatially-colored Gaussian and the (m, n)-th entry of the covariance matrix is given by R(m, n) = σ 2 n 0.85 |m−n| e j π|m−n| 16 [66,67]. The signal-to-noise ratio (SNR) is defined as 10 log 10 (σ 2 s /σ 2 n ). The accuracy of the DOA estimate is measured from 1000 Monte Carlo runs in terms of the root mean square error (RMSE) which is defined as is the estimate of θ i for the n-th trial, and N is the number of signals. We first assumed that four independent sources from [−23 • , −6 • , 14 • , 36 • ] impinging on a ten-element array with mutual coupling where P = 3 with c 1 = −0.1545 + 0.4755j, c 2 = 0.1618 − 0.1176j. From Figure 1, it can be seen that with the increase of the two variables, the SNR and the number of snapshots, the RMSE of the DOA estimates declines slowly for all three methods and then stabilizes for certain SNR or snapshot values due to the manifold error resulting from the unknown mutual coupling, whereas Liao's method has a relatively large and constant error at approximately 1.57 • . It is evident that our solution is superior to Liu's method and Liao's method, particularly at low SNRs as well as small snapshot sizes, since a larger array aperture has been utilized for in our algorithm. Between the other two approaches, one can tell that in this scenario, Liu's method generally outperforms Liao's method and achieves similar performance to our method in moderate conditions (i.e., SNR≥ −7dB or the number of snapshots is larger than 3500) at a cost of computation. This demonstrates that when the ULA has sufficient sensors and the degree of coupling is not high (i.e., P is small), Liu's method performs well. In contrast, Liao's method has a clear advantage over Liu's for SNR lower than −10dB or for less than 2000 snapshots.
In the second scenario, we further appraise the performance of the algorithm in this submission under stronger mutual coupling. More precisely, we consider the case of P = 4 with c 1 = −0.1545 + 0.4755j, c 2 = 0.1618 − 0.1176j, c 3 = 0.0211 + 0.0651j. The corresponding RMSEs versus SNR and snapshot size are illustrated in Figure 2. It can be seen that the proposed approach performs the best over the entire range of SNR values and the moderate numbers of snapshots. Nevertheless, a performance degradation appears for less than 2000 snapshots mainly because Equation (12) is guaranteed to hold under strong coupling effects only for sufficiently large samples. Contrary to the performance in the first scenario, Liu's method is strictly inferior to Liao's method for all SNRs. As analyzed in Section 4.1, the size of the middle sub-array adopted by Liu's method is M − 2P + 2 = 4, and the size of the noise subspace is 7 × 3, which is much smaller than the 200 × 196 noise subspace of the proposed method and the 24 × 4 signal subspace of Liao's method. Due to the fact that a consistent estimate of a small subspace requires a sufficiently large number of snapshots, Liu's method has relatively large errors for a moderate number of snapshots, even at high SNRs. Next, we change the number of sensors to M = 8. The other settings are kept the same as in the second scenario. As discussed in Section 4.1, Liu's approach fails to work in this example as the number of sensors in the middle sub-array becomes M − 2P + 2 = 2, and up to two DOAs can be identified, but actually N = 4 signals should be estimated herein. However, the developed algorithm in this submission and Liao's method can still work under this more challenging scenario. The resultant RMSEs versus SNR and the number of snapshots are depicted in Figure 3. One can observe that the proposed solution significantly outperforms Liao's method. This can be mainly attributed to the fact that Liao's method can only utilize the signal subspace with a size of 8 × 4, whereas the proposed method exploits the noise subspace with a larger size of 128 × 124 and achieves a lower RMSE for the DOA estimates. However, as in the second scenario, the performance of our method deteriorates quickly and becomes inferior to Liao's method when the number of snapshots is low. It has been shown in the literature that ESPRIT-like algorithms have clear advantages over MUSIC-like algorithms in terms of estimation errors and robustness to array manifold perturbations in the limited snapshots situation [68]. In our simulations, it was found that at least 2000 snapshots is required to sufficiently characterize the properties of FOC, especially when mutual coupling is relatively strong. It is thus reasonable that Liao's method (ESPRIT-like) offers better performance than the proposed method (MUSIC-like) when the number of snapshots is lower than 2000.  The number of times of FBSS is q = 2, and the SNR and number of snapshots are fixed at 15 dB and 5000. Figure 4 illustrates the DOA estimates from 50 independent trials by the developed method in such a harsh case. The solid lines mark the true DOAs. It can be seen that each DOA is correctly determined, and our solution performs satisfactorily.  In the last experiment, the computational complexity of the three algorithms versus the number of snapshots is examined. The searching step-size in the MUSIC-like methods is set as 0.1 • , and the number of snapshots varies from 1000 to 10, 000. Suppose that the number of sensors in the ULA M = 10, the coupling length P = 3, and the number of independent signals N = 2. It can be observed from Figure 5 that our method costs the highest computation, followed by Liao's method, while Liu's method is most computationally efficient, which supports the analysis in Section 4.2.

Conclusions
This paper has addressed the problem of DOA estimation of non-circular signals with ULAs under the coexistence of unknown mutual coupling and colored noise. A new MUSIC-like approach in the context of FOC is developed to offer superior performance over the existing techniques. The proposed algorithm makes efficient use of the extended array observations and enables us to deal with DOA estimation without a priori knowledge of mutual coupling between sensors. In addition, a solution to coherent signals with mutual coupling between sensors is introduced, by the middle sub-array to suppress the electromagnetic nuisance and blind separation of coherent signals to enhance DOFs. Compared with the previous work, our solution is able to handle more signals and achieve more accurate DOA estimates even with relatively stronger mutual coupling. A series of simulation results verify the validity and efficiency of the proposed method.
Author Contributions: All authors contributed extensively to the study presented in this manuscript. B.W. designed the main idea, methods, and experiments, interpreted the results, and wrote the paper. J.Z. carried out the experiments, edited the manuscript, and provided many valuable suggestions to this study. All authors have read and agreed to the published version of the manuscript.