Decomposition of functions into pairs of intrinsic mode functions

The intrinsic mode functions (IMFs) arise as basic modes from the application of the empirical mode decomposition (EMD) to functions or signals. In this procedure, instantaneous frequencies are subsequently extracted from the IMFs by the simple application of the Hilbert transform, thereby providing a multiscale analysis of the signal's nonlinear phases. The beauty of this redundant representation method is in its simplicity and extraordinary effectiveness in many important and diverse settings. A fundamental issue in the field is to better understand these demonstrated qualities of the EMD procedures and the elementary modes they produce. For example, it is easily observed that when an EMD procedure is applied to the sum of two arbitrary IMFs, the original modes are rarely reproduced in the generated collection of IMFs. An interesting question from a representation point of view may be stated as follows: for any given sufficiently smooth function and fixed n≥2, when is it possible to represent the function as a sum of (at most) n intrinsic modes? A more interesting question is whether such a decomposition is possible when the extracted modes are constructed from a common formulation of the intrinsic properties of the function being analysed. We provide an answer to these questions for a relaxed version of IMFs, called weak IMFs, which has been shown to be characterized in terms of eigenfunctions of Sturm–Liouville operators. The objective of this study is to further extend that analogy to the relationship between sums of weak IMFs and coupled Sturm–Liouville systems. The construction of this decomposition also provides a guide to an alternate characterization of the instantaneous frequency and bandwidth.


Introduction
The empirical mode decomposition (EMD) method was developed by Huang et al. (1998) to decompose functions into a superposition of natural modes, each of which could be easily analysed for their instantaneous frequencies and bandwidths. These natural modes, which were termed intrinsic mode functions (IMFs), were generated at each scale, going from fine to coarse, by an iterative procedure to locally isolate the modal behaviour. Application of EMD to real signals f(t) has the purpose of representing these signals as sums of simpler modes j, i.e.
f ðt Þ Z X M jZ1 j j ðt Þ; ð1:1Þ where each j j (see Cohen 1995) comes with a specific polar representation of the form j j ðt Þ Z r j ðt Þ sin q j ðt Þ: ð1:2Þ These modes generalize the standard real and imaginary parts of Fourier components, in which it is required that each r j ðtÞZf ðj Þ be constant and q j (t)Zjt be linear, but are more suitable for the study of non-stationary phase. A richer analysis of signals is provided when the amplitudes are allowed to vary and the phases are permitted nonlinear, or non-stationary, behaviour as in the representation (1.2). We wish to emphasize that even if one is provided a decomposition of a signal f in the form (1.1), for each given j j the particular polar representation (1.2) that should be used is ambiguous with many possible selections of reasonable pairs of amplitudes and phases. The objective of EMD is to extract the decompositions from among such highly redundant representations so that the amplitudes r(t) and the corresponding phases q(t) at each scale are both physically and mathematically meaningful. In the case that the signal j is causal, the polar representation may be recovered by the application of an appropriate (cf. Sharpley & Vatchev 2006) Hilbert transform. In this case each mode j, after a phase shift of p/2, is represented as the real part of a complex signal J Jðt Þ Z rðt Þ exp iqðt Þ: ð1:3Þ Obviously, the choice of amplitude-phase (r, q) in the representation (1.3) is equivalent to the selection of an imaginary part f since rðt Þ Z ffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi ffi j 2 ðt Þ C f 2 ðt Þ q ; qðtÞ Z arctan fðt Þ jðt Þ ; ð1:4Þ once some care is taken to handle the branch cut. An alternative to the Hilbert transform method for extracting instantaneous frequencies is the class of quadrature methods (see Cohen 1995). Whichever method is to be used should be governed by the context of the phenomena under study and the analysing procedure should produce, for each signal j within the signal population, a properly chosen companion f for the imaginary part. This companion function f should be unambiguously defined and should properly encode information about the component signal, which is in this case the IMF. Therefore, the collection of instantaneous phases present at a given instant (i.e. tZt 0 ) for a signal f(t) is heavily dependent upon both the decomposition (1.1) and the selection of representations (1.2) for each monocomponent. A discussion of the application of Hilbert transforms and alternative quadrature methods to IMFs to obtain instantaneous frequencies are provided in Sharpley & Vatchev (2006).
A partial sum of IMFs in the representation (1.1) can be considered as an approximation of a function f by IMFs. The stoppage criteria used in applying EMDs are numerous, but they could be considered as control of approximation errors in various metrics.
According to Huang et al. (1999, p. 423), the EMD method was motivated 'from the simple assumption that any data consists of different simple intrinsic mode oscillations'. Three methods of estimating the time scales of f at which these oscillations occur have been proposed are -the time between successive zero-crossings; -the time between successive extrema; and -the time between successive curvature extrema.
The use of a particular time scale is application dependent. Following the development in Huang et al. (1999), we define a specific class of signals with special properties that make them very well suited for phase and wave band analysis.
Definition 1.1. A function j(t) is defined to be an intrinsic mode function, or more briefly an IMF, of a real variable t, if it satisfies the following two characteristic properties: (i) j has exactly one zero between any two consecutive local extrema and (ii) j has zero 'local mean'.
A function that is required to satisfy only condition (i) will be called a weak IMF and W will be used to designate all such functions. Note that j2W precisely when the number of zeros and the number of extrema of j on I differ at most by 1.
In general, the local mean in condition (ii) in the EMD procedure is typically the pointwise average of the 'upper' and 'lower' envelopes as determined by the spline fits of knots consisting of, respectively, the local maxima and the local minima of j.
The concept of decomposing functions into oscillating modes is the essence of the Sturm-Liouville theory. In Sharpley & Vatchev (2006) we related the concept of monocomponents from phase analysis to the Sturm-Liouville theory by characterizing the weak IMFs as solutions of self-adjoint differential equations.
Theorem 1.2. Let j2C 2 (I ) be a weak IMF with simple zeros and extrema, then there exist positive continuously differentiable functions P and Q such that j is the solution of the initial-value problem for some t2I.
From this result, it follows that an instantaneous frequency and bandwidth pair ðq 0 ; r 0 =rÞ can be defined implicitly through the Prüfer substitution given by Pf 0 ðt Þ drðt Þ cos qðt Þ; f ðt Þ drðt Þ sin qðt Þ: ð1:6Þ The goal of this paper is to further extend the relationship between weak IMFs and eigenfunctions of self-adjoint differential operators by the representation of very general functions by superposition of weak IMFs from a common framework. Therefore, a natural question arises which must be answered if an EMD procedure, such as Hilbert-Huang transform (Huang et al. 1999), is to be realized in terms of some nonlinear minimization procedure in the context of structured dictionaries (e.g. the expository article of DeVore 1998).

