ISRN Signal ProcessingVolume 2013 (2013), Article ID 605035, 18 pageshttp://dx.doi.org/10.1155/2013/605035

Research Article

About a Partial Differential Equation-Based Interpolator for Signal Envelope Computing: Existence Results and Applications

Oumar Niang,1,2,3 Abdoulaye Thioune,2,4 Éric Deléchelle,² Mary Teuw Niane,³ and Jacques Lemoine²

¹Département Informatique et Télécommunications, Ecole Polytechnique de Thiès (EPT), Thiès BP A10, Senegal²Laboratoire Images, Signaux et Systèmes Intelligents (LISSI-E.A.3956), Université Paris-Est Créteil Val-de-Marne, Créteil, France³Laboratoire d’Analyse Numérique et d’Informatique (LANI), Université Gaston Berger (UGB), Saint-Louis BP 234, Senegal⁴Département de Mathématique et Informatique, Faculté des Sciences et Technique, Université Cheikh Anta Diop de Dakar, Dakar BP 5005, Senegal

Received 6 November 2012; Accepted 3 December 2012

Academic Editors: H. Hu and S. Li

Copyright © 2013 Oumar Niang et al. This is an open access article distributed under the Creative Commons Attribution License, which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited.

Abstract.

This paper models and solves the mathematical problem of interpolating characteristic points of signals by a partial differential Equation-(PDE-) based approach. The existence and uniqueness results are established in an appropriate space whose regularity is similar to cubic spline one. We show how this space is suitable for the empirical mode decomposition (EMD) sifting process. Numerical schemes and computing applications are also presented for signal envelopes calculation. The test results show the usefulness of the new PDE interpolator in some pathological cases like input class functions that are not so regular as in the cubic splines case. Some image filtering tests strengthen the demonstration of PDE interpolator performance.

1. Introduction

Interpolators are widely used in signal and processing or data analysis. In particular for the empirical mode decomposition (EMD) algorithm [1, 2], the iterative estimation of the signal trend is based on the computing of the envelopes obtained by the cubic spline interpolation of local extrema. The spline interpolation has been recognized as being very effective for EMD. But for signals that have no local extremum, the cubic spline interpolation fails. We proposed a PDE-based model which overcomes this limit of classical EMD implementation, in the computing of the envelopes for signals that have no local extremum [2–4]. Recently, the Spectral Intrinsic Decomposition (SID) method [5], based on the spectral decomposition of the PDE interpolator, provides a new application of our model. This PDE interpolator contribute, to the mathematical modeling of the EMD and has provided various applications in signal and image processing [4, 6]. In this paper, we describe the mathematical modeling of the new PDE interpolator by variational methods. The resolution of the variational problem leads to existence and uniqueness results in appropriate spaces. The paper is organized as follows. In Sections 2 and 2.1 recalls some PDE models in signal and image processing, and in Section 2.2 some mathematical preliminary notions are set out. In Section 3, the mathematical modeling is exposed. In Sections 4 and 5, the resolution of the variational problem is dealt. Subsequently, numerical implementation and applications are presented in Section 6. At last, we finish by conclusions.

2. Some Preliminaries

2.1. PDE Models in Signal and Image Processing and EMD Principle

2.1.1. Some Diffusion Equations

This part consists of a brief and nonexhaustive presentation of classical nonlinear diffusion filters for 1D and 2D signal processing, with more focus on the 2D case. In this purpose, we can recall the model for nonlinear diffusion in image filtering proposed by Catte et al. [7]. This filter is a modified version of the well-know Perona and Malik model [8]. The basic equation that governs nonlinear diffusion filtering is with , and where is a filtered version of the original image as the initial condition, with reflecting boundary conditions. In (1), is the conductivity (or diffusivity) function, which is dependent (in space and time) on the image gradient magnitude. Several forms of diffusivity were introduced in the original paper of Perona and Malik [8]. All forms of diffusivity are chosen to be a monotonically decreasing function of the signal gradient. Possible expressions for conductivity are Parameter is a threshold one, which influences the anisotropic smoothing process. The nonlinear equation (1) acts as a forward parabolic equation smoothing regions while preserving edges. Other methods based on high-order PDE are provided for image restoration like in [9–11]. Efficient numerical schemes were introduced in [12] based on additive operator splitting (AOS) scheme, or based on alternating direction implicit (ADI) scheme. See [12–15] for a review and extensions of these methods. Unlike these methods of high-order PDE that are specially developed for denoising, our model was constructed to interpolate the characteristic points of a signal.

