Matthieu Garcin

Dominique Guégan

Probability density of the wavelet coefficients of a noisy chaos

Keywords: wavelets, dynamical systems, chaos, noise, alpha-stable

We are interested in the random wavelet coefficients of a noisy signal when this signal is the unidimensional or multidimensional attractor of a chaos. More precisely we give an expression for the probability density of such coefficients. If the noise is a dynamic noise, then our expression is exact. If we face a measurement noise, then we propose two approximations, using Taylor expansion or Edgeworth expansion. We give some illustrations of these theoretical results for the logistic map, the tent map and the Hénon map, perturbed by a Gaussian or a Cauchy noise.

Introduction

Many methods have been proposed to recover a pure linear signal form noisy observations using wavelets. However, when the signal is non-linear, no satisfying solution has been proposed for the denoising problem. Because reality is often non-linear, we aim to propose a method to denoise those signals, in particular when they are chaotic systems.

We model the attractor of a chaotic system (xt) by the relation:

xt = z (xt-1) , (1)
where z is a real function with a bounded support, denoted Supp(z). Moreover, we suppose that this chaotic system is such that each small interval contained in its support is frequently visited by its trajectory. Thanks to that assumption, interpolation is not mandatory when integrating on a large quantity of observed data generated by the chaotic system. Such a chaotic system may be perturbed by some noise: in that case, the observable data is not xt but ut and the observed attractor is not z but a certain z ε defined by: ut = z ε (ut-1) .

The observable data, ut, may be defined by different approaches. Indeed, there exists several ways to model the noise influence. We will consider some kinds of noise, following [START_REF] Guégan | Les Chaos en finance : approche statistique[END_REF]: Measurement noise is a set of independent identically distributed random variables (εt) such that:

ut = xt + εt xt = z (xt-1) . (2)
This kind of noise is frequently used in the literature. It corresponds to a homogeneous perturbation of each observation of the chaotic system. Equation (2) can then be rewritten as: ut = z (ut-1 -εt-1) + εt.

Dynamic noise is a set of independent identically distributed random variables (εt) such that: ut = z (ut-1) + εt.

(3) This kind of noise corresponds to a propagation of small errors. The noisy perception of the image of z feeds the dynamical system.

In the present article, we consider that the noise can be described by a random variable belonging to a wide family of random variables, the symmetric alpha-stable variables, because it includes a lot of well-known random variables and forms a general framework. A symmetric alpha-stable variable X of parameters 0 < α ≤ 2, γ > 0, µ ∈ R, is characterized by its characteristic function:

E e itX = exp (iµt -γ|t| α ) , (4)
where t ∈ R and i is the imaginary unit [START_REF] Samorodnitsky | Stable non-Gaussian random Processes[END_REF] 1 . Since the observed attractor, z ε , is different from z, we are interested in applying a transformation on z ε in order to get a good estimation of z. That crucial step of de-noising can be achieved in different ways. One of them is based on the wavelet analysis and is largely used to de-noise simple signals for which the noise has a linear influence. In the framework of chaotic systems, this approach is more pioneering [START_REF] Guégan | Effect of noise filtering on predictions: on the routes of chaos[END_REF] [START_REF] Guégan | Denoising with wavelets method in chaotic time series: application in climatology, energy and finance in Noise and fluctuations in econophysics and finance[END_REF]. Let us now introduce the wavelet analysis and some of its properties.

We suppose that z ε is a squared-integrable function: z ε ∈ L 2 (R). It is the input signal. We want to decompose it in a discrete dyadic real wavelet basis. Such a wavelet basis is a countable subset of the function space L 2 (R). Each wavelet functions is defined by: ψ j,k : t ∈ R → 2 j/2 Ψ 2 j t -k , where Ψ ∈ L 2 (R) is the real mother wavelet, j ∈ Z is the resolution level and k ∈ Z is the translation parameter [START_REF] Mallat | Une exploration des signaux en ondelettes, Ellipses[END_REF] [START_REF] Rioul | Wavelets and signal processing[END_REF]. We give two simple examples of mother wavelet: the Haar wavelet, defined by:

Ψ Haar : t ∈ R → 1 [0,1/2) (t) -1 [1/2,1) (t);
the Daubechies wavelets extend the previous case and consist in families of wavelets defined on a compact support [START_REF] Daubechies | Ten lectures on wavelets[END_REF]. In general, there is no short formula for such wavelet functions, except for the Haar wavelet which is the simplest case of Daubechies wavelets. However, a simple cascade algorithm allows the construction of the wavelet functions by successive interpolations.

The discrete wavelet analysis consists in decomposing the function z ε in the discrete wavelet basis. For j, k ∈ Z, the projection of z ε on the subspace generated by the vector basis ψ j,k is the wavelet coefficient ẑε j,k of resolution level j and translation parameter k:

ẑε j,k = N n=1
z ε (un)ψ j,k (un)(un -un-1),

1 The alpha-stable distribution can be described by three parameters. If we add a fourth parameter, we can obtain the whole family of alpha-stable random variables and not only the symmetric ones. We are then able to describe a wide range of random variables. Some of them have a thin tail, others are heavy-tailed: it is therefore a much richer class than the only Gaussian random variable which is usually chosen in the literature to describe a noise. In another hand, we note that a rich class of random variables, in order to describe a noise, can be obtained with generalized hyperbolic random variables, which are described by five parameters [START_REF] Barndorff-Nielsen | Exponentially decreasing distributions for the logarithm of particle size[END_REF] [START_REF] Barndorff-Nielsen | Infinite divisibility of the hyperbolic and generalized inverse Gaussian distribution[END_REF]. Some articles make the same choice as us and they deal with a representation of the noise by alpha-stable random variables [START_REF] Georgiou | Alpha-stable modeling of noise and robust time-delay estimation in the presence of impulsive noise[END_REF] [START_REF] Ilow | Signal processing in alpha-stable noise environments: noise modeling, detection and estimation[END_REF] [17][18]. This family of random variables is convenient for the results that follow because, in the proofs, we use characteristic functions. We do not know the probability density function of many alpha-stable random variable but we can describe precisely at least two of them: the Gaussian variable and the Cauchy variable, which are respectively a 2-stable and a 1-stable random variable. where we have considered the discrete grid u1 < u2 < ... < uN of N observations and where u0 = min(inf(Supp(z ε ψ j,k )), u1)

is an additional discretization point.

The wavelet inverse transform consists in reconstructing the observed signal, z ε , according to each translation and resolution level. The reconstructed signal for a given resolution level j ∈ Z is called the detail signal:

Dj : t ∈ R → k∈Z ẑε j,k ψ j,k (t).
The whole reconstructed signal is the sum of all the details:

z ε = j∈Z Dj.
As our goal is to find an approximation for the true chaotic signal z using wavelet decomposition given by ( 5), we propose to filter the wavelet coefficients with a thresholding function, following [START_REF] Guégan | Denoising with wavelets method in chaotic time series: application in climatology, energy and finance in Noise and fluctuations in econophysics and finance[END_REF]. The aim of that thresholding function, Φ, consists in selecting the most significant vectors ψ j,k of the wavelet basis and to eliminate the too noisy terms of the wavelet decomposition. Therefore, for a wavelet coefficient w ∈ R, Φ is defined by: Φ(w) = w1 {|w|≥λ} , like in [START_REF] Donoho | Ideal spatial adaptation via wavelet shrinkage[END_REF], where we have set λ > 0 as the threshold. Then, we define z, the filtered wavelet inverse transform of the noisy attractor, by:

z : t ∈ R → j∈Z k∈Z Φ ẑε j,k ψ j,k (t). (7)
For a given λ > 0, the filtered function z, defined by equation [START_REF] Donoho | Ideal spatial adaptation via wavelet shrinkage[END_REF], is an approximation of z. The aim of the deconvolution with the wavelet method is to select the λ which allows to get the optimal approximation z of z, in the sense that the selected λ is the value which minimizes d(z, z), where d is a function determining the error made by our approximation:

d(z, z) = E R (z(u) -z(u)) 2 du . (8)
In [START_REF] Garcin | Optimal wavelet shrinkage of a noisy chaos[END_REF], we propose such an optimal threshold. In order to calculate it rigorously, we need to know the probability density function of each wavelet coefficient ẑε j,k of the noisy attractor. Indeed, such probability density functions are fundamental when we calculate the expected value d in [START_REF] Garcin | Extreme values of random or chaotic discretization steps and connected networks[END_REF].

Therefore, the present article aims to determine the probability density of each wavelet coefficients ẑε j,k . The importance of such a study for signal processing already arose for simpler frameworks. For example, in [START_REF] Schoukens | Modeling the noise influence on the Fourier coefficients after a discrete Fourier transform[END_REF], Schoukens and Renneboog studied the probability density of each Fourier coefficients of a linear noisy signal, where the noise is represented by Gaussian random variables.

The article is divided into two parts, describing successively the theory (Section 2) and some applications (Section 3).

Main results

In the present section we provide the probability density function of a wavelet coefficient ẑε j,k defined in equation [START_REF] Cornish | Moments and cumulants in the specification of distributions[END_REF] in the presence of both kinds of noise defined respectively in equations ( 3) and ( 2), namely a dynamic noise and a measurement noise. The proofs of the following theorems and propositions are reported at the end of the article.

We consider a sequence of N ∈ N noisy observations: u1, u2, ..., uN , that we rank: u1:N ≤ u2:N , ... ≤ uN:N . For convenience, in the following, we note directly: u1 ≤ u2, ... ≤ uN .

We introduce now some assumptions cited in the following theorems and propositions:

