Anastasia Volkova

Christoph Lauter

Thibault Hilaire

Reliable verification of digital implemented filters against frequency specifications

Reliable implementation of digital filters in finite-precision is based on accurate error analysis. However, a small error in the time domain does not guarantee that the implemented filter verifies the initial band specifications in the frequency domain. We propose a novel certified algorithm for the verification of a filter's transfer function, or of an existing finite-precision implementation. We show that this problem boils down to the verification of bounds on a rational function, and further to the positivity of a polynomial. Our algorithm has reasonable runtime efficiency to be used as a criterion in large implementation space explorations. We ensure that there are no false positives but false negative answers may occur. For negative answers we give a tight bound on the margin of acceptable specifications. We demonstrate application of our algorithm to the comparison of various finite-precision implementations of filters already fully designed.

I. INTRODUCTION

The great majority of signal processing or control algorithms are designed on our desktop computer, and then embedded in digital devices, such as general purposes processors, DSPs, FPGAs or ASICs.

Linear filters are basic bricks of such algorithms, and they are usually designed using software like Matlab 1 or SciPy2 which are used to design filters fulfilling a given frequency specification. These tools rely on algorithms based on double precision Floating-Point arithmetic. Then, the filter is implemented, often with Fixed-Point arithmetic, for cost or power consumption reasons. The filter algorithm used for the implementation may directly round and use the coefficients obtained at design stage, or deduce new coefficients to round and use, from the previous stage. In both cases, the finite precision used in design and implementation stages induces some uncontrolled errors that may make the implemented filter/controllers differ from the initial specification. Currently, engineers have no other choice than to add "some" design margin in their specification, trust these software tools and perform some simulations to somehow test that the implemented filter respects its specifications.

So this article deals with the reliable a posteriori verification of the filter implementation against its mathematical frequency specification. We understand by the term a posteriori verification that all filter design and implementation is already finished at that point. A filter comes with its (possibly quantized) coefficients and its filter order and it is not up to our algorithm to question these input data.

This problem cannot straightforwardly be solved in a rigorous way by computing the maximum and minimum frequency response and comparing it to a bound. Computing such an extremum value rigorously amounts to rigorously optimizing a function on a subset of the complex numbers, which is not an easy task, even using multiple-precision interval arithmetic [START_REF] Chevillard | Efficient and accurate computation of upper bounds of approximation errors[END_REF].

We shall rather consider verification only and decompose that verification problem into the following two parts: First, suppose, we have a transfer function with its coefficients expressed as fixed-or Floating-Point numbers at some precision (integers scaled by powers of 2) and a set of band specifications. Then, we need to verify whether this filter satisfies the specifications. Such a problem boils down to verification whether a rational function is between bounds. We propose an algorithm that guarantees that no false positive answer is ever given.

Secondly, suppose, we have an algorithm realizing the filter, defined as a set of equations or as a data-flow graph (like those in Matlab/Simulink). In general case, the algorithm's coefficients have been quantized and are hence not those of the transfer function. Therefore, to verify whether the filter satisfies band specifications, its transfer function must be recovered. Usually, this cannot be done exactly. We propose a rigorous approach which computes a multiple precision approximation on the filter's transfer function with a reliable error-bound. Then, we can check whether the implemented filter satisfies band specifications while taking into account the error of the transfer function re-computation, out of the actually used coefficients. Moreover, if the filter does not satisfy the band specifications, our algorithm computes a tight extension to the frequency bands that the filter does satisfy.

As a matter of course, filters do not only have to satisfy frequency specifications but also phase specifications. Our method might be amenable to this setting. We consider phase specifications to be beyond the scope of the paper, though.

The paper is organized as follows. Section II reminds some propositions about linear filters and defines the band specifications. In Section III, the reliable verification of a transfer function is exhibited. With the multiple precision approximation on the transfer function of any linear algorithm, Section IV checks whether an implemented filter satisfies its initial specification. Numerical examples that illustrate the method are given in Section V before conclusion.

Notation: throughout the article matrices are in uppercase boldface, vectors are in lowercase boldface, scalars are in lowercase. The complex unit is notated j (j 2 = -1). The conjugation operator on a complex number z is notated z .

II. PRE-REQUISITES

A. Digital filters and LTI systems

Digital filters are computational blocks that transform a signal u (i.e. a sequence u(k), where k ∈ Z stands for the step-time, according to the sampling period) into a signal y. In this article, we focus on discrete-time Linear Time Invariant (LTI) filters only, i.e. filters that are linear and for which a time shift (a delay) of the input sequence causes a corresponding shift in the output sequence [START_REF] Kailath | Linear Systems[END_REF], [START_REF] Oppenheim | Discrete-Time Signal Processing[END_REF]. They are standard basic bricks of digital filtering (and control) and include FIR (Finite Impulse Response) and IIR (Infinite Impulse Response) filters.

The output of a LTI filter can be expressed by

y(k) = n i=0 b i u(k -i) - n i=1 a i y(k -i), ∀k ∈ Z, (1)
where {a i } 1≤i≤n and {b i } 0≤i≤n are real numbers defining the filter (n is the order of the filter).

Instead of an input-output relationship in time (the output y(k) with respect to the previous and current inputs and previous outputs), such filters are more often described in the frequency-domain, by means of the Z-transform [START_REF] Oppenheim | Discrete-Time Signal Processing[END_REF]. For a given signal x, its Z-transform is defined as:

X(z) = +∞ k=-∞ x(k)z -k . ( 2
)
This transform is an equivalent of the Laplace transform but for discrete-time signals and systems. Applied to (1), we obtain a relationship between the Z-transform of the input and the output, denoted U and Y respectively:

Y (z) = n i=0 b i z -i U (z) - n i=1 a i z -i Y (z) (3)
and then Y (z) = H(z)U (z), with

H(z) n i=0 b i z -i 1 + n i=1 a i z -i . ( 4
)
H is called the transfer function of the filter, i.e. its input-output relationship in the z-domain. Its restriction to the unit circle ({z = e jω | ∀ω ∈ [0, 2π]}) is the Discrete-Time Fourier Transform (DTFT), so H(e jω ) gives the frequency response of the filter (if a sinusoid signal with normalized frequency ω 0 is applied as input of the filter, then the output will be a sinusoid with the same normalized frequency ω 0 , but amplified by H(e jω0 ) , and dephased by arg H(e jω0 ) ).

B. Filter specifications

In a classical signal processing flow, filters are designed from specifications in the frequency domain, in order to amplify (or preserve) signals in some frequency bands, and attenuate them in other bands.

A filter specification is then composed of several passbands (i.e. the gain of the filter for these frequencies should be bounded, often around 1) and stopbands (the gain should be lower than a given bound), formally described as:

β ≤ H(e jω ) ≤ β, ∀ω ∈ [ω 1 , ω 2 ]. (5)
The lower bound β is equal to 0 for stopbands, and usually β = 1δ and β = 1 + δ for passbands. For instance, Figure 1 exhibits a lowpass filter specification, mathematically described by

1 -δ p ≤ H(e jω ) ≤ 1 + δ p , ∀ω ∈ [0, ω p ] (passband) H(e jω ) ≤ δ s , ∀ω ∈ [ω s , π] (stopband)
For sake of generality, a filter specification will be described in this paper as a set of inequalities as in [START_REF] Roy | Basic algorithms in real algebraic geometry and their complexity: from Sturm's theorem to the existential theory of reals[END_REF]. This way of specifying filters is not the only conceivable way: instead of constant bounds β and β, one might consider upper and lower bounds be given as polynomials varying in the normalized frequency or even measures allowing spectral densities and partial violations of bounds to be taken into account. Similary, the filter specification only concerns here the magnitude, but not the phase. Considering these alternate ways of filter specifications shall be left to future work. We shall consider all frequency specification bounds as hard constraints. Our algorithm will return -in the first place-a boolean answer whether the specification is satisfied or not. Users of the algorithm might nevertheless decide to live with weaker bounds, allowing for margins that cancel out from a systems perspective. They may use our algorithm by adding their possible margins -possibly in a second run of the algorithm-to the specifications to allow for larger frequency specifications.

|H(e j! )| ! s ! p ! s 1 1 p 1 + p ⇡ passband stopband
Sometimes, filter designers prefer to give the bounds in decibels (x dB means 10

x 20 as a bound) and frequencies (normalized frequency ω and frequency f are linked by ω = 2π f Fs where F s is the sampling frequency of the filter). Due to Nyquist-Shannon theorem [START_REF] Oppenheim | Discrete-Time Signal Processing[END_REF], it is only necessary to specify the frequency up to Fs 2 (or the normalized frequency up to π instead of 2π). Due to lack of space, this article considers only filters with real coefficients, setting aside filters with complex coefficients.

