Nicolas Brisebarre
email: nicolas.brisebarre@ens-lyon.fr

Joldeş Mioara
email: mioara.joldes@ens-lyon.fr

Arénaire Lip

Chebyshev Interpolation Polynomial-based Tools for Rigorous Computing

Keywords: Rigorous Computing, Validated Numerics, Interpolation Polynomial, Chebyshev Polynomials, Taylor Models

Performing numerical computations, yet being able to provide rigorous mathematical statements about the obtained result, is required in many domains like global optimization, ODE solving or integration. Taylor models, which associate to a function a pair made of a Taylor approximation polynomial and a rigorous remainder bound, are a widely used rigorous computation tool. This approach benefits from the advantages of numerical methods, but also gives the ability to make reliable statements about the approximated function. Despite the fact that approximation polynomials based on interpolation at Chebyshev nodes offer a quasi-optimal approximation to a function, together with several other useful features, an analogous to Taylor models, based on such polynomials, has not been yet well-established in the field of validated numerics. This paper presents a preliminary work for obtaining such interpolation polynomials together with validated interval bounds for approximating univariate functions. We propose two methods that make practical the use of this: one is based on a representation in Newton basis and the other uses Chebyshev polynomial basis. We compare the quality of the obtained remainders and the performance of the approaches to the ones provided by Taylor models.

Introduction

Computers are used nowadays to quickly give numerical solutions to various global optimization, ODE solving or integration problems. However, traditional numeric methods usually provide only approximate values for the solution. Bounds for the approximation errors are only sometimes available, are not guaranteed to be accurate or are sometimes unreliable. In contrast, validated computing aims at providing rigorously verified information about solutions, in order to complete proofs, or to give rigorous mathematical statements about the obtained result.

Interval arithmetic [START_REF] Moore | Methods and Applications of Interval Analysis[END_REF] is a classical tool to perform validated computations with floatingpoint arithmetic. Intervals are well-suited to represent enclosures of real numbers on a machine. However, they propagate only information about function values, and fail to convey much information about the other properties about the function itself. In particular, when modeling functions with interval arithmetic, splitting the domain in subintervals is usually required. For some cases, known as "high dependency problems" [START_REF] Berz | Rigorous global search using Taylor models[END_REF][START_REF] Berz | Long-term stability of the tevatron by verified global optimization[END_REF][START_REF] Chevillard | Certified and fast computation of supremum norms of approximation errors[END_REF][START_REF] Chevillard | Efficient and accurate computation of upper bounds of approximation errors[END_REF] the number of necessary subintervals becomes unfeasible.

Taylor models [START_REF] Makino | Taylor models and other validated functional inclusion methods[END_REF][START_REF] Berz | Rigorous global search using Taylor models[END_REF][START_REF] Berz | Long-term stability of the tevatron by verified global optimization[END_REF], introduced by Berz and his group offer a remedy to this problem. They provide another way to rigorously manipulate and evaluate functions using floatingpoint arithmetic. They have been widely used for validated computing for global optimization and range bounding [START_REF] Makino | Rigorous Analysis of Nonlinear Motion in Particle Accelerators[END_REF][START_REF] Berz | Rigorous global search using Taylor models[END_REF][START_REF] Chevillard | Efficient and accurate computation of upper bounds of approximation errors[END_REF][START_REF] Berz | Long-term stability of the tevatron by verified global optimization[END_REF], solutions of ODEs [START_REF] Neher | On Taylor model based integration of ODEs[END_REF], quadrature [START_REF] Berz | New methods for high-dimensional verified quadrature[END_REF], etc.

A Taylor model (TM) of order n for a function f which is supposed to be n + 1 times continuously differentiable over an interval [a, b], is a rigorous polynomial approximation of f . More specifically, it is a couple (P, ∆) formed by a polynomial P of degree n, and an interval part ∆, such that f (x) -P (x) ∈ ∆, ∀x ∈ [a, b]. Roughly speaking, the polynomial can be seen as a Taylor expansion of the function at a given point. The interval ∆ provides the validation of the approximation, meaning that it provides an enclosure of all the approximation errors encountered (truncation, roundings).

A natural idea is to try to replace Taylor polynomials with better approximations such as minimax approximation, Chebyshev truncated series or interpolation polynomials (also called approximate Chebyshev truncated series when the points under consideration are Chebyshev nodes) for instance. The last two kind of approximations are of particular relevance for replacing Taylor polynomials since the series they define converge on domains better shaped for various usual applications than Taylor expansions (see, for instance, Section 2.7 of [START_REF] Boyd | Chebyshev and Fourier spectral methods[END_REF] for a more detailed account). Moreover, we can take advantage of numerous powerful techniques for computing these approximations. So far, the attempts for using these better approximations, in the context of rigorous computing, do not seem to have succeeded, see for example [START_REF] Makino | Taylor models and other validated functional inclusion methods[END_REF] for a comparison of existing techniques.

In this work we propose two approaches for computing models based on interpolation polynomials at Chebyshev nodes, what we call "Chebyshev interpolation models" (CM). The first method is based on Newton Basis and the second on Chebyshev polynomial basis. We believe that bringing a certified remainder to an approximate truncated Chebyshev series and providing effective tools for working with such models, opens the way to adapting to rigorous computing many numerical algorithms based on Chebyshev interpolation polynomials, for rigorous ODE solving, quadrature, etc.

The outline of the paper is the following. We first recall or prove various definitions and results required by our approaches in Section 2. In Section 3, we present in more details TMs and we discuss the use of better polynomial approximations. Then, we introduce the notion of "Chebyshev interpolation models" in Section 4. The CMs are implemented using multiple precision interval arithmetic in order to perform rigorous computing and yet, for the sake of clarity, we present their implementation in exact arithmetic in Section 5. We give some results and a comparison of our models with TMs in Section 6. We end with a brief conclusion and a mention of our future works on the subject.

Some preliminary statements about interpolation and Chebyshev polynomials

We first give a very short reminder on Chebyshev polynomials. A detailed presentation can be found in [START_REF] Borwein | Polynomials and Polynomial Inequalities[END_REF][START_REF] Rivlin | Chebyshev polynomials. From approximation theory to algebra[END_REF]. Then we state some interpolation results that we use in the sequel.

Some basic facts about Chebyshev polynomials

Over [-1, 1], Chebyshev polynomials can be defined as T n (x) = cos (n arccos x) , n ≥ 0. Since we consider functions over any interval I = [a, b], we define in the following the Chebyshev polynomials over I as

T [a,b] n (x) = T n 2x-b-a b-a . T [a,b] n+1 has n + 1 distinct real roots in [a, b],
called "Chebyshev nodes" since they are of utmost interest for interpolation:

x * i = a + b 2 + b -a 2 cos (i + 1/2) π n + 1 , i = 0, . . . , n. (1)
We now recall

Lemma 2.1 The polynomial W x * (x) = n i=0 (x -x * i )
, is the monic degree-n + 1 polynomial that minimizes the supremum norm over [a, b] of all monic polynomials in C[x] of degree at most n + 1. We have

W x * (x) = (b -a) n+1 2 2n+1 T [a,b] n+1 (x) and max x∈[a,b] |W x * (x)| = (b -a) n+1 2 2n+1 .

Brief overview of interpolation results used

Let I = [a, b] be an interval. Let f be a function that is at least n + 1 times continuously differentiable over I. Let {y i , i = 0, . . . , n} be a set of n + 1 points in I. There exists a unique polynomial P of degree ≤ n which interpolates f at these points [START_REF] Cheney | Introduction to Approximation Theory[END_REF]:

P (y i ) = f (y i ), i = 0, . . . , n, or if y i is repeated k times, P (j) (y i ) = f (j) (y i ), j = 0, . . . , k -1.
If all the points are distinct, this is called Lagrange interpolation. In the extreme case that all the y i are equal, P is just the Taylor polynomial of f at the considered point. Several algorithms and interpolation formulae in various basis exist for representing P , for example monomial basis, Lagrange, Newton, Barycentric Lagrange, Chebyshev basis [START_REF] Cheney | Introduction to Approximation Theory[END_REF][START_REF] Berrut | Barycentric Lagrange interpolation[END_REF][START_REF] Schatzman | Numerical Analysis, A Mathematical Introduction[END_REF]. The numerical properties (stability) of these formulas have been widely studied in the literature [START_REF] Higham | The numerical stability of barycentric Lagrange interpolation[END_REF].

Let us consider the polynomial P in Newton basis: P (x) = n i=0 c i N y,i (x), where N y,0 (x) = 1 and N y,i = i-1 j=0 (x-y j ), i = 1, . . . , n. The coefficients c i are the divided-differences f [y 0 , . . . , y i ] of f at the points y 0 , . . . , y i . As mentioned above, if k points coincide, it suffices to take the successive k -1 derivatives of f . Note that c i can be obtained thanks to the divideddifferences algorithm [START_REF] Gel | Calculus of finite differences[END_REF][START_REF] Schatzman | Numerical Analysis, A Mathematical Introduction[END_REF]. Moreover, the error between f and P is given [START_REF] Gel | Calculus of finite differences[END_REF][START_REF] Schatzman | Numerical Analysis, A Mathematical Introduction[END_REF] by:

∀x ∈ I, f (x) -P (x) = f [y 0 , . . . , y n , x]W y (x) (2)
with W y (x) = n i=0 (xy i ). By a repeated application of Rolle's theorem [START_REF] Cheney | Introduction to Approximation Theory[END_REF][START_REF] Gel | Calculus of finite differences[END_REF][START_REF] Schatzman | Numerical Analysis, A Mathematical Introduction[END_REF], we get:

∀x ∈ I, ∃ξ ∈ (a, b) s.t. f (x) -P (x) = 1 (n + 1)! f (n+1) (ξ)W y (x). (3)
In the sequel we denote the right member of this formula for the error by ∆ n (x, ξ). Lemma 2.1 suggests that a clever choice of interpolation points seems to be the Chebyshev nodes (1), which is indeed the case as mentioned in 3.2.2. We will take advantage of the following lemma which generalizes Lemma 5.12 of [START_REF] Zumkeller | Global Optimization in Type Theory[END_REF]:

Lemma 2.2 Under the assumptions on f and y i above, if f (n+1) is increasing (resp. decreasing) over I, then f [y 0 , . . . , y n , x] is increasing (resp. decreasing) over I.

