An improved envelope algorithm for eliminating undershoots

Digital Signal Processing 23 (2013) 401–411

Contents lists available at SciVerse ScienceDirect

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

An improved envelope algorithm for eliminating undershoots ✩ Lijun Yang a,c , Zhihua Yang b , Lihua Yang a,∗ , Ping Zhang d a

Guangdong Province Key Laboratory of Computational Science, School of Mathematics and Computational Science, Sun Yat-sen University, Guangzhou, 510275, China Information Science School, Guangdong University of Business Studies, China c College of Mathematics and Information Science, Henan University, China d Department of Mathematics and Computer Science, Alcorn State University, USA b

a r t i c l e

i n f o

Article history: Available online 31 August 2012 Keywords: Instantaneous frequency Envelope Undershoot Empirical Mode Decomposition (EMD)

a b s t r a c t According to the drawbacks of current envelope algorithms, we present an improved envelope algorithm based on cubic spline and monotone piecewise cubic polynomial interpolations. The new envelope can eliminate the undershoots and meanwhile keep the smoothness property. In addition, we show that the developed method is valid when it is applied to the Empirical Mode Decomposition for non-stationary signal processing. © 2012 Elsevier Inc. All rights reserved.

1. Introduction

Signal analysis can be conducted in time domain and/or frequency domain. The Fourier transform and its inverse transform have built a bridge between the analysis on these two domains. The Fourier transform decomposes a signal into the components of different scales in the frequency domains. But this presentation lacks local time information of the signal. Namely, for a certain frequency spectrum, it cannot be mapped to an exact time segment in which the frequency occurs. Essentially, the Fourier analysis is the best tool for analyzing and processing stationary signals. However, most signals that appear in practical applications are nonstationary. It is obviously not efficient and reasonable to describe non-stationary signals only on time domain or frequency domain. This motivates the research of time–frequency analysis [3]. The basic objective of time–frequency analysis is to devise a function, called time–frequency distribution, which can describe energy density of a signal simultaneously in time and frequency [1]. Since the frequencies of a non-stationary signal vary over time, we call these frequencies Instantaneous Frequencies (IF). The concept of IF has been extensively studied for more than a half century. But so far, IF is still far from well-defined. A widely used definition of IF is to employ the analytic signal [3].

✩ This research was partially supported by National Natural Science Foundation of China (Nos. 11071261, 91130009, 60873088, 60802061), the Computational Science Innovative Research Team program and Guangdong Province Key Laboratory of Computational Science at the Sun Yat-sen University. Corresponding author. Fax: +86 20 84111696. E-mail address: [email protected] (L. Yang).

*

1051-2004/$ – see front matter © 2012 Elsevier Inc. All rights reserved. http://dx.doi.org/10.1016/j.dsp.2012.08.011

Let us recall the relevant concepts on IF. Given a signal x(t ), its Hilbert transform is defined as the following Cauchy principal value integral







x(t − u )

H x(t ) = p. v.

u

du .

By taking x(t ) as the real part and H [x(t )] as the imaginary part we can obtain the analytic signal





Z (t ) = x(t ) + j H x(t ) = a(t )e j φ(t ) , where



a(t ) =

 

2

x2 (t ) + H x(t )

,

 H [x(t )] . φ(t ) = arctan x(t )

(1)

(2)

The functions a(t ) and φ(t ) in Eq. (2) are called analytic envelope and analytic phase of x(t ), respectively. The derivative of the analytic phase

ω(t ) =

dφ(t ) dt

(3)

is defined as the IF of the signal x(t ) [10]. The definition of IF in Eq. (3) can indeed describe the frequency physical-reasonably for many signals including those with form cos(ωt + θ0 ), where ω , θ0 are real constants. But for many other signals, it leads to some physical paradoxes [1]. Therefore, we usually consider the signals in a suitable scope whose IF can be reasonably computed. To use the definition IF in Eq. (3), the signals should be mono-component, that means the signal should have only one oscillation mode at any time [1]. However, the concept of mono-component is not rigorous, and it is hard to say whether a signal is mono-component or not. Therefore we can only consider the IF for special signal classes in the practical applications.

402

L. Yang et al. / Digital Signal Processing 23 (2013) 401–411

Fig. 1. The signals (black line) and their corresponding analytic envelopes (red line). (a) The signal s1 (t ). (b) The signal s2 (t ). (For interpretation of the references to color in this figure legend, the reader is referred to the web version of this article.)

Fig. 2. The signals (black line) and their corresponding CSEs (red line). (a) The signal s1 (t ). (b) The signal s2 (t ). (For interpretation of the references to color in this figure legend, the reader is referred to the web version of this article.)

