Fannie Castignoles
email: fannie.castignoles@univ-st-etienne.fr

Manuel Flury

Thierry Lépine

Thierry Lepine

Comparison of the efficiency, MTF and chromatic properties of four diffractive bifocal intraocular lens designs

come

Comparison of the efficiency, MTF and chromatic properties of four diffractive bifocal intraocular lens designs

Introduction

Multifocal intraocular lenses are nowadays widely used in the correction of presbyopia, offering the pseudophakic patient far and near vision with one implant [START_REF] Davison | History and development of the apodized diffractive intraocular lens[END_REF].

The profile of such bifocal diffractive lenses can be binary or continuous with various forms (parabolic, hyperbolic, sinusoidal…). This profile will determine the optical properties of the lens, notably in terms of efficiency, MTF (modulation transfer function) and chromatism.

In the present manuscript we will study four important multifocal diffractive profiles (a parabolic one, a parabolic profile with holes, a sinusoidal one and a binary one) with a common addition of + 3 dioptres ( + 3 D is the power of the IOL lens which corresponds to an equivalent of + 2.4D of glass addition), and we will compare their properties in a radar graph that highlights their respective characteristics.

The first criterion retained for this study is the diffractive efficiency at 550 nm, which will be obtained with a Hankel transform of the profile at stake. In a second part, a simplified model in Zemax® is employed to compare the MTF of the near and the far focus, at 550 nm. Eventually in a third part, the diffractive efficiency and the MTF at 440 nm and 640 nm are calculated, shedding light on the effects of chromatism due to the different profiles. All theses data are then gathered in a single figure in part 6.

Description of the lenses designs

Multifocal intraocular lenses used to correct presbyopia are made of a carrier lens which determines the power for the far vision of the bifocal, whereas the diffractive profile provides the added power (used to correct the near vision) [START_REF] Cohen | Diffractive bifocal lens designs[END_REF]. In this paper we will only study planar lenses in order to focus on the addition of power and the bifocality provided by the diffractive profile. This work can be easily transposed to a curved IOL, if the diffractive orders are moved on the power scale according to the refractive power added to the base of the lens.

All the lenses studied in this paper have a diameter of 3 mm, corresponding to the average diameter for eye pupils in photopic light. The choice of the addition is + 3D which corresponds to a focal length of 333 mm. The radii of the annular concentric zones depend on the focal length. In a parabolic profile, a parabolic profile with holes, and a sinusoidal one the rings are spaced according to the formula [START_REF] Moreno | High efficiency diffractive lenses: Deduction of kinoform profile[END_REF]:

( ) 2 R k k f λ = ⋅ ⋅ ⋅ (1)
where R (k) is the zone radii and k is the integer which corresponds to the k th zone. Otherwise in the binary profile which consists of full and empty zones, the radii of the full zones are equal to [START_REF] Moreno | High efficiency diffractive lenses: Deduction of kinoform profile[END_REF]:

( ) ' , R k k f k odd λ = ⋅ ⋅ (2)
Thus the parabolic, parabolic with holes, and sinusoidal lenses with a 3 mm diameter contain 6 annular concentric zones, and the binary profile contains 6 empty and 6 full zones. In order to have a bifocal lens which divides light fairly equally between the two images, the steps of the diffractive profile must delay light of ½ wavelength [START_REF] Cohen | Diffractive bifocal lens designs[END_REF][START_REF] Moreno | High efficiency diffractive lenses: Deduction of kinoform profile[END_REF]. The height of the steps is equal to:

max 0.5 h n λ ⋅ = ∆ ( 3
)
where ∆n is the difference in the refractive index between the lens and its environment, and λ is the design wavelength of the lens (here λ = 550 nm). However an exception is made to the sinusoidal profile, which needs a maximum height of:

' max 0.415 h n λ ⋅ = ∆ (4)
to become bifocal with the same efficiency in the near and far focus [START_REF] Cohen | Diffractive bifocal lens designs[END_REF]. The four different profiles examined in this paper are represented in Fig. 1. We have chosen to compare the results of a binary diffractive profile, a parabolic profile, a parabolic with holes, and a sinusoidal one.

The binary profile [Fig. 1(A)] is quite different from the other profiles, owing to the fact that its main diffractive orders are the -1st and 1st [START_REF] Kress | Digital Diffractive Optics: An introduction to planar diffractive optics and related technology[END_REF][START_REF] Baude | Caractérisation, optimisation et réalisation de composants optiques multifocaux diffractifs et réfractifs[END_REF] whereas the parabolic, parabolic with holes, and sinusoidal profiles diffract light into the zero (far focus) and first order (near focus) [START_REF] Kress | Digital Diffractive Optics: An introduction to planar diffractive optics and related technology[END_REF][START_REF] Baude | Caractérisation, optimisation et réalisation de composants optiques multifocaux diffractifs et réfractifs[END_REF]. The binary profile is described in Eq. [START_REF] Baude | Caractérisation, optimisation et réalisation de composants optiques multifocaux diffractifs et réfractifs[END_REF]. In all the following equations, the integer k corresponds to the k th ring, which is equal to:

r k ceil R   =     . ( ) max '( ) ' ( 1), 0 '( 1) '( 2)
,

h if R k r R k k even h r if R k r R k k even < < +  =  + < < +  (5)
The Fig. 1(B) represents the profile of a parabolic lens which is highly employed in ophthalmology [START_REF] Cohen | Diffractive bifocal lens designs[END_REF][START_REF] Moreno | High efficiency diffractive lenses: Deduction of kinoform profile[END_REF][START_REF] Cohen | Progressive intensity phase bifocal[END_REF]. Its profile is given by the Eq. ( 6):

( )

2 max max 2 , (1)
r h r h k h k R = ⋅ - ⋅ ∈ Ν (6)
The Fig. 1(C), corresponding to the Eq. ( 7), is also a parabolic profile but with full and empty steps as for the binary lens. In fact such profile can be useful for applications where holes are required in the lens.

( )

2 max max 2 1 ² '( ) '( 1), , ² 0.5 ² ² (1) 0 '( 1) '( 2), k r h h if R k r R k k even h r R if R k r R k k even α α α α  -+   ⋅ - ⋅ < < + =    = ⋅     + < < +  (7)
The last profile, Fig. 1(D) and Eq. ( 8) is sinusoidal, and can be employed in some applications instead of the parabolic one. In particular, a sinusoidal profile advantageously exhibits thinner and smoother steps [START_REF] Cohen | Diffractive bifocal lens designs[END_REF][START_REF] Cohen | Progressive intensity phase bifocal[END_REF].

( )

' 2 max 2 1 ( 1) cos( ) 2 (1) k h r h r R π    = ⋅ --⋅ ⋅      (8)

Diffractive efficiencies at 550 nm

Methods

One way to calculate the diffractive efficiencies is to use a Fourier transform [START_REF] Baude | Caractérisation, optimisation et réalisation de composants optiques multifocaux diffractifs et réfractifs[END_REF]. In fact when the transmission function t is periodic in r 2 , it can be expressed as a Fourier sum:

( )

² exp 2 . ( 1
)² p p Z r t r C i p R π ∈   = ⋅     ∑ (9)
where R( 1) is the period of the function t in the r 2 space and the p C are the coefficients given by:

(1)² 0 1 ( ) exp 2 (1)² (1)² R p u C t r i p du R R π   = - ⋅ ⋅     ∫ (10)
The diffractive efficiency p η of each order p will be the squared modulus of its coefficient in the Fourier series of the transmission function t:

² p p C η = (11)
However, even if the Fourier transform is a rapid and useful tool to calculate the efficiencies of a diffractive lens, it is limited to transmission functions periodic in the r 2 space, which is not always the case. For instance sawtooth profiles cannot be described with the Fourier transform. Moreover the pupil is not taken into account in this method.

That is why another method based on wave propagation is employed. This method is based on the theory of scalar diffraction and allows the efficiency of each kind of design to be calculated. Moreover the propagation of the field can be studied and represented, which is of interest to observe the variation of focal length of our systems.

This wave propagation is numerically implemented in Matlab, as is described in reference [START_REF] Guizar-Sicairos | Computation of quasi-discrete Hankel transforms of integer order for propagating optical wave fields[END_REF]. Quasi-discrete Hankel transform of the 0 order is used to adapt the circularly symmetric function of the different transmission functions.

The lens is normally illuminated with an incident plane wave. The efficiencies of the convergent orders (negatives orders in our convention) are calculated from the square modulus of the field. Thus for each order the energy contained in the first ring of the corresponding Airy pattern (at the focal plane and its fractions f/2, f/3 …) corresponds to 84% of the diffractive efficiency in this order.

The efficiencies of the divergent orders (positives ones), are calculated in the same way, applying the Hankel transform to the complementary phase of the diffractive lens.

Finally the calculation of the zero-order efficiency requires the use of an additional perfect refractive lens to determine the efficiency of the diffractive lens in its focal plane. In Fig. 2 are represented the distributions of the diffractive efficiencies between the orders -2 to + 2, for each of the four profiles described in Fig. 1. The total efficiencies plotted in Fig. 2 do not add to 100% because the orders greater than +/-2 have finite but small efficiencies too.

Results