Assume that f (n+1) is increasing. We know (see Chap. 4 of [START_REF] Schatzman | Numerical Analysis, A Mathematical Introduction[END_REF] for instance) that, for all x ∈ [a, b], we have

f [y 0 , . . . , y n , x] = 1 0 t 1 0 • • • tn 0 f (n+1) (y 0 + t 1 (y 1 -y 0 ) + • • • + t n+1 (x -y n ))dt 1 • • • dt n+1 . Let x, y ∈ [a, b], x ≤ y, let Z n denote y 0 + t 1 (y 1 -y 0 ) + • • • + t n (y n -y n-1 ) -t n+1 y n , we notice that f [y 0 , . . . , y n , y] -f [y 0 , . . . , y n , x] = 1 0 t 1 0 • • • tn 0 f (n+1) (Z n + t n+1 y) -f (n+1) (Z n + t n+1 x) dt 1 • • • dt n+1 . Since t n+1 ≥ 0, we have Z n + t n+1 y ≥ Z n + t n+1 x. As f (n+1
) is increasing, it follows that the integrand is nonnegative, which implies f [y 0 , . . . , y n , y] ≥ f [y 0 , . . . , y n , x].

3 Taylor models vs. using better approximations

Basic principles of Taylor models

Computing a TM consists in computing a polynomial together with an interval bound by applying simple rules recursively on the structure of the function f . In fact, for functions like trigonometric, exponential, logarithmic functions, as well as operations like 1/x or the power function, all referred to in this article as basic functions (or as intrinsics in [START_REF] Makino | Rigorous Analysis of Nonlinear Motion in Particle Accelerators[END_REF]), bounds for the remainders can be easily computed. For composite functions, TMs offer usually a much tighter bound than the one directly computed for the whole function, for example using automatic differentiation [1, [START_REF] Rall | The arithmetic of differentiation[END_REF][START_REF] Moore | Methods and Applications of Interval Analysis[END_REF]. A meaningful comparison for this phenomenon is provided in [START_REF] Chevillard | Efficient and accurate computation of upper bounds of approximation errors[END_REF][START_REF] Makino | Taylor models and other validated functional inclusion methods[END_REF]. Here we provide the reader with a quick overview of the situation. Let u be a basic function defined on I and v a basic function defined on J , an interval enclosing the image u(I). Consider the Taylor expansion of a composite function v • u about x 0 ∈ I.

The Taylor remainder is given by the Lagrange formula: for all x ∈ I, there exists ξ ∈ (a, b) such that

R(x) = (v • u) (n+1) (ξ) (n + 1)! (x -x 0 ) n+1 .
When bounding the remainder by means of automatic differentiation, an interval K enclosing the values of the n + 1-th derivative of this composite function (v • u) (n+1) (I) is obtained by performing many recursive operations involving enclosures of u (i) (I) and v (i) (J ) which finally may produce a considerable overestimation in the remainder [START_REF] Chevillard | Efficient and accurate computation of upper bounds of approximation errors[END_REF].

In contrast, TMs compute polynomials and interval bounds directly for v and u, and then use algebraic rules for computing a TM for the composition v •u. Since both v and u are basic functions, evaluating with interval arithmetic their n + 1-th derivative can be done in a fast way using simple formulae and does not lead to serious overestimation.

Morevoer, most of the functions we deal with have Taylor series whose coefficients decrease. In particular, when the functions are analytic over a sufficiently large domain (a disk of radius > 1 suffices), the magnitude of the coefficients of the underlying polynomial decreases exponentially (this is a consequence of Cauchy's integral formula). Hence when performing the composition of two such models, the intervals contributing to the final remainder become smaller for monomials of higher degree. This leads to a reduced overestimation in the computed remainder.

In conclusion, we highlight that, in practice, it is significantly more suitable to use a twostep procedure for handling composite functions: first step consists in computing models (P, ∆) for all basic functions; second, apply algebraic rules specifically designed for handling operations with these mixed models instead of operations with the corresponding functions.

Substituting tighter polynomial approximations to Taylor polynomials

It is well known that Taylor polynomials can be fairly poor approximations to functions, except perhaps on very small intervals. The idea of using better approximation polynomials, for example Chebyshev truncated series, in the context of validated computations was introduced in [START_REF] Epstein | Ultra-arithmetic i: Function data types[END_REF] under the name ultra-arithmetic, which we first briefly present. Then, we turn to a first introduction of the CMs and give the basic principles which they rely on. We finally address the case of minimax approximation. The case of truncated Chebyshev series is mentioned in the conclusion of this paper.

A previous attempt

As explained in [START_REF] Makino | Taylor models and other validated functional inclusion methods[END_REF], in the setting of ultra-arithmetic, the advantages of non-Taylor approximations cannot be explicitly maintained due to several drawbacks. It is noted in [START_REF] Makino | Taylor models and other validated functional inclusion methods[END_REF] that the Taylor representation is a special case, because for two functions f 1 and f 2 , the Taylor representation of order n, for the product f 1 • f 2 can be obtained merely from the Taylor expansions of f 1 and f 2 , simply by multiplying the polynomials and discarding the orders n + 1 to 2n. On the contrary, the Chebyshev truncated series of a product f 1 • f 2 can in general not be obtained from the Chebyshev representations of the factors f 1 and f 2 , and no operation analogous to TMs multiplication is given. Moreover, there is no systematic treatment of common basic functions. Finally, [START_REF] Makino | Taylor models and other validated functional inclusion methods[END_REF] explains that the methods developed in ultra-arithmetic can lead to an increase in the magnitude of the coefficients, which will increase both the computational errors and the difficulty of finding good interval enclosures of the polynomials involved.

Chebyshev interpolation polynomial

A natural idea is to consider an interpolation polynomial instead of a Taylor approximation: if the polynomial interpolates the function at Chebyshev nodes, such a polynomial, which we call Chebyshev interpolation polynomial, can usually provide a near-optimal approximation of f . See e.g. [START_REF] Powell | On the maximum errors of polynomial approximations defined by interpolation and by least squares criteria[END_REF] which states an effective measure of this property.

Trefethen [START_REF] Trefethen | Computing numerically with functions instead of numbers[END_REF] uses the idea of representing the functions by Chebyshev interpolation polynomials. He chooses their expansion length in such a way as to maintain the accuracy of the approximation obtained close to machine precision. Moreover, the idea of using basic functions and then consider an algebra on them is used. The developed software, Chebfun was successfully used for numerically solving differential equations, quadrature problems. However, the Chebfun system does not provide any validation of the results obtained. As mentioned in [START_REF] Trefethen | Computing numerically with functions instead of numbers[END_REF], the aim of this system are numerical computations and no formal proof or safeguards are yet implemented to guarantee a validated computation, although it seems to be a future project.

Towards Chebyshev interpolation models

The interpolation polynomial itself is easy to compute and a wide variety of methods and various basis representations exist (we refer the reader to any numerical analysis book, see [START_REF] Schatzman | Numerical Analysis, A Mathematical Introduction[END_REF] for instance, for its computation). Another advantage of interpolation polynomials is that a closed formula, hence an explicit bound, for the remainder exists, cf. ( 2) and (3). Lemma 2.1 and Formula (3) induces that for a Chebyshev interpolation polynomial of degree n, roughly speaking, compared to a Taylor remainder, the interpolation error will be scaled down by a factor of 2 -n . For bounding the remainder, the only remaining difficulty is to bound the term f (n+1) (ξ) for ξ ∈ I. However, as briefly discussed above and in [START_REF] Chevillard | Efficient and accurate computation of upper bounds of approximation errors[END_REF] the advantage of using directly this formula for the remainder can not be effectively maintained if the interval bound is obtained using automatic differentiation because of the overestimation. Hence in what follows we will try to adapt the "basic bricks" approach used by TMs to one using interpolation polynomials.

Minimax approximation

One might first try to directly use the minimax polynomial which is the polynomial which minimizes the supremum norm of the approximation error. In fact, such a polynomial can be obtained only through an iterative procedure, namely Remez algorithm, which can be considered too computationally expensive in our context. But the main drawback is linked to the computation of the remainder: either obtaining a certified remainder in such a procedure, raises a significant dependency problem as discussed in [START_REF] Chevillard | Certified and fast computation of supremum norms of approximation errors[END_REF][START_REF] Chevillard | Efficient and accurate computation of upper bounds of approximation errors[END_REF], or the closed formula, due to Bernstein [START_REF] Elliott | Error of truncated Chebyshev series and other near minimax polynomial approximations[END_REF], does not bring any improvement over Formula (3).

Chebyshev Interpolation Models

In the following we consider a Chebyshev Interpolation model (CM) of degree n for a function f over an interval I, as a couple (P, ∆), in the following sense: P is a polynomial of degree n which is "closely related" to the Chebyshev interpolation polynomial for f , ∆ is an interval enclosure for the remainder f -P . Specifically, from the mathematical point of view, for basic functions f , this polynomial P coincides with the Chebyshev interpolation polynomial and the computation of ∆ is detailed in 5.1.2. For composite functions, we will use some algebraic rules defined in Section 5 in a similar manner to the TMs arithmetic. We insist that our setting is multiple-precision interval arithmetic and all algorithms chosen are adapted for and implemented specifically for it. Since a thorough description of their implementation is tedious we indicate some implementation hints in 4.2. Before doing this, we explain why we chose to use Newton and Chebyshev bases for performing operations over CMs.

Choice of representation basis

When computing an interpolation polynomial P two key choices are: the interpolating nodes (which are Chebyshev nodes (1), in this paper), and the basis used for implementation.

In our case, however, a third essential requirement emerges: one has to be able to find suitable algebraic rules for operations with models (P, ∆) that will not lead to overestimation of the remainder in the resulted model. This implies that addition, multiplication or composition of models made out of polynomials P , together with an interval remainder bound, has to be an effective process not only in terms of performance and quality of the polynomial obtained, but also in terms of quality of remainder. We note, that in what conaddition, the representation of P in any basis will suffice, since this is a linear operator. For multiplication, as discussed in 3.2.1, it has not been obvious so far how to devise an efficient algorithm that, given two functions f 1 and f 2 and their models (P 1 , ∆ 1 ) and (P 2 , ∆ 2 ), where the polynomial part is obtained through an interpolation process, is able to efficiently compute a model (P, ∆) for f 1 • f 2 . Furthermore, composition rules f 1 • f 2 were even more difficult to obtain. On the contrary, these operations are straightforward for TMs, which is one of their incontestable advantages.

Consequently, we were led to searching for suitable basis representations for P such that all these requirements are successful. We eliminated a Lagrange basis representation for P not only because of its poor numerical properties [START_REF] Higham | The numerical stability of barycentric Lagrange interpolation[END_REF], but also because of the n -1 terms containing each n -1 products of xx * i . This leads to high overestimations on the interval remainder when a composition operation was tried.

In what concerns using barycentric Lagrange basis versus Newton basis, it is proved in [START_REF] Higham | The numerical stability of barycentric Lagrange interpolation[END_REF] that the first has better numerical properties. However, Newton basis also has good numerical properties for certain orderings of points [START_REF] Higham | The numerical stability of barycentric Lagrange interpolation[END_REF]. In our case, barycentric Lagrange has the disadvantage that the polynomial P is represented as

n i=0 w i x -x * i f (x * i ) / n i=0 w i x -x * i .
It seems cumbersome to implement multiplication and composition of two models using this representation for P . We note that when composition should be implemented, the reciprocal 1/(xx * i ) should be designed. This leads to considerable overestimations of the remainder. We note that Newton basis can be seen as an "immediate extension" to Taylor basis, so it is a natural try to make and leads to an adaptation of TMs algorithms. In what follows we consider the Newton basis at Chebyshev nodes taken in decreasing order, as given by Formula (1). This order has good numerical properties, see for example Section 5.3 [START_REF] Higham | Accuracy and Stability of Numerical Algorithms[END_REF]. In Section 5 we show how multiplications and composition operations can be implemented.

As for Chebyshev basis, we were led to this choice mainly because of the remarkable properties and success of numerical methods based on Chebyshev series expansions [START_REF] Boyd | Chebyshev and Fourier spectral methods[END_REF][START_REF] Berrut | Barycentric Lagrange interpolation[END_REF].

Moreover, as we will see in Section 5.4, for most of the functions we consider, the rapid decreasing of the coefficients in the representations in the last two bases allow for small overestimation when handling CMs.

We now give some basic operations that will be used during the operations on the models.

Operations with polynomials in Newton basis Consider two polynomials of degree n in

Newton basis:

P (x) = n i=0 p i N x * ,i (x) and Q(x) = n i=0 q i N x * ,i (x).
Addition. Adding two polynomials in Newton basis is straightforward:

P (x) + Q(x) = n i=0 (p i + q i )N x * ,i (x)
Multiplication. When multiplying two polynomials in Newton basis, we need that "the lower part" of the result be represented in this basis also, and the "upper part", should be represented such that it can be easily bounded with interval arithmetic. Hence, we chose to represent the product P Q in the basis N x * ,0 , . . . , N x * ,n , W x * N x * ,0 , . . . , W x * N x * ,n-1 . We have

P (x) • Q(x) = G(x) + W x * (x)H(x), where G(x) = n i=0 g i N x * ,i (x) and H(x) = n-1 i=0 h i N x * ,i (x).
One of the advantages of this representation is that we gave an exact interval bound for W x * in Section 2.1. Moreover, it can be easily shown that G(x) is the Chebyshev interpolation polynomial of P (x) • Q(x), i.e. of f g.

In order to compute the coefficients of G and H one can use several techniques. We mention that based on [START_REF] Bostan | Polynomial evaluation and interpolation on special sets of points[END_REF] a conversion back to monomial basis, followed by an interpolation seems to be the fastest asymptotically (needing O(M (n) log n) operations, where M (n) denotes the cost of multiplying univariate polynomials of degree less than n). However, in our current multiple precision interval arithmetic implementation, we use an O(n 2 ) algorithm based on applying the divided differences algorithm first for computing the coefficients of G and then for computing the coefficients of H = (P Q -G)/W x * .

Interval Range Bounding. We use a Horner-like algorithm [START_REF] Higham | Accuracy and Stability of Numerical Algorithms[END_REF] for bounding the range of a polynomial in Newton basis, which takes O(n) operations.

Operations with polynomials in Chebyshev basis Consider two polynomials of degree n

in Chebyshev basis:

P (x) = n i=0 p i T [a,b] i (x) and Q(x) = n i=0 q i T [a,b] i (x).
Addition. Adding two polynomials in Chebyshev basis is straightforward:

P (x)+Q(x) = n i=0 (p i + q i )T [a,b] i (x) and takes O(n) operations.
Multiplication. Their product can be expressed in Chebyshev basis as follows:

P (x) • Q(x) = 2n k=0 c k T [a,b] k (x), where c k = ( |i-j|=k p i • q j + i+j=k p i • q j )/2. This identity can be obtained noting [20] that T [a,b] i (x) • T [a,b] j (x) = (T [a,b] i+j + T [a,b]
|i-j| )/2. The cost using this simple identity is O(n 2 ) operations.

Interval Range Bounding. For our purpose, we will use the following identity for interval bounding range of a polynomial in Chebyshev basis, which takes

O(n) operations: ∀x ∈ [a, b], P (x) ∈ p 0 + n i=1 p i • [-1, 1].

Rigorous implementation of the interpolation process

In the following, we denote an interval x as a pair x = [x, x]. For a polynomial P with real number coefficients, we denote by P a polynomial obtained by replacing its coefficients with intervals which enclose them.

We implemented all the operations involved using multiple-precision interval arithmetic and all the algorithms we use are straightforwardly adaptable in this context. Hence, from the implementation point of view, the only notable change is that instead of polynomials with real number coefficients we have polynomials with tight interval coefficients. This means that in our implementation, for a function f we obtain a model (P , ∆), such that ∀x ∈ [a, b], f (x) ∈ P (x) + ∆. This design choice allows us to take into account all the rounding errors, because for each operation with interval arithmetic an outward rounding is performed. Moreover, it is proved in [START_REF] Revol | Newton's algorithm using multiple precision interval arithmetic[END_REF] that when evaluating a function ϕ over a point interval x = [x, x], the interval enclosure of ϕ(x) can be made arbitrarily tight by increasing the precision used for evaluation. In our case the computations needed for the coefficients of the polynomial are done with initially almost point intervals, so the overestimation of the coefficients can be made as small as possible by increasing the precision used. However, if needed, a certified floating-point approximation polynomial can be easily obtained from a CM (P , ∆). Let us consider P (x) = n i=0 c i β i (x) where β i is the ith element of the considered basis (Newton or Chebyshev). It suffices to take the middle points t i of c i as the coefficients of the approximation polynomial P (x) = n i=0 t i β i (x) and then to compute

a simple interval bound δ = n i=0 [c i -t i , c i -t i ] • β i (x)
, using the methods presented in Section 4.1. Then the error between the function f and its approximation polynomial P is bounded by δ + ∆.

Chebyshev Interpolation Models with exact arithmetic

We follow the two-step approach specific for TMs, as mentioned in Section 3. As detailed below, firstly, we compute models for basic functions; secondly we apply specifically designed algebraic operations on such models.

Basic functions

We consider as basic functions the trigonometric, exponential, logarithmic functions, 1/x, the power function, or any other function for which specific properties can be exploited. For such functions we compute directly a model (P, ∆) formed by an interpolation polynomial P and an interval bound for the remainder ∆.

Computation of the interpolation polynomial

We mentioned in 2.2 how to express an interpolation polynomial in Newton basis, using the divided differences procedure. The number of operations necessary to compute the coefficients of the interpolation polynomial is 3n 2 /2 [START_REF] Higham | Accuracy and Stability of Numerical Algorithms[END_REF].

We now need to know how to represent it using Chebyshev basis. The Chebyshev interpolation polynomial P can be expressed using the collocation method [START_REF] Mason | Chebyshev polynomials[END_REF] as follows:

P (x) = n i=0 p i T [a,b] i (x), with p i = n k=0 2 n+1 f (x * k )T [a,b] i (x * k ).
Using directly this formula, the computation cost is O(n 2 ) operations. It is known [START_REF] Platte | Chebfun: A new kind of numerical computing[END_REF] that the usage of Fast Fourier Transform, can speed-up this computation to O(n log n) operations. However, note that in this case, an interval arithmetic adaptation of FFT should be considered.

Computation of the remainder

We can compute an enclosure of the remainder f -P over [a, b] using Formula (3):

∆ n (x, ξ) = f (n+1) (ξ)W x * (x)/(n + 1)!, where x, ξ ∈ [a, b]. This reduces to bounding f (n+1) over [a, b],
which does not pose any problem for basic functions, since simple formulae are available for their derivatives. Moreover, W x * (x) can be bounded straightforwardly using Lemma 2.1.

We remark that, thanks to Lemma 2.2, an exact bound for the remainder can be computed when one can show that f (n+1) is increasing (resp. decreasing) over [a, b]. In such cases, one can use [START_REF] Benoit | Chebyshev expansions for solutions of linear differential equations[END_REF]

and bound f [x * 0 , . . . , x * n , x] by [ f [x * 0 , . . . , x * n , a], f [x * 0 , . . . , x * n , b] ] (resp. [ f [x * 0 , . . . , x * n , b], f [x * 0 , . . . , x * n , a] ]).
For basic functions, checking whether f (n+2) has constant sign over [a, b] is simple. This remark makes it possible to obtain smaller remainders and strengthens the effectiveness of the "basic bricks" approach.

Bounding the interpolation polynomial

When using this approach, one needs to compute bounds for the range of polynomials involved in such models. We denote by B(P ) an interval range bound for a polynomial P , over a given interval I. Several methods exist and a trade-off between their speed and the tightness of the bound is usually considered. For TMs, the fastest but "rough" method is a Horner-like interval evaluation. More complicated schemes exist, that usually give tighter bounds: LDB, QDB [START_REF] Berz | Rigorous global search using Taylor models[END_REF], or using a conversion to Bernstein basis [START_REF] Nataraj | Global optimization with higher order inclusion function forms part 1: A combined taylor-bernstein form[END_REF][START_REF] Zumkeller | Global Optimization in Type Theory[END_REF]. For univariate polynomials, slower techniques based on root isolation can be used for a very tight polynomial bounding [START_REF] Rouillier | Efficient isolation of polynomial's real roots[END_REF]. In our case, we focused on speed, and so, when considering interval range bounding for polynomials in Newton or Chebyshev basis we used two simple methods described in Section 4.1. Similarly to TMs, for a penalty in speed, more refined algorithms can also be plugged-in. However the simple techniques we used, proved effective in most of the examples we treated so far.

Addition and multiplication

In what follows we consider two CMs for two functions f 1 and f 2 , over I, of degree n: (P 1 , ∆ 1 ) and (P 2 , ∆ 2 ). Adding the two models is done by adding the two polynomials and the remainder bounds:

(P 1 , ∆ 1 ) + (P 2 , ∆ 2 ) = (P 1 + P 2 , ∆ 1 + ∆ 2 ).
Note that adding two polynomials in Newton or Chebyshev basis is straightforward and has a linear complexity, see Section 4.1.

For proving the correction we note that:

∀x ∈ I, ∃δ 1 ∈ ∆ 1 and δ 2 ∈ ∆ 2 s.t. f 1 (x) - P 1 (x) = δ 1 and f 2 (x)-P 2 (x) = δ 2 . Hence f 1 (x)+f 2 (x)-(P 1 (x) + P 2 (x)) = δ 1 +δ 2 ∈ ∆ 1 +∆ 2 .
For multiplication, we have

f 1 (x) • f 2 (x) = P 1 (x) • P 2 (x) + P 2 (x) • δ 1 + P 1 (x) • δ 2 + δ 1 • δ 2 .
We observe that P 1 • P 2 is a polynomial of degree 2n. Depending on the basis used, we split it into two parts: the polynomial consisting of the terms that "do not exceed n", (P 1 • P 2 ) 0...n and respectively the upper part (P 1 • P 2 ) n+1...2n , for the terms of the product P 1 • P 2 whose "order exceeds n". Now, a CM for f 1 • f 2 can be obtained by finding an interval bound ∆ for all the terms except P = (P 1 • P 2 ) 0...n :

∆ = B((P 1 • P 2 ) n+1...2n ) + B(P 2 ) • ∆ 1 + B(P 1 ) • ∆ 2 + ∆ 1 • ∆ 2 .
The interval bound for the polynomials involved can be computed as discussed in 5.1.3.

In our current setting, cf. Section 4.1, the number of operations necessary to multiply two such models is O(n 2 ).

For the sake of completeness, we mention that multiplying a CM with a constant, reduces to multiplying the polynomial and the remainder with the respective constant.

Composition

When the model for f 1 • f 2 is needed, we can consider (f 1 • f 2 )(x) as function f 1 evaluated at point y = f 2 (x). Hence, we have to take into account the additional constraint that the image of f 2 to be included in the definition range of f 1 . This can be checked by a simple interval bound computation of B(P 2 ) + ∆ 2 . Then we have:

(f 1 • f 2 )(x) -P 1 (f 2 (x)) ∈ ∆ 1 (4)
which implies that

(f 1 • f 2 )(x) -P 1 (P 2 (x) + ∆ 2 ) ∈ ∆ 1 (5)
In this formula, the only polynomial coefficients and remainders involved are those of the CMs of f 1 and f 2 which are basic functions. As we have seen above, fairly simple formulae exist for computing the coefficients and remainders of such functions. However, when using Formula (4), it is not obvious how to extract a polynomial and a final remainder bound from it. In fact, we have to reduce this extraction process to performing just multiplications and additions of CMs. A similar idea is used for composing TMs [START_REF] Makino | Taylor models and other validated functional inclusion methods[END_REF][START_REF] Zumkeller | Global Optimization in Type Theory[END_REF].

In our case, the difference is that P 1 and P 2 are polynomials represented in Newton/Chebyshev basis, and not in the monomial basis. In consequence, for computing the composition, we had to use a different algorithm. It is worth mentioning that a simple change of basis back and forth to monomial basis will not be successful. The problem is that the multiplications and additions used in such a composition process do not have to add too much overestimation to the final remainder. As we discussed in Section 3, for expansions of most of the functions we address, the size of the coefficients for the representation in the monomial basis is bounded by a decreasing sequence. Hence the contributions of the remainders in such a recursive algorithm are smaller and smaller. On the contrary, for interpolation polynomials, the coefficients represented in monomial basis oscillate too much and have poor numerical properties. Hence, a direct application of the principle of composing TMs will not be successful.

When using Newton basis, we perform the composition using a Horner-like algorithm [START_REF] Higham | Accuracy and Stability of Numerical Algorithms[END_REF]. It takes a linear number of operations between models.

Algorithm 1 Horner-like composition of CMs in Newton basis: Composing (P 1 , ∆ 1 ) with (P 2 , ∆ 2 ) /*We denote by P 1j the jth coefficient of P 1 */ /*and by x * j the jth interpolation point for P 1 */ (C n , R n ) := (P 1n , [0, 0])

For j = n -1, . . . , 0 do (C j , R j ) := ((P 2 , ∆ 2 ) -(x * j , [0, 0])) • (C j+1 , R j+1 ) + (P 1j , [0, 0]) ; Return (C 0 , R 0 + ∆ 1 )
When using Chebyshev basis, we perform the composition using an adaptation of Clenshaw algorithm [START_REF] Mason | Chebyshev polynomials[END_REF]. Algorithm 2 is used for efficient evaluation of a Chebyshev sum c i T i (x). It reduces evaluation of such polynomials to basic additions and multiplications and it is as efficient as Horner form for evaluating a polynomial as a sum of powers using nested multiplications. In our case, the variable x where the sum is to be evaluated is a CM, the multiplications and additions are operations between CMs. Moreover, using this algorithm, we perform a linear number of such operations between models.

Clenshaw-like composition of CMs in Chebyshev basis: Composing

(P 1 , ∆ 1 ) with (P 2 , ∆ 2 ) (C n+2 , R n+2 ) := (0, [0, 0]) /*CMs for 0*/ (C n+1 , R n+1 ) := (0, [0, 0]) /*We denote by P 1j the jth coefficient of P 1 */ For j = n, . . . , 1 do (C j , R j ) := 2 • (P 2 , ∆ 2 ) • (C j+1 , R j+1 ) -(C j+2 , R j+2 ) + (P 1j , [0, 0]) ; (C, R) := (P 2 , ∆ 2 ) • (C 1 , R 1 ) -(C 2 , R 2 ) + (P 10 , [0, 0]) Return (C, R + ∆ 1 )
The following examples (6 -8) are similar to some presented in [START_REF] Berrut | Barycentric Lagrange interpolation[END_REF]. There, the authors computed the minimax polynomials for these functions. Evidently, the polynomials obtained with CMs have a higher approximation error than the minimax, however, it is important to notice that in these tricky cases the remainder bound obtained for the CMs stays fairly good and is much better than the one obtained from TMs.

Examples 8 -9 present the case when the definition domain of the function is close to a singularity. As seen in these examples, when a direct interpolation process is used for a composite function, unfortunately, one can not apply Lemma 2.2 for bounding the remainder. Consequently, the bound obtained for the remainder is highly overestimated. However, when using the approach based on "basic bricks" and CMs benefit from it, yielding a much better remainder bound.

Example 10 deals with a function which is associated to the classical Runge phenomenon. Firstly, since the complex singularities of the function f defined by f (x) = 1/(1 + 4x 2 ) are close to the definition interval I, the Taylor polynomial is not a good approximation. Then, the interpolation method gives unfeasible interval remainder bounds due to the overestimation of the n + 1th derivative of the function f . TMs and CM1 fail both from the same cause: when computing a model for f = g • h one needs to compose the models for g(y) = 1/y and h(x) = 1 + 4x 2 . The model for h is simple to compute. However, one has to take into account that the model of g has to be computed over an interval range enclosure of the model of h. When such an enclosure is computed, using the presented simple methods for polynomial range bounding we have an overestimation of the image of h. This leads in fact to an interval that contains 0. Hence, the model for g which is not defined in 0 can not be computed. On the contrary, the overestimation in polynomial range bounding with the method of CM2 is smaller and we have a feasible remainder in this case.

We do not give specific timings, since for the moment our implementation is rather a "proof of concept" one. The algorithms necessary for TMs seem to be slightly simpler and hence more efficient than the ones needed by our approaches: we observed a factor between 2 and 3 of speed-up in favor of TMs, for cases when Taylor approximations are fairly good ones (tight intervals, analytic functions). However, we note that in some cases (when considering larger intervals or composed functions), in order to attain the same quality for the remainder, TMs need a much higher polynomial degree, more computation time, or they have to be applied over many subintervals, which favors the usage of CMs.

Conclusion and future work

We introduced two approaches for computing "Chebyshev interpolation models", a tool which is potentially useful in various rigorous computing applications. Currently, they always achieve smaller remainders than Taylor models but require nevertheless more computing time in some cases. This work is preliminary, in two senses. First, we did not address in this paper the opportunity to use Chebyshev truncated series instead of Chebyshev interpolation polynomials. Actually, this approach is a work in progress and seems very promising since the quality of the remainder should remain at least as good as the one provided by CM and the complexity of basic bricks computations should be lowered. This issue of complexity has to be addressed if we want CMs to replace TMs in most of univariate applications. The techniques developed in [START_REF] Rebillard | Étude théorique et algorithmique des séries de Chebyshev solutions d'équations différentielles holonomes[END_REF][START_REF] Benoit | Chebyshev expansions for solutions of linear differential equations[END_REF] could prove useful to achieve this goal. Secondly, we address the computation of such models for univariate functions only. A long-term project of ours is to extend these methods to the multivariate case.

The methods presented in this paper should be available in the months to come in the Sollya1 tool.

Table 1 :

 1 10 -14 1.19 • 10 -14 1.13 • 10 -14 1.19 • 10 -14 1.13 • 10 -14 1.22 • 10 -11 1.16 • 10 - 2 arctan(x), [-0.25, 0.25], 15 7.89 • 10 -15 7.89 • 10 -15 7.95 • 10 -17 7.89 • 10 -15 7.95 • 10 -17 2.58 • 10 -10 3.24 • 10 - Examples of bounds obtained by several methods

	No	f (x), I , n	CM1 bound	CMs CM2 bound	Exact bound	Interpolation Interp. bound Exact bound	TM bound	TMs	Exact bou
	1 1.19 • 3 sin(x), [3, 4], 10 arctan(x), [-0.9, 0.9], 15 5.10 • 10 -3 4 exp(1/ cos(x)), [0, 1], 14 6.69 • 10 -7 5 exp(x) log(2+x) cos(x) , [0, 1], 15 1.70 • 10 -8	5.10 • 10 -3 5.22 • 10 -7 9.11 • 10 -9	1.76 • 10 -8 4.95 • 10 -7 2.21 • 10 -9	5.10 • 10 -3 0.11 0.18	1.76 • 10 -8 6.10 • 10 -7 2.68 • 10 -9	1.67 • 10 2 9.06 • 10 -3 1.18 • 10 -3		5.70 • 10 2.59 • 10 3.38 • 10
	6 7 8 9 10	sin(exp(x)),[-1, 1], 10 tanh(x + 0.5) -tanh(x -0.5), [-1, 1], 10 √ x + 1.0001, [-1, 0], 10 √ x + 1.0001 • sin(x), [-1, 0], 10 1 1+4x 2 , [-1, 1], 10	4.10 • 10 -6 8.48 • 10 -3 3.64 • 10 -2 3.10 • 10 -2 +∞	9.47 • 10 -5 1.75 • 10 -3 3.64 • 10 -2 3.32 • 10 -2 1.13 • 10 -2	3.72 • 10 -6 4.88 • 10 -7 3.64 • 10 -2 3.08 • 10 -2 6.17 • 10 -3	4.42 • 10 -3 8.48 • 10 -3 3.64 • 10 -2 3.21 • 10 33 1.50 • 10 7	3.72 • 10 -6 4.88 • 10 -7 3.64 • 10 -2 3.08 • 10 -2 4.95 • 10 -3	2.96 • 10 -2 8.68 0.11 0.12 +∞		1.55 • 10 2.96 • 10 0.11 9.83 • 10 8.20 • 10

http://www.math.uni-wuppertal.de/~xsc/software/intpakX/

http://sollya.gforge.inria.fr/

Acknowledgements

The authors warmly thank the reviewers for their careful work. It made it possible to icantly improve the manuscript. They also would like to thank Martin Berz, Sylvain Chevillard, Kyoko Makino, Jean-Michel Muller and Markus Neher for many useful discussions.

Growth of the coefficients and overestimation

The overestimation does not grow too much during the recursion process. This is due to the nice convergence properties of the series expansions in Newton Basis or in Chebyshev polynomial basis. One can prove for instance that if the function under consideration is analytic over a "sufficiently large" domain of C which contains [a, b], the coefficients of the expansion in Newton [START_REF] Gel | Calculus of finite differences[END_REF] or in Chebyshev [START_REF] Boyd | Chebyshev and Fourier spectral methods[END_REF] bases decrease exponentially. As with TMs, when composing two such models, the intervals contributing to the final remainder become smaller for higher coefficients, which yields a reduced overestimation in the final remainder.

Experimental results

implemented the two methods in Maple, using the IntpakX 1 package. The Maple code can be downloaded from http://www.ens-lyon.fr/LIP/Arenaire/Ware/ChebModels. Table 1 shows the quality of some absolute error bounds obtained with the following presented methods: CMs, direct interpolation combined with AD, and TMs respectively. Each row of the table represents one example. The function f , the interval I and the degree n of the approximation polynomial are given in the first column.

We computed, as explained in Section 5, two CMs: one using polynomials in Newton basis CM 1 = (P 1 , ∆ 1 ) and the other CM 2 = (P 2 , ∆ 2 ) using Chebyshev basis. Let In order to observe the amount of overestimation we computed the exact error bounds ε 1 = sup x∈I {|f (x) -P 1 (x)|} and ε 2 = sup x∈I {|f (x) -P 2 (x)|} and we provide in the fourth column the minimum of the two:

We give in fifth and sixth columns the computed remainder bounds and the exact error obtained when an interpolating polynomial is directly used (directly means that the remainder is computed using (3) and automatic differentiation). Finally we present the computed remainder bound obtained using a TM and the corresponding exact error. The Taylor polynomial was developed in the middle of I and the necessary polynomials bounds were computed using a Horner scheme.

The first five examples were analyzed in [START_REF] Chevillard | Efficient and accurate computation of upper bounds of approximation errors[END_REF] for comparing the bounds obtained with interpolation vs. TMs. There, it was observed that in some cases the overestimation in the interpolation remainder is so big, that we can not benefit from using such a polynomial. We used them in order to highlight that CMs do not have this drawback and the remainders obtained with our methods have better quality than the TMs in all situations.

The first example presents a basic function which is analytic on the whole complex plane. There is almost no overestimation in this case, whatever method we use. The second is also a basic function. It has singularities in the complex plane (in π/2 + Zπ), but the interval I is relatively far from these singularities. All the methods present a relatively small overestimation. The third example is the same function but over a larger interval. In these case, the singularities are closer and Taylor polynomials are not very good approximations. The fourth and fifth examples are composite functions on larger intervals. The overestimation in the interpolation method becomes very large, rendering this method useless, while it stays reasonable with TMs and CMs.


Pascal Azerad
email: azerad@math.univ-montp2.fr

Afaf Bouharguane

Jean-Franc ¸ois Crouzet
email: crouzet@math.univ-montp2.fr

Simultaneous denoising and enhancement of signals by a fractal conservation law

Keywords: fractal operator, fractional calculus, Fourier transform, PDEs filters, denoising, enhancement, fast Fourier transform (FFT), finite difference scheme, Savitzky-Golay filter AMS Subject Classifications: 35R11, 60G35, 26A33

In this paper, a new filtering method is presented for simultaneous noise reduction and enhancement of signals using a fractal scalar conservation law which is simply the forward heat equation modified by a fractional anti-diffusive term of lower order. This kind of equation has been first introduced by physicists to describe morphodynamics of sand dunes. To evaluate the performance of this new filter, we perform a number of numerical tests on various signals. Numerical simulations are based on finite difference schemes or Fast Fourier Transform. We used two well-known measuring metrics in signal processing for the comparison. The results indicate that the proposed method outperforms the well-known Savitzky-Golay filter in signal denoising. Interesting multi-scale properties w.r.t. signal frequencies are exhibited allowing to control both denoising and contrast enhancement.

Introduction

Filtering is a process that removes some unwanted component from a signal. It is a very important task in signal processing, data analysis and communication systems. Many techniques have been proposed for this purpose. For instance, we can use simple averaging filters such as a Gaussian filter. It is wellknown that this filter can be realized by solving the heat equation. Other more general partial differential equations (PDE) have been used, with non-linear anisotropic diffusion [START_REF] Catte | Image selective smoothing and edge detection by nonlineal diffusion[END_REF][START_REF] Perona | Scale-space and edge detection using anistropic diffusion[END_REF], non-linear fractional diffusion [START_REF] Bai | Fractional-order anisotropic diffusion for image denoising[END_REF][START_REF] Guidotti | A new nonlocal nonlinear diffusion of image processing[END_REF] or fractional time-derivative [START_REF] Kirane | Anisotropic like approach to image denoising by means of generalized fractional time integrals[END_REF]. It has been proved that PDEs are suitable in signal denoising [START_REF] Nikpour | Using PDE's for Noise Reduction in Time Series[END_REF]. The denoising applications have to take into account two points. First, we want to obtain a clean and readily observable signal (improving signal-to-noise ratio SNR) and secondly, preserve the original shape characteristics (maxima, minima...) of the signal. This task is complex because it is very important that the denoising has no blurring effect on the images and does not change the location of image edges, [START_REF] Buadues | Neighborhood filters and PDE's[END_REF]. For some applications, it is interesting to amplify some features of the signal such as relative maxima or minima in order to enhance its contrast. Usually, these features are flattened by the denoising methods based on averaging techniques. Among the denoising methods which preserve characteristics of the initial signal is the Savitzky-Golay filter, with which we will compare our new method.

A basic and crude idea to enhance the contrast could be to solve the backward heat equation for a few time steps. Of course, this is an ill-posed PDE and we do not advocate this unsafe method but it illustrates the fact that enhancement and denoising are antagonistic operations. Our method is based on a linear PDE, with two antagonistic terms : a usual diffusion and a nonlocal fractional anti-diffusive term of lower order. The diffusion is used to reduce the noise whereas the nonlocal anti-diffusion is used to enhance the contrast. Let us emphasize that we perform at the same time noise reduction and contrast enhancement. Our method is based on the Cauchy problem of the following PDE :

∂ t u(t, x) -a ∂ 2 xx u(t, x) + b I λ [u(t, .)](x) = 0 t ∈ (0, T ), x ∈ R, u(0, x) = u 0 (x) x ∈ R, (1)
where T is any given positive time, u 0 ∈ L 2 (R), a, b are positive constants and I λ is a fractional operator defined as follows via the Fourier transform: for any Schwartz function ϕ ∈ S(R)

I λ [ϕ](x) := -F -1 (|.| λ F(ϕ))(x) (2)
where 0 < λ < 2 and F denote the Fourier transform defined by: for all ξ ∈ R Ff (ξ) := R e -2iπxξ f (x)dx.

Let us note that for λ ∈]1, 2[, we have an explicit nonlocal formula (see proposition 5)

I λ [ϕ](x) = α λ R ϕ ′′ (x -ξ) |ξ| λ-1 dξ (3)
where α λ is a suitable constant. Alternatively, we can also give a slightly different definition, inspired by fractional calculus (see remark below) :

I λ [ϕ](x) = α λ +∞ 0 ϕ ′′ (x -ξ) |ξ| λ-1 dξ. (4)
Remark 1. For causal signals (i.e ϕ(x) = 0 for x < 0) we have

1 Γ(2 -λ) +∞ 0 ϕ ′′ (x -ξ) |ξ| λ-1 dξ = d λ-2 dx λ-2 ϕ ′′ = d λ dx λ ϕ (5)
which is exactly the Riemann Liouville definition of the fractional derivative [START_REF] Podlubny | An introduction to fractional derivatives, fractional differential equations, to methods of their solution and some of their applications[END_REF].

Remark 2. Our model is closely related to a nonlocal conservation law first introduced to describe the morphodynamics of sand dunes and ripples sheared by a fluid flow. Namely, Fowler ([12, 13, 14]) introduced the following equation

∂ t u(t, x) + ∂ x u 2 2 (t, x) + I[u(t, .)](x) -∂ 2 xx u(t, x) = 0, (6)
where u = u(t, x) represents the dune height and I is a nonlocal operator defined as follows: for any Schwartz function ϕ ∈ S(R) and any x ∈ R,

I[ϕ](x) := +∞ 0 |ζ| -1 3 ϕ ′′ (x -ζ)dζ. (7)
Equation ( 6) is valid for a river flow over a erodible bottom u(t, x) with slow variation. The nonlocal term appears after a subtle modeling of the basal shear stress. See [START_REF] Alibaud | A non-monotone nonlocal conservation law for dune morphodynamics[END_REF][START_REF] Alvarez-Samaniego | Existence of travelling-wave and local well-posedness of the Fowler equation[END_REF] for theoretical results on this equation.

This nonlocal term appears also in the work of Kouakou & Lagree [START_REF] Kouakou | Evolution of a model dune in a shear flow[END_REF][START_REF] Lagree | Stability of an erodible bed in various shear flows[END_REF]. The operator I[u] is a weighted mean of second derivatives of u with the bad sign and has therefore an anti-diffusive effect and create instabilities which are controlled by the diffusive operator -∂ 2 xx . We find again this phenomenon for the equation [START_REF] Alibaud | A non-monotone nonlocal conservation law for dune morphodynamics[END_REF].

The solution of the linear PDE ( 1) is obtained by convolution with the kernel K of I λ -∂ 2

xx . Thereafter, this kernel will be our filter for denoising and enhancement of signals. The analysis of this kernel shows that the low frequencies are preserved, the medium frequencies are amplified and the high frequencies are eliminated. It is clear that this kernel depends on the parameters a, b, λ and that the choice of these coefficients will determine the quality of the noise reduction and of the enhancement. In this paper, we discuss the choice of these parameters. To evaluate the performance of our filter, we discretize the PDE by two methods : the fast Fourier transform and the finite difference method. Numerical studies of nonlocal equations are rather scarce: among them, we mention [START_REF] Droniou | A numerical method for fractal conservation laws[END_REF] which proves the convergence of a finite volume method to approximate the solutions of a fractal scalar conservation law, that is to say a conservation law regularized by a diffusive fractional power of the Laplace operator and [START_REF] Azerad | On the stability of a nonlocal conservation law[END_REF] which analyzes the stability of finite difference schemes for the solution of [START_REF] Catte | Image selective smoothing and edge detection by nonlineal diffusion[END_REF]. In these works, the discretization of the fractal operator is performed via an integral formula for I λ similar to the ones appearing in [START_REF] Alibaud | A non-monotone nonlocal conservation law for dune morphodynamics[END_REF][START_REF] Droniou | Global solution and smothing effect for a non-local regularization of an hyperbolic equation[END_REF]. Our method is then compared with the Savitzky-Golay filter. Results show that the PDE (1) is relevant and effective for denoising with preservation or enhancement of features of the signal.

The remaining of this paper is organized as follows: in the next section, we explicit the solution of problem (1) and we give some properties of the kernel. In section 3, we present explicit numerical schemes which approximate the fractal conservation law (1) and we give some numerical simulations. In this section, we also briefly present the Savitzky-Golay filter and compare it with our filter. Other numerical simulations based on the FFT are given in section 4. We also discuss the choice of parameters a, b, λ and we highlight both the ability of denoising and contrast enhancement of our model. Section 5 is devoted to the performance evaluation and metrics. Section 6 gives conclusions about this study.

Theoretical study of the PDEs

In this section, we verify that (1) is well-posed and in subsection 2.2 we analyze the properties of the kernel K of I λ -∂ 2 xx , for 0 < λ < 2.

Well-posedness of the problem

Using the Fourier transform, we see that any solution to (1) satisfies the formula [START_REF] Diethelm | Algorithms for the fractional calculus: a selection of numerical methods[END_REF].

Proposition 1. Let T > 0 and u 0 ∈ L 2 (R). The function u ∈ L ∞ ((0, T ); L 2 (R)
) is a solution of (1) if for any t ∈ (0, T ):

u(t, x) = K(t, •) * u 0 (x) (8)
where K(t, x) = F -1 e -tψ(•) (x) with ψ(ξ) := 4π 2 aξ 2 -b |ξ| λ is the kernel of the operator

I λ -∂ 2 xx .
Proposition 2. Let T > 0, u 0 ∈ L 2 (R). Then, the function

u : t ∈ (0, T ] → K(t, •) * u 0
is well-defined and belongs to C([0, T ]; L 2 (R)): u is extended at t = 0 by the value u(0, .) = u 0 .

Proof. We have

K(t, x) = R e 2iπxξ h(ξ) dξ,
where h(ξ) = e -tψ(ξ) . For 0 < λ < 2, it is easy to verify that h belongs to W 2,1 (R) where

W 2,1 (R) := v ∈ L 1 (R); ∂v ∂x , ∂ 2 v ∂x 2 ∈ L 1 (R) therefore we have that K(t, .) ∈ L 1 (R). Hence, ∀t > 0, K(t, •) * u 0 is in L 2 (R).
Let us prove the strong continuity ie

lim t→0 K (t, .) * u 0 = u 0 in L 2 (R) .
By Plancherel's formula,

||K(t, •) * u 0 -u 0 || 2 L 2 (R) = ||F(K(t, •) * u 0 ) -Fu 0 || 2 L 2 (R) = ||e -tψ Fu 0 -Fu 0 || 2 L 2 (R) = R |e -tψ -1| 2 |Fu 0 | 2 . ( 9
)
Since the function |e -tψ -1| 2 |Fu 0 | 2 converges pointwise to 0 on R, as t → 0 and as min ψ is finite then, by the dominated convergence theorem, the last term of (9) tends to 0 as t → 0.

Remark 3 (Regularity of the solution). It is easy to see

that u ∈ C ∞ ((0, T ] × R) because K is smooth.
The smoothness of K is an immediate consequence of the theorem of derivation under the integral sign applied to the definition of K by Fourier transform. To obtain the regularity at t = 0, we have to suppose that the initial condition u 0 is in C ∞ (R) and satisfies for all k ∈ N, u

(k) 0 ∈ L 2 (R).
The proof is similar with replacing u 0 by u

(k) 0 . Remark 4. From formula (8), it is straightforward that R u(t, x) dx = F (K(t, •)) (0) • R u 0 (x) dx = R u 0 (x) dx. ( 10
)
This is a "mass conservation" property.

Study of the kernel

In this subsection, we give some properties on the kernel

K of I λ -∂ 2 xx .
Proposition 3. The kernel K has a non-zero negative part.

Proof. Let us assume that K is nonnegative, then

|e -tψ(ξ) | ≤ ||F -1 (e -tψ )|| L 1 (R) = R |K(t, .)| = R K(t, .) = F(K(t, .))(0) = e -tψ(0) = 1 for all ξ ∈ R. But we have for 0 < |ξ| < b 4π 2 a 1 2-λ , |e -tψ(ξ) | = e -t(4π 2 aξ 2 -b|ξ| λ ) > 1 , this yields a contradiction.
The main consequence of the non-positivity of K is the failure of the maximum principle for the equation ( 1) [START_REF] Alibaud | A non-monotone nonlocal conservation law for dune morphodynamics[END_REF]. Thereby, u(x, t) is not forced to remain in the interval [min(u 0 ), max(u 0 )]. The signal can then be amplified when it is necessary. Using the proof of the previous proposition, we can deduce that the enhancement of the frequencies will be feasible only for the low/medium frequencies. Of course, this amplification will depend on the choice of the parameters a, b, λ and of the time t. We expose in figure 1 the evolution of K(t, .) for different times and the behaviour of ψ for a = 0.5, b = 2, λ = 1.5.

Finite difference schemes

In this section, we present a finite difference numerical scheme to directly approximate the solution to (1) for any λ ∈]0, 2[.

Integral representation of I λ

To approximate the nonlocal term, it is useful to give an integral formula for I λ : Proposition 4 ( [START_REF] Droniou | Fractal first-order partial differential equations[END_REF][START_REF] Imbert | A non-local regularization of first order Hamilton-Jacobi equations[END_REF]). If λ ∈]0, 2[, we have for all ϕ ∈ S(R), all x ∈ R and all r > 0,

I λ [ϕ](x) = C λ r -r ϕ(x + z) -ϕ(x) -ϕ ′ (x)z |z| 1+λ dz + R\(-r,r) ϕ(x + z) -ϕ(x) |z| 1+λ dz where C λ = λΓ( 1+λ 2 ) 2π 1 2 +λ Γ(1-λ 2 )
and Γ denote the Euler function.

1) If λ ∈]0, 1[ then I λ [ϕ](x) = C λ R ϕ(x + z) -ϕ(x) |z| 1+λ dz. ( 11
)
2

) If λ ∈]1, 2[ then I λ [ϕ](x) = C λ R ϕ(x + z) -ϕ(x) -ϕ ′ (x)z |z| 1+λ dz. (12)
From this proposition, we deduce the following useful result :

Proposition 5. For all ϕ ∈ S(R) and all x ∈ R we have for λ ∈]1, 2[ I λ [ϕ](x) = α λ R ϕ ′′ (x + z) |z| λ-1 dz = α λ R ϕ ′′ (x -z) |z| λ-1 dz (13)
where α λ = C λ λ(λ-1) . Proof. The regularity of ϕ ensures the validity of the following computations . Since

ϕ(x + z) -ϕ(x) -ϕ ′ (x)z = z 0 ϕ ′ (x + y) -ϕ ′ (x) dy = 1 0 ϕ ′ (x + tz) -ϕ ′ (x) z dt,
the last equality arises from the change of variable y = tz. Then, using Fubini's Theorem, we have

I λ (x) = C λ R 1 0 (ϕ ′ (x + tz) -ϕ ′ (x)) z |z| 1+λ dt dz = 1 0 R (ϕ ′ (x + tz) -ϕ ′ (x)) z |z| 1+λ dz dt = y=tz C λ 1 0 R (ϕ ′ (x + y) -ϕ ′ (x)) y |y| 1+λ dy t λ-1 dt = C λ 1 0 t λ-1 dt R (ϕ ′ (x + y) -ϕ ′ (x)) y |y| 1+λ dy = C λ λ R 1 0 ϕ ′′ (x + sy)y ds y |y| 1+λ dy = z=sy C λ λ R ϕ ′′ (x + z) |z| λ-1 dz 1 0 s λ-2 ds = C λ λ(λ -1) R ϕ ′′ (x + z) |z| λ-1 dz.

The numerical scheme

The spatial discretization is given by a set of points x i ; i = 1, ..., N and the discretization in time is represented by a sequence of times t 0 = 0 < ... < t n < ... < T . For the sake of simplicity we will assume constant step size δx and δt in space and time, respectively. The discrete solution at a point will be represented by

u n i ≈ u(t n , x i ).
In this section, we will present the behaviour of explicit numerical scheme which directly approximate the solution to [START_REF] Alibaud | A non-monotone nonlocal conservation law for dune morphodynamics[END_REF].

We discretize all terms of the equation using an explicit method. We consider for any n, j ∈ N

u n+1 j -u n j δt -a u n j+1 -2u n j + u n j-1 δx 2 + b I δx [u n ] j = 0 ( 14
)
where I δx is a discretization of the nonlocal term I λ . This scheme is explicit because the values of the solution at time t n+1 are obtained directly from the (known) values at time t n . For the Laplacian term, we use a standard centered finite difference approximation of second order. To discretize the fractal operator I λ , we consider the formulation (4), which for λ ∈]1, 2[ is a causal variant of [START_REF] Fowler | Evolution equations for dunes and drumlins[END_REF]. Next, we use a basic quadrature rule to approximate the integral and we use a finite difference approximation of the derivative:

I δx [v] j = δx -λ +∞ l=1 l 1-λ (v j-l+1 -2v j-l + v j-l-1 ); (15)
Note that we have absorbed the constant α λ in b.

Remark 5. The practical implementation of the schemes requires to make truncations. The main truncation concerns the integral operator for the nonlocal term I λ . We replace +∞ 0 with L 0 and the finite difference approximation becomes:

I δx [v] j = δx -λ A l=1 l 1-λ (v j-l+1 -2v j-l + v j-l-1 ); (16)
with L = A δx. Usually, we consider that functions are either compactly supported or constant near -∞ and +∞, so it is legitimate to consider a finite sum for the discretization of I λ . However the truncation parameter A has to be chosen judiciously. Indeed, when 0 < λ < 1 , the term l 1-λ in the discretization ( 16) is bigger and bigger for increasing l, hence A ought to be big enough. In contrast, whenever 1 < λ < 2, l 1-λ is negligible for large values and is important only for small values of l. Thus in this case, it is judicious to take A small enough. We see again here that the non local effect is all the more important than λ is small. These memory effects strongly depending on λ have been reported in [START_REF] Diethelm | Algorithms for the fractional calculus: a selection of numerical methods[END_REF]. To take this behaviour into account, we set the truncation parameter A := max(100, 10 λ ). Of course 100 is an arbitrary limit to avoid rounding effects for small λ.

Since the equation is implemented using an explicit method, this imposes a CFL condition on the time and space steps which ensures the numerical stability that is to say that the difference between the approximate solution and the exact solution remains bounded when n → +∞ for δx, δt given. The stability analysis of the nonlocal scheme ( 14) is done in [START_REF] Azerad | On the stability of a nonlocal conservation law[END_REF]. This requires a careful study, owing to the anti-diffusive behaviour of the nonlocal term. Indeed, the analysis of the PDE [START_REF] Alibaud | A non-monotone nonlocal conservation law for dune morphodynamics[END_REF] shows that for the continuous problem, low to medium frequencies are amplified by the nonlocal term, [START_REF] Azerad | On the stability of a nonlocal conservation law[END_REF]. We come back to this in section 4. Therefore, the standard Von Neumann definition of stability must be slightly modified. Denoting θ = k δx ∈ (0, 2π) and u n j = exp (ijθ), instead of imposing that the amplification factor G(θ) = u n+1 j u n j must fulfill |G(θ)| < 1 for any frequency θ, we only impose |G(θ)| < 1 for θ above a given threshold. We obtain the two following conditions on the time and space steps:

(1 -2 1-λ ) b δx λ < 2a δx 2 , ( 17
)
2a δt δx 2 + (2 -2 1-λ ) b δt δx λ < 1. (18)
The first stability condition [START_REF] Kouakou | Evolution of a model dune in a shear flow[END_REF] forces the mesh-size δx to be small enough in order that the diffusion term should dominate the nonlocal anti-diffusive term for high frequencies, whereas the second stability condition (18) looks like an usual CFL condition for explicit schemes and forces the time-step to be small enough. For more details, we refer the reader to [START_REF] Azerad | On the stability of a nonlocal conservation law[END_REF].

Numerical results

In the following numerical tests, we have to take care to choice δt and δx accuratly following the conditions ( 17) and [START_REF] Lagree | Stability of an erodible bed in various shear flows[END_REF]. Thus, the time and space steps depend on the choice of a, b and λ. We begin by consider an electrocardiogram (ECG) signal. Figure 2 illustrates a comparison between the PDE (1) implemented using a finite difference method (FDM) and the Savitzky-Golay smoothing filter. The Savitzky-Golay smoothing filter also called digital smoothing polynomial filter or least-squares smoothing filter was first described in 1964 by Abraham Savitzky and Marcel J. E. Golay. [START_REF] Savitzky | Smoothing and Differentiation of Data by Simplified Least Squares Procedures[END_REF]. Essentially, the Savitzky-Golay method performs a local polynomial regression on a series of values that is to say it replaces each value of the series with a new value which is obtained from a polynomial fit to neighboring points. The algorithm is based on the following equation

u j =   a 0 v j + n+1 2 i=1 a i (v j-i + v j+i )  
where n is the number of data points, a i are (positive) constants, v is the noisy signal and u defines the filtered signal. The main advantage of this approach is the preservation of features of the signal such as relative maxima, minima and width. Indeed, usually these characteristics are 'flattened' by classical averaging techniques such as moving averages. Figure 2 shows three plots: The noisy ECG signal, the smoothed signal (red) using the PDE filter (1) superimposed with the noiseless signal (blue) and the smoothed signal (red) using the Savitzky-Golay filter superimposed with the noiseless signal (blue).

As we can remark, for these parameters, the low/medium frequencies are not amplified but the relative maxima and minima are preserved. The denoising seems correct and the filtering with our numerical scheme and with Savitzky-Golay seems similar, see better. We evaluate our approach and we compare it with Savitzky-Golay in section 5. Figure 3 illustrates two filtered signals of different types. We start from the simplest possible example of nonlocal filter denoising applied to a signal u(x) consisting of a piecewise constant step function v(x) = -1 for 0 < x < 1 and v(x) = 1 for 1 < x < 2 corrupted by additive Gaussian white noise n(x) with standard deviation σ = 0.4: The result of the denoising is illustrated in figure 3(a). As we can notice, the noise is very well eliminated and we find again our original signal, the shape of the signal has been preserved. The result is better than filtering Savitzky-Golay. Moreover, unlike Discrete Fourier Transform where the contrast is low in the neighbourhood of the discontinuity (see figure 5(a)), the dicontinuity/shock is conserved. Therefore, the finite difference method is more suitable for this kind of signal. Figure 3(b) shows another example of good denoising.

u(x) = v(x) + n(x). (19)

Numerical results using Discrete Fourier Transform

In this section, we first discuss the choice of parameters a, b, λ and next, we give some numerical results based on fast Fourier transform (FFT) which is an efficient algorithm to compute the discrete Fourier transform.

Choice of parameters

We fix T = 1. So, we can rewrite the kernel in Fourier variable as follows:

K λ a,b (ξ) = e -ψ(ξ)
where ψ(ξ) = 4π . Therefore ξ 1 is the threshold above which dampening will occur. Now, we wish to control both the denoising and contrast enhancement. For that, we follow a simple strategy:

1. In a first step, we wish to control the frequency range to amplify. This one can be controlled by the following ratio:

ξ 1 ξ M = ( 2 λ ) 1/2-λ .
Indeed, this ratio is decreasing w.r.t. parameter λ and shrinks from +∞ to √ e when λ varies from 0 to 2. Hence, to choose a given amplified frequency range, it is enough to fix parameter λ and the ratio b a (which in turn determines

ξ 1 = b 4π 2 a 1 2-λ ).
2. In this step, we want to monitor the denoising, which will start for frequencies above the neutral one ξ 1 . One can easily check the following equality:

ψ(αξ 1 ) = 4π 2 a ξ 2 1 (α 2 -α λ )
where α is a any positive constant. Therefore, the greater a is, the smaller K λ a,b (αξ 1 ) will be, which means that the dampening rate will increase. Hence parameter a monitors the denoising intensity, this is coherent because it controls the Laplacian term.

3. To finish with, we wish to monitor the value M := K λ a,b (ξ M ) in order to control the contrast enhancement. Indeed, the higher M is, the better the constrast will be amplified. Parameters λ and b a being fixed by item 1 it is enough to adjust both parameters b and a while keeping the ratio b a within some bounds in order to preserve the amplified frequency range. The expression of M shows that to have a good amplification of medium frequencies, it is necessary to increase b as well as a. This is quite natural since the coefficient b controls the anti-diffusive term and has an opposite behaviour to the Laplacian term which tends to flatten the signal. In our method, contrast and denoising are no more antagonistic. We come back on the tuning of parameter b in section 5.

Denoising

In this subsection, we are only going to highlight the ability of noise reduction of our nonlocal equation (1).

To begin with, we consider an ECG signal to illustrate the denoising using the FFT to solve the fractal equation ( 1). The result is given in the figure 5(b). For this choice of parameters, we note that the denoising is suitable and that the relative maxima are preserved. However, the relative minima are not completely preserved, in fig. 5(b) bottom. At last, Figure 10(a) describes the denoising of a sinusoidal signal. We can see the performance of the denoising of nonlocal approach discretized with FFT, which works very well for this type of periodic signal.

Denoising & Enhancement of signals

In this part, we are going to highlight both the ability of noise reduction and contrast enhancement. Of course, to emphasize the amplification, we take into account the remarks done in subsection 4.1.

We start from a simple example of simulatenous denoising and enhancement applied to a strongly attenuated sinusoidal signal highly corrupted by a random noise. In figure 6 middle, we plot the original noiseless signal, the same signal amplified ×50 and the filtered signal, performed with our non local FFT method. As we can notice in figure 6, the noise is eliminated and the contrast is well amplified. Figure 7 illustrates the smoothing of an electrocardiogram signal by filtering the noise with Savitzky-Golay filter (third plot) and by denoising and enhancement with the fractal scalar conservation law (1) (second plot). In the second case, we can see that the relative maxima and minima are amplified. Of course, we can obtain even greater amplification of low/medium frequencies by suitably tuning the parameters but, in these conditions, we will obtain sizeable variations between each peak. Thus, when we wish to amplify the low frequencies, we have to be careful that we do not alter too much the shape of the noiseless signal.

To rate the performance of denoising and enhancement of our filter, we consider signals with different SNR, which corresponds to the power ratio between a signal and the background noise.

In figure 8, we took a small SNR. We can observe that shapes are amplified and that noise has been reduced significantly. But, we can also remark the undershoots just before and behind the shapes. This phenomenon has been highlighted in [START_REF] Alibaud | A non-monotone nonlocal conservation law for dune morphodynamics[END_REF] in the setting of dunes morphodynamics. It is a consequence of the mass conservation property, see equation [START_REF] Droniou | Fractal first-order partial differential equations[END_REF]. We obtain similar results in figure 9 with a large SNR. Indeed, as we can see, the proposed method of filtering yields both an interesting denoising and an amplification of low/medium frequencies. The third plot conveys the filtering using the Savitzky-Golay approach. As agreed upon, we obtain a preservation of relative maxima. Moreover, comparing the output obtained thanks to our nonlocal method with the one of Savitzky-Golay filter, we notice the better ability to increase local extrema (contrast enhancement) while keeping a good denoising. This statement is confirmed by Figure 10(b), where the signal has a medium SNR. We obtain a "good" smoothing and an interesting amplification of the low and medium frequencies.

Performance metrics

In this section, we wish to measure the denoising ability of our model. To evaluate our approach and compare it with the Savitzky-Golay filter, we use two measures: the Mean Square Error (MSE) and the Signal-to-Noise Ratio (SNR). These metrics are frequently used in signal processing. They are defined as follows:

M SE = 1 N N i=1 (u 0 (i) -u(i)) 2 SN R db = 10 log 10 N i=1 u 0 (i) 2 N i=1 (u 0 (i) -u(i))
2 where u 0 is the noiseless original signal, u is the filtered signal and N is the length of the filter. It is easy to see that a small MSE corresponds to a high noise reduction and that a large SNR indicates a good denoising. To compare the performance of filters, we consider a signal of trigonometric type and an ECG signal. Noise is added to these signals with SNR varying between 0 to 8. Results are plotted in figures 10(a) and 5(b). For each signal, we use a sample of 100 random noises. The performance of the two approaches is estimated using MSE and SNR criteria. We plot an average of the results on figures 11 and 12. Figures 11(a Remember that in this case, the algorithm used for the implementation of the equation is the FFT, which is most suitable for trigonometric signals. We may just notice in figure 12 that for high SNR, the Savitzky-Golay filter is better than the nonlocal filter approach. Still, when the SNR is low -and it is the critical case -the proposed method is more efficient than the Savitzky-Golay. Thus, the implementation of our PDE based on the FFT may not be effective for any type of signal, at least when the SNR is high. Alternatively, we may use the finite difference scheme. Results of the filtering are illustrated in figure 13. As we can see, the proposed model is always more efficient than the Savitzky-Golay approach, but for low SNR the FFT scheme remains best. Finally, regardless the type of signal considered, the fractal conservation law (1) is a good tool of denoising, provided that it is implemented with the right method.

Concluding remarks

Our first aim was to introduce a fractal conservation law for denoising and contrast enhancement of signals. This device permits to reduce considerably the noise and to increase contrasts simultaneously. The study showed that our filter eliminates the high frequencies and amplifies the low/medium frequencies. In this paper, we also discussed the choice of parameters a, b and λ. This equation has been implemented using both a finite difference scheme and the fast Fourier transform. Various well-known measuring metrics have been used to compare our method with the Savitzky-Golay filter. Results showed the good performance of our model. Moreover, the analysis highlighted that, depending on the considered signal, it may be more suitable to use the finite difference scheme or the FFT algorithm when the SNR is high. Obviously, for a sinusoidal type signal, it is preferable to use the FFT, whereas for a signal like step functions it is better to implement the fractal equation with finite difference method. Nevertheless, no matter the algorithm used, the fractal consersation law ( 1) is an interesting and natural method for denoising and contrast enhancement.

These satisfying properties for signal processing encourages us to implement it for image enhancement: this study is in progress.

  The kernel of I λ for t = 0.1 (red) and t = 0.5 s (blue)

Figure 1 :

 1 Figure 1: a = 0.5, b = 2 and λ = 1.5

  signal vs. Savitzky-Golay filter

Figure 2 :

 2 Figure 2: Top: a noisy electrocardiogram u; Middle: in red the noiseless signal, in blue the filtered signal with finite difference scheme (14) with a = 4, b = 0.5, λ = 1.7, δx = 1 and δt = 0.1; Bottom : Noiseless signal (blue) vs. filtered signal using Savitzky-Golay filter.

Filtered

  (a) a = 3.5, b = 0.2, λ = 0.1, δt = 10 -6 and δx = 0.001. signal using FDM (b) a = 3.5, b = 0.3, λ = 1.1, δt = 0.1 and δx = 1.

Figure 3 :

 3 Figure 3: (a) Top: Noisy signal (blue) vs. noiseless signal (red); Middle: Noiseless signal (blue) vs. filtered signal using FDM (red); Bottom: Noiseless signal (blue) vs. filtered signal using Savitzky-Golay filter (red). (b) Top: Noisy signal; Middle: Noiseless signal ; Bottom: filtered signal using FDM .

1 (

 1 4π 2 a = 0.01 and λ = 1b) 4π 2 a = 0.01, b = 0.05 and λ = 1.5

Figure 4 :λb 8π 2 a 1 2e 4π 2 a( λb 8π 2 a ) 2 2-λ ( 2 λ - 1 )b 4π 2 a 1 2

 412211 Figure 4: Evolution of the spatial Fourier transform of the kernel K for different parameters.

7 Figure 5 :

 75 Figure 5: Top: Noisy signal (blue) vs. filtered signal (red). Bottom: noiseless signal (blue) vs. filtered signal (red).

  with or without amplification) vs. filtered signal using FFT

Figure 6 :

 6 Figure 6: Top: Noisy signal (blue) vs. noiseless signal v(x) = 0.02 cos(x) (red); Middle: Noiseless signal (blue) vs. filtered signal (red) using FFT with 4π 2 a = 0.2, b = 1.25, λ = 0.5 superimposed with amplified noiseless signal cos(x) (green); Bottom: Noiseless signal (blue) vs. filtered signal (red) using Savitzky-Golay method.

  vs. Savitzky-Golay filter

Figure 7 :

 7 Figure 7: Top: a noisy electrocardiogram u(x); Middle: the noiseless signal v(x) (blue) and denoising & enhancement of the ECG signal (red) using the FFT on the fractal equation (1) with 4π 2 a = 0.005, b = 0.0225 and λ = 1.7. Bottom: filtering using Savitzky-Golay method.

Figure 8 :

 8 Figure 8: First plot (top) represents a noisy signal u(x) (blue). Second plot (bottom) illustrates the noiseless signal v(x) (blue) and denoising & enhancement of the signal (red) using the FFT on the fractal equation (1) with 4π 2 a = 0.1, b = 0.3 and λ = 1.5.

Figure 9 :

 9 Figure 9: Choice of parameters: 4π 2 a = 0.01, b = 0.05 and λ = 1.5. The red line depicts the filtered signal.

Figure 10 :

 10 Figure 10: The red line depicts the filtered signal. In this case, we took u 0 (x) = cos(5πx) + cos(20πx).

  ) value for input signal SNR (db) value for output signal Savitzky-Golay method Proposed method (a) SNR values in Savitzky-Golay method in comparison to proposed method ) value for input signal MSE value for output signal Savitzky-Golay method Proposed method (b) MSE values in Savitzky-Golay method in comparison to proposed method