(A0) z is a real function with a bounded support and j, k ∈ Z.

(A1) u1 ≤ u2 ≤ ... ≤ uN ∈ R and u0 is such that u0 = min(inf(Supp(z ε ψ j,k )), u1).

Dynamic noise

We begin with the case of a dynamic noise, introduced in equation ( 3). We introduce some assumptions for this dynamic noise:

(A2) ε1, ε2, ..., εN are N independent identically distributed symmetric alpha-stable random variables of parameters α, γ, µ, where 0 < α ≤ 2, γ > 0, µ ∈ R.

(A2g) ε1, ε2, ..., εN are N independent identically distributed centred Gaussian random variables of variance σ 2 , with σ > 0.

(A2c) ε1, ε2, ..., εN are N independent identically distributed Cauchy random variables with scale parameter γ and location parameter 0, with γ > 0.

(A3) The noisy and observed attractor is defined by: ∀n ∈ {1, ..., N }, z ε (un) = z(un) + εn.

We present, in Theorem 1, a general expression of the random wavelet coefficient of a chaos perturbed by alpha-stable dynamic noise. Then, we study two particular cases according to the nature of the alpha-stable noise: the Gaussian case is studied in Proposition 1 and the Cauchy noise in Proposition 2.

Theorem 1. Let assume (A0), (A1), (A2) and (A3). Then, the wavelet coefficient ẑε j,k of the noisy and observed attractor, z ε , is a symmetric alpha-stable random variable of parameters α, γ , µ , where:

γ = γ N n=1 |ψ j,k (un)(un -un-1)| α µ = ẑj,k + µ N n=1 ψ j,k (un)(un -un-1
). If we are facing a centred Gaussian noise, then every wavelet coefficient of the noisy signal is a Gaussian variable: Proposition 1. Let assume (A0), (A1), (A2g) and (A3). The wavelet coefficient ẑε j,k of the noisy observed attractor, z ε (un), is a Gaussian random variable with mean ẑj,k and variance (σ ) 2 , where:

σ = σ N n=1 |ψ j,k (un)(un -un-1)| 2 .
Moreover, if ψ j,k is an orthonormal wavelet, then σ ≤ σ when N → ∞.

If the noise is described by a Cauchy distribution, then the wavelet coefficients of the noisy signal are also Cauchy variables. This kind of noise is interesting because its fat tail is a good way to model the real perturbations of a signal. The Cauchy distribution example is developed in Proposition 2. We call respectively location parameter and scale parameter the numbers µ ∈ R and γ > 0 such that the Cauchy density function is:

δ Cauchy : x ∈ R → 1 πγ 1 + x-µ γ 2 .
For simplicity, in this article, we only consider Cauchy variables for which µ = 0.

Proposition 2. Let assume (A0), (A1), (A2c) and (A3). The wavelet coefficient ẑε j,k of the noisy observed attractor, z ε (un), is a Cauchy random variable with location parameter ẑj,k and scale parameter γ , where:

γ = γ N n=1 |ψ j,k (un)(un -un-1)|.
Proposition 2 is a direct consequence of Theorem 1 since any Cauchy random variable with scale parameter γ and location parameter 0 is also a symmetric 1-stable random variable with parameters µ = 0 and γ.

One may be surprised by the simplicity of the results obtained for a dynamic noise. The astonishment will grow when one will read in the next section that it is more complicated to obtain de probability density function of the wavelet coefficients in a measurement noise environment. Indeed, in the dynamic noise case, the noisy attractor and the pure one are linked by an affine relation, whereas this link is much intricate in the measurement noise case. Therefore, we remark in [START_REF] Garcin | Optimal wavelet shrinkage of a noisy chaos[END_REF] that it is easier to approximate the attractor by wavelet shrinkage when the noise is dynamic than when it is a measurement noise. On the contrary, it seems obvious that it is more complicated to predict the trajectory of the pure chaos when we are facing a dynamic noise than when we are facing a measurement one.

Measurement noise

We now assume that a measurement noise is perturbing the system, like in equation ( 2). We complete our list of assumptions with the following ones, specific to a measurement noise:

(A2 ) ε1, ε2, ..., εN , ε 1 , ε 2 , ..., ε N are 2N independent identically distributed symmetric alphastable random variables of parameters α, γ, µ, where 0 < α ≤ 2, γ > 0, µ ∈ R.

(A2 g) ε1, ε2, ..., εN , ε 1 , ε 2 , ..., ε N are 2N independent identically distributed centred Gaussian random variables of variance σ 2 , where σ > 0.

(A2 c) ε1, ε2, ..., εN , ε 1 , ε 2 , ..., ε N are 2N independent identically distributed Cauchy random variables with scale parameter γ and location parameter 0, where γ > 0.

(A4) The noise is a measurement noise described by equation (3). More precisely, the noisy and observed attractor is defined by: ∀n ∈ {1, ..., N }, z ε (un) = z(un -ε n ) + εn.

(A5) z is a continuous and piecewise continuously differentiable function on R: ∃M ∈ N, ∃χ0 < ... < χM ∈ R, such that Supp(z) ⊂ [χ0, χM ] and, ∀m ∈ {1, ..., M }, the restriction of z to the interval (χm-1, χm) is differentiable with continuous derivative and z| (χ m-1 ,χm) has a limit on its right in χm-1 and a limit on its left in χm. Moreover N is the set of the non-differentiable points of z and u1, ..., uN ∈ R \ N .

(A6) We only observe states of the dynamical system separated by, at least, two time steps.

(A7) max n∈{1,...,N } |ε n | is small 2 .
We have to pay attention to the fact that, for an observation time t, the observation ut depends both on εt and εt-1, whereas at the next observation time, t + 1, the observation ut+1 depends both on εt+1 and εt. Then, for the next observation time, t + 2, the observation ut+1 depends both on εt+2 and εt+1. Therefore, the independence of the noise implies the independence between ut and ut+2, but consecutive observations, such as ut and ut+1, are not independent. In order to avoid such a situation where some noises are correlated, we assume that we do not observe two consecutive states of the system 3 . But, for convenience, we simply 2 Theoretically, the adjective small is clarified when we describe, in Theorem 3, the error of the approximation of Theorem 2. More practically, the extreme value theory allows to link max n∈{1,...,N }

|ε n | to the scale parameter γ of the noise. For instance, in the example about the logistic map, the error, expressed as the value of the integral of the difference, in absolute value, between the empirical probability density function and the one obtained by the application of Proposition 3, is simply 4.7% for σ = 5%. If σ = 10%, the error becomes about 10.9%. 3 Formally, for an initial observation u ∈ R, we only observe z ν(n) 1≤n≤N -1 , where z n (u) denotes z n-1 (z(u)) and z 1 = z and where ν is a function N → N, such that:

∀n ∈ N, ν(n + 1) -ν(n) ≥ 2.
note u1 ≤ u2, ... ≤ uN these non-consecutive and ranked observations. That condition allows us to assume that ε1, ε2, ..., εN , ε 1 , ε 2 , ..., ε N are independent random variables.

We now present two different approximations of the probability density function of a wavelet coefficient of a chaos with measurement noise. The first one is based on the study of the sensitivity of such a coefficient to small noise. It is particularly suitable when the noise is small. The second method is based on Edgeworth expansion and provides even better results as the number of observation points increases.

Approximation of the noise influence

In this paragraph, we propose first an approximation of the wavelet coefficient distribution of the noisy attractor. This approximation is detailed in Theorem 2 and Propositions 3 and 4. Then, we give an estimation of the error of such an approximation in Theorem 3.

We consider that z is continuous and piecewise continuously differentiable on R, what is the assumption (A5). In fact, such an assumption provides for z a derivative function with a finite limit on the left and the right of each point. Therefore, a first-order Taylor expansion can be achieved because the sensitivity of z(u) to a small variation on u is well defined. That small variation is a consequence of the assumption (A7) made for the noise.

Theorem 2. Let assume (A0), (A1), (A2 ), (A4), (A5), (A6) and (A7). Then, the wavelet coefficient ẑε j,k of the noisy and observed attractor, z ε , can be approximated by a symmetric alpha-stable random variable of parameters α, γ , µ , where:

γ = γ N n=1 (1 + |z (un)| α )|ψ j,k (un)(un -un-1)| α µ = ẑj,k + µ N n=1 (1 -z (un))ψ j,k (un)(un -un-1).
In a practical use, we can have a good approximation thanks to Theorem 2 when we replace the condition u1, ..., uN ∈ R \ N by the condition: N is negligible for the Lebesgue measure, N is big, z is the attractor of an ergodic chaos and its invariant measure has no accumulation point in N . Indeed, the probability of observing a point un in N would be negligible. The same remark can be done for both following propositions, namely Proposition 3 and Proposition 4.

We assume, that the noise is small (A7). In practice, the random noise is unknown since we just observe the noisy data without viewing the clean signal. However, assuming we know the probability density function of the noise, we can reasonably use the results of Theorem 2 as an approximation, when such a probability density function is much localized around zero. It will be the case when the variance is small. But when the variance is not defined, other parameters may be considered, still as an approximation, like the scale parameter for a Cauchy distribution.

If we are facing a centred Gaussian noise, then the wavelet coefficient of the noisy signal is approximated by a Gaussian variable. This is the object of Proposition 3. Proposition 3. Let assume (A0), (A1), (A2 g), (A4), (A5), (A6) and (A7). Then, the limit of the wavelet coefficient ẑε j,k of the noisy and observed attractor, z ε , can be approximated by a Gaussian random variable with mean ẑj,k and variance (σ ) 2 , where:

σ = σ N n=1 (1 + |z (un)| 2 ) |ψ j,k (un)(un -un-1)| 2 .
Proposition 3 is a direct consequence of Theorem 2 since any centred Gaussian random variable of variance σ 2 is also a symmetric 2-stable random variable with parameters µ = 0 and γ = σ 2 /2.

If the noise is described by a Cauchy distribution, then the wavelet coefficient of the noisy signal is approximated by a Cauchy variable. This fat-tailed distribution case is the object of Proposition 4. Proposition 4. Let assume (A0), (A1), (A2 c), (A4), (A5), (A6) and (A7). Then, the limit of the wavelet coefficient ẑε j,k of the noisy and observed attractor, z ε , can be approximated by

γ = γ N n=1 (1 + |z (un)|) |ψ j,k (un)(un -un-1)|.
Proposition 4 is a direct consequence of Theorem 2 since any Cauchy random variable with scale parameter γ and location parameter 0 is also a symmetric 1-stable random variable with parameters µ = 0 and γ.

We see in Proposition 3 and Proposition 4 that the probability density functions of the resulting wavelet coefficients are more spread than in Proposition 1 and Proposition 2 around the wavelet coefficients of the signal without any perturbation by some noise. Indeed, a coefficient |z (un)| 2 or |z (un)| appears in the expression of the variance or of the scale parameter. This is intuitive since the noise only linearly affects z in Proposition 1 and Proposition 2, whereas, in addition to that linear effect, a non-linear noise influence is added in Proposition 3 and Proposition 4.

It can also be interesting to have an idea of the accuracy of the approximation presented in Theorem 2. The error made achieving a Taylor expansion can be controlled. This is the aim of the end of the present paragraph. In Theorem 3, we highlight an upper bound of any quantile of the error made by the approximation suggested in Theorem 2. We restrict here a little the framework to two times continuously differentiable attractors what allow us to have an upper bound of the error thanks to Taylor expansion. Nevertheless, that error given by Taylor expansion is a random variable. In order to get a deterministic upper bound of the error, we will then use extreme value theory. Indeed, the distribution function of the random variables considered is not bounded so that, theoretically, the error can be infinite. However, an upper bound of high quantiles gives a good idea of the error and extreme value theory can furnish such quantiles.

The error estimation, and more particularly the use of extreme value theory, incites us to add two assumptions:

(A8) gN is the cumulative distribution function of N n=1 |ε n | 2 . (A9) The cumulative distribution function of max n∈{1,...,N } |ε n | 2
when N is big exists and is denoted hN .

Theorem 3. Let assume (A0), (A1), (A2 ), (A4), (A5), (A6), (A7), (A8) and (A9). Moreover, z is assumed to be two times continuously differentiable and corresponding to the attractor of an ergodic chaos. We assume that N is Lebesgue measure-zero. Let 0 ≤ p ≤ 1. Let ∆U be the cumulative distribution function of the independent identically distributed observations (un) and let δU be the corresponding probability density function. Let Ẑε j,k be the approximation4 of the wavelet coefficient ẑj,k defined by:

Ẑε j,k = ẑj,k + N n=1 -z (un)ε n + εn ψ j,k (un)(un -un-1).
We define a random distance d j,k between Ẑε j,k and ẑε j,k by:

d j,k = ẑε j,k -Ẑε j,k .
We define D j,k (p), the quantiles of d j,k , by

P [d j,k ≤ D j,k (p)] = p. Then, there exists ΩN , Ω ∞ N : [0, 1] → R such that these quantiles D j,k (p) have an upper bound: D j,k (p) ≤ min (ΩN (p), Ω ∞ N (p)) ,
where

       ΩN (p) = 1 2 max x∈Supp(z) (|z (x)|) max x∈Supp(ψ j,k ) (|ψ j,k (x)|) f -1 N ( √ p)g -1 N ( √ p) Ω ∞ N (p) N →∞ ∼ 1 2 max x∈Supp(z) (|z (x)|) R |ψ j,k (x)| dx h -1 N (p)
and where

fN : v > 0 → max 0, 1 -N sup(Supp(δ U )) inf(Supp(δ U ))+v δU (u) [∆U (u -v) + 1 -∆U (u)] N -1 du .
Theorem 3 gives an upper bound for every quantile of the error made in the approximation of the wavelet coefficient of a chaos with a measurement noise as proposed in Theorem 2. It is therefore a useful tool in order to have a wise and prudent estimation of the accuracy of our approximation since it exaggerates the error. Indeed, the probability that d j,k is bellow ΩN (p) or Ω ∞ N (p) is greater than p. Moreover, it introduces two kinds of upper bounds, one when the number of observations N is finite and another one when that number tends towards infinity. From a practical point of view, Theorem 3 may be more explicit about gN and hN if we are facing a Gaussian noise. This is the aim of Proposition 5. We will further consider the case of a Cauchy noise in Proposition 6. Proposition 5. Let x ≥ 0. Let assume (A2 g), (A8) and (A9). Then:

     gN (x) = γ N 2 , σ 2 x 2 Γ( N 2 ) hN (x) N →∞ ∼ exp -exp α(N ) -x 2σ 2
, where α(N ) is such that:

α(N ) 1 2 e α(N ) = N √ π
and Γ and γ are respectively the Gamma function and the lower incomplete Gamma function 5 .

If the noise is a Cauchy random variable, then we can achieve the same work as in Proposition 5 in order to make some functions of Theorem 3 more explicit. This is the aim of Proposition 6. This result only highlights the case N → +∞. Proposition 6. Let x ≥ 0. Let assume (A2 c) and (A9). Then:

hN (x) N →∞ ∼ exp - 2γN π √ x .

Expansion of the probability density

In the approximation of the noise influence in Theorem 2, we limit our analysis to a local impact of the noise. Indeed, we perform a linear approximation of this influence at a given point of the support of the attractor function. Implicitly, we restrict ourself to a small noise. We could ameliorate our analysis by achieving a second order Taylor expansion, but such a result would still consist on a local analysis of the noise influence. In order to have a global view and to take into account a larger noise, we can realize an expansion of the probability density function. With that approach, the convergence of our estimator is much less dependent on the size of the noise but more on the size of the observation sample. Then, we use the well-known Edgeworth expansion as in [START_REF] Blacher | Multivariate quadratic forms of random vectors[END_REF][5] [START_REF] Lacoume | Statistiques d'ordres supérieurs pour le traitement du signal[END_REF] with the notations adapted to our framework. Let assume (A0), (A1), (A2 ), (A4), (A5), (A6). We assume moreover that all the moments of the 2N random variables representing the noise are well 5 For x, k ≥ 0: defined. Let φN be the probability density function of a random Gaussian variable with mean µN and variance σ 2 N , where:

Γ (k) = +∞ 0 t k-1 e -t dt γ (k, x) = x 0 t k-1 e -t dt.
σ 2 N = γ N n=1 (1 + |z (un)| α )|ψ j,k (un)(un -un-1)| α µN = ẑj,k + µ N n=1 (1 -z (un))ψ j,k (un)(un -un-1).
We note ΦN the corresponding cumulative distribution function. If σ 2 N < 2, then the probability density function, δ ẑε j,k , of the wavelet coefficient ẑε j,k of the noisy and observed attractor is such that, for x ∈ R:

δ ẑε j,k (x) = φN (x) 1 + +∞ m=1 E Hm ẑε j,k -µN √ 2σN Hm x -µN √ 2σN σ m N ( √ 2) m m! , (10)
where (Hm) m∈N are Hermite polynomials 6 . The corresponding cumulative distribution function, ∆ ẑε j,k , of the wavelet coefficient ẑε j,k is such that, for x ∈ R:

∆ ẑε j,k (x) = ΦN (x)+ +∞ m=1 E Hm ẑε j,k -µN √ 2σN Hm-1(0) √ 2σN -φN (x)Hm-1 x -µN √ 2σN σ m+1 N ( √ 2) m-1 m! .
(11) In a practical use of equation [START_REF] Georgiou | Alpha-stable modeling of noise and robust time-delay estimation in the presence of impulsive noise[END_REF], we restrict to truncations of the expansion. We have to pay attention to the fact that such truncated expressions do not strictly define a probability density function and a cumulative distribution function. Therefore, the result will be even more reliable for high-order truncations. It is then recommended to set a stopping criterion and to truncate the sum in equation ( 10) when new terms seem negligible.

In equation [START_REF] Georgiou | Alpha-stable modeling of noise and robust time-delay estimation in the presence of impulsive noise[END_REF], because of the Hermite polynomials, we implicitly need the expression of the moments of the wavelet coefficient of the noisy attractor. These moments cannot be deducted directly from the moments of the attractor value in un, even though we get the wavelet coefficient from a simple linear combination of z ε (u1), ..., z ε (uN ). Instead of that, we propose to calculate jointly the moments and the cumulants. Indeed, cumulants seem to be an appropriate tool 7 and they incite us to propose the following algorithm in order to get the value of E Hm

ẑε j,k -µ N √ 2σ N
in equation [START_REF] Georgiou | Alpha-stable modeling of noise and robust time-delay estimation in the presence of impulsive noise[END_REF] determining the probability density expansion of the wavelet coefficients. More precisely, we are going to define two algorithms: in Algorithm 1, we define a procedure which allows to calculate E Hm

ẑε j,k -µ N √ 2σ N from the cumulants of z ε (un);
considering that the cumulants of z ε (un) are not necessarily easy to get, we define in Algorithm 2 a way to get them from the cumulants of the noise in the case where z can be written as a convergent power series.

The reading of these algorithms is somewhat laborious, but their implementation is very easy.