One typical such signal class is the narrow-band signals, having the form of a(t ) cos(ωt + φ(t )), where the frequency of a(t ) is far smaller than ω and φ(t ) is a small perturbation [11]. We know that there are infinite pairs of functions (a(t ), φ(t )) to write a given signal x(t ) in the form of a(t ) cos(ωt + φ(t )), it is hence difficult for us to judge whether the signal is narrow-band or not. As a result, it cannot be asserted the efficiency of the method for computing the IF of a signal. The envelope plays an important role in studying the IF of non-stationary signals. Once the upper and lower envelopes are obtained and both are symmetric about the x-axis, we can easily calculate the IF of the signal. Rice proposed an envelope model for narrow-band signals of the form x(t ) = a(t ) cos[ωt + φ(t )] [10]. However, this envelope model is descriptive and the basic requirement for an envelope is to coincide with the physical intuition of human beings. In fact, Rice did not explain what “physical intuition” exactly means in his literature. It is usually a common sense that an envelope should meet the following two conditions: (1) an envelope must be smooth and it contains as little oscillation as possible; (2) the envelope wraps the signal as tightly as possible. The envelope should be situated above (below) the signal except for the points where the envelope and the signal are tangent. The envelopes may exist even if a signal is not mono-component or narrow-band, in which case we cannot compute IF. In [5], envelopes are used in the Empirical Mode Decomposition (EMD), which emerged as a new method for processing the non-stationary signals. The mathematical mean is computed using the upper and lower envelopes during the processing of EMD. The Intrinsic Mode Functions (IMFs) are obtained by subtracting the mean from the

signal iteratively. By using EMD, the Hilbert spectrum, an important new time–frequency distribution, can be obtained [5]. By calculating the cubic spline interpolation iteratively the envelope of each IMF can be obtained and consequently a new method to calculate the IF of the IMF, called AM/FM demodulation, is deduced in [6,7]. There are currently two main solutions to calculate envelopes, one is analytic envelope method based on the Hilbert transform [3,4,12–14]. The analytic envelope usually does not conform to the physical character of a signal. We take two examples to illustrate this point. The first example is a multi-component signal s1 (t ) = cos t + cos(10t )/4, whose analytic envelopes are shown in Fig. 1(a). As can be seen, the resulting envelopes make little physical sense. The second example is a narrow-band signal s2 (t ) = (3 + 2 cos t ) cos(6t + cos(3t )), and its analytic envelopes are shown in Fig. 1(b). It can be seen that the envelopes are not smooth, which does not conform to the basic physical intuition. The other envelope model is based on cubic spline interpolation [5], which is derived by fitting a cubic spline through the local maxima/minima and denoted by “CSE” for simplicity in this paper. The CSEs of s1 (t ) and s2 (t ) are shown in Fig. 2. It can be seen that the CSEs have much better intuitive character than the analytic envelopes. Although the CSEs perform well, there frequently appears “undershoot” near the maximum/minimum points of the signal. We take the signal x(t ) = 3 cos(π t ) + cos(4π t ) as an example, whose CSEs are shown in Fig. 3. From the enlarged drawing on the topright corner, it can be seen that the upper/lower envelope does not wrap the signal completely and there are apparent undershoots, such as the points “A” and “B”. Therefore the upper envelope

L. Yang et al. / Digital Signal Processing 23 (2013) 401–411

u (t ) < x(t ), and the lower envelope v (t ) > x(t ) at some locations. If the upper envelope and lower envelope are reflection of each other, i.e. v (t ) = −u (t ), it will lead to the contradiction: cos θ(t ) > 1 when computing IF by θ(t ) = cos−1 (x(t )/u (t )), ω(t ) = θ  (t ). In order to have |x(t )/u (t )|  1, Huang et al. proposed an empirical AM/FM demodulation method through a normalization of iterative modification [6]. This method needs many times of multiplication operations, which lead too much computational complexity. What is worse, the normalized phase signal may contain riding waves and cause negative IFs [16]. Qin et al. studied the problem of undershoot appearing in the CSEs [9]. They thought the reason of undershoot is due to the inadequate flexibility of the spline functions. Based on this idea, they proposed a new envelope model, which uses the piecewise power function to construct the envelope. Unfortunately, smoothness and flexibility cannot be met simultaneously in their method. That is, in order to have a better flexibility to eliminate the undershoot phenomenon, the envelope cannot have good smoothness. Furthermore, the authors neither gave the method to find the order of the power function for the envelope construction nor provided the convincing analysis for the efficiency of the method. We will show later that this method cannot effectively solve the problem of undershoot. Yang et al. presented an envelope operator model using the basic conditions for envelope and proved that the zero-mean

condition implies the zero-crossing condition in the definition of IMF [15]. However, a specifically feasible envelope still has not been given. In this paper, we propose a new algorithm for calculating the envelope, which can eliminate the undershoots and meanwhile keep the smoothness property. The improved envelope can also meet the physical intuition. Moreover, we will show that the developed method is valid when it is applied to the EMD. Because of no undershoot on the envelope, the instantaneous amplitude and instantaneous frequency of each IMF can be obtained directly without normalization as described in [7]. The rest of the paper is organized as follows: In Section 2, a new envelope algorithm is presented and the corresponding computational complexity issue is discussed. In Section 3 we use two typical signals to compare the new envelope algorithm with other methods. In Section 4, a new EMD is proposed based on the improved envelope, and consequently, in Section 5, a new time– frequency distribution is given based on the new envelope. Finally, Section 6 is the conclusion of the paper. 2. A new envelope algorithm In this section, we will propose a new envelope algorithm. We only discuss the upper envelope since the low envelope can be handled similarly. To begin with, an example signal is shown in Fig. 4(a). The blue curve in Fig. 4(a) represents an original signal x(t ) and the green-dash line is the CSE u 0 (t ) of x(t ). Although u 0 (t ) has good physical intuition in the nearly entire region, undershoots occur in a few locations. The enlarged drawing of the circle region in Fig. 4(a) is shown on the top-right corner of Fig. 4(a), from which we can clearly see that an undershoot appears between ta j and tb j . To solve this problem, we should modify the envelope u 0 (t ). The goal is not only to eliminate undershoots but also to: (i) maintain initial envelope smoothness as well as possible; (ii) do not change the initial envelope as far as possible in those regions without undershoots. For this we should determine the undershoot regions first. It can be implemented through the following straightforward procedure. Let

