Mamadou Mboup

Cédric Join
email: cedric.join@cran.uhp-nancy.fr

Michel Fliess
email: michel.fliess@polytechnique.edu

Numerical differentiation with annihilators in noisy environment

Keywords: Numerical differentiation, Operational calculus, Orthogonal polynomials, Linear filtering Mathematics Subject Classification (2000) 65D25, 44A40, 44A10

Numerical differentiation in noisy environment is revised through an algebraic approach. For each given order, an explicit formula yielding a pointwise derivative estimation is derived, using elementary differential algebraic operations. These expressions are composed of iterated integrals of the noisy observation signal. We show in particular that the introduction of delayed estimates affords significant improvement. An implementation in terms of a classical finite impulse response (FIR) digital filter is given. Several simulation results are presented.

Introduction

Numerical differentiation, i.e., the derivative estimation of noisy time signals, is a longstanding difficult ill-posed problem in numerical analysis and in signal processing and control. It has attracted a lot of attention due to its importance in many fields of engineering and applied mathematics. A number of different approaches have been proposed. Methods based on observer design may be found in the control literature [START_REF] Dabroom | Discrete-time implementation of high-gain observers for numerical differentiation[END_REF], [START_REF] Chitour | Time-varying high-gain observers for numerical differentiation[END_REF], [START_REF] Ibrir | Online exact differentiation and notion of asymptotic algebraic observers[END_REF]. See also [START_REF] Ibrir | Linear time-derivatives trackers[END_REF][START_REF] Levant | Higher-order sliding modes, differentiation and output-feedback control[END_REF][START_REF] Su | A simple improved velocity estimation for low-speed regions based on position measurements only[END_REF] for other approaches from the control litterature. In signal processing, it is very common to cast the problem in terms of frequency domain digital filter design [START_REF] Roberts | Digital signal processing[END_REF], [START_REF] Rader | Approximating noncausal IIR digital filters having arbitrary poles, including new Hilbert transformer designs, via forward/backward block recursion[END_REF], [START_REF] Chen | Design of high-order digital differentiators using L 1 error criteria[END_REF]. This is motivated by the observation that an ideal n th order differentiator has a frequency response of magnitude ω n . Another interesting approach, due to [START_REF] Al-Alaoui | A class of second-order integrators and low-pass differentiators[END_REF], consists in inverting the minimum phase transfer function of a properly designed numerical integrator to obtain a digital differentiator. Similar ideas have also been presented by [START_REF] Diop | On numerical differentiation algorithms for nonlinear estimation[END_REF], in the continuous-time context. All these interesting approaches have been developped as candidate alternatives to the very classical one, based on least-squares polynomial fitting or (spline) interpolation. When completed with a regularization step, this classical approach may be very efficient (see i.e. [START_REF] Ibrir | A numerical procedure for filtering and efficient high-order signal differentiation[END_REF], [START_REF] Duncan | Numerical differentiation and parameter estimation in higher-order linear stochastic systems[END_REF]) in off-line applications.

In this paper, the numerical differentiation problem is revised through the algebraic framework of parameter estimation initially presented in [START_REF] Fliess | An algebraic framework for linear identification[END_REF] (see also [START_REF] Fliess | Questioning some paradigms of signal processing via concrete examples[END_REF] and [START_REF] Mboup | Parameter estimation via differential algebra and operational calculus[END_REF]). Our purpose is to improve a new approach which started in [START_REF] Fliess | Compression différentielle de transitoires bruités[END_REF][START_REF] Fliess | Control via state estimations of some nonlinear systems[END_REF], and in [START_REF] Fliess | Estimation des dérivées d'un signal multidimensionnel avec applications aux images et aux vidéos[END_REF][START_REF] Fliess | Analyse et représentation de signaux transitoires : application à la compression, au débruitage et à la détection de ruptures[END_REF][START_REF] Fliess | Closed-loop fault-tolerant control for uncertain nonlinear systems[END_REF], for solving various questions in control and in signal and image processing. Given a smooth signal x and an order n, a key point of our approach is to consider d n x(t) dt n | t=τ , for each fixed τ 0, as a single parameter to be estimated from a noisy obervation of the signal. A pointwise derivative estimation therefore follows by varying τ . The general ideas leading to this estimation are exposed in section 2. As in the classical approaches, the starting point is an order N n truncation of the Taylor series expansion of the signal. One of the key feature of the proposed approach is to operate in the operational calculus domain [START_REF] Yosida | Operational Calculus -A Theory of Hyperfunctions[END_REF], [START_REF] Mikusiǹski | Operational Calculus[END_REF][START_REF] Mikusiǹski | Operational Calculus[END_REF]. Therein, an extensive use of differential elimination and a series of algebraic manipulations yield, back in the time domain, an explicit expression for the estimate of d n x (t) dt n | t=τ as an integral operator of the noisy observation within a short time interval [τ, τ +T ]. The differential algebraic manipulations may be chosen in such a way that the time domain integral operators specialise to iterated integrals. This corresponds to a family of pointwise derivative estimators introduced in section 3. This section contains the main contributions of the paper. Therein, we establish a direct link in between our algebraic numerical differentiators with the well known fact that a scalar product of the signal x(t) with a kernel having n vanishing moments essentially captures its n th order derivative, provided x is smooth. Note however that drawing an efficient numerical differentiator from this general principle is another history, especcially in a noisy setting. More precisely, recall that (see e.g. [START_REF] Alpay | Algorithme de Schur, espaces à noyau reproduisant et théorie des systèmes[END_REF]Theorem 2.3.5,p. 31]): If f (t) belongs to a reproducing kernel Hilbert space H, with reproducing kernel K H (τ, t), analytic in both t and τ * , then we have

f (t), ∂ ℓ ∂τ * ℓ K H (τ, t) = f (ℓ) (τ ). (1)
We mention that, the model space H and especially the evaluation point τ have to be judiciously chosen before the equation ( 1) yields some efficient derivative estimator. But, once such choices are met, high performance numerical differentiation method is obtained. Now, it is a happy fact that the proposed algebraic estimators inherently achieve these judicious choices. In particular, the estimators are as efficient as an appropriate time delay is admitted. We give an explicit expression for the corresconding delay. Section 4 tells us how to compute the estimates based on the samples of the noisy observation signal. For equally spaced samples, our numerical differentiators take the form of classical finite impulse response (FIR) digital filters. These filters are next used to present some simulation examples. To ease the presentation, all the proofs are deferred to an appendix.

Preliminaries

We consider the estimation of x (n) (t), the n th order derivative of a smooth signal x(t) defined on an interval I ⊂ R + . In a concrete application, the signal x(t) is not directly available but rather it is observed through a noise corruption. In all the sequel, ̟(t) will denote the noise and we set y(t) = x(t) + ̟(t) for the noisy observation. The problem is challenging since it is well known that differentiation tends to amplify the noise.

Operational domain

Let us ignore the noise for a moment. Assume that x(t) is analytic on I so that we may consider, without any loss of generality, the convergent Taylor expansion at t = 0

x(t) = i≥0 x (i) (0) t i i! . (2)
With N n, the truncated Taylor expansion

x N (t) = N i=0 x (i) (0) t i i! , (3)
satisfies the differential equation d N +1 dt N +1 x N (t) = 0. In the operational domain, this reads as [START_REF] Aronszajn | Theory of reproducing kernels[END_REF] where xN (s) stands for the operational analog of x N (t). In all the sequel, the operational analog of a signal u(t) will be denoted as û(s) and, to ease the notation to follow, we will drop the argument s and write û for short.

s N +1 xN (s) = s N x(0) + s N -1 ẋ(0) + • • • + s N -n x (n) (0) + • • • + x (N ) (0),
The basic step towards the estimation of x (n) (t), t 0 is the estimation of the coefficient x (n) (0) from the observation y(t). Indeed, assume that an algorithm is already designed for this task. Then, for any τ 0, it is clear that applying the same algorithm on the new signal observation y τ (t) △ = Heaviside(t)y(t + τ ), will yield an estimate for x (n) (τ ). We thus lay the focus on how to estimate x (n) (0) in (4).

Annihilators

Of course it is also possible to estimate all the coefficients x (i) (0), i = 0, . . . , N simultaneously. However, not only the coefficients x (i) (0), i = n are not necessary for the estimation of x (n) (t) as explaned above, but also simultaneous estimation is more sensitive to noise and numerical computation errors [START_REF] Mboup | A revised look at numerical differentiation with an application to nonlinear feedback control[END_REF]. These drawbacks are avoided in the proposed approach. All the terms s N -i x (i) (0) in ( 4) with i = n, are consequently considered as undesired perturbations which we proceed to annihilate. For this, it suffices to find a linear differential operator, i.e.

Π = finite ̺ ℓ (s) d ℓ ds ℓ , ̺ ℓ (s) ∈ C(s), (5)
satisfying

Π xN = ̺(s)x (n) (0), (6)
for some rational function ̺(s) ∈ C(s). Such a linear differential operator, subsequently called an annihilator for x (n) (0), obviously exists and is not unique.

It is also clear that to each annihilator Π, there is a unique ̺(s) ∈ C(s) such that (6) holds. We shall say that Π and ̺(s) are associated.

Integral operators

A linear differential operator is said to be proper (resp. strictly proper) if, and only if, each ̺ ℓ (s) in ( 5) is a proper (resp. strictly proper) rational function.

We recall that a rational function a b is said to be proper (resp. strictly proper) if, and only if, deg a deg b (resp. deg a < deg b).

Now, if Π is an annihilator associated with ̺(s), then left multiplying both sides of (6) by a rational function σ(s) immediately shows that σ(s)Π is also an annihilator, associated with σ(s)̺(s). It is therefore easy to see that Lemma 1 An annihilator Π and its associated rational function ̺ as in (6) may always be chosen strictly proper.

Henceforth, we shall consider only strictly proper annihilators with strictly proper associated rational function. This means that only integral operators will be involved, once (6) is translated back in the time domain.

Time domain

We now replace the model x N by the actual noisy observation y. With (5), equation ( 6) then becomes

̺(s)x (n) N (0) = finite ̺ ℓ (s) d ℓ ds ℓ ŷ. ( 7
)
This equation provides us with an operational estimator of x (n) (0). It now remains to express the estimate x(n) N (0), directly in the time domain. For this, remember that ̺(s) and ̺ ℓ (s) are strictly proper rational (transfer) functions. Let f (t) and h ℓ (t) be their corresponding impulse responses and recall that, by the classical rules of operational calculus, d ds corresponds to multiplication by -t. Then equation [START_REF] Dabroom | Discrete-time implementation of high-gain observers for numerical differentiation[END_REF] readily reads in the time domain as

f (t)x (n) N (0) = finite t 0 h ℓ (t -α)(-1) ℓ α ℓ y(α)dα. (8)
Here, the variable t represents the estimation time. This equation has therefore to be considered for fixed t, say t = T . Now, for each estimation time T > 0 such that f (T ) = 0, we obtain a numerical estimate x(n) N (0). Examine this estimate. To do this, write R N (t) = x(t) -x N (t) from ( 2) and (3) and decompose the observation y in y

(t) = x N (t) + R N (t) + ̟(t). Rewriting (8) in the form x(n) N (0) = T 0 h(α; T ){x N (α) + R N (α) + ̟(α)}dα = x (n) (0) + e RN (0) + e ̟ (0), (9)
shows that the estimate is corrupted by two sources of error:

• the mismodelling, stemming from the truncation of the Taylor expansion, which contributes to e RN (0) (the bias term) and • the noisy environment which produces the error e ̟ (0) (the variance term).

Needless to say, the contribution of each one of these two sources may vary markedly from one estimator to another. So, the game of the next section will be to devise an adequate annihilator Π, along with its associated rational function ̺. By adequateness, we mean the minimisation of the overall estimation error. However for N fixed, minimizing the mismodelling effect is, most often, in opposition with reducing the noise error contribution. As an example, note that the estimate x(n) N (0) will vary with T in general. Now, for N fixed, choosing a smaller estimation time, T , will tend to reduce |e R N (0)|. But at the same time, a large T is required to filter out the noise showing that the choice for T should obey a compromise.

Pointwise derivative estimation

Now that we have fixed the notations and the guidelines of our methodology, we investigate in this section some detailed properties and performance of a class of pointwise derivative estimators. These will be derived from a particular family of annihilators. And, as we shall shortly see, the Jacobi orthogonal polynomials [START_REF] Szegö | Orthogonal polynomials, 3rd edn[END_REF] are inherently connected with these estimators. A least squares interpretation then naturally follows [START_REF] Mboup | Parameter estimation via differential algebra and operational calculus[END_REF], [START_REF] Mboup | A revised look at numerical differentiation with an application to nonlinear feedback control[END_REF] and this leads to one of the main contribution of the paper: the numerical differentiation is as efficient as an appropriately chosen time delay is tolerated.

