Abstract
This paper proposes a new signal decomposition method that aims to decompose a multicomponent signal into monocomponent signal. The main procedure is to extract the components with frequencies higher than a given bisecting frequency by three steps: (1) the generalized demodulation is used to project the components with lower frequencies onto negative frequency domain, (2) the Hilbert transform is performed to eliminate the negative frequency components, and (3) the inverse generalized demodulation is used to obtain the signal which contains components with higher frequencies only. By running the procedure recursively, all monocomponent signals can be extracted efficiently. A comprehensive derivation of the decomposition method is provided. The validity of the proposed method has been demonstrated by extensive numerical analysis. The proposed method is also applied to decompose the dynamic strain signal of a cable-stayed bridge and the echolocation signal of a bat.
1. Introduction
Vibration and sound signals contain intrinsic information of dynamic systems. The famous Fourier analysis can be used to project the signal into the frequency domain and identify natural frequencies of linear time-invariant systems. However, the Fourier analysis fails to study time-varying or nonlinear systems due to the nonstationarity of the signals. Therefore, numerous time-frequency analysis methods have been proposed seeking to solve this problem. The time-frequency analysis methods can be roughly classified into two categories: energy distribution and signal decomposition.
As one of the most representative methods of the energy distribution category, the wavelet transform (WT) is essentially an adjustable window Fourier spectral analysis method. With the aid of the WT, Ruzzene et al. identified natural frequencies and damping with real world data from a bridge, and Wang et al. identified instantaneous frequency (IF) of time-varying structures [1, 2]. Although the WT method has many successful engineering applications, it is difficult to achieve high resolutions in time and frequency domains simultaneously due to the Heisenberg-Gabor uncertainty principle [3]. Notwithstanding, the WT is a powerful tool for nonstationary signals in the time-frequency domain and has motivated many analogous time-frequency energy distributions such as the transform, the chirplet transform, and the synchrosqueezed wavelet transforms [4–6]. The synchrosqueezed wavelet transforms developed by Daubechies et al. are a new time-frequency analysis tool with a special reassignment method [6]. It can offer better time-frequency resolution than many other methods, and its successful applications in dynamic signal reconstruction and gearbox fault diagnosis and so forth can be found in [7–9]. However, versatile as these energy distribution category methods are, the main problem is their nonadaptive nature, since these methods utilize a family of preselected oscillatory bases to represent signals. In spite of this, the WT and other methods of energy distribution category are still important for nonstationary signal processing. Therefore, we will use the WT method in this paper to preprocess the signal for subsequent decomposition.
Empirical mode decomposition (EMD) proposed by Huang et al. in 1998 has become a representative signal decomposition method [10]. The EMD can decompose a multicomponent signal into intrinsic mode functions whose amplitude and IF can be demodulated by Hilbert transform. Due to its adaptivity, the EMD has received increasing attentions in the field of signal processing and been applied in a broad domain such as vibration signal analysis, acoustic signal analysis, and geophysical studies [11–13]. Similar to EMD, the local mean decomposition (LMD) proposed by Smith decomposes signals into a set of functions, each of which is the product of an amplitude and a pure frequency modulation signal. The LMD method has been used for electroencephalogram (EEG) analyzing [14]. However, as semiempirical methods, EMD and LMD are heuristic in nature and lack solid mathematical foundation [8]. Huang and Wu also pointed out that the Hilbert transform of the intrinsic mode functions may contain error if Bedrosian’s theorem on Hilbert transform of product functions is not established [13, 15].
Feldman introduced a very simple signal decomposition method called Hilbert vibration decomposition (HVD), which decomposes an initial signal into a sum of components with slow varying instantaneous amplitudes and frequencies [16, 17]. Gianfelici et al. introduced an iterated Hilbert transform (IHT) method to obtain slow varying amplitude and its corresponding oscillatory signal by filtering and implement the method iteratively to the residue [18]. Qin et al. have successfully utilized the IHT method for mechanical fault diagnosis [19]. The idea to decompose a multicomponent signal into monotones is very useful and deserves further study.
More recently, Chen and Wang developed a new signal decomposition method named analytical mode decomposition (AMD) [20]. The AMD method is an efficient and accurate method that separates a signal into two parts below and above the bisecting frequency [21]. Wang et al. successfully applied the AMD method to many structural vibration signal decomposition cases for modal parameters identification [7, 22, 23]. However, an error that cannot be neglected arises when the AMD method is applied for processing discrete signals [24]. The reason of the error is that the AMD method involves multiplication of signal and makes the frequencies of some components of the signal exceed the Nyquist frequency [25]. An improved multistep AMD, or an interpolation of the discrete signal, can be adopted to reduce the error [24], but the computation cost is increased significantly.
In this study, we introduce a generalized demodulation and Hilbert transform (GDHT) based signal decomposition method, which possesses the capacity of AMD but avoids computational error. The generalized demodulation is first developed by Olhede and Walden aiming to track the time-dependent frequency content of each component in a multicomponent signal [26, 27]. Using generalized demodulation, monocomponent signals with curved IF profile can be converted to another analytical signal with a constant frequency, which is very useful for enhancing time-frequency representation [8, 9]. With this in mind, components with lower frequencies are projected onto negative frequency domain so that they can be eliminated by Hilbert transform. And an inverse generalized demodulation is conducted to restore the components with higher frequencies. This procedure operates like a high-pass signal filter and can be used to extract recursively all monocomponent signals in a multicomponent signal. In the next section, the generalized demodulation theory is introduced. In Section 3, a comprehensive derivation of the decomposition method is provided. Finally, the proposed method is validated by numerical analysis and applied to practical cases such as vibration signal filtering and echolocation signal decomposition.
2. Generalized Demodulation
Consider a monocomponent signal expressed aswhere and are the amplitude and the IF of , respectively. Define the quadrature signal of asWith this definition, a complex signal can be formed asThe generalized demodulation of the signal is achieved by multiplying it with a mapping function , which givesIf a proper phase makes the signal become a component with constant frequency , that is, , the IF of the original signal can be obtained byConversely, the inverse generalized demodulation recovers the original signal by multiplying the signal with the conjugate of the mapping function; that is, , which restores the original signalThe above six equations are exactly rigorous formulas so far. In practice, however, since the phase of the signal is unknown, the Hilbert transform is always used to obtain a substitution for the complex signal . The complex signal defined by Hilbert transform is given bywhere represents the Hilbert transform of the signal .
It should be noted that substituting by implies that the Bedrosian identity is established and is an analytic signal [15], so that the signal satisfiesThis condition can be well satisfied in signals where the amplitudes and the instantaneous frequencies (IFs) are slow varying functions. Otherwise, only approximated results will be obtained if the signals contain abruptly changes caused by sudden events (such as a brittle fracture of a structural component).
3. Signal Decomposition Method
In the following content, the multicomponent signal is investigated, which is defined bywhere and are the amplitude and the IF of the th component , respectively. In many practical applications, the amplitude and the IF of the signal components are always slow varying functions. The multicomponent signal is said to be well separated if the Fourier transform of each amplitude can be neglected for and the IFs satisfy This relationship of the th IF and the th IF is illustrated in Figure 1. Thus, the phase and bisecting frequency of the mapping function can be chosen asGiven the bisecting frequency the signal can be decomposed into two parts by 3 steps.
Step 1 (to project the components with lower frequencies onto negative frequency domain). According to the generalized demodulation theory, the original signal is first processed by the Hilbert transform to obtain the corresponding analytic signal; that is,It should be noted again that (12) implies that the monocomponents of satisfy the conditions of (8). Multiplying the complex signal by the mapping function with the phase , , we obtainwhereConsidering that for , the Fourier transform of varnishes for ; and considering for , the Fourier transform of varnishes for . Note that equalities similar to (8) are implied here; that is,
Step 2 (to eliminate the negative frequency components). In order to eliminate the slow varying term , a further Hilbert transform can be conducted to . Define an operator by is an altered version of the Hilbert transform that directly produces the analytical signal corresponding to the signal . It should be noted that the Hilbert transform of a complex signal, such as , contains two subtasks that simultaneously transform the real part and imaginary part of the signal. This operator doubles the spectral components with positive frequencies and eliminates the components with negative frequencies; that is,
Step 3 (inverse generalized demodulation). Finally, an inverse generalized demodulation is performed to restore the fast varying part of the signal ,Thus the GDHT method operates like an adaptive high-pass filter. The block diagram of the decomposition method is shown in Figure 2. With the above derivation, we can conclude brief formulas of the proposed GDHT method; that is,where
Furthermore, by taking as an updated signal to be decomposed and selecting a new mapping function with phase given by (11a), the th monocomponent of the original signal can be extracted by the proposed method; that is, . Similarly, with and , the th monocomponent can be extracted. By this way, the GDHT method can be used to extract recursively all monocomponent signals in a multicomponent signal. In the following sections, we are going to test the proposed method with numerical examples.
4. Performance Analysis
In this section, the proposed GDHT method is used to process synthetic multicomponent signals. The performance of the proposed method is compared with the AMD method developed by Chen and Wang [20, 21]. Signal decomposition with constant bisecting frequency and signal decomposition with time-varying bisecting frequency are discussed in Sections 4.1 and 4.2, respectively.
4.1. Signal Decomposition with Constant Bisecting Frequency
To investigate the frequency response characteristic of the GDHT method, a zero-mean white noise signal is decomposed with a constant bisecting frequency. The variance of the white noise is set to be . The sampling frequency = 20 Hz and total sampling points are used in the simulation.
A bisecting frequency = 1 Hz () is first chosen to decompose the white noise signal. Note that the GDHT method and the AMD method both decompose the original signal into two parts; that is, . Only the slow varying part is investigated here and the result for fast varying part can be obtained by a simple subtraction. It is expected that the slow varying part of the result contains components with frequencies below 1 Hz. The one-side Fourier amplitude spectrums of the original white noise signal and two decomposed results are depicted in Figure 3(a). The result given by AMD method contains high frequency error with frequency 9~10 Hz, and the result given by the proposed GDHT method performs as expected. The frequency response of the AMD and the GDHT method is shown in Figure 3(b), which illustrates that the GDHT method is a perfect signal decomposition method, but the AMD method retains and makes negative the high frequency error.
(a) Fourier amplitude spectrum
(b) Frequency response
The second simulation is conducted with higher bisecting frequency () for extracting components with frequencies below 6 Hz. Again, the one-side Fourier amplitude spectrums of the noise and results are depicted in Figure 4(a), and the frequency response of the AMD and the GDHT method with = 6 Hz is shown in Figure 4(b). The result given by AMD method contains high frequency error with frequency 6~10 Hz and eliminates the components with frequencies 4~6 Hz. The result given by the proposed GDHT method performs as expected, which illustrates that the GDHT method is valid as well.
(a) Fourier amplitude spectrum
(b) Frequency response
4.2. Signal Decomposition with Time-Varying Bisecting Frequency
The GDHT method can be used to decompose nonstationary signals with time-varying frequencies. To investigate the performance of the GDHT method, a signal with two frequency modulated components is considered:where , . Thus the IFs of the two components are and Hz. The sampling frequency = 20 Hz and a total sampling time = 30 s are used in the simulation. This signal is very similar to the “warblet,” which has been found to be very useful in analyzing actual radar data [28]. The radar signal returning from small ice fragments rises and falls in frequency in a periodic manner.
The purpose here is to retrieve these two components with overlapping frequencies. Firstly, the Fourier amplitude spectrum of the signal is shown in Figure 5(a), which gives no clue to select a bisecting frequency. This provides evidence that the Fourier transform is not suitable for nonstationary signal processing. Thus a continuous wavelet transform is performed to trace the time-frequency energy distribution of the signal, in which the complex Morlet wavelet is used [2]. The WT scalogram of the signal is shown in Figure 5(b), from which the fluctuations of the instantaneous frequency of the signal can be observed. The energy distribution in the scalogram coincides well with the IFs and . Although the WT scalogram cannot provide an unambiguous bisecting frequency for the decomposition method, a mapping function can be selected considering the variation tendency of the IFs.
(a) Fourier amplitude spectrum
(b) The WT scalogram of
To make the signal separable in its Fourier spectrum, a mapping function with phase function is adopted, which corresponds to the mapping frequency Hz. According to (4), the generalized demodulation of the signal is achieved by multiplying the mapping function with the analytical form of the original signalwhere the operator is defined by (16). Therefore, the IFs of the components are mapped into and Hz, respectively. The Fourier amplitude spectrum and the WT scalogram of the mapped signal are shown in Figures 6(a) and 6(b), respectively. Obviously, the two components of the mapped signal can be distinguished from each other by the Fourier spectrum or by the wavelet scalogram. There is a trough at frequency 1.55 Hz in the Fourier amplitude spectrum, which suggests that a proper bisecting frequency can be chosen as Hz. With this bisecting frequency, the signal can be decomposed into two parts and by the GDHT method.
(a) Fourier amplitude spectrum
(b) The WT scalogram of
As shown in Figure 7, the decomposed components and from the GDHT method are in excellent agreement with the exact components and , respectively. The IFs of the decomposed component are computed with Hilbert transform [29], and the results are compared with the exact IFs, as shown in Figure 8. The IFs of the decomposed components are very close to the exact ones except for the errors in two ends of the signal. The error is caused by the end effect of Hilbert transform and can be reduced by a simple mirror image technique [20]. Anyhow, with the decomposed components from the GDHT method, the IFs can be identified accurately at most of the time. Therefore the GDHT has application value in practice because the variation of frequencies of the signals always contains intrinsic information about the dynamic systems.
(a) Slow varying component
(b) Fast varying component
To further compare the GDHT with AMD method for the time-varying bisecting frequency, the WT is applied to analyze the decomposed components. The WT scalogram of the slow varying part and the fast varying part decomposed by the GDHT method is plotted in Figure 9. Although the time-frequency resolution of the WT is limited by Heisenberg uncertainty principle, it is obvious that the energy of the decomposed slow varying signal mainly distributes in the region below the bisecting frequency and, inversely, the energy of the decomposed fast varying signal mainly distributes in the region above the bisecting frequency . Two simple schematic diagrams are given in Figure 10 to illustrate the characteristics of the GDHT method. Figure 10 shows that the slow varying part decomposed by the GDHT method contains no signal component with frequency higher than the bisecting frequency, while the fast varying part contains no signal component with frequency lower than the bisecting frequency. This shows that the GDHT method is a perfect and adaptive filter for discrete signal.
(a) Slow varying part
(b) Fast varying part
(a) Slow varying part
(b) Fast varying part
As a comparison, the WTs of the decomposed components by the AMD method are also conducted and the wavelet scalograms are plotted in Figure 11. Obvious deviations from Figure 9 can be observed in the scalograms of Figure 11, which is caused by the discretization of the signal. The computed slow varying signal by the AMD method contains components with frequencies higher than the bisecting frequency, as shown in Figure 11(a). And the computed fast varying signal contains components with frequencies lower than the bisecting frequency, as shown in Figure 11(b).
(a) Slow varying part
(b) Fast varying part
The effect of the discretization for the AMD method with time-varying bisecting frequency is similar to the time-invariant scene given in Section 4.1. To illustrate this effect, two simple schematic diagrams are given in Figure 12 to explain the deviations observed in the wavelet scalograms. As shown in Figure 12(a), when , the decomposed slow varying signal retains and makes negative the signal component with frequency higher than ; and when , the decomposed slow varying signal retains and makes negative the signal component with frequency higher than and wrongly eliminates the signal component with frequency . The performance of AMD method to decompose the fast varying signal can be obtained by a simple subtraction, as shown in Figure 12(b). The results show that a sampling frequency 4 times higher than the bandwidth, or maximum component frequency, should be adopted for correct decomposition of the signal by the AMD algorithm, which doubles the AMD algorithm computational cost.
(a) Slow varying part
(b) Fast varying part
5. Case Study
5.1. Dynamic Strain Signal Filtering
The proposed GDHT signal decomposition is used to process the dynamic strain signal of the Tai-ping Lake Bridge. This bridge is a prestressed concrete cable-stayed bridge with a total span of 380 meters. The strain gauges are installed at the upper surface of the bottom plate of the box girder, and the sampling frequency is set to be 50 Hz. A typical dynamic strain signal for a period of 24 hours is selected and shown in Figure 13(a), which contains the slow varying components caused by the variation of environment temperature and the fast varying components caused by the vehicle load. The signal is decomposed into two parts by the GDHT method with a bisecting frequency of 0.001 Hz. The results are shown in Figures 13(b) and 13(c). The decomposed slow varying component contains no high frequency error and the fast varying component is free of slow varying excursion. The fast varying components are very useful for vehicle load statistics and fatigue analysis of the structure.
(a) Dynamic strain signal
(b) Slow varying component
(c) Fast varying component
The total sampling number of the dynamic stain signal is 4.32 × 106 and the computing time of GDHT is 3.75 sec (by a computer with 3.1 GHz processor, 4.0 GB of RAM). Given the huge number of the discrete sampling signals, the decomposition is relatively fast and is suitable for engineering applications.
5.2. Echolocation Signal Decomposition
The echolocation signal of a bat is decomposed in this subsection. It is well known that the bats judge distances and identify the objects through the echolocation signal. A typical echolocation signal of a bat is plotted in Figure 14. This signal has been studied by Yu and Zhou [30] and the data can be downloaded at [31]. It should be noted that the duration of the signal is 0.0028 seconds and the sample interval is 7 μs according to [27]. The WT of the signal is given in Figure 15, from which a set of bisecting frequencies can be easily determined for the GDHT method. The time-frequency domain is divided into five parts by the four bisecting frequencies shown in Figure 15.
The five decomposed components are shown in Figure 16. It should be noted that the amplitudes of the first component and the fifth component are very small. This means that the original signal can be well reconstructed by the three components C2, C3, and C4. The Hilbert transform is employed to compute the instantaneous frequencies of these five decomposed components. The results are shown in Figure 17, which provides a better time-frequency resolution than the WT. The amplitude is grey-coded in Figure 17, where white corresponds to the smallest values and black corresponds to the largest values. This time-frequency representation method is inspired by the Hilbert spectrum method proposed by Huang et al. [10].
6. Conclusions
This paper describes a novel generalized demodulation and Hilbert transform based signal decomposition method to separate a signal into two parts above and below a bisecting frequency. The bisecting frequency can be selected as a constant or a time-varying function. The generalized demodulation is first applied to project the signal components below the bisecting frequency onto negative frequency domain, and the Hilbert transform is then used to eliminate the negative frequency components. And an inverse generalized demodulation is performed to restore the components with frequencies higher than the bisecting frequency. The characteristic of the method is analyzed by theoretical derivation and numerical examples. The proposed method is finally applied to process a typical 24-hour dynamic strain signal and the echolocation signal of a bat for validating its effectiveness and high efficiency. The proposed method yields better results than the AMD method for discrete signals and provides a better time-frequency resolution than the WT.
Conflicts of Interest
The authors declare that they have no conflicts of interest.
Acknowledgments
The work described in this paper is supported by the National Natural Science Foundation of China (Project no. 51408177) and China Postdoctoral Science Foundation (Project no. 2014M551802). The authors would like to thank Fei-Yu Wang for modifying the manuscript.