x1 (t ) = x(t ) − u 0 (t ), Fig. 3. The signal x(t ) (blue line) and its CSEs (red-dotted line). (For interpretation of the references to color in this figure legend, the reader is referred to the web version of this article.)

403

(4)

which is shown in Fig. 4(b), and x1 (t ) > 0 in the undershoot intervals, such as the interval from ta j to t b j in Fig. 4(b). Consequently, the undershoot intervals can be represented by

Fig. 4. (a) The signal x(t ) = 3 cos(π t ) + cos(4π t ) (blue line), the CSE u 0 (t ) (green-dash line), (b) the difference signal x1 (t ) = x(t ) − u 0 (t ) (blue line). (For interpretation of the references to color in this figure legend, the reader is referred to the web version of this article.)

404

L. Yang et al. / Digital Signal Processing 23 (2013) 401–411

Fig. 5. (a) u 1 (t ) (red-dotted line) is the upper envelope of x1 (t ) (blue line) by the MPCI with all nonnegative maxima. (b) The difference signal x2 (t ) = x(t ) − [u 0 (t ) + u 1 (t )] (blue line). (For interpretation of the references to color in this figure legend, the reader is referred to the web version of this article.)



 (tai , tbi ): x1 (tai ) = x1 (tbi ) = 0, x1 (t ) > 0, t ∈ (tai , tbi ), i ∈ N , (5)

where N is the set of all the natural numbers. The second step is to modify the envelope near the undershoot regions, so as to make the final envelope more tightly wrap the original signal and meanwhile affect other regions as little as possible. Let us consider the nonnegative local maxima of x1 (t ) and divide them into two sets:













O = O i toi , x1 (toi ) , x1 (toi ) = 0, i = 1, 2, . . .



(6)

and



N = N j tn j , x1 (tn j ) , x1 (tn j ) > 0, j = 1, 2, . . . .

(7)

The points in O are also the maxima of x(t ) where x(t ) is tangent with u 0 (t ). As shown in Fig. 4(b), the point O i with coordinate (1, 0) is one of such points. It is not in the undershoot regions and should be affected as little as possible when modifying the envelope. The points in N locate in the undershoot regions, such as the point N j with x-coordinate t = 2.5 in Fig. 4(b), where the difference between x(t ) and u 0 (t ) arrives at maximum within the undershoot interval [ta j , tb j ]. Using monotone piecewise cubic interpolation (MPCI) [2] through the points in O ∪ N, we can obtain the upper envelope u 1 (t ) of x1 (t ), which is shown in Fig. 5(a). It obviously holds that u 1 (t )  0 at all t. Similarly by x2 (t ) we denote the residual between x1 (t ) and u 1 (t ), that is

x2 (t ) = x1 (t ) − u 1 (t ).

(8)

Now we have x1 (t )  u 1 (t ) at all t as shown in Fig. 5(b). Substituting (4) into (8) and introducing u (t ) = u 0 (t ) + u 1 (t ), we have





x2 (t ) = x(t ) − u 0 (t ) + u 1 (t ) = x(t ) − u (t )  0.

(9)

That is x(t )  u (t ). Therefore the undershoot intervals (tai , tbi ), i = 1, 2, . . . , have been eliminated if replacing u 0 (t ) with u (t ) as the new envelope. According to the properties of MPCI, the improved envelope u (t ) changes the initial envelope u 0 (t ) slightly in those regions without undershoot. It is easy to see that u (t ) is the first differentiable since u 0 (t ) and u 1 (t ) are the second and first differentiable, respectively. The new envelope u (t ) is shown in Fig. 6. However, u 1 (t ) may still have some undershoots for normal signals. As an example, let us consider the signal x(t ) = −2 sin t − sin(5t ). The difference x1 (t ) = x(t ) − u 0 (t ) is shown in Fig. 7(b). It can be easily seen that u 0 (t ) as the upper envelope has obvious undershoots. Using MPCI to interpolate the nonnegative maxima of x1 (t ), we obtain the upper envelope u 1 (t ) of x1 (t ), and the difference signal x2 (t ) = x(t ) − [u 0 (t ) + u 1 (t )], which is shown in Fig. 7(c). Since there still exist undershoots in x2 (t ),

Fig. 6. Signal and its new upper envelope u (t ) = u 0 (t ) + u 1 (t ).