Purely integral estimators

A linear differential operator Π is said to be in finite-integral form if, and only if, each ̺ ℓ (s) in ( 5) is in the form ̺ ℓ (s) = 1 s H ℓ ( 1 s ) for some polynomial H ℓ . We consider the family of annihilators defined below.

Lemma 2 For any κ, µ ∈ N, the differential operator

Π N,n κ,µ = 1 s N +µ+1 d n+κ ds n+κ • 1 s • d N -n ds N -n • s N +1 ( 10
)
is a finite-integral form annihilator for x (n) (0) and, it is associated with

ϕ κ,µ,N (s) = (-1) (n+κ) (n + κ)!(N -n)! s µ+κ+N +n+2 . (11)
This means that for each set of positive integers {κ, µ, N }, with N n 0, the equation ϕ κ,µ,N (s)x

(n)

N (0) = Π N,n κ,µ ŷ (12)
provides an operational estimator of the n th order derivative of x(t), observed by y(t). We thus obtain a family of estimators parametrized by κ, µ and N . In order to emphasis the dependence on these parameters, we introduce the notation x (n) (0; κ, µ; N ) for x

(n) N (0). If N = n, the differential operator Π N,n κ,µ reduces to Π n κ,µ = 1 s n+µ+1 d n+κ ds n+κ • s n . (13)
Accordingly, the simplified notations x(n) (0; κ, µ) and ϕ κ,µ (s) will be use instead of x(n) (0; κ, µ; n) and ϕ κ,µ,N (s) respectively. The next result will play a central rôle in the rest of the paper.

Theorem 1 Let N, n, κ and µ be four positive integers, with N n. Set q = N -n. Then the n th order derivative estimate x(n) (0; κ, µ; N ), obtained from a Taylor expansion of order N , is uniquely expressible in the form

x(n) (0; κ, µ; N ) = q ℓ=0 λ ℓ x(n) (0; κ ℓ , µ ℓ ), λ ℓ ∈ Q ( 14
)
where κ ℓ = κ+q+ℓ and µ ℓ = µ+ℓ. Moreover, if q n+κ, then the coefficients λ ℓ satisfy

q ℓ=0 λ ℓ = 1.
This theorem shows that an n th -order truncated Taylor expansion is appropriate for estimating the n th -order derivative. Also, it characterises x(n) (0; κ, µ; N ) as a point in the Q-affine hull of the set

S κ,µ,q = x(n) (0; κ + q, µ), • • • , x(n) (0; κ, µ + q) , (15)
provided q n + κ. Suppose now that, this affine hull is extended by allowing the barycentric coordinates λ ℓ in ( 14) to be real rather than rational. By doing so, it is clear that any point in this new larger set will represents an n th -order derivative estimate of x(t) at the origin, in some meaningful sense.

Characterising those points which minimise a given distance to x (n) (0) is an important question. This is the program of the following two subsections.

Minimizing the bias term

Jacobi orthogonal polynomials

We show in this subsection how the preceding special class of strictly proper differential operators are intimately connected to the Jacobi orthogonal polynomials. This connection, in turn, will allow us to attach a least squares interpretation to our estimators (see [START_REF] Mboup | Parameter estimation via differential algebra and operational calculus[END_REF] for more delails).

To begin let us recall that for a given signal û, and a positive integer α, the time domain analog of v = 1 s α d β ds β û is the iterated integral of order α, of (-1) β t β u(t). Using the well known Cauchy formula for repeated integration, this reduces to a single itegral

v(t) = 1 (α -1)! t 0 (t -τ ) α-1 (-1) β τ β u(τ )dτ.
By a direct application of this rule of operational calculus and by recalling ( 13), ( 12) and [START_REF] Fliess | Compression différentielle de transitoires bruités[END_REF], it becomes easy to see that, for any κ and µ,

x(n) (0; κ, µ) = (µ + κ + 2n + 1)! (µ + n)!(κ + n)! T 0 (T -τ ) µ+n τ κ+n T µ+κ+2n+1 y (n) (τ )dτ. ( 16
)
Note that even though the noise function may not be differentiable, the above expression still makes sense because the derivation under the integral sign is only formal; it disappears upon integrating by parts. Nonetheless, to avoid any possible misunderstanding, we rewite ( 16) by considering separately the contribution due to the signal x and that due to the noise ̟. For the former, we keep the form of ( 16) and for the latter, we integrate by parts, to obtain

x(n) (0; κ, µ) = γ κ,µ,n 1 0 τ κ+n (1 -τ ) µ+n x (n) (T τ )dτ +(-1) n γ κ,µ,n T n 1 0 d n dτ n τ κ+n (1 -τ ) µ+n ̟(T τ )dτ (17)
where γ κ,µ,n △ = (µ+κ+2n+1)! (µ+n)!(κ+n)! . Here we have also used a change of variable to reduce the estimation interval

I T 0 = [0, T ] to [0, 1]. Consider now P {κ,µ} i (t) i 0
, the Jacobi orthogonal polynomials (cf. [START_REF] Szegö | Orthogonal polynomials, 3rd edn[END_REF][START_REF]Handbook of mathematical functions[END_REF]) on the interval [0, 1], associated to the weight function

w κ,µ (t) = t κ+n (1 -t) µ+n . (18)
In all the sequel, • will denote the norm induced by the inner product defined, for two real valued functions f and g, by:

f, g {κ,µ} = 1 0 f (t)w κ,µ (t)g(t)dt. (19)
The short hand notation •, • will be used wherever there is no possible confusion. Upon noting that P {κ,µ} 0 (t) ≡ 1 and

P {κ,µ} 0 2 = (µ + n)!(κ + n)! (µ + κ + 2n + 1)! = 1 γ κ,µ,n , (20)
we have the following:

Proposition 1 Let x (n)
LS,1 (t) denote the first order least-squares polynomial approximation of

x (n) (t) in the interval [0, T ]. Then x(n) (0; κ, µ) is given by x(n) (0; κ, µ) = x (n) LS,1 (T ξ 1 ) + e ̟ (0; κ, µ), (21)
where

ξ 1 = κ + n + 1 µ + κ + 2(n + 1) (22)
is the root of P {κ,µ} 1

(t) and e ̟ (0; κ, µ) is the noise contribution as given by the second term in the right hand side of (17).

Time-delayed derivative estimation

For a given τ 0, let us substitute, in all the preceding developments, y(t) by [Heaviside(τ )y(τ + t)]. This simply amounts to moving the time origin from 0 to τ . As a result, we obtain x(n) (τ ; κ, µ), the estimate of x (n) (τ ). Observe however that, the corresponding estimation is anti-causal, since the estimate at time τ is based on the signal observation during the time window

I T τ+ = [τ, τ + T ].
Causal estimation is readily achieved by substituting y(t) by -[Heaviside(t)y(τ -t)], for τ T . Doing so, we still obtain in x(n) (τ ; κ, µ), the estimate of x (n) (τ ) but now from the signal observation on I T τ-= [τ -T, τ ]. These sliding windows thus induce a map t → x(n) (t; κ, µ), as an image of the map t → y(t) through all the preceding developments. The underlying operator will subsequently be called the (κ, µ)-algebraic numerical differentiator (AND) of order n. A (κ, µ)-AND will be termed minimal if it has been designed from the minimum truncation order of the Taylor expansion (3). It may be either causal or anti-causal, according to which estimation interval I T τ-or I T τ+ is used. Now, it is straightforward to infer from Proposition 1 that, the noise contribution e ̟ (t; κ, µ) apart,

x(n) (t; κ, µ) = x (n)
LS,1 (t -T ξ 1 ), t T, for the causal variant. As a consequence, we have:

Corollary 1 Let x(t) be observed through a noise corruption ̟(t) by y(t) = x(t) + ̟(t). The causal n th order minimal (κ, µ)-algebraic numerical differen- tiation of y(t) is time-delayed: x(n) (t; κ, µ) ≈ x (n) (t -τ 1 ), t T. (23)
The delay τ 1 reads as τ 1 = T ξ 1 where ξ 1 is given by [START_REF] Levant | Higher-order sliding modes, differentiation and output-feedback control[END_REF].

Similarly, we obviously have for the anti-causal variant,

x(n) (t; κ, µ) ≈ x (n) (t + τ 1 ), t 0. ( 24
)
Clearly, the above approximations become exact for polynomial signals of degree not exceeding n, in the noise free context. But since ξ 1 , in Proposition 1, is the root of P {κ,µ} 1

(t), we deduce:

Corollary 2
The approximations in [START_REF] Lorentz | Bernstein Polynomials[END_REF] and (24) are exact for polynomial signals of degree up to n + 1.

Delay and model complexity

In this paragraph, we investigate the non minimality effects on the algebraic numerical differentiators. We shall first establish that non minimality allows to compensate the previous delay. This is important because a delay may not be tolerable in some real time applications1 . We shall also show that for a fixed modelling order N , the abscence of delay undergoes a performance loss. Finally, the question of the "best" choice of the barycentric coordinates in ( 14) will be approached.

To begin, let us recall some well-known facts. Consider the subspace of L 2 ([0, 1]), defined by

H q = span P {κ,µ} 0 (t), P {κ,µ} 1 (t), • • • , P {κ,µ} q (t) . ( 25
)
Equipped with the inner product ( 19), H q is clearly a reproducing kernel Hilbert space [START_REF] Aronszajn | Theory of reproducing kernels[END_REF], [START_REF] Alpay | Algorithme de Schur, espaces à noyau reproduisant et théorie des systèmes[END_REF], [START_REF] Saitoh | Theory of reproducing kernels and its applications[END_REF], with the reproducing kernel

K q (τ, t) = q i=0 P {κ,µ} i (τ )P {κ,µ} i (t) P {κ,µ} i 2 . ( 26
)
The reproducing property implies that for any function f defined on [0, 1], we have

K q (τ, •), f (•) = f q (τ ),
where f q stands for the orthogonal projection of f on H q . In the rest of the paper, we will set q = N -n and assume that q 0. Theorem 2 Assume that q κ + n and let

x (n) LS,q (t) △ = q i=0 P {κ,µ} i (τ ), x (n) (T τ ) P {κ,µ} i 2 P {κ,µ} i (t/T ), (27)
denote the least-squares q th order polynomial approximation of

x (n) (•) in the interval I T 0-= [0, T ]. Then x(n) (0; κ, µ; N ) is given by x(n) (0; κ, µ; N ) = x (n) LS,q (0) + e ̟ (0; κ, µ; N ), (28)
where e ̟ (0; κ, µ; N ) represents the corresponding noise contribution.

Remark 1 This is for the anti-causal variant. To obtain the causal counterpart, just replace y(t) for t ∈ [0, T ] by -y(T -t).

We immediately deduce from the above theorem that the corresponding causal (κ, µ)-AND is delay-free:

x(n) (t; κ, µ; N ) ≈ x (n) (t), t T. (29)
Now, since x

LS,q (T t) is the orthogonal projection of x (n) (T t), t ∈ [0, 1] on H q , we may rewrite (28) as

x(n) (0; κ, µ; N ) = K q (0, t), x (n) (T t) + e ̟ (0; κ, µ; N ) = q ℓ=0 λ ℓ x(n) (0; κ + q -ℓ, µ + ℓ), (30)
where we recall (14) (for q κ + n). As already suggested, we henceforth allow the barycentric coordinates λ ℓ , ℓ = 0, . . . , q to belong to R instead of being restricted in Q. We obtain the next result.

Proposition 2 Let κ, µ, n and q be four positive integers as above. For any ξ ∈ [0, 1], there exists a unique set of real coordinates

λ ℓ (ξ) ∈ R, ℓ = 0, • • • , q such that q ℓ=0 λ ℓ (ξ) x(n) (0; κ ℓ , µ ℓ ) = K q (ξ, t), x (n) (T t) + e ξ ̟ (0; κ, µ; n + q), (31)
for some noise contribution e ξ ̟ (0; κ, µ; n + q). Moreover, these coordinates satisfy

q ℓ=0 λ ℓ (ξ) = 1
and, unless all the λ ℓ (•) are equal to zero except one, the following holds:

min ℓ λ ℓ (•) < 0.
For each ξ ∈ [0, 1], this proposition establishes a correspondance between x (n) LS,q (T ξ) and a point in the R-affine hull of S κ,µ,q in (15). Let