As it was explained previously, the diffractive efficiency of the binary profile [Fig. 2(D)] is shared between the -1st and 1st order with about 37% of efficiency. Moreover there is no efficiency for the even orders which is a characteristic of binary diffractive lenses. The diffractive efficiencies of the parabolic profile [Fig. 2(A)] are about 40% but affect the -1st and 0th orders. Furthermore its -2nd and 1st orders have a non negligible efficiency of about 5%. The parabolic profile with holes [Fig. 2(B)], has slightly lower diffractive efficiencies in the -1st and 0th order compare to the parabolic lens. It has also a non negligible efficiency in the -2nd and 1st orders (reaching 10%), which can lead to halos. The sinusoidal profile [Fig. 2(C)], keeps diffractive efficiencies close to 40% in the -1st and 0th orders, but exhibits also quite important efficiency in the first order.

Simulation of the MTF at 550 nm

A second performance characteristic of the different profiles is the modulation transfer function, which provides the spatial frequency response of the imaging system [START_REF] Lang | Interpreting multifocal intraocular lens modulation transfer functions[END_REF]. The MTFs are calculated with Zemax®, which is a software used to design optical systems and perform optical simulations. The four diffractive lenses are created from their radial coordinates inserted in a .TOB file in Zemax®. The resulting non sequential objects are then inserted in a very simplified model of the eye described in Fig. 3.This optical system allows to simulate the MTF at the focal plane of the eye, for both a near and a far object. The results are plotted in Fig. 4. The X-axis is in line pairs per millimeter. Hundred line pairs per millimeter correspond to a vision of 20/20.

The diffraction limit of the optical system is given on the figure as a reference.

The curves of the modulation transfer function of the four studied profiles are almost equivalent for the case of the near focus [Fig. 4(A)]. Those MTFs reach about 25% for 100 cycles per millimetres. Regarding the far focus [Fig. 4(B)], the MTFs of the parabolic and the sinusoidal profile are very similar and reach also 25% for 100 lp/mm. The MTF of the binary profile for the far focus is the MTF corresponding to the + 1st order. At 550 nm the optical quality of the binary profile is equivalent to those of the parabolic and sinusoidal ones, for both near and far vision. In return the drawback of the parabolic profile with holes is a weaker MTF for the far vision.

Effects of chromatism

Introduction

An important drawback of a diffractive lens is its wavelength dependency [START_REF] Baude | Caractérisation, optimisation et réalisation de composants optiques multifocaux diffractifs et réfractifs[END_REF][START_REF] Cohen | Diffraction bifocal with adjusted chromaticity[END_REF]. In fact the wavelength plays directly a part in the radius of the circular rings [Eq. ( 1), ( 2)], and in the height of those rings [Eq. ( 3), ( 4)]. Moreover the variation of λ implies also a variation of the refractive index n, which in its turn causes a modification of the radius and height of the rings. However, this dispersion of the material n (λ) is insignificant compared to the variation of λ itself. For example, with Polycarbonate (a very dispersive material with νd = 30) in air, and for λ varying from 440 nm (blue) to 640 nm (red),

d
- - - ∆ - = = = ∆ - - (13)
The variations of the radii and the height of the rings have an impact on the focal length and the optical power of the lens, on the diffractive efficiencies and on the MTFs. These produced modifications are compared in this paragraph for the four different profiles. The results are calculated at 440 nm (blue), 550 nm (design wavelength) and 640 nm (red), using the Hankel transform and Zemax® as it was described before.

Effects of chromatism on the focal length and the optical power

Given that the focal length is linked to the radii of the annular zones [Eq. ( 1), ( 2)], the variation of λ will have an effect on f and as a direct consequence on the optical power of the lens [START_REF] Cohen | Diffraction bifocal with adjusted chromaticity[END_REF]:

0 0 P P P P λ λ λ λ ⋅ ∆ = ⇒ ∆ = ⋅ (14)
According to the Eq. ( 14), the variation of the power of the lens between blue and red (Table 1) depends only on the value of the power for the calculated order. We can thus calculate ∆P for the first convergent order (-1st) which will be the same for the four profiles studied here: Thus the power deviation equals: ∆P = P red -P blue = 3.496-2.404 = 1.09 D, where P red is the power of the lens for red light and P blue the power for blue light. ∆P is positive which is characteristic of diffractive convergent components.

Besides, the variation in power for the higher convergent orders will increase proportionally to the power. Thus, ∆P corresponding to the -2nd order will be twice as high as ∆P for the -1st one. However in the case of the parabolic, parabolic with holes, and sinusoidal profiles studied in this paper, the near focus is given by the first negative order, and higher orders will thus not be considered.

Noteworthy, the far vision provided with the 0th order has no variation in its focal length (infinite) caused by chromatic effects.

In return the two orders (-1st and 1st) used for near and far visions with the binary profile suffer from chromatic effects. Furthermore, the positions of the red and blue foci for the far vision ( + 1st divergent order) will be the opposite of the near foci ones, so that the binary profile displays the disadvantage of being added to the chromatism of the human eye [START_REF] Cohen | Diffractive multifocal optical device[END_REF][START_REF] Thibos | The chromatic eye: a new reduced-eye model of ocular chromatic aberration in humans[END_REF].

Effects of chromatism on the distribution of the diffractive efficiencies

The change in the phase difference with λ will lead to a variation of the distribution of the efficiencies, and the MTFs.

Thus the variations of the diffractive efficiencies of the two interesting orders for each profile are plotted versus λ in Fig. 5 . They are calculated with the Hankel transform. The three parabolic, parabolic with holes, and sinusoidal profiles studied here behave the same way, so that the convergent -1st order decreases with λ, whereas the 0th order rises. Near 550 nm those two orders cross each other and reach about 40% of efficiency for the parabolic and the sinusoidal profile, and about 35% for the parabolic profile with holes. This latter has lower efficiencies throughout the curve. At 440 nm the three profiles proposed here are much more efficient for the near vision. On the contrary, at 640 nm, the dominant order is the far vision one.

The binary profile has a very different behavior, given that the -1st and + 1st orders have the same variation of efficiency with λ. Moreover this variation is much weaker than the other. Regarding the efficiencies' variations with λ, the binary profile is much more balanced than the others. This phenomenon can be intuitively understood since the complementary form of a binary lens is also a binary lens. This complementarity explains the balanced results between the positive and negative orders, whatever the wavelength.

Effects of chromatism on the MTFs

The impact of a wavelength change is also high in the MTFs (Fig. 6). First, the MTFs are now calculated in the focus corresponding to the blue and red foci. Moreover, the variation of the diffractive efficiencies has positive or negative impacts on the MTFs. In fact, the behavior of the MTFs is directly linked to the variations of efficiency. Thus the higher the diffractive efficiency is, the better the MTF curve will be. This is due to the fact that higher diffractive efficiency in a useful order means lower stray light, and thus less contrast lost.

Thus the phenomena observed with the variations of efficiency are similar for the MTFs. Indeed whereas the MTFs in blue and red for the binary profile are similar, it is obvious that the parabolic, parabolic with holes, and sinusoidal profile's MTFs are not balanced. Whereas in blue light the near vision becomes strongly better than for the binary profile, the situation is reversed in red light. On the contrary the MTF of the binary profile is higher than the three other profiles' MTFs for the far object in blue, and inversely in red light.

Synthesis of the different properties of each profile in a radar graph

In order to collect all the data previously obtained, a new way of illustrating the different properties of each profile was employed in a "radar-type" graph (Fig. 7). The diffractive efficiencies at 440 nm, 550 nm and 640 nm are plotted for the two useful orders. Besides punctual data on the MTF curves representative of the quality of the vision have to be selected. Thus the values of the MTF for near and far objects have been plotted at a frequency of 100 lp/mm, which corresponds to a visual acuity of 20/20. Furthermore it is interesting to represent the properties for the near vision in front of those of the far one. This underlines the symmetry of the results of the binary profile for near and far vision, contrary to the dissymmetry of the parabolic, parabolic with holes, and sinusoidal profiles behaviors. Advantages and drawbacks of each design can be assessed by the superimposition of those figures in a same graph (Fig. 8). The choice of the profile will then depends on the specific needs of the designer. For instance the use of a binary profile can be interesting for a bifocal lens requiring two balanced diffractive orders, and a small dependence with λ. Otherwise the superimposition of the merit figures for the parabolic and the sinusoidal designs highlights the similarity of their results. In return, the inclusion of the figure for the parabolic profile with holes inside the parabolic and sinusoidal profiles figures shows that only diffractive lenses that absolutely require holes can take advantage out of it.

Conclusion

This paper provides a study of the properties of four designs that can be used for multifocal diffractive lenses with a view to correcting presbyopia.

The binary profile, which diffracts light into the -1st and + 1st order, and the parabolic one, often used in diffractive IOL are well known in the art. The third one, a sinusoidal profile is sometimes used as a variant of the parabolic one. The study of the distribution of the diffractive efficiencies and the MTFs for several wavelengths shows that those two latter profiles exhibit very similar behaviors. On the contrary, the binary profile gives quite different results. Indeed its symmetric behavior leads to two main orders very balanced and less wavelength dependent. However other criterion like the sign of its chromatism must be considered before use. Besides, an intermediate design made of a parabolic profile with holes (like the binary one) has been compared to the previous ones. The results obtained with the Hankel transform and the Zemax® simulations show a slightly lower efficiency and MTF compared to the other profiles. It can still be used in a system where holes are needed.

