Nonlinear Time-Varying Spectral Analysis : HHT and MODWPT

Time-frequency distribution has received a growing utilization for analysis and interpretation of nonlinear and nonstationary processes in a variety of fields. Among them, two methods, such as, the empirical mode decomposition EMD with Hilbert transform HT which is termed as the Hilbert-Huang Transform HHT and the Hilbert spectrum based on maximal overlap discrete wavelet package transform MODWPT , are fairly noteworthy. Comparisons of HHT and MODWPT in analyzing several typical nonlinear systems and examinations of the effectiveness using these two methods are illustrated. This study demonstrates that HHT can provide comparatively more accurate identifications of nonlinear systems than MODWPT.


Introduction
Time-frequency TF analysis has experienced a number of qualitative and quantitative changes during the last three decades, and has gradually received growing attentions and further applications in a variety of fields such as radar, water waves 1 , fault diagnose, geophysics, and biological signals.However, most traditional signal processing methodologies, developed under rigorous mathematical rigor, are based on linear and stationary assumptions.As the data from the real world are generally neither linear nor stationary, the traditional data analysis methods aimed at linear and stationary signals and processes are becoming glaringly inadequate.In recent years, several new methods have been introduced to analyze nonlinear and nonstationary data.For instance, spectrogram and Wigner-Ville distribution 2 were designed for linear but nonstationary data.In addition, to accommodate for nonlinear but stationary and deterministic processes, Tong 1990 , Kantz, and Schreiber 1997 , and Diks 1999 raised various time-series-analysis methods 3-5 .

Mathematical Problems in Engineering
Whereas most signals and processes, either natural or artificial ones, are most probably be simultaneously nonlinear and nonstationary.Thus this makes finding a suitable approach to such kind of data a dread in a way.
Huang and Shen put forward that all the analyses in terms of a priori established basis have drained all the physics out of the analyzed results, because any a priori basis could not possibly fit all the variety of data from different driving mechanisms 6 .As the ultimate goal for data analysis is not only to find the mathematical properties of data, but also to excavate the physical insights and implications hidden in the data, the adaptivity becomes absolutely necessary for nonlinear and nonstationary data.Following with the two historical views of nonlinear mechanics of Fourier and of Poincaré, Huang et al. proposed the Hilbert view based on a new method, called empirical mode decomposition EMD and Hilbert spectral analysis, which is termed as the Hilbert-Huang Transform HHT 7 .The approach of combining EMD and HT differs from the Fourier transform and wavelet analysis.It provides a faithful representation for the nonlinear and nonstationary data analysis.The EMD is conceptualized as an alternative means to separate a multicomponent signal into its monocomponent constituents through a progressive sifting process to yield an empirical base consisting of intrinsic mode function IMF components.To ensure that these IMFs have the well-behaved Hilbert transform and conform to a narrowband condition, more definitive stopping criteria such as max sifting iteration are defined.In this way, the data are expanded in a basis derived from the data itself.
Recently, the EMD algorithm, acting as a manner with highly data-driven characteristic of data decomposition, plays a role of either nonlinear or nonstationary processes.With these nice properties, the EMD has been used to calculate the Hurst index of longrange dependence processes LRD and network traffic 8 , additionally, being an approach of synthesizing fractional processes 9 .In a number of studies, the combination of EMD and HT has been applied and advocated in a variety of problems covering geophysical 10 , biomedical engineering 11 , and also fluid mechanics as nonlinear water waves and turbulence data 12 .Although several problems still exist in the EMD and the associated Hilbert transform, such that the EMD method may produce mode mixing for some signals 13, 14 , the EMD HT owes strength of being data dependent and provides a potentially viable method and validity of both nonlinear and nonstationary data analyses.Furthermore, it offers a new sight for nonlinear and nonstationary signal processing, that is, individual component signal with physically meaningful instantaneous frequency can be obtained by appropriate signal decomposition method.
Wavelet transforms are one of the fast-evolving mathematical and signal processing tools 15, 16 .A wavelet transform is complete, orthogonal in the discrete form , and local 17 .A continuous wavelet transform CWT decomposes a function by band-pass filtering of the original signal at different bandwidths, while a discrete wavelet transform DWT is implemented by using quadrature mirror filter QMF banks 18 .The basic operation of wavelet transform consists of the procedures of dilation and translation 19, 20 , which lead to a multiscale analysis of a signal.All these vital advantages make wavelet transforms capable of analyzing nonlinear and nonstationary signals.Nevertheless, the deficiencies including the interference terms, border distortion, and energy leakage may generate a lot of undesired small spikes all over the frequency scales and make the results confusing and difficult to be interpreted.To deal with these shortcomings, Olhede and Walden developed another self-adaptive wavelet-based algorithm also via Hilbert transform, namely, maximal-overlap discrete wavelet packet transform MODWPT 21 .The ordinary DWT requires the sample size to be exactly a power of 2 for the full transform 22 .Besides, the scaling coefficient of DWT is not circularly shift equivariant.By avoiding down sampling, MODWPT overcomes these disadvantages of DWT.With optimum decomposition scale and disjoint dyadic decomposition, the complicated signal could be decomposed into a number of components with instantaneous frequencies physically meaningful at different levels.Furthermore, each single component obtained by MODWPT has desirable statistical characteristics, which are desirable properties to deal with nonstationary time series in practice 21 .
Duffing equation, Lorenz system, and R össler system are three of the typical nonlinear examples.The network traffic data is also one kind of practical nonlinear processes 23-30 .Several TF analyses and spectral analyses in network traffic have been used 31 .To discuss the applicability of nonlinear data analysis on HHT and on MODWPT, we will enumerate these several typical nonlinear systems and examine the effectiveness of these two methods in the following representation.The organization of this paper is given as follows.The rationale of HHT and MODWPT will be elaborated separately in Sections 2 and 3. Characteristics of three typical nonlinear systems using HHT and MODWPT will be discussed in Section 4 and the efficiency of each approach will also be demonstrated.In Section 5, we give the future works that are under investigation and exploration.Finally, we will offer the conclusion in Section 6.