III. VERIFYING BOUNDS ON A TRANSFER FUNCTION

The purpose of this Section is to detail our method that verifies that the modulus of a transfer function H stays between two bounds β and β for all z taken on a segment of the unit circle, corresponding to a certain frequency band, i.e. z = e jω for all ω ∈ Ω ⊆ [0, 2π]. In the case when the given bounds cannot be verified, our intention is to compute approximations to problematic frequencies for which the bounds are violated.

We proceed in three steps. In Section III-A, we show how we can reduce the given problem to showing that a rational function with real coefficients stays between two bounds for real arguments taken in a subset of [0, 1]. In Section III-B, we then further reduce the problem to showing that a polynomial stays non-negative over a subset of [0, 1]. In Section III-C, we briefly describe our approach to computing problematic frequencies in the case when the verification does not succeed.

A. Reducing the problem to a real rational function

We wish to verify that β ≤ |H(z)| ≤ β for all z = e jω with ω ∈ Ω ⊆ [0, 2π]. We suppose that H is given as a rational function H(z) = b(z) a(z) with real coefficients. Since we can suppose without lack of generality that β ≥ 0, this is equivalent to showing that

β 2 ≤ |H(z)| 2 ≤ β 2 , ∀z = e jω , ω ∈ Ω. (6)
Since z = e jω and the numerator and denominator polynomials a and b have real coefficients, conjugation of z yields z = 1/z and conjugation of the polynomials has no effect. So we have

|H(z)| 2 = b(z) b (z ) a(z) a (z ) = b(z) b(1/z) a(z) a(1/z) = v(z) w(z) , (7)
where v and w also are polynomials with real coefficients, obtained by simplifying the fraction b(z) b(1/z) a(z) a(1/z) . We have hence reduced the problem to verifying that

β 2 ≤ v(z) w(z) ≤ β 2 , ∀z = e jω , ω ∈ Ω. (8)
Taking now t = tan ω 2 , we can write z = e jω as

z = e jω = cos ω + j sin ω = 1 -t 2 1 + t 2 + j 2t 1 + t 2 . ( 9
)
By formally composing v and w with the expression z = 1-t 2 1+t 2 + j 2t 1+t 2 , for example by formal evaluation with Horner's scheme, and clearing numerators and denominators, we can hence obtain four polynomials r, s, ж, щ, all with real coefficients, such that

|H(z)| 2 = v(z) w(z) = v 1-t 2 1+t 2 + j 2t 1+t 2 w 1-t 2 1+t 2 + j 2t 1+t 2 = r(t) + j ж(t) s(t) + j щ(t) . ( 10
)
We can now observe that |H(z)| 2 is a real number and that the ratio r(t) s(t) is hence equal to the complex ratio r(t)+j ж(t) s(t)+j щ(t) . We may therefore drop ж and щ. We have now reduced the problem to verifying that

β 2 ≤ r(t) s(t) ≤ β 2 , ∀t = tan ω 2 , ω ∈ Ω ⊆ [0, 2π] , (11)
where the both polynomials r and s have real coefficients and all other quantities, β 2 , β 2 , ω and t are all real numbers. We must hence no longer deal with complex ratios and complex numbers and have reduced the problem to verifying the bounds of a real rational function over an interval, subset of the reals. Unfortunately, the mapping t = tan ω 2 maps the possible frequencies ω ∈ Ω ⊆ [0, 2π] onto to the whole real axis. In our experiments, we found this difficult to handle, partly because the tool we used, Sollya, has very little support for unbounded intervals and partly because having unbounded intervals meant searching for the zeros of certain functions over such unbounded intervals, which we found numerically unstable (see Section III-C for more details). We hence apply a second mapping: t = 1-2ξ ξ (1-ξ) . Still by formally composing the polynomials r and s with the expression t = 1-2ξ ξ (1-ξ) , we obtain two polynomials p and q with real coefficients such that

|H(z)| 2 = r(t) s(t) = p(ξ) q(ξ) . (12)
In the same step, we reduce the resulting rational function to its least terms to obtain p(ξ) q(ξ) . In order to do so, we extended the tool we used, Sollya, with an algorithm to compute the gcd of two polynomials [START_REF] Char | Gcdheu: Heuristic polynomial gcd algorithm based on integer gcd computation[END_REF].

As the inverse mapping ξ = t+2- √ t 2 +4 2t maps the reals onto the interval [0, 1], we have hence reduced our problem to verifying that

β 2 ≤ p(ξ) q(ξ) ≤ β 2 , ∀ξ ∈ Ξ ⊆ [0, 1] . (13)

B. Verifying the bounds of a rational function by showing the non-negativity of a polynomial

In order to verify an instance of (13), we can suppose that the interval Ξ the arguments ξ vary in is not reduced to a point and that β 2 = β 2 (otherwise a simple evaluation or a check whether p and q are constant polynomials suffice). We can hence reduce the problem further to obtain:

-1 ≤ 2p(ξ) -β 2 + β 2 q(ξ) β 2 -β 2 q(ξ) ≤ 1, ∀ξ ∈ Ξ. ( 14
) Let g(ξ) = 2p(ξ) -β 2 + β 2 q(ξ) and h(ξ) = β 2 -β 2 q(ξ).
It hence suffices to verify that

g(ξ) 2 h(ξ) 2 ≤ 1, ∀ξ ∈ Ξ (15)
which is equivalent to showing that

h(ξ) 2 -g(ξ) 2 ≥ 0, ∀ξ ∈ Ξ. ( 16
) Let f (ξ) = h(ξ) 2 -g(ξ) 2 .
Again f is a polynomial with real coefficients. We have reduced our problem to showing that the value of this polynomial f (ξ) stays non-negative over all ξ ∈ Ξ ⊆ [0, 1], where the interval Ξ is easily obtained from the original frequency domain

Ω = [ω 1 , ω 2 ].
Our approach to showing that f stays non-negative over Ξ is similar to the one set out in [START_REF] Chevillard | Efficient and accurate computation of upper bounds of approximation errors[END_REF]. We typically perform the following checks:

(i) Check whether f is positive at some (arbitrarily chosen) point ξ 1 ∈ Ξ by (interval arithmetic) evaluation of f at ξ 1 and that f has no zero over the whole interval Ξ. If so, f is non-negative over the whole interval Ξ. (i) Check whether f is positive at both endpoints of the interval Ξ by (interval) evaluation at these endpoints and that it has exactly one zero over whole interval Ξ, not counting multiplicities. The zero it has in the interval hence is of even multiplicity and the polynomial stays non-negative over the whole interval. (i) Check whether the interval Ξ can be split into subintervals such that one of the two aforementioned checks become satisfied. We test whether a polynomial (with real coefficients) has no, one or more zeros over an interval, bounded subset of the reals, utilizing Sturm's theorems on the Sturm sequence of the polynomial, similarly as done in [START_REF] Chevillard | Efficient and accurate computation of upper bounds of approximation errors[END_REF]. Sturm's theorem yields the number of real zeros of a real polynomial over a bounded interval, not counting multiplicities [START_REF] Roy | Basic algorithms in real algebraic geometry and their complexity: from Sturm's theorem to the existential theory of reals[END_REF]. The tool we used, Sollya, includes a fast but rigorous implementation of Sturm's technique [START_REF] Chevillard | Users' manual for the Sollya tool[END_REF].

C. Numerically computing problematic frequencies

In the case when our checks verifying if a given transfer function H(z) stays bounded in modulus by the two bounds β and β does not succeed, we numerically compute a list of problematic frequencies ωi , at which one of the bounds is violated. In contrast to the verification step which is completely rigorous in the sense that will never return a positive answer (i.e. the transfer function satisfies the given bounds) while the function actually does not, this numerical step is not fully rigorous. It may miss certain frequencies at which the bounds are violated. It is nevertheless pretty efficient with respect to speeding up the complete LTI filter verification algorithm we set out in Section IV-C, in particular concerning determining a reasonable verification margin (see below).

In our approach, we couple the verification process, described in Sections III-A and III-B, with the possibly needed step of computing problematic frequencies. These frequencies actually correspond, thru the different mappings ω → t → ξ, to points ξ at which the polynomial eventually obtained, f , takes negative values. We determine these points as the (negative) extremum points of f . We therefore differentiate f and compute approximations to the zeros of f by root isolation (still using Sturm's technique) and refinement with Newton-Raphson iterations. The tool we used, Sollya, offers all necessary basic bricks for these computations [START_REF] Chevillard | Users' manual for the Sollya tool[END_REF]. Once we obtain a list of points ξ i at which f becomes negative, we remap these values ξ i to a list of problematic frequencies, by following the inverse mappings ξ → t → ω.

¯ !

[ ] ] [ ] [ ] !1 !2 !3 ! 4
H(e j! )