x(n) ξ (0; κ, µ; N ) denotes such a point as given in [START_REF] Roberts | Digital signal processing[END_REF] and consider as before, the map t →

x(n) ξ (t; κ, µ; N ) :

x(n) ξ (0; κ, µ; N ) = K q (ξ, t), x (n) (T t) + e ξ ̟ (0; κ, µ; n + q). (32)
Then, following the previous paragraph we clearly see that the corresponding causal algebraic numerical differentiator is time delayed:

x(n) ξ (t; κ, µ; N ) ≈ x (n) (t -τ ), t T, (33)
where the delay is given by τ = T ξ.

In order to simplify some notations to come, we turn to the anti-causal variant (replace

x (n) (t -τ ), t T above by x (n) (t + τ ), t 0). Let now 0 < ξ {q+1} 1 < ξ {q+1} 2 < • • • < ξ {q+1} q+1
< 1 represent the (q + 1) zeros of P {κ,µ} q+1 . Then, the noise aside, one may readily check that for each i = 1, . . . , q + 1,

x(n) ξ {q+1} i (t; κ, µ; n + q) = K q+1 (ξ {q+1} i , u), x (n) (t + T u) = x (n) LS,q+1 (T ξ {q+1} i ). ( 34
)
In other words, in the noise-free case, the approximation (33) with ξ = ξ {q+1} i , 1 i q + 1 would be exact for all polynomial signals up to degree N + 1 although this approximation is based on a Taylor expansion of order N . This has an important implication as shown next.

Delay and precision

To continue the discussion, we establish in Proposition 3 below that the quality of the estimation may be significantly improved by admitting a delay. When relating this result with the last remark of the preceding subsection, we somehow recover the paradigm of model section [START_REF] Massart | Concentration inequalities and model selection[END_REF]: a simple and approximate model may outperform a precise and more complex one.

To begin, let us expend x (n) (t + T u), for each fixed t, as

x (n) (t + T u) = a i (t)P {κ,µ} i (u), u ∈ [0, 1], (35)
and set ξ {q} i , i = 1, . . . , q for the zeros of P {κ,µ} q

, in increasing order. Suppose that q = N -n 1. For some given i, consider first x(n)

ξ {q} i (t; κ, µ; N -1), the time-delayed τ = T ξ {q} i
estimate [START_REF] Su | A simple improved velocity estimation for low-speed regions based on position measurements only[END_REF], based on an (N -1)-order truncated Taylor expansion. Using [START_REF] Szegö | Orthogonal polynomials, 3rd edn[END_REF] and [START_REF] Su | A simple improved velocity estimation for low-speed regions based on position measurements only[END_REF], the estimation error due to the truncation reads, at each time t, as:

e ξ {q} i RN-1 (t; κ, µ; n + q -1) = a q+1 (t)P {κ,µ} q+1 (ξ {q} i ) + i 2 a q+i (t)P {κ,µ} q+i (ξ {q} i ). (36)
Next, suppose that 1) we do not want the estimation to be time-delayed and 2) we have a more precise model, namely a Taylor truncation of order N instead of N -1. Accordingly, the estimation error e RN (t; κ, µ; n + q) corresponding to the delay-free estimate x(n) (t; κ, µ; N ) in ( 28)-( 30) and stemming from the mismodelling is given by: e RN (t; κ, µ; n + q) = a q+1 (t)P {κ,µ}

q+1 (0) + i 2 a q+i (t)P {κ,µ} q+i (0). (37)
Let us now compare the magnitude of these two errors. For that, let us notice that for T small enough, it is reasonable to approximate the above error by the first term of the expansion (37):

e RN-1 (t; κ, µ; n + q) ≈ a q+1 (t)P {κ,µ}

q+1 (0) (38)
and likewise for e ξ {i} q RN-1 (t; κ, µ; n+q-1) in [START_REF] Yosida | Operational Calculus -A Theory of Hyperfunctions[END_REF]. Using the three terms recurrence relation of the Jacobi polynomials and the recurrence relation (61) (see proof of Theorem 2), we obtain the following: Proposition 3 Denote by z q (κ, µ) and Z q (κ, µ) the lower and upper bounds of the smallest and largest zeros of P {κ,µ} q , respectively. Then

e ξ {q} 1 RN-1 (t; κ, µ; n + q -1) |e RN (t; κ, µ; n + q)| < θ q (κ, µ) 1 - z q (κ, µ) Z q-1 (κ, µ) q-1 , (39)
where

θ q (κ, µ) = q(q + n + µ)(2(q + n) + µ + κ + 2) (q + n + κ + 1)(q + 2n + µ + κ + 1)(2q + 2n + µ + κ) < 1. ( 40
)
Explicit expressions for z q (κ, µ) and Z q (κ, µ) may be obtained from [21, Theorem 2]. Anyway, for i = 1 in [START_REF] Yosida | Operational Calculus -A Theory of Hyperfunctions[END_REF], which amounts to choosing the minimum delay2 τ = T ξ {q} 1 , a significant improvement is achieved over a delay-free setting, even with a more precise model.

Remark 2 For the causal variant, the minimum delay will be given by τ = T ξ {q} q and all the preceding arguments remain valid upon interchanging κ with µ, and τ with 1 -τ .

Finally, we have seen that a causal non minimal (κ, µ)-AND may provide either delay-free or time-delayed estimation. The performance obtained in the delay-free case are worst than for the time-delayed counterpart, even though a more precise signal model is used. Admitting a delay significantly improves the quality of an (κ, µ)-AND. Thus, we recover the paradigm of model selection [START_REF] Massart | Concentration inequalities and model selection[END_REF]. On the other hand, it appears that the question of a "good" choice for the barycentric coordinates reduces to the choice of the time-delay, provided the noise effect is ignored.

Reducing the noise effect

In the preceding analysis, we have confined our attention to minimizing the mismodelling error e RN . Here, we investigate the effect of the noise. The analysis relies on the very classical properties of linear time-invariant filtering. Basically, it amounts to adjust the free parameters of a filter in order to maximize the output signal-to-noise ratio (SNR). Clearly, this approach is in contrast with the previous local analysis.

Implementation issues

As a matter of fact, a causal (κ, µ)-AND is equivalent to a stable and causal linear time-invariant filter, with finite impulse response. To see this, rewrite the equation [START_REF] Haykin | Signals and Systems[END_REF] for the causal variant,

x(n) (t; κ, µ) = (-1) n γ κ,µ,n T n 1 0 d n dτ n τ κ+n (1 -τ ) µ+n y(t -T τ )dτ , ( 41
)
where t stands for the present time. Setting y T (τ ) = y(T τ ) = x(T τ ) + ̟(T τ ) and

h κ,µ (τ ) △ = (-1) n γ κ,µ,n T n rect(τ ) d n dτ n w {κ,µ} (τ ), (42)
where rect(•) is the characteristic function of the interval [0, 1], one may immediately rewrite (41) as the convolution product:

x(n) (T t; κ, µ) = ∞ 0 h κ,µ (τ )y T (t -τ )dτ. ( 43
)

Minimizing the output noise variance

Given κ, µ and n, we consider that the estimation time T is chosen so that the mismodelling error is small. We are now interested in describing the points in the R-affine hull of the set S κ,µ,q (t) = {x (n) (t; κ + q -ℓ, µ -ℓ), ℓ = 0, . . . , q} which minimise the noise effect. To proceed, we suppose that the noise is a wide-sense stationary random process. Consider

x(n) (T t; κ, µ; N ) = ∞ 0 q ℓ=0 λ ℓ h κ+q-ℓ,µ+ℓ (τ )y T (t -τ )dτ, (44)
a point, with barycentric coordinates λ ℓ , ℓ = 0, . . . , q, in this affine hull. Set λ = [λ 0 • • • λ q ] t where the superscript t denotes transposition. Then one may readily check that the associated output noise variance is

var(e ̟ ) = q ℓ,m λ ℓ λ m 1 0 1 0 h κ ℓ ,µ ℓ (τ ) C(T (t -τ )) h κm,µm (τ )dτ dt, = λ t Rλ, (45)
where

κ i = κ + q -i, µ i = µ + i, C(τ ) △ = E[̟(t)̟(t -τ )
] stands for the noise convariance function and R is the (q + 1)-square matrix with entries

R i,j = 1 0 1 0 h κ+q-i,µ+i (τ )C(T (t -τ ))h κ+q-j,µ+j (τ )dτ dt.
Define by u the (q + 1)-dimensional real vector, with 1 in each entry. Now, the barycentric coordinates λ ̟ minimizing this variance are immediately given by

λ ̟ = R -1 u u t R -1 u . ( 46
)
By proposition 2, we know that x(n) (t; κ, µ; N ) associated to λ ̟ is in the form x(n) (t; κ, µ; N ) = K q (ξ, τ ), (-1) n x (n) (t -T τ ) + e ξ ̟ (t; κ, µ; N ), for some ξ ∈ [0, 1] if, and only if, at least one entry of λ ̟ is strictly negative. However, the barycentric coordinates which minimise the noise variance (45) do not fulfill this requirement, in general. If ̟(t) is a (bandlimited) white noise for instance, then we directly verify that (46) yields

λ ̟ = 1 q + 1 u, (47)
showing that the minimum output mean square error is achieved by the centroid of the points x(n) (t; κ + q -ℓ, µ -ℓ), ℓ = 0, . . . , q.

If, more generally, we consider the (κ, µ)-AND associated now to the points in the convex hull of S κ,µ,q (t), then we obtain a class of differentiators different from the one described so far. The difference is, as stated by proposition 2, that these (κ, µ)-AND do not have the orthogonal projection principle attached to them.

We postpone the detailed analysis of this class to a forthcoming paper and go directly to the simulation studies.

Numerical experiments

A digital FIR filter

Before all, let us remind that the numerical estimates are obtained from (44) which we rewrite as

x(n) (T t; κ, µ; N ) = 1 0 g(τ )y T (t -τ )dτ, ( 48
)
with the simplified notation g(t) = q ℓ=0 λ ℓ h κ+q-ℓ,µ+ℓ (τ ). This is a very classical input-output linear time-invariant filtering relation, the implementation of which is easy. To further simplify the implementation, we consider that only discrete-time samples of the observation y are available. We assume that the samples are regularly spaced3 , with a sampling period T s and we use the common notation y m = y(mT s ). For an estimation time T = M T s , we set g m = g(mT s /T ), m = 0, . . . , M . These numbers g m are computed once for all, after all the parameters λ ℓ , κ, µ, n, q have been selected.

Let now W m and t m = mT s /T, m = 0, . . . , M be the weights and abscissas associated with a given numerical integration method as in

1 0 f (t)dt ≈ M m=0 W m f (t m ).
For each sampling time t ℓ = ℓT s , we compute

x(n) (t ℓ ; κ, µ; N ) ≈ x(n) ℓ = M m=0 W m g m y ℓ-m . (49)
The n th order derivative estimate is thus obtained from the output of a classical FIR digital filter with impulse response c m = W m g m , m = 0, . . . M , where the input is the noisy observation samples signal. In all the subsequent simulation, we use the work horse trapezoidal method, corresponding to

W 0 = W M = T s 2
, and W m = T s , m = 1, . . . , M -1.

Simulation results

Let y(t) = x(t) + ̟(t), 0 ≤ t ≤ 5 (without unit), be the noisy measurement of the signal x(t) = tanh(t -1) + e -t/1.2 sin(6t + π).

We consider a sampling period (given without unit) of T s = 1/200 which corresponds to 1000 samples in the interval [0, 5]. The samples of the noise are simulated from a zero-mean white Gaussian iid sequence. The variance is adjusted in such a way that the signal to noise ratio in dB 4 , i.e., SN R = 10 log 10 P |yi| 2 P |̟i| 2 , corresponds to SN R = 25 dB (see figure 1). Begin with the first order derivative. The estimates x(1) (t; 0, 0), obtained with the minimal (0, 0)-AND are displayed in figure 2 (dotted line). Each estimation window I T t-= [t -T, t] has a width T = 60T s . We also depict (dashed line) the exact derivative of the noise-free signal in order to gauge the estimation accuracy. Compare with the solid line curve, obtained by a left shift of the dotted line curve with an amount corresponding to the delay ξ 1 T in equation [START_REF] Levant | Higher-order sliding modes, differentiation and output-feedback control[END_REF]. With N = 1, the linear signal model results in poor estimates on the intervals where the signal's dynamic is strong. For high signal-to-noise ratio, the estimates may be improved by reducing the estimation time T . 0.0 0.5 Alternatively, one may consider a richer signal model e.g. with N = 2. This is the case of the next experiment. We consider the non minimal estimator from ( 10)- [START_REF] Fliess | Analyse et représentation de signaux transitoires : application à la compression, au débruitage et à la détection de ruptures[END_REF], with sliding windows of width T = 110T s . The estimates x(1) (t; 0, 0; 2) are plotted (solid line) in figure 3 below. There is no estimation delay, as expected. However, the performance significantly degrades as compared to the preceding results although the signal model is more precise. If we now relax the delay-free constraint, it becomes possible to take advantage of the more flexible second order model for the signal. This is illustrated in the following simulation (see Figure 4), where we keep the same settings for T , but with κ = µ = 1 now. The solid line curve in Figure 4 represents the estimates obtained as in Proposition 2, with the barycentric coordinates given by equation ( 68) where ξ = ξ {2} 2 is the largest root of P {κ,µ} 2