Rationale of HHT
The development of HHT provides an alternative view of the time-frequency-energy paradigm of nonlinear and nonstationary data.To examine data from real-world nonlinear and nonstationary processes, the detailed dynamics in the processes from the data need to be determined because the intrawave frequency modulation, which indicates the instantaneous frequency changes within one oscillation cycle, is one of the typical characteristics of a nonlinear system.As Huang et al. 7 pointed out, the intrafrequency variation is the hallmark of nonlinear systems.One way to express the nonlinearity is to find instantaneous frequency, which reveals the intrawave frequency modulations.But actually, the detailed frequency representation cannot be obtained from a priori approach hampered by a collection of the endless harmonics 6 .Thus, as an easier approach, the Hilbert Transform is used, which is defined as The purpose is to separate function x t into a set of nearly monocomponent signals called IMFs.An IMF is a single frequency component within the length of the signal.The way of extracting instantaneous frequencies is called EMD 7 .

Empirical Mode Decomposition
Physically speaking, the necessary conditions to define a meaningful instantaneous frequency are that the signal must be symmetric concerning the local zero mean, and have the same numbers of zero crossings and extrema.This means that, in an IMF function, the number of extrema and the number of zero crossings must be either equal or different at most by one in the whole data set, and the mean value of the envelope defined by the local maxima and the envelope defined by the local minima is zero at every point.All these conditions are so strict that the determined IMF may not satisfy them precisely.Consequently, the resultant IMF is nearly a monocomponent function.
The EMD is developed based on the assumption that any signal consists of a set of different IMFs.The procedures to decompose signal x t can be enumerated as following steps.
a Find all the local maxima from x t and connect them with the cubic spline to form the upper envelope denoted by x up t .
b Find all the local minima from x t and connect them with the cubic spline to form the lower envelope denoted by and it is designated as: g The criterion suggested by Huang for stopping the sifting process is and the SD value is regularly 0.2∼0.3.
h Remove c 1 t from the rest of the signal by r 1 t x t − c 1 t and treat r 1 t as a new signal and repeating step a to f as described above.Thus, we can obtain a series of IMFs c i i 1, 2, ..., n and the final residue r n t .
Summing up all the IMFs and the final residue, we should be able to reconstruct the original signal x t by 7 Then, the HT of c i t yields a i t e j ω i t dt r n t .