with the same process we can obtain u 2 (t ) and the difference x3 (t ) = x(t ) − [u 0 (t ) + u 1 (t ) + u 2 (t )]. Fig. 7(d) shows that the undershoots have been completely eliminated. Hence the final envelope is obtained as u (t ) = u 0 (t ) + u 1 (t ) + u 2 (t ), which is shown with the red curve in Fig. 7(a). We can see that the improved envelope maintains good smoothness and meets the physical intuition. It is a natural question to ask whether all the undershoots can be eliminated after finite procedures. The stop criterion is xk+1 (t ) = xk (t ) − uk (t )  0. Since uk (t ) is determined by the nonnegative maxima of xk (t ), it can be written as the action of some operator T k on xk (t ): uk (t ) = T k xk (t ) (k = 1, 2, . . .). Therefore we have

xk+1 (t ) = xk (t ) − uk (t ) = ( I − T k )xk (t )

= · · · = ( I − T k )( I − T k−1 ) · · · ( I − T 1 )x1 (t ),

(10)

where I is the identity operator. Since the location of maxima of each xk (t ) is different from those of xk−1 (t ), the operator { T k } is nonlinear and therefore the convergence of {xk (t )} is hard to prove rigorously in mathematics. In spite of this, it is easily seen that uk (t )  0, which implies that {xk (t )} is decreasing. Besides, since u 0 (t ) is the CSE of x(t ), there holds u 0 (t ) < x(t ) only in the undershoot regions and the difference x1 (t ) = x(t ) − u 0 (t ) above x-axis is very small. In addition, the region of xk+1 (t ) = xk (t ) − uk (t ) > 0 shrinks rapidly and meanwhile the corresponding positive maxima decrease rapidly. This phenomenon is not difficult to be understood intuitively from Fig. 7. In fact, lots of experiments show that, for general signals, all the undershoots can be eliminated after a few times procedures as above, usually 2–3 times. The final envelope is

u (t ) = u 0 (t ) + · · · + uk−1 (t )  x(t ).

(11)

L. Yang et al. / Digital Signal Processing 23 (2013) 401–411

405

Fig. 7. (a) x(t ) = −2 sin(t ) − sin(5t ) (blue line) and its upper envelope u (t ) = u 0 (t ) + u 1 (t ) + u 2 (t ) (red line); (b) the difference signal x1 (t ) = x(t ) − u 0 (t ); (c) the difference signal x2 (t ) = x(t ) − [u 0 (t ) + u 1 (t )]; (d) the difference signal x3 (t ) = x(t ) − [u 0 (t ) + u 1 (t ) + u 2 (t )]. (For interpretation of the references to color in this figure legend, the reader is referred to the web version of this article.)

Let us consider the possible cases that undershoots appear. Suppose t i −1 , t i and t i +1 are three consecutive maximum points of x(t ). It can be concluded that there is no undershoot at t i if x(t i )  min{x(t i −1 ), x(t i +1 )}. Otherwise undershoots may exist on the left or/and right side(s) of t i . A particular case is that there are two continual undershoots appearing on both sides of t i , which is called “twin-undershoot” briefly. Fig. 8 is an example of such case, in which the blue-solid line presents the original signal x(t ), the green-dash line shows the CSE u 0 (t ). From the enlarged drawing at the top-right corner, we can see that a twin-undershoot appears at tb j . The point t b j is also the minimum of x1 (t ) and satisfies x1 (tb j ) = 0 which is remarked with “B” in Fig. 8. Using MPCI, the envelope u 1 (t ) of x1 (t ) is obtained and shown with the red-dash line in this figure. Finally, u (t ) = u 0 (t ) + u 1 (t ) is the improved envelope shown in the top chart of Fig. 8. The improved envelope no longer has undershoot. At the point of twin-undershoot, the new envelope cannot tightly wrap the original signal, but arise up a little (see u (t ) near t = 2 in Fig. 8). We can modify the above algorithm to make the envelope wrap the original signal tightly. We describe the modification with an example. In Fig. 9, x(t ), u 0 (t ) and x1 (t ) are the same as those in Fig. 8. Firstly, we search for the zero-value minima of x1 (t ), for instance, the point “B” in Fig. 9. For each such minimum, find the nearest two neighboring maxima on left and right around it, such as the points “A” and “C” in Fig. 9. The interval [ A , C ] is called twin-undershoot interval. Secondly, in order that the envelope wrap the signal tightly on [ A , C ], we define u 1 (t ) := x1 (t ) for t ∈ [ A , C ]. For the other intervals, we still use the MPCI to calculate u 1 (t ) as discussed above. Finally, we can obtain the upper

envelope u 1 (t ) of x1 (t ) by linking the curves, which is the reddotted line shown in Fig. 9. From the properties of MPCI, u 1 (t ) is differentiable in the splicing position, so the revised envelope u (t ) = u 0 (t ) + u 1 (t ) is smooth. Since u (t ) = x(t ) on [ A , C ], it looks like that the envelope “presses closely” against the original signal on this interval. This technique for twin-undershoot is called “splicing processing”. Each type of envelope we mentioned above, as shown in Figs. 8 and 9, has its rationality, they both satisfy the physical intuition of envelopes. In practical applications, one can choose whether or not to apply splicing processing for twin-undershoot envelope according to the requirements. For a given x(t ), the new envelope algorithm can be summarized as follows: Algorithm 1. Suppose x(t ) is a real signal,