The major problem of nonlinear diffusion-based process is that it is generally difficult to correctly separate the high frequency components from the low frequency ones. In case of denoising applications, the objective of this process is to use the diffusivity function as a guide to retain useful data and suppress noise.

Numerous authors have proposed fourth-order PDEs for image smoothing and denoising with the hope that these methods would perform better than their second-order analogues [16–18]. Indeed, there are good reasons to consider fourth-order equations. The first reason that can be retained is the fact that fourth-order linear diffusion damps oscillations at high frequencies (e.g., noise) much faster than second-order diffusion. On the other hand, the theory of fourth order nonlinear PDEs is far less developed than the second order one. Also such equations often do not satisfy a maximum principle or comparison principle, and implementation of the equations could thus introduce artificial singularities or other undesirable behavior. In recent studies, Tumblin [19], Tumblin and Turk [16], and Wei [17] proposed the following form: where as in (2) and is some measurement of . In [16], (3) is called a “low curvature image simplifier” (LCIS), and a good choice for is defined as to enforce isotropic diffusion [19]. These PDE tools for digital signal and image processing make more reachable the 2D extension of the 1D PDE-based method for characteristic points interpolation presented in this paper.

2.1.2. The Empirical Mode Decomposition Principle

The EMD [1] method decomposes iteratively a signal into amplitude modulation-frequency modulation (AM-FM) type components called intrinsic mode functions (IMF). The underlying principle of this decomposition is to locally identify in the signal, the most rapid oscillations defined as the waveform interpolating interwoven local maxima and minima. To do this, local extrema points are interpolated with a cubic spline, to yield the upper and lower envelopes. The mean envelope (half sum of upper and lower envelopes) is then subtracted from the initial signal, and the same interpolation scheme is reiterated on the remainder. The so-called sifting process stops when the mean envelope is reasonably zero everywhere, and the resulting signal is called the first IMF. The higher-order IMFs are iteratively extracted applying the same procedure to the initial signal after the previous IMFs have been removed.

2.2. Some Useful Mathematical Concepts

Throughout the paper, denotes an open-bounded subset of ,    or  2. In the sequel, we will need the following definitions and results.

Definition 1. we define one the spaces and as follows

Definition 2. For a given function , if and denote, respectively, the right and left derivative of at ; the minmod derivatives of is define by A function being absolutely continuous admits right and left derivatives, then has obviously left and right derivatives, so that we can validate numerically computing of the diffusivity function defined in Section 3.

Definition 3. Let be a sequence of elements in a vectorial normed space ; it is said that converge, weakly [20] in , and noted by , if exists an element such that ,  , where denotes the set of continuous linear forms on .

Definition 4 (the Green formula). Let , then where denotes the normal derivative of on the boundary of .

Definition 5 (the Poincaré-Wirtinger inequality). Let be an open-bounded set, and let , then there exists a constant such that the norm of in and the norm of in are linked by the following inequality: where denotes the length of . The best constant in the Poincaré-Wirtinger inequality is , for example, the inverse of the first positive eigenvalue of the Laplace operator with homogenous Neumann boundaries conditions. In our case, , and is the diameter of .

Theorem 6   (regularity theorem). Let then there exist two strictly positive reals and such that firstly, and secondly,

The main existence and uniqueness result of the solution of our henceforth problem (29) is due to the application of the following Lions theorem [21].

Theorem 7 (the Lions theorem). Let and be two Hilbert spaces with . One considers a bilinear application on and an operator on . Under the following conditions: (1) and such that ,(2) such that . For , the problem such that and , has a unique solution. Furthermore, .

3. Mathematical Modeling of the New PDE-Interpolator

Recall that the PDE model aimed to contribute to the mathematical modeling of the EMD, as it is in [2]. In a first step, to be in line with the classical EMD, our goal is to model the upper or lower envelope as the asymptotic solution of a PDE system whose initial condition is the input signal that we want to interpolate local extrema (or more generally characteristic points). Initially, this envelope is obtained by cubic spline interpolation. Let be the open domain of an given finite energy signal . We construct an operator of appropriate domain as follows:

The need of asymptotic solution existence denoted by implies , for example, according to (13), we have . Moreover, the requirement of the same regularity between the solution and a cubic spline leads to . Thus, in the first analysis, we can choose . To leave invariant features points during the diffusion process, simply multiply in the expression the term by function depending on spatial variable and vanishing at characteristic points. That gives Another reminiscent form of long-range diffusion [22] is the following: Function can be interpreted as a diffusivity function whose role is to control the diffusion process. We take it necessarily positive to have a direct diffusion and the existence of solution. Other forms of operators are possible [3] we present some ones as follows: for (14) or for (15), or accomplete diffusion involving a flow where is the tension factor which controls the regularity of the solution.