The question then becomes
Under what conditions can functions be decomposed into a superposition of at most n IMFs that are inherently extracted from the function?
Here, n is a fixed number that is set in advance. We prove that any smooth function can be decomposed into two (or fewer) weak IMFs (see theorem 3.1). This is not surprising since the two IMFs are unrelated to f other than the one which should be a sinusoidal with amplitude larger than the uniform norm of f and frequency higher than any local frequency of f. In this case, there are no requirements coupling the two corresponding Sturm-Liouville operators for the j j . However, with more care, weak IMF pairs j 1 , j 2 can be constructed from f so that they are solutions corresponding to different eigenvalues of a Sturm-Liouville operator (SLO; with coefficients P, Q) and f Zj 1 Cj 2 . This is the content of theorem 3.2. The analogy between a decomposition such as EMD and coupled Sturm-Liouville systems with the same primary coefficients P and Q is provided in theorem 2.1. We prove that for any sufficiently smooth function f, there exists a coupled linear system of two differential equations, such that f represents one of the components of the solution. These coupled equations can be formulated in terms of a mechanical system of a pair of variable masses and springs and f then represents the displacement from equilibrium of the terminal mass of the system. This formulation gives an indication of a possible deeper relationship between decompositions of signals into superpositions of IMFs and concepts from system identification.
We begin by defining a general SLO that acts on sufficiently smooth functions f and is defined by Lf dK 1 Q ðPf 0 Þ 0 ; ð1:7Þ where P and Q are positive continuous functions on some interval IZ [a,b]. The spectrum of L consists of all u for which the equation LjZuj has a solution on I that satisfies certain boundary conditions. In a more general situation the Sturm-Liouville (SL) system is defined as the solutions of the boundary-value problem (BVP) ðPj 0 Þ 0 C ðuQKrÞj Z 0; jðaÞ cos a C PðaÞj 0 ðaÞ sin a Z 0; for positive functions P and Q and real constants u, a and b. It is well known (see Birkhoff & Rota 1989) that (1.8) has a solution for a discrete set of eigenvalues, such that 0!u 1 !u 2 !/ and lim n/N u n ZN. The corresponding eigenfunctions j u n form an orthogonal basis in the weighted Hilbert space L 2 ðI ÞðQÞ dff : Ð I jf ðt Þj 2 Qðt Þ dt!Ng with inner product Ð I fgQ dt, and hence any function f2L 2 (I )(Q) has a decomposition with convergence in L 2 (I )(Q). For a fixed SL system, there exists a unique representation of any f2L 2 in the form of (1.9) but, in general, requiring infinitely many terms. This should be contrasted to theorem 3.2.