Finally a radar graph (Fig. 8) is proposed to better comprehend the different parameters involved in the choice of a specific profile. This graph takes into account the distribution of efficiencies of the two main diffractive orders at 440 nm, 550 nm, and 640 nm and draws a parallel with the corresponding MTFs which illustrates the variation in image quality. However although this radar graph consists in a useful tool, the choice of a profile should also take into account the fabrication technology. For instance the choice of a binary lithography technology will be more appropriated for a binary profile, whereas the use of diamond tool or grey-tone lithography enables to design continuous profiles. Other parameters like geometrical constraints of the lens may also play a role in the choice of the profile.

Fig. 1 .

 1 Fig. 1. the four diffractive profiles studied in this paper. (A) Binary profile; (B) Parabolic profile; (C) Parabolic profile with holes; (D) Sinusoidal profile.

Fig. 2 .

 2 Fig. 2. Distribution of the diffractive efficiencies at 550 nm for the four different lenses. (A) Parabolic profile; (B) Parabolic profile with holes; (C) Sinusoidal profile; (D) Binary profile

Fig. 3 .

 3 Fig. 3. Optical system used in the simulation of the MTF with Zemax® for near and far objects.

Fig. 4 .

 4 Fig. 4. MTF of the four profiles at 550 nm, for near (A) and far object (B).

Fig. 5 .

 5 Fig. 5. Variations of the diffractive efficiencies versus λ

Fig. 6 .

 6 Fig. 6. MTFs for the near and far vision, at 440 nm and 640 nm.

Fig. 7 .

 7 Fig. 7. radar graphs for the four different profiles.

Fig. 8 .

 8 Fig. 8. Superimposition of the radar graphs for the four different profiles.

Table 1 . Variation of the convergent focal length with wavelength.

 1

	λ (nm)	Focal length (mm)	Optical power (dioptres)
		of the -1st order
	440	416	2.404
	550	333	3.0
	640	286	3.496

Acknowledgements

The authors wish to thank G. Cohen and P. Chavel for their helpful discussions and the National Association for Research and Technology (ANRT) for its financial support.


Dorothée Salvatore Monaco
email: salvatore.monaco@uniroma1.it

Mattia Mattioni
email: mattia.mattioni@l2s.centralesupelec.fr

D Normand-Cyrot

Stabilization of feedforward discrete-time dynamics through immersion and invariance

Keywords: Nonlinear output feedback, Stability of nonlinear systems, Lyapunov methods

The paper deals with the problem of stabilizing discrete-time feedforward dynamics through Immersion and Invariance. Closed loop stabilization of the equilibrium is achieved making use of a passivity-based controller combined with a domination argument. A simulated example illustrates the performances.

I. INTRODUCTION

Forwarding is a stabilizing approach developed in continuous time (see [START_REF] Kristic | Feedback linearizability and explicit integrator forwarding controllers for classes of feedforward systems[END_REF], [START_REF] Sepulchre | Integrator forwarding: a new recursive nonlinear robust control design[END_REF]) for dynamics admitting a particular cascaded (or triangular) structure. It provides a systematic bottom up recursive Lyapunov-based design procedure which can be interpreted as the dual of the celebrated back-stepping one [START_REF] Kokotović | Constructive nonlinear control: a historical perspective[END_REF]: instead of assuming a state component as a virtual control and controlling through forwarding, stabilization is achieved by iteratively adding a state component which "integrates" the other ones. Such an approach has been developed in discrete time in [START_REF] Mazenc | Forwarding in discrete-time nonlinear systems[END_REF] to stabilize classes of nonlinear dynamics under Lyapunov-based bounded feedback. In the present work we consider discrete-time strictfeedforward dynamics of the form x jk+1 = x jk + F j (x 1k , x 2k , . . . , x j-1k , u k ); j = 2, . . . , n

x 1k+1 = F 1 (x 1k , u k ) (1)
where x 1 ∈ R p and the x j ∈ R for j = 2, . . . , n and u ∈ U ⊆ R m , m = 1. The proposed control strategy extends the design introduced in [START_REF] Monaco | Forwarding stabilization of nonlinear discrete-time dynamics[END_REF] for systems in strict-feedforward form (when setting x 1k+1 = u k , m = 1). Here, the general non linear x 1 -dynamics is known to be stabilizable under a suitable state feedback but the knowledge of a control Lyapunov Function is not assumed. The problem is presently set in the context of Immersion and Invariance -I&I. I&I was proposed in [START_REF] Astolfi | Immersion and invariance: a new tool for stabilization and adaptive control of nonlinear systems[END_REF]- [START_REF] Astolfi | Nonlinear and adaptive control with applications[END_REF] for stabilizing continuous-time systems and reformulated in [START_REF] Yalcin | Discrete time immersion and invariance adaptive control for systems in strict feedback form[END_REF]- [START_REF] Mattioni | Digital stabilization of strict feedback dynamics through invariance and immersion[END_REF] in discrete-time; the overall design results to be less demanding in such a context.

With reference to (1), a preliminary controller ensuring global asymptotic stability -GAS -of the origin of the x 1dynamics is computed so defining the target system. Then the design is completed for ensuring attractiveness of a certain invariant set associated with the target x 1 -dynamics. Boundedness of the state trajectories guarantees global asymptotic stabilization of the closed-loop equilibrium. Attractiveness of the invariant set is achieved thanks to passivity arguments and negative output feedback as proposed in [START_REF]Nonlinear avarage passivity and stabilizing controller in discretetime[END_REF]. As usual in a discrete-time context, the negative output feedback is only implicitly defined; for this reason, domination arguments as developed in [START_REF] Mazenc | Adding integrations, saturated controls and stabilization for feedforward systems[END_REF]- [START_REF] Mazenc | Forwarding in discrete-time nonlinear systems[END_REF] are used to provide an explicit bounded solution.

The paper is organized as follows. Some preliminaries on discrete-time state space representations and I&I stabilizability are in Section II. A brief motivating discussion is in Section III. The proposed control design is developed for the elementary feedforward two block dynamics in Section IV. An example concludes the paper in Section V.

II. RECALLS AND PRELIMINARIES

A. The Differential Difference Representation

Following [START_REF] Monaco | A unifying representation for nonlinear discrete-time and sampled dynamics[END_REF], nonlinear discrete-time dynamics in the form of a map x k+1 = F (x k , u k ) with F (•, u) a R n -valued smooth map, smootly parameterized by u ∈ U, can be represented as a couple of a difference and a differential equation

x + = F 0 (x) (2) ∂ x + (u) ∂ u = G (x + (u), u) (3)
where x + (u) represents a curve in R n parameterized by u, F 0 (x) := F (x, 0) := x + (0), defines the initial condition of the differential equation [START_REF] Kokotović | Constructive nonlinear control: a historical perspective[END_REF]. G (•, u) on R n , parameterized by u, is computed to satisfy the equality

G (F (x, u), u) = ∂ F (x, u) ∂ u .
Given F (•, u), the existence of G (•, u) is ensured by the reversibility of F 0 (•) as a function of x, so uniquely defining and computing G (•, u) for u sufficiently small as

G (x, u) := ∂ F (•, u) ∂ u x=F -1 (x,u) (4)
where F -1 (x, u) denotes the reverse function, i.e., satisfying

F (F -1 (x, u), u) = x.
For any x, completeness of the vector field G (•, u) for all u ∈ U ensures integrability of (3) so recovering the usual representation in the form of a map

x + (u) = F (x, u) = x + (0) + u 0 G (x + (v), v)dv
with x + (0) := F 0 (x). Consequently, given any smooth function H : R n → R, its variation with respect to u around H (x + (0)) admits the integral form representation

H (F (x, u)) -H (F 0 (x)) = u 0 L G (•,v) H (x + (v))dv where L G (•,•) H (•) represents the Lie derivative of H (•) along the vector field G (•, •); i.e. L G (•,v) H (x + (v)) = ∂ H (•) ∂ x G (•, v) x + (v)
.

(2) and (3) define the (F 0 , G ) representation of a discretetime dynamics.

B. (F 0 , G ) representation of feedforward dynamics

According to these definitions it is a matter of computations to verify that any nonlinear feedforward dynamics of the form (1), with F 1 (x, u) reversible in x 1 for all u ∈ U, admits the following (F 0 , G ) representation

x + j = x j + F j0 (x j-1 , . . . , x 1 ); j = 2, . . . , n x + 1 = F 10 (x 1 ) (5)
∂ x + j (u) ∂ u = G j (x + j-1 (u), . . . , x + 1 (u), u); j = 2, . . . , n ∂ x + 1 (u) ∂ u = G 1 (x + 1 (u), u) (6)
with