Both forms allow more freedom on the regularity of the solution . To access the mathematical properties of the solution, we chose the second form in (17) with and the corresponding equation (15).

At the local scale, between two consecutive characteristic points—two local maxima, for example—the diffusion induces a smoothing phenomenon that deletes the local minimum. A simple form for to calculate the envelope is given by a positive piecewise function lower than that is constant between two characteristic points of and zeroed only at these points. Characteristic points are often being defined by their values and the signs of first, second, and third local derivatives of . We characterize the function as depending on , and .

For the purpose of existence and regularity of the solution, we are led to work with the regularized version of : where is a regularization coefficient.

In the following, we define for extrema detection, for turning points detection, and for maximum and minimum curvature points detection: (+) for maximum curvature points and (−) for minimum curvature points. All these functions are of the form   with at characteristic points.

So, we have if . This property formally allows us to cancel and of the form (15) at the points where is null. The factor is a normalization term and verifies . A more general diffusivity function could be envisaged.

Let us consider the following differential operators and ( for upper and for inflexion or curvature): or

If , as it is the case in the classical implementation EMD, there is no problem for interpolation by cubic spline method. But it is indeed a restriction on the regularity of which is due to the choice of the interpolation technique that requires a good detection of local extrema for the envelopes calculation. The basic EMD principle must apply for input functions that are not regular as functions in , for example, for signal , as it is the case in reality [23]. But numerically, without even a regularization function , if the derivatives are taken in the sense of minmod called flux limiter or slope, the function still has a meaning.

3.1. Interpretation of the Diffusivity Action in the Diffusion Process

Between two consecutive maxima and , the function is piecewise constant and can be written as follows:

We have a diffusion in the sense that unhook to the maximum of the curve . The reason is that the diffusion is more pronounced in local minimum because the smoothing effect tends to regularize curves. Thus, the upper envelope of is less oscillating than is. The same behavior occurs for the mean envelope. We would then try to interpret locally the relationship by evoking the maximum principle [24]. But, this principle is not immediately checked for equations of the type (15). The justification lies in the fact that smoothing implies .

In summary, we have just built a PDE system (1)having the input signal to decompose as an initial value condition, (2)such that the result of a diffusion process preserves characteristic points of , (3)with a regularity similar to a cubic spline.

3.2. Formulation of the Mathematical Problem

Let be one of the above constructed interpolation operators, and be a diffusion time sufficiently long. To fix ideas, we solve the problem of upper envelope calculus. Indeed, for the other kinds of characteristic points, the same problem arises with analogous mathematical resolution methods.

The problem posed by our mathematical modeling is to find a solution of the system Then, the mean envelope (interpolating turning points) calculus is stated as follows.

For given , find such that

3.3. Comments on the Degeneracy in Problem (24)

To prevent the zero degeneracy case in the system (25) whose solution is not directly accessible by variational methods, we are led to solve an intermediate problem with That is the nondegenerate problem defined hereafter.

3.4. Formulation of the Zero Degeneracy Induced by Diffusivity Function

The diffusivity function in (26) is very close to zero at the characteristic points while refraining to cancel and to solve the system (25), now, just to move on to the limit when tends to zero to retrieve (25) and its solution.

Thus, in particular, we can define the function which verifies the conditions of (26).

The new system that we call nondegenerate problem is fomulated as follows.

Find such that

3.5. Formulation of the Nondegenerate Problem

For sufficiently large, is a close approximation of the envelope , interpolating the characteristic points.

By a density technique, we can demonstrate that this sequence converges to the solution of the system (25). By posing for the convenience of notation , the non-degenerate problem is reworded as.

Find such that

Equation (29) is an operational differential equations type. The techniques for the resolution of this kind of problem are numerous we mainly refer to the variational formulation to find weak solutions [25] and the method of the approximation of evolutionary operators. The last allows to work on a new elliptic problem. The solution is obtained by passing to the limit of the approximate solution nstead of the equation.

In the next section, we solve (29) by a variational approach, according to the resolution methods for parabolic problems described in [7, 20].

4. Existence and Uniqueness of Solutions for PDE Interpolator

4.1. The Guideline for the Resolution of the Mathematical Problem

Now, we summarize the general procedure to solve the problem. (1)First step: resolution of the problem (29). This step includes (i)the variational formulation in Section 4.3, (ii)the resolution of the variational problem in appropriate functional space,(iii)the converse showing that the solution of the variational problem solves the departure problem in Section 4.3.4. (2)Second step: in Section 5, construction of the sequence of solutions of problem (29), where is the solution obtained with the diffusivity. ,  , see Section 5. Finally, demonstrate that the sequence converges to a limit , where is the solution of problem (25), in Section 5.

