An oblique-extrema-based approach for empirical mode decomposition

Digital Signal Processing 20 (2010) 699–714

Contents lists available at ScienceDirect

Digital Signal Processing www.elsevier.com/locate/dsp

An oblique-extrema-based approach for empirical mode decomposition ✩ Zhijing Yang a , Lihua Yang a,∗ , Chunmei Qing b a b

School of Mathematics and Computing Science, Sun Yat-sen University, Guangzhou 510275, China School of Informatics, University of Bradford, Bradford BD7 1DP, UK

a r t i c l e

i n f o

a b s t r a c t

Article history: Available online 12 August 2009 Keywords: Empirical mode decomposition (EMD) Intrinsic mode function (IMF) Inflection point Oblique extremum point Oblique-extrema IMF (OIMF) Oblique-extrema empirical mode decomposition (OEMD)

It has been found that envelopes established by extrema in the empirical mode decomposition cannot always depict the local characteristics of a signal very well. This is due in part to the slight oscillations characterized as hidden scales which are almost left untreated during the sifting process. When involving hidden scales, the intrinsic mode function usually contains at a given instance multiple oscillation modes. In view of this, based on inflection points this paper presents a new decomposition algorithm called ‘oblique-extrema empirical mode decomposition’ to settle these problems. With this algorithm, any signal can be decomposed into a finite number of ‘oblique-extrema intrinsic mode functions’ which may possess better-behaved Hilbert transforms and produce more accurate instantaneous frequencies. It can suppress the effect of hidden scales and gets one step further in extracting finer scales. Experimental results demonstrate good performances of this new method. © 2009 Elsevier Inc. All rights reserved.

1. Introduction Nowadays, nonlinear and non-stationary signal processing is one of the most important research areas in information science. In the past, a non-stationary signal was usually assumed to be piecewise stationary and its processing was based on the classical theories and techniques for stationary signals. Windowed Fourier transform was a typical example based on such an assumption. Only in recent years have new methods been introduced to analyze nonlinear and non-stationary data. For example, wavelet analysis and the Wigner–Ville distribution are designed for non-stationary but linear data. Meanwhile, various nonlinear time series analysis methods (see, for example, [1–3]) are designed for nonlinear but stationary systems. However, data from natural phenomena are usually both nonlinear and non-stationary. To characterize these data, the concept of instantaneous frequency (IF) and its various definitions are introduced. Although none of the existing definitions for IF can give physically meaningful interpretations for all signals, the one based on the analytic signal (AS) generated through the Hilbert transform is commonly acceptable. Furthermore, Vakman showed in [4,5] that the AS is the unique mathematical model to produce physically meaningful IFs. For a multicomponent signal, there may be more than one IF at a time, thus it should be made up of several so-called ‘monocomponents’ whose IFs defined by the AS method provide physically meaningful characterizations of the original signal [6]. In 1998, Huang et al. proposed a novel method called Hilbert–Huang transform (HHT) for analyzing nonlinear and nonstationary signals [7]. It consists of two procedures: the empirical mode decomposition (EMD) and the Hilbert spectrum. With EMD, any complicated data set can be decomposed into a finite and often small number of intrinsic mode functions (IMFs), where an IMF is a function that satisfies the following two conditions:

✩ This work was supported in part by NSFC (Nos. 60873088, 10631080, 60475042, 10926139), NSFGD (No. 9451027501002552), and the Young Faculty Career Start Program under grant 3171921. Corresponding author. Fax: +86 20 84111696. E-mail addresses: [email protected] (Z. Yang), [email protected] (L. Yang), [email protected] (C. Qing).

*

1051-2004/$ – see front matter doi:10.1016/j.dsp.2009.08.013

©

2009 Elsevier Inc. All rights reserved.

700

Z. Yang et al. / Digital Signal Processing 20 (2010) 699–714

Fig. 1. x1 (t ) (solid line) and its inflection points (asterisk).

(1◦ ) In the whole data set, the number of extrema and the number of zero crossings must either equal or differ at most by one; (2◦ ) At any point, the mean value of the envelope defined by the local maxima and the envelope defined by the local minima is zero. For the convenience of depiction, a function which is required to only satisfy the condition (1◦ ) is called a weak-IMF. Applying the Hilbert transform to the IMFs, Huang et al. constructed the energy–frequency–time distribution, designated as the Hilbert spectrum. Recently, the HHT has received more attention in terms of interpretations [8–11] and applications in many disciplines such as ocean science [12], biomedicine [13–15], speech signal processing [16], image processing [17], pattern recognition [18–20] and so on. Though EMD is effective for detecting large and apparent oscillations, it may miss the gentle humps characterized as hidden scales [12,21]. If there is more than one inflection point between two consecutive extrema where the waveforms have unwanted fluctuations, it is reasonable to think that hidden scales exist. In EMD, the characteristic time scale is defined by the time lapse between the extrema. The local mean of the envelopes defined by the extrema is used to characterize the local symmetry of a signal [7]. If the extrema can reveal the changes of oscillation modes, the envelopes defined by them indeed can depict the local characteristics of a signal very well. Unfortunately, not all the changes of oscillation modes produce extrema. If a signal is distorted without introducing any additional extremum, its envelopes defined by the extrema remain the same. This is the main reason why EMD misses the gentle humps and some IMFs contain hidden scales. To show this problem more intuitively, let us consider the following typical signal:

x1 (t ) = t + 0.9 sin(t ),

t ∈ [0, 63].

(1)

The signal is plotted in Fig. 1. Since x1 (t ) is strictly monotonic increasing everywhere, it will be treated as a residue by EMD. However, from either Eq. (1) or Fig. 1 we know that it obviously consists of a harmonic component and a trend r (t ) = t. The amplitude of the harmonic component is so small that x1 (t ) does not produce any extremum. So the slight oscillation mode of that harmonic component is missed by the EMD algorithm. This problem was considered by Huang and his coworkers in [12,21]. They suggested to find the extremum points of the curvature

k(t ) =

x (t )

(1 + (x (t ))2 )3/2

,

and get one step further for EMD in extracting the fine-scale signal from a given data. However they thought this curvaturebased sifting was not necessarily better and recommended to use it only for special cases when the EMD result suggested too many hidden scales. The curvature-based sifting need take the second order derivative of the data. It is well known in mathematics that the numerical computation of high order differential is typically ill-posed [22]. Besides, we observe that there may be some technical problems in employing the extrema of the curvature. Fig. 2(a) shows a monotonic increasing signal, designated as x2 (t ), whose expression is as follows



x2 (t ) =

1 − (t − 2k − 1)2 + 2k, t ∈ (2k, 2k + 1]  − 1 − (t − 2k − 1)2 + 2(k + 1), t ∈ (2k + 1, 2k + 2]