F 0 (x) = col[x n + F n0 (x n-1 , . . . , x ), . . . ., x 2 + F 20 (x 1 ), F 10 (x 1 ] F 0 (0) = 0 and G (•, u) := col[G n (•, u), . . . , G 1 (•, u)].
For, it is sufficient to verify that reversibility of F 1 (x 1 , u) in x 1 is sufficient to imply reversibility of F 0 (x) in x. More precisely, the reverse dynamics can be iteratively computed so getting

x 1 = F -1 1 (x + 1 (u), u) x 2 = x + 2 (u) -F 2 (x 1 , u) = x + 2 (u) -F 2 F -1 1 (x + 1 (u), u), u x 3 = x + 3 (u) -F 3 x + 2 (u) - F 2 (F -1 1 (x + 1 (u), u), u), F -1 1 (x + 1 (u), u), u • • • .
Then, according to (4), one computes for j = 2, . . . , n, the control vector fields G j (x + j-1 (u), . . . , x + 1 (u), u) as

G j (x j-1 , . . . , x 1 , u) = ∂ F j (x j-1 , . . . , x 1 , u) ∂ u x=F -1 (x,u) G 1 (x 1 , u) = ∂ F 1 (x 1 , u) ∂ u x 1 =F -1 1 (x 1 ,u)
which maintain the required triangular form.

In the sequel, the design is instrumentally developed with reference to strict-feedforward dynamics which admit the (F 0 , G ) representation [START_REF] Monaco | Forwarding stabilization of nonlinear discrete-time dynamics[END_REF][START_REF] Astolfi | Immersion and invariance: a new tool for stabilization and adaptive control of nonlinear systems[END_REF], but the solution can be applied to a discrete-time system of the form [START_REF] Kristic | Feedback linearizability and explicit integrator forwarding controllers for classes of feedforward systems[END_REF].

C. The discrete-time I&I stabilizability conditions

Following [START_REF] Astolfi | Immersion and invariance: a new tool for stabilization and adaptive control of nonlinear systems[END_REF], let us preliminarily formulate I&I stabilization for generally nonlinear difference equations [START_REF] Mattioni | Digital stabilization of strict feedback dynamics through invariance and immersion[END_REF]. Theorem 2.1: Consider the nonlinear discrete-time dynamics

x k+1 = F (x k , u k ) (7)
with state x ∈ R n , control u ∈ R and let x * ∈ R n the equilibrium to be stabilized. Let p < n and assume that we can find mappings

α(•) : R p → R p ; π(•) : R p → R n φ (•) : R n → R n-p ; ψ(•, •) : R n×(n-p) → R such that the following hold. H 1d (Target System) -The dynamics with state ξ ∈ R p ξ k+1 = α(ξ k ) (8)
has a globally asymptotically stable equilibrium at ξ * ∈ R p and x * = π(ξ * ). H 2d (Immersion and invariance condition) -For all ξ ∈ R p , there exists c(•) : R p → R such that

F (π(ξ ), c(ξ )) = π(α(ξ )). (9)
H 3d (Implicit manifold) -The following identity between sets holds

{x ∈ R n |φ (x) = 0} = {x ∈ R n |x = π(ξ ) for ξ ∈ R p } (10)
H 4d (Manifold attractivity and trajectory boundedness)

-All the trajectories of the system

x k+1 = F (x k , ψ(x k , z k )) (11a) z k+1 = φ (F (x k , ψ(x k , z k ))) (11b)
with z ∈ R n-p , z 0 = φ (x 0 ), are bounded for all k ≥ 0 and satisfy lim k→∞ z k = 0 and ψ(•, 0)

x=π(ξ ) = c(ξ ).
Then, x * is a globally asymptotically stable equilibrium of the closed loop dynamics

x k+1 = F (x k , ψ(x k , φ (x k ))).
Definition 2.1: The discrete-time nonlinear dynamics ( 7) is said to be I&I-stabilizable with target dynamics

ξ k+1 = α(ξ k ) when H 1d -H 2d -H 3d -H 4d in Theorem 2.1 are satisfied.
We show in this paper that I&I stabilization of feedforward dynamics reduces the control design to the iterative stabilization of suitably defined one-dimensional dynamics.

III. A MOTIVATING DISCUSSION

To briefly discuss the idea behind forwarding, let the elementary discrete-time cascade-connected dynamics

y k+1 = y k + H(x k ); x k+1 = F(x k )
where x ∈ R p , y is scalar, H and F are continuous functions and the origin of the x-subsystem is assumed asymptotically stable; i.e. there exists a positive definite Lyapunov function V (x) such that V (x k+1 ) -V (x k ) < 0 for all x k = 0. The existence of an invariant stable set described by the graph of a function y = Φ(x), implies the stability of the cascade.

In fact, if such a function Φ(x) (with Φ(0) = 0) exists, it satisfies the invariance implication

(y 0 , x 0 ) ∈ Ω = {(y, x)s.t.y = Φ(x)} ⇒ (y k , x k ) ∈ Ω; ∀k > 0
for any (y k , x k ) solutions to the cascade system. Φ(•) can be computed by solving the equality

Φ(F(x)) -Φ(x) = H(x) (12)
with initial condition Φ(0) = 0. On these bases, it is easily verified that a Lyapunov function for the overall system is given by W (x, y) = V (x) + 1 2 (y -Φ(x)) 2 with first order increment

∆ k W = W (x k+1 , y k+1 ) -W (x k , y k ) negative definite equal to ∆ k V = V (F(x k )) -V (x k ).
Forwarding relies on this basic idea to stabilize controlled cascade-connected dynamics of the form

y k+1 = y k + H(x k ), x k+1 = F(x k , u k ).
In this case, according to the (F 0 , G) representation of the x-dynamics

x + k (u k ) := F(x k , u k ) = F 0 (x k ) + u k 0 L G(•,v) (x + (v))dv
with F 0 (x) := F(x, 0), one has 13)

∆ k W = V (F(x k )) -V (x k ) + u k 0 L G(•,v) V (x + (v))dv + 1 2 ( u k 0 L G(•,v) Φ(x + (v))dv) 2 (
-(y k -Φ(x k )) u k 0 L G(•,v) Φ(x + (v))dv.
It clearly comes out that choosing u k to render negative ∆ k W in (13) is a difficult task since it involves solving an implicit inequality in u k . An equivalent solution, which recalls the structure of the continuous-time one, can be obtained by rewriting ∆ k W as

∆ k W = V (F(x k )) -V (x k ) + u k 1 0 L G(•,su k ) V (x + (su k ))ds + u 2 k 2 1 0 L G(•,su k ) Φ(x + (su k ))ds 2 -u k (y k -Φ(x k )) 1 0 L G(•,su k ) Φ(x + (su k ))ds
and solving the implicit equality below

u = -Γ -1 1 0 [L G(•,su) V x + (su) -(y -Φ(x))L G(•,su) Φ x + (su)
]ds (14

)
with positive gain function

Γ = [1 + [ 1 0 L G(•,su) Φ x + (su) ds] 2 2 ]. (15)
Remark 3.1: The control solution (14) recalls the continuous-time stabilizing control

u = -[L g V (x) -(y -φ (x))L g φ (x)] (16)
which ensures global asymptotic stabilization of the system

ẏ = h(x), ẋ = f (x) + ug(x) (17)
when the origin of the x-subsystem is assumed asymptotically stable (there exists a positive definite Lyapunov function V (x) such that L f V < 0 for all x = 0) and when the function φ (x) is computed to satisfy L f φ = h with φ (0) = 0.

It will be shown in the sequel that the concept of I&I stabilization combined with a domination argument makes the design constructive in discrete-time too.

IV. I&I FORWARDING STABILIZATION Consider the following elementary feedforward dynamics over R p+1

x 2k+1 = x 2k + F 2 (x 1k , u k ), x 1k+1 = F 1 (x 1k , u k ) (18)
and assume that: A1) -the origin of the x 1 -subsystem is asymptotically stable, i.e. there exists a positive definite Lyapunov function V (•) : R p → R such that V (F 1 (x 1k , 0)) -V (x 1k ) < 0 for all x 1k = 0; A2) -there exists a function Φ 1 (x 1 ) with Φ 1 (0) = 0 satisfying the following equality

Φ 1 (F 10 (x 1 )) -Φ 1 (x 1 ) = F 20 (x 1 ) (19)
with F 10 (x 1 ) = F 1 (x 1 , 0) and F 20 (x 1 ) = F 20 (x 1 , 0). The (F 0 , G) representation of (18) takes the form:

x + 2 =x 2 + F 20 (x 1 ), ∂ x + 2 (u) ∂ u = G 2 (x + 1 (u), u) x + 1 =F 10 (x 1 ), ∂ x + 1 (u) ∂ u = G 1 (x + 1 (u), u) (20)
with

G 1 (F 1 (x 1 , u), u) = F 1u (x, u); G 2 (F 1 (x 1 , u), u) = F 2u (x 1 , u).
Setting

z 2 = x 2 -Φ 1 (x 1 )
and

G z (•, u) := G 2 (•, u) -L G 1 (•,u) Φ 1 (•) (21) with L G 1 (•,u) Φ 1 (F 1 (x 1 , u)) := ∂ Φ 1 ∂ x 1 (F 1 (x 1 , u))F 1u (x 1 , u)
, the following result holds true. Proposition 4.1: -Given the discrete-time feedforward dynamics (18) satisfying A1) A2), then the equilibrium is I&I stabilizable with target dynamics ξ k+1 = F 10 (ξ 1k ).