4.2. The Main Results Formulation

The main result which is going to be solves is the following.

Theorem 8. For , exists a unique solution of the variational problem such that and . Reversely, let be a solution of the variational problem (30), then .

And the complete follows result.

Theorem 9. For , exists a unique solution verifying

4.3. Variational Formulation

Let , and (to make sense to integrals) multiply the equation of the system (29) by the test function after we integrate on , then it comes out that By using the Green formula (see Definition 4) and integration by parts, we get firstly and secondly, As we are in one-dimension case, the normal derivative and the gradient on the edge in the use of the Green formula are the same, that is to say, In the following, we adopt the notation .

Considering such that it comes out that Let us consider the following bilinear forms The addition of boundaries conditions transforms our initial problem into the system

The variational formulation of our problem consequently, subject to additional conditions that would cause a change of space, is the following where and are like in Theorem 7.

This kind of problem is studied in [21]. However, before continuing the resolution, we explore some properties of the space .

4.3.1. On the Quotient Space

The space is defined by , where and are the domains of and . Let us consider the quotient space and the new space   defined by Let , and let be a representative of class, then is defined by The nullity of means ,  constant such that  ,  with    one representative of  

So for all ,  we have  . Accordingly, any element of the quotient space is, on one hand, in with normal derivative on the boundaries of null and, on the other hand, is zero average on . Thus, the quotient space is identified with the space . We can work later in the space with zero mean. The final variational problem to solve is the following

4.3.2. Application of the Lions Theorem

To apply the Lions theorem, we need to identify the applications ,  , and and the spaces and . And finally, we must prove that the norm defined in , which derives from the scalar product given by the bilinear application , is equivalent to the norm of .

Let us define as the identity, and let us define and by , and , obtained from (39), (40); the spaces and are define as follows.

(a) The linear operator : in the variational problem (11) in the Lions theorem, denotes a scalar product on . Indeed, we define .

The operator , , is the identity which is continuous (verify the Lions conditions Theorem 7), and the term is equal to In this identification, is a well scalar product in the subspace of , formed by functions with zero mean. Indeed, according to Poincare-Wirtinger inequality (see Definition 5) in , where is the first positive eigenvalue of Laplace operator in . But as ,  then  , consequently there exists such that: that is to say the norm is equivalent to the gradient norm in . Hence, one can replace the scalar product of the gradient by scalar product.

Hereinafter, we complete the proof that the precedent choices allow the application of the Lions theorem in our context.

(b) is an Hilbert space with equivalent norm to norm.

Theorem 10. The space , resulting from (11) and analysis, defined by is an Hilbert space with the a norm equivalent to norm.

Proof. At first is a closed vector subspace of . Note that the condition of the nullity of the mean in may be replaced by , where is a nonempty open subset of . In this case becomes rather From the fact that the norms built with the two bilinear applications and are equivalent to the norm constructed by the linear form, noted by with , and given by Indeed, is as well a norm in . There is no difficulty to formally verify the triangular inequality. If , then and . Consequently, is necessarily constant almost everywhere in , for example, is equal to a constant real . But as   implies , so . We have thus built a standard scalar product on from bilinear forms of the variational formulation. is a Hilbert space with the scalar product associated with the norm (54). In addition, this norm is equivalent to the norm, as it shown in Lemma A.2. Now, we just need to prove that the correspondence is continuous and coercive on the Hilbert space .

4.3.3. The Application Is Continuous and Coercive

For , hence the continuity of is on . In order to complete the application conditions of Theorem 7, we only check the first condition of the same theorem.

For hence nevertheless Thus, We have found thus such that All the conditions of Theorem 7 are acquired; now, the main existence result for solution of system (41) can be enunciated as follows.

Theorem 11. For , exists a unique solution of the problem

4.3.4. Equivalence with the Initial Problem

Reversely, let be a solution of the variational problem As we must consider a test function

Let us take a function test and such that

We pose the following Neumann problem:

Its variational formulation gives a unique solution in , and the regularity of implies the regularity of the solution which is in , consequently .

Returning to (62) with , solution of the problem (65), the first term can be written as However, , subsequently The second term of (62) is written as follows: which gives in the distribution sense from which We deduce that We verify now that . Returning to (70) with and according to the Green formula, we have But as , it gives according to (71) , and consequently on boundaries of , or as ,    on boundaries of .