(k = 0, ±1, ±2, . . .).

It is easy to verify that the absolute value of the curvature of x2 (t ) is 1 everywhere. Therefore we cannot get the extrema of the signal and extract the gentle oscillation modes. The following signal is another example,

Z. Yang et al. / Digital Signal Processing 20 (2010) 699–714

701

Fig. 2. (a) A monotonic increasing signal whose absolute value of the curvature is 1 everywhere. (b) A monotonic increasing signal whose curvature has three extrema between two consecutive inflection points, maximum points of the curvature are marked by asterisks and minimum points of the curvature are marked by circles.

 x3 (t ) = 2t +

f (t ), t ∈ (4kπ , 4kπ + 2π ] − f (t ), t ∈ (4kπ + 2π , 4kπ + 4π ]

(k = 0, ±1, ±2, . . .),

where f (t ) is t + sin(t ) on (0, π ] and 2π − t − sin(t ) on (π , 2π ]. It can be easy to verify that it is also monotonic increasing and its curvature has three extrema between any two consecutive inflection points, as shown in Fig. 2(b) where the maximum points of the curvature are marked by asterisks and minimum points of the curvature are marked by circles. In this case the cubic spline interpolation passing all the curvature maxima or minima cannot give a physically meaningful envelope. In this paper, we study this problem seriously. It is known that an inflection point is a separation between the concave and convex of a signal and even a gentle distortion may generate a new inflection point. With inflection points, we can depict the local characteristics of a signal better. However, if we define the local mean by inflection points directly, we cannot guarantee the local symmetry of an IMF. Motivated by the extremum point and based on inflection points, we propose the concept of ‘oblique extremum point’. It is proved in this paper that there exists a unique ‘oblique extremum point’ between any two consecutive inflection points. So the oblique-extrema-based envelopes can characterize the properties of a signal more locally. Based on ‘oblique extremum points’, this paper introduces a new kind of IMF, designated as ‘oblique-extrema IMF’ (OIMF). It is proved that an OIMF is a weak-IMF, i.e. it does not contain any riding wave. What is more, there exists a unique inflection point between its any two consecutive extrema, so it can avoid hidden scales existing. Then, a new algorithm called ‘oblique-extrema empirical mode decomposition’ (OEMD) is presented to decompose any signal into a finite number of OIMFs. With this new decomposition algorithm, we can separate different oscillation modes even hidden scales from one another. Therefore OEMD gets one step further in extracting the finer-scale signal from a given data compared with EMD and OIMFs may produce more accurate instantaneous frequencies. We will show that OIMFs can also admit intrawave modulations and the intra-wave frequency modulation admitted in an OIMF is weaker than that in an IMF. In addition, in order to identify inflection points in the OEMD algorithm stably, we utilize the convolution technique to calculate the numerical differentiation. Encouraging results and new insights are obtained by applying the OEMD to various data: from the typical simulation signals to data representing natural phenomena. The rest of this paper is organized as follows. The definition of OIMF and its properties are presented in Section 2. In Section 3, we review the EMD algorithm and propose a new algorithm called OEMD. Experiments and discussions are conducted to demonstrate the performance of our method in Section 4. Finally, conclusions of this paper are summarized in Section 5. 2. A new kind of IMF: Oblique-extrema IMF 2.1. The definition of oblique-extrema IMF As stated in [7], asymmetric waveforms will induce unwanted fluctuations of the IF. So Huang et al. utilize the condition (2◦ ) to guarantee the local symmetry of an IMF. This condition utilizes the local mean of the envelopes defined by extrema to force the local symmetry. Nevertheless, the envelopes defined by extrema cannot always depict the local characteristics of a signal very well. If an IMF is distorted without introducing any additional extremum, it will still satisfy the condition (2◦ ). But its waveforms may no longer be symmetric. That is the reason why an IMF may contain hidden scales.

702

Z. Yang et al. / Digital Signal Processing 20 (2010) 699–714

Fig. 3. The asterisk mark ‘∗’ represents a local maximum point of f l (t ) and the circle mark ‘◦’ represents an oblique local maximum point of f (t ).

Fig. 4. (a) The oblique extremum point ξ1 (circle) is an extremum point since f (ta ) = f (tb ); (b) the oblique extremum point ξ3 (circle) is different from the extremum point ξ2 (asterisk) since f (ta ) = f (tb ); (c) ξ4 (circle) is an oblique extremum point but there is no extremum point on (ta , tb ).

So we should choose other more powerful characteristic points to define the envelopes. As we know that an inflection point is where the signal changes the concavity (see [23] for the precise definition). It is a separation between the concave and convex of a signal and the joint point between various oscillation modes. With inflection points, we can depict the characteristics of a signal more locally, since even a gentle oscillation may generate a new inflection point. However, if we define the local mean by inflection points directly, we cannot guarantee the local symmetry. Motivated by the extremum point and based on inflection points, we propose the concept of ‘oblique extremum point’ whose precise definition is given as follows. Definition 1 (Oblique extremum point). Let (ta , f (ta )) and (t b , f (tb )) (ta < tb ) be two consecutive inflection points of a curve f (t ) and l(t ) be the straight line connecting these two inflection points. Denote

f l (t ) = f (t ) − l(t ),

t ∈ [ta , tb ],

(2)

then the local maximum (or minimum) point ξ of f l (t ) is called an oblique local maximum (or minimum) point of f (t ). A typical oblique extremum point is shown in Fig. 3, where (ξ, f l (ξ )) (asterisk) is a local maximum point of f l (t ) and (ξ, f (ξ )) (circle) is an oblique local maximum point of f (t ). From Definition 1, it is easy to check that the oblique extremum point is an extremum point if and only if f (ta ) = f (tb ). An illustration is given in Fig. 4. In Fig. 4(a) the oblique extremum point ξ1 is just an extremum point due to the fact that f (ta ) = f (tb ). But if f (ta ) = f (tb ) the oblique extremum point need not be an extremum point even if it exists nearby, as shown in Fig. 4(b) or not, as shown in Fig. 4(c). The latter case attracts our particular attention, because in the monotonic segment from ta to tb of the wave a gentle oscillation mode is included which cannot be extracted by EMD. To extract finer scale of a signal, oblique extremum points will be employed to improve the EMD algorithm instead. It will be shown that there is a unique oblique extremum point between any two consecutive inflection points. Before this, let us clarify two concepts on extremum points and inflection points, which are practically useful in EMD and our algorithm.

Z. Yang et al. / Digital Signal Processing 20 (2010) 699–714

703

b d Fig. 5. (a) Only a+ is considered as an extremum point on [a, b]; (b) the original signal (solid line) where ( c + , f ( c +2 d )) is considered as an inflection 2 2 point; (c) the second order derivative of the signal shown in Fig. 5(b) where it is zero on [c , d].