Figure 2: If needed our algorithm returns the problematic frequencies as small intervals ωi .

IV. COMPLETE ALGORITHM FOR A LTI REALIZATION VERIFICATION

The above-described algorithm can be used to verify whether a transfer function satisfies a band specification. However, it is the quantization of filter's coefficients that can drastically impact the behavior of the implemented filter, leading to violation of the initial band specifications. In this section we propose an approach on verifying an implemented filter against user-given band specifications.

To apply our algorithm from Section III, we need to first recover the transfer function that corresponds to the implemented filter with quantized coefficients. Just like polynomials, a rational transfer function can be evaluated in many different ways. Numerous filter structures (filter algorithms) have been developed over time. With some exceptions (Direct Forms) filter structures do not directly use the coefficients of the transfer function but some modification of them. Deducing a transfer function from a structure differs from one structure to another [START_REF] Varga | Computation of transfer-function matrices of generalized state-space models[END_REF], [START_REF] Constantinides | Wordlength optimization for linear digital signal processing[END_REF]. Moreover, if they involve approximate computation, the error analysis must be done on a case-by-case basis.

In order to unify the approach of transfer function determination for any filter structure, we propose to use the Specialized Implicit Form (SIF) [START_REF] Hilaire | A unifying framework for finite wordlength realizations[END_REF]. In this paper we restrict ourselves to the Single-Input Single-Output (SISO) case, but an extension to the Multiple-Input Multiple-Output (MIMO) case is straightforward. In Section IV-A we present a quick overview of this representation. In Section IV-B we propose our approach on the computation of the filter's transfer function with the corresponding bound on the approximation error. Finally, in Section IV-C we propose an algorithm for the verification of any implemented filter against some band specifications.

A. Specialized Implicit Form

The Specialized Implicit Form (SIF) has been proposed as an extension of the state-space form:

  J 0 0 -K In x 0 -l 0 1     t(k + 1) x(k + 1) y(k)   =   0 M n 0 P q 0 r s     t(k) x(k) u(k)   (17)
where u(k) represents the input, and y(k) the output, x(k +1) is the n x stored states. Vector t(k +1) holds the n t intermediate variables in the calculations on step k. These calculations are topologically ordered in the matrix J , which is lower-triangular with 1s on the diagonal [START_REF] Hilaire | Reliable fixed-point implementation of linear data-flows[END_REF]. The diagonal matrix on the left side of the implicit equation ( 17) allows us to describe the sequence of computations within a filter. For example, y ← m 2 (M 1 x) is computed as t ← M 1 x and then y ← m 2 t. So this sequence is described as

I 0 -m 2 1 t y = M 1 0 x. (18)
A similar approach to (18) can be used to transform any linear signal processing algorithm or any input/output data-flow graph with delays, multiplications by constants and additions into a SIF. A systematic algorithm for the conversion of any linear data-flow graph represented in Simulink format to the SIF is presented in [START_REF] Hilaire | Reliable fixed-point implementation of linear data-flows[END_REF].

For example, a 3 rd order linear digital filter implemented with the ρ-Direct Form II transposed (ρDFIIt) algorithm [START_REF] Zhao | Roundoff noise analysis of two efficient digital filter structures[END_REF], [START_REF] Li | On the generalized DFIIt structure and its state-space realization in digital filter implementation[END_REF], will be described in Simulink as the data-flow graph given in Figure 3. We can observe, that this structure uses ten non-trivial coefficients. The corresponding SIF is the following 7 × 7 sparse matrix Z:

Z   -J M n K P q l r s   =    
where blue and pink rectangles are '-1' and '1', and circles are non-trivial coefficients.

Once rigorous filter analysis and implementation techniques are developed for the SIF formalism, they can be applied upon any linear realization (classical structures like direct forms, ρ-forms, lattice, etc. have already been directly converted to SIF [START_REF] Hilaire | A unifying framework for finite wordlength realizations[END_REF], [START_REF] Hilaire | Sensitivity-based pole and input-output errors of linear filters as indicators of the implementation deterioration in fixed-point context[END_REF], [START_REF] Volkova | Fixed-point implementation of lattice wave digital filter: Comparison and error analysis[END_REF]; others can be obtained from their input-output data flow graph [START_REF] Hilaire | Reliable fixed-point implementation of linear data-flows[END_REF]).

B. Reliable computation of the Transfer Function of a filter

The direct way to obtain the transfer function of a given SIF is to first convert it to a discrete-time state-space (dSS) representation S (A, b, c, d)

S x(k + 1) = Ax(k) + bu(k) y(k) = cx(k) + du(k) (19)
A = KJ -1 M + P , b = KJ -1 n + q (20a) c = lJ -1 M + R, d = lJ -1 n + s. ( 20b
)
The coefficients A, b, c and d can be computed exactly from a SIF, since matrix J is lower-triangular with 1s on the main diagonal.

Then, applying the Z-transform of ( 2) to [START_REF] Gevers | Parametrizations in Control, Estimation and Filtering Probems[END_REF] we have the classical formula for the transfer function of a state-space [START_REF] Oppenheim | Discrete-Time Signal Processing[END_REF]:

H(z) = c(zI -A) -1 b + d. ( 21
)
By considering the eigendecomposition V EV -1 of the matrix A (with {λ i } 1≤i≤nx its associated eigenvalues, i.e. the diagonal elements of E), we have

H(z) = c(zI -V EV -1 ) -1 b + d (22) = cV (zI -E) -1 V -1 b + d (23) = cV    1 z-λ1 . . . 1 z-λn    V -1 b + d, (24)
Therefore, the transfer function can be written as a rational function

H(z) = b(z) a(z) with a(z) = nx j=1 (z -λ j ) (25) b(z) = nx i=1 (cV ) i (V -1 b) i j =i (z -λ j ) + d a(z) (26)
Then, the accuracy of the computation of the transfer function relies on the accuracy of the computation of its eigendecomposition and on the computation of the polynomial coefficients. Exhibiting a bound on the accuracy of the eigendecomposition is a non-trivial task. Moreover, as it will be shown below, we may need to iteratively increase the accuracy of the computed transfer function. Thus, we propose to use Multiple Precision Floating-Point arithmetic to compute an approximation H(z) using formulas ( 25) and ( 26).

Suppose we have available a multiple precision eigensolver such that its error monotonically decreases when its working precision increases, and multiple precision basic bricks for matrix arithmetic. Then, increasing the precision of the computation of H decreases its error H -H for any transfer function norm . . In order to bound this absolute error we propose the following:

Step 1: since we cannot exactly compute H (at least not for any structure), we first compute the approximation H on the transfer function of the dSS S = (A, b, c, d) using an eigendecomposition and Multiple Precision arithmetic (equations ( 25) and ( 26));

Step 2: we compute the system S which exactly corresponds to the approximation H(z) using the controllable canonical state-space [START_REF] Oppenheim | Discrete-Time Signal Processing[END_REF]: if { b i } 0≤i≤n and { a i } 1≤i≤n are the coefficients of the polynomial defining H, then S is given by ( A, b, c, d) with:

A =       -a 1 1 . . . . . . . . . 1 -a n 0 . . . 0       , b =       1 -a 1 b 0 . . . . . . b n -a n b 0       c = 1 0 • • • 0 , d = b 0 (27)
Step 3: we exactly compute the state-space difference ∆S = S -S, which is the difference between outputs of the two state-spaces, using the classical formula [START_REF] Kailath | Linear Systems[END_REF]:

∆A = A 0 0 A , ∆b = b b (28a) ∆c = c -c , ∆d = d -d (28b)
Step 4: then, we can compute a bound on the approximation error H -H. By definition of the state-space ∆S, we have H = H + ∆H where ∆H is the transfer function of ∆S. Using the relations between the transfer function and impulse response norms of a filter [START_REF] Boyd | Comparison of peak and rms gains for discrete-time systems[END_REF], we have ∀ω ∈ [0, 2π], ∆H(e jω ) ≤ sup ω ∆H(e jω ) ≤ ∆h 1 where ∆h(k) is the impulse response of the system (response to the filter ∆S with transfer function ∆H to an impulse signal), and ∆h 1 its 1 -norm, i.e. ∞ k=0 |∆h(k)|. As shown in [START_REF] Balakrishnan | On computing the worst-case peak gain of linear systems[END_REF], [START_REF] Volkova | Reliable evaluation of the worst-case peak gain matrix in multiple precision[END_REF], the norm ∆h 1 can also be computed using the Worst-Case Peak-Gain (WCPG) measure of ∆S, denoted ∆S , and defined by:

∆S ∞ k=0 ∆c(∆A) k ∆b + |∆d| . ( 29
)
Using the result given in [START_REF] Volkova | Reliable evaluation of the worst-case peak gain matrix in multiple precision[END_REF], we can compute ∆dSS with arbitrary precision. Finally, the approximation error H(z) -H(z) is bounded, for any z in the unit circle:

∀ω ∈ [0, 2π], H(e jω ) -H(e jω ) ≤ Θ
where Θ = ∆S + ε is the Worst-Case Peak Gain of the system ∆S computed with arbitrary small absolute error bounded by ε > 0.

C. Reliable verification of the implemented filter

Now, given an implemented filter, we represent it in SIF and then compute, with some precision, an approximation H on its transfer function along with a rigorous upper bound Θ on the corresponding ∆S . We choose the precision of the computations based on an heuristic which increases the compute precision in a loop and is, hence, reasonably fast and accurate but does not provide any guarantee that the precision will eventually be enough. The result will always be reliable, though.

Then, to verify whether the exact magnitude response evaluated on the unit circle is in the interval If the verification is not successful, we obtain a list of problematic frequencies with the algorithm from Section III. Thus, we can compute the maximum excess of the bounds, enlarge the band by this amount and repeat the process. For example, on Figure 6i we suppose that the approximation H(e jω ) is too close to the bound or Θ is too large. Then, we enlarge the band margin up to some new bound β while decreasing Θ, and perform the verification again with the updated band. This algorithm may yield a false negative answer if the initial precision of the transfer function computation was not large enough, which we never observed on the numerical examples we ran the algorithm on.

|H(e j! )| ! ¯ ⇥ ¯ (i) |H(e j! )| ! ¯ 0 ¯ 0 ⇥ ¯ (ii)
The detailed description of the heuristics that we used can be found in the long version of the paper that will be published along with our software implementation.

V. NUMERICAL EXAMPLES

The algorithm presented in Section III was implemented using a modified Sollya tool 3 and further bound with our automatic filter generator using pythonSollya 4 . We use an implementation of the algorithm for the WCPG5 computation in arbitrary precision [START_REF] Volkova | Reliable evaluation of the worst-case peak gain matrix in multiple precision[END_REF], which is written in C, using GNU MPFR version 3.1.12, GNU MPFI version 1.5.1 and CLAPACK version 3.2.1. Experiments were done on a laptop computer with an Intel Core i5 processor running at 2.8 GHz and 16 GB of RAM. We present three different use cases of our algorithm illustrated with examples.

Example 1: Our tool can be used to certify an already existing filter implementation. Consider the implementation represented as the Simulink data-flow graph of Figure 7, with gains given as 8-bit constants g 1 = 0.34765625, g 2 = 0.3359375 and g 3 = 0.08496094. The task is to verify whether the given implementation is a lowpass filter which satisfies following normalized frequency constraints: passband ω p = 0.1 with passband amplitude between 1 dB and 3 dB, stopband starting ω s = 0.3 with minimum attenuation 20 dB. First, we convert the Simulink graph to the SIF (for which we can compute the transfer function with arbitrary precision). Then, we apply our algorithm from Section IV-C and obtain a positive result in 1.9s. Thus, we obtain a guarantee that the given filter implementation satisfies the desired frequency response requirements. Remark Figure 7 represents a Lattice Wave Digital Filter, the coefficients of which are usually derived from the specifications using direct formulas without actually computing the filter's transfer function [START_REF] Gazsi | Explicit formulas for lattice wave digital filters[END_REF].

g1 -1 -1 g2 -1 -1 -1 g3 -1 -1 1 u(k) 1 y(k) Z -1 Z -1 Z -1 0.5
Example 2: Suppose that we need to choose among four different filter realizations: a Direct Form II transposed, a Direct Form II transposed with optimized ρ operator [START_REF] Li | On the generalized DFIIt structure and its state-space realization in digital filter implementation[END_REF], a balanced State-Space [START_REF] Gevers | Parametrizations in Control, Estimation and Filtering Probems[END_REF] and a Lattice Wave Digital Filter [START_REF] Gazsi | Explicit formulas for lattice wave digital filters[END_REF]. These realizations have different number of coefficients and were designed using different approaches that are beyond the scope of our article. Their floating-point coefficients are results of various approximations and optimizations specific for each structure and can be found in the Appendix.

Our goal is to verify whether after quantization of the coefficients to different formats these realizations will satisfy following frequency response constraints: sampling frequency F s = 48 kHz, passband up to 2.4 kHz with amplitude in [-0.5dB; 0.5dB]; stopband starting 7.2 kHz with minimum attenuation 80 dB.

The results of verification of the above realizations to 32, 16 and 8 bits are listed in Table I. We can observe that for certain filters a specification is fulfilled only if a safe margin is added; for certain cases that margin is as large as 6.71e -1 dB. Exceptions similar to the balanced State-Space and Lattice Wave structures may also occur: with 32 bits coefficients, their transfer functions do not pass the verification but the coefficient quantization effects cancel out for 16 bits and verification "luckily" passes.

If 16-bit quantization is the target format, we see that except DFIIt, all structures verify the specifications and now the designer may concentrate on choosing the best realization according to other criteria (for example, number of coefficients or datapath delay). If 8-bit quantization is the target format, the safe margins computed with our tool can be used to redesign the realizations and repeat the verification.

On the other hand, if the realizations are impossible to redesign, we can determine what is the maximum quantization of the coefficients for which our tool gives a positive answer.

Finally, for the cases when the transfer function respects the band specifications our algorithm gives an answer quickly. Most time is spent on computing the safe margin, especially when an overflow is on the edge of the band. Example 3: Finally, our algorithm can be used to certify the result of different design methods even before structure choice or quantization. For example, we may compare the quality of different transfer function design methods in Matlab. Some of the widely used algorithms for the transfer function design are Butterworth [START_REF] Butterworth | On the theory of filter amplifiers[END_REF], Elliptic [START_REF] Zolotarev | Application of elliptic functions to problems about functions with least and greatest deviation from zero[END_REF], [START_REF] Cauer | Synthesis of linear communication networks[END_REF] and Chebyshev [START_REF] Daniels | Approximation Methods for Electronic Filter Design[END_REF].

We propose to verify whether transfer functions (of order suggested to be minimal by Matlab) satisfy their initial frequency constraints. They always should. For a simple test, we consider the frequency specifications used in Example 2. To take into account that the computations in Matlab are done in double precision, we find to reasonable to apply a design margin 10 -15 dB before using our algorithm from Section III. From Table II we see that only Butterworth design method instantly satisfies the specifications. Despite the initial design margin, Chebyshev and Elliptic methods still may require up to 5.65e-12 dB of margin to verify the desired specifications. On the other hand, even these margins may be acceptable depending on filter designer's needs.

Thus, our algorithm can be used to certify that the filter design method is rigorous, or give the designer a perception of the sufficient correction of the design margin. Lattice Wave Digital Filter: This realization has the smallest number of coefficients that is equal to the order of the filter. This structure consists of two parallel branches (top and bottom) that realize all-pass filters. Coefficients γ t of the top branch of the realization:

Figure 1 :

 1 Figure 1: A lowpass filter specification.

Figure 3 :

 3 Figure 3: Simulink data-flow of a ρDFIIt form.

Figure 4 Figure 4 :

 44 Figure 4: Dotted line: inexact transformation. Straight line: exact transformation.

Figure 5 :

 5 Figure 5: Verify whether H(e jω ) in [β + Θ; β -Θ].

Figure 6 :

 6 Figure 6: If verification fails (i), we enlarge the bound (ii).

Figure 7 :

 7 Figure 7: Implemented digital filter.

  γ t = (0.309545756846, 0.426986744076, 0.0648857726338, 0.115929866936, 0.0648857724342) (32) Coefficients γ b of the bottom branch of the realization: γ b = (0.499598692531, 0.0648857723853, 0.30101467827, 0.0648857724933) (33)

Table I :

 I Checking quantized realizations.

		wordlength	32	16	8
	DFIIt	margin time	12.49s	unstable -	unstable -
	ρ DFIIt	margin time	13.12s	4.19s	4.68e-3 dB 104.01s
	State-Space	margin	6.16e-10 dB		6.71e-1 dB
	Balanced	time	12.27s	18.18s	92.05s
	Lattice Wave	margin time	3.80e-10 dB 920.88s	4.58s	1.73e-2 dB 200.83s

http://www.mathworks.com/

https://www.scipy.org/

https://scm.gforge.inria.fr/anonscm/git/sollya/sollya.git

https://scm.gforge.inria.fr/anonscm/git/metalibm/pythonsollya.git