In conclusion is solution of the initial equation (41), then we enunciate the following theorem.

Theorem 12. For , exists a unique solution , verifying

4.3.5. Some Remarks on the Modeling Space

The nullity condition of the average for functions in has been decisive in the construction of the norm on this space. This space is prima facie clearly not natural for EMD because of entry function that is not necessarily zero mean. But at the cost of working with an initial condition , such that , we can consider as functions with null local mean.

The Nullity of the Mean Envelope: A Characteristic of and EMD Algorithm. This condition, far from being superfluous, was even expected because the space is the space of the mean envelopes which must be zero mean in EMD stifting process for mode extraction. If it were to calculate the upper or lower envelopes, it should be noted that the condition of the mean nullity raised a query. In EMD algorithm, this nullity condition is affected on the mean envelope and not on upper or lower envelopes. This question does not arise in the model giving the envelope interpolating inflection points of a signal. The elements of are not zero mean on , but on a non-empty subset of .

The nullity of the mean of the signal on a is not immediately visible and warranty to decompose any function. We show in Lemma A.1 that each relevant function to be decomposed by EMD is an integrable function possessing at least three extrema belonging to . The modeling space of the extrema interpolation is defined in (52).

Let be a function admitting at least three extrema we can consider that passes at least once a zero at a point in its definition domain. Otherwise, we just work with the translatory (as our model is invariant by translation of the input function) vector which is equals to the half amplitude of .

5. Convergence of the Sequence of Solutions of (41) to a Solution of the Called Degenerate Problem

In the following sections, we demonstrate that the sequence of solutions of the non degenerate problem is bounded. Next, we prove that there exists a subsequence of which converges weakly to an element and, finally, that this element is solution of the degenerate initial problem.

5.1. Some Estimations

From (28) we have in the distributions sense (multiplying by ): which is equivalent to After integration by parts, it comes out that then by integrating with respect to the variable , it comes out that which implies that But as , we have or from which Thus, on one hand, But where comes from inequality (49), which implies that Then, the sequence   is bounded in and in . And on the other hand which gives But as therefore Thus, which implies according to (84) that We finally deduced the second estimate Then, the sequence is bounded in .

5.2. Weak Convergence of the Subsequence of

The sequence is bounded in the space . Then there exists a subsequence that converge weakly to in with .

Leting be a test function, we have the following.

Firstly. As weakly in , therefore, weakly in .

Thus, in the distributions sense, by integrating on , we have By integrating by parts again the last equation, we have This expression converges to Therefore,converges when   tends to infinity to Consequently, weakly in .

Secondly. As , then, converges to zero when tends to infinity.

On the other hand, from (98), we have which converges to zero in according to weak convergence (see Definition 3) in   of    to  . Thus, we have shown that in the distributions sense of the one hand converges to and on the other hand, converges to Consequently, for each , converges in the distributions sense to that is to say, converges when   tends to infinity to . It remains to show that the limit is a solution of the degenerate problem. But , which implies in the distributions sense . In other words, the weak limit of the sequence is weak solution of the initial degenerate problem. Moreover, this solution is unique and is in . Better, using the compact inclusion of (which behaves like because of the equivalence of the norms) in , we conclude that the sequence of solutions of non-degenerate problem (approximated problem) converges strongly in to the solution of the initial degenerate problem. The space is, moreover, the space of strong solutions of the initial system (25).

6. Numerical Implementation and Applications

Deliberately, we present the numerical implementation in the case of dimension . The reason is that we give examples both on 1D and 2D signals.

In the sequel, we adopt the following notation. The mesh of the domain is given by

We note by the approximate value of the solution at the point .

The discrete solution at the iteration is obtained at the discrete time and is denoted by . An approximation of the temporal derivative is at this same time is then .

6.1. Numerical Schemes for the PDE Resolution

The discrete model calculating the mean envelope can be written as follows: where denotes the derivative of order relatively to the spatial coordinate .

To take into account the discontinuity of the diffusivity function, we use the harmonic mean of defined by Recall that is strictly positive, which allows its inversion in the previous system.

We define also the numerical derivatives of by and finally we pose . For positive integer, and , the numerical approximation , of in (108) is given by We present the first natural explicit scheme.

6.1.1. An Explicit Scheme

An explicit scheme to solve the PDE is the following: or . Let us pose .

Thus, If is the image or the signal size, the matrix is a sparse matrix of order .

