`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

1Département Informatique et Télécommunications, Ecole Polytechnique de Thiès (EPT), Thiès BP A10, Senegal
2Laboratoire Images, Signaux et Systèmes Intelligents (LISSI-E.A.3956), Université Paris-Est Créteil Val-de-Marne, Créteil, France
3Laboratoire d’Analyse Numérique et d’Informatique (LANI), Université Gaston Berger (UGB), Saint-Louis BP 234, Senegal
4Dé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 [24]. 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 [911]. Efficient numerical schemes were introduced in [12] based on additive operator splitting (AOS) scheme, or based on alternating direction implicit (ADI) scheme. See [1215] 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 [1618]. 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