6 ∀(m, x) ∈ N × R, Hm(x) = (-1) m e x 2 /2 d m
dx m e -x 2 /2 . 7 Indeed, it is easy to get moments from cumulants and cumulants from moments thanks to the fact that C 1 (X) = M 1 (X) and to the recursive formula, for r ∈ {1, 2, ...}:

Cr(X) = Mr(X) - r-1 s=1 r s Cs(X)M r-s (X),
where C denotes any finite cumulant, M any finite moment and X a random variable for which cumulants and moments are well defined. Moreover, two properties concerning cumulants are very interesting:

The cumulant of any sum of independent random variables is the sum of the cumulants of each random variable.

The cumulants are homogeneous: the r-th cumulant of aX, for any a ∈ R, is a r Cr(X).

Algorithm 1. For r ∈ {1, 2, ...}, Cr(X) and Mr(X) are respectively the r-th cumulant and the r-th moment of the random variable X for which cumulants and moments are well defined. We suppose that those cumulants are known for the random variables z ε (un), where n ∈ {1, ..., N }.

Let M ∈ N and µ, σ ∈ R. Then, the following algorithm allows to calculate E Hm

ẑε j,k -µ √ 2σ
for every m ∈ {1, ..., M }, where Hm is the m-th Hermite polynomial:

1. We begin with r = 1. At each iteration, we replace r by r + 1.

2. For each j and k indexing the considered wavelet coefficients, we calculate Cr ẑε j,k as:

Cr ẑε j,k = N n=1
[ψ j,k (un)(un -un-1)] r Cr (z ε (un)).

3. For each j and k indexing the considered wavelet coefficients, we calculate Cr

ẑε j,k -µ √ 2σ
as:

Cr ẑε j,k -µ √ 2σ = Cr ẑε j,k √ 2σ r - Cr (µ) √ 2σ r .
4. For each j and k indexing the considered wavelet coefficients, we calculate Mr

ẑε j,k -µ √ 2σ
as:

   Mr ẑε j,k -µ √ 2σ = Cr ẑε j,k -µ √ 2σ if r = 1 Mr ẑε j,k -µ √ 2σ = Cr ẑε j,k -µ √ 2σ + r-1 s=1 r s Cs ẑε j,k -µ √ 2σ
Mr-s

ẑε j,k -µ √ 2σ
else.

5. We go back to the step 1 of the algorithm until r = M is reached.

6. For each m ∈ {0, ..., M }, we deduct E Hm

ẑε j,k -µ √ 2σ
from the previously calculated moments, say:

               E H0 ẑε j,k -µ √ 2σ = 1 E H1 ẑε j,k -µ √ 2σ = M1 ẑε j,k -µ √ 2σ E H2 ẑε j,k -µ √ 2σ = M2 ẑε j,k -µ √ 2σ -1 E H3 ẑε j,k -µ √ 2σ = M3 ẑε j,k -µ √ 2σ -3M1 ẑε j,k -µ √ 2σ
and so on until the M -th polynomial.

The degree of the m-th order Hermite polynomial is precisely m. Therefore, we only need the cumulants until the order m in order to calculate E Hm

ẑε j,k -µ √ 2σ
. That consideration justify the fifth step of Algorithm 1. We also note that the algorithmic complexity of Algorithm 1, when the size of the sample grows towards infinity, is dominated8 by M (N + ξ), where ξ is the number of wavelet coefficients considered, depending of course on the number of observed data, N .

Algorithm 2 calculates the cumulants of z ε (un) from the cumulants of the noise in the case where z can be written as a convergent power series. Algorithm 2. Let z be a function on R which can be written as a convergent power series around a center c on its support: for each x ∈ Supp(z), ∃(a k ) k∈N such that z(x) = +∞ k=0 a k (x -c) k . Let K ∈ N be the truncation order retained in the computation of z: z(x) ≈ K k=0 a k (x -c) k . Let ε 1 , ε 2 , ..., ε N , ε1, ε2, ..., εN be 2N independent identically distributed random variables with all their cumulants well defined. For r ∈ {1, 2, ...}, Cr(X) and Mr(X) are respectively the r-th cumulant and the r-th moment of the random variable X for which cumulants and moments are well defined. Let R ∈ N. Then, the following algorithm allows to calculate the R first cumulants of z ε (un), where n ∈ {1, ..., N } and z ε (un) = z(un -ε n ) + εn:

1. For each n ∈ {1, ..., N } and each r ∈ {1, ..., KR}, we calculate Mr (un -c -ε n ) as:

Mr (un -c -ε n ) = Cr (un -c)-Cr (ε n )+ r-1 s=1 r s [Cs (un -c) -Cs (ε n )] Mr-s (un -c -ε n ).
2. For each n ∈ {1, ..., N }, each r ∈ {1, ..., R} and each k ∈ {1, ..., K} we calculate Cr (un -c -ε n ) k as:

Cr (un -c -ε n ) k = M kr (un-c-ε n )- r-1 s=1 r s Cs((un -c -ε n ) k )M k(r-s) (un -c -ε n ).
3. For each n ∈ {1, ..., N } and r ∈ {1, ..., R}, we calculate the cumulant of z ε (un) as:

Cr (z ε (un)) ≈ Cr (εn) + K k=0 a r k Cr (un -c -ε n ) k ,
where the approximation stands for the truncation at the order K in the power series.

We can then use successively Algorithm 2 and Algorithm 1 in order to estimate E Hm

ẑε j,k -µ N √ 2σ N
in the case of a convergent power series. The algorithmic complexity of Algorithm 2, when the size of the sample grows towards infinity, is dominated9 by (KR) 2 N . However, that algorithmic cost may be strongly reduced when the noise is a Gaussian variable, while remaining linear in N . Indeed, only two cumulants are then different from zero.

We note that such an approach does not work for every kind of symmetric stable noise because cumulants have to be defined. This is the case, for example, for Gaussian noise but not for Cauchy noise.

We also note that few maps are convergent power series. However, if z is simply a continuous function, then, according to Weierstrass Theorem, it can be uniformly approximated by a polynomial function and therefore by a power series. In practice, one can sometimes get an accurate polynomial function by Lagrange interpolation: the interpolated polynomial is equal to z at pre-specified nodes. If those nodes are equally spaced, one can face Runge's phenomenon: the interpolation is not always better when the number of interpolation nodes increases. To avoid such a polynomial wiggle, one may use Chebyshev nodes. Then, such interpolation polynomials often converge towards the interpolated function when then number of nodes increases [START_REF] Gil | Numerical methods for special functions[END_REF][21] [START_REF] Trefethen | Approximation theory and approximation practice[END_REF]. Therefore, Algorithm 1 and Algorithm 2 allow to apply equation [START_REF] Georgiou | Alpha-stable modeling of noise and robust time-delay estimation in the presence of impulsive noise[END_REF] and equation [START_REF] Gil | Numerical methods for special functions[END_REF] in many practical cases, where z is simply continuous.

Extension to multidimensional chaos

Until this section, we have considered unidimensional chaos. We used therefore unidimensional wavelet basis. In practice we can often face multidimensional chaos. This implies the use of multidimensional wavelet functions. In this new framework, we denote an ergodic chaotic system by the vector Xt = (x1,t, ..., xp,t), where p > 1 is the dimension of the chaos. The attractor of X is a function Z : R p → R p with a bounded support, such that: Xt = Z(Xt-1) = (z1(x1,t-1, ..., xp,t-1), ..., zp(x1,t-1, ..., xp,t-1)) .

Subsequently, we define the noisy attractor, Z ε , of the noisy chaos, Ut:

Ut = Z ε (Ut-1) = z ε
1 (u1,t-1, ..., up,t-1), ..., z ε p (u1,t-1, ..., up,t-1) .

After having described the multidimensional chaos, we have to define the multidimensional wavelet function. We consider that φ is a scaling function and that Ψ is the corresponding mother wavelet generating, by translation and dilatation, an orthonormal basis of the function space L 2 (R). Let θ 0 = φ and θ 1 = Ψ. Then, for each integer η ∈ {1, ..., 2 p -1} written as a binary expression η = η1...ηp, we define the family of 2 p -1 multidimensional mother wavelets by: Ψ η : T = (t1, ..., tp) ∈ R p → θ η 1 (t1) × ... × θ ηp (tp).

By dilating and translating these mother wavelets, with all the parameters j ∈ Z, k = (k1, ..., kp) ∈ Z p and η ∈ {1, ..., 2 p -1}, we get an orthonormal basis of L 2 (R p ) [START_REF] Mallat | Une exploration des signaux en ondelettes, Ellipses[END_REF]. The wavelet functions thus obtained are:

ψ η j,k : T = (t1, ..., tp) ∈ R p → 2 -pj/2 Ψ η 2 j t1 -k1, ..., 2 j tp -kp .
Finally, the discrete wavelet coefficient, of resolution level j ∈ Z and translation parameter k ∈ Z p , of the multidimensional noisy attractor, Z ε , is a vector of R p :

Ẑε η,j,k = N n=1 z ε 1 (Un)ψ η j,k (Un) V (Un) dx, ..., N n=1 z ε p (Un)ψ η j,k (Un) V (Un) dx ,
where, for n ∈ {1, ..., N },

V (Un) = {X ∈ Supp(Z) |∀m = n, X -Un R p ≤ X -Um R p } is
the Voronoi region of one of the N ∈ N observations, . R p being the Euclidean norm in R p . The Voronoi region is a straight way to discretize the integral. Since we had another simple expression for unidimensional chaos, we have not used Voronoi regions in dimension 1. However, the chaos being ergodic, when N tends towards infinity both discretization techniques lead to the same result. We can directly adapt the results developed in the unidimensional environment to multidimensional chaos. Particularly, Theorem 1 and Theorem 2 are slightly transformed to become Theorem 4 and Theorem 5.