. For comparison 4 See, e.g., [START_REF] Haykin | Signals and Systems[END_REF] for this well known concept in signal processing

Time

Fig. 3 Non minimal, delay-free derivative estimation.

mation of the higher order derivatives. As an illustration, we give in figure 5 and figure 6 below the estimation results corresponding to the second order derivative from the noisy observation depicted in figure 1. We represent the transient components of the corresponding filter's output by the dotted parts of the curves. The dash-dotted line curves are the delayed estimates. The predicted values of the respective delays match the experiment as can be seen by comparing the exact second order derivative (dashed line curves) with the shifted versions of the estimates (solid line curves).

Observe that for q = 1 and for all n, the filter corresponding to the minimal (κ, µ)-AND coincides with that which minimises the output noise variance for the same κ and µ. More generally, we may verify that any minimal (κ, µ)-AND, with small |κ-µ|, is close to the corresponding minimum variance (κ, µ)-AND, for q = 1. This stems from the identity

h M V κ,µ (τ ) = κ + µ + 2n + 2 2(µ + n + 1) h M in κ,µ (τ ) + mu -κ 2(κ + µ + 2n + 2) h M in κ+1,µ (τ )
where the superscripts "MV" and "Min" refer to the "Minimum Variance" and "Minimal" estimators respectively. An illustration is given in the figure 7 below, with the settings of the two preceding simulation. Finally, let us notice that all the preceding simulation results are based on one sample realization of the noisy signal. Now, in order to validate the presented results, we investigate the bias and variance in each simulation above. Recall that the signal derivatives estimates are obtained as output of FIR filters with a noisy input. So, with the zero mean white noise assumption, the associated output noise variance is readily computed by taking the ℓ 2 sequence norm of the corresponding impulse response, {c m }, multiplied by the variance of the input noise. The output signal to noise ratio (SNRout) is therefore given by

SNRout = 1 L-M L ℓ=M x (n) ℓ 2 σ 2 ̟ M -1 m=0 |c m | 2 .
Concerning the bias, we use 100 Monte-Carlo simulation runs to compute the sample means E s {x (n) (ℓ)}. The relative sample mean squared error, due to the mismodelling (the bias term) and to the numerical approximation of the integrals, will be called the signal to bias ratio. It is computed as:

SBR = L ℓ=M x (n) ℓ 2 L ℓ=M x (n) ℓ -E s {x (n) (ℓ)} 2 .
Table 4.2 displays the signal to bias ratio (SBR) and the output signal to noise ratio (SNRout). The values are expressed in decibel (dB). Here, L = 1000 is the signal duration, that is: LT s = 5. Each column, except the last one, represents the results obtained in the context of the simulation of the figure referenced in the first row of the table. The last column displays the results obtained with the second order minimum variance derivative estimator filter show in figure 7.

Proofs

Proof (Proof of Lemma 2) A direct inspection of [START_REF] Fliess | Closed-loop fault-tolerant control for uncertain nonlinear systems[END_REF] shows that the maximun power of s, appearing in the coefficients of the differential operator Π N,n κ,µ , is s -µ-1 . That Π N,n κ,µ is in finite-integral form is thus clear since µ 0.

Next apply this differential operator to xN . The (N -n) th order derivative of the right hand side of (4) results in a polynomial of degree n, the constant term of which is equal to (N -n)!x (n) (0). If we divide this polynomial by s and take the (n + κ) th order derivative of the result, then we immediately see that there is only one remaining term which, when divided by s N +µ+1 , yields ϕκ,µ(s)x (n) (0), where ϕκ,µ(s) is given by [START_REF] Fliess | Compression différentielle de transitoires bruités[END_REF]. ds q s q+1 s n in the form

D = P q i=0 `q i
´(q+1)! (q+1-i)! d p ds p s q-i d q-i ds q-i s n = P q i=0 P min (p+i,q) j=i a i,j s q-j d p+q-j ds p+q-j s n ,

where a i,j = " q i "" p j -i " (q + 1)! (q + 1 -i)(q -j)! .

Dividing D by s µ+N +1 and rearranging the terms in the summation above allows one to get:

Π N,n κ,µ = q X j=0 0 @ j X i=m(j) a i,j 1 A Π n κ+q-j,µ+j ,
where, if q p then m(j) ≡ 0 and else m(j) = max(0, j -p). Using the expression [START_REF] Fliess | Compression différentielle de transitoires bruités[END_REF] and the relation [START_REF] Fliess | Analyse et représentation de signaux transitoires : application à la compression, au débruitage et à la détection de ruptures[END_REF], we get the linear combination in [START_REF] Fliess | Questioning some paradigms of signal processing via concrete examples[END_REF] with coefficients given by:

λ j = (-1) q-j (n + κ + q -j)! p!q! j X i=m(j)
a i,j .

It now remains to show that, ( 12) is a Q-affine combination. That the coefficients λ j are rational numbers is indeed plain. By rearranging the factors of the numerator and of denominator, we obtain the simpler expression:

λ j = (-1) q-j `p+q-j p ´Pj i=0 `p j-i ´`q+1 i = (-1) q-j `p+q-j p ´`p+q+1 j ´ (51)
for q p. Define I(q) = P q j=0 λ j . Then, with the help of the identity `m+1 ℓ ´= `m ℓ ´+ `m ℓ-1 ´, ℓ > 0, we have the recurrence relation

I(q + 1) = I(q) + P q+1 j=0 (-1) j `p+j j ´`p+q+1 p+j = I(q) + `p+q+1 p ´Pq+1 j=0 (-1) j `q+1 j ´= I(q),
showing that I(q) is constant. Therefore, we get P q j=0 λ j = 1 upon noting that I(1) = 1. ⊓ ⊔ Proof (Proof of Proposition 1) Observe that x

(n) LS,1 (T t) is defined by x (n) LS,1 (T t) = P {κ,µ} 0 (τ ), x (n) (T τ ) P {κ,µ} 0 2 + P {κ,µ} 1 (τ ), x (n) (T τ ) P {κ,µ} 1 2
P {κ,µ} 1 (t).

For t = ξ 1 as in [START_REF] Levant | Higher-order sliding modes, differentiation and output-feedback control[END_REF], i.e., the zero of P {κ,µ} 1 (t), the second term in the right hand side of the equality above vanishes. Hence, we recognize the contribution due to the signal in equation [START_REF] Haykin | Signals and Systems[END_REF]. The contribution due to the noise is denoted by the term e̟(0; κ, µ) in [START_REF] Ismail | Bound on the extreme zeros of orthogonal polynomials[END_REF].

⊓ ⊔ Proof (Proof of Theorem 2) To simplify the expressions to come, we shall write x(n) (κ, µ; N ) in place of x(n) (0; κ, µ; N ). Begin with the following lemma.

Lemma 3 The recurrence relation

x(n) (κ, µ; N ) = λx (n) (κ, µ + 1; N -1) + (1 -λ)x (n) (κ + 1, µ; N -1), ( 52
)
where λ is defined by λ = N +κ+1 N -n , holds for any κ, µ, n, N positive, with N > n.

The following miscellaneous recurrence relations, which are adapted from [1, p.782], will be useful:

(

P {κ,µ} m-1 (0) = (m + α + 1)(m + α + β + 1) (2m + α + β + 1)(2m + α + β + 2) (0)
With these relations, the equation (59) becomes

x(n) (κ, µ; n + q + 1) = q+1 X m=0 ζm P {κ,µ} m , x (n) , (63)
with

ζ 0 = λa 0 (κ, µ)η 0 (κ, µ + 1) + (1 -λ)a 0 (µ, κ)η 0 (κ + 1, µ), ζ q+1 = -λη q (κ, µ + 1) + (1 -λ)η q (κ + 1, µ), ζ m = λ[am(κ, µ)ηm(κ, µ + 1) -η m-1 (κ, µ + 1)] +(1 -λ)[a m (µ, κ)η m (κ + 1, µ) -η m-1 (κ + 1, µ)] 1 m q,
Now, it is straightforward to verify that the following identities hold:

ζ 0 = 1 P {κ,µ} 0 2
,

ζ q+1 = η q+1 (κ, µ), η m(κ, µ) = am(κ, µ)ηm(κ, µ + 1) -η m-1 (κ, µ + 1) = a m(µ, κ)ηm(κ + 1, µ) -η m-1 (κ + 1, µ) = ζm, m = 1, . . . , q
This shows that if the relation ( 28) is valid for q > 1, then it must also be valid for q + 1.

Since it has been established for q = 1, we conclude that (28) holds for all q 1. ⊓ ⊔ Proof (Proof of Proposition 2) As for the equation (58) in the previous proof, it is easy to see that the linear combination P q ℓ=0 λ ℓ x(n) (0; κ ℓ , µ ℓ ) is given by

q X ℓ=0 λ ℓ x(n) (0; κ ℓ , µ ℓ ) = Z 1 0 B(τ )w {κ,µ} (τ )y (n) (T τ )dτ = B(τ ), y (n) (T τ ) . (64)
with

B(τ ) = q X ℓ=0 λ ℓ γκ ℓ ,µ ℓ ,n τ q-ℓ (1 -τ ) ℓ = q X ℓ=0 λ ℓ φ ℓ b ℓ,q (τ ), (65)
where b ℓ,q (τ )

△ = " q ℓ " τ q-ℓ (1 -τ ) ℓ , ℓ = 0, . . . , q,
are the Bernstein basis polynomials [START_REF] Lorentz | Bernstein Polynomials[END_REF] of degree q and φ ℓ

△ = γκ ℓ ,µ ℓ ,n " q ℓ " .
Since the set ˘bℓ,q (τ ), ℓ = 0, . . . , q ¯also forms a basis for the subspace Hq (see ( 25)), we may equate

Kq(ξ, τ ) = q X ℓ=0 q X m=0 bm,q(ξ) (B -1 ) m,ℓ b ℓ,q (τ ). (66)
Here, the matrix B stands for the Gramian of the Bernstein polynomials b ℓ,q (•), with respect to the inner product of Hq. Now, observe that the first part of the proposition claims that, for each ξ ∈ [0, 1], the equation b 0,q (τ ) . . . b q,q (τ )

B(τ ) = K q (ξ, τ ) (67
It now remains to establish the second part of the proposition. To proceed, observe that the entries of the Gramian matrix B are readily identified as:

(B) i,j = " q i " ' ' 'P

{κ+2q-(i+j),µ+(i+j)} 0 ' ' '

2 " q j " .

Next, using [START_REF] Ibrir | A numerical procedure for filtering and efficient high-order signal differentiation[END_REF] and some straightforward algebra allows one to verify that (B) i,j φ j = `µ+n+i+j µ+n+j