Coupled differential systems for smooth functions
In this section, we prove the main result of this paper for representing a function as a component of a solution of a linear system of differential equations. In the proofs of the theorems, we require two technical lemmas. In order not to disturb the main flow of ideas, we refer the reader to appendix A for their proofs.
Throughout this section assume f is defined on a finite interval IZ[a,b] and has only simple inflections and set XZ fx 2 I : f 00 ðxÞZ 0g has M!N elements.
Theorem 2.1. For a function f2C 2 (I ) with simple extrema and for any positive numbers k 1 and k 2 , there exist continuous and positive functions P and Q, depending upon f, and twice differentiable function x such that f and x satisfy the following system of differential equations: First assume that a 0 2J 0 . The existence of P and Q is shown by inductively constructing them on I 0 and J 0 and repeating the same procedure on I j , J j for jO0.
For fixed k 1 , k 2 positive, set Since f, f 0 and f 00 are bounded on I, then for any given real number u, it follows that C 1 ZjðbKaÞujCðjl 1 j=jl 2 jÞðjf 0 ðaÞðbKaÞjCjf ðaÞjCkf k N Þ and C 2 ZððbKaÞ 2 Þ= jl 2 jkff 00 k N are finite constants.
(We note that these are the primary constants and parameters in lemma A.2 of the appendix. Although somewhat redundant, introducing these here aids in motivating the proof of alternating the successive application of the two lemmas.) On the interval J 0 the function f does not have an inflection, and hence jf 00 jO hO 0. By applying lemma A.2 with T ZKsignðf 00 ða 0 ÞÞT and f 0 ða 0 ÞZ u, we obtain a twice differentiable solution F that does not change sign on J 0 . Furthermore, since f 00 does not change its sign on J 0 and owing to the choice of the initial condition for F(a 0 ), it follows that the function Q 0 dKðP 0 f 00 =FÞ is positive and continuous on J 0 for any fixed function P 0 O0. Fix P 0 a positive constant and define xZfK(1/k 1 )F, then it follows that P 0 f 00 ZKk 1 Q 0 ðf KxÞ: ð2:4Þ Twice differentiating the integral equation in lemma A.2, we get F 00 ZKðl 1 =l 2 Þf 00 Kð1=l 2 Þðff 00 =FÞ and after substituting in x 00 Z f 00 Kð1=k 1 ÞF 00 , we obtain the following differential equation for x: By combining (2.4) and (2.5), we have constructed the desired system (2.1) for J 0 . Next we construct the functions P and Q on the interval I 0 , but here the function f has an inflection point and lemma A.2 cannot be applied. Instead, on this interval this will be the role of lemma A.1. Note that the system (2.1) is equivalent to the system where pd(P 0 /P) and qd(Q/P). Let q be the linear function connecting the points ðb 0 ; ðQ 0 ðb 0 Þ=P 0 ÞÞ and ða 1 ; ðjf 00 ða 1 Þj=T ÞÞ. Since f 0 R hO 0 on I 0 , then the second equation can be solved for p and upon substituting in the first equation of (2.6), we obtain the differential equation in lemma A.2 for x. The initial conditions naturally come from the already determined values y 0 Z xðb 0 Þ; y 1 Z x 0 ðb 0 Þ. On J 0 we have that FZKð1=qÞf 00 and jFjR ðT=2Þ including the right endpoint b 0 , and hence the estimate jqðb 0 ÞjZ ðjf 00 ðb 0 Þj=jFðb 0 ÞjÞ% ð2jf 00 ðb 0 Þj=TÞ% 2kf 00 k N holds true. From the definition of q it follows that jqða 1 Þj% ð2kf 00 k N =TÞ and since q is linear on I 0 we have that the last estimate holds on over all of I 0 . All the conditions in lemma A.1 are met, and hence there exists a function x that is a solution to the differential equation (A 2). Substituting x in (2.7), we get a continuous function p on J 0 gI 0 . The functions P 1 and Q 1 are determined from the system and have solutions P 1 ðt ÞZ P 0 expð Ð t b 0 pðvÞ dvÞ; Q 1 ðtÞZ P 0 qðtÞ expð Ð t b 0 pðvÞ dvÞ, where P 0 ZP 1 (b 0 ) and Q 0 (b 0 )ZQ 1 (b 0 ). We have now constructed the continuous functions P and Q on the interval J 0 gI 0 . Inductively, we can continue the construction on J 1 with the same constants l 1 and l 2 as on J 0 and initial conditions fða 1 ÞZKsign ðf 00 ða 1 ÞÞð1=qða 1 ÞÞ; f 0 ða 1 ÞZ u, u arbitrary. Note that jfða 1 ÞjZ jTj and hence lemma A.1 provides continuous positive P 2 and Q 2 on J 1 . Repeating the same procedure on I j and J j , jR1, we construct continuous functions P, Q and x, which satisfy the system (2.1). This completes the proof if the partition of IZ[a,b] begins with J 0 .
In the case a 0 2I 0 , we then start by instead using lemma A.1 with appropriate choice for the initial conditions and the function q. & In the previous theorem, we showed the existence of continuous positive coefficients P and Q assuming that the function f has a continuous second derivative. If f has fourth-order continuous derivatives, i.e. f2C 4 (I ), the next result shows that we can modify our construction to guarantee in addition that the coefficient P is continuously differentiable, so that the equations (2.1) hold in the classical sense.
Theorem 2.2. Let f be a real-valued function in C 4 (I ). If f has only simple extrema, then for any two positive numbers k 1 and k 2 there exist positive functions P and Q such that P, P 0 and Q are continuous on I, and f is a solution to (2.1).
Proof. Following the proof of theorem 2.1, it is clear that the only points where P 0 could have breaks are the points a j , jZ1, 2, ., M. We show that around these points the construction of P can be adjusted in a way that provides continuous P 0 .
Let k 1 and k 2 be fixed and l 1 , l 2 andT be as in theorem 2.1, and where mZ(P/Q) and sZ(P 0 /Q). The system (2.1) is equivalent to the fourthorder equation where L 2 fZL(Lf ). After substituting (2.8) into (2.9) and isolating s 00 in the resulting equation we obtain, for the given f, the second-order linear equation for s s 00 Z Fðt; s; s 0 Þ; ð2:10Þ where In the proof of theorem 2.1, the functions P and Q were constructed continuous and bounded on I, and hence the corresponding p and q are continuous and bounded on I. Clearly P 0 Z(s/m)P, and hence P, P 0 , and Q are continuous if and only if s and m are continuous. On the other hand, qZ(1/m) and pZ(s/m). By using the existence theorem for a solution of BVP (appendix lemma A.1), we alter p and q in one-sided neighbourhoods of the points a j in a way that p and q are continuous on the entire interval I. For a fixed a j , there exists a neighbourhood in which the function f does not have extremal points and inflections, and hence the function (f 00 /f 0 ) has a constant sign on that neighbourhood, say eZ sign t!a j ðf 00 ðtÞ=f 0 ðt ÞÞ. We modify q, obtained from theorem 2.1, on the interval I j preserving all necessary conditions required for the proof of theorem 2.1. Let eZ ðT=jf 00 ða j ÞjÞ and rZ(e/4R), where RO0, then on the interval I j,r Z[a j Kr,a j ] the straight line mðxÞZK2eRðx K a j ÞC e is strictly positive and ð1=mÞ% ð2 max t2I jf 00 j=TÞ. If necessary, R might be increased in order that I j,r be contained entirely in I j . By using the previously constructed P and Q from theorem 2.1, the new q is modified only on the interval I j . Namely, on I jK1 \I j,r we define q as the linear function that connects the points ðb j ; ðQ j ðb j Þ=P j ÞÞ and ða j Kr; 1=ðmða j KrÞÞÞ and on I j,r we define q(t)Z1/m(t). Similar to the construction in theorem 2.1, we construct the function p on I jK1 \I j,r and let s 0 Z sða j KrÞZ pða j KrÞmða j KrÞ. By using the already defined p on J j , let s r Z sða j ÞZ pða j Þmða j ÞZ pða j Þe. We show that on I j,r , the BVP (2.10) with boundary conditions sða j KrÞZ s 0 ; sða j ÞZ s r has a solution sðt ÞZ pðt Þmðt Þ, and hence the newly defined q and p are continuous on I jK1 g J j .
Without loss of generality, we can consider I r Z I j; r Z ½0; r. Let RO0 be large enough such that js 0 j; js r j% R. Since Fðt; s; s 0 Þ is linear with respect to s 0 on E(r, R) (see appendix A 1), then for a fixed R there exists constants g and CO0 depending only on f and R such that gR!1 and jF j% gðs 0 Þ 2 C C on E(r, R). Next we show that the conditions for Fðt;GR; 0Þ hold in limit for R/N. From the choice of r, it follows that jmj% ð3=2Þe independent of R on I r . Let m be a constant with a value either 1 or K1. Since f and its derivatives do not depend on R and m 00 Z0 on I r , it follows that sign lim holds uniformly on I r . Hence we can pick large enough R for which Fðt; R; 0ÞR 0 and Fðt;KR; 0Þ% 0. With that choice of R, all the necessary conditions in lemma A.1 are met and the solution to (2.10) with sða j KrÞZ s 0 ; sða j ÞZ s r provides a solution p. By repeating the same considerations for all of the a j 's, we can construct positive coefficients P and Q for the system (2.1) such that P, P 0 and Q are continuous functions on the entire I. Remark 2.3. In the previous discussion, we proved that for an appropriate function f there exists companion function x such that the pair (f, x) is a solution of a system of type (2.1). The function x was constructed first on neighbourhoods around the inflections of f and then, starting from the beginning of the interval I, the pieces were smoothly connected. In order to construct a particular function x, we need to specify the initial values x(a) and x 0 (a). On the other hand, we can start from any interior pointt 2 I with initial values xðt Þ and x 0 ðt Þ and construct the function on the left and on the right fromt .
Furthermore, for functions that satisfy the conditions in theorem 2.2, we considered a method to smoothly connect two disjoint segments of x. Summing up, we can construct the function x to satisfy boundary-initial conditions (BVP) by initializing the 'left' branch of x at a and the 'right' branch of x at b, followed by connecting smoothly the two branches on an appropriate subinterval of I.
In the next section, the relation between the EMD and the representation of f as a solution of system of differential equations is discussed.