However, this explicit scheme requires for its stability very small. To overcome this drawback, other numerical schemes can be used such as Du Fort and Frankel [22], which is unconditionally stable, but conditionally consisting of time steps not too great, to satisfy the condition of numerical stability of Courant-Friedrich-Levy [26] schemes. Another very effective scheme is proposed by Vogel and Oman in [27] it is based on the principle of the fixed point. It allows us to address the problem stationary with the result of the existence of the asymptotic solution. With this scheme, the solution is an eigenvector of the operator represented by the matrix .

6.1.2. The Additive Operator Splitting Scheme

We can use the additive operator splitting scheme (AOS) which is written as follows: This scheme is more stable, but requires a matrix inversion, even if the matrix is pentadiagonal. This inversion may be affected by the Thomas [28] algorithm modified. In (113), the reverse of the matrix is independent of the iteration and is carried out once only for .

6.1.3. The Alternate Direction Implicit Scheme (ADI)

Simplification which reduces the complexity of the computation of the solution can be made in the implementation of the numerical equation (108) combining an implicit form for half the computation time and an explicit form for the remainder.

It is more difficult to implement this scheme, ADI, for the fourth-order PDE. This is due to cross derivatives in (108). Witelski and Bowen [29] suggest a readjustment of the ADI in which a mixed derivative is calculated explicitly in (108) by the following: Besides the approximation errors inherent in all numerical scheme, the fundamental difference between theses methods remains the computation time.

6.2. The Boundary Conditions

The boundary conditions used are those of system (41). In fact, the solution is a polynomial of order 3 and can be written as . The nullity of its derivatives of orders and on the edges leads to ,   on the boundary. Hence, , which finally gives . We have taken as boundary conditions that outside, the domain, the signal is constant. This means into practice the cancellation of all the derivatives of order less than or equal to . For diffusivity function at the edges, it is constant. When this constant is null, we join a nonnatural manner to care of edge effects, which was to take the signal ends at the same time as local maximum and local minimum. We do not this here.

6.3. Test Results of PDE Interpolator on Signal

In our numerical tests, we used the implicit scheme of Crank-Nikelson which gives good results with that can take very large values. In the following, we give some examples of calculation of upper and lower envelopes and for the mean envelope interpolating the inflection points of the signal. This scheme in consistent of order both in space and in time.

6.3.1. Some Examples

First, we present in Figure 1 a pathological case of a signal without local extremum, when we cannot calculate the envelopes by spline interpolation of the extrema. Our PDE interpolator works well with the max-min curvature points detection in the diffusivity function.

Figure 1: PDE interpolator for signal without local extrema. The input signal is in blue unbroken line, and the envelopes are in dotted lines.

Figure 2 represents the evolution of the PDE solution with the form given by (14). Upon iteration , the solution almost evolves addition; it is potentially an envelope.

(a)   Input signal is in blue color. In red and green colors are the envelopes upon iteration 4000

(b) Input signal is in blue color. In red and green colors are the envelopes upon iteration 40000

(c) Input signal is in blue color. In red and green colors are the envelopes upon iteration 400000

Figure 2: In (a) input signal is in blue color. In red and green colors are the envelopes calculus upon iteration . In (b) the envelopes calculus upon iteration . In (c) the envelopes calculus upon iteration .

In Figure 3, we represent the envelope calculation of the same function with the form given by (16) for different tension factor . We can see that the regularity of envelopes is no longer for values of close to zero.

(a) Input signal is in blue color. In red and green colors are the envelopes at the convergence for

(b) Input signal is in blue color. In red and green colors are the envelopes at the convergence for

(c) Input signal is in blue color. In red and green colors are the envelopes at the convergence for

Figure 3: In (a, b, and c) input signal is in blue color. In red and green colors are the envelopes at the convergence of the evolutive PDE.

In Figure 4, we represent the envelope calculation with the form given by (17) for different tension factor .

(a)   Input signal in blue color. In red and green colors the envelopes upon iteration 4000

(b) Input signal in blue color. In red and green colors the envelopes upon iteration 40000

(c) Input signal in blue color. In red and green colors the envelopes upon iteration 400000

Figure 4: Envelope calculus by (17) with . In (a), at this stage, the two envelopes overlap by location. In (b), less overlap and practically the asymptotic solution are almost reached. In (c), the stability of the evolution leading to the solution.

In Figure 5, we make a comparison between our PDE interpolator and the cubic spline interpolation for an input signal that equal to the superposition of two sinusoidal components with different amplitudes and frequencies, and . When it comes to interpolate the classical local extrema, we have almost the same result, the two not being satisfactory, which certainly requires a bit more iterations in the sifting process. In contrast by, our interpolator which passes through the points of maximum and minimum local curvature gives the correct envelopes quite properly.