Proof: For, it is sufficient to show that the H 1d to H 3d conditions in Theorem 2.1 are satisfied. First, setting x 1 = ξ , the target dynamics is defined by the x 1 -dynamics in free evolution because of assumption A1) (i.e. ξ k+1 = F 10 (ξ 1k )). Then, the immersion mapping π(ξ ) is defined from R p to R p+1 as π(ξ ) = col(Φ 1 (ξ ), ξ ) so immediately verifying the invariance condition (9) because of assumption A2) with c(ξ ) = 0. Setting now z 2 = x 2 -Φ 1 (x 1 ), one expresses with (21) the z 2 -dynamics, driftless by construction of Φ 1 ; i.e.

z + 2 =z 2 ∂ z + 2 (u) ∂ u =G z (x + 1 (u), u). ( 22
)
By construction, the p-dimensional set described by z 2 = 0 is invariant and the I&I design reduces to find a control that makes such set attractive while guaranteeing boundedness of the closed loop trajectories. The control design is discussed below.

rewritten as in ( 29) is bounded and guarantees negativity of uH + av (x 1 , z 2 , u); i.e. for all (x 1 , z 2 ), one has |u(x 1 , z 2 )| ≤ µ and

u 0 z + 2 (v)G z (x + 1 (v), v)dv ≤ - 1 2 ε(x 1 )λ (x 1 , z 2 )|z 2 G z (F 10 (x 1 ))| 2 .
When z 2 = 0, u(x 1 , 0) = 0. Boundedness of the whole state trajectories follows from Proposition 4.3 with sufficiently small ε(x 1 ) > 0.

Remark 4.1: Theorem 4.1 can be used repeatedly to deduce that the feedforward system (1) is globally asymptotically stabilizable under A1) and successive conditions of the type A2). At each step of the design, the closed loop col(x 1 , z 2 )-dynamics defines a new stable x1 -dynamics in R p+1 which defines a new target dynamics and one adds an integrating variable z 3 = x 3 -Φ( x1 ) which satisfies the invariance condition Φ( F10 ( x1 )) = F30 ( x1 ). This is repeated for the n blocks.

Remark 4.2:

The choice of the gain λ (x 1 , z 2 ) in ( 27) can be made (within the suitable interval) according to several strategies. In particular, it can be chosen so as to obtain robustness of the closed-loop system under the nominal control with respect to parameters uncertainty and external disturbances. In this case, suitable Lipschitz-like assumptions should be introduced for the mappings F i (i = 1, 2).

Remark 4.3:

The same approach can be pursued when relaxing assumption A1) to the existence of a preliminary stabilizing feedback c(x 1 ) for the p-dimensional x 1 -dynamics. Then, the same result holds when substituting the dynamics F 10 (x 1 ) with F 1 (x 1 , c(x 1 )) and F 20 (x 1 ) with F 2 (x 1 , c(x 1 )).

= f 2 (x 1 ) + ug 2 (x 1 ); ẋ1 = f 1 (x 1 ) + ug 1 (x 1 ) (30)
with the origin of the x 1 -subsystem asymptotically stable (there exists a positive definite Lyapunov function V (x 1 ) such that L f 1 V (x 1 ) < 0 for all x 1 = 0) and assume the existence of the function φ 1 (x 1 ) solving L f 1 φ 1 = f 2 with φ 1 (0) = 0. Under piecewise constant control over time intervals of length δ , the sampled-data equivalent model takes the form (2)

x 2k+1 = x 2k + F δ 2 (x 1k , u k ); x 1k+1 = F δ 1 (x 1k , u k ) with F δ 2 (x 1k , u k ) = δ 0 e τ( f 1 +u k g 1 ) ( f 2 + u k g 2 )(x 1 )dτ x 1k F δ 1 (x 1k , u k ) = e δ ( f 1 +u k g 1 ) x 1 x 1k
.

It is a matter of computations to show that the two conditions A1)-A2) hold with Φ 1 (•) = φ 1 (•); which proves the existence of a piecewise constant solution from sampled state measures.

V. EXAMPLE

Let the academic example on the plane

ẋ2 (t) =x 1 (t) + x 3 1 (t); ẋ1 (t) = -x 1 (t) + u(t). (31)
It is easily verified that continuous-time I&I stabilization is achieved by the control u = -Kz 2 (1 + x 2 1 ) with K > 0, z 2 = x 2φ 1 (x 1 ) and φ 1 (x 1 ) = -x 1 -1 3 x 3 1 . Furthermore, the dynamics is zero state detectable with respect to the "output" map z 2 (1 + x 2 1 ). Consider now the sampled-data equivalent to (31) which is exactly computable and given by

x 2k+1 =x 2k + (1 -e -δ )x 1k + x 3 1k δ 0 e -3τ dτ +u k [e -δ + δ -1 + 3x 2 1k δ 0 e -2τ (1 -e -τ )dτ] +3x 1k u 2 k δ 0 e -τ (1 -e -τ ) 2 dτ + u 3 k δ 0 (1 -e -τ ) 3 dτ x 1k+1 =e -δ x 1k + (1 -e -δ )u k . ( 32
)
The origin is still GAS for the discrete-time x 1 -dynamics. Hence, one defines the discrete-time target as ξ k+1 = e -δ ξ k . Then, setting z 2 = x 2φ 1 (x 1 ) with the same φ 1 (x 1 ) as in the continuous-time case, one verifies that z 2k+1 = z 2k under u k ≡ 0. Consequently, the immersion mapping is described, as in the continuous-time case, by π(ξ ) = col(-ξ -1 3 ξ 3 , ξ ). The z 2 -dynamics rewrites in the (F 0 , G) form as

z + 2 = z 2 ∂ z + 2 (u) ∂ u = δ + 3[e δ x + 1 (u) + (1 -e δ )u] 2 δ 0 e -2τ (1 -e -τ )dτ + [e δ x + 1 (u) + (1 -e δ )u] 2 e -2δ (1 -e -δ ) + 6u[δ e δ x + 1 (u) + (1 -e δ )u][ δ 0 e -τ (1 -e -τ ) 2 dτ + e -δ (1 -e -δ ) 3 ] + 3u 2 [ δ 0 (1 -e -τ ) 3 dτ + (1 -e -δ ) 3 3 ] (33)
with z + (u) = z 2k+1 and x + 1 (u) = x 1k+1 . The problem results in defining a digital I&I control law which makes z 2 → 0 as k → ∞ preserving boundedness of the whole state trajectories (32)-(33). According to Section IV.2, one defines the control Lyapunov function V (z 2 ) = z 2 2 and the output

H(x 1 , z 2 , u) = L G z (x 1 ,u) V z (z 2 ).
It is a matter of computation to verify that zero-state detectability is inherited from the continuous-time one. At this point, one notices that in order to find u = -1 0 H(x + 1 (su), z + 2 (su), su)ds one has to solve a fifth degree equation in u, which is hard. Hence, one looks for the bounded I&I negative output feedback defined as in Theorem 4.1. In particular, one computes

G z (F 0 (x 1 ), 0) =δ + 3x 2 1 δ 0 e -2τ (1 -e -τ )dτ + x 2 1 e -2δ (1 -e -δ )
and sets u = -λ (x 1 , z 2 )ε(x 1 )z 2 G z (F 0 (x 1 ), 0) choosing any µ > 0 and ε(x 1 )λ (x 1 , z 2 ) so as to satisfy Theorem 4.1.

A. Simulations

Simulations are referred to the continuous-time I&I feedback (red); its emulation (dotted red) which corresponds to hold constant over the sampling time the continuous-time controller; the dead beat approach (blue) and the bounded I&I output feedback (dotted blue). The invariant set is plotted in black. The bounded feedback is implemented by choosing λ (•) as the upper-bound of (28). Simulations are reported for increasing values of the sampling period δ (0.1, 0.4 and 0.6 seconds). We can see that while the proposed controllers (dead-beat and average-passivity based) succeed in preserving the I&I properties and the stability of the closed-loop equilibrium, the emulated-based control yields degradate performances for δ = 0.4s and instability for δ = 0.6s. We also note that, as the sampling period increases, the sampled-data controllers require a lower control effort than the continuous-time and emulated-based ones.

VI. CONCLUSIONS

With respect to usual forwarding, the proposed I&I approach enables us to relax the knowledge of the x 1 -control Lyapunov function V . Secondly, the bounded design relaxes the necessity to solve an implicit algebraic equality thanks to a suitable dynamic gain. Regarding sampled-data dynamics, it results that I&I stabilization of continuous-time feedforward dynamics implies I&I stabilization of its equivalent sampled-data dynamics in the case of two interconnected blocks.

Remark 4 . 4 :

 44 Consider the continuous-time cascade ẋ2

Fig. 1 .Fig. 2 .

 12 Fig. 1. δ = 0.1 seconds.

Fig. 3 .

 3 Fig. 3. δ = 0.6 seconds.

This work was partially supported by a CNRS-ST2I International Scientific Project -PICS -for cooperation between France and Italy.

S.Monaco is with

A. I&I dead-beat stabilizing control

A first dead-beat solution can be computed by solving in u k for all z 2k , the implicit equality z 2k+1 = 0; i.e.

1 (su k ))ds = 0 so bringing in one step z 2k+1 to zero (equivalently x 2k+1 to Φ(x 1k+1 )). Then u k+i = 0 for i ≥ 1 guarantees z 2k+i+1 = 0 so that the trajectory lays on the stable set.

B. I&I negative output feedback

