Abstract
We focus on the decomposition problem for nonstationary multicomponent signals involving Big Data. We propose the kernel sparse learning (KSL), developed for the T-F reassignment algorithm by the path penalty function, to decompose the instantaneous frequencies (IFs) ridges of the overlapped multicomponent from a time-frequency representation (TFR). The main objective of KSL is to minimize the error of the prediction process while minimizing the amount of training samples used and thus to cut the costs interrelated with the training sample collection. The IFs first extraction is decided using the framework of the intrinsic mode polynomial chirp transform (IMPCT), which obtains a brief local orthogonal TFR of signals. Then, the IFs curves of the multicomponent signal can be easily reconstructed by the T-F reassignment. After the IFs are extracted, component decomposition is performed through KSL. Finally, the performance of the method is compared when applied to several simulated micro-Doppler signals, which shows its effectiveness in various applications.
1. Introduction
In practical applications such as oceanic investigation [1], radar [2], biomedical application [3], and mechanical fault diagnosis [4], we need to represent and process nonstationary signals. While Big Data can be explicitly regarded as a good fortune, great challenges also appear with extensive datasets. In general, Big Data introduce challenges in which properties, such as time domain, space-time domain, computational complexity, and information energy, are mingled in complicated ways with data resources, as shown in Figure 1. Particularly, a multicomponent signal (MCS) is a superposition of monocomponent signals whose may overlap (or cross) in time-frequency (T-F) domains. The fundamental problem in processing signals is to extract the useful information stored in the amplitude-modulated (AM) and/or the frequency-modulated (FM) signals. The time-frequency analysis (TFA) is a central tool to capture the time-varying features of multicomponent chirp signals. To analyze further characteristics of these signals, we are more eager to decompose them into monocomponent signals.
1.1. Related Work
In the past few decades, many signal decomposition approaches have been developed to capture the accurate monocomponent signals. Commonly, the typical methods can be summarized as time-frequency distributions [5], the empirical mode decomposition [6] and its multivariate extensions [7], the local mean decomposition [8], the reassignment method [9, 10] such as the synchrosqueezing transform (SST) [9], the sparsification methods [11, 12], and so on.
However, most of the decomposition methods are signal-dependent and presume that the signal components meet strict divisible conditions in the T-F domain [9], and thus they are affected by additive noise and close or overlapped instantaneous frequencies (IFs) of signal components. Recently, the intrinsic chirp component decomposition (ICCD) method proposed in [13] and used in [14] revealed the incapable of separating crossed or overlapped chirp signals by using a joint-estimation scheme. Essentially, the ICCD based on the short-time Fourier transform (STFT) can accurately reconstruct overlapped components. Unfortunately, according to the Heisenberg uncertainty principle, the STFT, as the liner TFA, cannot achieve an arbitrarily high T-F resolution. The overlapped components are very close in the T-F plane and the mask with low T-F resolution may cover all components and the unwanted cross terms, which limits the applications of the ICCD.
To analyze overlapped signal components, TFA is an important issue which deserves to be mentioned for the extraction of feature information from the components. Currently, the extraction accuracy relies on the energy concentration of TFR generated by TFA methods. There exist a variety of types of TFA methods, and they can be normally divided into two categories: the parametric TFA (PTFA) methods and the nonparametric TFA (NPTFA) methods. PTFA methods, such as polynomial [1, 15], spline-kernelled chirplet transform (SCT) [16], and sinusoidal models [17], often involve the high-dimensional search of the IFs, which is very time consuming. Moreover, the predesigned parametric models may be only suitable for special applications. In this paper, we are more concerned about the NPTFA methods which are known for their more adaptability in practical applications. The traditional NPTFA methods by short-time Fourier transform (STFT), continuous wavelet transform (CWT), and Wigner-Ville distribution (WVD) are subjected to poor energy concentration. Recently, a class of TFA techniques, referred to as T-F reassignment method, proposes to improve the readability of TFR clearly. The synchrosqueezing transform (SST), belonging to this class, is a postprocessing technique based on the CWT to obtain better localized TFR of nonstationary signals [18, 19]. The SST provides a higher resolution that moderates the limitations of linear projection based traditional NPTFA methods, such as STFT, CWT, and SWT [19, 20]. Moreover, the SST reconstructs the energies of these transforms, so that the resulting energies of coefficients are condensed around the IF curves of the components. However, the interference of serious noise on TFA has attracted a major attention, since it brings about a dominant error source when it appears in the contaminated signal. Errors on account of the serious noise induces false maxima (maxima other than the autoterms) in the T-F plane. In [21], an ant colony optimization (ACO) algorithm and TFR by WVD is adopted to estimate IFs of the components under serious noise. This method is applicable to represent overlapped multicomponent signals. Nevertheless, we find that these algorithms mentioned above cannot further decompose the nonstationary overlapped components, which may be often seen in various signal processing of practical applications.
1.2. Proposed Method
Notice that despite all those efforts decomposing overlapped MCSs is still a challenging task. In this paper, we develop the intrinsic mode polynomial chirp transform (IMPCT) which takes into account IFs characteristics and a novel KSL method is proposed in order to obtain a robust decomposition overlapped components in high noise. For a nonstationary signal, instantaneous amplitudes (IAs) and IFs are two of the foremost properties. In [22], the IAs of the real signal are approximated by a third-order or more orders polynomial. To represent complex chirp signals, the IFs and IAs are both modeled as discrete chirp series in the IMPCT, which enables energy concentration and refines the most important coefficients in high noise. The coefficients for this model are distributed according to the time-varying band-pass filter in the T-F plane and not projected from either side. In addition, the distributed properties of the IMPCT method and carefully designed series make the representation more robust to swift IFs variations.
For the practical implementation of a decomposition method, a very interesting solution is the kernel sparse learning (KSL) approach, which has not yet been developed for the multicomponent decomposition. While we propose the KSL, as an active learning, similar to kernel dictionary learning, which has been successfully applied to classification problems [23–25], this model works differently. In fact, the purpose of the algorithm is to maximize the accuracy to minimize the amount of training samples. For such aim, smart strategies are adopted to pick out the most prominent training samples. An advantage of our approach is that it separates closer crossing components by KSL in parallel.
The remainder of this paper is organized as follows. In Section 2, the theory of IMPCT is described in detail. The IFs estimation by T-F reassignment by the path penalty function and the decomposition model by kernel sparse learning are given in Section 3. Section 4 shows simulation results of decomposing some real and micro-Doppler signals are provided, and conclusions are followed in Section 5.
2. IMPCT for TFR
2.1. The Motivation of IMPCT
First of all, a numerical example of a noisy environment is applied to demonstrate the motivation of our method proposed. A simulated signal is considered as follows: with which is sampled at a sampling frequency of 2000 Hz and the signal is polluted by a noise with SNR=-2dB. The T-F representations shown in Figures 2(a) and 2(b) are generated by SCT [16] and SWT [19], respectively.
(a) SWT
(b) SCT
(c) Iteration=0 IMPCT
(d) Iteration=50 IMPCT
The SWT, the continuous wavelet transform- (CWT-) based SST, sharpens the time-frequency representation of a multicomponent signal but suffers from a lower T-F resolution, as shown in Figure 2(a). The T-F coefficients of the IFs trajectories are obtained by SWT only in the frequency domains. Thus, it is no wonder that the SWT cannot accurately estimate time-varying IFs of an overlapped signal as well. In Figure 2(b), the SCT is able to generate a TFR with a poor energy concentration for the contaminated signal with nonlinearly time-varying IFs. If the chirp rate of SCT deviates from the true IFs, the T-F representation will smear heavily.
An intuitive solution stems from a combination of high-energy concentrations that partially result from the different modulation frequencies. The new T-F representations of the IMPCT will have better performance than the SCT, as shown in Figures\ 2(c) and 2(d) with 50 iterations.
Therefore, it is necessary to develop alternative approaches to deal with the challenging cases that there exist components overlapping in the T-F domain. We provide the following definition.
Definition 1. The two components and in model (3) are said to cross (or overlap) in the T-F domain. If there exist a time instant , a bound and a step shift up and down such that their IFs satisfy
Equation (3) illustrates a strict definition of these overlapped components for our algorithm to identify the case that some components are tangent to each other (). Figure 3 shows the strict overlapping and tangent of the two components. After shifting up and down, if there is a point such that at other intersecting time , the components are strictly overlapping. It is worth noticing that the IFs of the arbitrary two components have different curvature at the intersecting time .
(a)
(b)
A signal can be considered an intrinsic mode chirp function, with accuracy and and meeting the following conditions [14]:
Here, to overcome these limitations, a robust technology architecture for overlapped nonstationary signals is proposed, as shown in Figure 4. Initially, the input signal is passed by the -point IMPCT stage. The stage brings about the IMPCT coefficients shown as . The inverse IMPCT operation for a signal sequence of length is performed. These characterized polynomial coefficients are concatenated and rearranged to the IFs estimation for each of the input nonstationary signals. Then, to decompose overlapped components in practical applications, we need to deal with a very significant but challenging issue of how to estimate IFs of overlapped components for the IMPCT. The second problem is how to determine the criterion which can separate the signals under micro-Doppler effect and noise interference. In the following sections, the solution for these two problems will be provided about detailed theoretical analysis.
2.2. IMPCT Representation
The Fourier series analysis brings to light the sinusoidal property, so it may be utilized for the reconstruction of some signals with linear phase [14]. However, when processing nonlinear polynomial phase signals (NPPS) often involved in the practical application, a signal representations can be obtained by using the IMPCT. Namely, the IMPCT can be adopted to demodulating components of interest, reducing them to the sparse sinusoidal components. In this section, the Fourier series is stretched into the complex polynomial Fourier signal model as the IMPCT. Consider a multicomponent polynomial phase signal contaminated by a white Gaussian noise aswhere is the amplitude coefficient function and indicates the phase of the -th component; shows the noise. In general, the IMPCT of can be formulated as follows: For all points in the -dimensional frequency domain that correspond to the positions of signal components characterized by , we can acquire That is to say, at the positions of signal components: .Thus, can be deemed to a sparse -component representation. In addition, for has limited value, negligible when compared to (7).
Let us inspect the signal vector with elements , superimposed of a sum of NPPS: where the polynomial coefficients are supposed to be bounded integers and is the total number of discrete samples. To offer a sparse representation of nonlinear signals, the discrete form of the IMPCT is adopted, which may be expressed as When the set of IMPCT parameters is selected to obtain the adaptive coefficients of the NPPS , thus we can obtain the -th signal component is demodulated and the sinusoid is paramount in the IMPCT spectrum. Under this circumstances, the energy of the spectrum is extremely concentrated at . On the contrary, while for all the spectrum consisted of some dispersed phase. Under high noise environment, we suppose that the minimum spectral component having amplitude coefficient can be also primary in the spectrum, meaning that it is above the aggregate value of all other components (in the worst case above their accumulation): , where is the minimum of all , for . Therefore, a search is performed for sets of parameters so that holds, and is strongly concentrated. All in all, multicomponent set of parameters are calculated to obtain the optimal concentrated vector for all . The maximized component will lie in the frequency .
The IMPCT representation explains the motivation for obtaining an initial IFs estimate in the algorithm. The IFs estimate will be based on the value of . In the next section, the T-F reassignment algorithm by the path penalty function is proposed in order to more accurately estimate IFs with the serious noise cases.
3. Multicomponent Decomposition
3.1. IFs Estimation by the T-F Reassignment
There are several methods to implement the IFs estimation. Here, we will adopt a reassignment by the path penalty function to further reduce search space. The set of the IFs estimation is detected in the selected points with the best path penalty function gain. The algorithm is performed recursively until we are no longer able to search out a new point that can reduce values of the path penalty function.
In T-F plane, the frequencies and time sequence can be contained. The set of paths between two ending points is , which aids in searching out the optimal path impossible. Fortunately, the algorithm can be implemented recursively as an example about the generalized Viterbi algorithm [26]. The optimal paths can be obtained as which join the moment and entire points to the moment , , for . It can be defined aswhere the set is obtained in (9) and contains all paths between the moment and entire points to the moment . Simultaneously, a sum of the path penalty functions is revealed by for each ridge . In the Viterbi algorithm, a chirplet path (11) is known as the partial optimal path. Currently, the IFs estimation, within the time interval , can be expressed as At the next moment , the partial optimal paths will be represented as the cascade of (11) with the points at the new moment , , .
By selecting paths with shorter chirplets, one can search out chirplets with an increasing function of the distance to further decrease search space. A new minimum value at the next moment is represented as for each . It is worth noting that is positive constant to capture solutions of the constrained partial optimal path with different values of length. To illustrate algorithm, an example about two overlapping signals is shown in Figure 5. In the considered case, an optimal path (IFs estimation) connects points (2,1) and (1,1).
To form a best chirplet path using the path penalty functions, two key constraints in the connections stage from each point should be satisfied as where stand for the estimated polynomial chirp phase in (9), respectively. The parameters and are determined by the boundaries of the frequency modulation rates of the considered signal. In Figure 3, it can be easily seen that two best chirplet paths are synthesized by connecting the partial optimal path to points (2,1) and (1,1), as well as (1,2) and (2,2), respectively. Note that the new IFs estimation can perform an update of the previous value.
Furthermore, we only need to solve the 1D optimization problem in partial optimal path to estimate IFs ridges. Thus, our method can be efficiently used in various real applications and a numerical overlapped signal mentioned in Figure 2(a) is analyzed in the high noise. The relative result of the IFs estimation is shown in Figure 6. Each mode occupies difference boundary region of the T-F domain as shown in Figure 6(a), and the corresponding energy is concentrated around the IFs ridges in Figure 6(b).
(a)
(b)
In [27], the grid search method is adopted to solving the optimization problem. Especially assuming that the chirp phase search space is , the set of the discrete polynomial chirp phase can adjust the grid resolution. Considering the training set in the kernel sparse learning, our method searches out proper path from a dictionary which is reconstructed by possible partial optimal paths in the set (13). Each subset is generated from a different starting point of the sample, implemented independently to link the each other and used to update a different pool of optimal paths. Therefore, a parallel optimization calculation can be adopted in our learning method.
3.2. Decomposition by Kernel Sparse Learning
One can rewrite equation (5) where ; ; with where and indicate complex polynomial Fourier coefficients which should be identified. For each signal component , its inversion formula corresponding to (7) is written as where shows a column vector whose elements are the signal values and indicates the Hermitian transpose of . The signal is supposed to be time limited within . We have obtained the IMPCT of this -th component, according to (8). Through accumulating the above relationship for , the following is available:
Considering formula (9), for the signals meeting the proposed model, this relation can be rewritten as By indicating and the decomposition of the kernel matrix with the elements , we obtain where will be equal to the -th signal components.
To recover the expected components, we need to deal with nonlinear equations in (15). Since this inverse problem is commonly ill-posed, the kernel sparse learning is employed to process real data. In traditional kernel dictionary learning the dictionary is expressed as a linear combination of a nonlinear representation of the data in [23]. This is defined aswhere is defined as a series of coefficients and is nonlinear function. is the sparse weight matrix combining fixed basis in order to engender the dictionary atoms.
In kernel sparse learning, the fixed basis is replaced by a nonlinear combination of the input. In (22), a kernel could only be used while solving the inverse problem. This allows expressing (22) as where the kernelized data matrix is formulated instead of the data matrix. One can give the kernel upfront , referred to as . We design the learning as
For kernel sparse learning, one can notice that the gradients for different terms and updating the coefficients in (25) are easy to calculate. This is given by
When the amount of samples is much bigger than the dimension of the data, storing and processing the kernel matrix are great difficulty in terms of memory. Consequently, a more efficient scheme needs to be implemented. Using the eigenvalue decomposition, the kernel matrix can be denoted as where shows the eigenvector matrix and indicates the diagonal matrix of eigenvalues sorted in descending order. Alternating minimization of (25) brings about
Solving (28) is straightforwardly updated along with the update of IMPCT, identically; plays the role of data. In this study, we denote how (29) can be efficiently addressed by adopting the separation variable method.
An agent variable is introduced, which can ideally ensure that the agents and variables are equivalent in each iteration. In reality, the equality is only enforced at convergence. Therefore, under the slight relaxation of constraints, the augmented Lagrangian method [28] can be applied to improve scalability and availability of KSL. One can rewrite (29) as In the above description, the equivalence between the agents and the variables can be guaranteed by the hyperparameter . By regulating the level of the value, the strength of the constraint will be controlled to enforce the equivalence at convergence. Initially, the algorithm sets a low value of and gradually enhances the value of the hyperparameter to efficiently deal with the optimization problem.
In the high noise environment, one may choose more complex methods to make sure the number of components is accurate, such as by the split Bregman approach [29, 30]. To further relax the constraint, a Bregman variable is used in the following mode: A relatively large is conducive to make up for IFs estimation errors. One starts a high value of the Bregman variable but it is automatically updated in each iteration. Therefore, an effective approach is to automatically adapt the equality between the agents and the variables to the noise level such as the methods in [29].
In order to achieve parallel computing, through adopting the alternating directions method of multipliers, (31) can be decomposed into the following two easier subproblems.
It can be easily seen that the first subproblem (32) is an easy least squares problem forming a convergence in a mode of the pseudo inverse problem. In the subproblem (33), a hard threshold needs to be taken to enforce the equality between the agents and the variables at convergence. The model is efficient to deal with two subproblems in parallel. In addition, the KSL has the ability of select optimal training samples in order to minimize the error of the decomposition process while minimizing the number of learning samples to label and thus to cut down the costs interrelated with the training sample collection.
To test the result of the training samples, we need to generate the features for the sample signal . Simultaneously, we obtain a learned dictionary as shown in Figure 7, where signals are merged only to provide a more valuable time-frequency representation from the IMPCT. It is worth noting that, in all the following simulation experiments, the initial training samples can be obtained by the IMPCT based on different feature modalities, such as noise, micro-Doppler effect, and scattering. Feature modalities quantity and the number of the generated samples corresponding to each feature modality are determined according to the complexity of the signal to be decomposed.
The question is how to determine the dictionary size hyperparameter . Note that although the same initial value is adopted for all proposed numerical results, a more appropriate choice of this parameter may bring about better estimation results. The choice is an empirical balance between the MSE acquired on the training features, the sparsity of the kernel dictionary, and the interpret-ability of the final dictionary (decision criteria that rely on the specific application can also be considered).
Having estimated the IFs of the signal () in the previous section, the decomposition of signal with severely overlapped components is executed by KSL as shown in Figure 8. It indicates that all the components are successfully separated and the SNR is obviously improved. The detected IFs ridges from the TFR of the signal are generated in Figure 8(a). Figures 8(b) and 8(c) denote the corresponding individual components. The improvement of the superimposed of all components () is shown in Figure 8(d).
(a)
(b)
(c)
(d)
4. Validation
It is known that micro-Doppler effect has been diffusely applied to target detection and recognition in sonar systems [31, 32]. In general, such signals from scattering points of the transmitted signal (or from the time delay and Doppler shift of the received echo) may seriously overlap in the T-F domain. To accurately characterize the TFR of micro-Doppler signals, it is key to decompose overlapped signal components, first and foremost.
In this section, the presented approach will be used to decompose both simulated micro-Doppler signals to verify the performance of the proposed method. The maximum number of iterations is set to 80 and we adopt a more suitable . To quantify the noise level of the inspected signal or the accuracy of IFs estimation and signal reconstruction, the SNR (unit: dB) of the signal is defined as where indicates original signal or true IFs; is reconstructed signal or estimation IFs.
In all simulated examples, a zero-mean Gaussian white noise with variance and with i.i.d. real and imaginary parts is artificially added to the resulting micro-Doppler signal to analyze the effect of complicated AWGN.
Example 2. Consider an FM signalinduced by some typical moving targets [2] (such as rotating, tumbling, and coning targets); the signal is contaminated with noise with . Herein, the sampling frequency is 100 Hz and the time duration of the signal is set to 15 s. The IMPCT of the signal, as the starting point of the presented algorithm, is displayed in Figure 9(a).
The experimental results are observed after every 10 iterations, until the 80th iteration, are shown in Figures 9(b)–9(i). It can be easily seen that after the 10th iteration a large amount of elements has concentrated the spectrum energy at many T-F points. However, due to the nonlinear polynomial phase effect in (6), the spectrum energy has already started to concentrate at the autoterm, as shown in Figure 9(b). With the iteration, the stronger the spectral energy is enhanced, in contrast, the weaker the noise interference is. The optimal TFR with a superior energy concentration from the 80th iteration is shown in Figure 9(i) due to the energy positioning update mechanism (7) and a sparse representation (8) avoiding the noise interference, used in (9).
The detected IFs ridges from the TFR of the signal are revealed in Figures 10(a)–10(c). It reveals the results of IFs ridges with 0th iteration, 10th iteration, and 80th iteration, respectively. Before the 10th iteration, the T-F reassignment detector by the path penalty function cannot track precise IFs for the micro-Doppler signal with highly overlapped components. On the other hand, the proposed T-F reassignment can effectively distinguish these crossed components by considering the partial optimal path of the ridges. To abate the noise interference, one can further smooth the estimated IFs by nonlinear spline adaptive fitting [33] with a polynomial Fourier model, as illustrated in Figure 10(c).
Having estimated the IFs ridges, the separation of the multicomponent signal is executed by the KSL, as shown in Figures 10(d) and 10(e). It displays that all the components are successfully discomposed and the SNR is availably improved.
(a)
(b)
(c)
(d)
(e)
(f)
(g)
(h)
(i)
(a)
(b)
(c)
(d)
(e)
Example 3. Here, we consider again a poly-harmonic FM signal, expressed as with where a sampling frequency is set to 100 Hz and the its time duration is 15 s. The signal is contaminated by high white noise and the SNR is 0 dB. The TFRs acquired by the presented IMPCT are compared with the results acquired by SWT and SCT, as shown in Figure 11. In addition, IFs ridges of the noise signal are estimated by the proposed T-F reassignment (see Figures 11(c), 11(f), and 11(i)).
The TFRs of nonnoisy and noise signal obtained by SWT, shown in Figures 11(a) and 11(b), respectively, can barely show the inherent T-F pattern of the signal. Specifically, the TFRs only reveals the clear IFs trajectory for the autoterms of signal due to an easy approximation by a linear function. However, the TFRs of the cross terms are too blur to accurately estimate the IFs ridges of the overlapping components. In Figures 11(d) and 11(e), it can be noticed that the TFRs obtained by SCT scatter the energy around the IFs due to their rough frequency resolution induced by the high noise, and it is inaccurate to extract the IFs. In Figures 11(g) and 11(h), it is revealed that, due to the proposed optimal concentrated vector, the IMPCT surpasses SWT and SCT as it clearly shows the true T-F pattern of the signal. The IMPCT generates the TFRs with an outstanding concentration, based on which the IFs can be extracted accurately (see Figure 11(i)). However, in Figures 11(c) and 11(f), it can be seen that partial IFs ridges of cross terms corresponding to some points are truncated in disorder. The simulation results display that both methods suffer from high noise and cross terms. Obviously, one cannot accurately decompose the overlapped components by these two methods.
For all the following tests, the MSE versus variable SNR will be calculated (from -10 to 10 dB below, with a 1 dB step). The MSE is expressed as where indicates the estimated IFs and the true IFs index.
We consider the IFs estimation from the TFRs of IMPCT, SWT, and SCT, using the T-F reassignment based on the path penalty function. The obtained MSEs are displayed in Figure 12. Signals with fast IFs variations have been chosen to determine good performance of the presented method in such scenarios. Because of the distributed property of the optimization, the IMPCT is robust to characterize the TFRs, i.e., to the nonlinearity IFs ridges. Analyzing the results in Figure 12, we demonstrate that, through slightly promoting the IFs variations or adding the noise, the better improvement level of IMPCT over SWT and SCT is obtained.
(a)
(b)
(c)
(d)
(e)
(f)
(g)
(h)
(i)
5. Conclusion
In this paper, we presented a novel approach which combines the KSL algorithm with the T-F reassignment to decompose seriously overlapped components. The IMPCT is first proposed for the TFRs of a multicomponent nonstationary signal in high noise and can adapt to some practical applications. Comparing with traditional the decomposition approaches, the KSL based on the T-F reassignment can decompose signals with overlapped components. The modified T-F reassignment algorithm by the best path penalty function estimates the main IFs ridges of all the components and then reconstructs the obtained IFs ridges according to the partial optimal path. To estimate performance of the presented method, we analyze the numerical experiments including a real-world signal, which clearly illustrates that it outperforms the IMPCT based on the T-F reassignment and KSL in high noise environments. In addition, it outperforms the state-of-the-art TFRs methods when polynomial phase signal is contaminated by a white Gaussian noise. It should be emphasized that, in essence, the KSL algorithm involves postprocessing of ridges to optimize the IFs estimation of components contaminated by noise or induced by moving targets. Further, the applications in practical scenarios will be the part of our research development.
Data Availability
The data used to support the findings of this study are available from the corresponding author upon request.
Conflicts of Interest
The authors declare that they have no conflicts of interest.
Acknowledgment
This work was supported by the National Natural Science Foundation of China (61471309 and 61671394) and the Fundamental Research Funds for the Central Universities (20720170044).