´`κ+n+2q-(i+j) κ+n+q-j μ+κ+2n+2q+1 q

´. Now, with the help of the classical identity

X m " m ζ "" i -m j -ζ " = " i + 1 j + 1 " ,
which holds independently of ζ, one may readily check that q X j=0 (B) i,j φ j = 1, ∀i.

We have therefore shown that u t = [1 1 . . . 1] is a left eigenvector of the matrix BΦ, and the associated eigenvalue is 1. Last, since the Bernstein polynomials b ℓ,q (τ ), ℓ = 0, . . . , q form a partition of unity [START_REF] Lorentz | Bernstein Polynomials[END_REF], we get: q X ℓ=0 λ ℓ (ξ) = u t λ(ξ) = u t {BΦ} -1 bq(ξ) = u t bq(ξ) ≡ 1.

Finally, note that since the Bernstein basis polynomials are strictly positive in (0, 1), the definition of B(•) in (65) and the equation (67) imply that the condition λ ℓ (ξ) 0 for all ℓ is equivalent to We now proceed to bound the right hand side of the above expression. For that, recall first the three terms recurrence relations generating the family of the Jacobi polynomials in [0, 1], associated with the weight defined in (18): (0) = aq(µ, κ)a q-1 (µ, κ)P {κ,µ} q-1 (0). Using all these relations and owing that the polynomials P = dq • 1 aq(µ, κ)a q-1 (µ, κ) we end up with the expression (40) by a direct computation. Also, it is readily checked that θq(κ, µ) < 1. On the other hand, notice that the interlacing property 0 < ξ {q} 1 < ξ {q-1} 1 < ξ {q-1} j , j = 2, . . . , q -1 implies that each term 1 -ξ {q} 1 ξ {q-1} j of the product in (72) is positive and strictly less than 1. This product is therefore much lower than one and particular, we have

q-1 Y j=1 0 @ 1 - ξ {q} 1 ξ {q-1} j 1 A < "
1 -zq(κ, µ) Z q-1 (κ, µ) « q-1 , zq(κ, µ) and Zq(κ, µ) the lower and upper bounds of the smallest and largest zeros of P {κ,µ} q , respectively.

⊓ ⊔

TimeFig. 1

 1 Fig. 1 Noisy observation signal, SN R = 25dB.

TimeFig. 2

 2 Fig.2Estimation of the signal derivative: minimal (0, 0)-AND.

TimeFig. 4 Fig. 5

 45 Fig. 4 Non minimal, delayed derivative estimation.

Fig. 6

 6 Fig.6Non minimal, delayed second order derivative estimation.

1 Fig. 7

 17 Fig. 7 Minimum variance vs Minimal second order derivative filters.

  of Theorem1) Set q = N -n and p = n + κ. The first part follows by direct inspection, upon writing D = d p

  ) admits a unique solution for (λ 0 , . . . , λq). With the notations

  and with (65) and (66), this equation becomes equivalent to Φλ = B -1 bq(ξ). It is clear that it admits the unique solution:λ(ξ) = Φ -1 B -1 bq(ξ).

1 R N - 1 ( 1 R N - 1 (

 1111 τ ∈ (0, 1). Now, it is clear that this cannot hold. The above equality is the classical Christoffel-Darboux formula.⊓ ⊔Proof (Proof of Proposition 3) From the first order approximations of ˛eR N (t; κ, µ; n + q) ąs in (38) and of ˛eξ {q} t; κ, µ; n + q -1) ˛, we may write ˛eξ {q} t; κ, µ; n + q -1) ęR N (t; κ, µ; n + q)

P

2 =

 2 m(m + β)(m + α + β)(m + α) (2m + α + β -1)(2m + α + β) 2 (2m + α + β + 1),with α = µ + n and β = κ + n. For t = ξ{q} 1 , we get P On the other hand, setting t = 0 in (61) provides us with the recurrence P {κ,µ} q+1 (0) = -aq(µ, κ)P {κ,µ} q

Table 1

 1 7 Bias and variance

	SBR in dB	23.2	12	33.7	18.1	28.3	17.9
	SNRout in dB	26.5	12	26.6	22.6	24.1	22.9

  1 -t)P

		{κ,µ+1} m	(t) = a m (κ, µ)P m {κ,µ}	(t) -P m+1 (t) {κ,µ}	(60)
		tP m {κ+1,µ}	(t) = am(µ, κ)P m {κ,µ}	(t) + P m+1 (t) {κ,µ}	(61)
	where	a m(κ, µ) = -	P m {κ,µ}

For instance, it is known that introducing delayed signals in a control loop of a system tends to destabilize it[START_REF] Richard | Time delay systems: An overview of some recent advances and open problems[END_REF],[START_REF] Seuret | Robust exponential stabilization for systems with time-varying delays[END_REF].

We still use the term delay although one should rather talk about advance as long as the anti-causal variant is in concern.

Observe that non uniform sampling can also be handled in exactly the same way.

Acknowledgements The authors are gratful to the referees for the comments and suggestions which improve the clarity of this paper.

Proof For the sake of simplicity, we ignore the noise. This is without any loss of generality in the present proof. Set q = N -n, ν = µ + N + 1 and

Also, define x T (t) △ = x(T t). A direct application of the rules of operational calculus shows that:

Integrate by parts for the second integral. Then, we may rewrite the above expression as:

Let B N -1 (κ, µ) be defined by the second term of the right hand side of the previous equation:

). Now reconsider (54), and integrate by parts again but this time, with respect to the first integral. This will result in:

Combining ( 55) and (56), we arrive at the recurrence relation:

Finally, the relation (52) follows by considering the definition (53).

⊓ ⊔

We proceed by induction. For q = 1, equation (52) coincides with ( 14), with λ 1 = n + κ + 2 = λ and λ 0 = 1 -λ. Using [START_REF] Haykin | Signals and Systems[END_REF], it can be written as

Set α = µ + n and β = κ + n. Now, it is immediate to verify that the bracketed term under the integral above has the form g(τ )w {κ,µ} (τ ) where g is a polynomial of degree 1. The following identity η m(κ, µ)

may be verified directly from [START_REF]Handbook of mathematical functions[END_REF]. Then it is a matter of straightforward algebra to check that the polynomial g expresses as:

where we recall that P {κ,µ} 1

(τ ) = τ -ξ 1 , with ξ 1 given by [START_REF] Levant | Higher-order sliding modes, differentiation and output-feedback control[END_REF]. Plugging this expression of g into (58) immediately shows that (28) is valid for q = 1. Suppose (28) holds for a given q > 1. Then, we obtain from (52),

x(n) (κ, µ; n + q + 1) = λ (59)


Bertrand Pourroy

Christophe Curti

Edouard Lamy

Jauhel N'gbesso

Innocent Krasse

Patrice Vanelle

Availability of Morphine Oral Solution for Childhood Cancer Patients in Low-Income Countries: Compounding and Stability Study in a Cote d'Ivoire University Teaching Hospital

HAL is

Availability of Morphine Oral Solution for Childhood Cancer Patients in Low-Income Countries: Compounding and Stability Study in a Cote d'Ivoire University Teaching Hospital

To the Editor Cancer is an important cause of morbidity and mortality worldwide. In Africa, focusing on childhood cancer showed a constant increase in the annual number of new cases and a very low overall survival rate. [START_REF] Couitchere | Childhood cancer: epidemiology in the pediatric oncology department of a university hospital center in Abidjan, Ivory Coast[END_REF] Cancer treatment is often associated with pain management. The World Health Organization's guidelines recommend oral morphine (hydrochloride or sulfate) for severe chronic pain treatment in children. A previous work evidenced the good feasibility of oral morphine solution in Uganda. [START_REF] Merriman | Pain control in the African context: the Ugandan introduction of affordable morphine to relieve suffering at the end of life[END_REF] However, the authors highlighted several limitations to successfully replicate this pain management in other African countries. In low-income and middle-income countries, such as Côte d'Ivoire, limitations include both expense and restricted access to active pharmaceutical ingredients. This lack of access can prevent pain relief and a peaceful end of life for children suffering from cancer. [START_REF] Harding | How can we improve palliative care patient outcomes in low-and middle-income countries? Successful outcomes research in sub-Saharan Africa[END_REF] Low concentrations of morphine solutions for oral administration are not commercially available in several countries. Various compounding formulations were developed, but only a few have been reported in the literature. 4e7 The purpose of our study was to develop an affordable and easily available oral morphine formulation for pediatric patients with cancer pain. We devised an easy-to-make and low-cost morphine syrup formulation using local African (Côte d'Ivoire) raw materials and commercial injectable morphine. In this report, we describe the results of a study to assess the required quality controls and stability of this formulation.

Methods

Reagents and Instrumentation

Ampoules of morphine reference standard (Cerilliant Corporation, Round Rock, TX) at 1.00 mg/mL in methanol were used for method validation, and morphine hydrochloride trihydrate (Aguettant, Lyon, France) of pharmaceutical grade was used for the forced degradation study. The high-performance liquid chromatography (HPLC) mobile phases were prepared using ultrapure water (HiPerSolv Chromanorm Ò ; VWR International, Radnor, PA), methanol (HiPerSolv Chromanorm Ò ; VWR International), and acetonitrile (HiPerSolv Chromanorm Ò ; VWR International) of HPLC grade. Sodium lauryl sulfate buffer was prepared from sodium lauryl sulfate (0.65 g/L; Rectapur Ò ; Prolabo, Sion, Switzerland) and sodium sulfate (1.625 g/L; Merck, Darmstadt, Germany) in ultrapure water adjusted to pH 2.5 with 0.1 M of sulfuric acid. The mobile phases were filtered using 0.45 mm of Millipore cellulose filters. Volumes were aliquoted with a precision pipette (Thermo Scientific Finnpipette Ò ; F2 500 mL, Thermo Scientific, Waltham, MA). pH was determined with a Thermo Scientific Orion 4 Star Ò pH meter, calibrated with Radiometer Analytical standard etalons Ò (pH 4.005, 7.000, and 10.012).

The chromatographic method was carried out on an automatic HPLC Dionex Ultimate 3000 Ò with a ultraviolet diode array detector. The apparatus was connected to an HP 1702 computer equipped with chromatographic data processing software (Chromeleon Ò Chromatography Management System, Version 6.80 SRH Biold 3161; 1994e2011 Dionex Corporation, Sunnyvale, CA). Morphine separation was achieved using a Stability Ò C18 column with a particle size of 5 mm (150 Â 4.6 mm).

Viscosity was determined with a viscometer (Brookfield Ò LVDV-1M), equipped with a thermostatically controlled low-volume adaptor and its mobile (SC4-31). Measures were done at 100 rpm, 25 C.

Evaluation of the suitability of microbiological analyses was realized with reference strains (Biom erieux Bioball Ò ; Biom erieux, Marcy L'Etoile, France): Staphylococcus aureus NCTC10788, Bacillus subtilis NCTC10400, Pseudomonas aeruginosa NCTC12924, Candida albicans NCPF3179, and Aspergillus braziliensis NCPF2275. Microbiological analyses were done on Tryptic Soy Agar (Biom erieux) and Sabouraud chloramphenicol gentamicin agar (Biom erieux) under a microbiological safety cabinet (Herasafe Ò KS; Thermo Scientific). Nine milliliter tubes of pharmacopeia diluent with Tween 80 (Dominique Dutcher, Brumath, France) and incubators (Heratherm Ò and Heraeus Ò Thermo Scientific) were also used.

Stability and Microbiological Methods

Morphine quantification by HPLC was based on a methodology previously described, [START_REF] Hoang | Stabilit e du chlorhydrate de morphine en solution dans les syst emes PCA portables a usage unique: interaction contenant-contenu [Stability of morphine hydrochloride solution in single use PCA systems: container-content interaction[END_REF] and forced degradation studies were conducted to prove the stability-indicating capability of our chromatographic method. Degradation products did not interfere with morphine determination. About 1.0 mg/mL of morphine syrup samples were treated as follows: 200 mL of syrup were mixed and vortexed with 1800 mL of ultrapure water to obtain a solution of theoretically 100 mg/mL of morphine hydrochloride trihydrate (approximately 75.92 mg/mL of morphine base form).

The stability of the morphine syrups was evaluated for six months, under several storage conditions: classical ambient temperature (25 C), hot ambient temperature (40 C), and refrigeration (2 C/8 C). We studied four parameters on three independent batches: morphine concentration (triplicate), pH (one batch analyzed), viscosity (triplicate), and microbiological contamination (triplicate for initial and final points and one batch analyzed for other points).

Microbiological analyses of morphine syrups were done by surface-spread method. Suitability of the method was done according to European Pharmacopeia standards with five reference strains (Staphylococcus aureus, Bacillus subtilis, Pseudomonas aeruginosa, Candida albicans, and Aspergillus braziliensis).

Results

Formulation

Morphine syrup was compounded from commercial parenteral morphine. This choice was supported by the fact that morphine sulfate and morphine hydrochloride active pharmaceutical ingredients are not only less readily available in developing countries (anecdotally, also in France) than parenteral marketed morphine but also because the stability and storage parameters of commercial glass vials can be easily managed. To ensure an easy-to-made formulation, only one excipient was chosen to have a single dilution step of parenteral morphine. Because microbiological stability could be problematic, a simple syrup was selected for the formulation. The cost and availability of the pharmaceutical commercial form of a simple syrup were important limitations in Côte d'Ivoire. A simple syrup was therefore extemporaneously compounded from local mineral water (C eleste Ò ; Ciprem-Ci, Bonoua, Côte d'Ivoire, calcium 59.2 mg/ L, magnesium 6.90 mg/L, bicarbonates 202.50 mg/ L, and chloride 14.2 mg/L) and local commercial sugars (Batch 1: Princesse Tatie Ò ; SUCAF, Ferkess edougou, Côte d'Ivoire, and Batch 2: Sucrivoire Ò ; SIFCA, Abidjan, Côte d'Ivoire). Briefly, 320 g of sugar were dissolved in 200 mL of water and heated at 105 C. Two batches were prepared from two different commercial sugars. From these two batches of simple syrup, three batches of 1 mg/mL of morphine syrups were compounded (from simple syrup Batch 1, from simple syrup Batch 2, and from a mixture in equal parts of the two simple syrup batches).

Glycerol and citric acid were added to the formulation to avoid morphine syrup hydration and ensure a low pH, respectively. Specifically, 54 g of glycerol (Fagron, Thiais, France) and 3.6 g of citric acid (Cooper, Melun, France) were mixed in 36 mL of morphine hydrochloride 10 mg/mL (Aguettant). Thereafter, simple syrup was added at 360 mL, and the morphine syrup was finally aliquoted in 30-mL glass bottles (Coop eration Pharmaceutique Franc ¸aise, Melun, France) (20 mL per bottle).

Stability and Microbiological Study

The initial values of the morphine hydrochloride formulation were of concentration 1.01 AE 0.03 mg/ mL, pH 2.38, and viscosity 162 AE 18 Cp. The International Conference on Harmonization defined significant changes as a 5% change in assay from its initial value for stability studies conducted in the pharmaceutical industry. This 5% limit is often increasing to a 10% change in assay from its initial value for hospital compounded preparations done in a smaller scale. [START_REF] Curti | Stability studies of five anti-infectious eye drops under exhaustive storage conditions[END_REF] In the present study, 10% limit were used to define significant changes. These batches were individually analyzed by quality control laboratory with regard to European Pharmacopeia specifications.

Under all the studied conditions, pH and viscosity remained stable, and morphine concentration was also found to be stable for four months under refrigerated temperature (2 C/8 C) and for one month at ambient (25 C) and hot ambient temperature (40 C). There were no decreases of morphine content in syrup and no significative increase of detectable degradation product for these time points investigated. Results are summarized in Table 1.

Total aerobic microbial count was found to be lower than 200 CFU/mL, and total combined yeasts and mold count was found to be lower than 20 CFU/mL. At all periods of the stability study and for each storage condition, microbiological contamination was found to be lower than the method detection limit (total aerobic microbial count and total combined yeasts and mold count lower than 10 CFU/mL).

Comment

This study demonstrates that a new formulation of oral morphine 1.00 mg/mL compounded with simple syrup made from Côte d'Ivoire local ingredients had no contamination and was stable for 120 days at þ4 C and 30 days at þ25 C and þ40 C. These results are in accordance with the literature. Several studies proved the long-term stability of morphine hydrochloride compounded as oral solutions. 4e7 In low-income countries, opioids are mainly available in solid dosage forms (capsule and tablet), but their cost is an important limitation to their use. [START_REF] Lima | Crosssectional pilot study to monitor the availability, dispensed prices, and affordability of opioids around the globe[END_REF] Drinkable oral morphine formulations are more suitable for pediatric uses. Several oral morphine formulations have been described, but they contain excipients that could be problematic in countries like Côte d'Ivoire. Problems include the unavailability of excipients and aqueous solutions at risk of microbiological contamination. Our formulation is a solution to these problems because a syrup form inhibits microbiological growth, and its formulation made from local ingredients decreases costs.

Pain management is a key element for patient suffering from cancer, including children. Drinkable morphine solutions have been developed for palliative care in Africa, 2 but there is a critical lack of oral formulations that can be used by children. Preliminary studies showed the necessity not only to train health care workers but also to develop morphine supply chain. Our formulation can reduce undertreatment of pain and delay or prevent the transition to parenteral administration. It can also improve acceptance by children because of its sugar taste and reduce nurse miscalculations by the use of a fixed concentration. It also eliminates the need for nurses to prepare oral morphine, which is time consuming and at risk of error.

Table 1

 1 pH, Viscosity, and Morphine a Content of Morphine Syrups AE 0.03 0.98 AE 0.05 1.02 AE 0.07 1.02 AE 0.04 1.04 AE 0.01

	Parameter	T0	1 Day	7 Days	15 Days	30 Days	45 Days	60 Days	90 Days	120 Days
	þ2/þ8 C
	pH	2.38	2.38	2.38	2.46	2.18	2.28	2.36	2.29	2.32
	Viscosity (Cp) 162 AE 18 Morphine a 1.1 AE 0.3	163 AE 20 1.3 AE 0.2	157 AE 19 1.6 AE 0.1	174 AE 19 1.7 AE 0.1	159 AE 17 1.3 AE 0.6	164 AE 17 160 AE 13 161 AE 12 1.3 AE 0.4 1.8 AE 0.1 1.0 AE 0.5 1.01 AE 0.05 157 AE 9
	(mg/mL)
	25 C
	pH	2.38	2.42	2.39	2.46	2.25	d	d	d	d
	Viscosity (Cp) 162 AE 18 Morphine a 1.01 AE 0.03 1.03 AE 0.02 1.03 AE 0.09 1.01 AE 0.08 1.08 AE 0.01 163 AE 17 158 AE 22 166 AE 16 157 AE 16	d	d	d	d
	(mg/mL)
	40 C
	pH	2.38	2.39	2.37	2.42	2.24	d	d	d	d
	Viscosity (Cp) 162 AE 18 Morphine a 1.01	162 AE 15	158 AE 19	163 AE 16	157 AE 16	d	d	d	d
	(mg/mL)

a Expressed in morphine hydrochloride trihydrate.

Bertrand Pourroy, PhD, PharmD Oncopharma Unit Pharmacy Department University Teaching Hospital of La Timone Marseille, France Groupe Franco-Africain d'Oncologie P ediatrique Institut Gustave Roussy Villejuif, France Christophe Curti, PhD, PharmD Pharmacy Department Service central de la qualit e et de l'information pharmaceutiques (SCQIP) AP-HM Marseille, France E-mail: christophe.curti@univ-amu.fr Aix-Marseille Univ CNRS Institut de Chimie Radicalaire ICR UMR 7273 Equipe de Pharmaco-Chimie Radicalaire Marseille, France Edouard Lamy, PhD, PharmD Pharmacy Department Service central de la qualit e et de l'information pharmaceutiques (SCQIP) AP-HM Marseille, France Aix-Marseille Univ UMR 7287 CNRS Institut des Sciences du Mouvement ISM Facult e des Sciences du Sport Marseille, France Jauhel N'Gbesso, CPhT Innocent Krasse, PharmD Pharmacy Department University Teaching Hospital of Treichville Abidjan, Cote d'Ivoire Patrice Vanelle, PhD, PharmD Pharmacy Department Service central de la qualit e et de l'information pharmaceutiques (SCQIP) AP-HM Marseille, France Aix-Marseille Univ CNRS Institut de Chimie Radicalaire ICR UMR 7273 Equipe de Pharmaco-Chimie Radicalaire Marseille, France

Disclosures and Acknowledgments

This research received no specific funding/grant from any funding agency in the public, commercial, or not-for-profit sectors. The authors declare no conflicts of interest.


Alexandre Baud

Xavier Flecher

Richard Alexandre Rochwerger

Jean-Camille Mattei

Jean Noël Argenson

Comparing the outcomes of the induced membrane technique between the tibia and femur: Retrospective single-center study of 33 patients

Keywords:

Introduction: Bone defects are challenging to treat surgically. The primary objective of our study was to compare the union rate and time to union between the tibia and femur when using the induced membrane technique. The secondary objective was to document how failures were managed. Material and Methods: This retrospective, single-center study involved 33 patients (23 men, 10 women) who were older than 18 years of age. They were treated surgically for a leg fracture or long bone nonunion (22 tibia, 11 femur) using the induced membrane technique between January 2011 and December 2016 and had a complete follow-up. The minimum follow-up was 1 year for fractures and 2 years for non-union cases. Bone union was defined as the presence of at least two cortices with bridging on two radiographic views and return to full weight bearing. Results: The mean patient age was 38.3 ± 15.5 years (18-72). The mean bone defect size was 7.9 ± 5.0 cm (2.3-18.0). The mean follow-up was 3.3 ± 1.8 years (1-7.2). The union rate was 61% (20 patients). The mean time to union was 10 ± 6.4 months . The time to union was significantly longer in the tibia (11.6 ± 6.9 months ) than in the femur (6.3 ± 2.9 months [3.4-10.3]) (p = 0.025). The failure rate did not differ between the tibia and femur. Nine of the 13 patients (69%) in which the treatment failed were reoperated; 7 of them underwent nonunion treatment (78%) and 2 underwent amputation (22%). The other 4 patients were waiting for an infection to resolve before being reoperated.

Conclusion:

The induced membrane technique is an effective surgical procedure for large bone defects in both the tibia and femur. However, the time to union was shorter in the femur than the tibia in our cohort. Level of evidence: IV, retrospective study.

Introduction

Limb preservation in the context of large bone defects has always been a priority for surgeons, although it may be controversial. Nevertheless, the long-term outcomes are very encouraging [START_REF] Masquelet | Very long-term results of post-traumatic bone defect reconstruction by the induced membrane technique[END_REF]. Reconstruction is a challenging surgery [START_REF] Obremskey | Current practice in the management of open fractures among orthopaedic trauma surgeons. Part A: initial management. A survey of orthopaedic trauma surgeons[END_REF][START_REF] Obremskey | Current practice in the management of open fractures among orthopaedic trauma surgeons. Part B: Management of Segmental Long Bone Defects. A Survey of Orthopaedic Trauma Association Members[END_REF] and requires multidisciplinary [START_REF] Nauth | Managing bone defects[END_REF], early and effective treatment of infections [START_REF] Patzakis | Factors influencing infection rate in open fracture wounds[END_REF][START_REF] Yusof | Outcomes of infected grade IIIB open tibial fractures[END_REF][START_REF] Gustilo | Prevention of infection in the treatment of one thousand and twenty-five open fractures of long bones: retrospective and prospective analyses[END_REF] and tissue lesions [START_REF] Masquelet | Muscle reconstruction in reconstructive surgery: soft tissue repair and long bone reconstruction[END_REF]. However, there is no consensus on the management of large bone defects [START_REF] Hake | Difficulties and challenges to diagnose and treat post-traumatic long bone osteomyelitis[END_REF][START_REF] Ashman | Treatment of non-unions with bone defects: which option and why?[END_REF]. Several surgical techniques have been described in the literature, such as Ilizarov distractionosteogenesis [START_REF] Cattaneo | The treatment of infected nonunions and segmental defects of the tibia by the methods of Ilizarov[END_REF], bone transport [START_REF] Paley | Ilizarov bone transport treatment for tibial defects[END_REF], transfer of vascularized bone (with the fibula being the most used graft) [START_REF] Wood | Free vascularized fibular grafting-25 years' experience: tips, techniques, and pearls[END_REF], and segmental banked bone allograft [START_REF] Jaffe | Massive bone allografts for traumatic skeletal defects[END_REF].

In the 1980s, Masquelet described a novel bone reconstruction techniques in which polymethylmethacrylate (PMMA) was placed in the defect to induce a foreign body reaction and formation of an induced membrane [START_REF] Masquelet | Muscle reconstruction in reconstructive surgery: soft tissue repair and long bone reconstruction[END_REF][START_REF] Masquelet | Reconstruction of the long bones by the induced membrane and spongy autograft[END_REF][START_REF] Masquelet | The concept of induced membrane for reconstruction of long bone defects[END_REF] having osteogenic properties: secretion of growth factors, recruitment and osteogenic differentiation of stem cells, and production of angiogenic factors contributing to vascularization and corticalization of the bone graft [START_REF] Pelissier | Induced membranes secrete growth factors including vascular and osteoinductive factors and could stimulate bone regeneration[END_REF][START_REF] Liu | Histological characteristics of induced membranes in subcutaneous, intramuscular sites and bone defect[END_REF][START_REF] Gruber | Osteogenic, stem cell and molecular characterisation of the human induced membrane from extremity bone defects[END_REF][START_REF] Aho | The mechanism of action of induced membranes in bone repair[END_REF]. The surgical induced membrane technique is one of the most effective strategies for filling segmental bone defects with reported union rates of 82% to 100% [START_REF] Taylor | Treatment of bone loss with the induced membrane technique: techniques and outcomes[END_REF][START_REF] Mauffrey | Reconstruction of long bone infections using the induced membrane technique: tips and tricks[END_REF][START_REF] Wong | Masquelet technique for treatment of posttraumatic bone defects[END_REF].

Recent publications have focused on the biological and mechanical mechanisms that contribute to bone regeneration. Based on the work of Giannoudis et al. and the "diamond concept", five essential factors have been identified: osteogenic stem cells, osteoconductive scaffold, growth factors, stable mechanical environment and good vascularization [START_REF] Giannoudis | Fracture healing: the diamond concept[END_REF][START_REF] Giannoudis | The diamond conceptopen questions[END_REF]. Bone reconstruction using the induced membrane technique that also complies with the diamond concept will optimize the bone regeneration opportunity [START_REF] Yee | Scientific understanding of the induced membrane technique: current status and future directions[END_REF][START_REF] Aurégan | Induced membrane for treatment of critical sized bone defect: a review of experimental and clinical experiences[END_REF]. Some authors have pointed out that the outcomes differ depending on where the defect is located on a given bone [START_REF] Aurégan | Induced membrane for treatment of critical sized bone defect: a review of experimental and clinical experiences[END_REF]. To our knowledge, no published study has documented how the outcomes of this surgical technique differ between two different bones. We hypothesized that the effectiveness of the induced membrane technique (time to union and union rate) will vary depending on the location of the bone defect. The aim of this study was to compare the time to union and union rate between tibia and femur locations and to describe how failures are managed.

Materials and methods

Cohort description

We carried out a retrospective study of all the patients with extensive posttraumatic bone loss due to fracture or non-union who were seen at the assistance publique hôpitaux de Marseille (AP-HM) between January 2011 and December 2016. Two groups were compared: tibia versus femur. All the patients had undergone surgical treatment using the induced membrane technique [START_REF] Masquelet | Reconstruction of the long bones by the induced membrane and spongy autograft[END_REF]. The minimum patient age was 18 years. A comprehensive followup was done throughout the treatment period. The minimum follow-up was 1 year for fracture cases and 2 years for non-union cases. Exclusion criteria were the use of a vascularized fibular graft during the second surgical phase, pathological fractures and bone defects related to gunshot wounds. The primary outcome measures were the time to union and union rate; the secondary outcome measure was an analysis of the management of failed bone grafts in these two locations.

Induced membrane surgical technique

The 1 st stage consisted of wide excision of damaged or infected tissues, with addition of a musculocutaneous coverage flap in some cases. Next, the fracture site was stabilized with an external fixator or internal fixation. The bone defect was filled with a Gentamycinloaded PMMA cement spacer (Palacos ® R + G, Zimmer-Biomet, Warsaw, USA). Several microbiological samples were harvested to determine which bacteria were present in view of treating surgical site infections. The 2 nd phase was performed after infections were resolved and the soft tissues had healed completely, approximately 6-8 weeks later. This procedure consisted of removing the cement spacer and filling the gap with bone graft and in some cases, additional bone substitutes. New bone fixation was performed at this point. The options for fixation consisted of an external fixator, internal fixation using a screw plate or intramedullary locked nail, and Sarmiento ® type rigid splints.

Two types of bone grafts were used in our cases: bone allograft from a tissue bank and autograft obtained by either reaming cancellous bone (RIA ® , DePuy-Synthes, Switzerland) in the ipsilateral or contralateral limb, or by harvesting corticocancellous bone from the iliac crest.

Data collection

After obtaining the patients' consent, information to support our data analysis was collected (Table 1). The patients were reviewed in person at last 1 year after fracture treatment and at least 2 years after nonunion treatment. All the data were stored in a computerized medical records system (Axigate ® ) at the AP-HM, namely the information collected during the postoperative follow-up visits at 6 weeks, 3 months, 6 months and then every year. This study was approved by our hospital's research ethics board and data processing was approved by the hospitals' data privacy officer. The data collected were compared between the tibia and femur groups (Tables 2 and3).

Study population

The data collected consisted of age at the time of the accident, sex, smoking habits, body mass index (BMI), bone involved (tibia or femur), surgical indication (fracture, nonunion), size of bone defect using the Karger classification [START_REF] Karger | French Society of Orthopaedic Surgery and Traumatology (SoFCOT). Treatment of posttraumatic bone defects by the induced membrane technique[END_REF] -group 1 (> 2 cm), group 2 (2-5 cm), group 3 (5-10 cm) and group 4 (> 10 cm)-determined using Centricity ® software, whether or not the fracture site was open using the Gustilo and Anderson classification [START_REF] Thakore | The Gustilo-Anderson classification system as predictor of nonunion and infection in open tibia fractures[END_REF], risk factors for nonunion using the Non-Union Scoring System (NUSS) [START_REF] Calori | Classification of non-union: need for a new scoring system?[END_REF][START_REF] Calori | Validation of the Non-Union Scoring System in 300 long bone non-unions[END_REF], and the preoperative infection rate based on the presence of at least one bacterium at the surgical site identified by microbiology and treated with antibiotics before the surgical procedure.

Follow-up period

During the follow-up period, we recorded the time elapsed between the 1 st stage and the fracture, time elapsed between the cement spacer implantation and bone grafting, type of fixation (external fixator, internal fixation or rigid splint) used during each stage, whether a coverage flap was needed at any stage, postoperative infection rate defined by the presence of at least one microorganism found at the surgical site during microbiology sampling and treated by postoperative antibiotics, type of bone graft used during the second surgical stage (RIA ® , iliac crest, allograft) and lastly, the length of hospital stay.

Primary endpoint

Bone union was defined using radiological and clinical criteria: presence of at least two continuous bone bridges on two radiographic views or CT scan and return to full weight bearing. The time to union was defined as the time elapsed between the application of the bone graft and the date bone union was achieved. The union rate was documented in cases that required reoperation after the second stage of surgery along with the time to union and was compared to cases that did not require reoperation.

Secondary endpoint

The main complication of the induced membrane technique is nonunion, which is defined as lack of bone healing on radiographs and the patient being unable to resume walking with full weight bearing at least 1 year after the bone graft was introduced. Among these patients, we determined the reoperation rate: new treatment of nonunion with bone graft or amputation; union rate after revision surgery; number of patients still waiting for treatment because of an ongoing infection (which makes it impossible to perform surgery).

Statistical methods

A univariate analysis was carried out using Fisher's exact test for qualitative variables and Student's t-test with Welch's correction for quantitative variables using GraphPad Prism ® software. The significance threshold was set at p < 0.05. (Gustilo 2) that was treated at another hospital. External fixator applied, then plate fixation + bone graft (iliac crest). After 2 years: infected nonunion with Enterococcus faecalis identified. Patient transferred to our hospital for continuation of care; b: AP and lateral view of the knee after the 1 st stage of the Masquelet procedure showing an approximately 8 cm defect which was filled with a cement spacer and stabilized with a blade-plate; c: AP and lateral views of the femur after the second stage (6 week interval) with minimal bone margin resection (especially anterior cortex) followed by plate fixation and cancellous bone grafting (RIA ® contralateral); d: AP and lateral views of the knee after hardware removed 1½ years after grafting. Bone union was achieved 3½ months after grafting.

Results

Study cohort

The study cohort consisted of 33 patients (23 men, 10 women) who were treated with the induced membrane technique (Tables 2 and3). The procedure was performed at the tibia in 22 cases and the femur in 11 cases. There were more men in the tibia group than the femur group: 19 patients (86.4%) versus 4 patients (36.4%), p = 0.006. The mean age at the time of the accident was 38.3 ± 15.5 years . The reason for the surgical treatment was non-union in 15 patients and fracture in 18 patients. The mean defect size was 7.9 cm ± 5.0 cm (2.3-18.0). An open fracture was present in 28 patients (84.8%).

Before the first surgical stage, the tibia and femur groups were comparable in terms of age at the time of the accident, cardiovascular risk factors, reason for the surgery (fracture or nonunion), size of bone defect, percentage of open fractures, preoperative infection rate. The NUSS did not differ between groups: 37.8 ± 12.7 (20-60) for the tibia versus 39.1 ± 13.8 (20-68) for the femur, p = 0.7. However, there were more open Gustilo 3b fractures in the tibia group than in the femur group: 0 patients versus 8 patients, p = 0.03.

The mean follow-up was 3.3 ± 1.8 years (1-7.2) (Table 1). It was 2.7 ± 0.6 years (2-3.3) in the femur group and 3.5 ± 2.1 years (1-7.25) in the tibia group. This difference between groups was not significant.

The tibia and femur groups were comparable in terms of the surgery characteristics except for plate fixation being used more often during the 1 st stage in the femur group than the tibia group: 4 patients versus 0 patients, p = 0.008. The infection rate after the 1 st surgical stage was significantly higher in the tibia group than in the femur group: 13 patients (59%) versus 2 patients (18%), p = 0.03. The infection rate after the 2 nd surgical stage did not differ between groups.

Primary endpoint

The initial union rate was 61% (20 of 33 patients) (Table 4). The union rate was similar between the tibia and femur: 14 patients (64%) versus 6 patients (54%).

The mean time to union was 10.0 ± 6.4 months (3-23) (Table 1). The time to union was significantly higher in the tibia group than in the femur: 11.6 ± 6.9 months (3-23) versus 6.3 ± 2.9 months (3.4-10.3), p = 0.025.

Secondary endpoint

The union rate after surgical revision for nonunion was 76% (25 of 33 patients). The union rate was similar between the tibia and the femur: 18 patients (84%) versus 7 patients (64%). The mean time to union after surgical revision for nonunion was 13.6 ± 9.6 months (3-35) (Table 1). The time to union was significantly higher in the tibia group than in the femur: 15.8 ± 10.2 months (3-35) versus 7.4 ± 4.9 months (3.4-17), p = 0.015 (Table 4).

The non-union rate after the induced membrane treatment was 39% (13 patients). The non-union rate did not differ between the tibia and femur sites. Among the patients who had a non-union, 9 patients (69.2%) were reoperated: 7 underwent another non-union treatment (77.8%) and 2 underwent amputation (22.2%). Of the 7 patients who were reoperated with a new bone graft, 5 achieved bone union. Four patients (31%) were still waiting for an infection to resolve as of our last review point. There was no difference in the management of failures between the tibia and femur locations (Table 5) (Figs. 123).

Discussion

The induced membrane technique pioneered by Masquelet has between successfully used for many years for bone reconstruction when a single-stage procedure is not feasible [START_REF] Masquelet | Very long-term results of post-traumatic bone defect reconstruction by the induced membrane technique[END_REF][START_REF] Yee | Scientific understanding of the induced membrane technique: current status and future directions[END_REF][START_REF] Karger | French Society of Orthopaedic Surgery and Traumatology (SoFCOT). Treatment of posttraumatic bone defects by the induced membrane technique[END_REF]. The primary aim of our study was to compare the outcomes of this technique in two different bones: tibia and femur. There was no difference in the union rate between the two sites. However, the time to union was significantly longer in the tibia than the femur.

The focus of our study was on the union rate, which was 61% (20 of 33 patients) (Tables 4 and5). Our results are consistent with published studies on this topic reporting union rates from 41.6% to 100% [START_REF] Morelli | Masquelet technique: myth or reality? A systematic review and meta-analysis[END_REF]. Karger et al. [START_REF] Karger | French Society of Orthopaedic Surgery and Traumatology (SoFCOT). Treatment of posttraumatic bone defects by the induced membrane technique[END_REF] reported a 90% union rate in 84 patients with a mean union time of 13.2 months (57 weeks). Giannoudis et al. [START_REF] Giannoudis | Restoration of long bone defects treated with the induced membrane technique: protocol and outcomes[END_REF] found a 98% union rate in a 42-patient cohort with a mean time of 5.3 months (23 weeks). However, other teams have reported lower union rates. In a recent study of 12 patients, Morris et al. [START_REF] Morris | Induced membrane technique for treating tibial defects gives mixed results[END_REF] reported a 42% union rate. The failures were due to the high number of preoperative infections (11 of 12 patients). In our study, there was no difference in the preoperative infection rate between the tibia and femur: 45% (10 patients) versus 54% (6 patients), p = 0.72.

The non-union rate after the induced membrane technique was 39% (13 patients) in our study. Masquelet and his team were the first to describe the complications related to his surgical technique: 5 cases of infected non-union after the second stage among 35 patients [START_REF] Masquelet | Reconstruction of the long bones by the induced membrane and spongy autograft[END_REF]. Azi et al. confirmed that infection was the main reason for failure of this surgical technique, as they found a 20% delayed infection rate (7 of 34 patients) [START_REF] Azi | Membrane induced osteogenesis in the management of posttraumatic bone defects[END_REF]. Apard et al. reported a 42% infection rate after the second phase [START_REF] Apard | Two-stage reconstruction of post-traumatic segmental tibia bone loss with nailing[END_REF]. El-Alfy et al. reported a 30% non-union rate [START_REF] El-Alfy | Management of segmental skeletal defects by the induced membrane technique[END_REF], while it was 17% in a study by Taylor et al. [START_REF] Calori | Validation of the Non-Union Scoring System in 300 long bone non-unions[END_REF]. The latter identified three risk factors for non-union: post-graft infection (relative risk = 11.8, 95% [4.4-31.9]), wound dehiscence (relative risk = 9.7, 95% [4.3-22.1]) and presence of an infection before the first surgical phase (relative risk = 2.9, 95% [1.1-7.9]). We could not determine which factors contributed to the non-healing in our study since the small sample size precluded a multivariate analysis. However, there was a significant difference between the tibia and femur groups in terms of the number of open fractures and infection rate after the first surgery phase. There were more Gustilo 3b open fractures in the tibia group than in the femur group (8 patients versus 0 patients, p = 0.03) while the infection rate after the 1 st phase was higher in the tibia group than in the femur group (59% versus 18%, p = 0.03).

Among the 13 patients whose bone defect did not heal after the procedure, 2 underwent an amputation due to chronic osteomyelitis despite antibiotic therapy and several surgical revision procedures. In their study, Karger et al. [START_REF] Karger | French Society of Orthopaedic Surgery and Traumatology (SoFCOT). Treatment of posttraumatic bone defects by the induced membrane technique[END_REF] reported 6 amputations out of the 8 patients who failed treatment. In their meta-analysis, Morelli et al. [START_REF] Morelli | Masquelet technique: myth or reality? A systematic review and meta-analysis[END_REF] found that at least 1 amputation was reported in 11 of 17 publications about the Masquelet procedure. Our study confirms that the amputation rate is not insignificant after failed treatment with the induced membrane technique.

Conclusion

In summary, the induced membrane technique described by Masquelet is an effective treatment solution for managing large bone defects in the lower limbs of adults, with a satisfactory union rate, despite several postoperative complications, especially infections. We found that this technique can be used effectively to treat both fractures and midshaft non-union of the tibia and femur. The time to union appeared to be less in the femur than the tibia. This new information gives surgeons a better idea of the time to union for different bone locations, which they can communicate to their patients to prepare them for the long treatment course.

But we have asked ourselves whether a different surgical technique would be more effective than the induced membrane technique? In our study, this procedure was less effective at the tibia because the time to union was longer than in the femur. In these cases, what is the role of other surgical techniques such as the Ilizarov, bone transport or vascularized bone graft? These questions are good avenues for future research into optimizing treatment of large bone defects.

Disclosure of interest

The authors declare that they have no competing interest.

Fig. 1 .

 1 Fig. 1. Male patient, 55 years old, who was treated for infected nonunion of the tibia. a: AP and lateral radiographs of the tibia after application of external fixator in the emergency room for an open (Gustilo 3a) fracture of the distal end of the tibia; progression to infected nonunion; b: photograph of the 1 st stage of the Masquelet technique after resection of margins (left) and after implanting the cement spacer (right) 1 year after the fracture event. Bone defect was 4.5 cm long and 2 cm wide; c: AP and lateral radiographs after the 1 st stage of surgery with plate fixation and cement spacer; d: Radiographs made 16 months after the second phase (ORIF: plate and intramedullary nailing; bone graft: RIA ® ) showed bone union.

Fig. 2 .

 2 Fig. 2. Female patient, 42 years old, who was treated because of an infected nonunion of a prior femur fracture. a: preoperative AP and lateral radiographs of the knee after motor vehicle accident showing an open distal femur fracture (AO type C3)(Gustilo 2) that was treated at another hospital. External fixator applied, then plate fixation + bone graft (iliac crest). After 2 years: infected nonunion with Enterococcus faecalis identified. Patient transferred to our hospital for continuation of care; b: AP and lateral view of the knee after the 1 st stage of the Masquelet procedure showing an approximately 8 cm defect which was filled with a cement spacer and stabilized with a blade-plate; c: AP and lateral views of the femur after the second stage (6 week interval) with minimal bone margin resection (especially anterior cortex) followed by plate fixation and cancellous bone grafting (RIA ® contralateral); d: AP and lateral views of the knee after hardware removed 1½ years after grafting. Bone union was achieved 3½ months after grafting.

Fig. 3 .

 3 Fig. 3. Male patient, 20 years old, who was treated for a tibial metaphysis fracture. a: postoperative AP radiograph taken after motor vehicle accident showing an open (Gustilo 2) proximal tibial shaft fracture with an approximate 40 mm bone defect. First stage of the induced membrane procedure done in the emergency room: implantation of cement spacer and stabilization with external fixator; b: AP radiograph of the proximal tibia after the 2 nd stage of surgery showing the corticocancellous graft from the iliac crest and cancellous graft harvested with RIA ® along with locking plate fixation 3 months after the injury event; c, d: AP and lateral views of the tibia showing bone union. The patient could walk without pain or crutches and had full weightbearing 12 months after the accident.

Table 1

 1 Study population.

		Number of patients	%	Mean	Median	Standard Deviation	Min	Max
	Sex
	Male	23	69.7
	Female	10	30.3
	Age at time of accident (years)			38.3	37	15.5	18	72
	Cardiovascular risk factors
	Smoking	20	60.6
	Body mass index (kg/m 2 )			25.8	25	6.07	12.5	36.8
	Bone
	Tibia	22	66.7
	Femur	11	33.3
	Reason for procedure
	Nonunion	15	45.5
	Fracture	18	54.5
	Size of bone defect (cm)			7.9	6.2	5	2.3	18
	Classification of bone defect
	Karger 1 (< 2 cm)	0	0.0
	Karger 2 (2-5 cm)	12	36.4
	Karger 3 (5-10 cm)	11	33.3
	Karger 4 (> 10 cm)	10	30.3
	Number of open fractures	28	84.8
	Open fracture classification system
	Gustillo 1	2	6.1
	Gustillo 2	11	33.3
	Gustillo 3a	7	21.2
	Gustillo 3b	8	24.2
	Preoperative infection rate	16	48.5
	Non-Union Scoring System (NUSS)			38.7	36	13.24	20	68
	Time interval: Fracture-1st stage of surgery (months)			9.4	1.8	13.90	0	65.6
	Type of fixation
	External fixator	25	75.8
	Plate	4	12.1
	Rigid splint	4	12.1
	Additional surgical procedures	16	48.5
	Type of plastic surgery
	Gastrocnemius flap	7	21.2
	Latissimus dorsi flap	4	12.1
	Free thigh flap	3	9.1
	Skin graft	2	6.1
	Postoperative infection rate	15	45.5
	Length of hospital stay (days)			25.2	13.5	29.4	3	107
	Time interval between two stages of surgery (months)			4.1	3.2	3.2	1.3	18.8
	Type of fixation
	Intramedullary nailing	10	30.3
	Plate	22	66.7
	External fixator	5	15.2
	Rigid splint	2	6.1
	Bone graft
	RIA	20	60.6
	Iliac crest	18	54.5
	Allograft	8	24.2
	Postoperative infection rate	9	27.3
	Length of hospital stay (days)			8.2	8	3.1	4	15
	Union rate
	Without surgical revision	20	60.6
	After surgical revision	25	75.7
	Time interval 2 nd stage-Union (months)
	Without surgical revision			10	8.9	6.4	3	23
	After surgical revision			13.6	10.3	9.6	3	35
	Follow-up (years)			3.3	2.6	1.8	1.00	7.2

Table 2

 2 Preoperative data for the femur and tibia groups.

		FEMUR (n = 11)					TIBIA (n = 22)
		Number of	%	Mean Median Standard	Min Max Number of	%	Mean Median Standard	Min Max p
		patients				deviation		patients				deviation
	Sex
	Male	4	36.4				19	86.4			0.006
	Female	7	63.6				3	14			0.006
	Age at time of accident (years)			37.2	36	17.6	18	72		39.0	39	14.81	19	67
	Cardiovascular risk factors
	Smoking	8	72.7				12	54.5
	Body mass index (kg/m 2 )			26.3	26	8.4	13	36.8		25.4	24.9	4.03	18	33.5
	Reason for procedure
	Nonunion	6	54.5				9	40.9
	Fracture	5	45.5				13	59
	Size of bone defect (cm)			9.9	8	5.1	4.3 18		6.9	5.1	4.80	2.3 18
	Classification of bone defect
	Karger 1 (< 2 cm)	0	0.0					0	0.0
	Karger 2 (2-5 cm)	2	18.2				10	45.5
	Karger 3 (5-10 cm)	4	36.4				7	31.8
	Karger 4 (> 10 cm)	5	45.5				5	22.7
	Number of open fractures	8	72.7				20	90.9
	Open fracture classification system
	Gustillo 1	1	9.1					1	4.5
	Gustillo 2	4	36.4				7	31.8
	Gustillo 3a	3	27.3				4	18.2
	Gustillo 3b	0	0.0					8	36.4			0.03
	Preoperative infection rate	6	54.5				10	45.5
	Non-Union Scoring System (NUSS)		37.8	38	12.70	20	60		39.1	36	13.80	20	68
	Table 3
	Intraoperative data for the femur and tibia groups.
		FEMUR (n = 11)					TIBIA (n = 22)
		Number of	%	Mean Median Standard	Min Max Number of	%	Mean Median Standard	Min Max p
		patients				deviation		patients				deviation
	Time interval: Fracture-1st			15.5	12.6	19.7	0	65.7		6.3	1.7	8.90	0	30.2
	stage of surgery (months)
	Type of fixation
	External fixator	7	63.6					18	81.8
	Plate	4	36.4					0	0.0			0.008
	Rigid splint	0	0.0					4	18.2
	Additional surgical procedures	3	27.3					13	59.1
	Type of plastic surgery
	Gastrocnemius flap	0	0.0					7	31.8
	Latissimus dorsi flap	1	9.1					3	13.6
	Free thigh flap	0	0.0					3	13.6
	Skin graft	2	18.2					0	0.0
	Postoperative infection rate	2	18.2					13	59.1			0.03
	Length of hospital stay (days)			19.8	15	16.3	4	59		28	12	34	3	107
	Time interval between two			3.4	3.3	1.5	1.5	6.7		4.5	3.2	3.7	1.3	18.8
	stages of surgery (months)
	Type of fixation
	Intramedullary nailing	2	18.2					8	36.4
	Plate	10	90.9					12	54.5
	External fixator	0	0.0					5	22.7
	Rigid splint	0	0.0					2	9.1
	Bone graft
	RIA	8	72.7					12	54.5
	Iliac crest	4	36.4					14	63.6
	Allograft	3	27.3					5	22.7
	Postoperative infection rate	2	18.2					7	31.8
	Length of hospital stay (days)			8.3	9	2.7	4	13		8.2	8.0	3.4	4.0	15.0

Table 4

 4 Union rate and time to union after the 2nd phase of surgery in the patients with successful and failed induced membrane treatment.

		FEMUR (n = 11)					TIBIA (n = 22)
		Number of	%	Mean Median Standard	Min Max Number of	%	Mean Median Standard	Min Max p
		patients				deviation		patients				deviation
	Union rate
	Without surgical revision	6	54.5					14	63.6
	After surgical revision	7	63.6					18	84
	Time interval: 2 nd stage-Union (months)
	Without surgical revision			6.3	5.3	2.9	3.4 10.3		11.6 10.3	6.9	3	23	0.025
	After surgical revision			7.4	6.4	4.9	3.4 17		15.8 16	10.2	3	35	0.015
	Follow-up (years)			2.7	2.6	0.6	2	3.3		3.5	3.3	2.1	1	7.25

Table 5

 5 Outcomes of failed cases of induced membrane treatment.

		Population (n = 33)		Femur (n = 11)		Tibia (n = 22)
		Number of patients	%	Number of patients	%	Number of patients	%	P
	Nonunion rate	13	39.4	5	45.5	8	36.4
	Number of reoperations	9	27.3	2	40	7	87.5
	Nonunion treatment	7	77.8	2	100	5	71.4
	Amputation	2	22.2	0	0	2	28.6
	Union rate after revision	5	55.5	1	50	4	57.1
	Patients waiting for revision surgery	4	12.1	3	60	1	12.5

Funding

None.

Author contributions

A. Baud: data collection, writing of article. X. Flecher: study coordinator; helped to write article. R.A. Rochwerger: critical review of article. J.C. Mattéi: final critical review of article. J.N. Argenson: final approval of article.