• Step 1: interpolate the maxima of x(t ) with cubic spline to obtain u 0 (t ). Let u (t ) = u 0 (t ) and x1 (t ) = x(t ) − u (t ). If x1 (t )  0 for all points, go to Step 4; else, go to Step 2.

• Step 2: interpolate the maxima of x1 (t ) satisfying x1 (t )  0 (for the twin-undershoot envelopes one can choose whether or not to apply splicing processing according to practical requirements) with MPCI to obtain u 1 (t ). Let u (t ) = u (t ) + u 1 (t ), x1 (t ) = x(t ) − u (t ). • Step 3: if x1 (t )  0 for all points, go to Step 4; else, go to Step 2. • Step 4: u (t ) is the final upper envelope. The program ends.

406

L. Yang et al. / Digital Signal Processing 23 (2013) 401–411

Fig. 8. The top chart contains the original signal x(t ) (blue), the cubic spline envelope u 0 (t ) (green) and the final envelope u (t ) (red). The bottom chart contains the difference signal x1 (t ) and its MPCI envelope u 1 (t ). A twin-undershoot appears in the interval [taj , t c j ], in which “B” is a minimum point of x1 (t ) whose value is zero. (For interpretation of the references to color in this figure legend, the reader is referred to the web version of this article.)

Fig. 9. The top chart contains the original signal x(t ) (blue), the cubic spline envelope u 0 (t ) (green) and the envelope u (t ) (red) obtained by Algorithm 1 in which the splicing processing is implemented. The bottom chart contains the difference signal x1 (t ) and its envelope u 1 (t ). (For interpretation of the references to color in this figure legend, the reader is referred to the web version of this article.)

We can obtain the lower envelope similarly, or equivalently, find it as follows: let y (t ) = −x(t ), find the upper envelope u (t ) of y (t ) by Algorithm 1, then v (t ) = −u (t ) is the lower envelope of x(t ).

As an improvement of the CSE, the proposed envelope needs higher computational cost. To find the upper envelope of a signal containing N maxima, CSE needs to solve a linear system of N − 1 equations which is well-conditioned: symmetric, tridiagonal

L. Yang et al. / Digital Signal Processing 23 (2013) 401–411

407

Fig. 10. x(t ) = 2 sin(t ) + 0.6 sin(5t ), t ∈ [0, 5π ] (blue line), the PE (black line); the CSE (red-dotted line). (For interpretation of the references to color in this figure legend, the reader is referred to the web version of this article.)

and diagonally dominant [8]. Thus the computational complexity of CSE is O ( N ). For the improved envelope obtained by Algorithm 1, the first step is to compute the CSE, whose computational complexity is O ( N ). The rest steps use MPCI to interpolate the nonnegative maxima of the difference signal iteratively. The computational complexity of the k-th MPCI is O ( N k ), where N k is the number of positive maxima of the difference signal for the k-th iteration. There holds N > N 1 > · · · > N k → 0, but it is a difficult theoretic problem how fast N k decrease to 0 as the number of iterations k goes. However, a lot of experiments show that it usually takes K = 2–3 times MPCI iteration to eliminate all the undershoots. Thus, the Kcomputational complexity of the improved envelope is O ( N ) + k=1 O ( N k ), that is still O ( N ) for bounded K . 3. Examples compared with other methods In this section, we show the efficiency of Algorithm 1 and compare it with other envelope methods through several examples. In the following, the proposed envelope based on Algorithm 1 will be denoted by “PE” for simplicity. Example 1. Consider the signal x(t ) = 2 sin(t ) + 0.6 sin(5t ), t ∈ [0, 5π ], Fig. 10 compares the two different envelopes obtained by PE and CSE. In the figure, the blue-solid line, the black-solid line and the red-dotted line represent x(t ), the PE of x(t ) and the CSE of x(t ), respectively. From the enlarged drawing of the circled part on the top-right corner of Fig. 10, it is obviously that the CSE has some undershoots, but the PE does not. Example 2. We still use the same signal to compare Algorithm 1 with the method in [9]. Qin et al. proposed an envelope based on the segment power function instead of cubic spline, which they think having the better flexibility. However, they did not give convincing argument to support their ideas. In the following we will show that the power function envelope cannot remove undershoots completely. For simplicity we denote the envelope obtained by the method in [9] by SPFE. From the enlarged drawings (a) and (b) of the circles “A” and “B” in Fig. 11, it is easy to see that the SPFE does not remove the undershoot regardless of whether β = 2.0, β = 2.5 or β = 2.75, where β represents the order of the segment power function for calculating SPFE and β = 2.5 is recommended in [9]. Finally, let us employ a practical signal to compare PE and CSE. Example 3. The solid line in Fig. 12(a) is a section of brain wave intercepted from the SMNI electroencephalogram database and the time length is 0.5 s. The black-dash line is the PE and the reddotted line is the CSE. Fig. 12(b) is the enlarged drawing of the