2.8
In the polar coordinates system, x t is expressed by Practically, the residue r n t can be ignored.Let a i ω, t be the combination of the amplitude a i t and the instantaneous frequency ω i t of the ith IMF.The HHT of x t is given by HHT ω, t n i 1 a i ω, t . 2.10

Hilbert Spectrum via MODWPT
Assume that we sample a continuous-time signal at intervals Δt 1 to a sequence of the observation X X 0 , X 1 , . . ., X N−1 and N is a power of 2. For the class of discrete compactly supported Daubechies wavelets, we denote the scaling low-pass filter by {g l : l 0, . . ., L − 1} and the wavelet high-pass filter by {h l : l 0, . . ., L − 1}.These even-length filters satisfy 3.1 for all nonzero integers n, and are related by being quadrature mirror filters:

Mathematical Problems in Engineering
For t 0, . . ., N − 1, the jth level wavelet and scaling coefficients are given by where mod means modulus after division.The maximal overlap discrete wavelet transform MODWT can be considered as a revised version of the DWT 21 .As previously mentioned, the DWT of level j restricts the sample size to a power of 2. However, the MODWT of level j is well defined for any sample size.To conserve energy, we define Thus, 3.1 becomes and the quadrature mirror filters are defined likewise The MODWT creates new filters at each stage by inserting 2 j−1 − 1 zeros between the elements of { g l } and { h l } to avoid downsampling.The MODWT pyramid algorithm generates the MODWT wavelet coefficients {W M j,t } and the scaling coefficients {V M j,t }, respectively by

3.7
The coefficients at level j and frequency-index n can be expressed as W j,n {W j,n,t , t 0, . . ., N − 1}, and then we produce {W j,n,t } using when n mod 4 0 or 3, f n,l g l ; when n mod 4 1 or 2, f n,l h l .For any signal, the analytic form can be represented as s t W j,n t jH W j,n t .

3.9
Then the instantaneous amplitude is denoted by a j,n t W 2 j,n t H 2 W j,n t , 3.10 and the instantaneous phase function is φ j,n t tg −1 H W j,n t W j,n t .

3.11
Accordingly, the instantaneous frequency is f j,n t 1 2π φ j,n t . 3.12

Nonlinear System
Generally speaking, the instantaneous frequency changes within one oscillation cycle for nonlinear systems, and can be used to describe intrawave frequency modulation.To discuss the characteristics of data from a nonlinear system by EMD HT and MODWPT, we will take several typical examples to aid the discussions.

Duffing System
The Duffing oscillator under harmonic excitation described by a second-order differential equation is one of the well-known nonlinear examples where ε, γ constants, and ω harmonic forcing frequency.If we rewrite the equation in the following form the term in parenthesis could be interpreted as nonlinearity in the stiffness of the oscillator, which will lead to a frequency that is ever changing with amplitude from location to location, and time to time, even within a single period.
Embodied by a Duffing oscillator with ε −1 and γ 0.1, the system is allowed to freely vibrate and with the initial conditions of x 0 1 and ẋ 0 1, as shown in Figure 1.After subjected to the EMD, the numerical result yields the IMF components and the corresponding Fourier spectrum of each IMFs as shown in Figure 2, where the first IMF, indicating a concentration of energy near 0.12 Hz, represents the intrinsic frequency of the system.The second IMF, identifying a weak concentration of energy around 0.04 Hz, represents the forcing function, and a low-frequency component, with low energy at 0.017 Hz, represents the very low-intensity subharmonics.Figure 3 illustrates the TF outcome of the Duffing equation using EMD and MODWPT.Apparently, the HHT result reveals the intrinsic frequency clearly, which shows strong introwave frequency modulation, presented as a variable frequency oscillation between 0.06∼0.18Hz.The forcing function is also perfect shown at 0.04 Hz.The low frequency and low amplitude at 0.02 Hz are unexpected but explicable, for they represent the slow aperiodic wobbling of the phase.The Hilbert marginal spectrum of HHT shows the intrawave frequency modulation from 0.06 to 0.18 Hz and the forcing function near 0.04 Hz likewise.Compared to the HHT, while the TF spectrum of MODWPT indicates the frequency oscillation near 0.1 Hz, the forcing function at 0.04 Hz is not displayed, just similar to the Hilbert marginal spectrum on right.