( 1) The first concept is on identifying extremum points in the EMD and our algorithm. Mathematically t 0 is called an extremum point of a function f (t ) if there exists a neighborhood of t 0 , U (t 0 ), such that f (t 0 )  f (t ) or f (t 0 )  f (t )× (t ∈ U (t 0 )). In practice, for a special case that f (t ) is a constant on an interval [a, b] and there exists a constant δ b such that f (t ) is strictly monotonic on the intervals [a − δ, a] and [b, b + δ], we only consider the midpoint a+ as an 2 extremum point on [a, b] if it has distinct monotonicity on [a − δ, a] and [b, b + δ]; otherwise, we think that there is no b is considered as an extremum point on extremum point on [a, b]. (An example is shown in Fig. 5(a), where only a+ 2 [a, b].) ( 2) The second concept is on identifying inflection points in our algorithm. We extend the definition of the inflection point to ensure that there definitely exists an inflection point between the concave and convex of a function. In mathematics for a twice differentiable function f (t ), (t 1 , f (t 1 )) is called an inflection point if there exists a constant δ1 such that f  (t ) has distinct signs on the intervals (t 1 − δ1 , t 1 ) and (t 1 , t 1 + δ1 ). In practice, for a special case that f  (t ) is zero d , f ( c+2 d )) is considered as an inflection point in our algorithm if there exists a constant δ2 on an interval [c , d], ( c + 2  such that f (t ) has distinct signs on the intervals (c − δ2 , c ) and (d, d + δ2 ) (an example is shown in Fig. 5(b) where ( c+2 d , f ( c+2 d )) is considered as an inflection point and the second order derivative of the signal is shown in Fig. 5(c)). 1), it is obvious that the extremum points of f  (t ) are all the inflection points of f (t ). By ( With these in mind, we will prove an important theorem as follows. Theorem 1. Suppose (ta , f (ta )) and (t b , f (tb )) (ta < tb ) are two consecutive inflection points of a twice differentiable function f (t ). Then there exists a unique oblique extremum point ξ of f (t ) on the open interval (ta , tb ). In addition, ξ is an oblique maximum (or minimum) point if and only if f  (t ) is decreasing (or increasing) on (ta , tb ). Proof. Since f l (ta ) = f l (tb ) = 0, where f l (t ) is defined as Eq. (2), by Rolle’s mean-value theorem, there is at least one point ξ belonging to (ta , tb ) such that f l (ξ ) = 0. By ( 2) and the hypothesis that (ta , f (ta )) and (t b , f (tb )) are two consecutive inflection points of f (t ), we conclude that f  (t ) is monotone on (ta , tb ). So f l (t ) is also monotone on (ta , tb ). By ( 1), ξ is the unique extremum point of f l (t ) on (ta , tb ), that is, ξ is the unique oblique extremum point of f (t ) on (ta , tb ). It is obvious that ξ is an maximum (or minimum) point if and only if f l (t ) is decreasing (or increasing) on (ta , tb ). Since f l (t ) and f  (t ) have the same monotonicity, the proof is complete. 2 From Theorem 1, we know that there exists a unique oblique extremum point between two consecutive inflection points. It also implies that there exists a unique inflection point between two consecutive oblique extremum points. In other words, any change of oscillation modes even a gentle one will produce an oblique extremum point. So the envelopes defined by the oblique extremum points can depict the characteristics of a signal more locally. Furthermore, the local mean of envelops defined in such a way can force better local symmetry and smoothness. This is the motivation for us to present a new kind of IMF, designated as ‘oblique-extrema IMF’ with the following formal definition. Definition 2 (Oblique-extrema IMF (OIMF)). An oblique-extrema IMF (OIMF) is a function that satisfies the following two conditions: (1) All the oblique local maxima are positive and all the oblique local minima are negative. (2) At any point, the mean value of the envelope defined by the oblique local maxima and the envelope defined by the oblique local minima is zero.

704

Z. Yang et al. / Digital Signal Processing 20 (2010) 699–714

As will be proved in the next subsection that a signal satisfying the condition (1) of an OIMF is a weak-IMF and there is a unique inflection point between its two consecutive extrema. So the condition (1) can ensure that no riding waves and hidden scales are allowed in an OIMF. The condition (2) is used to replace the condition (2◦ ) to force the local symmetry and smoothness so that the IF will not have unwanted fluctuations induced by asymmetric waveforms. Furthermore, OIMFs may admit better-behaved Hilbert transforms and produce more accurate IFs by the AS method compared with IMFs. 2.2. The properties of OIMF In this subsection, we will investigate the properties of OIMF. Specifically, we are interested in the relation between OIMF and IMF. It can be proved that a signal satisfying the condition (1) of an OIMF is a weak-IMF and its first order derivative is also a weak-IMF. This can guarantee no riding waves in the signal and the singleness of the oscillation mode at a given instance. The result is proved as follows. Theorem 2. If a twice differentiable function f satisfies the condition (1) of an OIMF, then it is a weak-IMF and its first order derivative f  is also a weak-IMF. Proof. Firstly, we prove that f (t ) is a weak-IMF, which means that all the local maxima of f (t ) are positive and all the local minima of f (t ) are negative. We prove the fact anagogically: suppose there exists a local minimum point t 0 of f (t ) such that f (t 0 )  0. Let (ta , f (ta )) and (t b , f (tb )) be the left and right inflection point of f (t ) closest to t 0 respectively. Then obviously, f (t )  0 on the interval (ta , tb ) and consequently the unique oblique extremum point ξ of f (t ) on (ta , tb ) (according to Theorem 1) satisfies f (ξ )  0. Since f (t ) arrives at local minimum at t 0 , we conclude that f  (t ) is monotone increasing on (ta , tb ). Therefore ξ is an oblique local minimum of f (t ). The fact that f (ξ )  0 contradicts the condition (1) of OIMF. Therefore f (t 0 ) < 0. It can be shown similarly that all the local maxima of f (t ) are positive. Hence f (t ) is a weak-IMF. Secondly, let us prove that f  (t ) is also a weak-IMF. Suppose there exists a local minimum point t 0 of f  (t ) such that f  (t 0 )  0. Let ta and tb be the left and right maximum point of f  (t ) closest to t 0 respectively. Obviously f  (t )  0 and f (t ) is monotonic increasing on the interval [ta , tb ]. By ( 2), we know that ta , t 0 , tb are the abscissas of three inflection points of f (t ). Note that f  (t ) is monotonic decreasing on the interval [ta , t 0 ], by Theorem 1, we conclude that f has a unique oblique local maximum point ξ1 on (ta , t 0 ). Similarly f has a unique oblique local minimum point ξ2 on (t 0 , tb ). Since f (t ) is monotonic increasing on [ta , tb ], it is impossible to ensure f (ξ1 ) > 0 and f (ξ2 ) < 0. It contradicts the condition (1) of OIMF. Therefore f  (t 0 ) < 0. Similarly we can prove that all the local maximum of f  (t ) must be positive. Hence f  (t ) is a weak-IMF. 2 In a signal, there may be more than one inflection point between two consecutive extremum points, which may be just the case of hidden scales. Since it has been shown that there must exist a unique inflection point between two consecutive oblique extremum points, oblique extremum points can characterize the oscillations better than extremum points. Based on oblique extremum points, we will introduce a new algorithm called ‘oblique-extrema empirical mode decomposition’ to decompose any data set into the needed OIMFs in the next section. 3. The oblique-extrema empirical mode decomposition 3.1. The EMD algorithm and discussions Since the ‘oblique-extrema empirical mode decomposition’ is motivated by EMD, let us review the EMD algorithm briefly here. Given a signal x(t ), let s(t ) = x(t ), then the first IMF c 1 (t ) can be extracted by the following iteration. 1) Identify all the extrema of s(t ). 2) Connect all the local maxima of s(t ) by a cubic spline curve as the upper envelope, designated as e + (t ). Repeat the procedure for the local minima to produce the lower envelope e − (t ). 3) Compute the mean m(t ) = (e + (t ) + e − (t ))/2. 4) Sifting: h(t ) = s(t ) − m(t ). 5) If h(t ) is not an IMF, let s(t ) = h(t ) and go back to 1); otherwise, h(t ) is an IMF, then let c 1 (t ) = h(t ) and end the iteration. The residue r1 (t ) = x(t ) − c 1 (t ) is treated as the new data, i.e., s(t ) = x(t ) − c 1 (t ), and subjected to the same iteration as described above. The sifting process is stopped by any of the following criteria: after extracting n-IMFs, the residue, rn (t ) is either an IMF or a monotonic function. Finally, we obtain

x(t ) =

n  i =1

c i (t ) + rn (t ).

(3)

Z. Yang et al. / Digital Signal Processing 20 (2010) 699–714

705

Fig. 6. x1 (t ) (solid line), its inflection points (asterisk) and oblique extremum points (circle).

As analyzed in Section 1, the EMD sifting process based on the envelops defined by the extrema may miss gentle oscillations characterized as hidden scales. Let us consider the signal x1 (t ) discussed in Section 1 again. As shown in Fig. 6, there exists one oblique extremum point (circle) between two consecutive inflection points (asterisk). If decompose x1 (t ) based on oblique extremum points, we can extract the gentle oscillation mode. 3.2. The oblique-extrema empirical mode decomposition Knowing the nice properties of OIMF is only the starting point. Since most of the data are not OIMFs, it is necessary to decompose them into sums of OIMFs. Based on the previous discussions, we know that the decomposition based on oblique extrema is more reasonable. In view of this, we propose a new decomposition algorithm called ‘oblique-extrema empirical mode decomposition’. With this new decomposition algorithm, any complicated data set can be decomposed into a finite number of OIMFs. First of all, it is necessary to design an algorithm to identify oblique extremum points in experiments. According to the definition of the oblique extremum point, we must identify the inflection points of a signal. In practice, we determine the inflection points of a signal by the extrema of its first order derivative. The algorithm of extracting oblique extremum points in experiments is presented as follows. Algorithm 1 (Algorithm of extracting oblique extremum points). Let x(t ) be the original signal. Step 1: Compute the first order derivative of x(t ), x (t ), and identify all the extremum points of x (t ), i.e. the abscissas of all the inflection points of x(t ), designated as a monotonic increasing sequence g 1 , g 2 , . . . , gm . Step 2: Connect all the inflection points of x(t ) by a piecewise line as the trend of x(t ), designated as l(t ). Step 3: Subtract the trend from the original signal: xl (t ) = x(t ) − l(t ) on the closed interval [ g 1 , gm ]. Step 4: Get all the extremum points of xl (t ): X 1 , X 2 , . . . , X m−1 . Then they are all the oblique extremum points of x(t ). The program ends. Remarks. In Step 1, we need to compute the first order derivative of a signal. As we know that the numerical calculation of derivative is ill-posed and unstable. In order to identify inflection points stably, we utilize a stable method to calculate the numerical differentiation in Section 4.4. Now let us present the ‘oblique-extrema empirical mode decomposition’ (OEMD) as follows. Algorithm 2 (OEMD). Given a signal x(t ), let s(t ) = x(t ), then the first OIMF s1 (t ) can be extracted by the following iteration. Step 1: Identify all the oblique extrema of s(t ) by Algorithm 1. Step 2: Connect all the oblique local maxima of s(t ) by a cubic spline curve as the upper envelope, designated as e + (t ). Repeat the procedure for the oblique local minima to produce the lower envelope e − (t ). Step 3: Compute the mean m(t ) = (e + (t ) + e − (t ))/2. Step 4: Sifting: h(t ) = s(t ) − m(t ). Step 5: If h(t ) is not an OIMF, let s(t ) = h(t ) and go back to Step 1; otherwise, h(t ) is an OIMF, then let s1 (t ) = h(t ) and end the iteration.

706

Z. Yang et al. / Digital Signal Processing 20 (2010) 699–714

Fig. 7. (a) The decomposition result of EMD for the signal y 1 (t ). EMD treats y 1 (t ) as a residue. (b) The decomposition result of OEMD for the signal y 1 (t ). Top: the original signal y 1 (t ); middle: an OIMF; bottom: the residue.

Thereafter, the residue r1 (t ) = x(t ) − s1 (t ) is treated as the new data, i.e. s(t ) = r1 (t ), and subjected to the same iteration as described above. The sifting process can be stopped by any of the following criteria: after extracting n-OIMFs, either the residue rn (t ) is an OIMF or the number of oblique extrema of rn (t ) is less than two. Finally, we obtain

x(t ) =

n 

si (t ) + rn (t ).

(4)

i =1

Thus, the signal x(t ) is decomposed into n-OIMFs and a residue. Remarks.

• In this paper, we use the standard deviation (SD) proposed in [24] as follows

N

SD = tN=1

|m(t )|2

t =1 |h (t )|

2

.

SD  10−4 for analog signals and SD  10−3 for noisy signals. • We adopt the boundary conditions proposed in [25] here to treat the end effects of the spline fitting. The sifting process of OEMD has two effects: (a) to eliminate riding waves and hidden scales; and (b) to smooth uneven amplitudes. Furthermore, similar to EMD, the OEMD is intuitive, direct and adaptive, with the basis of the decomposition derived from the data. An OIMF produced by OEMD does not contain hidden scales and its waveforms will be symmetric and smooth. As will be shown in our examples, OIMFs may admit better-behaved Hilbert transforms and produce more accurate instantaneous frequencies. 4. Experiments and discussions In this section, we will present some typical experiments and their analysis to show the effectiveness and performance of OEMD. 4.1. Extraction of hidden scales Rather than the EMD algorithm, OEMD can extract gentle oscillations characterized as hidden scales. Let us consider the signal discussed in Section 1:

y 1 (t ) = t + 0.9 sin(t ),

t ∈ [0, 63].

Since y 1 (t ) is strictly monotonic increasing, EMD treats it as a residue without any oscillation, as shown in Fig. 7(a). However, OEMD can extract the harmonic component from the original signal correctly, as shown in Fig. 7(b).

Z. Yang et al. / Digital Signal Processing 20 (2010) 699–714

707

Fig. 8. (a) The decomposition result of EMD for the signal y 2 (t ; 100, 0.1). EMD considers it as an IMF. (b) The decomposition result of OEMD for the signal y 2 (t ; 100, 0.1). Top: the original signal; middle and bottom: two OIMFs.

Fig. 9. (a) The decomposition result of EMD for the signal y 2 (t ; 100, 0.12). Top: the original signal; middle: two IMFs; bottom: the residue. (b) The decomposition result of OEMD for the signal y 2 (t ; 100, 0.12). Top: the original signal; middle and bottom: two OIMFs.

Since y 1 (t ) does not have any extremum, it may be too special for the EMD algorithm. Therefore, let us consider the following simple model discussed by Rilling and Flandrin in [10]:

y 2 (t ; a, f ) = cos(2π t ) + a cos(2π f t ),

t ∈ R,

where a > 0 and 0  f  1. It consists of two harmonic components, where the cos(2π t ) term is the higher frequency component and a cos(2π f t ) is the lower frequency component. We did a lot of tests on this model and found out that OEMD behaved better than EMD especially in extracting gentle oscillations. Because of the limitation of length, we just give two examples here to illustrate the advantage of our method. When f = 0.1 and a = 100, EMD considers y 2 (t ; 100, 0.1) as an IMF, yet OEMD can separate the two components correctly, as shown in Fig. 8. When f = 0.12 and a = 100, the decomposition results of EMD and OEMD are given in Fig. 9. It can be seen that both EMD and OEMD mainly obtain two components. OEMD separates the two components correctly, while EMD does not. Moreover, OEMD can also extract the gentle oscillations for x2 (t ) and x3 (t ) on the interval [0, 20π ] discussed in Section 1 very well, as shown in Figs. 10(a) and 10(b) respectively, while EMD treats them as residues without any oscillation. Though the above examples are simple, they still can validate the advantage of our method in extracting gentle oscillations. Moreover, we have great confidence to believe that this method can be applied to some practical problems.

708

Z. Yang et al. / Digital Signal Processing 20 (2010) 699–714

Fig. 10. (a) The decomposition result of OEMD for the signal x2 (t ) discussed in Section 1. (b) The decomposition result of OEMD for the signal x3 (t ) discussed in Section 1 on the interval [0, 20π ].

Fig. 11. (a) The signal y 3 (t ) (solid line) and its inflection points (asterisk); (b) y 3 (t ) (solid line) and its inflection points (asterisk) on the interval [0, 6].

4.2. Better-behaved Hilbert transforms of OIMFs It is known that most IMFs admit well-behaved Hilbert transforms and their IFs obtained through the AS are coincident with its essential physics [7]. However, it is also observed that for some IMFs the IF computed by the AS will contain negative values. The first theoretic example was given by Sharpley and Vatchev in [26], whose expression is given by





y 3 (t ) = e 2.97 cos(t ) sin 2.97 sin(t ) ,

t ∈ [−60, 60].

It is proved in [26] that y 3 (t ) is an IMF. Clearly y 3 (t ) is 2π -periodic, so we only plot the original signal on the interval [−4π , 4π ] as shown in Fig. 11(a). It should be noted that hidden scales exist on some intervals, such as the interval [0, 6]. The magnified graph of Fig. 11(a) on [0, 6] is shown in Fig. 11(b). With the Hilbert transform, the IF of y 3 (t ) is calculated as shown in Fig. 12(a). y 3 (t ) and its IF are plotted together in Fig. 12(b) for convenient observations. It can be seen that if there is only one inflection point between any two consecutive extrema of y 3 (t ), i.e. the waveforms do not fluctuate, the IF values of y 3 (t ) are always positive. Nevertheless, if there is more than one inflection point between any two consecutive extrema of y 3 (t ), i.e. hidden scales exist, the IF of y 3 (t ) will have many negative values. To clarify further, we magnify the graph of Fig. 12(b) on the interval [0, 6] as shown in Fig. 12(c). Hence, we guess that the negative frequency may be caused by the hidden oscillation. To support our conjecture, we use OEMD to decompose y 3 (t ) into two OIMFs, whose amplitudes are constant as shown in Fig. 13(a). The first OIMF is a high-frequency oscillation, while the second is a low-frequency one. With the Hilbert transform, we obtain the IFs of the two OIMFs as shown in Fig. 13(b). It can be seen from Fig. 13(b) that the two OIMFs admit better-behaved Hilbert

Z. Yang et al. / Digital Signal Processing 20 (2010) 699–714

709

Fig. 12. (a) The IF of y 3 (t ) (solid line) and x-axis (dashed line); (b) y 3 (t ) (solid line), the IF of y 3 (t ) (dash-dotted line) and the inflection points of y 3 (t ) (asterisk); (c) a magnified display of Fig. 12(b) on the interval [0, 6].

Fig. 13. (a) The decomposition result of OEMD for the signal y 3 (t ). Top: the original signal y 3 (t ); middle and bottom: two OIMFs. (b) The IF of the first OIMF of y 3 (t ) (solid line) and the IF of the second OIMF of y 3 (t ) (dash-dotted line).

transforms without any negative IF value. Besides, the IF of the first OIMF is greater than that of the second one at any time. This gives a physically meaningful interpretation to y 3 (t ). From this example, we can see that hidden scales may be harmful for producing physically meaningful IFs. Therefore OIMFs may admit better-behaved Hilbert transform and produce more accurate IFs compared with IMFs. 4.3. Discussions on intra-wave frequency modulation As stated in [7], there are actually two types of frequency modulations: the inter-wave and the intra-wave modulations. The first type is familiar to us; the frequency of the oscillation is gradually changing with the waves in a dispersive system. The second type is less familiar, but it is also a common phenomenon and essential for nonlinear signals: if the frequency may change from time to time within a wave, its profile can no longer be described by simple trigonometric functions. EMD views most wave profile deformations from the simple sinusoidal form as intra-wave frequency modulation which means that IMFs admit intra-wave modulations. This gives us new physical insight in nonlinear and non-stationary phenomena. In this subsection we will show that OIMFs can also admit intra-wave modulations. Since OEMD can extract gentle oscillations characterized as hidden scales, the waveforms of OIMFs are smoother than IMFs. That is to say the intra-wave frequency modulation admitted in an OIMF is weaker than that in an IMF. An IMF considers the wave profile deformations without introducing any additional extremum as intra-wave frequency modulations, while an OIMF considers such deformations without introducing any inflection point as intra-wave frequency modulations. To show this more intuitionally, let us consider the following mathematical model discussed in [7]:



f ε (t ) = cos



ωt + ε sin(2ωt ) ,

(5)

710

Z. Yang et al. / Digital Signal Processing 20 (2010) 699–714

Fig. 14. (a) The decomposition result of OEMD for the signal f ε (t ) when ε = 0.4 and ω = 1. Top: the original signal; middle and bottom: two OIMFs. (b) Top: the first OIMF of signal f ε (t ) when ε = 0.4 and ω = 1 (solid line) and the component 12 ε cos(3ωt ) (dotted line); bottom: the second OIMF of signal f ε (t ) when

ε = 0.4 and ω = 1 (solid line) and the component (1 − 12 ε ) cos(ωt ) (dotted line). The agreement is excellent for this case.

with ω as a constant frequency and ε as a small parameter. According to the classic wave theory, this expression is a clear case of intra-wave frequency modulation. It is easy to see that when ε  0.5, EMD considers f ε (t ) as an IMF for any ω > 0 [7]. On the other hand, if let ε0 = sup{ε  0.5 | f ε (t )  0, t ∈ (0, 2πω )}, then we can prove that when 0  ε  ε0 OEMD also considers it as an OIMF for any ω > 0. When ε0 < ε  0.5, f ε (t ) is no longer an OIMF for any ω > 0 (the proof is given in Appendix A Proposition 1). For ε0 < ε  0.5, OEMD mainly decomposes f ε (t ) into two OIMFs. As an example, the decomposition result for ε = 0.4 and ω = 1 is given in Fig. 14(a). As stated in [7], the second-order approximation of Eq. (5) is



f ε (t ) = cos



ωt + ε sin(2ωt )

  = cos(ωt ) cos ε sin(2ωt ) − sin(ωt ) sin ε sin(2ωt ) 

1 1 ≈ 1 − ε cos(ωt ) + ε cos(3ωt ) + · · · . 2

(6)

2

It can be seen from Fig. 14(b) that the first OIMF is corresponding to the component 1 2

1 2

ε cos(3ωt ) and the second OIMF is

corresponding to the component (1 − ε ) cos(ωt ). The decomposition result of OEMD is consistent with this mathematical expansion. When ε > 0.5, f ε (t ) is neither an IMF nor an OIMF. The decomposition results of EMD and OEMD are almost the same. As an example, choose ε = 0.8 and ω = 1. The decomposition results of EMD and OEMD are shown in Figs. 15(a) and 15(b) respectively. The two IMFs and OIMFs are almost the same and also very similar to the two components in Eq. (6). To sum up, both IMFs and OIMFs admit intra-wave frequency modulations. The intra-wave frequency modulation admitted in an OIMF is weaker than that in an IMF. The above discussions on intra-wave frequency modulation also show that OEMD gets one step further in extracting the finer-scale signal from the given data compared with EMD. 4.4. Discussions on robustness We know that OEMD needs to compute the first order derivative of a signal, so it will make the method less robust than EMD. To improve the robustness of OEMD, we use a stable +∞method to calculate the numerical differentiation as follows: we smooth the signal first with a convolution f ∗ K  (t ) := −∞ f (x) K  (t − x) dx, where K  (t ) = 1 K ( t ) is a kernel function sat-

∞

isfying −∞ K (x) dx = 1. It is well known in mathematics that f ∗ K  (t ) is as smooth as K (t ) and f ∗ K  (t ) can approximate f (t ) as close as possible when  → 0. Thus, the first order derivative of f ∗ K  (t ) can be employed as an approximation of f  (t ). For practical considering, we derive the discrete algorithm of f ∗ K  (t ) as follows: Given a discrete signal f (n) (n = 1, . . . , N ), assume its values outside the domain are zeros, then 

( f ∗ K  ) (n) =

1



∞

2

−∞

f (n − x) K 

 x



dx

(n = 1, 2, . . . , N ).

(7)

It is easy to see that the numerical differentiation is replaced by numerical integration if the derivative of the kernel function K  (x) is known. Furthermore, using the approximation f (n − x) ≈ f (n − j ) for all x ∈ [ j , j + 1), we have

Z. Yang et al. / Digital Signal Processing 20 (2010) 699–714

711

Fig. 15. (a) The decomposition result of EMD for the signal f ε (t ) when ε = 0.8 and ω = 1. Top: the original signal; middle and bottom: two IMFs. (b) The decomposition result of OEMD for the signal f ε (t ) when ε = 0.8 and ω = 1. Top: the original signal; middle and bottom: two OIMFs.



( f ∗ K  ) (n) ≈

=



2

j =− N

K

f (n − j )

j =−∞

+N 1 



j +1 

+∞ 

1

x



dx

j





 j+1 j f (n − j ) K −K (n = 1, 2, . . . , N ).





√

(8)

In this paper, we choose the Gauss kernel function, i.e. K (t ) = π e −π t . It should be pointed out that the numerical derivative method as Eq. (8) may smooth out some hidden scales. If the smoothing parameter  is large, OEMD will be more robust but less sensitive to hidden scales. Therefore if the data contain low-level noise we choose small  , while if the data contain high-level noise we choose large  . It should be pointed out that both EMD and OEMD algorithms, because of their dependence on the extrema or oblique extrema, are not stable to noise. Indeed OEMD may be less stable to noise than EMD due to the numerical differentiation. However, the numerical differentiation method as Eq. (8) can overcome the drawback somewhat as discussed above. Below we consider an experiment of OEMD and EMD for a noisy signal: Let s(t ) be the signal of sin(t ) plus 45 dB white Gaussian noise. With EMD and OEMD to s(t ) the decomposition results are shown in Figs. 16(a) and 16(b) respectively. It is seen that both EMD and OEMD are stable to noise in this level. If the SNR of the signal decreases OEMD may be less stable than EMD. In conclusion, though OEMD is a little less robust than EMD, it is also stable for certain level noise. Moreover, our experiments show that a smoothness before the decomposition is very helpful to improve the stability to noise for both EMD and OEMD. 2 2