Decomposition into pairs of weak IMFs
In the current section, we prove that any function f2C 2 (I ) with simple zeros and extrema can be decomposed into a sum of two or fewer weak IMFs. The decomposition is constructive and it is not obtained by applying EMD. Two types of decompositions are considered. The first one is a more general result about the existence of only two weak IMFs. In general, the two weak IMFs are not strongly related to physical properties of the function. The second method for decomposition can be associated with the displacement from equilibrium of one of the two masses in a simple mechanical system with time-dependent masses and forces. Details are given later in this section. In either case, the decompositions are non-unique. We denote the interval by IZ[a,b].
Theorem 3.1. Let f2C 1 (I ), then there exist two weak IMFs j 1 and j 2 , such that fZj 1 Cj 2 on I.
Proof. Let j 1 be any weak IMF such that j 0 1 is also a weak IMF and the extrema of j 1 do not coincide with the extrema of f. For example, j 1 (t) could be sin(mtCz) for appropriate choices of m and z. Define g a Z f Kaj 1 for a2R. Since f is bounded on I, then for any point t2I we have that lim a/N sign ðg a ðt ÞÞZ sign ðj 1 ðt ÞÞ and lim a/N sign ðg 0 a ðtÞÞZ sign ðj 0 1 ðtÞÞ. Since f, f 0 , j 1 and j 0 1 are continuous on I, there exists large enough a 1 for which j 2 Z g a 1 is a weak IMF and f Z a 1 j 1 C j 2 on I. Since a 1 j 1 is a weak IMF the proof is complete. & From the proof, it is clear that j 1 and j 2 may not be related in general to the local properties of the function and further analysis based on that representation will not reveal any useful local characteristics of f. Another type of decomposition into weak IMFs, which does admit physical interpretation for a sufficiently smooth function f, can be obtained from theorem 2.1. That is the content of the next result.
Theorem 3.2. Let f2C 2 (I ), then there exist weak IMFs j 1 and j 2 , such that f Z a 1 j 1 C a 2 j 2 on I, and j 1 and j 2 are solutions of the self-adjoint differential equations ðPj 0 j Þ 0 ZKu j Q j j , for some positive continuous functions P and Q and positive u j ; j Z 1; 2.
Proof. Fix k 1 , k 2 O0. By applying theorem 2.1 to f, it follows that there exist functions P, Q and x such that the vector is a solution of the system of linear differential equations ðPEX 0 Þ 0 Z QKX, where Direct calculations show that the numbers Ku 1 and Ku 2 , where u 1 Z ð2k 1 C k 2 K ffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi ffi 4k 2 1 C k 2 2 p Þ=2O 0 and u 2 Z ð2k 1 C k 2 C ffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi ffi 4k 2 1 C k 2 2 p Þ=2O 0 are the two eigenvalues of the matrix K, and hence there exists an invertible matrix of real numbers A such that The functions j 1 and j 2 are two weak IMFs such that f Z a 1 j 1 C a 2 j 2 for some real a's. Indeed, we have From (3.1) it follows that the functions j 1 and j 2 are solutions to the differential equations ðPj 0 Þ 0 ZKuQj for uZu 1 and uZu 2 , respectively. Hence they are weak IMFs generated by one and the same self-adjoint operator. Adding the initial conditions from f, and xðaÞZ x a ; x 0 ðaÞZ x 0 a , where x a ; x 0 a are determined from the proof of theorem 2.1, we can determine constants a 1 and a 2 such that f Z a 1 j 1 C a 2 j 2 . & In general, the functions j 1 and j 2 do not satisfy the same boundary conditions, and hence they are not part of one and the same Sturm-Liouville system. In the case f2C 4 (I ), we can use remark 2.3 to construct j 1 and j 2 as two eigenfunctions of one and the same Sturm-Liouville system.
Corollary 3.1. Let f2C 4 (I ) and have simple zeros, extrema and inflections. If u 1 and u 2 are as in theorem 3.2, then there exists a Sturm-Liouville system of type (1.8), such that f is a linear combination of two, or fewer, of its eigenfunctions.
Proof. We modify P and Q in such a way that the functions j 1 and j 2 constructed in theorem 3.2 satisfy the same boundary condition of type (1.8). From the proof of theorem 3.2, it follows that where A K1 is a non-singular matrix, and hence j 1 and j 2 satisfy the same boundary conditions if and only if f and x satisfy the same boundary conditions.
Since the function f is given, the parameters in the boundary condition should be chosen appropriately. The function P is coupled with the function x but a and b could be arbitrary.
Let x l denote a function x that is constructed in theorem 2.2 by starting from the end point a. We show that x l can be constructed to satisfy (3.2). Depending on the values of f and f 0 at a, we consider the following three possible cases.
-Case (i). The function f does not have both a zero and an extrema at a. Then, (3.2) is equivalent to f ðaÞKf ðaÞ=f 0 ðaÞf 0 ðaÞZ 0, and so x must satisfy x l ðaÞKf 0 ðaÞ=f ðaÞx 0 l ðaÞZ 0. In the case a2J 0 , from theorem 2.1, it follows that x l constructed as x l ðtÞ Z f ðt ÞKð1=k 1 ÞFðtÞ satisfies x l ðaÞKðf ðaÞ=f 0 ðaÞÞx 0 l ðaÞZ 0, if F is the solution of the initial-value integral equation (A 11) with initial condition F(a)ZT and F 0 ðaÞZ ðf 0 ðaÞ=f ðaÞÞT.
In the case a2I 0 , from theorem 2.1, it follows that the function x l , constructed as the solution of the IVP (A 2) with initial conditions x l (a)Zf (a) and x 0 l Zf 0 (a), satisfies (3.2). -Case (ii ). The function f does not have zero at a, but f 0 (a)Z0. By choosing aZp/2 in (3.2), the initial condition for f becomes f 0 (a)Z0 and the same considerations as in case (i) lead to a construction of x l with x 0 l Z0 and any choice of x l (a).
-Case (iii). The function f does not have extrema at a, but f (a)Z0. In that case, we construct P with zero at a and positive elsewhere. Choosing aZp/2 in (3.2), we have a singular Sturm-Liouville system. Since f has simple zeros and extrema, it follows that f 0 (a)s0 and a2I 0 . Furthermore, by repeating the same construction of p and q as in theorem 2.2 with boundary values s(a)Z0 and s(r)Zv, for an appropriate r and v, we obtain a continuously differentiable non-negative function P with its only zero value at the point a.
The corresponding Q is a continuous positive function.
In a similar way, starting from b, a function x r can be constructed that satisfies the boundary conditions (3.3). Finally, x l and x r can be connected smoothly on an appropriate interval as in theorem 2.2 to construct a function x that satisfies (3.1) and (3.2). The resulting weak IMFs j 1 and j 2 are eigenfunctions of a Sturm-Liouville system of type (1.8) with corresponding eigenvalues u 1 and u 2 . One immediate consequence is that j 1 and j 2 are orthogonal on I with a weight function Q. & The main purpose of the EMD method is to decompose a function into a sequence of IMFs in order to extract important physical characteristics from each component by the application of the Hilbert transform. Next we discuss possible physical interpretations of the decompositions considered in theorems 3.1 and 3.2.
Let SL(P, Q) denote the set of all the eigenfunctions of Sturm-Liouville operators defined in (1.7). It is clear that SL(P, Q) contains the eigenfunctions of any Sturm-Liouville system. For any fixed boundary conditions, the system of corresponding eigenfunctions is complete in L 2 (I, Q) and it follows that the set of functions SL(P, Q) is complete and redundant in L 2 (I, Q). Furthermore, since the weight function Q is strictly positive, continuous on the finite interval I, then the space L 2 (I, Q) is norm equivalent to the space L 2 (I ). Varying P and Q over all continuous positive functions on I, we observe from theorem 1.2 that W Z fj : j 2 SLðP; QÞ; P; Q positive continuous functions on I g.
Self-adjoint differential equations describe many physical processes and have embedded physical characteristics. A typical example is the periodic BVP j 00 C u 2 j Z 0; jð0Þ Z jð2pÞ; j 0 ð0Þ Z j 0 ð2pÞ; that defines the trigonometric system fcos nt; sin ntg N nZ0 on the interval [0,2p). This differential equation models the frictionless displacement from equilibrium of a unit mass attached to a spring with a spring constant u 2 . If the 'periodic' boundary condition is replaced with an initial condition, then the SL operator LfZKf 00 generates trigonometric functions cos ut, sin ut that still obey the physical relation with mass-spring system, but is a redundant system. One method to decompose a function into a linear combination of elements from a redundant system is by using 'greedy'-type algorithms (for reference see DeVore 1998 andTemlyakov 1999). Generally, the resulting representation consists of infinitely many terms.
In analogy with the trigonometric system, the characteristic equation for a weak IMF ðPj 0 Þ 0 ZKQj can be interpreted as the physical model according to which a variable mass P(t) is attached to a spring with a variable spring constant Q(t) and vibrates frictionlessly around its equilibrium position j(0). The solution j(t) is the displacement of the mass at the instant t. An equivalent form of the characteristic differential equation is j 00 C pj 0 C qj Z 0; ð3:4Þ where qO0 and p is an arbitrary function. In that form the mass-spring interpretation could be that j(t) is the displacement of a unit mass attached to a spring with a variable spring constant q if the motion is subjected to a frictional force Kpj 0 . In both interpretations, ffiffiffiffi Q p or ffiffi ffi q p is considered the frequency of the motion. In the second case, p can be considered as the instantaneous bandwidth of j. The instantaneous quantities of j can be defined by using the Prüfer substitution (1.6).
Summarizing the results of this section, theorems 3.1 and 3.2 not only show the existence of the decomposition of a function f into finitely many weak IMFs, but they also provide an additional analytical resource to the Hilbert transform method in order to define instantaneous frequency and bandwidth by using the physical interpretation from a model system of differential equations.
Further, select positive c m such that J is then set to the interval centred at x with length d dð1=c 2 m Þ. If necessary, c m can be increased to ensure d!1 and J3J 0 . We show that J is the interval whose existence is stated in the lemma.
Write the second-order equation (A 2) as its equivalent first-order system in integral form, where G(t,x,z) is as in expression (A 1). A standard technique to show the existence of solutions of differential equations is the method of successive approximations. In particular, we may, without loss of generality, assume that y 0 Zy 1 Z0 and let x 0 Z0, z 0 Z0, then inductively construct the sequence pairs x nC1 ðtÞ Z y 0 C ð t a z n ðtÞ dt; Gðt; x n ðtÞ; z n ðtÞÞ dt: If we show jx n j%c m and jz n j%c m for nZ0,1, 2,., we can use a standard argument to then establish that the limit functions x and z exist and the convergence is uniform. Using the inductive assumption jx n j%c m and jz n j%c m , we show that jx nC1 j%c m and jz nC1 j%c m . For x nC1 , we have jx nC1 j% jz n jd% jz n j% c m , since d!1. For z nC1 , using the bound jqj% 2kf 00 k N , we have jz nC1 j% dðm 0 C m 1 jx n j C m 2 jz n j C m 3 jx n jjz n jÞ: The polynomial function bðx; zÞZ m 0 C m 1 x C m 2 z C m 3 xz has only one local extremum and the value at that extremum is mZ jm 0 K3ðm 1 m 2 =m 3 Þj! c 3 m , from the choice of c m . The global maxima of b(x,z) on the square jxj% c m ; jzj% c m are either the value at the extremum or a value from the boundary. It is easy to see that the maximum value on the boundary is bðc m ; c m ÞZ m 0 C ðm 1 C m 2 Þc m C m 3 c 2 m ! c 3 m , from the choice of c m . Hence jz nC1 j% d max jbðjxj; jzjÞj% ð1=c 2 m Þc 3 m Z c m and the limit functions x Z lim n/N x n , and z Z lim n/N z n is the solution pair of (A 9). Hence x is a solution of the IVP (A 2) and (A 3) on the interval [a,aCd] that contains the point x. &