Fig. 11. x(t ) = 2 sin(t ) + 0.6 sin(5t ) (blue line); the PE (black line); the SPFE with β = 2.0, 2.5, 2.75 by dash line, dot line and dash–dot line, respectively. (a) The enlarged drawing of the circle “A”; (b) the enlarged drawing of the circle “B”. (For interpretation of the references to color in this figure legend, the reader is referred to the web version of this article.)

dash line square of Fig. 12(a), from which it is clearly shown that the CSE contains many undershoots, but the PE eliminates them completely. 4. EMD based on the new envelope The EMD algorithm iteratively obtains IMFs by subtracting the mean of the upper and lower envelopes from the input signal [5]. Each IMF is mono-component and can be written as the form imfk (t ) = ρk (t ) cos θk (t ), where ρk (t ) and θk (t ) are the amplitude and phase of imfk (t ), respectively. Therefore the signal x(t ) can be represented as the form

x(t ) =

ρk (t ) cos θk (t ) + r (t ),

k

where r (t ) is the residual. Because the amplitude ρk (t ) and the phase θk (t ) are time-varying, the EMD can be regarded as an extension of the classical Fourier decomposition. EMD was designed to decompose a multi-component signal into IMFs, which are not predefined but produced adaptively according to the intrinsic physical property of the signal. Each IMF can efficiently describe the frequency attribute of the input signal locally on time domain [5]. In EMD, CSE is employed to compute the mean in each iteration. It inevitably produces undershoots, and even worse it possibly produces pseudo waves because of the “elasticity” of cubic spline [16]. As an example, let us consider the signal x(t ) = exp(20 cos t ) sin(20 sin t ), |t |  12.55, which is shown in the first row in Fig. 13(a), the curves in other rows represent the IMFs and the residual by EMD. The oscillations of the original signal mainly concentrate on several very short time intervals and vanish elsewhere. However, from the results of EMD, we can see that the first component imf1 oscillates with big amplitude in the entire duration, including those intervals where the original signal does not oscillate. Those are typical pseudo waves. Similarly, pseudo waves occur in the second component imf2 and the third component imf3 to counteract and amend those in imf1 . This phenomenon shows that traditional EMD based on CSE has the insufficiency in time–frequency localization of signals. If we use the new envelopes proposed in this paper in the EMD, the results can be improved significantly, which is shown in Fig. 13(b). Let us consider the undershoot problem again. The blue line in Fig. 14(a) is the imf1 in Fig. 13(a), the red and violet lines

408

L. Yang et al. / Digital Signal Processing 23 (2013) 401–411

Fig. 12. (a) A section of brain wave (solid line), “PE” (black-dash line) and “CSE” (red-dotted line); (b) the enlarged drawing of the dash line square in (a). (For interpretation of the references to color in this figure legend, the reader is referred to the web version of this article.)

Fig. 13. The EMD of the signal x(t ) = exp(20 cos t ) sin(20 sin t ), |t |  12.55. (a) EMD based on the cubic spline envelope. (b) EMD based on Algorithm 1.

Fig. 14. (a) The imf1 in Fig. 13(a) and its CSE. (b) The imf1 in Fig. 13(b) and its PE. (For interpretation of the references to color in this figure, the reader is referred to the web version of this article.)

are its upper and lower CSEs. The right figure corresponds to the first IMF in Fig. 13(b) and its envelopes obtained with the proposed algorithm. We can see that the envelope in the left figure has undershoots in many locations, on the contrary, the right side envelope has no undershoot. 5. Time–frequency analysis based on the new envelope Any IMF derived from the EMD can be seen as a monocomponent and then be demodulated to find its instantaneous frequency and amplitude. Furthermore the Hilbert spectrum of the signal can be computed [5]. The most classical method to compute

IF is through the analytic signal based on the Hilbert transform. As shown in Section 1, it is not always consistent with essential physics of the signal. Therefore, Huang et al. proposed a normalization scheme in [7] to compute the IF of each IMF, called empirical AM/FM demodulation. The basic idea is to calculate the CSEs of the IMFs, consequently, regarding the envelope as the signal amplitude, then the cosine of the phase function can be obtained by dividing the signal by the amplitude. However because of the possible undershoots from CSEs, the following normalization scheme is applied to eliminate undershoots in [7]. Suppose that x(t ) is an IMF. Find all the local maxima of |x(t )| and connect all these points with a cubic spline curve to produce

L. Yang et al. / Digital Signal Processing 23 (2013) 401–411

409

Fig. 15. The time–frequency distributions of x(t ) = exp(20 cos t ) sin(20 sin t ), |t |  12.55. (a) TFDold . (b) TFDnew .

the empirical envelope of |x(t )|, denoted by e 1 (t ). Then this envelope can be used to normalize x(t ) by

y 1 (t ) =

x(t ) e 1 (t )

(12)

.

Ideally | y 1 (t )|  1 since e 1 (t ) is the upper envelope of |x(t )|. However, because of possible undershoots in e 1 (t ), there may exist some local extreme points t i (i = 1, 2, . . .), such that | y 1 (t )| > 1 near these points. Repeat the above normalization procedure:

y 2 (t ) =

y 1 (t ) e 2 (t )

,

...,

yn (t ) =

yn−1 (t ) en (t )

,

(13)