Theorem 4. Let 0 < α ≤ 2, γ > 0, µ ∈ R, j ∈ Z, k ∈ Z p , η ∈ {1, ..., 2 p -1}, ( U1
, ..., UN ) ∈ R N . Let ε1,1, ε1,2, ..., ε1,N , ..., εp,1, εp,2, ..., εp,N be pN independent identically distributed symmetric alpha-stable random variables of parameters α, γ, µ. The wavelet coefficient Ẑε η,j,k of the noisy and observed attractor, ∀n ∈ {1, ..., N }, Z ε (Un) = Z(Un) + (ε1,n, ..., εp,n), is a vector whose coordinate indexed by m ∈ {1, ..., p} is a symmetric alpha-stable random variable of parameters α, γm, µm, where:

γm = γ N n=1 |ψ η j,k (Un) V (Un) dx| α µm = Ẑη,j,k m + µ N n=1 ψ η j,k (Un) V (Un) dx.
Theorem 5. Let Z be a continuous and piecewise continuously differentiable function on R p . Let N be the set of the non-differentiable points of

Z. Let 0 < α ≤ 2, γ > 0, µ ∈ R, j ∈ Z, k ∈ Z p , η ∈ {1, ..., 2 p -1}, (U1, ..., UN ) ∈ R N \ N . Let ε1,
γm = γ N n=1 1 + p l=1 (D l Z(Un)) m α |ψ η j,k (Un) V (Un) dx| α µm = Ẑη,j,k m + µ N n=1 1 -p l=1 (D l Z(Un)) m ψ η j,k (Un) V (Un) dx,
where Dm is the differential on the m-th coordinate.

Examples for diverse chaos and noises

The aim of this section is to present the different results for specific cases. We will consider two simple chaos described by the logistic map and the tent map. For each case, we present different approximations of the probability density function of a wavelet coefficient of the noisy chaos and we estimate their accuracy.

The logistic map

The logistic map of parameter α is defined by the recurrence relation:

un+1 = z(un) = αun(1 -un).
If α = 4, then (un) n∈N is a chaotic system. It is a popular example of chaos due to its simple expression [4][12]. We begin with a Gaussian noise and we represent the results in the Figure 1. We have chosen a Daubechies mother wavelet with 5 vanishing moments and the wavelet coefficient ẑε 3,4 , whose wavelet function ψ3,4 has its support roughly contained in the support of z. For a centred Gaussian noise with standard deviation σ = 5% and a sample size N = 129, we obtain a negligible wavelet coefficient for the pure signal: ẑ3,4 = -3.75 × 10 -13 . The impact of the noise on ẑε 3,4 is therefore preponderant as we can observe in the Figure 1. 3,4 of the logistic map of parameter 4 for a centred Gaussian dynamic noise with standard deviation 0.05. If the dynamic noise is replaced by a centred Gaussian measurement noise with standard deviation 0.05, we get the black lines. In both cases, the solid thick line is the empirical density obtained by 1, 000 simulations and represented with a Gaussian kernel, the dotted line is the exact density obtained by Proposition 1 for the dynamic noise and the approximation obtained by Proposition 3 for the measurement noise. The solid thin black line is the approximation of the probability density function obtained by expansion as in equation ( 10) truncated at the 10th order. The mother wavelet is a Daubechies wavelet with 5 vanishing moments. Moreover, the number of observations, N , is 129.

We can also see that both the proposed estimators are very good for the measurement noise as they are close to the empirical distribution. Moreover, in that case, the approximation obtained by Proposition 3, corresponding to the approximation of the noise influence, is also close to the approximation of the probability density function obtained by expansion as in equation [START_REF] Georgiou | Alpha-stable modeling of noise and robust time-delay estimation in the presence of impulsive noise[END_REF]. However, the second approximation is supposed to be more accurate than the first one since it is an expansion from that first approximation. In fact, there is a slight We can also compare, in the Figure 1, the probability density functions for the measurement noise and the dynamic noise. The first one has heavier tails than the second one, what is consistent with the fact that the noise influence is greater in the measurement noise framework.

Gaussian noise has the disadvantage to present thin tails in its probability density function, what is poorly realistic. However, the results developed in the present article allow us to deal with many kinds of noise. In particular, we present in the Figure 3 the probability density function of the wavelet coefficient ẑε 3,4 when we face a Cauchy noise with scale parameter γ = 1%.

In addition to the apparent consistency between theoretic and empirical results in the Figure 3, it seems difficult to build an empirical probability density function in a small time for Cauchy variables. Therefore, Proposition 2 and Proposition 4 provide efficient tools in order to get the exact or an approximated probability density function of the wavelet coefficients of a noisy chaos.

The tent map

The tent map of parameter α is defined by the recurrence relation:

un+1 = z(un) =    αun if 0 < un ≤ 1 2 α (1 -un) if 1 2 < un ≤ 1 0 else.
We consider that α = 2. Then (un) n∈N is a chaotic system [START_REF] Guégan | Les Chaos en finance : approche statistique[END_REF]. We consider a Gaussian noise and we represent the results in the Figure 4. We have chosen the same Daubechies mother wavelet with 5 vanishing moments and the same wavelet coefficient ẑε 3,4 than for the logistic map. For a centred Gaussian noise with standard deviation σ = 5% and a sample size N = 129, we obtain a non-negligible wavelet coefficient for the pure signal: ẑ3,4 = -1.24 × 10 -2 . The impact of the noise on ẑε 3,4 is nevertheless considerable as we can observe in the Figure 1.

For the dynamic noise environment, we have a good fit between the empirical distribution of the wavelet coefficient ẑε 3,4 of the logistic map of parameter 4 for a Cauchy dynamic noise with scale parameter 0.01 and location parameter 0. If the dynamic noise is replaced by a Cauchy measurement noise with scale parameter 0.01 and location parameter 0, we get the black lines. In both cases, the solid thick line is the empirical density obtained by 5, 000 simulations and represented with a Gaussian kernel, the dotted line is the exact density obtained by Proposition 2 for the dynamic noise and the approximation obtained by Proposition 4 for the measurement noise. The mother wavelet is a Daubechies wavelet with 5 vanishing moments. Moreover, the number of observations, N , is 129. measurement noise case is more complex. A translation appears between the approximation obtained by Proposition 3 and the empirical distribution 10 . Thus, we must ameliorate that approximation by proposing an other method. Can it be achieved by a density expansion?

For the expansion as in equation ( 10), we face the problem that the tent map is not a polynomial. Unfortunately, the approximation with Chebyshev interpolation [START_REF] Mason | Chebyshev polynomials[END_REF] seems to fail as it gives unsatisfactory results 11 strongly depending on the choice of parameters of the interpolation. In particular, we observe in Figure 4 a translation of the wavelet coefficient distribution in the measurement noise environment compared with the dynamic noise. This effect (which may be different for other wavelet coefficients) is caused by the non-differentiable point. Since the proposed polynomial interpolation does not approximate accurately z around that point, the translation effect in the approximation of the wavelet coefficient does not appear. We can therefore not use equation [START_REF] Georgiou | Alpha-stable modeling of noise and robust time-delay estimation in the presence of impulsive noise[END_REF] in order to get the probability density function of the wavelet coefficient ẑε 3,4 . However, we are able, for the tent map, to calculate directly the mean of any wavelet coefficient of the map with a measurement noise. Indeed, if the noise is a centred Gaussian variable of variance σ 2 , then, for j, k ∈ Z, φσ the Gaussian density function with variance σ 2 and Φσ the Gaussian cumulative distribution function with the same variance σ 2 :

E ẑε j,k = E N n=1 [z(un -ε n ) + εn] ψ j,k (un)(un -un-1) = N n=1 E [z(un -ε n )] ψ j,k (un)(un -un-1) = N n=1 R z(un -x)φσ(x)dx ψ j,k (un)(un -un-1) = N n=1 un-1 2 un-1 [2 -2(un -x)]φσ(x)dx + un un-1 2 2(un -x)φσ(x)dx ψ j,k (un)(un -un-1) = N n=1 [2 -2un] Φσ(un - 1
2 ) -Φσ(un -1) + 2un Φσ(un) -Φσ(un -1 2 ) +2σ 2 φσ(un) -2φσ(un - 1 2 ) + φσ(un -1) ψ j,k (un)(un -un-1). 10 We remark that the only observation of that empirical distribution, which is almost centred on 0, can lead to the wrong intuition that the wavelet coefficient of the pure signal is negligible. 11 We get a polynomial written like a truncated power series. 3,4 of the tent map of parameter 2 for a centred Gaussian dynamic noise with standard deviation 0.05. If the dynamic noise is replaced by a centred Gaussian measurement noise with standard deviation 0.05, we get the black lines. In both cases, the solid thick line is the empirical density obtained by 1, 000 simulations and represented with a Gaussian kernel, the dotted line is the exact density obtained by Proposition 1 for the dynamic noise and the approximation obtained by Proposition 3 for the measurement noise. That approximation is ameliorated by the approximation represented with the solid thin black which is a Gaussian density with the same standard deviation as the one calculated thanks to Proposition 3 but with the exact mean E ẑε j,k . The mother wavelet is a Daubechies wavelet with 5 vanishing moments. Moreover, the number of observations, N , is 129.

