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

Authors: Oumar Niang; Abdoulaye Thioune; Éric Deléchelle; Mary Teuw Niane; Jacques Lemoine
Journal: ISRN Signal Processing
 (2013)

---

## 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.


---

## Body

## 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
(1)ut(x,t)=div(g(|∇u(x,t)|2)∇u(x,t)),
with x=(x1,x2), and where u(x,t) is a filtered version of the original image u(x,t)=u0(x) as the initial condition, with reflecting boundary conditions. In (1), g(·) 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
(2)g1(x,t)=11+(|∇u(x,t)|/β)2,g2(x,t)=exp(-(|∇u(x,t)|/β)2).
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:
(3)ut(x,t)=-div(g(m(u))∇Δu(x,t)),
where g(·)=g1(·) as in (2) and m is some measurement of u(x,t). In [16], (3) is called a “low curvature image simplifier” (LCIS), and a good choice for m is defined as m=Δu 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 ℝn,  n=1  or  2. In the sequel, we will need the following definitions and results.Definition 1.
we define one the spacesV and 𝒱 as follows
(4)V={v∈H3(Ω)∣∂v∂ν|∂Ω=0},𝒱={v∈H3(Ω)∣∂v∂ν|∂Ω=0,∫Ωvdx=0},𝒱e={∫ΓΩ‍vdx=0andv∈H3(Ω)suchthat∃ΓΩ⊂Ω,openandverifying∫ΓΩvdx=0and∂v∂ν|∂Ω=0}.Definition 2.
For a given functionv∈H1(Ω), if Dd(v)(x) and Dg(v)(x) denote, respectively, the right and left derivative of v at x; the minmod derivatives of v is define by
(5)∂v∂x(x)=minmod(Dg(v)(x),Dd(v)(x))whereminmod(p,q)={sign(p)min(|p|,|q|),ifpq>0,0,ifpq≤0.
A H1(Ω) function being absolutely continuous admits right and left derivatives, then u0∈H1(Ω) has obviously left and right derivatives, so that we can validate numerically computing of the diffusivity function g defined in Section 3.Definition 3.
Let(xn)n∈ℕ be a sequence of elements in a vectorial normed space (E,∥·∥E); it is said that (xn)n∈ℕ converge, weakly [20] in E, and noted by xn⇀x, if exists an element x∈E such that ∀f∈E′,  limn→∞f(xn)=f(x), where E′ denotes the set of continuous linear forms on E.Definition 4 (the Green formula).
Letu,v∈C2(Ω), then
(6)∫Ω‍Δudx=∫∂Ω‍∂u∂nds,∫Ω‍∇u∇vdx=-∫Ω‍uΔvdx+∫∂Ω‍∂v∂nds,∫Ω‍(uΔv-vΔu)dx=∫∂Ω‍u∂v∂n-v∂u∂ndsGreenformula,
where ∂u/∂n denotes the normal derivative of u on the boundary ∂Ω of Ω.Definition 5 (the Poincaré-Wirtinger inequality).
LetΩ be an open-bounded set, and let u∈H1(Ω), then there exists a constant C>0 such that the norm of u in H1(Ω) and the norm of u in L2(Ω) are linked by the following inequality:
(7)∥u-1|Ω|∫Ω‍uds∥H1(Ω)≤C∥∇u∥L2(Ω),
where |Ω| denotes the length of Ω. The best constant C in the Poincaré-Wirtinger inequality is 1/λ, 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(8)u∈{u∈H3(Ω)suchthat∂u∂n|∂Ω=0,∫Ω‍udx=0},
then there exist two strictly positive reals C and K such that firstly,(9)∥∇u∥L2(Ω)2≤C∥Δu∥L2(Ω)2,
and secondly,
(10)∥u∥H3(Ω)2≤K∥Δu∥H1(Ω)2.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).
LetH and V be two Hilbert spaces with V⊂H. One considers a bilinear application a on V×V and B an operator on V. Under the following conditions:  (1)
∃μ and ρ>0 such that |a(u,v)|≥ρ∥v∥V2-μ|v|H2,(2)
∃c2>0 such that |(v,Bv)V|≤c2|v|V2.
Foru0=u(0)∈H, the problem
(11)u∈L2(0,T;V)∩C(0,T;H),
such that
(12)∀v∈V,(dudt,Bv)H+a(u,v)=0
and u0=u(0), has a unique solution.
Furthermore,du/dt∈L2(0,T;H′).


## 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
(1)ut(x,t)=div(g(|∇u(x,t)|2)∇u(x,t)),
with x=(x1,x2), and where u(x,t) is a filtered version of the original image u(x,t)=u0(x) as the initial condition, with reflecting boundary conditions. In (1), g(·) 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
(2)g1(x,t)=11+(|∇u(x,t)|/β)2,g2(x,t)=exp(-(|∇u(x,t)|/β)2).
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:
(3)ut(x,t)=-div(g(m(u))∇Δu(x,t)),
where g(·)=g1(·) as in (2) and m is some measurement of u(x,t). In [16], (3) is called a “low curvature image simplifier” (LCIS), and a good choice for m is defined as m=Δu 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.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
(1)ut(x,t)=div(g(|∇u(x,t)|2)∇u(x,t)),
with x=(x1,x2), and where u(x,t) is a filtered version of the original image u(x,t)=u0(x) as the initial condition, with reflecting boundary conditions. In (1), g(·) 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
(2)g1(x,t)=11+(|∇u(x,t)|/β)2,g2(x,t)=exp(-(|∇u(x,t)|/β)2).
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:
(3)ut(x,t)=-div(g(m(u))∇Δu(x,t)),
where g(·)=g1(·) as in (2) and m is some measurement of u(x,t). In [16], (3) is called a “low curvature image simplifier” (LCIS), and a good choice for m is defined as m=Δu 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 ℝn,  n=1  or  2. In the sequel, we will need the following definitions and results.Definition 1.
we define one the spacesV and 𝒱 as follows
(4)V={v∈H3(Ω)∣∂v∂ν|∂Ω=0},𝒱={v∈H3(Ω)∣∂v∂ν|∂Ω=0,∫Ωvdx=0},𝒱e={∫ΓΩ‍vdx=0andv∈H3(Ω)suchthat∃ΓΩ⊂Ω,openandverifying∫ΓΩvdx=0and∂v∂ν|∂Ω=0}.Definition 2.
For a given functionv∈H1(Ω), if Dd(v)(x) and Dg(v)(x) denote, respectively, the right and left derivative of v at x; the minmod derivatives of v is define by
(5)∂v∂x(x)=minmod(Dg(v)(x),Dd(v)(x))whereminmod(p,q)={sign(p)min(|p|,|q|),ifpq>0,0,ifpq≤0.
A H1(Ω) function being absolutely continuous admits right and left derivatives, then u0∈H1(Ω) has obviously left and right derivatives, so that we can validate numerically computing of the diffusivity function g defined in Section 3.Definition 3.
Let(xn)n∈ℕ be a sequence of elements in a vectorial normed space (E,∥·∥E); it is said that (xn)n∈ℕ converge, weakly [20] in E, and noted by xn⇀x, if exists an element x∈E such that ∀f∈E′,  limn→∞f(xn)=f(x), where E′ denotes the set of continuous linear forms on E.Definition 4 (the Green formula).
Letu,v∈C2(Ω), then
(6)∫Ω‍Δudx=∫∂Ω‍∂u∂nds,∫Ω‍∇u∇vdx=-∫Ω‍uΔvdx+∫∂Ω‍∂v∂nds,∫Ω‍(uΔv-vΔu)dx=∫∂Ω‍u∂v∂n-v∂u∂ndsGreenformula,
where ∂u/∂n denotes the normal derivative of u on the boundary ∂Ω of Ω.Definition 5 (the Poincaré-Wirtinger inequality).
LetΩ be an open-bounded set, and let u∈H1(Ω), then there exists a constant C>0 such that the norm of u in H1(Ω) and the norm of u in L2(Ω) are linked by the following inequality:
(7)∥u-1|Ω|∫Ω‍uds∥H1(Ω)≤C∥∇u∥L2(Ω),
where |Ω| denotes the length of Ω. The best constant C in the Poincaré-Wirtinger inequality is 1/λ, 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(8)u∈{u∈H3(Ω)suchthat∂u∂n|∂Ω=0,∫Ω‍udx=0},
then there exist two strictly positive reals C and K such that firstly,(9)∥∇u∥L2(Ω)2≤C∥Δu∥L2(Ω)2,
and secondly,
(10)∥u∥H3(Ω)2≤K∥Δu∥H1(Ω)2.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).
LetH and V be two Hilbert spaces with V⊂H. One considers a bilinear application a on V×V and B an operator on V. Under the following conditions:  (1)
∃μ and ρ>0 such that |a(u,v)|≥ρ∥v∥V2-μ|v|H2,(2)
∃c2>0 such that |(v,Bv)V|≤c2|v|V2.
Foru0=u(0)∈H, the problem
(11)u∈L2(0,T;V)∩C(0,T;H),
such that
(12)∀v∈V,(dudt,Bv)H+a(u,v)=0
and u0=u(0), has a unique solution.
Furthermore,du/dt∈L2(0,T;H′).


## 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 u0∈L2(Ω). We construct an operator A of appropriate domain D(A) as follows:
(13)∂u∂t+Au=0in[0,T]×Ω,u(0,x)=u0inΩ,plusboundariesconditions.The need of asymptotic solution existence denoted byu∞ implies (∂u/∂t)(∞,x)=0, for example, according to (13), we have Au(∞,x)=0. Moreover, the requirement of the same regularity between the solution and a cubic spline leads to ∂4u/∂x4(∞,x)=0. Thus, in the first analysis, we can choose A=∂4u/∂x4. To leave invariant features points during the diffusion process, simply multiply in the expression Au the term ∂4u/∂x4 by function g depending on spatial variable x and vanishing at characteristic points. That gives
(14)(Au)(x)=g(x)∂4u∂x4.
Another reminiscent form of long-range diffusion [22] is the following:
(15)(Au)(x)=∂∂x(g(x)∂3u∂x3).
Function g 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)
(16)(Au)(x)=g(x)[-θ∂2u∂x2+(1-θ)∂4u∂x4],
or for (15),
(17)(Au)(x)=∂∂x[g(x)(-θ∂u∂x+(1-θ)∂3u∂x3)],
or accomplete diffusion involving a flow
(18)(Au)(x)=∂∂x[-θg(x)∂u∂x+(1-θ)∂∂x(g(x)∂2u∂x2)],
where 0≤θ≤1 is the tension factor which controls the regularity of the solution.Both forms allow more freedom on the regularity of the solutionu. To access the mathematical properties of the solution, we chose the second form in (17) with θ=0 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 forg to calculate the envelope is given by a positive piecewise function lower than 1 that is constant between two characteristic points of u0 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 u0. We characterize the function g as depending on sign(∂u0/∂x),sign(∂2u0/∂x2), and sign(∂3u0/∂x3).For the purpose of existence and regularity of the solution, we are led to work with the regularized version ofsign:
(19)signϑ(z)=2πarctan(πzϑ),
where ϑ is a regularization coefficient.In the following, we definege± for extrema detection, gf for turning points detection, and gmc± for maximum and minimum curvature points detection:
(20)ge±(x)=19[|sign(∂u0∂x)|±sign(∂2u0∂x2)+1]2(+)formaximaand(-)forminima,gf(x)=[signϑ(∂2u0∂x2)]2forturningpoints,gmc±(x)=19[|sign(∂3u0∂x3)|±sign(∂2u0∂x2)+1]2,
(+) for maximum curvature points and (−) for minimum curvature points. All these functions are of the form g(x)=[h(x)]2  with h(x)=0 at characteristic points.So, we haveg′(x)=0if g(x)=0. This property formally allows us to cancel Au and ∂u/∂t of the form (15) at the points where g is null. The factor 1/9 is a normalization term and verifies 0≤g(x)≤1. A more general diffusivity function g could be envisaged.Let us consider the following differential operatorsAe and Af (e for upper and f for inflexion or curvature):
(21)∀u∈D(Ae),Aeu=∂∂x(ge±(x)∂3u∂x3),
or
(22)∀u∈D(Af),Afu=∂∂x(gf(x)∂3u∂x3).Ifu0∈C2(Ω), 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 u0 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 C2(Ω), for example, for signal u0∈H1(Ω), as it is the case in reality [23]. But numerically, without even a regularization function sign, if the derivatives are taken in the sense of minmod called flux limiter or slope, the function g(x) still has a meaning.