until there exists some yn (t ) satisfying | yn (t )|  1. Finally we have

x(t ) = e 1 (t )e 2 (t ) · · · en (t ) yn (t ) = a(t ) cos φ(t ),

(14)

where, a(t ) = e 1 (t )e 2 (t ) · · · en (t ) is assigned as the empirical envelope of x(t ), and yn (t ) = cos φ(t ) is the carrier wave. The key of the above iteration is to ensure the final yn (t ) satisfying | yn (t )|  1, so that yn (t ) can be represented as yn (t ) = cos φ(t ), where φ(t ) is the phase function of x(t ), and the derivative φ  (t ) is the IF. For each imf j we compute the empirical envelope a j (t ) and IF ω j (t ) through the above method, then a time–frequency distribution based on EMD and Empirical Demodulation is obtained and denoted as follows:



0



TFDold (ω, t ) =

if J ω,t is an empty set, j ∈ J ω,t a j (t ),

otherwise,

(15)

where J ω,t = { j | 0  j  n satisfying φ j (t ) = ω}. For each IMF x(t ) obtained with our new EMD described in Section 4, let u (t ) and v (t ) be its upper and lower envelopes by Algorithm 1. By the symmetry of the upper and lower envelopes, we have that v (t ) = −u (t ) (this equality holds approximately in practical applications). Since the new algorithm eliminates undershoots, which implies that v (t )  x(t )  u (t ), i.e. |x(t )|  u (t ), we can use the following equalities:

cos θ(t ) =

x(t ) u (t )

,

θ(t ) = cos−1

x(t ) u (t )

(16)

to compute the carrier wave and the phase function θ(t ). Then the derivative θ  (t ) is the IF. It should be pointed out that we do not need to do any normalization here. Suppose that imf1 , . . . , imfm are the IMFs generated through the new EMD. For i = 1, 2, . . . , m, let u i (t ) be the corresponding envelope of imfi , and θi (t ) be the IF computed by (16). Then we can

obtain the following new time–frequency distribution based on the proposed new envelope algorithm:



TFDnew (ω, t ) =

0



if K ω,t is an empty set, i ∈ I ω,t

u i (t ),

otherwise,

(17)

where I ω,t = {i | 0  i  m satisfying θi (t ) = ω}. Fig. 15 shows the time–frequency distributions TFDold (a) and TFDnew (b) of the signal x(t ) = exp(20 cos t ) sin(20 sin t ), |t |  12.55, computed by (15) and (17), respectively. Obviously, in some intervals where the original signal does not contain any frequency component, TFDold contains several horizontal frequency lines of about 15–25 Hz, but TFDnew hardly contains any frequency, which coincide with the frequency structures of x(t ) well. It shows that TFDnew can provide better characterization of time–frequency localization than TFDold . Let us compare the computational complexity of the demodulation for IMFs obtained with the classical EMD and the new EMD. For each IMF x(t ) obtained with the classical EMD, the demodulation needs to compute the CSE ek (t ) of | yk−1 (t )| and the quotient yk−1 (t )/ek (t ) at each iteration of the normalization as described by Eqs. (12) and (13). If x(t ) is a discrete signal of length n with N extrema (N < n), the computational complexity for CSE is O ( N ) and that for the division is O (n). Therefore the computational complexity for each iteration of the normalization is O (n). If it takes M i iterations to meet | y M i (t )|  1 for the i-th IMF, then the computational complexity for the demodulation of this IMF is M i O (n). Accordingly, the total r computational complexity for the computation of TFDold is ( i =1 M i )O (n) provided that the classical EMD produces r IMFs. For the computation of TFDnew , even though the new EMD spends more computation to generate the IMFs, as discussed above, we can demodulate each IMF directly without normalization, which save much computational complexity. At the end of this section, let us discuss about the effectiveness of the new EMD and TFDnew for noisy signals. Fig. 16 shows the experimental result of EMD and TFDnew of a noisy signal. The first row of Fig. 16(a) is the signal x(t ) = exp(20 cos t ) sin(20 sin t ), |t |  12.55, plus a white Gaussian noise whose signal-to-noise ratio (SNR) is 8.5 in decibels. The rest nine rows are its eight IMFs and the residual decomposed by the new EMD. It is easy to see that noise is successfully separated as the first IMF imf1 , as shown on the second row of Fig. 16(a). The other seven IMFs, imf2 , . . . , imf8 and the residual mainly come from the original signal x(t ). Fig. 16(b) depicts the corresponding time–frequency distribution TFDnew . It contains oscillations with frequency 0–50 Hz

410

L. Yang et al. / Digital Signal Processing 23 (2013) 401–411

Fig. 16. The EMD and time–frequency distribution of signal x(t ) = exp(20 cos t ) sin(20 sin t ), |t |  12.55, plus a white Gaussian noise whose signal-to-noise ratio is 8.5 in decibels. (a) The new EMD. (b) TFDnew .