Duffing System
The Lorenz system has also been widely studied and is described by where σ, r, and b positive constants, assumed to be 10, 20, and 3, with initial position of 10, 0, 0 .The numerical result is shown in Figure 4.The IMFs are displayed in Figure 5 with the corresponding Fourier spectrum on right.The comparisons of TF resolutions using HHT and MODWPT are provided in Figure 6.The sharp peak that appears at 1.4 Hz in the Fourier spectrum represents the main oscillating frequency.
In the diagram of HHT, the transient nature of both components is perfectly located, with the main component being intrawave modulated and a fairly clear indication of the nonlinear effect of the oscillation.Comparatively speaking, the result of MODWPT is not so satisfactory.Only some blurry frequency components can be recognized near 1.4 Hz.The oscillation of the frequency, supposed to be demonstrating the nonlinearity, is not represented here in the spectrum of the MODWPT.where μ is a constant parameter which stands for the famous period-doubling event and we let it be μ 3.5.Also we assume the initial condition to be −4,4,0 .Wave form of the x-component is displayed in Figure 7.

Duffing System
The result of the Fourier spectrum of x-component shows many harmonics at 0.27 Hz, 0.35 Hz, and 0.42 Hz, separately 3 times, 5 times, and 7 times of 0.07 Hz.The HHT spectrum and the MODWPT distribution are given in Figure 8 on the left column.The corresponding Hilbert marginal spectrum of each method is on the right column.The marginal spectrum gives a low-frequency peak and a broad bimodal distribution near the main peak frequency, which indicates a typical distribution of a periodic variable.None of the peaks in the marginal spectrum agree with the main peak in the Fourier spectrum in Figure 7. Comparatively, the peak near 0.07 Hz in marginal spectrum of HHT does not appear in the marginal spectrum of MODWPT.Similarly, the high bimodal frequency distribution around 0.18 Hz generated by the intrawave frequency modulation in marginal spectrum of MODWPT agrees with these in marginal spectrum of HHT, which indicates the two time scales involved in the period doubling.As these results showed, the EMD decomposed the data into two nonlinear components more successfully than the MODWPT.

Future Work
This paper focuses on the performance of nonlinear processes using the HHT and the MODWPT with the aim of bringing a better choice of time-frequency decomposition in nonlinear data analysis.Compared with the HHT, the MODWPT shows some weakness in nonlinear data analysis and the study of identifying the computational burden required by these two methods is still at the exploratory stage.

Conclusion
HHT, which decomposes data through EMD, offers a potentially viable method for nonlinear and nonstationary data analysis.Verified by three typical nonlinear systems, the HHT can not only perform and locate main frequency components but also force function frequency details.Furthermore, the intrawave modulation, which is the important characteristic of nonlinear system, can as well be obtained in the distribution of HHT.Compared with the MODWPT that decomposes data into a number of components alike, the EMD gives results much sharper and more supportable to nonlinear system identification.

d 2 x dt 2 Figure 1 :
Figure 1: a Numerical solution of the Duffing equation 200s .b Power spectral density estimate.

bFigure 2 :
Figure 2: The IMF components of the Duffing equation from EMD in a and the corresponding power spectrum in b .

Figure 3 :
Figure 3: Diagrams at left represent the HHT result and the MODWPT result at level 3 of the Duffing equation separately; the corresponding Hilbert Marginal spectrum results of EMD and MODWPT are on the right side.

Figure 4 :
Figure 4: The numerical solution of the Lorenz equation a and the Fourier spectrum of the given xcomponent b .

Figure 5 :
Figure 5: The IMF components of the Lorenz equation from EMD in a and the corresponding power spectrum in b .

Figure 6 :Figure 7 :Figure 8 :
Figure 6: The HHT spectrum a and the MODWPT spectrum at level 3 b for the Lorenz equation solution.
1where P is the Cauchy principal value of the singular integral and in which y t is the Hilbert transform of the function x t .The analytic signal is defined as