That example of the tent map is a good illustration of the complexity of a chaotic signal with measurement noise: its filtering is more complex than for a signal with a linear noise influence. Therefore, a signal with non-linear noise influence has not to be filtered with the same methods.

The Hénon map

In the following example, we propose a chaos in higher dimension, the Hénon map, which is described by:

xn+1 = yn + 1 -ax 2 n yn+1 = bxn.
If a = 1.4 and b = 0.3, then (xn, yn) n∈N is a chaotic system [START_REF] Guégan | Les Chaos en finance : approche statistique[END_REF], whose phase plot is represented in the Figure 5.

We consider a Gaussian noise and we represent the results in the Figure 6. We have chosen the two-dimensional Daubechies mother wavelet with 5 vanishing moments and the wavelet coefficient j = 3, k = (4, 4) and η = 2: ẑε 2,3,(4,4) . For a centred Gaussian noise with standard deviation σ = 5% and a sample size N = 200, we obtain the following wavelet coefficient for the pure signal12 : ẑ2,3,(4,4) = 2.00 × 10 -3 . We can observe the impact of the noise on ẑ2,3,(4,4) in the Figure 6.

For the dynamic or the measurement noise environment, we have a good fit between the empirical distribution of the wavelet coefficient ẑε 2,3,(4,4) and its exact distribution obtained by Theorem 4 and Theorem 5. Thus the multidimensional case seems well handled by the approximations we proposed. We presented several methods to get the probability distribution of each wavelet coefficients of a noisy chaos. If we face a dynamical noise, Theorem 1 provide an exact expression. For measurement noise, two approximations are proposed in Theorem 2 and Algorithms 1 and 2. Those approximations are of course more accurate when the noise is small and we are able to quantify what small means: the variance of the wavelet coefficient, which depends on the chaos, must be lower than 2. For au Gaussian noise, for the mother wavelet considered and the resolution and translation parameters considered in the article, this corresponds approximately to a variance of the noise equal to 48 for a tent map and equal to 155 for a logistic map.

Moreover, using the multidimensional wavelets, we generalize these results to multidimensional chaos in Theorem 4 and Theorem 5. This allows us to present a quite wide set of examples, from the logistic map to the Hénon map.

These results are very useful to determine the optimal threshold to use when filtering a noisy chaotic attractor to recover the pure chaotic attractor. This is the aim of another article [START_REF] Garcin | Optimal wavelet shrinkage of a noisy chaos[END_REF].

A Proof of Theorem 1

Proof. The wavelet coefficient ẑε j,k of z ε is:

ẑε j,k = N n=1 z ε (un)ψ j,k (un)(un -un-1) = N n=1 [z(un) + εn] ψ j,k (un)(un -un-1) = ẑj,k + N n=1 εnψ j,k (un)(un -un-1). Let t ∈ R. The characteristic function of ẑε j,k is then: E e itẑ ε j,k = exp (itẑ j,k ) E N n=1 exp (itεnψ j,k (un)(un -un-1)) = exp (itẑ j,k ) N n=1 E [exp (itεnψ j,k (un)(un -un-1))],
since the variables ε1, ε2, ..., εN are independent. Then, using the equation ( 4), we get: where:

E e itẑ ε j,k = exp (itẑ j,k ) N n=1 exp (iµtψ j,k (un)(un -un-1) -γ|tψ j,k (un)(un -un-1)| α ) = exp itẑ j,k + N n=1 iµtψ j,k (un)(un -un-1) -N n=1 γ|tψ j,k (un)(un -un-1)| α = exp (iµ t -γ |t| α ) ,
γ = γ N n=1 |ψ j,k (un)(un -un-1)| α µ = ẑj,k + µ N n=1 ψ j,k (
un)(un -un-1). We recognize the characteristic function of a symmetric alpha-stable random variable of parameters α, γ , µ , like in equation ( 4).

B Proof of Proposition 1

Proof. This is a direct consequence of Theorem 1 since any centred Gaussian random variable of variance σ 2 is also a symmetric 2-stable random variable with parameters µ = 0 and γ = σ 2 /2. We also note that, when N → ∞, the maximum discretization step size is decreasing towards 0 because of the ergodic condition on the chaos whose attractor is z. A more precise expression of that maximum step size may be found in [START_REF] Garcin | Extreme values of random or chaotic discretization steps and connected networks[END_REF]. Nevertheless, when N → ∞:

σ N →∞ -→ lim N →∞ σ N n=1 un u n-1 ψ j,k (u)du 2 ≤ lim N →∞ σ N n=1 un u n-1 |ψ j,k (u)| 2 du
thanks to Jensen's inequality. Then:

σ ≤ lim N →∞ σ u N u 0 |ψ j,k (u)| 2 du ≤ σ R |ψ j,k (u)| 2 du = σ ψ j,k L 2 (R) .
Finally, if ψ j,k is an orthonormal wavelet we simply get σ ≤ σ.

C Proof of Theorem 2

Proof. From the equation ( 9), we can achieve a first-order Taylor expansion on z, which is continuous and is continuously differentiable on each segment included in R \ N and which is therefore continuous and continuously differentiable in each un, because they are assumed to be in R \ N : ∀n ∈ {1, ..., N }, z ε (un)

ε n →0 ∼ z(un) -z (un)ε n + εn. (12)
In the equation ( 12), we used the notation a(x) 12), the wavelet coefficient ẑε j,k of z ε is:

ẑε j,k = N n=1 z ε (un)ψ j,k (un)(un -un-1) max n∈{1,...,N } |ε n |→0 ∼ N n=1 [z(un) -z (un)ε n + εn] ψ j,k (un)(un -un-1) max n∈{1,...,N } |ε n |→0 ∼ ẑj,k + N n=1 [-z (un)ε n + εn] ψ j,k (
un)(un -un-1).

Then, the limit of the random variable ẑε j,k when max n∈{1,...,N }

|ε n | → 0, denoted Ẑε j,k , is simply:

Ẑε j,k = ẑj,k + N n=1
-z (un)ε n + εn ψ j,k (un)(un -un-1).

Let t ∈ R. The characteristic function of Ẑε j,k is then: since the variables ε1, ε2, ..., εN , ε 1 , ε 2 , ..., ε N are independent. Then, using the equation ( 4) in the equation ( 14), we get: where:

E e it Ẑε
E e it Ẑε
γ = γ N n=1 (1 + |z (un)| α )|ψ j,k (un)(un -un-1)| α µ = ẑj,k + µ N n=1 (1 -z (un))ψ j,k (
un)(un -un-1). We recognize the characteristic function of a symmetric alpha-stable random variable of parameters α, γ , µ , like in equation ( 4).

D Proof of Theorem 3

Proof. Let r be the remainder of the Taylor-Lagrange expansion of z ε :

∀n ∈ {1, ..., N }, z(un -ε n ) = z(un) -z (un)ε n + r(n),
where

|r(n)| ≤ |ε n | 2 max x∈Supp(z) (|z (x)|) 2 .
Then,

z ε (un) = z(un) -z (un)ε n + εn + r(n)
and we get:

d j,k = ẑε j,k -Ẑε j,k = N n=1 r(n)ψ j,k (
From the equation ( 15), two possible bounds arise. The first one leads to an approximation of the distribution of the error when N tends towards infinity. The second one works particularly well when the number of observations is limited.

1. We can keep the wavelet function in the sum in the equation [START_REF] De Haan | Extreme value theory: an introduction[END_REF] and extract the noise from that sum:

d j,k ≤ 1 2 max x∈Supp(z) z (x) max n∈{1,...,N } |ε n | 2 N n=1 |ψ j,k (un)| (un -un-1). (16)
For an ergodic chaos, when N is big, the biggest step size (un -un-1) is small, so that we can write:

N n=1 |ψ j,k (un)| (un -un-1) N →∞ -→ u N u 0 |ψ j,k (x)| dx ≤ R |ψ j,k (x)| dx. (17)
Then, using together equation ( 16) and equation ( 17), we get:

∀ > 0, ∃N ∈ N such that when N ≥ N d j,k ≤ 1 2 max x∈Supp(z) z (x) M |ε | 2 R |ψ j,k (x)| dx + , (18)
where

M |ε | 2 = lim N →∞ max n∈{1,...,N } |ε n | 2
is a random variable whose distribution may be explicitly known thanks to the extreme value theory.

21

Documents de Travail du Centre d'Economie de la Sorbonne -2013.15

2. On the contrary, we can keep the noise in the sum of the equation ( 15) and extract the wavelet function from that sum:

d j,k ≤ 1 2 max x∈Supp(z) z (x) max x∈Supp(ψ j,k ) (|ψ j,k (x)|) max n∈{1,...,N } (un -un-1)S N,|ε | 2 , (19)
where

S N,|ε | 2 = N n=1 |ε n | 2 .
In equation [START_REF] Lacoume | Statistiques d'ordres supérieurs pour le traitement du signal[END_REF], depending on the cumulative distribution function of the observations, ∆U , we can find the cumulative distribution function of max n∈{1,...,N } (un -un-1). In particular, even when N does not tends towards infinity, we can use the bound presented in [START_REF] Garcin | Extreme values of random or chaotic discretization steps and connected networks[END_REF]:

∆ max n∈{1,...,N } (un-u n-1 ) ≥ max 0, 1 -N sup(Supp(δ U )) inf(Supp(δ U ))+v δU (u) [∆U (u -v) + 1 -∆U (u)] N -1 du , ( 20
)
where δU is the probability density function corresponding to the cumulative distribution function ∆U . In addition to that, we note that if X and Y are positive independent random variables, and x, y > 0, then X ≤ x and Y ≤ y ⇒ XY ≤ xy which leads to:

P(XY ≤ xy) ≥ P(X ≤ x, Y ≤ y) = P(X ≤ x)P(Y ≤ y). (21)
In particular, in equation ( 19), using equation ( 21), we have:

∆ max n∈{1,...,N } (un-u n-1 )S N,|ε | 2 (xy) ≥ ∆ max n∈{1,...,N } (un-u n-1 ) (x)∆S N,|ε | 2 (y). (22)
On one hand thanks to equation [START_REF] Kalluri | Adaptative weighted median filter algorithms for robust signal processing in alpha-stable noise environments[END_REF] where can be arbitrarily small and on the other hand thanks to equations ( 19), ( 20) and ( 22), we can prove the result of Theorem 3, namely:

D j,k (p) ≤ Ω ∞ N (p) D j,k (p) ≤ ΩN (p).

E Proof of Proposition 5

Proof. gN is the cumulative distribution function of

N n=1

ε n σ 2 is a chi-squared random variable with N degrees of freedom. Thus:

gN : x ≥ 0 → γ N 2 , σ 2 x 2 Γ N 2 .
2. For a random variable X, we denote ∆X its cumulative distribution function and δX the probability density function. Since ε 1 is Gaussian, then |ε 1 /σ| 2 is a chi-squared random variable with only one degree of freedom. We denote Y a random variable with the same probability density function than |ε 1 | 2 . Let x > 0. Then: In order to apply some theorem from the extreme value theory, more particularly the von Mises' condition [START_REF] De Haan | Extreme value theory: an introduction[END_REF], we need to find the limit, when x → ∞, of the derivative of the function D, defined by:

δY (x) = 1 σ √ 2πx exp - x 2σ 2 . ( 23
D(x) = 1 -∆Y (x) δY (x) . ( 24
)
Its derivative is:

D (x) = -δY (x) 2 -(1 -∆Y (x))δ Y (x) δY (x) 2 . ( 25
)
From the equation ( 23), we calculate that:

δ Y (x) = -1 2σ √ 2πx 1 x + 1 σ 2 exp -x 2σ 2 x→∞ ∼ - 1 2σ 3 √ 2πx exp -x 2σ 2 . (26)
We now have to find an equivalent, when x → ∞, of 1 -∆Y (x). First, let y > x and I the function defined by:

I(x, y) = y x δY (s)ds. ( 27
)
Integrating the equation ( 27) by parts, we get:

I(x, y) = σ 2 πx exp - x 2σ 2 -σ 2 πy exp - y 2σ 2 - y x 2σ 2 2s 1 σ √ 2πs exp - s 2σ 2 ds, (28)
so that, making y → ∞ in the equation (28):

1 -∆Y (x) = σ 2 πx exp - x 2σ 2 - ∞ x σ 2 s 1 σ √ 2πs exp - s 2σ 2 ds. (29)
Since x, σ > 0, we note that:

∞ x σ 2 s 1 σ √ 2πs exp - s 2σ 2 ds ≤ σ 2 x (1 -∆Y (x)), (30)
so that, using the equations ( 29) and (30), we get the equivalence relation:

1 -∆Y (x) x→∞ ∼ σ 2 πx exp - x 2σ 2 . (31)
Finally, from the equations ( 23), ( 25), ( 26) and (31) we get:

D (x) x→∞ ∼ -1 - σ 2 πx exp -x 2σ 2 - 1 2σ 3 √ 2πx exp -x 2σ 2 1 σ √ 2πx exp -x 2σ 2 2 = 0. ( 32
)
Then, using the equation (32) and von Mises' condition, we can assert that ∆Y is in the max-domain of attraction of a GEV of parameter 0, namely a Gumbel distribution: there exists a real normalization series (aN , bN ) such that the distribution function of

: x -→ exp -e -x . (33)
The parameters of the affine transformation will be determined accurately next in this proof. Moreover, von Mises' condition allows us to write that:

aN N →∞ ∼ ∆ -1 Y 1 -1 N N →∞ -→ +∞ bN N →∞ ∼ D(aN ) N →∞ -→ 2σ 2 , (34)
thanks to the equations ( 23), ( 24) and (31). Let α be a function defined by:

α(N ) 1 2 e α(N ) = N √ π .
One can numerically easily get α but also a proof of its existence and uniqueness, the function t → t 1 2 e t being a bijection of [0, ∞). Then, thanks to the equation (31), we remark that:

1

-∆Y 2σ 2 α(N ) N →∞ ∼ 1 N ,
so that, ∆Y being an increasing function, we can choose to define aN by:

aN = 2σ 2 α(N ). ( 35
)
The limit of the translated and scaled maximum of the set of random variables can now be written. Since bN x = N →∞ o (aN ), then:

P max n∈{1,...,N } |ε n | 2 -a N b N ≤ x = [∆Y (xbN + aN )] N N →∞ ∼ 1 -(1 -∆Y (aN )) exp -xb N 2σ 2 N N →∞ ∼ 1 -1 N exp (-x) N N →∞ -→ G0(x). (36)
Finally, using the equations ( 33), (34), ( 35) and (36), we get:

hN (x) N →∞ ∼ G0 x-a N b N N →∞ ∼ exp -exp α(N ) -x 2σ 2
.

We remark that successive Taylor expansions using the equation (31) can lead to build α as the limit of the sequence of functions (α k ) k∈N , defined by:

     α0(N ) = log N √ π α k+1 (N ) = log N √ πα k (N ) . ( 37
)
Then α is the limit of α k , when k → ∞. Indeed, when we replace α k+1 and α k by α in the equation (37), we find that α is such that:

α(N ) 1 2 e α(N ) = N √ π ,
what is consistent with its definition in Proposition 5.

F Proof of Proposition 6

Proof. Let x ≥ 0. We denote Y a random variable with the same probability density function than |ε 1 | 2 . We remark that:

∆Y (x) = P (|ε 1 | ≤ √ x) = ∆ ε 1 ( √ x) -∆ ε 1 (- √ x) = 2 π arctan √ x γ
, because the cumulative distribution function of a Cauchy random variable with scale parameter γ is:

∆ ε 1 (x) = 1 2 + 1 π arctan x γ .
24 Documents de Travail du Centre d'Economie de la Sorbonne -2013.15

Therefore, the density probability function of Y is:

δY (x) = ∂∆ ε 1 ( √ x) ∂x - ∂∆ ε 1 (- √ x) ∂x = ∂ ∂x 2 π arctan √ x γ = 1 γπ √ x 1+ x γ 2 x→∞ ∼ γ πx 3/2 , (38)
Like for Proposition 5, we are going to use the von Mises' condition: we look for the limit, when x → ∞, of the derivative of the function D, defined by:

D(x) = 1 -∆Y (x) δY (x) .
Its derivative is:

D (x) = -δY (x) 2 -(1 -∆Y (x))δ Y (x) δY (x) 2 . ( 39
)
From the equation (38), we calculate that:

δ Y (x) = -
G 1/2 : x -→ exp -x -1/2 . ( 42
)
The parameters of the affine transformation will be determined accurately next in this proof. Moreover, von Mises' condition allows us to write that:

aN N →∞ ∼ ∆ -1 Y 1 -1 N bN N →∞ ∼ D(aN ). ( 43
)
The equation (41) incites us to choose:

aN = 4γ 2 N 2 π 2 , (44)
from which we get, using the equations ( 38) and ( 43 The limit of the translated and scaled maximum of the set of random variables can now be written, using the equations (41), ( 44) and ( 45 (46)

Finally, using the equations ( 42), ( 44), ( 45) and (46), we get:

hN (x) N →∞ ∼ G 1/2 2 x-a N b N + 1 N →∞ ∼ exp -2γN π √
x .

G Proof of Theorem 4

Proof. The wavelet coefficient Ẑε η,j,k of Z ε is: + µ N n=1 ψ η j,k (Un) V (Un) dx. We recognize the characteristic function of a symmetric alpha-stable random variable of parameters α, γm, µm, like in equation ( 4).

H Proof of Theorem 5

Proof. Like in the proof of Theorem 2, we achieve a first-order Taylor expansion on Z: ∀n ∈ {1, ..., N }, Z ε (Un) We recognize the characteristic function of a symmetric alpha-stable random variable of parameters α, γm, µm, like in equation ( 4).

I Chebyshev interpolation polynomials

Chebyshev polynomials of the first kind are described by:

Tn(x) =

2

  Documents de Travail du Centre d'Economie de la Sorbonne -2013.15

8

  Documents de Travail du Centre d'Economie de la Sorbonne -2013.15

Figure 1 :

 1 Figure 1: In grey, probability density function of the wavelet coefficient ẑε

13

  Documents de Travail du Centre d'Economie de laSorbonne -2013.15 difference between both densities, represented in Figure2: in the first approximation, the tail are slightly underestimated.

Figure 2 :

 2 Figure 2: Difference between the approximations of the probability density functions for the measurement noise framework represented in Figure 1: the one obtained by Proposition 3 minus the one obtained by equation (10).

3 , 4

 34 and its exact distribution obtained by Proposition 1. The 14 Documents de Travail du Centre d'Economie de la Sorbonne -2013.15

Figure 3 :

 3 Figure 3: In grey, probability density function of the wavelet coefficient ẑε

See appendix I 15 Documents

 15 de Travail du Centre d'Economie de la Sorbonne -2013.15

Figure 4 :

 4 Figure 4: In grey, probability density function of the wavelet coefficient ẑε

Figure 5 :

 5 Figure 5: Phase plot of the chaotic Hénon map with a = 1.4 and b = 0.3.

Figure 6 :

 6 Figure 6: In grey, probability density function of the wavelet coefficient ẑε

2, 3 ,

 3 [START_REF] Breeden | Noise in the modeling and control of dynamical systems[END_REF][START_REF] Breeden | Noise in the modeling and control of dynamical systems[END_REF] of the Hénon map of parameters a = 1.4 and b = 0.3 for a centred Gaussian dynamic noise with standard deviation 0.05. If the dynamic noise is replaced by a centred Gaussian measurement noise with standard deviation 0.05, we get the black lines. In both cases, the solid thick line is the empirical density obtained by 1, 000 simulations and represented with a Gaussian kernel, the dotted line is the exact density obtained by Theorem 4 for the dynamic noise and the approximation obtained by Proposition 5 for the measurement noise. The mother wavelet is a Daubechies wavelet with 5 vanishing moments. Moreover, the number of observations, N , is 200.

19

  Documents de Travail du Centre d'Economie de la Sorbonne -2013.15

x→0∼

  b(x), for functions a and b. This means that, when lim x→0 b(x) = 0, then lim x→0 a(x) b(x) = 1. Using the equation (

  j,k = exp (itẑ j,k ) E N n=1 exp (itεnψ j,k (un)(un -un-1)) exp (-itz (un)ε n ψ j,k (un)(un -un-1)) = exp (itẑ j,k ) N n=1 E [exp (itεnψ j,k (un)(un -un-1))] E [exp (-itz (un)ε n ψ j,k (un)(un -un-1))],(14) 20 Documents de Travail du Centre d'Economie de la Sorbonne -2013.15

  itẑ j,k ) N n=1 exp (iµtψ j,k (un)(un -un-1) -γ|tψ j,k (un)(un -un-1)| α ) N n=1 exp (-iµtz (un)ψ j,k (un)(un -un-1) -γ|z (un)| α |tψ j,k (un)(un -un-1)| α ) = exp itẑ j,k + N n=1 iµt(1 -z (un))ψ j,k (un)(un -un-1) -N n=1 γ(1 + |z (un)| α )|tψ j,k (un)(un -un-1)| α = exp (iµ t -γ |t| α ) ,

  un)(un -un-1) x)|) N n=1 |ε n | 2 |ψ j,k (un)(un -un-1)|.

N n=1 |ε n | 2

 2 and hN is the cumulative distribution function of limN →∞ max n∈{1,...,N } |ε n | 2 .1.

) 22 Documents

 22 de Travail du Centre d'Economie de la Sorbonne -2013.15

the random variable b - 1 N

 1 max n∈{1,...,N } |ε n | 2 -aN converges towards the distribution function of an affine transformation of a Gumbel random variable, whose distribution function is: G0

2 = 2 .

 22 the equations (38), (39), (40) and (41) we get an equivalent of the derivative of D: According to von Mises' condition, we can assert that ∆Y is in the max-domain of attraction of a GEV of parameter 2, namely a Fréchet distribution: there exists a real normalization series (aN , bN ) such that the distribution function of the random variable b -1 N max n∈{1,...,N } |ε n | 2 -aN converges towards the distribution function of an affine transformation of a Fréchet random variable of parameter 1/2, whose distribution function is:

  ): bN = 8γ 2 N 2 π 2 . (45) 25 Documents de Travail du Centre d'Economie de la Sorbonne -2013.15

  ..,N } |ε n | 2 -a N

  ε (Un)ψ η j,k (Un) V (Un) dx = N n=1 [Z(Un) + (ε1,n, ..., εp,n)] ψ η j,k (Un) V (Un) dx = Ẑη,j,k + N n=1 ε1,nψ η j,k (Un) V (Un) dx, ..., N n=1 ε1,nψ η j,k (Un) V (Un) dx .Let t ∈ R. The characteristic function of the m-th coordinate of Ẑε η,j,k is then:E exp it Ẑε η,j,k m = exp it Ẑη,j,k m E N n=1 exp itεm,nψ η j,k (Un) V (Un) dx = exp it Ẑη,j,k m N n=1 E exp itεm,nψ η j,k (Un) V (Un) dx ,since the variables εm,1, εm,2, ..., εm,N are independent. Then, using the equation (4), we get:E exp it Ẑε η,j,k m = exp it Ẑη,j,k m N n=1 exp iµtψ η j,k (Un) V (Un) dx -γ|tψ η j,k (Un) V (Un) dx| α = exp it Ẑη,j,k m + N n=1 iµtψ η j,k (Un) V (Un) dx -N n=1 γ|tψ η j,k (Un) V (Un) dx| α = exp (iµmt -γm|t| α ) ,where: γm = γ N n=1 |ψ η j,k (Un) V (Un) dx| α µm = Ẑη,j,k m

ε 1 ,D

 1 n ,...,ε p,n →0 ∼ Z(Un) -p m=1 DmZ(Un)ε m,n + (ε1,n, ..., εp,n) . (47) 26 Documents de Travail du Centre d'Economie de la Sorbonne -2013.15Using the equation (47), the wavelet coefficient Ẑε η,j,k of Z ε is:Ẑε η,j,k = N n=1 Z ε (Un)ψ η j,k (Un) V (Un) dx max n∈{1,...,N } |ε 1,n |+...+|ε p,n |→0 ∼ N n=1 Z(Un) -p m=1 DmZ(Un)ε m,n + (ε1,n, ..., εp,n) ψ η j,k (Un) V (Un) dx max n∈{1,...,N } |ε 1,n |+...+|ε p,n |→0 ∼ Ẑη,j,k + N n=1 -p m=1 DmZ(Un)ε m,n + (ε1,n, ..., εp,n) ψ η j,k (Un) V (Un) dx.Then, the limit of the random variable Ẑε η,j,k when maxn∈{1,...,N } |ε 1,n | + ... + |ε p,n | → 0, denoted Ẑε,lim η,j,k , is simply: Ẑε,lim η,j,k = Ẑη,j,k + l Z(Un)ε l,n + (ε1,n, ..., εp,n) ψ η j,k (Un) V (Un)dx.Let t ∈ R. The characteristic function of the m-th coordinate of Ẑε,lim η,j,k is then:E exp it Ẑε,lim η,j,k m = exp it Ẑη,j,k m N n=1 E exp itεm,nψ η j,k (Un) V (Un) dx N n=1 E exp -it p l=1 (D l Z(Un)) m ε l,n ψ η j,k (Un) V (Un) dx = exp it Ẑη,j,k m N n=1 E exp itεm,nψ η j,k (Un) V (Un) dx N n=1 p l=1 E exp -it (D l Z(Un)) m ε l,n ψ η j,k (Un) V (Un) dx ,(48)since the variables ε 1,1 , ε 1,2 , ..., ε 1,N , ..., ε p,1 , ε p,2 , ..., ε p,N , εm,1, εm,2, ..., εm,N are independent. Then, using the equation (4) in the equation (48), we get:E exp it Ẑε,lim η,j,k m = exp it Ẑη,j,k m N n=1 exp iµtψ η j,k (Un) V (Un) dx -γ|tψ η j,k (Un) V (Un) dx| α N n=1 p l=1 exp -iµt (D l Z(Un)) m ψ η j,k (Un) V (Un) dx N n=1 p l=1 exp -γ| (D l Z(Un)) m | α |tψ η j,k (Un) V (Un) dx| α = exp it Ẑη,j,k m + N n=1 iµt 1 -p l=1 (D l Z(Un)) m ψ η j,k (Un) V (Un) dx exp -N n=1 γ 1 + p l=1 (D l Z(Un)) m α |tψ η j,k(Un) V (Un) dx| α = exp (iµmt -γm|t| α ) , where: γm = γ N n=1 1 + p l=1 (D l Z(Un)) m α |ψ η j,k (Un) V (Un) dx| α µm = Ẑη,j,k m + µ N n=1 1 -p l=1 (D l Z(Un)) m ψ η j,k (Un) V (Un) dx.

n 2 n 2 k=0 (- 1 )

 21 k (n-k-1)! k!(n-2k)! (2x) n-2k if n ≥ 1 1 else,where x ∈ R. The roots of Tn+1 are noted x0, ..., xn and they are defined by:xj = cos j + 1 2 π n + 1.

Documents de Travail du Centre d'Economie de la Sorbonne -2013.15

It is the approximation appearing in equation (13) in the proof of Theorem 2 7 Documents de Travail du Centre d'Economie de la Sorbonne -2013.15

If M = 10, N = 1000 and ξ < N , the execution time is about 0.5 second for a classical interpreted language running on an Intel Core 2 Duo 3GHz. 10 Documents de Travail du Centre d'Economie de la Sorbonne -2013.15

If K =

10, R = 10 and N = 1000, the execution time is about 45 seconds for a classical interpreted language running on an Intel Core 2 Duo 3GHz.

This empirical value is obtained by taking arbitrarily the same Voronoi region size for all the observations.

We want to interpolate the function z, knowing its value at some nodes which are the x0, ..., xn translated and scaled. More precisely, the nodes are taken on an interval (a, b) wider than the support of z. We note them ζ(x0), ..., ζ(xn), where ζ :

x is an affine function. Because of the discrete orthogonality of the Chebyshev polynomials, [START_REF] Mason | Chebyshev polynomials[END_REF], one gets Pn, the Chebyshev interpolation polynomial of z:

and

else.

Then

Therefore

what may also be written with the power series form with a center c = a+b 2 :

n-1 2 k=j (-1) k-j 2 4j+1 (b-a) 2j+1 c 2k+1 (2k + 1) (k+j)! (k-j)!(2j+1)! (x -c) 2j+1 .