at five locations, which are mainly contributed by the original signal x(t ). The other frequency spectra distributed on the whole time–frequency plane come from the white Gaussian noise. This time–frequency distribution is consistent with the physical frequency structures of the signal. 6. Conclusion Although some empirical algorithms of computing envelope were widely applied in engineering problems, they have many drawbacks. Envelopes obtained from them may contradict basic physical intuition, therefore leading to inaccurate even wrong results. In this paper, we proposed an improved envelope algorithm which can agree with physical intuition well and meanwhile keep the smoothness property. Compared with the traditional methods, the new envelope completely eliminates undershoots. Moreover, the classical EMD and Hilbert spectrum are improved consequently based on the new envelope. This paper also discusses the computational complexity and effectiveness for noisy signals of the new method. Experimental results presented in this paper confirm that the proposed algorithm performs well. Acknowledgments The authors would like to acknowledge the anonymous referees for their valuable remarks and comments to this work. References [1] L. Cohen, Time–Frequency Analysis, Prentice–Hall, Englewood Cliffs, NJ, 1995. [2] F.N. Fritsch, R.E. Carlson, Monotone piecewise cubic interpolation, SIAM J. Numer. Anal. 17 (2) (1980) 238–246. [3] D. Gabor, Theory of communication, J. Inst. Electr. Eng. 93 (1946) 429–457. [4] J. Huang, L. Yang, Vakman’s analysis in L 2 (R), Int. J. Comput. Math. 88 (3) (2011) 545–554. [5] N.E. Huang, Z. Shen, S.R. Long, et al., 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. [6] N.E. Huang, Z. Wu, S.R. Long, On instantaneous frequency, in: Workshop on the Recent Developments of the Hilbert–Huang Transform Methodology and Its Application, Taipei, China, March 15th–17th, 2006. [7] N.E. Huang, Z. Wu, S.R. Long, K.C. Arnold, X. Chen, K. Blank, On instantaneous frequency, Adv. Adapt. Data Anal. 2 (2009) 177–229. [8] D. Kincaid, W. Cheney, Numerical Analysis: Mathematics of Scientific Computing, third edition, Brooks/Cole Publishing Co., 2002.

[9] S.R. Qin, Y.M. Zhong, A new envelope algorithm of Hilbert–Huang Transform, Mech. Syst. Signal Process. 20 (2006) 1941–1952. [10] S.O. Rice, Envelope of narrow-band signal, Proc. IEEE 70 (7) (1982) 692–699. [11] M. Schwartz, W.R. Bennett, S. Stein, Communications Systems and Techniques, John Wiley and Sons, New York, 1995. [12] 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. [13] D.E. Vakman, On the definition of concepts of amplitude, phase and instantaneous frequency of a signal, Radio Eng. Electron. Phys. 17 (5) (1972) 754–759. [14] D.E. Vakman, L.A. Vainshtein, Amplitude, phase, frequency—fundamental concepts of oscillation theory, Sov. Phys. Usp. 20 (12) (1977) 1002–1016. [15] Z. Yang, L. Yang, A new definition of the intrinsic mode function, in: Proceedings of World Academy of Science, Engineering and Technology, vol. 60, 2009, pp. 822–825. [16] Z. Yang, L. Yang, C. Qing, D. Huang, A method to eliminate riding waves appearing in the empirical AM/FM demodulation, Digital Signal Process. 18 (2008) 488–504.

Lijun Yang was born in Henan province, China, in 1979. She received the B.S. degree in the College of Mathematics and Information Science, Henan University, in 2002. Now she is pursuing the Ph.D. degree at the mathematics and computing science, Sun Yat-sen University. Her current research interests include wavelet and signal processing and time– frequency analysis. Zhihua Yang was born in Hunan province, China, in 1964. He received his M.S. degree in automation from Hunan University and Ph.D. degree in computer science from Sun Yat-sen University respectively in 1995 and 2005. Now he is a professor at Information Science School, Guangdong University of Business Studies. His research interests include Signal Analysis, Pattern Recognition and Image Processing. Lihua Yang received the B.S. degree from Hunan Normal University, China, in 1984, the M.S. degree from Beijing Normal University, China, in 1987, and the Ph.D. degree from 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 Hong Kong Baptist University from 1997 to 1999, in Concordia University, Canada, from July of 2001 to January of 2002, in Syracuse University, USA from August to November, 2006, in Universite Pierre et Marie Curie, Paris, France, from September to December, 2008, in Michigan State University, USA, from June to September, 2010, and in City University of Hong Kong from January to March, 2009 and 2010. He is now a Professor and the director of the Institute of Computing Science and Computer Application, School of Mathematics and Computing Science, Sun Yat-sen University, China. His current research interests include wavelet and time–frequency

L. Yang et al. / Digital Signal Processing 23 (2013) 401–411

analysis, signal processing and pattern recognition. He has published more than 60 papers, 1 book and 2 translations. He was the supervisor of 10 Ph.D. candidates. Ping Zhang (Ph.D.) is an Assistant Professor with the Department of Mathematics and Computer Science, Alcorn State University, USA. As a se-

411

nior IEEE member, he has been working in academy and industry for more than twenty years. He has published more than 40 academic papers and book chapters. He has reviewed more than 60 research journal papers for Pattern Recognition, Pattern Recognition Letters, and IEEE Transactions on SMC, NN, etc. His research interests include Pattern Recognition, OCR, Image Processing and Computer Vision.