4.5. Discussions and application to a real-world signal The examples presented so far are constructed so that we can analyze and verify the solutions of different algorithms. We now describe the application of OEMD to a real-world signal, the length-of-the-day data, discussed by Huang et al. in [12]. The data note the length of every day covering the period from 1978 to 1988. The original data are shown in Fig. 17(a), where the abscissa represents the days and the ordinate represents the lengths of the corresponding days. Clearly, the data are quite complicated. The decomposition result of OEMD is shown in Fig. 17(b). From Fig. 17(a) it can be seen that an approximatively yearly cycle imbeds in the data. We can immediately see the efficiency of OEMD, since the sixth OIMF clearly characterizes the yearly cycle. 5. Conclusions HHT has found many applications in a variety of problems covering signal processing, pattern recognition, biomedicine engineering and ocean science since its advent. Although it often proves remarkably effective, it has been found that when involving hidden scales, the output of EMD, namely, the IMF, usually contains multiple oscillation modes at a given instance. In view of this, based on inflection points we present a novel decomposition algorithm called OEMD to settle these problems. With this new decomposition algorithm, any signal can be decomposed into a sum of proposed OIMFs. It is proved in this paper that an OIMF is necessarily a weak-IMF, which means it does not have any riding wave. Moreover, there exists a

712

Z. Yang et al. / Digital Signal Processing 20 (2010) 699–714

Fig. 16. (a) The decomposition result of EMD for the resultant noisy signal s(t ). (b) The decomposition result of OEMD for the resultant noisy signal s(t ).

Fig. 17. (a) The length-of-the-day data covering the period from 1978 to 1988. (b) The decomposition result of OEMD for the length-of-the-day data.

unique inflection point between its any two consecutive extremum points, which can avoid hidden scales existing and may produce more accurate IFs. Another benefit of the proposed method is its relative high accuracy and stability by using a convolution kernel to substitute numerical integration for numerical differentiation. Experimental results demonstrate good performances of this new method. We recommend the OEMD as a candidate for further investigation and applications. Acknowledgments The authors would like to thank the anonymous referees for their thorough reviews and constructive comments. Appendix A Proposition 1. Let ε0 = sup{ε  0.5 | f ε (t )  0, t ∈ (0, 2πω )}. For any ω > 0, if 0  ε  ε0 , then f ε (t ) = cos(ωt + ε sin(2ωt )) is an OIMF; if ε0 < ε  0.5, then f ε (t ) is not an OIMF. Proof. It is easy to calculate its first and second order derivatives as follows:

f  (t ) = −ω sin







ωt + ε sin (2ωt ) 1 + 2ε cos(2ωt ) ,

 2  f  (t ) = ω 4ε sin(2ωt ) sin ωt + ε sin (2ωt ) − cos ωt + ε sin(2ωt ) 1 + 2ε cos(2ωt ) .  2







(9) (10)

Z. Yang et al. / Digital Signal Processing 20 (2010) 699–714

From f  (t ) = 0 and

ε  0.5, we know that the extremum points of f (t ) are t = kπ

713

kπ

2π ω (k is an arbitrary integer). Since ω is

the period of f (t ) and f (t ) is symmetrical about t = ω , we only need to consider the inflection point between these two π  π  π  consecutive extremum points: t = 0 and t = π ω . It is easy to check that f ( 2ω ) = 0 and f ( ω − t ) = − f (t ), so ξ0 = 2ω is the abscissa of an inflection point of f (t ). Since f (ξ0 ) = 0, f (t ) is an OIMF if and only if (ξ0 , f (ξ0 )) is the unique inflection point on the interval (0, π ω ). π Since f  (0) < 0 and f  ( π ω ) > 0, (ξ0 , f (ξ0 )) is the unique inflection point on (0, ω ) if and only if the following two inequalities hold:

f  (t )  0





t ∈ 0,





f (t )  0



π 2ω



π π t∈ , 2ω ω

,

(11)

 .

(12)

 Since f  ( π ω − t ) = − f (t ), inequality (11) is equivalent to inequality (12). Therefore (ξ0 , f (ξ0 )) is the unique inflection point ) if and only if inequality (11) holds. on (0, π ω Let ε0 = sup{ε  0.5 | f ε (t )  0, t ∈ (0, 2πω )}. Since for any ω > 0, when ε = 0, f ε (t ) satisfies inequality (11), so ε0  0. For t ∈ (0, 4πω ], if 0  ε < 0.1, it is easy to check that the following two inequalities hold:

0  4ε sin(2ωt )  0.4 < cos 0  sin





π 4

  + 0.1  cos ωt + ε sin(2ωt ) ,



2

ωt + ε sin (2ωt )  1 + 2ε cos(2ωt ) .

That is

f ε (t ) < 0





0  ε < 0.1, t ∈ 0,

π 4ω

 (13)

.

If 0  ε  0.5 and t ∈ ( 4πω , 2πω ), then it is easy to check that sin(2ωt ) > 0, cos(2ωt ) < 0, (ωt + ε sin(2ωt )) ∈ (0, π2 ) and cos(ωt + ε sin(2ωt )) > 0. Therefore for fixed t, sin(ωt + ε sin (2ωt )) is strictly monotone increasing and cos(ωt + ε sin(2ωt ))(1 + 2ε cos(2ωt ))2 is strictly monotone decreasing respect to ε . That is to say for fixed t, f ε (t ) is strictly monoπ tone increasing respect to ε . Furthermore it is easy to check that f 0.077 ( 0.296 ω ) > 0, so 0  ε0 < 0.077. From the definition of ε0 , it is easy to find that if ε0 < ε  0.5 then f ε (t ) dissatisfies inequality (11), i.e. f ε (t ) is not an OIMF. Therefore we just need to prove if 0  ε  ε0 , f ε (t ) satisfies inequality (11), i.e. f ε (t ) is an OIMF. By 0  ε0 < 0.077 and inequality (13), the following inequality holds

f ε (t ) < 0





0  ε  ε0 , t ∈ 0,

π 4ω

 .

(14)

Since f ε (t ) is continuous respect to ε and for fixed t ∈ ( 4πω , 2πω ), f ε (t ) is strictly monotone increasing respect to the definition of ε0 it is easy to check that

f ε0 (t )  0,

t∈

ε , from



π π , . 4ω 2ω

(15)

Therefore

f ε (t )  0



0  ε  ε0 , t ∈

π π , 4ω 2ω

 (16)

.

Combining inequality (14) and inequality (16), the following inequality holds 

f ε (t )  0





π 0  ε  ε0 , t ∈ 0, 2ω



That is if 0  ε  ε0 , f ε (t ) satisfies inequality (11). In conclusion, for any ω > 0, if 0  ε  ε0 , f ε (t ) is an OIMF; if It can be found that

ε0 ≈ 0.075.

(17)

.

ε0 < ε  0.5, f ε (t ) is not an OIMF. 2

714

Z. Yang et al. / Digital Signal Processing 20 (2010) 699–714

References [1] [2] [3] [4] [5] [6] [7] [8] [9] [10] [11] [12] [13] [14] [15]

[16] [17] [18] [19] [20] [21] [22] [23] [24] [25] [26]

C. Diks, Nonlinear Time Series Analysis, World Scientific, Singapore, 1997. H. Kantz, T. Schreiber, Nonlinear Time Series Analysis, Cambridge University Press, Cambridge, 1997. H. Tong, Nonlinear Time Series Analysis, Oxford University Press, Oxford, 1997. D. Vakman, On the definition of concepts of amplitude, phase and instantaneous frequency of a signal, Radio Eng. Electron. Phys. 17 (5) (1972) 754–759. D. Vakman, On the analytic signal, the Teager–Kaiser energy algorithm, and other methods for defining amplitude and frequency, IEEE Trans. Signal Process. 44 (4) (1996) 791–797. L. Cohen, Time-Frequency Analysis: Theory and Applications, Prentice–Hall, Upper Saddle River, NJ, 1995. N.E. Huang, Z. Shen, S.R. Long, M.C. Wu, H.H. Shih, Q. Zheng, N.-C. Yen, C.C. Tung, H.H. Liu, The empirical mode decomposition and the Hilbert spectrum for nonlinear and non-stationary time series analysis, Proc. R. Soc. Lond. A 454 (1998) 903–995. P. Flandrin, G. Rilling, P. Goncalves, Empirical mode decomposition as a filter bank, IEEE Signal Process. Lett. 11 (2) (2004) 112–114. Q. Chen, N. Huang, S. Riemenschneider, Y. Xu, A B-spline approach for empirical mode decompositions, Adv. Comput. Math. 24 (2006) 171–195. G. Rilling, P. Flandrin, One or two frequencies? The empirical mode decomposition answers, IEEE Trans. Signal Process. 56 (1) (2007) 85–95. S. Meignen, V. Perrier, A new formulation for empirical mode decomposition based on constrained optimization, IEEE Signal Process. Lett. 14 (12) (2007) 932–935. N.E. Huang, Z. Shen, S.R. Long, A new view of nonlinear water waves: The Hilbert spectrum, Annu. Rev. Fluid Mech. 31 (1999) 417–457. W. Huang, Z. Shen, N.E. Huang, Y.C. Fung, Engineering analysis of biological variables: An example of blood pressure over 1 day, Proc. Natl. Acad. Sci. USA 95 (1998) 4816–4821. S.C. Phillips, R.J. Gledhill, J.W. Essex, C.M. Edge, Application of the Hilbert–Huang transform to the analysis of molecular dynamic simulations, J. Phys. Chem. A 107 (2003) 4869–4876. Z. Yang, L. Yang, D. Qi, Detection of spindles in sleep EEGs using a novel algorithm based on the Hilbert–Huang transform, in: Proceeding of the 4th International Conference on Wavelet Analysis and Its Applications, Book on Applied and Numerical Harmonic Analysis, Macau, SRA, China, Nov. 30– Dec. 2, 2005. Z. Yang, D. Huang, L. Yang, A novel pitch period detection algorithm based on Hilbert–Huang transform, Lect. Notes Comput. Sci. 3338 (2004) 586–593. N. Bi, Q. Sun, D. Huang, Z. Yang, J. Huang, Robust image watermarking based on multiband wavelets and empirical mode decomposition, IEEE Trans. Image Process. 16 (8) (2007) 1956–1966. J.C. Nunes, S. Guyot, E. Delechelle, Texture analysis based on local analysis of the bidimensional empirical mode decomposition, Mach. Vision Appl. 16 (3) (2005) 177–188. Z. Yang, D. Qi, L. Yang, Signal period analysis based on Hilbert–Huang transform and its application to texture analysis, in: Proceedings of the Third International Conference on Image and Graphics, Hong Kong, China, 18–20 Dec. 2004, pp. 430–433. Z. Yang, L. Yang, D. Qi, C.Y. Suen, An EMD-based recognition method for Chinese fonts and styles, Pattern Recogn. Lett. 27 (2006) 1692–1701. N.E. Huang, M.-L.C. Wu, S.R. Long, S.S.P. Shen, W. Qu, P. Gloersen, K.L. Fan, A confidence limit for empirical mode decomposition and Hilbert spectral analysis, Proc. R. Soc. Lond. A 459 (2003) 2317–2345. A.N. Tikhonov, V.Y. Arsenin, Solutions of Ill-posed Problems, Wiley, New York, 1977. S.M. Nikolsky, A Course of Mathematical Analysis, vol. 1, Mir, Moscow, 1977. N.E. Huang, S.S.P. Shen, The Hilbert–Huang Transform and Its Applications, World Scientific, New Jersey, 2005. G. Rilling, P. Flandrin, P. Gongalves, On empirical mode decomposition and its algorithms, in: IEEE-EURASIP Workshop on Nonlinear Signal and Image Processing NSIP-03, Grado (I), 2003. R.C. Sharpley, V. Vatchev, Analysis of the intrinsic mode functions, Constr. Approx. 24 (1) (2006) 17–47.

Zhijing Yang was born in Guangdong, China, in 1980. He received the B.S. degree and Ph.D. degree in mathematics and computing science from Sun Yat-sen University, China, in 2003 and 2008, respectively. Now he is a lecturer at the School of Mathematics and Computing Science, Sun Yat-sen University, China. His current research interests include signal processing, time-frequency analysis and pattern recognition. Lihua Yang received the B.S. degree from the Mathematics Department, Hunan Normal University, China, in 1984, the M.S. degree from the Mathematics Department, Beijing Normal University, China, in 1987, and the Ph.D. degree from the Department of Scientific Computing and Computer Application, Sun Yat-sen University, China, in 1995. From 1996 to 1998, he was a postdoctoral fellow in the Institute of Mathematics, Academia Sinica, China. As a visiting scholar he worked in the Department of Computer Science of Hong Kong Baptist University from 1997 to 1999, in the Center of Pattern Recognition and Machine Intelligence, Concordia University, Canada, from July of 2001 to January of 2002, and in the Department of Mathematics, Syracuse University, USA, from August to November 2006. He has published more than 40 papers and 1 book on approximation theory, wavelets, image processing, and pattern recognition. He is now a professor at the School of Mathematics and Computing Science, Sun Yat-sen University, China. His current research interests include wavelet and time-frequency analysis, signal processing, and pattern recognition. Chunmei Qing was born in Sichuan, China, in 1980. She received the B.S. degree in mathematics and computing science from Sun Yat-sen University, China, in 2003. She is now pursuing the Ph.D. degree at the School of Informatics, University of Bradford, Bradford, UK. Her current research interests include image/video processing, time-frequency analysis, and pattern recognition.