https://scm.gforge.inria.fr/anonscm/git/metalibm/wcpg.git

   8.09616443886441e -8 7.28654809378781e -7 2.91461913803914e -6 6.80077821613168e -6 1.02011670932711e -5 1.02011672709068e -5 6.80077808468127e -6 2.91461921175795e -6 7.28654797166328e -7 8.09616444719108e -8  -1 1 0 0 0 0 0 0 0 0 8.09616444e -8 6.91501406e + 0 0 1 0 0 0 0 0 0 0 7.28654799e -7 -2.14466416e + 1 0 0 1 0 0 0 0 0 0 2.91461920e -6

3.91222980e + 1 0 0 0 1 0 0 0 0 0 6.80077813e -6 -4.62252163e + 1 0 0 0 0 1 0 0 0 0 1.02011672e -5

3.66656397e + 1 0 0 0 0 0 1 0 0 0 1.02011672e -5 -1.95139227e + 1 0 0 0 0 0 0 1 0 0 6.80077813e -6 6.71659680e + 0 0 0 0 0 0 0 0 1 0 2.91461920e -6 -1.35615063e + 0 0 0 0 0 0 0 0 0 1 7.28654799e -7

1.22341213e -1 0 0 0 0 0 0 0 0 0 8.09616444e -

9.98078553860207e -4 2.15657523951793e -4 2.76696614924035e -5 3.32572318297719e -3 6.88163327582342e -4 8.89436597243942e -5 7.35371936546042e -3 1.58919858005338e -3 2.03894253732018e -4 -1.80048266605967e -2 -3.65628640741023e -3 -4.74107807447374e -4 4.36139632080109e -2 1.01333758827138e -2 1.28486543968479e -3 -2.02071915732453e -1 -3.20997699114184e -2 -4.34923676306156e -3 6.5167003556539e -1 -1.75022739029844e -1 -1.86034616397461e -2 1.75022739029834e -1 5.96527629913795e -1 -1.29541395028042e -1 -1.86034616397113e -2 1.2954139502792e -1 5.43595944182089e -1

Butterworth

Chebyshev Elliptic Matlab OK/margin 5.65e-12 dB 5.65e-12 dB time 0.81s 11.47s 2.36s

Table II: Verification of transfer function design methods.

VI. CONCLUSION AND PERSPECTIVES

With this paper, a rigorous method to verify a transfer function of LTI digital filters against band specifications in the frequency domain has been developed. It relies on translating the problem of verifying bounds on magnitude response evaluated on a unit circle to the verification of positivity of a real polynomial. Our algorithm guarantees that no false positives occur. In the case of unsuccessful verification a list of problematic frequencies is provided, for which we compute the maximum overflow over the band specification. We propose an implementation using a combination of interval and rational arithmetic in Sollya tool.

We applied this method to develop an approach on the verification of any implemented filter. It relies on the multiple precision computation of a transfer function corresponding to the filter realization with quantized coefficients. Then, we bound the approximation error using the Worst-Case Peak Gain measure in arbitrary precision.

This approach opens various possibilities on the verification of digital filters. As it is automatic, our approach allows for easy integration into automatic filter implementation tools.

Moreover, our verification algorithm can be applied to compare different filter structures with various Fixed-Point settings. Such a comparison offers a filter designer an overview of the implementation possibilities. On top of that, the information on the safe margin can be used as a margin in design process. For instance, we can narrow down the initial band by the safe margin and re-design the filter with the new band specifications. Then, the re-computed filter with quantized coefficients should respect the initial conditions. Finally, using the list of problematic frequencies we can probably improve the rational Remez [START_REF] Selesnick | Magnitude squared design of recursive filters with the chebyshev norm using a constrained rational remez algorithm[END_REF] algorithm, that is usually behind the transfer function design.

However, overall time-efficiency of our implementation can be improved. Passing too much time on the computation of the safe margin can be a significant drawback were the algorithm to be used during the exploration of a large design space.

The approach sets out this paper is concerned with the static quantization of the filters coefficients only. The dynamic roundoff error observed when the implemented filter runs is taken into account by other approaches [START_REF] Menard | Analytical fixed-point accuracy evaluation in linear time-invariant systems[END_REF], [START_REF] Volkova | Determining fixed-point formats for a digital filter implementation using the worst-case peak gain measure[END_REF], [START_REF] Constantinides | Roundoff-noise shaping in filter design[END_REF]. Combining the analysis of both error sources is future work.

APPENDIX

A. Coefficients of filters from Example 2

In Example 2 from Section V we consider four realizations that realize a 9 th order IIR filter. Their coefficients have been subject to different optimizations specific to each structure.

Direct Form II transposed: This realization has the same coefficients as the transfer function of the filter. Its numerator b and denumerator a are:


Mathilde Plard

Body in motion: experience of walking amongst elders immigrant women in Montréal

As part of an ongoing study on "older immigrant women and their relationship to aging" conducted in Montreal, the experiences of aging in regards to the self and body are particularly questioned. The data presented in this paper are from 15 focus groups, 90 women from different backgrounds and cultural communities with diverse migration routes. From this data, the lives of these immigrant women are presented through a common practice: walking. A comparative analysis of transcripts especially has helped to understand the importance granted to movement. Notably walking attests to the unanimous commitment these women reveal; that of remaining autonomous and independent. This is the result found in all focus groups. Indeed, whatever the person's age, autonomy, cultural background, religious practices, their migratory route or family status, walking is a way to experience and live their environment -their lived space. Beyond the vision that describes aging in terms of limitations and disabilities, we mobilize the prism of movement to show how these women adjust their body to stay as much as possible consistent with their environment. We highlight in this section if getting older and having physical limitations may represent obstacles in every day movement, however not necessarily in a definitive way. Furthermore, we will show what this attitude towards walking reveals about the person practicing it. Travel and daily movements, for example, enables social relations.

INTRODUCTION

Relational dynamics of the subject: "to the self, to others and to the environment" This paper presents a discussion based on data from a research that explores the experience of aging in older immigrant women in Montreal. It aims to show how the practice of walking daily is an activity in which subjects encountered build relationships at different levels. These connections are illustrated in the figure (see opposite) and link the "subject" to the self and body (B), to others (O) and to the environment (E). The articulation of these three relational dynamics stemming from the subject also supports the theoretical notion of plasticity of these relationships in order to reach a balance. Thus, the socioplasty evoked allows for the sharing of the subject's creative dimensions that mobilizes different relationships according to biographical experiences. In this way, the experiences of aging of older immigrant women met in Montreal reveal how these relationships are arranged. Walking is also a relevant entry to grasp how a spatial practice is a tool of expression and identity construction [START_REF] Kusenbach | Street Phenomenology: The Go-Along as Ethnographic Research Tool[END_REF].

The framework of the study offers, firstly, general elements in order to situate the results (1). The results are then presented in the form of relational typologies in order to understand the meaning behind walking for women interviewed through the different relationships expressed in regards to their body, to others and to the environment (2). The discussion reflects on how the practice of walking ultimately shapes their meaning. This last part suggests that these "relationships" stem from the "meaning" (3). The conclusion addresses three important steps in the process of relationship building and opens onto a general notion of plasticity.

Context of the study

This is a presentation of some of the results of a qualitative research aimed to analyze the personal, social and family-oriented experiences of aging. This research1 was conducted in Montreal, Quebec (Canada) with 83 women aged 65 years and older. Here, I present some of the findings concerning migration patterns of women and the changes and adjustments that immigration brings to their lives.

An accelerated aging of the population of Quebec

In 2011, 15.9% of Quebec's population is 65 years old and over.

(Gouvernement du Québec, 2012. Les aînés du Québec. Quelques données récentes)

Methodology

Qualitative research through interviews 18 focus groups each comprising 3-6 elderly women (n = 83) of the same ethno-cultural background and speaking the same language, but having diverse characteristics such as marital status, education level, income and migration routes. Individual interviews (ongoing analysis) (n = 15)

Sample

Our sample included elderly immigrant women with diverse variables such as: age, ethnocultural background (Arab, African, Haitian, Japanese, Chinese, Portuguese, Romanian, etc.); their migration routes or timing of migration (during youth, adulthood or retirement) and immigration status. We formed 18 focus groups; each consisted of 3-6 women sharing the same ethno-cultural background and speaking the same language, but have varied characteristics in terms of marital status, education, income and migration patterns. The 83 elders which we will present are from three continents and a variety of countries: Ex-Yugoslavia, Bosnia, Serbia, Croatia (10), Romania (6), Mexico (2), Guatemala (1), El Salvador (3), 5 Colombia (5), Algeria (3), Egypt (6), Lebanon (4), China (9), Japan (9), Congo (5) , Haiti (10), Jamaica (3) and Portugal (7).

RESULTS. "Walking in order to shape relationships"