Figure 5: Comparison between PDE interpolator and cubic spline. The input signal is in blue color, and in red and black colors are the envelopes.

6.3.2. The Specific Case of the PDE Interpolator of the Turning Points

In the case of the PDE interpolator of inflection points, we have the advantage of not only calculating asymptotic solution of PDE but also dividing the computing time almost in half. Only, we must numerically have a good estimate of the envelope to guard against a sampling of the input signal which must also be quite dense. This model asks more numerically regularity and in some cases an over-sampling.

In Figures 6 and 7, we proceed to calculate the solution of the PDE interpolator passing through the inflection points, for different sampling for the same composite signal . The convergence is slower if the number of sampling points is higher, but the solution is a good estimate of the mean envelope. In Figure 6, with a more dense sampling, convergence is slower, and we get a mean envelope judged more suitable than 128 sampling points. In Figure 7, with adequate sampling, we presented the mean envelope expected trend, that is the component of the signal . To capture the inflection points, a sampling of four once of the smallest local period of the signal is generally eligible for a suitable to estimate the envelope. This is a consequence of the theory of Shannon sampling.

(a)   Envelope after 40 iterations

(b) Envelope after 400 iterations

Figure 6: The envelope of the signal test sampling with 256 points.

(a)   Envelope after 40000 iterations

(b) Envelope after 400000 iterations

Figure 7: The envelope of the signal test sampling with 1024 points.

6.3.3. Test on 2D Signal

Finally, in Figure 8, we use the 2D version of the PDE interpolator for image restoration. The corrupted image taken as an input 2D signal in our model which proceeds by edge propagating provides an interpolated version of the data that is clean enough.

(a)

(b)

(c)

Figure 8: PDE interpolator for image restoration. An original image (a) is corrupted (b) and restored (c) by our 2D PDE-interpolator. Original image.

7. Conclusion

In this paper, we model a new PDE interpolator which permits to calculate the envelopes of a signal. The existence and uniqueness of the solution to the mathematical problem are determined by a variational approach. This theoretical framework contributes to the mathematical modeling of EMD algorithm. We have shown, in particular, how the space reminds the extraction conditions of the EMD modes, and we prove that the set of eligible functions to be decomposed by EMD are in . The PDE interpolator is not based on any prior knowledge on the noise level as opposed to the total variation method in [11]. The tests in both signal and image processing demonstrate the effectiveness of the new interpolator and provide an insight into opportunities in multiscale analysis of multidimensional signal.

Appendix

Lemma A.1. Let having at least three local extrema and null at least at one point of . One can assume that vanishes by changing sign. Then, there exist and such that the local mean at is null, for example,

Proof. Let and such that ,  , and . Proceed by contradiction as the local mean is continuous, and not null implies that it keeps a constant sign on . Suppose that ,  ,  . In particular, the positive sequence converges to . This is absurd because . Thus, the proof of Lemma A.1 is performed.

From this lemma, any eligible function for EMD is in .

Lemma A.2. There exist two constants such that :

Proof of Lemma A.2  Step  1  (second part of (A.2)). Rewriting first the norms, , that is to say, and secondly, we have immediately posing .Step  2  (first part of (A.2)). Firstly, and by integration by parts, As ,  we have and, therefore, the use of Schwarz inegality gives then, we deduce More  , then according to Poincaré-Wirtinger theorem [20], it comes out that .Therefore, Thus, If we pose we have The inegality (A.2) is established and the lemma’s proof is completed.

References

N. E. Huang, Z. Shen, S. R. Long et al., “The empirical mode decomposition and Hilbert spectrum for nonlinear and non-stationary time series analysis,” Proceedings of the Royal Society A, vol. 545, no. 1971, pp. 903–995, 1998.

E. Deléchelle, J. Lemoine, and O. Niang, “Empirical mode decomposition: an analytical approach for sifting process,” IEEE Signal Processing Letters, vol. 12, no. 11, pp. 764–767, 2005.

O. Niang, Empirical mode decomposition: contribution à la modélisation mathématique et application en traitement du signal et l'image [Ph.D. thesis], University Paris-Est Créteil, Paris, France, 2007.

O. Niang, E. Deléchelle, and J. Lemoine, “A spectral approach for sifting process in empirical mode decomposition,” IEEE Transactions on Signal Processing, vol. 58, no. 11, pp. 5612–5623, 2010.

O. Niang, A. Thioune, E. Delechelle, and J. Lemoine, “Spectral intrinsic decomposition method for adaptive signal representation,” ISRN Signal Processing, vol. 2012, Article ID 457152, 10 pages, 2012.