(b ) Integral equation
One additional technical lemma is required for the proof of theorem 2.1. The result establishes the existence of positive solutions of an integral equation on intervals that avoid inflections of f and will be used on subintervals complementary to those handled by lemma A.1.
Lemma A.2. Let f 2 C 2 ðI Þ with simple zeros and extrema. Lett 2 I . Then for any non-zero l 1 and l 2 and any u, there exists a positive real number T 0 , depending only on I and f, such that for each T, jT jR T 0 , the integral equation, has a solution F, such that jFjO ðjTj=2Þ on I, with fðt ÞZ T, f 0 ðt ÞZ u.
Proof. We prove the lemma in the casetZ a, where IZ[a,b]. The casẽ t 2 ða; b is similar. Let c 1 Z jðbKaÞujC ðjl 1 j=jl 2 jÞðjf 0 ðaÞðbKaÞjC jf ðaÞjC kf k N Þ and c 2 Z ðjI j 2 =jl 2 jÞkff 00 k N . Since f, f 0 and f 00 are bounded on I, it follows that c 1 and c 2 are finite constants. Then pick T 0 positive such that jT jRT 0 implies jT jO c 1 C ffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi c 2 1 C 4c 2 q : ðA 12Þ Substituting fCT for f in (A 11) and applying the method of successive approximations, we consider a sequence of functions f n that satisfy the initial conditions By construction jf 0 jZ 0! ðjT j=2Þ. Assume jf n j! ðjT j=2Þ, then jf n C TjR ðjT j=2Þ. Estimating the r.h.s. of (A 13), it follows from the selection of T in (A 12) that jf nC1 j% c 1 C 2c 2 jTj ! jTj 2 : ðA 14Þ Hence the limit function FZ T C lim n/N f n satisfies (A 11) with Fðt ÞZ T, F 0 ðt ÞZ u and jFjR ðjTj=2Þ on I. & Finally, we include a necessary condition for the existence of a solution of a second-order BVP, the details of which may be found in Hartman (1964). Proposition A.3. Suppose F(t, x, x) is continuous on the tube domain Eðr; RÞZ fðt; x; xÞ : 0% t% r; jxj% R; x realg and that F satisfies the condition that both F(t, R, 0) is non-negative and F(t, KR, 0) is non-positive for all 0%t%r. Assume further that there exist positive numbers C, g with g!1/R so that jFj% gx 2 C C, then there exists a solution to the BVP x 00 Z Fðt; x; x 0 Þ xð0Þ Z x 0 ; xðrÞ Z x r ; jx 0 j; jx r j% R: ) ðA 15Þ