Figure 11 :

 11 Figure 11: Choice of parameters: 4π 2 a = 0.01, b = 0.03 and λ = 1.5. The initial signal is u 0 (x) = cos(5πx) + cos(20πx).

Figure 12 :

 12 Figure 12: Choice of parameters: 4π 2 a = 0.005, b = 0.015 and λ = 1.7 . The initial signal u 0 (x) is an electrocardiogram (ECG) signal.

  ) value for input signal SNR (db) value for output signal Savitzky-Golay method Proposed method (a) SNR values in Savitzky-Golay method in comparison to proposed method using finite difference scheme ) value for input signal MSE value for output signal Savitzky-Golay filter Proposed method (b) MSE values in Savitzky-Golay method in comparison to proposed method

Figure 13 :Figure 14 :Remark 6 .

 13146 Figure 13: Choice of parameters: a = 0.5, b = 0.3, λ = 1/3. The initial signal u 0 (x) is an electrocardiogram (ECG) signal.

  ) and12(a) show SNR values for these two methods applied to trigonometric

							Noisy signal vs.filtered signal
	2
	1
	0
	-1	0	0.1	0.2	0.3	0.4	0.5	0.6	0.7	0.8	0.9	1
						Noiseless signal vs.filtered signal
	2
	1
	0
	-1	0	0.1	0.2	0.3	0.4	0.5	0.6	0.7	0.8	0.9	1
					Noiseless signal vs. Savitzky-Golay filter
	2
	1
	0
	-1	0	0.1	0.2	0.3	0.4	0.5	0.6	0.7	0.8	0.9	1

Acknowledgements

We thank Bijan Mohammadi for advice on numerical schemes and for helpful comments. P. Azerad and A. Bouharguane are supported by the ANR MATHOCEAN ANR-08-BLAN-0301-02.