It is shown in this section that the concept of u-average passivity, introduced by the authors in [START_REF]Nonlinear avarage passivity and stabilizing controller in discretetime[END_REF], can be fruitfully used to get asymptotic stabilization of the origin of the z 2dynamics. With this in mind, let us associate to (20) in the coordinates (x 1 , z 2 ) the output

ant its average [10]

The next result is an immediate consequence of the driftless property of (22).

Proposition 4.2: The feedforward dynamics (18) satisfying A1)-A2) with output (23) is u-average lossless with storage function V z (z 2 ) = 1 2 z 2 2 ; i.e. it satisfies

In the present case, one easily computes with G z (•, u) in ( 21)

From [START_REF]Nonlinear avarage passivity and stabilizing controller in discretetime[END_REF], any feedback law making (23) negative renders the origin of the z 2 -dynamics GAS, provided the so defined output H(x 1 , z 2 , u) is zero state detectable -ZSD -; i.e.

ZSD: no solution of the uncontrolled dynamics (18) can stay in the set

), 0)} other than solutions converging asymptotically to the zero equilibrium.

Proposition 4.3 (I&I Negative u-average output feedback):

Given the feedforward dynamics (18) satisfying A1) A2) with output mapping z + 2 (u)G z (x + 1 (u), u) assumed zero state detectable then, for all (x 1 , z 2 ) ∈ R p+1 , the control law u solution of the algebraic equation

with sufficiently small ε(x 1 ) > 0, ensures global asymptotic stabilization of the origin of (18).

Proof: From the lossless equality (23) rewritten as

and according to Proposition 4.2, the control law solution of the algebraic equality

with sufficiently small ε(x 1 ) > 0 achieves in closed loop

Asymptotic stability to the origin of the z 2 -dynamics follows under the requested zero state detectability condition of the mapping z 2 G z (F 10 (x 1 ), 0). Boundedness of the trajectories is ensured by sufficiently small ε(x 1 ). Then, I&I stabilization of the origin of the feedforward dynamics (18) follows.

C. A constructive bounded solution

The solution proposed in ( 24) is implicitly defined and thus cannot be exactly computed in general. Setting s = 0 in the right hand side of (24), one gets a computable approximation of the solution in the form

). Lemma 4.1 below, recalled from [START_REF] Mazenc | Adding integrations, saturated controls and stabilization for feedforward systems[END_REF], [START_REF] Mazenc | Forwarding in discrete-time nonlinear systems[END_REF] is used to provide a bounded feedback of the form (26) with suitably chosen gain.

Lemma 4.1: [START_REF] Mazenc | Adding integrations, saturated controls and stabilization for feedforward systems[END_REF] Let K(ξ , u) be a continuous function. For any strictly positive real number ū, there exists a function λ (ξ ), as smooth as K(ξ , u) is, such that if u(ξ ) = -λ (ξ )K(ξ , u) then, for all ξ , we have

Moreover, if K(ξ , u) is C 1 , then λ is strictly positive on any compact set.

On the basis of Proposition 4.3 and Lemma 4.1 we can now prove the following result Theorem 4.1: Bounded I&I negative u-average output feedback -Given the feedforward dynamics (18) satisfying A1)-A2) then, for any bound µ > 0, the feedback

where λ (x 1 , z 2 ) is any function that satisfies


F Aman

Michel Vacher
email: michel.vacher@imag.fr

François Portet
email: francois.portet@imag.fr

W Duclot

Benjamin Lecouteux

CirdoX: an On/Off-line Multisource Speech and Sound Analysis Software

Keywords: audio and speech processing, natural language and multimodal interactions, Ambient Assisted Living (AAL)

HAL is

Introduction

Voice User Interfaces (VUI) are becoming more and more popular in ambient intelligence environments whether it is for smartphones or for smart homes (e.g., S.A.H.R.A.1 ). For instance, in domestic environments, VUI recently gained interest in the speech processing community as exemplified by the rising number of smart home projects that consider Automatic Speech Recognition (ASR) in their design [START_REF] Istrate | Embedded implementation of distress situation identification through sound analysis[END_REF][START_REF] Badii | CompanionAble -integrated cognitive assistive & domotic companion robotic systems for ability & security. In 1st Congres of the Société Française des Technologies pour l[END_REF][START_REF] Hamill | Development of an automated speech recognition interface for personal emergency response systems[END_REF][START_REF] Filho | From science fiction to science fact: a smart-house interface using speech technology and a photorealistic avatar[END_REF][START_REF] Gemmeke | Self-taught assistive vocal interfaces: an overview of the aladin project[END_REF][START_REF] Christensen | homeservice: Voice-enabled assistive technology in the home using cloud-based automatic speech recognition[END_REF][START_REF] Cristoforetti | The DIRHA simulated corpus[END_REF][START_REF] Vacher | Evaluation of a context-aware voice interface for Ambient Assisted Living: qualitative user study vs. quantitative system evaluation[END_REF]. However, the speech technology, though mature is still difficult to master for non-expert and to tune for a particular usage (Vacher et al., 2015a). Although many good quality commercial systems were made available during this decade, these are neither open source nor customisable by the user. Moreover, for some of them, respect to user privacy may be questioned. Hence, it becomes important for research in natural language interaction to access opensource and highly configurable real-time VUI software. Many toolboxes and libraries exist for audio processing such as Planet CCRMA at Home (CCRMA, 2016), ChucK [START_REF]ChucK : Strongly-timed, Concurrent, and On-the-fly Music Programming Language[END_REF], etc. But these are dedicated for music processing and not for speech interaction. General purpose recorder and editor such as SoX [START_REF] Sox | SoX -Sound eXchange[END_REF] and Audacity [START_REF] Audacity | [END_REF]) make it possible to acquire audio signals in real-time and edit them, but these tools are mainly designed to perform off-line processing. Regarding speech interaction, there are a high number of ASR (Automatic Speech Recognition) systems2 but again, using them in an on-line hand free setting requires a high expertise in signal and speech processing. Ideally, a complete hand free VUI system would be composed of five main components:

1. an multisource acquisition stage; A few tools have been built to perform some of these stages such as AuditHIS [START_REF] Vacher | Complete Sound and Speech Recognition System for Health Smart Homes: Application to the Recognition of Activities of Daily Living[END_REF] or more recently PATSH [START_REF] Vacher | Experimental Evaluation of Speech Recognition Technologies for Voice-based Home Automation Control in a Smart Home[END_REF]. They were built to be used in continuous sensing uncontrolled ambient system using a pipeline architecture of audio acquisition, sound detection, classification and ASR. However, these systems are not open-source and difficult to access. In this paper, we present an open-source VUI software, called CIR-DOX which performs multisource audio analysis to extract speech occurrences and handle the speech recognition. Moreover, the system performs non-speech audio event recognition (cough, shock, etc.) so that interaction can also be performed using sound. The software is designed to be modular and to process on-line the audio stream (as soon as audio sample chunks arrive). In this way, inclusion of new modules is designed to be easy.

CirdoX: an Open Source Audio Analysis Software

The architecture of the presented system is depicted Figure 1. It has been entirely developed in C and it is currently running on Unix-like POSIX systems only. We plan to diffuse the software under the GNU GPL license. The source code will be available on https://forge.imag.fr/projects/cirdox/.

The backbone of CIRDOX is its decomposition into three processes running in parallel:

1. The Acquisition and Detection Process which handles the continuous audio sensing and the extraction of audio events of sufficient energy from multichannel microphone or record (1 to 8 channels); 2. The Classification Process, which performs first the speech/non-speech discrimination and secondly the identification of each non-speech event. Thirdly, for speech events, a call to an external ASR system is made by sending speech signal and receiving hypothesis;

3. and, the Communication Process, which performs subsequent processing and formatting to communicate, the interpreted audio events to external devices (home automation system, logger, alarm system, etc.).

The CIRDOX application is designed to be modular and configurable, with independent modules so that each module can be chosen from among several plugins correspond-ing to different ways of achieving a task (for instance, classification using either Gaussian Mixture Model or decision trees). The plugins are loaded using dynamic library linking so that any user can build her own plugins to perform some of the module tasks. Table 1 shows the different plugins implemented and that we plan to implement.

The information exchange is mainly performed using a sound object data structure which contains the full information necessary to characterize a detected audio event (time of begin/end, signal to noise ratio) and whose fields are filled in along the processing chain (type of sound: speech or non-speech, ASR hypothesis, every day sound class, etc.). These sound objects are exchanged between processes using IPC (Inter-Process Communication) sig- CIRDOX can be executed in two modes:

1. online, meaning the input sources are the different microphones, 2. off-line, meaning that the input is actually extracted from audio files.

The second mode is useful for the development of new modules in order to test them in reproducible conditions. Moreover, the off-line mode can be executed either in realtime (chunks of data are processed at the sampling frequency of the audio record) or in stream mode (chunks of data are processed as fast as possible).

Acquisition/Detection Process

In online mode, the audio stream is captured by a set of microphones and recorded by the Acquisition module. At this time of writing, this module can be executed through the following plugins:

1. the PortAudio library 3 ;

3 www.portaudio.com is used by many open-source applications including Audacity to access the sound card 2. The Kinect API version 04 ; 3. the driver of a National Instrument card 5 .