### 3.1. Interpretation of the Diffusivity Action in the Diffusion Process

Between two consecutive maximaXmaxi and Xmaxi+1, the function ge+ is piecewise constant and can be written as follows:
(23)g(x)={0,atXmaxi,Xmaxi+1,49,atXinfi,Xinfi+1,Xmini,1,on]Xinfi,Xmini[∪]Xmini,Xinfi+1[,19,sur]Xmaxi,Xinfi[∪]Xinfi+1,Xmaxi+1[.We have a diffusion in the sense thatut=∂u/∂t unhook to the maximum of the curve u0. The reason is that the diffusion is more pronounced in local minimum because the smoothing effect tends to regularize curves. Thus, the upper envelope u(t≅∞,·) of u0 is less oscillating than u0 is. The same behavior occurs for the mean envelope. We would then try to interpret locally the relationship u≥u0 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 ∥∇u(t,·)∥L2(Ω)≤∥∇u0∥L2(Ω).In summary, we have just built a PDE system(1)
having the input signalu0 to decompose as an initial value condition,(2)
such that the result of a diffusion process preserves characteristic points ofu0,(3)
with a regularity similar to a cubic spline.


### 3.2. Formulation of the Mathematical Problem

LetA be one of the above constructed interpolation operators, and T>0 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 solutionu of the system
(24)∂u∂t+Au=0in[0,T]×Ω,u(0,x)=u0inΩ,plusboundariesconditions.
Then, the mean envelope (interpolating turning points) calculus is stated as follows.For givenu0, find u such that
(25)∂u∂t+∂∂x(gf(x)∂3u∂x3)=0in[0,T]×Ω,u(0,x)=u0inΩ,plusboundariesconditions.


### 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
(26)0<α≤gf,α(x)≤β≤1.
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(27)gf,n(x)=gf(x)+1n,n∈ℕ*
which verifies the conditions of (26).The new system that we call nondegenerate problem is fomulated as follows.Findun such that
(28)∂un∂t+∂∂x(gn(x)∂3un∂x3)=0in[0,T]×Ω,un(0,x)=u0inΩ,plusboundariesconditions.


### 3.5. Formulation of the Nondegenerate Problem

Forn sufficiently large, un is a close approximation of the envelope u, interpolating the characteristic points.By a density technique, we can demonstrate that this sequence(un)n∈ℕ converges to the solution of the system (25). By posing for the convenience of notation 0<α≤gf,n=h≤β≤1, the non-degenerate problem is reworded as.Findu such that
(29)∂u∂t+∂∂x(h(x)∂3u∂x3)=0in[0,T]×Ω,u(0,x)=u0inΩ,plusboundariesconditions.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].


## 3.1. Interpretation of the Diffusivity Action in the Diffusion Process

