Download Empirical Wavelet Transform: A Novel Approach for Signal Decomposition and more Thesis Computer-Aided Power System Analysis in PDF only on Docsity!
IEEE TRANSACTIONS ON SIGNAL PROCESSING, VOL. 61, NO. 16, AUGUST 15, 2013 3999
Empirical Wavelet Transform
Jérôme Gilles
Abstract— Some recent methods, like the empirical mode decom- position (EMD), propose to decompose a signal accordingly to its contained information. Even though its adaptability seems useful for many applications, the main issue with this approach is its lack of theory. This paper presents a new approach to build adaptive wavelets. The main idea is to extract the different modes of a signal by designing an appropriate wavelet filter bank. This construction leads us to a new wavelet transform, called the empirical wavelet transform. Many experiments are presented showing the useful- ness of this method compared to the classic EMD.
Index Terms— Adaptive filtering, empirical mode decomposi- tion, wavelet.
I. I NTRODUCTION
A
DAPTIVE methods to analyze a signal is of great in- terest to find sparse representations in the context of com- pressive sensing. “Rigid” methods, like the Fourier or wavelets transforms, correspond to the use of some basis (or frame) de- signed independently of the processed signal. The aim of adap- tive methods is to construct such a basis directly based on the in- formation contained in the signal. A well known way to build an adaptive representations is the basis pursuit approach which is used in the wavelet packets transform. Even though the wavelet packets have shown interesting results for practical applications, they still are based on a prescribed subdivision scheme. A com- pletely different approach to build an adaptive representation is the algorithm called “Empirical Mode Decomposition” (EMD) proposed by Huang et al. [9]. The purpose of this method is to detect the principal “modes” which represent the signal (roughly speaking, a mode corresponds to a signal which have a com- pactly supported Fourier spectrum). This method has gained a lot of interest in signal analysis this last decade, mainly because it is able to separate stationary and non-stationary components from a signal. However, the main issue of the EMD approach is its lack of mathematical theory. Indeed, it is an algorithmic ap- proach and, due to its non-linearity, is difficult to model. Never- theless, some experiments [5]–[7] show that EMD behaves like an adaptive filter bank. Some recent works attempt to model EMD in a variational framework. In [4], the authors proposed to model a mode as an amplitude modulated-frequency modu- lated (AM-FM) signal and then use the properties of such sig- nals to build a functional to represent the whole signal. Then they are able to retrieve the different modes by minimizing this
Manuscript received October 15, 2012; revised February 10, 2013; accepted May 21, 2013. Date of publication May 30, 2013; date of current version July 19, 2013. The associate editor coordinating the review of this manuscript and approving it for publication was Prof. Trac D. Tran. This work was partially founded by the following grants NSF DMS-0914856, NSF DMS-1118971, ONR N00014-08-1-1119, ONR N0014-09-1-0360, ONR MURI USC, the UC Lab Fees Research and the Keck Foundation. The author is with the Department of Mathematics, University of California, Los Angeles (UCLA), Los Angeles, CA 90024 USA (e-mail: jegilles@math. ucla.edu). Digital Object Identifier 10.1109/TSP.2013.
functional. Another proposed variational approach is the work of Hou et al. [8] where the authors also use the AM-FM for- malism. They propose to minimize a functional which is build on some regularity assumptions about the different components and uses higher-order total variation priors. In this paper, we propose a new approach to build adaptive wavelets capable of extracting AM-FM components of a signal. The key idea is that such AM-FM components have a com- pact support Fourier spectrum. Separating the different modes is equivalent to segment the Fourier spectrum and to apply some filtering corresponding to each detected support. We will show that it is possible to adapt the wavelet formalism by considering distinct Fourier supports and then build a set of functions which form an orthonormal basis. Based on this construction, we pro- pose an empirical wavelet transform (and its inverse) to analyze a signal. The remainder of the paper is organized as follows. Section II has two distinct subsections: in Section II-A, we recall the prin- ciple of the EMD algorithm and the AM-FM model; while in Section II-B, we recall some wavelet formalism which will be useful in our own construction and we discuss some of the ex- isting adaptive wavelet methods. In Section III, we build the proposed empirical wavelets and give some of their properties, then the empirical wavelet transform and its inverse are intro- duced. Section IV show many experiments based on simulated and real signals. The time-frequency representation based on the Hilbert transform is introduced in Section V. In Section VI, we address the question of the estimation of the number of modes. An extension to 2D signals (images) is presented in Section VII. Finally, we conclude and give some perspectives in Section VIII.
II. E XISTING A PPROACHES
A. Empirical Mode Decomposition In 1998, Huang et al. [9] proposed an original method called Empirical Mode Decomposition (EMD) to decompose a signal into specific modes (we define the meaning of “mode” here- after). Its particularity is that it does not use any prescribed function basis but it is self adapting accordingly to the analyzed signal. In this paper, as we will use the Fourier formalism in Section III, we adopt the description used in [4] which is slightly different from the original used in [9]. EMD aims to decompose a signal as a (finite) sum of Intrinsic Mode Functions (IMF) such that
An IMF is an amplitude modulated-frequency modulated func- tion which can be written in the form
1053-587X/$31.00 © 2013 IEEE
4000 IEEE TRANSACTIONS ON SIGNAL PROCESSING, VOL. 61, NO. 16, AUGUST 15, 2013
Fig. 1. EMD: basic IMF detection. Envelopes detection on top (thin contin- uous: , dashed: and and thick continuous: ). On bottom: the first IMF candidate.
The main assumption is that and vary much slower than
. The IMF behaves as a harmonic component. Originally, the method of Huang et al. [9] to extract such IMFs is a pure algorithmic method. Candidates for an IMF are extracted by first computing the upper, , and lower, , envelopes via a cubic spline interpolation from the maxima and minima of . Then the mean envelope is obtained by computing and finally the candidate by (see Fig. 1). Generally, does not fulfill the properties of an IMF. A good candidate can be reached by iterating the same process to and the subsequent. The final retained IMF is. Then the next IMF is obtained by the same algorithm applied on. The remaining IMFs can be computed by repeating this algorithm on the successive residues. The interesting fact about this algorithm is that it is highly adaptable and is able to extract the non-stationary part of the original function. However, its main problem is that it is based on an ad-hoc process which is mathematically difficult to model. Consequently it is difficult to really understand what the EMD provides. For example, some problems appear when some noise is present in the signal. To deal with this problem, an Ensemble EMD (EEMD) was proposed in [17]. The authors propose to compute several EMD decompositions of the original signal corrupted by different artificial noises. Then the final EEMD is the average of each EMD. This approach seems to stabilize the obtained decomposition but it increases the computational cost. Another EMD approach is proposed in [8]. The authors pro- posed to minimize a functional which looks for a sparse repre- sentation of in a dictionary of IMFs. This variational method provides similar results as the original EMD algorithm. How- ever, this functional is based on a scheme which uses higher order total variation terms, this makes the method sensitive to the presence of noise and some filtering must be added to the method.
Fig. 2. On top: dyadic wavelet tiling of the frequency line. On bottom: a wavelet packet like tiling.
B. Wavelets Approaches Nowadays, wavelet analysis is one of the most used tool in signal analysis. Let us fix some notations and recall the very basics about wavelet theory. For further details, we refer the reader to the extensive literature about the wavelet theory, see for example [3], [10], [12], [14]. The Fourier transform and its inverse are denoted and , respectively. In the temporal do- main, a wavelet dictionary is defined as the dilated, with a parameter , and translated by of a mother wavelet (of zero-mean) as
(3)
Then the wavelet transform of is obtained by computing the inner products. If is a continuous vari- able then is called the continuous wavelet transform while if then is called the discrete wavelet transform. A useful property of the wavelet transform is that it can be viewed as the application of a filter bank (each filter corresponds to one scale). In practice, the most used case is the dyadic case,. It can be shown that such a case cor- responds to tile the time-frequency plane like in top of Fig. 2. As we are interested in developing adaptive representations, we recall some existing tentatives of adaptive wavelets con- struction. As far as we know there are a very few attempts in the literature. Probably the most known method is the wavelet packets in a basis pursuit framework based on successive scale refinements of the expansion. It provides an adaptive time-fre- quency plane tiling like in bottom of Fig. 2. Even though the wavelet packets are useful in many applications, they use a con- stant prescribed ratio in the subdivision scheme, which limits their adaptability. Another approach, called the Malvar-Wilson wavelets [10], [13], tries to build an adaptive representation by segmenting the temporal signal itself in order to separate the time intervals con- taining different spectral information. While the original idea is interesting, it turns out that the temporal segmentation is a dif- ficult task to perform efficiently. In [15], the authors propose a method, called the brush- lets, which aims to build an adaptive filter bank directly in the Fourier domain. Basically, it uses the idea of the Malvar-Wilson
4002 IEEE TRANSACTIONS ON SIGNAL PROCESSING, VOL. 61, NO. 16, AUGUST 15, 2013
assume that the number of segment, , is given (in Section VI we propose a method to estimate the number of bands). This implies that we need a total of boundaries, but 0 an are always used in our definition and consequently we need to find extra boundaries. To find such boundaries, we first detect the local maxima in the spectrum and sort them in decreasing order (0 and are excluded). Let us assume that the algorithm found maxima. Two cases can appear:
- : the algorithm found enough maxima to define the wanted number of segments, then we keep only the first maxima,
- : the signal has less modes than expected, then we keep all the detected maxima and reset to the appropriate value. Now, equipped with this set of maxima plus 0 and , we define the boundaries of each segment as the center between two consecutive maxima.
B. Frame
The following proposition shows that, by properly choosing the parameter , we can obtain a tight frame. Proposition 1: If , then the set is a tight frame of. Proof: We follow the idea behind the construction of Meyer’s wavelet. The set is a tight frame if
Accordingly to the periodicity (see Fig. 5), it is enough to focus on the interval. Following the previous notations, we can write
where is a copy of but centered at instead
. First, it is easy to see that for we have
Then, it remains to look at the transition areas. Because of prop- erties of , this result also holds in if consecutive do not overlap:
(13) (14) (15)
This condition must be true for all which is equivalent to said that condition (15) must be true for the smallest and finally
we get the result if.
Fig. 5. Periodicity of the filter bank.
Fig. 6. Example of a fourier partitioning of an empirical filter bank (see text for details).
Fig. 6 gives an empirical filter bank example based on the set with (the theory tells us that ).
C. Empirical Wavelet Transform From the previous section, we know how to build a tight frame set of empirical wavelets. We can now define the Em- pirical Wavelet Transform (EWT), , in the same way as for the classic wavelet transform. The detail coefficients are given by the inner products with the empirical wavelets:
and the approximation coefficients (we adopt the convention to denote them) by the inner product with the scaling function:
where and are defined by (9) and (8), respectively. The reconstruction is obtained by
Following this formalism, the empirical mode , as defined in Section II-A, is given by
(22) (23)
GILLES: EMPIRICAL WAVELET TRANSFORM 4003
Fig. 7. test signal: signal on top. From second to last row: the different components constituting Sig.
IV. E XPERIMENTS
A. Test Signals
We propose to test the Empirical Wavelet Transform on four different signals: three are artificial signals taken from [8] and one real electrocardiogram signal. Their description is given hereafter. Simulated : The first test signal, Sig1 , is made with the sum of three distinct components (for ) (see Fig. 7):
(24) (25) (26)
and then
(27)
Simulated : The second test signal, Sig2 , is made with the sum of three distinct components (for ) (see Fig. 8):
(28)
Fig. 8. test signal: signal on top. From second to last row: the different components constituting Sig.
if otherwise
and then
(31)
Simulated : The third test signal, Sig3 , is made with three distinct components (for ) (see Fig. 9):
(32)
and then
(35)
GILLES: EMPIRICAL WAVELET TRANSFORM 4005
Fig. 14. Modes extracted by the Empirical Wavelet Transform for.
Fig. 15. Modes extracted by the Empirical Wavelet Transform for.
Fig. 16. Sum of the two first EWT modes of.
surprising as those modes have significant individual energy and can be considered as independent modes. About , another phenomenon appears. The initial part is decomposed as the sum of two modes (this sum is given in Fig. 16). If we look more closely, the detected boundaries, shown in Fig. 12, we can
Fig. 17. Modes extracted by the Empirical Wavelet Transform for.
see that the scaling function is built upon a very small support while the next support contains the main information of. If we ask the EWT to consider only (as it supposed to be from the construction of ), the first two supports shown in Fig. 12 are merged except that its boundary is moved closer to the zero frequency. This change in the boundary position has the consequence that some information of is now on the second support and then appear in the second mode instead the first one. In some sense, if we ask to use , we get a perturbed de- composition. This issue comes from the method we use to detect the boundaries of the Fourier support and suggests more inves- tigations. Experiments on the real ECG signal seem to give the advan- tage to the EWT because the EMD provides too many modes. Typically, the EMD modes six to nine are really difficult to in- terpret as such behavior is clearly not visible in the signal itself. A contrary, the EWT focuses on the oscillating patterns we can observe in. Of course it will be very interesting to have the opinion of a cardiologist about the medical interpretation of such components. The different modes for the seismic signal are
4006 IEEE TRANSACTIONS ON SIGNAL PROCESSING, VOL. 61, NO. 16, AUGUST 15, 2013
Fig. 18. Empirical Mode Decomposition of.
not provided here but we will present its time-frequency repre- sentation in Section V as it more relevant for analysis purposes.
V. T IME -F REQUENCY R EPRESENTATION The time-frequency representation is useful to have the infor- mation of all components summarized in a single domain. In this paper, we follow the idea used in the Hilbert-Huang transform [9]. First, let us recall the definition of the Hilbert transform of a function :
where the integral is defined by using the Cauchy principal value [11]. The Hilbert transform can be used to derive the ana- lytical form of :. The Hilbert trans- form has the interesting property that in the case of AM-FM signals , it provides which permits to extract the instantaneous amplitude and frequency. The Hilbert-Huang transform consists in ap- plying the Hilbert transform to each extracted IMF and then plotting each curve in the time-frequency plane where the intensity of the plot is given by. For a time-frequency rep- resentation based on the EWT, we can apply the same process by evaluating the Hilbert transform of each filter output. For instance, Figs. 22 and 23 give the time-frequency representa- tion of and based on the EMD and EWT, respec- tively. We can notice that in the case of the Hilbert-Huang trans-
Fig. 19. Empirical Mode Decomposition of.
form some artifacts appear in the low frequencies while repre- sentations based on the EWT avoid such problems. These ar- tifacts are mainly due to the fact that the EMD decomposition forces the extraction of IMFs even if the initial components are not. Fig. 24 gives the time-frequency representations obtained from the EMD and the EWT of the seismic signal given in Fig. 11. For such signals the EWT seems to provide a more consistent representation during the earthquake while the EMD shows a sparse representation which is difficult to interpret.
VI. A UTOMATIC D ETECTION OF THE N UMBER OF M ODES
If for a given class of signals it should be possible to guess the best number of mode , this is not the case in general where no a priori information are available. In such cases, it should be inter- esting to estimate the appropriate number of modes. In general such task is difficult, hereafter we present a simple way to esti- mate but a deeper analysis is necessary to provide a more ro- bust method. We follow the idea that the most important maxima
4008 IEEE TRANSACTIONS ON SIGNAL PROCESSING, VOL. 61, NO. 16, AUGUST 15, 2013
Fig. 22. Time-frequency representation based on the EMD (Hilbert-Huang transform) for signals (top) and (bottom).
the supports corresponding to the same band number can be very far away from each other. The direct consequence is that in the output image corresponding to a specific band, very different frequency modes can be mixed providing a “strange” representation. In order to avoid these issues, a solution consists to use the same filters for each rows (or columns). A simple way to build such solution is to consider a mean spectrum for the rows (or the columns), then perform the detection of the Fourier supports based on this mean spectrum and finally use the same filters for all rows (or columns). The complete algo- rithm corresponds to perform the following steps (we denote the input image):
- Take the 1D FFT of each rows of ; ; and compute the mean row spectrum magnitude: ,
- Take the 1D FFT of each columns of ; ; and compute the mean row spectrum magnitude: ,
- Perform the boundaries detection on and build the corresponding filter bank ,
Fig. 23. Time-frequency representation based on the EWT for signals (top) and (bottom).
- Perform the boundaries detection on and build the corresponding filter bank ,
- Filter along the rows with which provides output images,
- Filter each previous output image along the columns with , this provides at the end subband images. In order to illustrate this 2D version of the EWT, we per- form the previously described algorithm on the synthetic image given in left of Fig. 25. In order to see the adaptability of the 2D-EWT, this image is composed of objects and different 2D AM-FM components as shown in the Fourier domain (right of Fig. 25). We fix the maximum number of band per direction to two. The resulting subband images are given in Fig. 26 and the corresponding empirical tilling of the Fourier domain in Fig. 27. We can see that the algorithm is capable to isolate each AM-FM component in different Fourier domains as we expected.
VIII. C ONCLUSION In this paper, we proposed a novel approach to build wavelets adapted to represent the processed signal. The key idea is to
GILLES: EMPIRICAL WAVELET TRANSFORM 4009
Fig. 24. Time-frequency representation based on the EMD (on top) and the EWT (on bottom) for the seismic signal.
TABLE I N UMBER OF B ANDS A UTOMATICALLY D ETECTED FOR D IFFERENT FOR E ACH S IGNAL
Fig. 25. Test image on left and the logarithm of its Fourier transform on right.
Fig. 26. Output of the 2D EWT of the test image.
Fig. 27. Resulting partitioning of the 2D Fourier domain.
build a wavelet filter bank based on Fourier supports detected from the information contained in the processed signal spec- trum. We showed that it is possible to build a tight frame set of wavelets and then defined the Empirical Wavelet Transform (EWT) and its inverse. Like the classic wavelets, the proposed empirical wavelets correspond, in the temporal domain, to di- lated versions of a single mother wavelet. The new insight lies in the fact the corresponding dilation factors don’t follow a pre- scribed scheme but are detected empirically. Some experiments made on simulated and real signals showed that this approach is able to separate different AM-FM modes composing the input signal. A comparison with the Empirical Mode Decomposition (EMD) showed that the EWT gives a more consistent decom- position while, generally, the EMD exhibits too much modes, which are sometime really difficult to interpret. Another advan- tage of the EWT compared to the EMD is that we can adapt the classic wavelet formalism to understand it. We also proposed a 2D extension based on the tensor product idea to process im- ages. In this case, the 2D-EWT provides an adaptive tilling of the 2D Fourier plane. However, in some cases, the proposed method can fail (like any other wavelet approaches) compared to the EMD. For in-