It is possible to structure into three categories the types of relationships that walking establishes. For most of the women interviewed, discussing and defining their experiences of aging implies identifying what reduces or increases that capacity. What increases or decreases physical and functional capacities (2.1.), what increases or decreases cognitive, social and interpersonal skills (2.2.) or what increases or decreases the capacity to experience their own living space (2.3.).

2.1.

Walking in order to feel the relationships towards the self, towards the body It seems as though aging is first felt, physically, through the body. The analysis of transcripts shows how women relate the experience of aging to their own changing body. For them it is "feeling their body", to set themselves "in motion" and "to be mobile and independent". These elements reveal the anguish of the dependence particularly present in Western societies. In a social and cultural context in which the experience of aging is associated with a vision articulating dependence VS autonomy, it is conceivable to ask how the subjects encountered construct themselves. The practice of walking daily is an activity that first reveals the experience of aging by physically making the subject feel that the "body is changing" and then presented as a tool to act upon the self and "keep fit".

Feeling your body age: sensations of an aging body: "Until I was 65, I did not feel like an old person inside, I felt like a young person. But at 65, I was operated on my knees and since that time, I started... It was like going down the stairs, we stop at every step, then decrease, and decrease" (Portuguese, p.25). "Personally, I started to feel that I was getting old when I had back pains. When I had back pains I could not walk so I said, "it starts" (Congolese, p.11).

Walking, an activity to stay fit and a normative imperative: "Even if my health is not good I walk anyway. Today I moved. Someone called me at home and told me: "C. you should not stay at home", "No, God talks to me, he told me not to stay at home" I have to walk; I have to move, move a lot" (Haitian, p.10).

Walking towards others

Walking as a meeting tool: destinations of choice in order to get to others. The common itinerary usage of "to" or "by" certain places reveals clearly structured strategies that develop social ties (church, coffee, etc.). Destinations where women meet friend(s), relatives, acquaintances, etc.: "Then, I want to go out; I cannot walk much because my legs do not let me. But I like to go to Mass, I always get to talk to some ladies, they also have the same problems as me, when I get home, I am relieved but the days are long ... "(Portuguese, p.32). Walking with others: "We walk together. I always walk with my husband "(Ex-Yugoslavian, p.9).

Walking to identify and develop a sense of space

Just as it is possible to connect the experience of aging with the body, the relationship to space runs parallel to this dialectic. Therefore, we must grasp what limits or increases women's practices in their local area -direct geographical environment, lived space. Specific habitual routes, neighborhood practices, strategies in relation to the weather, etc.: In regards to snow, "it is slippery, I'm afraid of falling and hurting my wrist again. I have an arm that has no strength, if I fall and hurt my other arm, I could do nothing" (Chinese, p.39). In regards to walking for one's own self, one can note the importance of temporality and spatiality. Indeed, there is seasonality in the practice of this activity. It is not the same in the snow in winter than it is in summer. The itineraries vary in shape and duration. The discussion initiated by Bernard Duperrein is relevant to this discussion as it brings up the notion of how it is possible to "make sense of space" [START_REF] Duperrein | La maison de retraite, dernière installation[END_REF].

DISCUSSION "bonds born through meaning"

3.1. Walking activates processes of plural domestication Indeed, "the domestication of the world is Man's primary concern in society. The purpose of this process is none other than the production of meanings by humans to the space they are in. In this way, individuals domesticating their environment establish the ontological foundations of their world, making it livable and practical. The meaning given and produced provides points of references (orthonormal) necessary for individual practices. The methods for making sense of space and time are also directly related to the social, political, economic and cultural environment in which each individual evolves. Also, forms of domestication results from both the subject that produces them, and the situation in which the individual is a part of" [START_REF] Plard | Les lieux du vieillir[END_REF].

Walking and the experience of aging: daily updates of the self

Objective: To remain consistent with one's environment in terms of one's physical, functional, psychological capacities, etc. In other words, it is to adjust the terms of interactions in terms of one's own aging experience. The personal and intimate relationship "of" and "in" the aging body expresses the dialectic between body and identity. In this way, we find the dynamics that bind the women interviewed in both the process of change and continuity in everyday life. Indeed, we observed interactions between transformations /variations (important transitioning elements towards the end of life) and constancy /stability /sustainability of self-identity (which is not simply what the signs of time passing leave as traces on the body nor what society -and others in general -determine as stigmas.

In other words, it is to find a balance in the constantly changing body as expressed by this Chinese woman in regards to her relationship to her appearance. "Every day, there is a slight difference. It is worse than before. Our 'sifu' (Buddhist masters) that tell us about the lessons given by the Buddha told us that our appearance changes every day, but we do not see it "(Chinese, p.30). The aim is to deal with one's aging body without forgetting its identity, the heart that does not change -or "narrative identity" according to Paul Ricoeur. "I'm not old! [Speaking loudly] I am still very young! Especially do not say I'm old [laughing]. The person is old, but the heart is not old [expression: meaning that the heart remains young despite appearances] "(Chinese, p.30).

Walking, a metaphysical activity: translation of a vital impulse

Practice to embody, feel and live "Aging and walking: counting on movement"