Between two consecutive maximaXmaxi and Xmaxi+1, the function ge+ is piecewise constant and can be written as follows:
(23)g(x)={0,atXmaxi,Xmaxi+1,49,atXinfi,Xinfi+1,Xmini,1,on]Xinfi,Xmini[∪]Xmini,Xinfi+1[,19,sur]Xmaxi,Xinfi[∪]Xinfi+1,Xmaxi+1[.We have a diffusion in the sense thatut=∂u/∂t unhook to the maximum of the curve u0. The reason is that the diffusion is more pronounced in local minimum because the smoothing effect tends to regularize curves. Thus, the upper envelope u(t≅∞,·) of u0 is less oscillating than u0 is. The same behavior occurs for the mean envelope. We would then try to interpret locally the relationship u≥u0 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 ∥∇u(t,·)∥L2(Ω)≤∥∇u0∥L2(Ω).In summary, we have just built a PDE system(1)
having the input signalu0 to decompose as an initial value condition,(2)
such that the result of a diffusion process preserves characteristic points ofu0,(3)
with a regularity similar to a cubic spline.


## 3.2. Formulation of the Mathematical Problem

LetA be one of the above constructed interpolation operators, and T>0 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 solutionu of the system
(24)∂u∂t+Au=0in[0,T]×Ω,u(0,x)=u0inΩ,plusboundariesconditions.
Then, the mean envelope (interpolating turning points) calculus is stated as follows.For givenu0, find u such that
(25)∂u∂t+∂∂x(gf(x)∂3u∂x3)=0in[0,T]×Ω,u(0,x)=u0inΩ,plusboundariesconditions.


## 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
(26)0<α≤gf,α(x)≤β≤1.
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(27)gf,n(x)=gf(x)+1n,n∈ℕ*
which verifies the conditions of (26).The new system that we call nondegenerate problem is fomulated as follows.Findun such that
(28)∂un∂t+∂∂x(gn(x)∂3un∂x3)=0in[0,T]×Ω,un(0,x)=u0inΩ,plusboundariesconditions.


## 3.5. Formulation of the Nondegenerate Problem

Forn sufficiently large, un is a close approximation of the envelope u, interpolating the characteristic points.By a density technique, we can demonstrate that this sequence(un)n∈ℕ converges to the solution of the system (25). By posing for the convenience of notation 0<α≤gf,n=h≤β≤1, the non-degenerate problem is reworded as.Findu such that
(29)∂u∂t+∂∂x(h(x)∂3u∂x3)=0in[0,T]×Ω,u(0,x)=u0inΩ,plusboundariesconditions.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 Section4.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 Section4.3.4.(2)
Second step: in Section5, construction of the sequence (un)n∈ℕ of solutions of problem (29), where un is the solution obtained with the diffusivity. hn(x)=h(x)+(1/n),  n∈ℕ*, see Section 5. Finally, demonstrate that the sequence (un)n∈ℕ converges to a limit u, where u 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.
Foru0∈H1(Ω), exists a unique solution u∈L2(0,T;𝒱)∩C([0,T];H1(Ω)) of the variational problem
(30)ddt∫Ω‍∂u∂x∂v∂xdx+∫Ω‍h(x)∂3u∂x3∂3v∂x3dx=0,foranyv∈𝒱
such that u(0)=u0 and du/dt∈L2(0,T;H).  Reversely, let u be a solution of the variational problem (30), then (∂u/∂t)+(∂/∂x)(h(x)(∂3u/∂x3))=0.And the complete follows result.Theorem 9.
Foru0∈H1(Ω), exists a unique solution
(31)u∈L2(0,T;𝒱)∩C([0,T];H1(Ω)),
verifying
(32)∂u∂t+∂∂x(h(x)∂3u∂x3)=0in[0,T]×Ω,∂3u∂x3|∂Ω=0(onboundariesofΩ),∂u∂n|∂Ω=0(onboundariesofΩ),u(0,x)=u0inΩ.


### 4.3. Variational Formulation

Letv∈H3(Ω), and (to make sense to integrals) multiply the equation
(33)∂u∂t+∂∂x(h(x)∂3u∂x3)=0
of the system (29) by the test function ∂2v/∂x2 after we integrate on Ω, then it comes out that
(34)∫Ω‍dudt∂2v∂x2dx+∫Ω‍∂∂x(h(x)∂3u∂x3)∂2v∂x2dx=0.
By using the Green formula (see Definition 4) and integration by parts, we get firstly
(35)∫Ω‍∂∂x(h(x)∂3u∂x3)∂2v∂x2dx=[h(x)∂3u∂x3∂2v∂x2]∂Ω-∫Ω‍h(x)∂3u∂x3∂3v∂x3dx,
and secondly,
(36)∫Ω‍dudt∂2v∂x2dx=ddt∫Ω‍u∂2v∂x2dx,∫Ω‍u∂2v∂x2dx=∫∂Ω‍u∂v∂ndσ-∫Ω‍∂u∂x∂v∂xdx.
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,
(37)∫∂Ω‍u∂v∂ndσ=∫∂Ω‍u∂v∂xdσ.
In the following, we adopt the notation ∂v/∂n.Consideringu∈H3(Ω) such that
(38)∂3u∂x3|∂Ω=0,v∈V={v∈H3(Ω)∣∂v∂n|∂Ω=0}
it comes out that
(39)ddt∫Ω‍∂u∂x∂v∂xdx+∫Ω‍h(x)∂3u∂x3∂3v∂x3dx=0,∀v∈V.
Let us consider the following bilinear forms
(40)a(u,v)=∫Ω‍h(x)∂3u∂x3∂3v∂x3dx,b(u,v)=∫Ω‍∂u∂x∂v∂xdx.
The addition of boundaries conditions transforms our initial problem into the system
(41)∂u∂t+∂∂x(h(x)∂3u∂x3)=0in[0,T]×Ω,∂3u∂x3|∂Ω=0(ontheedgeofΩ),∂u∂n|∂Ω=0(ontheedgeofΩ),u(0,x)=u0inΩ.The variational formulation of our problem consequently, subject to additional conditions that would cause a change of space, is the following(42)Findu∈Vsuchthatddt(b(u,v))+a(u,v)=0,∀v∈V,
where a and b 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 V.


#### 4.3.1. On the Quotient SpaceV/ℝ

The spaceV is defined by D(a)∩D(b), where D(a) and D(b) are the domains of a and b. Let us consider the quotient space
(43)Vℝ={∂v∂n|{v+c,c∈ℝ}suchthatv∈H3(Ω),and∂v∂n|∂Ω=0}
and the new space 𝒱  defined by
(44)𝒱={v∈H3(Ω)∣∂v∂n|∂Ω=0,∫Ωvdx=0}.
Let v~∈V/ℝ, and let u be a representative of v~ class, then (∂v~/∂n)|∂Ω is defined by
(45)∂v~∂n|∂Ω=∂(u+c)∂n|∂Ω=∂u∂n|∂Ω=0withc∈ℝ.
The nullity of ∫Ω‍v~dx means ∃c∈ℝ,  constant such that  ∫Ω‍(u+c)dx=0,  with  u  one representative of  v~(46)ifc=-1|Ω|∫Ω‍udx,then∫Ω‍v~dx=∫Ω‍(u+c)dx=∫Ω‍(u-1|Ω|∫Ω‍udx)dx=0.So for allv~∈V/ℝ,  we have  ∫Ω‍v~dx=0. Accordingly, any element of the quotient space V/ℝ is, on one hand, in H3(Ω) with normal derivative on the boundaries of Ω null and, on the other hand, is zero average on Ω. Thus, the quotient space V/ℝ is identified with the space 𝒱. We can work later in the space 𝒱 with zero mean. The final variational problem to solve is the following
(47)findu∈𝒱suchthatddt(b(u,v))+a(u,v)=0,∀v∈𝒱.


#### 4.3.2. Application of the Lions Theorem

To apply the Lions theorem, we need to identify the applicationsa,  b, and B and the spaces H and V. And finally, we must prove that the norm defined in 𝒱, which derives from the scalar product given by the bilinear application a, is equivalent to the norm of H3(Ω).Let us defineB as the identity, and let us define a and b by |a(u,v)|=(h(∂3u/∂x3),h(∂3v/∂x3))L2(Ω), and b(u,v)=∫Ω‍(∂u/∂x)(∂v/∂x)dx, obtained from (39), (40); the spaces V=𝒱 and H are define as follows.(a)The linear operatorB: in the variational problem (11) in the Lions theorem, (du/dt,Bv)H denotes a scalar product on H. Indeed, we define H={v∈H1(Ω)suchthat(∂v/∂n)|∂Ω=0,∫Ω‍vdx=0}.The operatorB=I, Bv=v, is the identity which is continuous (verify the Lions conditions Theorem 7), and the term (du/dt,Bv)H=(d/dt)(b(u,v)) is equal to
(48)∫Ω‍∂∂x(dudt)∂v∂xdx=∫Ω‍ddt(∂u∂x)∂v∂xdx=(∂∂x(dudt),∂v∂x)L2(Ω).
In this identification, ((∂/∂x)(du/dt),∂v/∂x)L2(Ω) is a well scalar product in the subspace H of H1(Ω), formed by functions with zero mean. Indeed, according to Poincare-Wirtinger inequality (see Definition 5) in H1(Ω),
(49)∥u-1|Ω|∫Ω‍udx∥H1(Ω)≤1λ∥∂u∂x∥L2(Ω),
where λ is the first positive eigenvalue of Laplace operator in Ω. But as u,v∈𝒱,  then  ∫Ω‍udx=0, consequently there exists C2>0 such that:
(50)∥u∥H1(Ω)≤C2∥∂u∂x∥L2(Ω),
that is to say the H1(Ω) norm is equivalent to the gradient norm in 𝒱. Hence, one can replace the scalar product of the gradient by H1(Ω) 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 toH3(Ω)norm.Theorem 10.
The space𝒱, resulting from (11) and 𝒱 analysis, defined by
(51)𝒱={v∈H3(Ω)∣∂v∂n|∂Ω=0,∫Ω‍vdx=0},
is an Hilbert space with the a norm equivalent to H3(Ω) norm.Proof.
At first𝒱={v∈H3(Ω)∣(∂v/∂n)∣∂Ω=0,∫Ω‍vdx=0} is a closed vector subspace of H3(Ω). Note that the condition of the nullity of the mean in 𝒱 may be replaced by ∫ΓΩ‍udx=0, where ΓΩ is a nonempty open subset of Ω. In this case 𝒱 becomes rather
(52)𝒱e={∫ΓΩvdx=0,and∂v∂n|∂Ωv∈H3(Ω)suchthat∃ΓΩ⊂Ωverifying∫ΓΩvdx=0,and∂v∂n|∂Ω=0}.
From the fact that(53)0<α≤h(x)≤β<1,∀x∈Ω,
the norms built with the two bilinear applications a and b are equivalent to the norm constructed by the linear form, noted by  |∥∥|𝒱2 with v∈𝒱, and given by
(54)|∥v∥|𝒱2=∥∂v∂x∥L2(Ω)2+∥∂3v∂x3∥L2(Ω)2.Indeed,|∥∥|𝒱2is as well a norm in𝒱. There is no difficulty to formally verify the triangular inequality. If |∥v∥|𝒱2=0, then ∥∂v/∂x∥L2(Ω)2=0 and ∥∂3v/∂x3∥L2(Ω)2=0. Consequently, v is necessarily constant almost everywhere in Ω, for example, v is equal to a constant real K. But as  ∫Ω‍vdx=0 implies K|Ω|=0, so K=v=0. 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 H3(Ω) norm, as it shown in Lemma A.2. Now, we just need to prove that the correspondence (u,v)∈𝒱×𝒱↦a(u,v) is continuous and coercive on the Hilbert space 𝒱.


#### 4.3.3. The Applicationa(·,·) Is Continuous and Coercive

Foru,v∈𝒱,
(55)|a(u,v)|=|(h∂3u∂x3,h∂3v∂x3)L2(Ω)|≤β∥∂3u∂x3∥L2(Ω)·∥∂3v∂x3∥L2(Ω)≤β|∥u∥|𝒱·|∥v∥|𝒱,
hence the continuity of a(·,·) is on 𝒱×𝒱. In order to complete the application conditions of Theorem 7, we only check the first condition of the same theorem.For(56)v∈𝒱,-∥v∥H1(Ω)2≤-∥∂v∂x∥L2(Ω)2,
hence
(57)|∥v∥|𝒱2-∥v∥H1(Ω)2≤|∥v∥|𝒱2-∥∂v∂x∥L2(Ω)2,
nevertheless
(58)a(v,v)=∥h∂3v∂x3∥L2(Ω)2≥α∥∂3v∂x3∥L2(Ω)2=α[|∥v∥|𝒱2-∥∂v∂x∥L2(Ω)2].
Thus,
(59)a(v,v)≥α[|∥v∥|𝒱2-∥v∥H1(Ω)2],∀v∈𝒱.
We have found thus (ρ,μ)=(α,α)such that
(60)a(v,v)≥ρ|∥v∥|𝒱2-μ∥v∥H1(Ω)2,∀v∈𝒱.
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.
Foru0∈H1(Ω), exists a unique solution u∈L2(0,T;𝒱)∩C([0,T];H1(Ω)) of the problem
(61)ddt∫Ω‍∂u∂x∂v∂xdx+∫Ω‍h(x)∂3u∂x3∂3v∂x3dx=0,foranyv∈𝒱suchthatu(0)=u0,dudt∈L2(0,T;H).


#### 4.3.4.  Equivalence with the Initial Problem

Reversely, letu be a solution of the variational problem
(62)ddt∫Ω‍∂u∂x∂v∂xdx+∫Ω‍h(x)∂3u∂x3∂3v∂x3dx=0∀v∈𝒱.
As we must consider a test function
(63)v∈H3(Ω)suchthat∂2v∂x2stillinD(Ω),∂v∂ν|∂Ω=0.Let us take a function testϕ∈D(Ω) and K∈ℝ such that
(64)∫Ω‍(ϕ+K)dx=0.We pose the following Neumann problem:(65)∂2v∂x2=ϕ+K,∂v∂n|∂Ω=0.Its variational formulation gives a unique solution inH={w∈H1(Ω)∣∫Ω‍wdx=0}, and the regularity of ϕ∈H2(Ω) implies the regularity of the solution which is in v∈H3(Ω), consequently v∈𝒱.Returning to (62) with v, solution of the problem (65), the first term can be written as
(66)ddt∫Ω‍∂u∂x∂v∂xdx=ddt(∫Ω‍u∂2v∂x2dx-∫∂Ω‍u∂v∂ndσ),as∂2v∂x2=ϕ+K,∂v∂n|∂Ω=0,anditcomesoutthatddt∫Ω‍∂u∂x∂v∂xdx=ddt∫Ω‍u(ϕ+K)dx=ddt(∫Ω‍uϕdx+K∫Ω‍udx).
However, ∫Ω‍udx=0, subsequently
(67)ddt∫Ω‍∂u∂x∂v∂xdx=ddt∫Ω‍uϕdx=∫Ω‍dudtϕdx.
The second term of (62) is written as follows:
(68)∫Ω‍h(x)∂3u∂x3∂3v∂x3dx=∫Ω‍h(x)∂3u∂x3∂ϕ∂xdx,
which gives in the distribution sense
(69)∫Ω‍h(x)∂3u∂x3∂3v∂x3dx=∫Ω‍∂∂x(h(x)∂3u∂x3)ϕdx,
from which
(70)∫Ω‍dudtϕdx+∫Ω‍∂∂x(h(x)∂3u∂x3)ϕdx=0,∀ϕ∈D(Ω).
We deduce that
(71)dudt+∂∂x(h(x)∂3u∂x3)=0.
We verify now that (∂3u/∂x3)|∂Ω=0. Returning to (70) with ϕ=∂2v/∂x2 and according to the Green formula, we have
(72)∫Ω‍[dudt+∂∂x(h(x)∂3u∂x3)]∂2v∂x2dx-∫∂Ω‍h(x)∂3u∂x3∂∂n(∂v∂x)·ndσ=0∀v∈𝒱.
But as (∂/∂n)(∂v/∂x)·n=(∂/∂x)(∂v/∂x)=∂2v/∂x2, it gives according to (71) ∫∂Ω‍h(x)(∂3u/∂x3)(∂2v/∂x2)dσ=0∀v∈𝒱, and consequently h(x)(∂3u/∂x3)=0 on boundaries of Ω, or as h(x)≠0,  ∂3u/∂x3=0  on boundaries of Ω.In conclusionu is solution of the initial equation (41), then we enunciate the following theorem.Theorem 12.
Foru0∈H1(Ω), exists a unique solution u∈L2(0,T;𝒱)∩C([0,T];H1(Ω)), verifying
(73)∂u∂t+∂∂x(h(x)∂3u∂x3)=0in[0,T]×Ω,∂3u∂x3|∂Ω=0(onboundariesofΩ),∂u∂n|∂Ω=0(onboundariesofΩ),u(0,x)=u0inΩ.


#### 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 u0 that is not necessarily zero mean. But at the cost of working with an initial condition v0∈H1(Ω), such that v0=u0-(1/|Ω|)∫Ω‍u0dx, we can consider H as H1(Ω) 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 𝒱e 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 𝒱e. The modeling space of the extrema interpolation is 𝒱e defined in (52).Letf be a function admitting at least three extrema we can consider that f passes at least once a zero at a point xz 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 f.


## 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 Section4.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 Section4.3.4.(2)
Second step: in Section5, construction of the sequence (un)n∈ℕ of solutions of problem (29), where un is the solution obtained with the diffusivity. hn(x)=h(x)+(1/n),  n∈ℕ*, see Section 5. Finally, demonstrate that the sequence (un)n∈ℕ converges to a limit u, where u 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.
Foru0∈H1(Ω), exists a unique solution u∈L2(0,T;𝒱)∩C([0,T];H1(Ω)) of the variational problem
(30)ddt∫Ω‍∂u∂x∂v∂xdx+∫Ω‍h(x)∂3u∂x3∂3v∂x3dx=0,foranyv∈𝒱
such that u(0)=u0 and du/dt∈L2(0,T;H).  Reversely, let u be a solution of the variational problem (30), then (∂u/∂t)+(∂/∂x)(h(x)(∂3u/∂x3))=0.And the complete follows result.Theorem 9.
Foru0∈H1(Ω), exists a unique solution
(31)u∈L2(0,T;𝒱)∩C([0,T];H1(Ω)),
verifying
(32)∂u∂t+∂∂x(h(x)∂3u∂x3)=0in[0,T]×Ω,∂3u∂x3|∂Ω=0(onboundariesofΩ),∂u∂n|∂Ω=0(onboundariesofΩ),u(0,x)=u0inΩ.


## 4.3. Variational Formulation

Letv∈H3(Ω), and (to make sense to integrals) multiply the equation
(33)∂u∂t+∂∂x(h(x)∂3u∂x3)=0
of the system (29) by the test function ∂2v/∂x2 after we integrate on Ω, then it comes out that
(34)∫Ω‍dudt∂2v∂x2dx+∫Ω‍∂∂x(h(x)∂3u∂x3)∂2v∂x2dx=0.
By using the Green formula (see Definition 4) and integration by parts, we get firstly
(35)∫Ω‍∂∂x(h(x)∂3u∂x3)∂2v∂x2dx=[h(x)∂3u∂x3∂2v∂x2]∂Ω-∫Ω‍h(x)∂3u∂x3∂3v∂x3dx,
and secondly,
(36)∫Ω‍dudt∂2v∂x2dx=ddt∫Ω‍u∂2v∂x2dx,∫Ω‍u∂2v∂x2dx=∫∂Ω‍u∂v∂ndσ-∫Ω‍∂u∂x∂v∂xdx.
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,
(37)∫∂Ω‍u∂v∂ndσ=∫∂Ω‍u∂v∂xdσ.
In the following, we adopt the notation ∂v/∂n.Consideringu∈H3(Ω) such that
(38)∂3u∂x3|∂Ω=0,v∈V={v∈H3(Ω)∣∂v∂n|∂Ω=0}
it comes out that
(39)ddt∫Ω‍∂u∂x∂v∂xdx+∫Ω‍h(x)∂3u∂x3∂3v∂x3dx=0,∀v∈V.
Let us consider the following bilinear forms
(40)a(u,v)=∫Ω‍h(x)∂3u∂x3∂3v∂x3dx,b(u,v)=∫Ω‍∂u∂x∂v∂xdx.
The addition of boundaries conditions transforms our initial problem into the system
(41)∂u∂t+∂∂x(h(x)∂3u∂x3)=0in[0,T]×Ω,∂3u∂x3|∂Ω=0(ontheedgeofΩ),∂u∂n|∂Ω=0(ontheedgeofΩ),u(0,x)=u0inΩ.The variational formulation of our problem consequently, subject to additional conditions that would cause a change of space, is the following(42)Findu∈Vsuchthatddt(b(u,v))+a(u,v)=0,∀v∈V,
where a and b 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 V.


### 4.3.1. On the Quotient SpaceV/ℝ

The spaceV is defined by D(a)∩D(b), where D(a) and D(b) are the domains of a and b. Let us consider the quotient space
(43)Vℝ={∂v∂n|{v+c,c∈ℝ}suchthatv∈H3(Ω),and∂v∂n|∂Ω=0}
and the new space 𝒱  defined by
(44)𝒱={v∈H3(Ω)∣∂v∂n|∂Ω=0,∫Ωvdx=0}.
Let v~∈V/ℝ, and let u be a representative of v~ class, then (∂v~/∂n)|∂Ω is defined by
(45)∂v~∂n|∂Ω=∂(u+c)∂n|∂Ω=∂u∂n|∂Ω=0withc∈ℝ.
The nullity of ∫Ω‍v~dx means ∃c∈ℝ,  constant such that  ∫Ω‍(u+c)dx=0,  with  u  one representative of  v~(46)ifc=-1|Ω|∫Ω‍udx,then∫Ω‍v~dx=∫Ω‍(u+c)dx=∫Ω‍(u-1|Ω|∫Ω‍udx)dx=0.So for allv~∈V/ℝ,  we have  ∫Ω‍v~dx=0. Accordingly, any element of the quotient space V/ℝ is, on one hand, in H3(Ω) with normal derivative on the boundaries of Ω null and, on the other hand, is zero average on Ω. Thus, the quotient space V/ℝ is identified with the space 𝒱. We can work later in the space 𝒱 with zero mean. The final variational problem to solve is the following
(47)findu∈𝒱suchthatddt(b(u,v))+a(u,v)=0,∀v∈𝒱.


### 4.3.2. Application of the Lions Theorem

To apply the Lions theorem, we need to identify the applicationsa,  b, and B and the spaces H and V. And finally, we must prove that the norm defined in 𝒱, which derives from the scalar product given by the bilinear application a, is equivalent to the norm of H3(Ω).Let us defineB as the identity, and let us define a and b by |a(u,v)|=(h(∂3u/∂x3),h(∂3v/∂x3))L2(Ω), and b(u,v)=∫Ω‍(∂u/∂x)(∂v/∂x)dx, obtained from (39), (40); the spaces V=𝒱 and H are define as follows.(a)The linear operatorB: in the variational problem (11) in the Lions theorem, (du/dt,Bv)H denotes a scalar product on H. Indeed, we define H={v∈H1(Ω)suchthat(∂v/∂n)|∂Ω=0,∫Ω‍vdx=0}.The operatorB=I, Bv=v, is the identity which is continuous (verify the Lions conditions Theorem 7), and the term (du/dt,Bv)H=(d/dt)(b(u,v)) is equal to
(48)∫Ω‍∂∂x(dudt)∂v∂xdx=∫Ω‍ddt(∂u∂x)∂v∂xdx=(∂∂x(dudt),∂v∂x)L2(Ω).
In this identification, ((∂/∂x)(du/dt),∂v/∂x)L2(Ω) is a well scalar product in the subspace H of H1(Ω), formed by functions with zero mean. Indeed, according to Poincare-Wirtinger inequality (see Definition 5) in H1(Ω),
(49)∥u-1|Ω|∫Ω‍udx∥H1(Ω)≤1λ∥∂u∂x∥L2(Ω),
where λ is the first positive eigenvalue of Laplace operator in Ω. But as u,v∈𝒱,  then  ∫Ω‍udx=0, consequently there exists C2>0 such that:
(50)∥u∥H1(Ω)≤C2∥∂u∂x∥L2(Ω),
that is to say the H1(Ω) norm is equivalent to the gradient norm in 𝒱. Hence, one can replace the scalar product of the gradient by H1(Ω) 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 toH3(Ω)norm.Theorem 10.
The space𝒱, resulting from (11) and 𝒱 analysis, defined by
(51)𝒱={v∈H3(Ω)∣∂v∂n|∂Ω=0,∫Ω‍vdx=0},
is an Hilbert space with the a norm equivalent to H3(Ω) norm.Proof.
At first𝒱={v∈H3(Ω)∣(∂v/∂n)∣∂Ω=0,∫Ω‍vdx=0} is a closed vector subspace of H3(Ω). Note that the condition of the nullity of the mean in 𝒱 may be replaced by ∫ΓΩ‍udx=0, where ΓΩ is a nonempty open subset of Ω. In this case 𝒱 becomes rather
(52)𝒱e={∫ΓΩvdx=0,and∂v∂n|∂Ωv∈H3(Ω)suchthat∃ΓΩ⊂Ωverifying∫ΓΩvdx=0,and∂v∂n|∂Ω=0}.
From the fact that(53)0<α≤h(x)≤β<1,∀x∈Ω,
the norms built with the two bilinear applications a and b are equivalent to the norm constructed by the linear form, noted by  |∥∥|𝒱2 with v∈𝒱, and given by
(54)|∥v∥|𝒱2=∥∂v∂x∥L2(Ω)2+∥∂3v∂x3∥L2(Ω)2.Indeed,|∥∥|𝒱2is as well a norm in𝒱. There is no difficulty to formally verify the triangular inequality. If |∥v∥|𝒱2=0, then ∥∂v/∂x∥L2(Ω)2=0 and ∥∂3v/∂x3∥L2(Ω)2=0. Consequently, v is necessarily constant almost everywhere in Ω, for example, v is equal to a constant real K. But as  ∫Ω‍vdx=0 implies K|Ω|=0, so K=v=0. 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 H3(Ω) norm, as it shown in Lemma A.2. Now, we just need to prove that the correspondence (u,v)∈𝒱×𝒱↦a(u,v) is continuous and coercive on the Hilbert space 𝒱.


### 4.3.3. The Applicationa(·,·) Is Continuous and Coercive

Foru,v∈𝒱,
(55)|a(u,v)|=|(h∂3u∂x3,h∂3v∂x3)L2(Ω)|≤β∥∂3u∂x3∥L2(Ω)·∥∂3v∂x3∥L2(Ω)≤β|∥u∥|𝒱·|∥v∥|𝒱,
hence the continuity of a(·,·) is on 𝒱×𝒱. In order to complete the application conditions of Theorem 7, we only check the first condition of the same theorem.For(56)v∈𝒱,-∥v∥H1(Ω)2≤-∥∂v∂x∥L2(Ω)2,
hence
(57)|∥v∥|𝒱2-∥v∥H1(Ω)2≤|∥v∥|𝒱2-∥∂v∂x∥L2(Ω)2,
nevertheless
(58)a(v,v)=∥h∂3v∂x3∥L2(Ω)2≥α∥∂3v∂x3∥L2(Ω)2=α[|∥v∥|𝒱2-∥∂v∂x∥L2(Ω)2].
Thus,
(59)a(v,v)≥α[|∥v∥|𝒱2-∥v∥H1(Ω)2],∀v∈𝒱.
We have found thus (ρ,μ)=(α,α)such that
(60)a(v,v)≥ρ|∥v∥|𝒱2-μ∥v∥H1(Ω)2,∀v∈𝒱.
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.
Foru0∈H1(Ω), exists a unique solution u∈L2(0,T;𝒱)∩C([0,T];H1(Ω)) of the problem
(61)ddt∫Ω‍∂u∂x∂v∂xdx+∫Ω‍h(x)∂3u∂x3∂3v∂x3dx=0,foranyv∈𝒱suchthatu(0)=u0,dudt∈L2(0,T;H).


### 4.3.4.  Equivalence with the Initial Problem

Reversely, letu be a solution of the variational problem
(62)ddt∫Ω‍∂u∂x∂v∂xdx+∫Ω‍h(x)∂3u∂x3∂3v∂x3dx=0∀v∈𝒱.
As we must consider a test function
(63)v∈H3(Ω)suchthat∂2v∂x2stillinD(Ω),∂v∂ν|∂Ω=0.Let us take a function testϕ∈D(Ω) and K∈ℝ such that
(64)∫Ω‍(ϕ+K)dx=0.We pose the following Neumann problem:(65)∂2v∂x2=ϕ+K,∂v∂n|∂Ω=0.Its variational formulation gives a unique solution inH={w∈H1(Ω)∣∫Ω‍wdx=0}, and the regularity of ϕ∈H2(Ω) implies the regularity of the solution which is in v∈H3(Ω), consequently v∈𝒱.Returning to (62) with v, solution of the problem (65), the first term can be written as
(66)ddt∫Ω‍∂u∂x∂v∂xdx=ddt(∫Ω‍u∂2v∂x2dx-∫∂Ω‍u∂v∂ndσ),as∂2v∂x2=ϕ+K,∂v∂n|∂Ω=0,anditcomesoutthatddt∫Ω‍∂u∂x∂v∂xdx=ddt∫Ω‍u(ϕ+K)dx=ddt(∫Ω‍uϕdx+K∫Ω‍udx).
However, ∫Ω‍udx=0, subsequently
(67)ddt∫Ω‍∂u∂x∂v∂xdx=ddt∫Ω‍uϕdx=∫Ω‍dudtϕdx.
The second term of (62) is written as follows:
(68)∫Ω‍h(x)∂3u∂x3∂3v∂x3dx=∫Ω‍h(x)∂3u∂x3∂ϕ∂xdx,
which gives in the distribution sense
(69)∫Ω‍h(x)∂3u∂x3∂3v∂x3dx=∫Ω‍∂∂x(h(x)∂3u∂x3)ϕdx,
from which
(70)∫Ω‍dudtϕdx+∫Ω‍∂∂x(h(x)∂3u∂x3)ϕdx=0,∀ϕ∈D(Ω).
We deduce that
(71)dudt+∂∂x(h(x)∂3u∂x3)=0.
We verify now that (∂3u/∂x3)|∂Ω=0. Returning to (70) with ϕ=∂2v/∂x2 and according to the Green formula, we have
(72)∫Ω‍[dudt+∂∂x(h(x)∂3u∂x3)]∂2v∂x2dx-∫∂Ω‍h(x)∂3u∂x3∂∂n(∂v∂x)·ndσ=0∀v∈𝒱.
But as (∂/∂n)(∂v/∂x)·n=(∂/∂x)(∂v/∂x)=∂2v/∂x2, it gives according to (71) ∫∂Ω‍h(x)(∂3u/∂x3)(∂2v/∂x2)dσ=0∀v∈𝒱, and consequently h(x)(∂3u/∂x3)=0 on boundaries of Ω, or as h(x)≠0,  ∂3u/∂x3=0  on boundaries of Ω.In conclusionu is solution of the initial equation (41), then we enunciate the following theorem.Theorem 12.
Foru0∈H1(Ω), exists a unique solution u∈L2(0,T;𝒱)∩C([0,T];H1(Ω)), verifying
(73)∂u∂t+∂∂x(h(x)∂3u∂x3)=0in[0,T]×Ω,∂3u∂x3|∂Ω=0(onboundariesofΩ),∂u∂n|∂Ω=0(onboundariesofΩ),u(0,x)=u0inΩ.


### 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 u0 that is not necessarily zero mean. But at the cost of working with an initial condition v0∈H1(Ω), such that v0=u0-(1/|Ω|)∫Ω‍u0dx, we can consider H as H1(Ω) 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 𝒱e 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 𝒱e. The modeling space of the extrema interpolation is 𝒱e defined in (52).Letf be a function admitting at least three extrema we can consider that f passes at least once a zero at a point xz 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 f.


## 4.3.1. On the Quotient SpaceV/ℝ

The spaceV is defined by D(a)∩D(b), where D(a) and D(b) are the domains of a and b. Let us consider the quotient space
(43)Vℝ={∂v∂n|{v+c,c∈ℝ}suchthatv∈H3(Ω),and∂v∂n|∂Ω=0}
and the new space 𝒱  defined by
(44)𝒱={v∈H3(Ω)∣∂v∂n|∂Ω=0,∫Ωvdx=0}.
Let v~∈V/ℝ, and let u be a representative of v~ class, then (∂v~/∂n)|∂Ω is defined by
(45)∂v~∂n|∂Ω=∂(u+c)∂n|∂Ω=∂u∂n|∂Ω=0withc∈ℝ.
The nullity of ∫Ω‍v~dx means ∃c∈ℝ,  constant such that  ∫Ω‍(u+c)dx=0,  with  u  one representative of  v~(46)ifc=-1|Ω|∫Ω‍udx,then∫Ω‍v~dx=∫Ω‍(u+c)dx=∫Ω‍(u-1|Ω|∫Ω‍udx)dx=0.So for allv~∈V/ℝ,  we have  ∫Ω‍v~dx=0. Accordingly, any element of the quotient space V/ℝ is, on one hand, in H3(Ω) with normal derivative on the boundaries of Ω null and, on the other hand, is zero average on Ω. Thus, the quotient space V/ℝ is identified with the space 𝒱. We can work later in the space 𝒱 with zero mean. The final variational problem to solve is the following
(47)findu∈𝒱suchthatddt(b(u,v))+a(u,v)=0,∀v∈𝒱.


## 4.3.2. Application of the Lions Theorem

To apply the Lions theorem, we need to identify the applicationsa,  b, and B and the spaces H and V. And finally, we must prove that the norm defined in 𝒱, which derives from the scalar product given by the bilinear application a, is equivalent to the norm of H3(Ω).Let us defineB as the identity, and let us define a and b by |a(u,v)|=(h(∂3u/∂x3),h(∂3v/∂x3))L2(Ω), and b(u,v)=∫Ω‍(∂u/∂x)(∂v/∂x)dx, obtained from (39), (40); the spaces V=𝒱 and H are define as follows.(a)The linear operatorB: in the variational problem (11) in the Lions theorem, (du/dt,Bv)H denotes a scalar product on H. Indeed, we define H={v∈H1(Ω)suchthat(∂v/∂n)|∂Ω=0,∫Ω‍vdx=0}.The operatorB=I, Bv=v, is the identity which is continuous (verify the Lions conditions Theorem 7), and the term (du/dt,Bv)H=(d/dt)(b(u,v)) is equal to
(48)∫Ω‍∂∂x(dudt)∂v∂xdx=∫Ω‍ddt(∂u∂x)∂v∂xdx=(∂∂x(dudt),∂v∂x)L2(Ω).
In this identification, ((∂/∂x)(du/dt),∂v/∂x)L2(Ω) is a well scalar product in the subspace H of H1(Ω), formed by functions with zero mean. Indeed, according to Poincare-Wirtinger inequality (see Definition 5) in H1(Ω),
(49)∥u-1|Ω|∫Ω‍udx∥H1(Ω)≤1λ∥∂u∂x∥L2(Ω),
where λ is the first positive eigenvalue of Laplace operator in Ω. But as u,v∈𝒱,  then  ∫Ω‍udx=0, consequently there exists C2>0 such that:
(50)∥u∥H1(Ω)≤C2∥∂u∂x∥L2(Ω),
that is to say the H1(Ω) norm is equivalent to the gradient norm in 𝒱. Hence, one can replace the scalar product of the gradient by H1(Ω) 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 toH3(Ω)norm.Theorem 10.
The space𝒱, resulting from (11) and 𝒱 analysis, defined by
(51)𝒱={v∈H3(Ω)∣∂v∂n|∂Ω=0,∫Ω‍vdx=0},
is an Hilbert space with the a norm equivalent to H3(Ω) norm.Proof.
At first𝒱={v∈H3(Ω)∣(∂v/∂n)∣∂Ω=0,∫Ω‍vdx=0} is a closed vector subspace of H3(Ω). Note that the condition of the nullity of the mean in 𝒱 may be replaced by ∫ΓΩ‍udx=0, where ΓΩ is a nonempty open subset of Ω. In this case 𝒱 becomes rather
(52)𝒱e={∫ΓΩvdx=0,and∂v∂n|∂Ωv∈H3(Ω)suchthat∃ΓΩ⊂Ωverifying∫ΓΩvdx=0,and∂v∂n|∂Ω=0}.
From the fact that(53)0<α≤h(x)≤β<1,∀x∈Ω,
the norms built with the two bilinear applications a and b are equivalent to the norm constructed by the linear form, noted by  |∥∥|𝒱2 with v∈𝒱, and given by
(54)|∥v∥|𝒱2=∥∂v∂x∥L2(Ω)2+∥∂3v∂x3∥L2(Ω)2.Indeed,|∥∥|𝒱2is as well a norm in𝒱. There is no difficulty to formally verify the triangular inequality. If |∥v∥|𝒱2=0, then ∥∂v/∂x∥L2(Ω)2=0 and ∥∂3v/∂x3∥L2(Ω)2=0. Consequently, v is necessarily constant almost everywhere in Ω, for example, v is equal to a constant real K. But as  ∫Ω‍vdx=0 implies K|Ω|=0, so K=v=0. 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 H3(Ω) norm, as it shown in Lemma A.2. Now, we just need to prove that the correspondence (u,v)∈𝒱×𝒱↦a(u,v) is continuous and coercive on the Hilbert space 𝒱.


## 4.3.3. The Applicationa(·,·) Is Continuous and Coercive

Foru,v∈𝒱,
(55)|a(u,v)|=|(h∂3u∂x3,h∂3v∂x3)L2(Ω)|≤β∥∂3u∂x3∥L2(Ω)·∥∂3v∂x3∥L2(Ω)≤β|∥u∥|𝒱·|∥v∥|𝒱,
hence the continuity of a(·,·) is on 𝒱×𝒱. In order to complete the application conditions of Theorem 7, we only check the first condition of the same theorem.For(56)v∈𝒱,-∥v∥H1(Ω)2≤-∥∂v∂x∥L2(Ω)2,
hence
(57)|∥v∥|𝒱2-∥v∥H1(Ω)2≤|∥v∥|𝒱2-∥∂v∂x∥L2(Ω)2,
nevertheless
(58)a(v,v)=∥h∂3v∂x3∥L2(Ω)2≥α∥∂3v∂x3∥L2(Ω)2=α[|∥v∥|𝒱2-∥∂v∂x∥L2(Ω)2].
Thus,
(59)a(v,v)≥α[|∥v∥|𝒱2-∥v∥H1(Ω)2],∀v∈𝒱.
We have found thus (ρ,μ)=(α,α)such that
(60)a(v,v)≥ρ|∥v∥|𝒱2-μ∥v∥H1(Ω)2,∀v∈𝒱.
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.
Foru0∈H1(Ω), exists a unique solution u∈L2(0,T;𝒱)∩C([0,T];H1(Ω)) of the problem
(61)ddt∫Ω‍∂u∂x∂v∂xdx+∫Ω‍h(x)∂3u∂x3∂3v∂x3dx=0,foranyv∈𝒱suchthatu(0)=u0,dudt∈L2(0,T;H).


## 4.3.4.  Equivalence with the Initial Problem

Reversely, letu be a solution of the variational problem
(62)ddt∫Ω‍∂u∂x∂v∂xdx+∫Ω‍h(x)∂3u∂x3∂3v∂x3dx=0∀v∈𝒱.
As we must consider a test function
(63)v∈H3(Ω)suchthat∂2v∂x2stillinD(Ω),∂v∂ν|∂Ω=0.Let us take a function testϕ∈D(Ω) and K∈ℝ such that
(64)∫Ω‍(ϕ+K)dx=0.We pose the following Neumann problem:(65)∂2v∂x2=ϕ+K,∂v∂n|∂Ω=0.Its variational formulation gives a unique solution inH={w∈H1(Ω)∣∫Ω‍wdx=0}, and the regularity of ϕ∈H2(Ω) implies the regularity of the solution which is in v∈H3(Ω), consequently v∈𝒱.Returning to (62) with v, solution of the problem (65), the first term can be written as
(66)ddt∫Ω‍∂u∂x∂v∂xdx=ddt(∫Ω‍u∂2v∂x2dx-∫∂Ω‍u∂v∂ndσ),as∂2v∂x2=ϕ+K,∂v∂n|∂Ω=0,anditcomesoutthatddt∫Ω‍∂u∂x∂v∂xdx=ddt∫Ω‍u(ϕ+K)dx=ddt(∫Ω‍uϕdx+K∫Ω‍udx).
However, ∫Ω‍udx=0, subsequently
(67)ddt∫Ω‍∂u∂x∂v∂xdx=ddt∫Ω‍uϕdx=∫Ω‍dudtϕdx.
The second term of (62) is written as follows:
(68)∫Ω‍h(x)∂3u∂x3∂3v∂x3dx=∫Ω‍h(x)∂3u∂x3∂ϕ∂xdx,
which gives in the distribution sense
(69)∫Ω‍h(x)∂3u∂x3∂3v∂x3dx=∫Ω‍∂∂x(h(x)∂3u∂x3)ϕdx,
from which
(70)∫Ω‍dudtϕdx+∫Ω‍∂∂x(h(x)∂3u∂x3)ϕdx=0,∀ϕ∈D(Ω).
We deduce that
(71)dudt+∂∂x(h(x)∂3u∂x3)=0.
We verify now that (∂3u/∂x3)|∂Ω=0. Returning to (70) with ϕ=∂2v/∂x2 and according to the Green formula, we have
(72)∫Ω‍[dudt+∂∂x(h(x)∂3u∂x3)]∂2v∂x2dx-∫∂Ω‍h(x)∂3u∂x3∂∂n(∂v∂x)·ndσ=0∀v∈𝒱.
But as (∂/∂n)(∂v/∂x)·n=(∂/∂x)(∂v/∂x)=∂2v/∂x2, it gives according to (71) ∫∂Ω‍h(x)(∂3u/∂x3)(∂2v/∂x2)dσ=0∀v∈𝒱, and consequently h(x)(∂3u/∂x3)=0 on boundaries of Ω, or as h(x)≠0,  ∂3u/∂x3=0  on boundaries of Ω.In conclusionu is solution of the initial equation (41), then we enunciate the following theorem.Theorem 12.
Foru0∈H1(Ω), exists a unique solution u∈L2(0,T;𝒱)∩C([0,T];H1(Ω)), verifying
(73)∂u∂t+∂∂x(h(x)∂3u∂x3)=0in[0,T]×Ω,∂3u∂x3|∂Ω=0(onboundariesofΩ),∂u∂n|∂Ω=0(onboundariesofΩ),u(0,x)=u0inΩ.


## 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 u0 that is not necessarily zero mean. But at the cost of working with an initial condition v0∈H1(Ω), such that v0=u0-(1/|Ω|)∫Ω‍u0dx, we can consider H as H1(Ω) 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 𝒱e 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 𝒱e. The modeling space of the extrema interpolation is 𝒱e defined in (52).Letf be a function admitting at least three extrema we can consider that f passes at least once a zero at a point xz 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 f.


## 5.  Convergence of the Sequence(un)n∈ℕ of Solutions of (41) to a Solution of the Called Degenerate Problem

In the following sections, we demonstrate that the sequence(un)n∈ℕ of solutions of the non degenerate problem is bounded. Next, we prove that there exists a subsequence of (un)n∈ℕ which converges weakly to an element u∈𝒱 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 ∂2un/∂x2):
(74)∫Ω‍∂un∂t∂2un∂x2dx+∫Ω‍∂∂x(hn(x)∂3un∂x3)∂2un∂x2dx=0,
which is equivalent to
(75)〈∂un∂t,∂2un∂x2〉V′V+〈∂∂x(hn(x)∂3un∂x3),∂2un∂x2〉V′V=0.
After integration by parts, it comes out that
(76)12∫Ω‍∂∂t(|∂un∂x|2)dx+∫Ω‍hn(x)(∂3un∂x3)2dx=0,
then by integrating with respect to the variable t∈]0,T[, it comes out that
(77)12∥∂un(t,·)∂x∥L2(Ω)2-12∥∂u0∂x∥L2(Ω)2+∫0t‍∫Ω‍hn(x)(∂3un∂x3)2dxds=0,
which implies that
(78)∥∂un(t,·)∂x∥L2(Ω)2+2∫0t‍∫Ω‍hn(x)(∂3un∂x3)2dxds=∥∂u0∂x∥L2(Ω)2.
But as hn(x)≥α, we have
(79)∥∂un(t,·)∂x∥L2(Ω)2+2∫0t‍∫Ω‍hn(x)(∂3un∂x3)2dxds≥∥∂un(t,·)∂x∥L2(Ω)2+2α∫0t‍∫Ω‍(∂3un∂x3)2dxds
or
(80)∥∂un(t,·)∂x∥L2(Ω)2+2∫0t‍∫Ω‍hn(x)(∂3un∂x3)2dxds≥∥∂un(t,·)∂x∥L2(Ω)2+2α∫0t‍∥∂3un∂x3∥L2(Ω)2ds,
from which
(81)∥∂u0∂x∥L2(Ω)2≥∥∂un(t,·)∂x∥L2(Ω)2+2α∫0t‍[|∥un∥|𝒱2-∥∂un(s,·)∂x∥L2(Ω)2]ds.
Thus, on one hand,
(82)∥∂u0∂x∥L2(Ω)2≥∥∂un(t,·)∂x∥L2(Ω).2
But
(83)∥∂un(t,·)∂x∥L2(Ω)2≥λ2∥un∥H1(Ω)2,
where λ comes from inequality (49), which implies that
(84)∥un∥L∞(0,T;H1(Ω))2≤1λ2∥∂u0∂x∥L2(Ω)2,(85)∥un∥L2(0,T;H1(Ω))2≤Tλ2∥∂u0∂x∥L2(Ω)2.Then, the sequence(un)n∈ℕ  is bounded inL∞(0,T;H1(Ω)) and in L2(0,T;H1(Ω)). And on the other hand
(86)2α∫0T‍[|∥un∥|𝒱2-∥∂un(s,·)∂x∥L2(Ω)2]ds≤∥∂u0∂x∥L2(Ω)2,
which gives
(87)∫0T‍|∥un∥|𝒱2ds≤∫0T‍∥∂un(s,·)∂x∥L2(Ω)2ds+12α∥∂u0∂x∥L2(Ω)2.
But as
(88)∥∂un(s,·)∂x∥L2(Ω)2≤∥un(s,·)∥H1(Ω)2,
therefore
(89)∫0T‍|∥un∥|𝒱2ds≤∫0T‍∥un(s,·)∥H1(Ω)2ds+12α∥∂u0∂x∥L2(Ω)2.
Thus,
(90)∫0T‍|∥un∥|𝒱2ds≤T∥un∥L∞(0,T;H1(Ω))2+12α∥∂u0∂x∥L2(Ω)2,
which implies according to (84) that
(91)∫0T‍|∥un∥|𝒱2ds≤Tλ2∥∂u0∂x∥L2(Ω)2+12α∥∂u0∂x∥L2(Ω)2.
We finally deduced the second estimate
(92)∥un∥L2(0,T;𝒱)2=∥un∥L2(0,T;𝒱e)2≤(Tλ2+12α)∥∂u0∂x∥L2(Ω)2.Then, the sequence(un)n∈ℕis bounded inL2(0,T;𝒱).


### 5.2. Weak Convergence of  the Subsequence of(un)n∈ℕ

The sequence(un)n∈ℕ is bounded in the space L2(0,T;𝒱). Then there exists a subsequence (unk)k∈ℕ that converge weakly to u in L2(0,T;𝒱) with du/dt∈L2(0,T;V′).Letingv∈𝒱 be a test function, we have the following.Firstly. As unk→u weakly in L2(0,T;V), therefore, unk→u weakly in L2(0,T;H1(Ω)).Thus, in the distributions sense,(93)∫Ω‍dunkdt∂2v∂x2dx=ddt∫Ω‍∂unk∂x∂v∂xdx,
by integrating on [0,t], we have
(94)∫0t‍∫Ω‍dunkds∂2v∂x2dxds=-∫0t‍dds∫Ω‍∂unk∂x∂v∂xdxds=-∫0t‍〈dunkds,∂2v∂x2〉V′Vds=-∫Ω‍(unk(t)-unk(0))∂2v∂x2dx.
By integrating by parts again the last equation, we have
(95)∫Ω‍(unk(t)-unk(0))∂2v∂x2dx=∫Ω‍∂unk(t)∂x∂v∂xdx-∫Ω‍∂unk(0)∂x∂v∂xdx.
This expression converges to
(96)∫Ω‍∂u(t)∂x∂v∂xdx-∫Ω‍∂u0∂x∂v∂xdx=∫Ω‍[∂u(t)∂x-∂u0∂x]∂v∂xdx=∫Ω‍[u(t)-u0]∂2v∂x2dx=∫0t‍〈duds,∂2v∂x2〉V′Vds.
Therefore,∫0t‍∫Ω‍(dunk/ds)(∂2v/∂x2)dxdsconverges when nk  tends to infinity to
(97)∫0t‍∫Ω‍duds∂2v∂x2dxds.
Consequently, ∂unk/∂t→∂u/∂t weakly in L2(0,T;𝒱′).Secondly. As hn(x)=h(x)+(1/n),(98)wehave∂∂x(hnk∂3unk∂x3-h(x)∂3u∂x3)=∂∂x(h(x)∂3∂x3(unk-u))+1nk∂4unk∂x4.However,∫Ω‍∂4unk∂x4vdx=∫∂Ω‍∂3u∂x3v·ndσ-∫Ω‍∂v∂x∂3unk∂x3dx.Andbecause∂3u∂x3|∂Ω=0,∂3unk∂x3convergesto∂3u∂x3inL2(0,T;L2(Ω)),
then,
(99)1nk|∫Ω‍∂4unk∂x4vdx|=1nk|∫Ω‍∂3unk∂x3∂v∂xdx|≤1nk∥∂v∂x∥L2(Ω)∥∂3unk∂x3∥L2(Ω)
converges to zero when nk tends to infinity.On the other hand, from (98), we have
(100)|∫Ω‍∂∂x(h(x)∂3∂x3(unk-u))vdx|≤β|∫Ω‍∂3∂x3(unk-u)∂v∂xdx|
which converges to zero in L2(0,T;L2(Ω)) according to weak convergence (see Definition 3) in H3(Ω)  of  ∂3unk/∂x3  to  ∂3u/∂x3. Thus, we have shown that in the distributions sense of the one hand
(101)∫Ω‍dunkdt∂2v∂x2dx
converges to
(102)∫Ω‍dudt∂2v∂x2dx,
and on the other hand,
(103)∫Ω‍∂∂x(hnk∂3unk∂x3)vdx
converges to
(104)∫Ω‍∂∂x(h(x)∂3u∂x3)vdx.
Consequently, for each v∈L2(0,T;V′),
(105)∫Ω‍dunkdtvdx+∫Ω‍∂∂x(hnk∂3unk∂x3)vdx
converges in the distributions sense to
(106)∫Ω‍dudtvdx+∫Ω‍∂∂x(h(x)∂3u∂x3)vdx,
that is to say, 〈(dunk/dt),v〉V′V+〈(∂/∂x)(hnk(x)(∂3unk/∂x3)),v〉V′V converges when nk  tends to infinity to 〈(du/dt),v〉V′V+〈(∂/∂x)(h(x)(∂3u/∂x3)),v〉V′V. It remains to show that the limit u is a solution of the degenerate problem. But ∫Ω‍(dunk/dt)vdx+∫Ω‍(∂/∂x)(hnk(∂3unk/∂x3))vdx=0, which implies in the distributions sense ∫Ω‍(du/dt)vdx+∫Ω‍(∂/∂x)(h(x)(∂3u/∂x3))vdx=0. In other words, the weak limit of the sequence (unk)k∈ℕ is weak solution of the initial degenerate problem. Moreover, this solution is unique and is in L2(0,T;𝒱)∩C([0,T];H1(Ω)). Better, using the compact inclusion of 𝒱 (which behaves like H3(Ω) because of the equivalence of the norms) in H2(Ω), we conclude that the sequence of solutions of non-degenerate problem (approximated problem) converges strongly in H2(Ω) to the solution of the initial degenerate problem. The space H2(Ω) is, moreover, the space of strong solutions of the initial system (25).


## 5.1.  Some Estimations

From (28) we have in the distributions sense (multiplying by ∂2un/∂x2):
(74)∫Ω‍∂un∂t∂2un∂x2dx+∫Ω‍∂∂x(hn(x)∂3un∂x3)∂2un∂x2dx=0,
which is equivalent to
(75)〈∂un∂t,∂2un∂x2〉V′V+〈∂∂x(hn(x)∂3un∂x3),∂2un∂x2〉V′V=0.
After integration by parts, it comes out that
(76)12∫Ω‍∂∂t(|∂un∂x|2)dx+∫Ω‍hn(x)(∂3un∂x3)2dx=0,
then by integrating with respect to the variable t∈]0,T[, it comes out that
(77)12∥∂un(t,·)∂x∥L2(Ω)2-12∥∂u0∂x∥L2(Ω)2+∫0t‍∫Ω‍hn(x)(∂3un∂x3)2dxds=0,
which implies that
(78)∥∂un(t,·)∂x∥L2(Ω)2+2∫0t‍∫Ω‍hn(x)(∂3un∂x3)2dxds=∥∂u0∂x∥L2(Ω)2.
But as hn(x)≥α, we have
(79)∥∂un(t,·)∂x∥L2(Ω)2+2∫0t‍∫Ω‍hn(x)(∂3un∂x3)2dxds≥∥∂un(t,·)∂x∥L2(Ω)2+2α∫0t‍∫Ω‍(∂3un∂x3)2dxds
or
(80)∥∂un(t,·)∂x∥L2(Ω)2+2∫0t‍∫Ω‍hn(x)(∂3un∂x3)2dxds≥∥∂un(t,·)∂x∥L2(Ω)2+2α∫0t‍∥∂3un∂x3∥L2(Ω)2ds,
from which
(81)∥∂u0∂x∥L2(Ω)2≥∥∂un(t,·)∂x∥L2(Ω)2+2α∫0t‍[|∥un∥|𝒱2-∥∂un(s,·)∂x∥L2(Ω)2]ds.
Thus, on one hand,
(82)∥∂u0∂x∥L2(Ω)2≥∥∂un(t,·)∂x∥L2(Ω).2
But
(83)∥∂un(t,·)∂x∥L2(Ω)2≥λ2∥un∥H1(Ω)2,
where λ comes from inequality (49), which implies that
(84)∥un∥L∞(0,T;H1(Ω))2≤1λ2∥∂u0∂x∥L2(Ω)2,(85)∥un∥L2(0,T;H1(Ω))2≤Tλ2∥∂u0∂x∥L2(Ω)2.Then, the sequence(un)n∈ℕ  is bounded inL∞(0,T;H1(Ω)) and in L2(0,T;H1(Ω)). And on the other hand
(86)2α∫0T‍[|∥un∥|𝒱2-∥∂un(s,·)∂x∥L2(Ω)2]ds≤∥∂u0∂x∥L2(Ω)2,
which gives
(87)∫0T‍|∥un∥|𝒱2ds≤∫0T‍∥∂un(s,·)∂x∥L2(Ω)2ds+12α∥∂u0∂x∥L2(Ω)2.
But as
(88)∥∂un(s,·)∂x∥L2(Ω)2≤∥un(s,·)∥H1(Ω)2,
therefore
(89)∫0T‍|∥un∥|𝒱2ds≤∫0T‍∥un(s,·)∥H1(Ω)2ds+12α∥∂u0∂x∥L2(Ω)2.
Thus,
(90)∫0T‍|∥un∥|𝒱2ds≤T∥un∥L∞(0,T;H1(Ω))2+12α∥∂u0∂x∥L2(Ω)2,
which implies according to (84) that
(91)∫0T‍|∥un∥|𝒱2ds≤Tλ2∥∂u0∂x∥L2(Ω)2+12α∥∂u0∂x∥L2(Ω)2.
We finally deduced the second estimate
(92)∥un∥L2(0,T;𝒱)2=∥un∥L2(0,T;𝒱e)2≤(Tλ2+12α)∥∂u0∂x∥L2(Ω)2.Then, the sequence(un)n∈ℕis bounded inL2(0,T;𝒱).


## 5.2. Weak Convergence of  the Subsequence of(un)n∈ℕ

The sequence(un)n∈ℕ is bounded in the space L2(0,T;𝒱). Then there exists a subsequence (unk)k∈ℕ that converge weakly to u in L2(0,T;𝒱) with du/dt∈L2(0,T;V′).Letingv∈𝒱 be a test function, we have the following.Firstly. As unk→u weakly in L2(0,T;V), therefore, unk→u weakly in L2(0,T;H1(Ω)).Thus, in the distributions sense,(93)∫Ω‍dunkdt∂2v∂x2dx=ddt∫Ω‍∂unk∂x∂v∂xdx,
by integrating on [0,t], we have
(94)∫0t‍∫Ω‍dunkds∂2v∂x2dxds=-∫0t‍dds∫Ω‍∂unk∂x∂v∂xdxds=-∫0t‍〈dunkds,∂2v∂x2〉V′Vds=-∫Ω‍(unk(t)-unk(0))∂2v∂x2dx.
By integrating by parts again the last equation, we have
(95)∫Ω‍(unk(t)-unk(0))∂2v∂x2dx=∫Ω‍∂unk(t)∂x∂v∂xdx-∫Ω‍∂unk(0)∂x∂v∂xdx.
This expression converges to
(96)∫Ω‍∂u(t)∂x∂v∂xdx-∫Ω‍∂u0∂x∂v∂xdx=∫Ω‍[∂u(t)∂x-∂u0∂x]∂v∂xdx=∫Ω‍[u(t)-u0]∂2v∂x2dx=∫0t‍〈duds,∂2v∂x2〉V′Vds.
Therefore,∫0t‍∫Ω‍(dunk/ds)(∂2v/∂x2)dxdsconverges when nk  tends to infinity to
(97)∫0t‍∫Ω‍duds∂2v∂x2dxds.
Consequently, ∂unk/∂t→∂u/∂t weakly in L2(0,T;𝒱′).Secondly. As hn(x)=h(x)+(1/n),(98)wehave∂∂x(hnk∂3unk∂x3-h(x)∂3u∂x3)=∂∂x(h(x)∂3∂x3(unk-u))+1nk∂4unk∂x4.However,∫Ω‍∂4unk∂x4vdx=∫∂Ω‍∂3u∂x3v·ndσ-∫Ω‍∂v∂x∂3unk∂x3dx.Andbecause∂3u∂x3|∂Ω=0,∂3unk∂x3convergesto∂3u∂x3inL2(0,T;L2(Ω)),
then,
(99)1nk|∫Ω‍∂4unk∂x4vdx|=1nk|∫Ω‍∂3unk∂x3∂v∂xdx|≤1nk∥∂v∂x∥L2(Ω)∥∂3unk∂x3∥L2(Ω)
converges to zero when nk tends to infinity.On the other hand, from (98), we have
(100)|∫Ω‍∂∂x(h(x)∂3∂x3(unk-u))vdx|≤β|∫Ω‍∂3∂x3(unk-u)∂v∂xdx|
which converges to zero in L2(0,T;L2(Ω)) according to weak convergence (see Definition 3) in H3(Ω)  of  ∂3unk/∂x3  to  ∂3u/∂x3. Thus, we have shown that in the distributions sense of the one hand
(101)∫Ω‍dunkdt∂2v∂x2dx
converges to
(102)∫Ω‍dudt∂2v∂x2dx,
and on the other hand,
(103)∫Ω‍∂∂x(hnk∂3unk∂x3)vdx
converges to
(104)∫Ω‍∂∂x(h(x)∂3u∂x3)vdx.
Consequently, for each v∈L2(0,T;V′),
(105)∫Ω‍dunkdtvdx+∫Ω‍∂∂x(hnk∂3unk∂x3)vdx
converges in the distributions sense to
(106)∫Ω‍dudtvdx+∫Ω‍∂∂x(h(x)∂3u∂x3)vdx,
that is to say, 〈(dunk/dt),v〉V′V+〈(∂/∂x)(hnk(x)(∂3unk/∂x3)),v〉V′V converges when nk  tends to infinity to 〈(du/dt),v〉V′V+〈(∂/∂x)(h(x)(∂3u/∂x3)),v〉V′V. It remains to show that the limit u is a solution of the degenerate problem. But ∫Ω‍(dunk/dt)vdx+∫Ω‍(∂/∂x)(hnk(∂3unk/∂x3))vdx=0, which implies in the distributions sense ∫Ω‍(du/dt)vdx+∫Ω‍(∂/∂x)(h(x)(∂3u/∂x3))vdx=0. In other words, the weak limit of the sequence (unk)k∈ℕ is weak solution of the initial degenerate problem. Moreover, this solution is unique and is in L2(0,T;𝒱)∩C([0,T];H1(Ω)). Better, using the compact inclusion of 𝒱 (which behaves like H3(Ω) because of the equivalence of the norms) in H2(Ω), we conclude that the sequence of solutions of non-degenerate problem (approximated problem) converges strongly in H2(Ω) to the solution of the initial degenerate problem. The space H2(Ω) 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 dimensionn=2. 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Ω×[0,T] is given by
(107)(tk,xi)=(kΔt,iΔx)fork≥0,i∈ℕ.We note byUik the approximate value of the solution u(t,x) at the point (kΔt,iΔx)=(kΔt,i).The discrete solution at the iterationk is obtained at the discrete time tk and is denoted by Uk. An approximation of the temporal derivative is ut=∂u/∂t at this same time tk=kΔt is then (Uk+1-Uk)/Δt.


### 6.1. Numerical Schemes for the PDE Resolution

The discrete model calculating the mean envelope can be written as follows:(108)ut=-∑i,j=12‍∂xi1(gi,j∂xi1∂xj2u),
where ∂xiq denotes the derivative of order q∈{1,2} relatively to the spatial coordinate xi.To take into account the discontinuity of the diffusivity function, we use the harmonic mean ofgi,j defined by
(109)gi+1/2,j=[gi,j-1+gi+1,j-12]-1,gi-1/2,j=[gi-1,j-1+gi+1,j-12]-1,gi,j+1/2=[gi,j-1+gi,j+1-12]-1,gi,j-1/2=[gi,j-1+gi,j-1-12]-1.
Recall that g is strictly positive, which allows its inversion in the previous system.We define also the numerical derivatives ofv by
(110)Di+v(i,j)=v(i+1,j)-v(i,j)forwarddifferenceonthexi-direction,Di-v(i,j)=v(i,j)-v(i-1,j)backwarddifferenceonthexi-direction,Di0v(i,j)=v(i+12,j)-v(i-12,j)centraldifferenceonthexi-direction,Dj+v(i,j)=v(i,j+1)-v(i,j)forwarddifferenceonthexj-direction,Dj-v(i,j)=v(i,j-1)-v(i,j)backwarddifferenceonthexj-direction,
and finally we pose wiju(i,j)=Di+Dj-Dj+u(i,j). For positive integer, i and j, the numerical approximation Liju, of gi,j∂xi1∂xj2u in (108) is given by
(111)Lij=gi+1/2,jDi+wij-gi-1/2,jDi-wij+gi,j+1/2Dj+wji-gi,j-1/2Dj-wji.
We present the first natural explicit scheme.


#### 6.1.1.  An Explicit Scheme

An explicit scheme to solve the PDE is the following:(Uk+1-Uk)/Δt=-∑i,j=12‍LijUk or Uk+1=(I-Δt∑i,j=12‍Lij)Uk. Let us pose 𝔇=I-Δt∑i,j=12‍Lij.Thus,(112)Uk+1=𝔇UkwithU0=u0.
If N is the image or the signal size, the matrix 𝒟 is a sparse matrix of order N2.However, this explicit scheme requires for its stabilityΔt 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 Δt 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 theadditive operator splitting scheme (AOS) which is written as follows:
(113)Uk+1=12∑p=12‍(I+2ΔtLpp)-1(I-Δt∑i=12‍∑j≠i‍Lij)Uk.
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 k and is carried out once only for k=0.


#### 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:
(114)Uk+1=(∏p=12‍(I-ΔtLpp))-1(I+Δt∑i=12‍∑j≠i‍Lij)Uk.
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 u(x)=a3x3+a2x2+a1x+a0. The nullity of its derivatives of orders 3 and 1 on the edges leads to a3=0,  ∂u/∂x=2a2x+a1=0 on the boundary. Hence, a2=a1=0, which finally gives u(x)=a0. 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 3. 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Δt 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 1 both in space and in time.


#### 6.3.1. Some Examples

First, we present in Figure1 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.Figure2 represents the evolution of the PDE solution with the form given by (14). Upon iteration 40000, the solution almost evolves addition; it is potentially an envelope.In (a) input signal is in blue color. In red and green colors are the envelopes calculus upon iteration4000. In (b) the envelopes calculus upon iteration 40000. In (c) the envelopes calculus upon iteration 400000.
(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 400000In Figure3, 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.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.
(a)
Input signal is in blue color. In red and green colors are the envelopes at the convergence forθ=0.05
(b)
Input signal is in blue color. In red and green colors are the envelopes at the convergence forθ=0.5
(c)
Input signal is in blue color. In red and green colors are the envelopes at the convergence forθ=1In Figure4, we represent the envelope calculation with the form given by (17) for different tension factor θ=0.Envelope calculus by (17) with θ=0. 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.
(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 400000In Figure5, we make a comparison between our PDE interpolator and the cubic spline interpolation for an input signal s that equal to the superposition of two sinusoidal components with different amplitudes and frequencies, s1 and s1. 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 Figures6 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 s=s1+s2. 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 s1 of the signal s. 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.The envelope of the signal tests=s1+s2 sampling with 256 points.
(a)
Envelope after 40 iterations
(b)
Envelope after 400 iterationsThe envelope of the signal tests=s1+s2 sampling with 1024 points.
(a)
Envelope after 40000 iterations
(b)
Envelope after 400000 iterations


#### 6.3.3. Test on 2D Signal

Finally, in Figure8, 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.PDE interpolator for image restoration. An original image (a) is corrupted (b) and restored (c) by our 2D PDE-interpolator. Original image.
(a)
(b)
(c)


## 6.1. Numerical Schemes for the PDE Resolution

The discrete model calculating the mean envelope can be written as follows:(108)ut=-∑i,j=12‍∂xi1(gi,j∂xi1∂xj2u),
where ∂xiq denotes the derivative of order q∈{1,2} relatively to the spatial coordinate xi.To take into account the discontinuity of the diffusivity function, we use the harmonic mean ofgi,j defined by
(109)gi+1/2,j=[gi,j-1+gi+1,j-12]-1,gi-1/2,j=[gi-1,j-1+gi+1,j-12]-1,gi,j+1/2=[gi,j-1+gi,j+1-12]-1,gi,j-1/2=[gi,j-1+gi,j-1-12]-1.
Recall that g is strictly positive, which allows its inversion in the previous system.We define also the numerical derivatives ofv by
(110)Di+v(i,j)=v(i+1,j)-v(i,j)forwarddifferenceonthexi-direction,Di-v(i,j)=v(i,j)-v(i-1,j)backwarddifferenceonthexi-direction,Di0v(i,j)=v(i+12,j)-v(i-12,j)centraldifferenceonthexi-direction,Dj+v(i,j)=v(i,j+1)-v(i,j)forwarddifferenceonthexj-direction,Dj-v(i,j)=v(i,j-1)-v(i,j)backwarddifferenceonthexj-direction,
and finally we pose wiju(i,j)=Di+Dj-Dj+u(i,j). For positive integer, i and j, the numerical approximation Liju, of gi,j∂xi1∂xj2u in (108) is given by
(111)Lij=gi+1/2,jDi+wij-gi-1/2,jDi-wij+gi,j+1/2Dj+wji-gi,j-1/2Dj-wji.
We present the first natural explicit scheme.


### 6.1.1.  An Explicit Scheme

An explicit scheme to solve the PDE is the following:(Uk+1-Uk)/Δt=-∑i,j=12‍LijUk or Uk+1=(I-Δt∑i,j=12‍Lij)Uk. Let us pose 𝔇=I-Δt∑i,j=12‍Lij.Thus,(112)Uk+1=𝔇UkwithU0=u0.
If N is the image or the signal size, the matrix 𝒟 is a sparse matrix of order N2.However, this explicit scheme requires for its stabilityΔt 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 Δt 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 theadditive operator splitting scheme (AOS) which is written as follows:
(113)Uk+1=12∑p=12‍(I+2ΔtLpp)-1(I-Δt∑i=12‍∑j≠i‍Lij)Uk.
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 k and is carried out once only for k=0.


### 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:
(114)Uk+1=(∏p=12‍(I-ΔtLpp))-1(I+Δt∑i=12‍∑j≠i‍Lij)Uk.
Besides the approximation errors inherent in all numerical scheme, the fundamental difference between theses methods remains the computation time.


## 6.1.1.  An Explicit Scheme

An explicit scheme to solve the PDE is the following:(Uk+1-Uk)/Δt=-∑i,j=12‍LijUk or Uk+1=(I-Δt∑i,j=12‍Lij)Uk. Let us pose 𝔇=I-Δt∑i,j=12‍Lij.Thus,(112)Uk+1=𝔇UkwithU0=u0.
If N is the image or the signal size, the matrix 𝒟 is a sparse matrix of order N2.However, this explicit scheme requires for its stabilityΔt 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 Δt 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 theadditive operator splitting scheme (AOS) which is written as follows:
(113)Uk+1=12∑p=12‍(I+2ΔtLpp)-1(I-Δt∑i=12‍∑j≠i‍Lij)Uk.
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 k and is carried out once only for k=0.


## 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:
(114)Uk+1=(∏p=12‍(I-ΔtLpp))-1(I+Δt∑i=12‍∑j≠i‍Lij)Uk.
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 u(x)=a3x3+a2x2+a1x+a0. The nullity of its derivatives of orders 3 and 1 on the edges leads to a3=0,  ∂u/∂x=2a2x+a1=0 on the boundary. Hence, a2=a1=0, which finally gives u(x)=a0. 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 3. 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Δt 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 1 both in space and in time.


### 6.3.1. Some Examples

First, we present in Figure1 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.Figure2 represents the evolution of the PDE solution with the form given by (14). Upon iteration 40000, the solution almost evolves addition; it is potentially an envelope.In (a) input signal is in blue color. In red and green colors are the envelopes calculus upon iteration4000. In (b) the envelopes calculus upon iteration 40000. In (c) the envelopes calculus upon iteration 400000.
(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 400000In Figure3, 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.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.
(a)
Input signal is in blue color. In red and green colors are the envelopes at the convergence forθ=0.05
(b)
Input signal is in blue color. In red and green colors are the envelopes at the convergence forθ=0.5
(c)
Input signal is in blue color. In red and green colors are the envelopes at the convergence forθ=1In Figure4, we represent the envelope calculation with the form given by (17) for different tension factor θ=0.Envelope calculus by (17) with θ=0. 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.
(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 400000In Figure5, we make a comparison between our PDE interpolator and the cubic spline interpolation for an input signal s that equal to the superposition of two sinusoidal components with different amplitudes and frequencies, s1 and s1. 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 Figures6 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 s=s1+s2. 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 s1 of the signal s. 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.The envelope of the signal tests=s1+s2 sampling with 256 points.
(a)
Envelope after 40 iterations
(b)
Envelope after 400 iterationsThe envelope of the signal tests=s1+s2 sampling with 1024 points.
(a)
Envelope after 40000 iterations
(b)
Envelope after 400000 iterations


### 6.3.3. Test on 2D Signal

Finally, in Figure8, 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.PDE interpolator for image restoration. An original image (a) is corrupted (b) and restored (c) by our 2D PDE-interpolator. Original image.
(a)
(b)
(c)


## 6.3.1. Some Examples

First, we present in Figure1 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.Figure2 represents the evolution of the PDE solution with the form given by (14). Upon iteration 40000, the solution almost evolves addition; it is potentially an envelope.In (a) input signal is in blue color. In red and green colors are the envelopes calculus upon iteration4000. In (b) the envelopes calculus upon iteration 40000. In (c) the envelopes calculus upon iteration 400000.
(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 400000In Figure3, 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.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.
(a)
Input signal is in blue color. In red and green colors are the envelopes at the convergence forθ=0.05
(b)
Input signal is in blue color. In red and green colors are the envelopes at the convergence forθ=0.5
(c)
Input signal is in blue color. In red and green colors are the envelopes at the convergence forθ=1In Figure4, we represent the envelope calculation with the form given by (17) for different tension factor θ=0.Envelope calculus by (17) with θ=0. 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.
(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 400000In Figure5, we make a comparison between our PDE interpolator and the cubic spline interpolation for an input signal s that equal to the superposition of two sinusoidal components with different amplitudes and frequencies, s1 and s1. 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 Figures6 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 s=s1+s2. 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 s1 of the signal s. 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.The envelope of the signal tests=s1+s2 sampling with 256 points.
(a)
Envelope after 40 iterations
(b)
Envelope after 400 iterationsThe envelope of the signal tests=s1+s2 sampling with 1024 points.
(a)
Envelope after 40000 iterations
(b)
Envelope after 400000 iterations


## 6.3.3. Test on 2D Signal

Finally, in Figure8, 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.PDE interpolator for image restoration. An original image (a) is corrupted (b) and restored (c) by our 2D PDE-interpolator. Original image.
(a)
(b)
(c)


## 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 H1(Ω). 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.



---

*Source: 605035-2013-03-07.xml*