O. Niang, A. Thioune, M. C. El Gueirea, E. Deléchelle, and J. Lemoine, “Partial differential equation-based approach for empirical mode decomposition: application on image analysis,” IEEE Transactions on Image Processing, vol. 21, no. 9, pp. 3991–4001, 2012.

F. Catte, P. L. Lions, J. M. Morel, and T. Coll, “Image selective smoothing and edge detection by nonlinear diffusion,” SIAM Journal on Numerical Analysis, vol. 29, no. 1, pp. 182–193, 1992.

P. Perona and J. Malik, “Scale-space and edge detection using anisotropic diffusion,” IEEE Transactions on Pattern Analysis and Machine Intelligence, vol. 12, no. 7, pp. 629–639, 1990.

Y. Meyer, Oscillating Patterns in Image Processing and Nonlinear Evolution Equations, vol. 22 of University Lecture Series, AMS, Providence, RI, USA, 2002.

S. Osher, A. Sole, and L. Vese, “Image decomposition and restoration using total variation minimization and the H1 norm,” Multiscale Modeling and Simulation: A SIAM Interdisciplinary Journal, vol. 1, no. 3, pp. 349–3370, 2003.

M. Lysaker, A. Lundervold, and X. C. Tai, “Noise removal using fourth-order partial differential equation with applications to medical magnetic resonance images in space and time,” IEEE Transactions on Image Processing, vol. 12, no. 12, pp. 1579–1589, 2003.

J. Weickert, B. M. ter Haar Romeny, and M. A. Viergever, “Efficient and reliable schemes for nonlinear diffusion filtering,” IEEE Transactions on Image Processing, vol. 7, no. 3, pp. 398–410, 1998.

J. Weickert, “A review of nonlinear diffusion filtering,” in Scale-Space Theory for Computer Vision, B. H. Romeny, Ed., vol. 1252 of Lecture Notes in Computer Science, pp. 3–28, Springer, New York, NY, USA, 1997.

D. W. Peaceman and H. H. Rachford, “The numerical solution of parabolic and elliptic differential equations,” Journal of the Society For Industrial and Applied Mathematics, vol. 3, no. 1, pp. 28–41, 1955.

D. Barash and R. Kimmel, An Accurate Operator Splitting Scheme for Nonlinear Diffusion Filter, HP Company, 2000.

J. Tumblin and G. Turk, “LCIS: a boundary hierarchy for detail-preserving contrast reduction,” in Proceedings of the SIGGRAPH 1999 Annual Conference on Computer Graphics, Los Angeles, Calif, USA, 1999.

G. W. Wei, “Generalized Perona-Malik equation for image restoration,” IEEE Signal Processing Letters, vol. 6, no. 7, pp. 165–167, 1999.

Y. L. You and M. Kaveh, “Fourth-order partial differential equations for noise removal,” IEEE Transactions on Image Processing, vol. 9, no. 10, pp. 1723–1730, 2000.

J. Tumblin, Private Communication, 2003.

H. Brezis, Analyse fonctionnelle. Théorie et Applications, Masson, Paris, France, 1983.

J. L. Lions, Équations Différentielles Opérationnelles et Problèmes aux Limites, Springer, Berlin, Germany, 1961.

J. D. Murray, Mathematical Biology, Springer, Berlin, Germany, 1993.

P. Flandrin, G. Rilling, and P. Gonçalvés, “Empirical mode decomposition as a filter bank,” IEEE Signal Processing Letters, vol. 11, no. 2, pp. 112–114, 2004.

R. Sperb, Maximum Principle and Their Applications, vol. 157 of Mathematics in Science and Engineering Series, Academic Press, New York, NY, USA, 1981.

H. Brezis, Operateurs Maximaux monotones et semigropes de contraction dans les espaces de Hilbert, North-Holland, Amsterdam, The Netherlands, 1972.

G. Aubert and P. Kornprobst, Mathematical Problems in Images Processing, Partial Differential Equations and the Calculus of Variations, vol. 147 of Applied Mathematical Sciences Series, Springer, Berlin, Germany, 2002.

C. Vogel and M. Oman, “Iteration methods for total variation denoising,” SIAM Journal on Scientific Computing, vol. 17, no. 1, pp. 227–238, 1996.

G. Engeln-Muellges and F. Uhlig, Numerical Algorithms with C, chapter 4, Springer, Berlin, Germany, 1996.

T. P. Witelski and M. Bowen, “ADI schemes for higher-order nonlinear diffusion equations,” Applied Numerical Mathematics, vol. 45, no. 2-3, pp. 331–351, 2003.