The acquisition can also store the full recorded signal for further analysis as well as simulate an acquisition from a audio file for development and debugging purpose. Once the acquisition has reached a certain amount of data, a callback function triggers the Detection module. The Detection module detects the occurrences of audio events (speech or non-speech). It keeps the state of the previous detection between the successive calls by the acquisition module. Rather than multi-threading the detection for each channel, the detection is performed successively on each channel. This is made possible by the synchronous multichannel acquisition. At the time of writing only one plugin as been written to perform the module task. This plugin is based on the change of energy level of the three highest frequency coefficients of the Discrete Wavelet Transform (DWT) in a floating window frame (last 2048 samples without overlapping) [START_REF] Vacher | Sound detection and classification through transient models using wavelet coefficient trees[END_REF]. Each time the energy on a channel goes beyond a self-adaptive threshold, an audio event is detected until the energy decreases below this level for at least an imposed duration. Hence, this plugin adapts itself to changes in background noise level. The module also identifies simultaneous events using simple temporal overlapping 6 . This make it possible to perform multisource processing or data fusion. Too lengthy audio events can also be discarded from the process 7 . During detection, the Signal to Noise Ratio (SNR) of the event is calculated for each channel. Once a detection is complete, event data are saved in an audio file in a local repository. In the case of simultaneous events, the saved data are those from the event with the best SNR. A sound object with a link to the audio file is created and sent by IPC to the next process.

Classification Process

This process is in charge of interpreting the audio events. A first module, called the Discrimination module, determines whether the audio event is a speech or non-speech event. This module can be executed through plugins using Gaussian Mixture Models (GMM) or Decision Tree (DT) (the DT plugin is not yet implemented). Two plugins with GMM are implemented and use the ALIZE library [START_REF] Bonastre | Alize, a free toolkit for speaker recognition[END_REF]: the first uses features in the SPRO3 format and integrates ALIZE API functions, and the second uses features in the SPRO4 format and calls LIA_RAL tools. In both cases, signal discrimination is achieved through a Gaussian Mixture Model (GMM) classifier, trained with a everyday life sound corpus, and a speech corpus recorded in the LIG laboratory [START_REF] Istrate | Information Extraction from Sound for Medical Telemonitoring[END_REF]. Acoustical features can be Linear-Frequency Cepstral Coefficients (LFCC) or Mel-Frequency Cepstral Coefficients (MFCC). The classifier uses 24 Gaussian models, on a logarithmic Mel scale. In the case of non-speech events, the Classification module attributes a sound class (e.g., cough, shock, etc.) to the event. Again, some plugins based on GMMs, Decision Tree (DT) are or will be available. For the moment, only the GMM plugin, with a classifier using the ALIZE API functions, is implemented. In the case of speech, the signal can be sent to an external ASR (Automatic Speech Recognition) system so that a transcription of the speech can be obtained. For the ASR synchronisation task, several plugins can be used, corresponding to different ASR systems: Sphinx3, Google Speech API and KALDI. 6 An event might be captured by several microphones, for instance when they are set in the same room.

7 Depending on the application, sound events might have a minimum and maximum duration. This make it possible to filter out too short events and installed background noise.

Regarding the plugin connecting to Sphinx3, ASR is managed with a Tcl/Tk script that calls Sphinx3 tools8 for signal to MFCC conversion and decoding. The synchronization between the Classification process and the Tcl/Tk script is done by lock files. A text file with the hypothesis is generated by Sphinx3. Regarding the plugin connecting to Google Speech API, the wav file is converted into the FLAC format and sent to the Google Speech API server9 through a POST request. The Google server return a JSON response with the hypothesis. Thanks to the plugin connecting to KALDI, the wav file is sent to an ASR server, based on the Kaldi toolkit [START_REF] Povey | The kaldi speech recognition toolkit[END_REF] and the GStreamer framework, and implemented in Python10 . Like the Google Speech API plugin, a JSON response with the hypothesis is returned. Sound classification results or ASR hypothesis are included in the sound object for processing by the next module.

Communication Process

The communication process handles the interaction with the exterior. It also logs all the sound events. This part is highly dependant on the application to communicate with. To take this dependence into account a filter module can be activated to communicate only the event to communicate to the exterior (e.g., only speech of a certain kind) and/or to perform more processing on the event. The formatting module takes care only of the communication protocol.

CirdoX in Ambient intelligence Setting

CIRDOX has been used at different occasions to capture and analyse audio events in different smart spaces including the smart class room of the MSTIC Fab-Lab and the Domus smart home of the LIG. We present below the results of a setting for voice command recognition.

Voice command recognition for Home Automation

Figure 2 shows the connections of CIRDOX with the input (wireless microphones set in rooms of a home) and the output (home automation system) that has be employed both in the CIRDO project [START_REF] Bouakaz | CIRDO: Smart companion for helping elderly to live at home for longer[END_REF] and the sweet-Home project [START_REF] Vacher | Evaluation of a context-aware voice interface for Ambient Assisted Living: qualitative user study vs. quantitative system evaluation[END_REF]. The aim is to constantly capture the audio signal to seek for specific voice commands and send them to the controller of the home automation system in order to perform the commands when applicable. In the case of a generic smart home equipped with a home automation system, actuators are used to activate lighting, stores, multimedia devices, etc. In the particular case of calls for help, the only action is alarm sending to a specialized service or to a relative. In all cases, the speaker is distant from the microphone (i.e., no worn microphone) and the interaction is hand-free (no push button to start the recording).

To address voice command recognition realistically, 3 kinds of speech were considered:

• vocal command (e.g., "home turn on the light"),

• distress calls (e.g., "home call for help")

• or colloquial speech (e.g., "I ate spinach at lunch").

In this last case, no action must be operated and the audio event must be discarded. Therefore, it is necessary to check if a speech audio event corresponds to a voice command or a distress call. To do so, the exhaustive list of all possible voice commands and distress calls was generated wrt a predefined grammar. Since a word-by-word comparison would be too crude, each speech transcription hypothesis was first phonetized and then a distance was computed against the list of phonetized voice commands. Formally, for each phonetized ASR output T , every voice commands H is aligned to T using Levenshtein distance [START_REF] Levenshtein | Binary codes capable of correcting deletions, insertions and reversals[END_REF]. The deletion, insertion and substitution costs were computed empirically while the cumulative distance γ(i, j) between H j and T i is given by Equation 1.

γ(i, j) = d(T i , H j )+ min{γ(i -1, j -1), γ(i -1, j), γ(i, j -1)} (1)
The sentence with the aligned symbols score is then selected for decision according a detection threshold. This approach takes into account some recognition errors such as word endings or light variations. Moreover, in a lot of cases, a miss-decoded word is phonetically close to the true one (due to the close pronunciation). At the end of the process, the sound object is received by the Data Formatting module that sends a socket containing all the relevant data to the controller system.

Experiment in Living lab

CIRDOX has been used in the framework of the CIRDO project with participants who played scenarios in the Living lab DOMUS, including falls down on the floor and calls for help uttered in French. Since it is difficult to find elderly people capable and accepting to play such scenarios, some recruited people were under 60. Overall 17 French participants were recruited (9 men and 8 women) and 4 were elderly people [START_REF] Vacher | The CIRDO corpus: comprehensive audio/video database of domestic falls of elderly people[END_REF].

The younger participants were instructed to wear an equipment i.e. old age simulator. Old age simulator hampered mobility, reduced vision and hearing. Figure 3 shows a young participant during the training phase before playing a fall scenario on the carpet. Scenarios included call for help sentences.

CIRDOX Configuration

CIRDOX was used in order to process data and to extract online the sentences uttered during the experiment. GMM models were trained using the ALIZE software in order to adapt the Speech/Non-speech Discrimination module. Training data were made of distant speech records from SWEET-HOME corpus [START_REF] Vacher | The Sweet-Home speech and multimodal corpus for home automation interaction[END_REF] recorded in the Living lab DOMUS: 9,147 speech files (home automation orders or distress calls) were used to train the speech model (i.e. 4 hours 53 minutes), 13,448 sound files (i.e. 5 hours 10 minutes) were used to train the non-speech model. ASR was operated using Sphinx3 with an acoustic model trained on the BREF120 corpus [START_REF] Lamel | BREF, a large vocabulary spoken corpus for french[END_REF] and adapted to elderly voice, to distant speech and to expressive speech conditions. BREF120 is a large corpus of oral read texts, containing more than 100 hours of speech produced in French by 120 speakers (65 women and 55 men), all recorded texts were extracted from the newspaper Le Monde. The adaptation was obtained thanks to the MLLR method with the "User Specific Interaction" subset of the SWEET-HOME corpus [START_REF] Vacher | The Sweet-Home speech and multimodal corpus for home automation interaction[END_REF] which is made of 337 sentences in French (9min 30s of signal) uttered in the DOMUS flat by elderly or visually impaired people (average age: 72 years), and with 742 expressive speech sen-Type of sound Classified Classified as speech as non-speech Speech 237 85 Non-speech 7 1621

Table 3: Discrimination between "Speech" and "Nonspeech".

tences of the "Voix détresse" corpus recorded by the LIG laboratory.

A language model specific to sentences to recognize was estimated from distress sentences and calls to caregivers of AD80 corpus (Vacher et al., 2015a). This is a trigram model containing 88 unigrams, 193 bigrams and 223 trigrams. A general language model has been learned from the French Gigaword corpus (Linguistic Data Consortium, 2011) which is a collection of newspaper articles that have been acquired over several years by the Linguistic Data Consortium (LDC) from the University of Pennsylvania. This model is unigram and contains 11018 words. The final model is the combination of the specific model with the general model, giving more weight (90%) on the probabilities of specific model.

Audio Analysis Results

The Detection module extracted 1,950 sound events. The confusion matrix of the results of the Speech/Non-speech Discrimination module are given in Table This module classified 244 of the events as speech, 237 were really speech, 7 were non-speech events classified as speech. 85 speech events were badly classified as sound. At the output of the Speech/Non-speech Discrimination module, the specificity is high, Sp = 0.996, this parameter represents the classification rate of non-speech events classified as non-speech. However, the sensitivity is low, Se = 0.736, this parameter represents the classification rate of speech events as speech.

According to results of previous stages of CIRDOX, 244 audio events were decoded, 204 are distress calls, 33 colloquial sentences and 7 non-speech signal. Table 4 presents the decoding results. WER is very high for colloquial sentences, this is due to the language model which is non adapted to this task and corresponds to a desired effect: this kind of sentences is not to be well recognized for privacy reason. Regarding distress calls, the Word Error Rate (WER) is quite high with 38.52%. This result is far from perfect but they were obtained under harsh conditions. Indeed, the participants played scenarios which included falls on the floor: they generated a lot of noise sounds which were often mixed with speech, and they were far from the microphone.

At the output of the filtering stage, 29 colloquial sentences and sounds were rejected (True Negative), 35 distress call were rejected by mistake (False Negative), 11 colloquial sentences and sounds were kept by mistake (False Positive), and 169 distress calls were recognized (True Positive). Among the 204 distress calls, 82.84% of the calls are well detected. The presented ASR results were obtained using Sphinx3. We obtained better results using a more sophisticated ASR approach by modellling acoustic using SGMM and on manually annotated data from the same corpus. Acoustic models were based on subspace GMM (sGMM) with the KALDI toolkit [START_REF] Vacher | Recognition of Distress Calls in Distant Speech Setting: a Preliminary Experiment in a Smart Home[END_REF].

In the experiments, the acoustic models were contextdependent classical three-state left-right HMMs. Acoustic features were based on Mel-frequency cepstral coefficients, 13 MFCC-features coefficients were first extracted and then expanded with delta and double delta features and energy (40 features). Acoustic models were composed of 11,000 context-dependent states and 150,000 Gaussians. The state tying is performed using a decision tree based on a treeclustering of the phones.

In addition, off-line fMLLR linear transformation acoustic adaptation was performed. The acoustic models were trained on 500 hours of transcribed French speech composed of the ESTER 1&2 (broadcast news and conversational speech recorded on the radio) and REPERE (TV news and talk-shows) challenges as well as from 7 hours of transcribed French speech of the SWEET-HOME corpus [START_REF] Vacher | Evaluation of a context-aware voice interface for Ambient Assisted Living: qualitative user study vs. quantitative system evaluation[END_REF] eliciting a distress emotion. WER was 34% and 72% of the calls were detected using Kaldi and SGMM instead of 61% using Sphinx3 .

Conclusion

Audioprocessing in the Smart Homes, despite recent developments, have led to few experimentations in daily living conditions as far as audio analysis is concerned. This paper presents a modular audio analysis software able to analyse online the audio environment from multichannel, to extract the uttered sentences and then to process them thanks to an ASR stage. Moreover this system perfoms non-speech audio event identification; in this case, specific models must be trained. The software is designed to be modular and to process on-line the audio stream. This software was used for call for help recognition in a smart home. The records were made in distant speech conditions, moreover the conditions were very harsh because the scenarios included falls of the person on the carpet. Audio events were correctly processed along the different stages. Results must be improved by using new models. We plan to diffuse the software under the GNU GPL license.

  2. a signal enhancement stage; speech_recognition_software 3. a Voice Activity Detection (VAD) stage; 4. an Automatic Speech Recognition (ASR) stage; 5. a Natural Language Understanding (NLU) stage, when required; 6. a decision stage; 7. and, a communication stage.

Figure 1 :

 1 Figure 1: Architecture of the CirdoX software.

Figure 2 :

 2 Figure 2: Home Automation system using CirdoX.

Figure 3 :

 3 Figure 3: A participant during the training phase of a fall scenario (Bouakaz et al., 2014).

Figure 4 :

 4 Figure 4: Global view of result experimentation in the DOMUS Living lab (Distress calls: Colloquial sentence: C, Non-Speech: NS).

Table 1 :

 1 Functionalities and Implementations of CIRDOX modules.

	Process	Functionality & modules	Implemented Library & model
	Acquisition & Detection
		Acquisition: NI on-line	Not yet	NI-DAQmx driver
		Acquisition: NI off-line	Not yet	NI-DAQmx driver
		Acquisition: Portaudio on-line	Yes	libportaudio2 v19
				portaudio19dev v19
		Acquisition: Portaudio off-line	Yes	libportaudio2 v19
				portaudio19dev v19
		Acquisition: Kinect on-line	Yes	libfreenect v0.2
		Acquisition: Kinect off-line	Not yet	-
		Detection: DWT	Yes	(Vacher et al., 2004)
	Classification
		Speech/non-speech discrimination: GMM	Yes	LIA_RAL v3.0-to be adapted
				to home audio environment
		Speech/non-speech discrimination: GMM	Yes	ALIZE API v3.0 -to be adapted
				to home audio environment
		Non-speech classification: GMM	Yes	ALIZE API v3.0 -to be adapted
				to home audio environment
		Non-speech classification: DT	Not yet	-
		ASR synchronisation: Sphinx3	Yes	to be adapted to the application
		ASR synchronisation: Google API	Yes	different languages supported
		ASR synchronisation: KALDI python server Yes	to be adapted to the application
	Communication
		Filtering: Levenshtein Distance	Yes	list of sentences to furnish
		Formatting	Yes	to be adapted to the application

nals and a common local repository which contain the raw audio data (written by the Acquisition and Detection process and read by the Classification).

  Je peux pas me relever Appelle quelqu'un e-lio Qu'est-ce qu'il m'arrive e-lio, appelle quelqu'un Aïe ! J'ai mal e-lio appelle ma fille Oh là ! Je saigne ! Je me suis blessé e-lio appelle le SAMU

	Distress Sentence	AD80 Sentence
	Aïe aïe aïe	Aidez-moi
	Oh là	Au secours
	Merde	e-lio, appelle du secours
	Je suis tombé	e-lio appelle ma fille

Table 2 gives some examples of these sentences. The size of audio records is 645.54 seconds.

Table 2 :

 2 Some examples of French sentences identified for inclusion in the scenarios

Table 5

 5

	Type	WER (%)
	Calls for help	38.52
	Other	115.31
	(colloquial and non-speech)

summarize the performances of this stage. The False Alarm Rate (FAR) is 6.11%.

Table 4 :

 4 WER for the DOMUS recorded data.

	Performance measure	Value (%)
	Se	82.84%
	Sp	72.50%
	FAR	6.11%

Table 5 :

 5 Sensitivity, specificity and false alarm rate of the filtering stage from automatically detected speech Results obtained by the Filtering stage and the previous stages are summarized Figure4. Among the 1950 audio events, 277 are distress calls, 45 are colloquial sentences and 1628 are non-speech. Considering the global performances, i.e. all the steps (discrimination, ASR and filtering modules), 61% of the distress calls uttered by the speakers lead to send an alert to the concerned service. This value is still relatively too low to consider a real use in realworld applications. One explanation could be the presence of noise in the speech. Thus 95 speech events were noisy, in most cases due to the superposition of the noise made by the person during the fall to the speech. Only 43.2% of speech events were classified as speech by the Discrimination module when they were noisy versus 86.3% for clean speech. It is also not surprising, all depend on the proportion of noise related to speech in term of energy or time.

	Moreover, for clean speech classified as speech, WER is
	33.28%, versus 77.91% for noisy speech.

  which consists of records of 60 speakers interacting in the smart home and from 28 minutes of the Voix-détresse corpus recorded by the LIG laboratory[START_REF] Aman | Analysing the Performance of Automatic Speech Recognition for Ageing Voice: Does it Correlate with Dependency Level? In 4th Workshop on Speech and Language Processing for Assistive Technologies[END_REF] which is made of records of speakers

		Audio: 1950 events	D: 277, C: 45, NS: 1628
	Speech/Non-speech discrimination			Classified as?
	D: 204, C: 33, NS: 7	Speech: 244		Non-speech: 1706	D: 73, C: 12, NS: 1621
	Filter: n(M1)/n(M2)
	169	29		73	1621 + 12
		11	35
	True Positive	False Positive	False Negative	True Negative

http://sarah.encausse.net

https://en.wikipedia.org/wiki/List_of_

http://openkinect.org/

http://www.ni.com/product-documentation/ 11787/en/

http ://cmusphinx.sourceforge.net

http://www.google.com/speech-api/v1/recognize

https://github.com/alumae/kaldi-gstreamer-server

Acknowledgements

This work was supported by the French funding agencies ANR and CNSA through CIRDO project (ANR-2010-TECS-012).