Walking -vitality (practice of one's lived space, "moving forward", "If I could walk, it would be ideal, I feel as though I have a creative power in me!" (Romanian p.2). Continuing walking while recognizing the aging process: "I like to walk, I feel like I have a lot of vitality. I walk every day. I do not use transportation because I like to walk. Yes, I look in the mirror and I know I'm older, I realize that my body is not the same every year ... "(Colombian, p.5). "Aging is like making a stop during the long race we call youth. Aging is like stopping a little bit and starting to walk more slowly in the path of life "(Colombian, p.2).

CONCLUSION: Walking, a tool to create and rebuild relationships

In order to further this presentation, there should be complementary studies that explore more in detail the evolution of the relationship that women have with themselves, others and their environment. It would be possible to explore how spatial practices of the environment and walking in particular change over time. What were their spatial and social practices at the time of their arrival in Montreal? Were these practices more developed at that time because of their young age and broader functional capabilities? Or on the contrary, is their lived space, as of today-while aging-more developed, better built and more organized than it was a decade ago since their functional capacities are now reduced due to their advanced age? Does knowledge of the environment and the relationships built during a lifetime help stabilize a lived space set in motion by international migration? In view of life events experienced, trajectories and life histories, subjects would develop and accumulate a capacity to cope -coping capacities [START_REF] Schröder-Butterfill | A framework for understanding oldage vulnerabilities[END_REF] variables. This idea is also supported by a dynamic vision of biographical trajectories developed by Andrew Delano Abbott in the book Time matters: on theory and method [START_REF] Abbott | Time matters : on theory and method[END_REF]. On the figure below you can read the different steps through which people who have migrated pass in order to "domesticate" a new environment, Montreal, in our study. Walking is one tool amongst others for building relationships that migration has initially broken.

Research on the development of partnerships funded by the Social Science Research Council of Canada (Charpentier

[2011][2012][2013][2014].


Guillaume Damiand
email: damiand@lirmm.fr

Patrick Resch
email: resch@lirmm.fr

Topological Map Based Algorithms for 3D Image Segmentation

One of the most commonly used approach to segment a 2D image is the split and merge approach. In this paper, we are defining these two operations in 3D within the topological maps framework. This mathematic model of regions segmented image representation allows us to define these algorithms in a local and generic way. Moreover, we are defining a new operation, the corefining, which allows to treat big images. They are cut into small units, treated separately, then the result of each of them are combined to reconstruct the final representation. These three operations let us view efficient 3D segmentation algorithms, which is a difficult problem due to the size of data to treat.

Introduction

The region segmentation is a difficult problem that was studied in many different works in 2 dimensions. It consists in making a partition of an image into connected sets of pixels verifying an homogeneity criterion, and that we are calling regions. The main approaches for the region segmentation are the split and merge methods.

The top-down approach [START_REF] Ohlander | Picture segmentation using a recursive region splitting method[END_REF][START_REF] Lee | Recursive region splitting at hierarchical scope views[END_REF] consists in taking big regions and cutting into smaller and smaller regions. The bottom-up approach [START_REF] Brice | Scene analysis using regions[END_REF][START_REF] Fiorio | Two linear time union-find strategies for image processing[END_REF] is the opposite approach, which begins with many small regions that are merged into bigger and bigger regions. At last the mixed approach [START_REF] Horowitz | Picture segmentation by a directed split-andmerge procedure[END_REF][START_REF] Pietikainen | Split and link algorithms for image segmentation[END_REF] consists in mixing the two previous ones, possibly making again the process until for example, result stabilization.

These approaches require a "good" model of images representation. Many works in 2 dimensions [START_REF] Domenger | Conception et implémentation du noyeau graphique d'un environnement 2D1/2 d'édition d'images discrètes[END_REF][START_REF] Brun | Segmentation d'images couleur à base Topologique[END_REF][START_REF] Brun | A new split and merge algorithm with topological maps and inter-pixel boundaries[END_REF][START_REF] Braquelaire | Image segmentation with topological maps and inter-pixel representation[END_REF] have shown that topological maps are a model that allows to introduce these segmentation algorithms in an efficient way. Moreover, thanks to that model we can also define many processing algorithms allowing to access or to modify the result of this segmentation.

Recently, the topological maps have been extended in dimension 3 [START_REF] Braquelaire | A topological structuring for aggregates of 3d discrete objects[END_REF][START_REF] Bertrand | Topological map: Minimal encoding of 3d segmented images[END_REF]. Indeed, more and more domains need to work in dimension 3, as medical imagery, geology, or industry. In order to carry out segmentation algorithms as in dimension 2, we need to define basic operations on topological maps. This problem turns out to be more difficult than in dimension 2, because first now we have an additional dimension, that sets down new problems, but also because it is much more difficult to represent and to make visual objects in 3 dimensions. These visualization problems act as a brake to the comprehension and the development of new algorithms. Moreover, there are complexity constraints, in memory space as well as in execution time, that are much more important than in dimension 2. Indeed, the treated data quantity is much more important in dimension 3. It requires us to define very efficient algorithms to use them on big images.

In this paper, we present the two algorithms of merge and split on the 3D topological maps. These two algorithms are the basic operations for the regions segmentation. The use of the topological maps let us define these algorithms in a generic way, because they work on any configuration, and in a local way because they treat the map element by element, only looking the direct neighborhood of the current element. These two properties make these algorithms more simple to understand but more efficient in complexity too. We also present the algorithm of corefining that allows us to treat big 3D images in parallel. We are going to cut this image into several small units, we are going to segment each unit in parallel, then we will reconstruct all the image with this units, using the corefining operation. Because of a lack of space and to not lose ourselves into technical details we are just presenting the principle of the algorithms and the main ideas of these three operations. For more details about these operations, we can report to [START_REF] Resch | Algorithmes pour la manipulation des cartes topologiques en 2 et 3 dimensions[END_REF][START_REF] Resch | Algorithmes pour la manipulation des cartes topologiques en 2 et 3 dimensions[END_REF].

We are first presenting Section 2 the combinatorial maps then the topological maps that are combinatorial maps verifying specific properties. Then we are presenting our three algorithms : merge Section 3, split Section 4 and corefining Section 5. We are describing their principle and the different cases we have met. At last, we are concluding and presenting some perspectives Section 6.

Topological Maps Recall

Topological maps allow to represent the nD regions segmented images. They encode at the same time the topology and the geometry of images. Topological maps are combinatorial maps with particular properties. So we are beginning by recalling the notion of combinatorial map. This is just a short reminder; a more detailed description can be found in [START_REF] Braquelaire | A topological structuring for aggregates of 3d discrete objects[END_REF][START_REF] Bertrand | Topological encoding of 3d segmented images[END_REF].

Combinatorial Maps

Combinatorial maps are a mathematical model of representation of space subdivisions in any dimension. They were introduced in the sixties by [START_REF] Cori | Un code pour les graphes planaires et ses applications[END_REF], at first as a planar graph representation model, and extended by [START_REF] Lienhardt | Subdivision of n-dimensional spaces and n-dimensional generalized maps[END_REF] in dimension n to represent orientable or not-orientable quasi-manifold. Combinatorial maps encode space subdivisions and all the incidency relations. They are made of abstract elements, called darts, on which are defined application, called β i . We are giving here the combinatorial map definition in n dimensions, that we can find for example in [START_REF] Lienhardt | Topological models for boundary representation: a comparison with n-dimensional generalized maps[END_REF]. Definition 1 (combinatorial maps). Let n ≥ 0. A n combinatorial map, (or n-map) is an (n + 1)-uplet M = (B, β 1 , . . . , β n ) where :

1. B is a finite set of darts; 2. β 1 is a permutation on B; 3. ∀i, 2 ≤ i ≤ n, β i is an involution on B; 4. ∀i, 1 ≤ i ≤ n -2, ∀j, i + 2 ≤ j ≤ n, β i • β j is an involution.
In this definition, there is an application β i for each space dimension which puts in relation two i-dimensional cells. When two darts are linked with β i , we say that they are β i -sewed. Each space cell is implicitly represented by a set of darts. We can see figure 1.a an example of an image and figure 1.b the corresponding combinatorial map. Each dart is represented by a segment, the β 1 relation by light grey arrows and the β 2 relation by dark grey arrows. β 1 put in relation a dart and the next dart of the same face. For example, the light grey face of the image is represented by four β 1 -sewed darts in the map. The adjacency between this face and the dark grey face is represented by two darts β 2 -sewed together. We are using the simplified representation (figure 1.c) that does not represent explicitly the applications, because it is more understandable.

The Topological Maps

In the combinatorial maps framework, several representations of a single object exist. We want a unique characterization of objects, for example to make easier isomorphisms. That is the main goal of the topological maps. They are mathematical model of 3D segmented images representation, which encode all incidency and adjacency relations. They represent interpixel elements composing the edge of boundary faces of an 3D image. Moreover, they are minimal, stable for rigid transformations, and they characterize the image's objects with their topology. We remind that a boundary face is a surface between two neighbouring regions. We can see on figure 2.a an example of a 3D image, and on figure 2.b its boundary faces. The construction of a topological map is progressive : at the beginning all image's voxels are encoded with a combinatorial map, then we simplify this map with successive mergings. The corresponding topological map of our example is shown on figure 2.c. We can see that each face of the map encodes a boundary face of the image. Inclusions of volumes are represented by a inclusion tree where each node corresponds to a region and its son nodes correspond to included regions. The root is R 0 , the infinite region which rounds up the image. When there is a hole on a face, the connexion between exterior and interior borders of this face is represented with a virtual edge. This is a one degree edge (that is adjacent twice to the same face). Such an edge should have been removed during the construction of the topological map, but have been retained to conserve the map connected, and each face homeomorphic to a topological disk. These edges are also useful to represent closed faces (without border), as we can see for the torus example shown figure 3. This torus is represented with a single closed boundary face. An intermediate step of the topological map's construction is shown on figure 3.b, where we can see two virtual edges (in grey) which keep the upper and lower faces connected. The topological map shown figure 3.c is only composed of virtual edges. We obtain the classic minimal representation of the torus composed of one face, one vertex and two edges.

Combinatorial maps encode only the topological part of our model. To encode the geometry of the corresponding image, we are using a geometrical model, which links a geometrical face to each topological face of the map . We are calling embedding this geometrical model. This distinction between topology and embedding allows the differentiation of treatments and sometimes enables a hierarchization of these treatments. Indeed, some operations only work on the topological model, others only on the geometrical model, and some on both. We can see figure 4.a the topological map of two adjacent objects and figure 4.a the embedding of these objects. Each dart of the topological is linked with the border of an embedding surface. For our example, the dart named a of the topological map is linked with the surface named 1, the two darts named b and c are linked with the same surface 2 because they belong to the same topological face, and the dart named d is linked with the surface 3.

Merge

The merging operation consists, starting from two adjacent regions R 1 and R 2 , to gather them into a single region R', union of the two first regions. Algorithm 1 make this operation on a topological map. It is local because it treats independently each dart of R 1 . There are two different cases depending if the current dart belongs to a virtual edge or not.

In the first case, we are testing if the deletion of this edge disconnects the map into several connected components. If it happends, we are modifying the inclusion tree this way. In the second case, we are just sewing the two faces adjacent to the currently treated edge together, and we are deleting this edge. This operation is made for all the darts of the boundary faces, so the two regions are merged and boundary faces are finally destroyed. At last, the topological map is simplified, because operations could have made the map incoherent or not minimal. For that, first we are removing degree 2 vertex and edges. Then, as virtual edges have just a topological existence, they can be moved to reduce the edge number. These simplification step guarantee the minimality of the topological map. Moreover, during this step, we are using the Euler characteristic to keep invariable all topological characteristics of the map. We can see figure 5 an example of the merging operation. Figure 5.a shows a topological map that represent three adjacent objects, before the merging of R 1 and R 2 . Figure 5.b represents the map obtained after the merging. We have destroyed the face between R 1 and R 2 and β 2 -sewn each other darts incident to this face. We can see on this figure that this map is not a topological map because it is not the minimal representation. Indeed, we have two edges incident to only two faces. These edges can be removed without lost of topological information. This is done by the simplification step of our algorithm. After the removal of these edges, we obtain some vertices incident to only two different edges. They can be also removed without lost of information. Finally, we obtain the map shown figure 5.c that is the topological map representing two adjacent objects. We can note that this simplification step can be performed during the merging, but this leads to a less understandable algorithm, but also with a more efficient complexity. We chose to present here the more simple algorithm, even though its complexity is a little higher, to make easier its comprehension.

Split

This operation is the opposite of the merging one. It consists in spliting a region by a separation face. This face is composed of a list of vertices which represent the intersections of a plane and the edges of the map. Constraints have been fixed on this separation face to simplify the algorithm and to limit the number of different cases to treat. A simple solution to do complex splits is to combine several simple splits and merge operation to obtain the wanted result. The principle of the algorithm is first to insert a new boundary face into the embedding then the corresponding boundary face in the topological map which represents the region.

There are three different cases :

1. If the separation face of R 1 cuts the edge of an already existing boundary face, like the example shown figure 6.a. New vertices will be inserted into this boundary face then they have to be separated with inserted edges between these vertices. At last the new boundary face can be inserted in the middle of the region. It is sewed on this new edge. We can see on figure 6.c the resulting topological map of the operation of spliting.

R 1 1 R R 1 β 2 β 3 β 2 1 R R 2
a. Two adjacents cubes. b. The corresponding map. c. After the split. Insert the 4 faces between F and the face which was β2-sewed to F ; Simplify the map;

R 1 R 1 R 1 ' β 3 β 2 β 2 β 2 β 3 β 2
a. Two adjacents cubes. b. After the split. 3. If the separation face cuts a closed face. A closed face is only represented with virtual edges in the topological map. There is an example on figure 8.b with its embedding in a. To obtain the minimal map, the new inserted boundary faces are only composed with one dart and one of this faces has to be directly β 1 -sewed to the closed face. The result of the split for a torus is shown on figure 8.c. We can verify with the Euler formula that the characterization of the torus is preserved. There is now 1 vertex, 3 edges and 2 faces, which gives a genius of 1.

Corefining

One method to segment a big 3D image consists cutting it into several parts and segmenting each part in parallel. During the reconstruction of the image, we have to glue back the different segmented parts. The operation which performs this glue back is the corefining. To design the two faces of each volume which are going to be put in contact, we are passing two embedding darts to the algorithm. The faces which contain these darts gives the coordinates of the corresponding planes (as we works with interpixel). In a first time, we build the new boundary faces. For that, we extract the embeddings of the two faces belonging to the two planes, because we consider that the initial image was cut in regular parts. We are building the map of these new border faces A and B. This construction can be different if the new boundary face cuts a closed face or the border of an existing boundary face, in a similar way that the three different cases explained for the split. We can see figure 9.a an example of this first step.

In a second step, we are inserting the intersection vertices of the two faces A and B in the embedding of the two maps of these two faces. These two maps have now exactly the same vertices, topological as well as geometrical. We are β 3 -unsewing and destroying faces belonging to R 0 which are previously β 3 -sewn to A and B. Then, for each dart of A, we are looking for a similar dart in B (a dart which connects the same vertices). If such a dart does not exist, we are adding in B a copy of this dart. And we are exactly doing the same operation for each dart of B, with eventually adding some copies in A. Then we are obtaining two faces A and B totally similar, as we can see on figure 9.b. We just have to β 3 -sew these two faces. For that, we are using the same coverage of the darts of A and B, and β 3 -sew each couple of darts. We can see on figure 9.c the result of this final operation.

Algorithm 3 Corefining 3D

Data: Two maps M1 and M2 and two darts d1 and d2 of their embeddings.

Result: The maps M1 and M2 are β3-sewn by the faces containig d1 and d2.

Step

Conclusion

In this paper, we have presented the basic operations for the 3D segmentation : merge and split plus an interesting operation in a parallelization goal : corefining. These operations are defined on the 3D topological map, a mathematical model of 3D images representation. The two first operations were already presented in [START_REF] Braquelaire | 3d split and merge with 3-maps[END_REF], but our approach differs from this solution by our aims to define local algorithms. The last operation is very interesting, because it allows us to segment big images by working in parallel onto different small parts of the image. Such a segmentation is very difficult to perform if we want to do in a direct way.

We have totally used all the topological map properties to obtain algorithms the most generic as possible. We are thus obtaining operations working on every possible case. Moreover, these algorithms operate in a local way, which simplifies a lot the different treatments. Indeed, we are treating each element without particular order, and we just have to look at the direct neighborhood of the current element. These different advantages allow us to obtain algorithms relatively simple and understandable, but also to keep a good complexity. Now, we have to implement these three operations in our computer software, to obtain finally some 3D segmentation algorithms. We are currently working to use these algorithms to perform segmentation refinement. Our approach consists in starting from a first segmentation, performed in a classical way on the voxel matrix with some algorithms used in the signal treatment research. Then, we are computing the topological map corresponding to this first segmentation. We can thus perform some treatments onto this map to refine this segmentation : automatic treatments, for example to remove small regions or particular configurations, or interactifve treatments made by an expert. These results are under development to obtain a software for cerebral tumor diagnostic. Moreover, we also work on some other operations, that could be combination of these three basic operations, but also some specific ones as the chamfering or boolean operations.

Fig. 1 .

 1 Fig. 1. An object and the corresponding combinatorial map represented by two different way.

  boundary faces. c. The topological map.

Fig. 2 .

 2 Fig. 2. A 3D image and the corresponding topological map.

  a. A torus.b. Intermediary step. c. The topological map.

Fig. 3 .

 3 Fig. 3. The topological map of a torus.

Fig. 4 .

 4 Fig. 4. Embedding example of a topological map.

Algorithm 1

 1 Merge 3D Data: Two adjacent regions R1 and R2. Result: The two regions are merged into R2. foreach dart b of the region R1 do if b belongs to the boundary between R1 and R2 then if b belongs to a virtual edge then foreach dart t of the orbit < β1, β2 > of β1(b) and β13(b) do if all regions of β3(t) are included into R1 ∪ R2 then Daughter(R2)=Daughter(R2)+regions of the connected component of β3(t) Destroy the virtual edge incident to b; else β2-sew(β2(b), β32(b)); Destroy β3(b) and b; Simplification of R2;

  topological map. b. After the merging. c. After the simplification.

Fig. 5 .

 5 Fig. 5. An example of Merge 3D

Fig. 6 .Algorithm 2

 62 Fig. 6. Case 1 of Split 3D

Fig. 7 .

 7 Fig. 7. Case 2 of Split 3D

  corresponding map. c. After the split.

Fig. 8 .

 8 Fig. 8. Case 3 of Split 3D

Fig. 9 .

 9 Fig. 9. Corefining principle.


Laurent Lapasset

Clément Peyronne

Benoit Hermeline

Conflict-free trajectory optimization using B-Splines and Genetic Algorithm

come   L'archive ouverte pluridisciplinaire

Conflict-free trajectory optimization using B-splines and Genetic Algorithm

Contacts : Clément Peyronne clement.peyronne@capgemini.com -Ingénieur R&D ATM Laurent Lapasset laurent.lapasset@capgemini.com -Ingénieur R&D ATM Benoît Hermeline benoit.hermeline@capgemini.com -Manager ATM Unit

Trajectory model and chromosome encoding

• B-splines are calculated from the control points designed by the GA

• One trajectory is determined by its control points Optimization method: Genetic algorithm A population of solution (chromosome) evolves using evolution concepts Crossover Mutation

A chromosome contains the control points for all the aircrafts involved

The roundabout test problem resolution Conflict detection and fitness calculation :

• Fitness expression:

• Conflict detection using a space discretization:

Main Idea :

Our methodology aims to generate optimal conflict-free trajectories for an en-route conflict situation. The genetic algorithm (GA) pilots the B-splines through their control points. Consequently, trajectories are directly monitored by the GA which handles both deconfliction and optimization.

About Capgemini

Capgemini, one of the world's foremost providers of consulting, technology and outsourcing services, enables its clients to transform and perform through technologies. Capgemini provides its clients with insights and capabilities that boost their freedom to achieve superior results through a unique way of working, the Collaborative Business Experience TM . The Group relies on its global delivery model called Rightshore ® , which aims to get the right balance of the best talent from multiple locations, working as one team to create and deliver the optimum solution for clients. Present in more than 35 countries, Capgemini reported 2009 global revenues of EUR 8.4 billion and employs over 100,000 people worldwide. More information is available at www.capgemini.com.

  1. Fill each cell with entry and exit time for each airplane 2. Check for time conflict the eight neighboring cells for each cell between a stochastic optimization method and a smooth trajectory model: Genetic Algorithm (GA) B-splines ✓ B-splines are controlled by GA via their control points Our methodology: