Rene Abib
email: renee.abib@univ-rouen.fr

GEODESICS IN NILPOTENT LIE GROUPS

Keywords: geometry, nonlinear control systems, distributions

We study the sub-Riemannian geodesic problem in the 2n + 1-Heisenberg group. An independent set of 2n left-invariant vector fields determines a SR structure. We are interested in sub-Riemannian length minimizers for this structure. We use the hamiltonian formalism and apply Pontryagin maximum principle to write the equations for the extremals and describe some properties of the extremal curves. We consider geodesics in the group with a left-invariant Riemannian metric. We obtain their equations and prove that the set of directions of all rays issuing form zero is the sphere of unit horizontal directions.

INTRODUCTION

Let L be the 2n+1-dimensional nilpotent Lie algebra with the following commutation rules in some basis X 1 , Y 1 , X 2 , Y 2 , . . . , X n , Y n , T :

[X i , Y j ] = 2δ ij T, [X i , X j ] = [Y i , Y j ] = ad T = 0,
and G be the connected, simply connected Lie group with the Lie algebra L. We consider X i , Y i , T as left-invariant vector fields on G. The set {X 1 , Y 1 , • • • , X n , Y n } determines a leftinvariant sub-Riemannian structure on G:

D = span (X 1 , Y 1 , • • • , X n , Y n ), X i , Y j = X i , X j = Y i , Y j = δ ij i, j = 1, 2, • • • , n,
called the flat (2n, 2n + 1) sub-Riemannian structure. Such structure is unique, up to isomorphism of Lie groups.

We are interested in sub-Riemannian length minimizers for this sub-Riemannian structure, i.e., in solutions to the fol-lowing optimal control problem:

γ = n i=1 (u i X i (γ) + v i Y i (γ)) (1) γ ∈ G, u i , v i ∈ R, γ(0) = q 0 , γ(T ) = q 1 fixed, (2)
l = T 0 n i=1 (u 2 i + v 2 i ) → min . (3)
Problem (1)-(3) was considered by [START_REF] Brockett | Control theory and singular Riemannian geometry[END_REF], and [START_REF] Liu | Shortest paths for sub-Riemannian metrics on rank-two distributions[END_REF] for n = 1. G is a 2-step nilpotent Lie group and

L = D ⊕ [L, L].
In this note we consider also geodesic lines in the Lie group G with a leftinvariant Riemannian metric.

The flat (2n, 2n + 1) sub-Riemannian structure gives a local nilpotent approximation for an arbitrary sub-Riemannian structure with growth vector (2n, 2n + 1), see [START_REF] Montgomery | A tour of Sub-Riemannian geometries, their geodesics and applications[END_REF]. The dynamics of a classical electric charge in the plan un-der the influence of a perpendicular magnetic field can be described by means of flat sub-Riemannian structure.

In general, if G is a 2-step nilpotent Lie group and (D, , ) is a left-invariant sub-Riemannian structure such that L = D ⊕ [L, L], any minimizer is normal in some sub-group of G and hence any minimizer is smooth, see [START_REF] Abib | Sub-Riemannian geodesics on Lie groups[END_REF]. This extends known results on the so-called Gaveau-Brockett problem, Brockett (1981) and [START_REF] Gaveau | Principe de moindre action, propagation de la chaleur et estimes sous elliptiques sur certains groupes nilpotents[END_REF].

MODEL

We choose the following model for the flat (2n, 2n + 1) sub-Riemannian structure:

G = R 2n+1 with coordinates q = (x, z) where x = (x 1 , y 1 , . . . . . . , x n , y n ), z ∈ R, and group law defined as follows:

(x, z)(x ′ , z ′ ) = (x + x ′ , z + z ′ + + n i=1 x i y ′ i - n i=1 x ′ i y i ) X i = ∂ ∂x i -y i ∂ ∂z Y i = ∂ ∂y i + x i ∂ ∂z D = Ker θ θ = dz + n i=1 (y i dx i -x i dy i )
Control system (1) has the full rank and the state space G is connected, thus the systems is globally controllable on G. A standard existence theorem from optimal control theory, see [START_REF] Liu | Shortest paths for sub-Riemannian metrics on rank-two distributions[END_REF], implies that any point can be connected with q 0 ∈ G by an admissible length minimizer.

EXTREMALS

A admissible curve for D is an absolutely continuous curve on G which is tangent to D almost everywhere; applying the Pontryagin maximum principle each minimizer parametrized by arclength is a normal extremal or abnormal extremal, see [START_REF] Liu | Shortest paths for sub-Riemannian metrics on rank-two distributions[END_REF].

Abnormal extremals

An abnormal extremal is a admissible curve for D which is the projection onto G of characteristic curve (in the symplectic sense) of the annihilator D ⊥ of D in T * G. Introduce the linear Hamiltonians corresponding to the basis fields:

h i (λ) = λ(X i ), k i (λ) = λ(Y i ), h(λ) = λ(T ), λ ∈ T * G. D ⊥ is defined by equations h i = k i = 0, i = 1, • • • , n. Assume Γ(t) ∈ T * G -{0} caracteristic curve of D ⊥ , then Γ(t) ∈ span { h i (Γ(t)), k i (Γ(t) ; i = 1, 2, • • • , n} and h i (Γ(t)) = k i (Γ(t)) = 0.
Abnormal extremals are exactly the constant curves.

Normal extremals

Normal biextremals are trajectories of the Hamiltonian system

λ = H(λ), λ ∈ T * G (4) with the sub-Riemannian Hamiltonian H = 1 2 n i=1 (k 2 i + k 2 i ). In the coordinates (q, λ) = (x 1 , • • • , x n , y 1 , • • • , y n , z, ξ 1 , • • • , ξ n , η 1 , • • • , η n , ζ) of T * G system (4) reads ẋi = ξ i + y i ζ (5) ẏi = η i -x i ζ (6) ż = n i=1 ( ẋi y i -ẏi x i ) (7) ξi = ẏi ζ (8) ηi = ẋi ζ (9) ζ = 0 (10)
If the constant ζ is 0, the curve (x(t), y(t)) is the straight line, along which the motion takes place with velocity constant. The function z(t) is then computed by solving (7) and the curve (x(t),

y(t), z(t)) is a straight line in G. If ζ is non zero constant, the equations (8)-(9) imply ẍi -2 ẏi ζ = 0, ÿi + 2 ẋi ζ = 0. So ẋi = A i cos 2ζt + B i sin 2ζt, ẏi = -A i sin 2ζt + B i cos 2ζt, A i , B i
are constant. The condition that our curve q(t) = (x(t), y(t), z(t)) is parametrized by arc-length then says that

n i=1 (A 2 i + B 2 i ) = 1. We obtain x i (t) = r i 2ζ (cos θ i -cos(2ζt + θ i )), y i (t) = r i 2ζ (sin(2ζt + θ i ) -sin θ i ) with A i = r i sin θ i , B i = r i cos θ i . The equation (7) imply z = ζ n i=1 d dt (x 2 i + y 2 i ) and then z(t) = n i=1 r 2 i ( t 2ζ - sin 2ζt (2ζ) 2 ).
The Pontryagin maximum principle, applied to this case, implies that all the minimizers curve parametrized by arc-length are normal extremal q(t).

CONTACT TRANSFORMATIONS

For n = 1, the Lie algebra of symmetries of the flat (2,3) distribution D is parametrized by arbitrary smooth functions of three variables and the Lie algebra of symmetries of the flat (2,3) sub-Riemannian structure is 4-dimensional, see [START_REF] Sachkov | Symmetries of flat rank two distributions and sub-Riemannian structures[END_REF].

For n ≥ 1, let C be the sheaf of germs of infinitesimal automorphisms X such that L X θ = f θ, the function f depending on X. C is the sheaf of contact transformations.

Let X ∈ C vanishing at 0, and φ t be the one parameter group generated by X. Then φ t (0) = 0 and φ * t θ = f t θ. Let V 1 be the subspace of L defined by θ 0 = 0. Then φ t leaves V 1 invariant. Furthermore φ * t dθ = df t ∧ θ + f t dθ. This implies that the form dθ restricted to V 1 is preserved up to scalar factor. If we pass to the infinitesimal we see that the linear isotropy algebra of C is a subalgebra of the algebra g defined as follows:

let u ∈ L not lying in V 1 and g = K + N + M where K is the set of all A ∈ End(L) with Au = 0 and A/V 1 ∈ sp(V 1 ), N is the of all B of the form Bv 1 = 0 for v 1 ∈ V 1 and Bu ∈ V 1 , M is the set of all multiples of the linear trans- formation C where Cv 1 = v 1 for v 1 ∈ V 1 and Cu = 2u.
The space g (1) ⊂ Hom (g, L) called the first prolongation of g, is the set of all T ∈ Hom (g, L) which satisfy (T (ω))(v) = (T (v))(ω) for all ω, v ∈ L. In the case, g (1) = sp (V 1 ) (1) + g and g (k) = sp (V 1 ) (k) + g (k-1) for all k. The algebra sp (V 1 ) is of infinite type because sp (V 1 ) (k) can be identified with S k+2 (V * 1 ); then g is of infinite type.

Corresponding to the coordinates z, x 1 , • • • , x n , y 1 , • • • , y n , let us choose the basis u = e 0 , e 1 , • • • , e n , e n+1 , • • • , e 2 of L. Then it is easy to verify that the following vector fields are all infinitesimal contact transformations (where

θ = dz + n i=1 (y i dx i -x i dy i )): α) ∂ ∂x i + y i ∂ ∂z , ∂ ∂y i -x i ∂ ∂z , ∂ ∂z β) n i,j=1 A ij x i ∂ ∂x j + A i+n,j y i ∂ ∂x j + A i,j+n x i ∂ ∂y j +A i+n,j+n y i ∂ ∂y j where 2n r,s=1 A rs e s ⊗ e * r ∈ sp(V 1 ), 2z ∂ ∂z + n j=1 x j ∂ ∂x j + n j=1 y j ∂ ∂y j , z ∂ ∂x i + y i n j=1 x j ∂ ∂x j + y j ∂ ∂y j + z ∂ ∂z , z ∂ ∂y i + x i n j=1 x j ∂ ∂x j + y j ∂ ∂y j + z ∂ ∂z γ) z n j=1 x j ∂ ∂x j + y j ∂ ∂y j + z ∂ ∂z ¿From α) we see that C is transitive. From β)
we see that the linear isotropy algebra of C is indeed g. A direct computation shows that the vector fields α), β) and γ) form a Lie algebra. C is an non-flat, infinite, transitive LAS (Lie algebra sheaf).

RIEMANNIAN CASE

We define the left-invariant metric on G by taking X i , Y i , T as the orthonormal frame.

Let ∇ the Riemannian connection of leftinvariant metric an {e

h , h = 1, 2, • • • , 2n + 1} = {X 1 , Y 1 , • • • , X n , Y n , T }. Then ∇ ei e j = 1 2 2n k=1 (c k ij + c j ki -c i jk )e k where [e i , e j ] = 2n k=1 c k ij e k .
For G the matrice (∇ ej ei ) i,j is:

         0 T 0 • • • 0 -Y 1 -T 0 0 • • • 0 X 1 . . . . . . . . . 0 • • • • • • 0 T -Y n 0 • • • 0 -T 0 X n -Y 1 X 1 • • • -Y n X n 0          (11) Let ċ(t) = n i=1 (u i (t)X i (t) + v i (t)Y i (t)) + γ(t)T geodesics issuing from 0 ∈ G. Then ∇ ċ(t) ċ(t) = 0 and the table (11) give ui + 2γv i = 0, vi -2γu i = 0, γ = 0.
Therefore, because the parameter t is natural we have n i=1 (u 2 i + v 2 i ) + γ 2 = 1 and we could take γ(t) = γ where the constant γ ∈ [-1, 1] is the cosine of the angle between ċ(0) and the T -axe.

For γ = 0, u i (t) = r i cos(2γt + θ i ) and v i (t) = r i cos(2γt + θ i ) where r 2 i = u 2 i + v 2 i .

In coordinates (x 1 , y 1 , • • • , x n , y n , z) the equations for geodesics c(t) = (x 1 (t), • • • , y n (t), z(t)) are:

ẋi = u i , ẏi = v i , ż = γ + n i=1 v i (t)x i (t) - n i=1 u i (t)y i (t) (12)
We have

x i (t) = - r i 2γ cos 2γt + a i y i (t) = r i 2γ sin 2γt + b i z(t) = 3γ 2 -1 2γ t + b i 2γ r i cos 2γt + a i 2γ r i sin 2γt + c
for some numbers a i , b i , c which could be defined from the initial condition c(0) = 0. If γ = 0, then they are horizontal and satisfy to the following:

x i (t) = u i (0)t, y i (t) = v i (0)t, z(t) = 0
The non-horizontal geodesics (i.e., with γ = 0) are not rays. Indeed. Let |γ| = 1. For every geodesic c(t) their x i , y i coordinates are periodic functions. Their projection on z = 0 hyperspace is a circle. Take two points on c(t) with difference of their parameters equals π γ . Then, besides the geodesics c(t) their is also the vertical geodesic z(s) = (0, • • • , 0, s) connecting them. Its length is equal to the difference of z coordinates of considering points. But, as easy to check this is always strictly less than π γ , i.e., the length of the interval of the c(t) between them, because t is natural one by definition. Therefore, every geodesic c(t) with |γ| = 1 is not minimal. Vertical geodesic z(s) (with |γ| = 1) also is not minimal, because that the vector field Z(s) = X 1 cos s + Y 1 sin is parallel along it and the sectional curvatures of G in 2-dimensional direction generated by Z(s) and ż(s) equals 1. This means that the index of every interval of z(s) with a length greater than π is positive, and z is not minimal. Now let us check that all horizontal geodesics (with γ = 0) are rays. Suppose that the horizontal geodesic a(t) intersects in some point with another geodesic c(s) issuing the same point 0. Because horizontal geodesics do not intersect, c(s) is non-horizontal one. For c(t) from ( 12) we see that the length l(c) of the projection curve c(t) = (c 1 (t), • • • , c 2n (t), 0), measured in euclidean coordinates, is strictly less than t, i.e., the length of the curve c in the left-invariant metric. Than in the submanifold z = 0 considered as euclidean space the length of every chord, connecting some points of c is less than l(c). If the geodesic c issuing from 0 intersects some horizontal geodesic a(t) also issuing from 0, then a is the chord of the corresponding c. For horizontal geodesics their lengths in the left-invariant metric coincide with the usual euclidean length in z = 0. The argument above means that the length of a is strictly less then that of c. Therefore, we have the following statement: the set of directions of all rays issuing from 0 in the group G is the sphere of unit horizontal directions in the point 0 i.e., {w = (w 1 , • • • , w 2n , 0); w = 1}.


M Choo-Kun
email: marlene.choo-kun@suez-env.com

D Conteau

S Martin Ruel

P Camacho

R Gourdon

P Buffière

Comparison of the Anaerobic Digestion of Biosorption i.e. A Sludge, Conventional Primary Sludge and Conventional Mixed Sludge

Keywords: Anaerobic Digestion, Biosorption sludge, A sludge, Primary sludge, Mixed sludge, Energy selfsufficiency

 with its two stage of activated sludge leads to reduce the energy demand from the aeration on the water line and increase the methane production on the sludge line which are the key elements to reach energy self-sufficiency for wastewater treatment plants. The present work studies the anaerobic digestion of Biosorption i.e. A sludge from the first stage of the AB process and compares it to the digestion of conventional primary and mixed sludge. The operation of 5-Litres digesters for 100 days shows that A sludge and conventional primary sludge are equivalent in terms of sludge characterisation and behaviour through the anaerobic digestion process. They both produce about 27% more methane per VS fed than the conventional mixed sludge. Considering that A sludge represents a larger portion (85%) of the mixed AB sludge that the primary sludge (65%) does for conventional mixed sludge, the AB process does have a higher potential of energy production than conventional wastewater treatment systems.

INTRODUCTION

The AB (Adsorption-Biosorption) Process for wastewater treatment was developed by [START_REF] Böhnke | Das Adsorptions-Belebungsverfahren[END_REF]. It consists in a two-stage activated sludge: A and B stages, each composed of a bioreactor and a clarifier. The A stage hosts a highly efficient biosorption process for early entrapment of the organic matter, which can be directly transferred to the digesters, maximizing the recovery of the biodegradable carbon contained in the wastewater. Around the world, a few wastewater treatment plants use this process and show a much more efficient energy balance, sometimes even positive, than conventional plants also equipped with digesters. The present study compares the anaerobic digestion process of Biosorption i.e. A sludge with conventional primary sludge and conventional mixed sludge, within semi-continuous 5-Liter reactors under mesophilic conditions (37°C).

MATERIALS AND METHODS

Experimental Setup

Three double jacketed glass reactors of 5 Litres working volume are equipped with a water bath circulator controlling the temperature at 37°C. They are fed on a daily basis thanks to a peristaltic pump. The conventional primary and mixed sludge are sampled from the conventional Water Reclamation Plant of Kranji in Singapore. The A sludge comes from a 50m 3 /d AB process pilot plant on the same WWTP and treating the same raw wastewater. The digesters were inoculated with the digestate of a 3m 3 pilot plant digester at mesophilic conditions fed with Kranji WRP mixed sludge.

The operating conditions are the same for the three reactors with a Sludge Retention Time of 16 days and an Organic Loading Rate of 1.5 kgVS/m 3 /d. The three different feed sludge are thickened through centrifugation to meet the common OLR.

Follow-up and measurement

The biogas production is monitored continuously with milligas counters (Ritter) and the biogas composition is analysed twice a week with a 490 Micro Gas Chromatograph (Agilent). Total Solids (TS) content is calculated as the amount of solids remaining after oven-drying the samples at 105°C during 24 h and Volatile Solids (VS) content is calculated as the loss of weight after igniting the oven-dried residue at 550°C for 5 h. Total COD and soluble COD concentrations, Volatile Fatty Acids (VFAs), Total Alkalinity (TA), Total Ammonia Nitrogen (TAN) are measured with kits (Hach). All these parameters are analysed from both feed and digested sludge three times a week. Biochemical Methane Potential tests were run to characterise the bioreactivity of the three substrates.

RESULTS AND DISCUSSION

Characterisation of the three feed sludge

The table 1 presents a characterisation of the A sludge from the 50m 3 /d AB process pilot plant, the conventional primary and the mixed sludge from Kranji Water Reclamation Plant.

Table1: Characterisation of the A sludge, Primary sludge and Mixed sludge.

A Sludge Reactor -CH4 / VM in

The A sludge and Primary sludge appear to be quite equivalent in terms of behaviour during anaerobic digestion. The Biochemical Methane Potential tests show about the same potential for this two sludge, with a slightly lower potential for the primary, which can be attributed to operational issues of the plant during the sampling (longer SRT of the primary clarifier affecting significantly the quality of the sludge as the clarifier temperature reaches 30°C). The Methane Potential of A sludge is 335 L STP CH 4 / kg VS introduced . The kinetics of the A sludge BMP test is slightly higher than the primary sludge test.

Start-up and steady state operation

The acclimation period lasted for 3 SRTs with a first phase of 1 SRT at half Loading Rate. The methane production steadily increased without any peaks of VFA (VFA<500mg/L, VFA/Alkalinity<0.15) nor Soluble COD (<2g/L). At 2.5 SRTs, the steady state operation was reached for all the reactors with a constant methane production, VS and COD total content.

The following results were obtained during the steady state operation. The primary sludge reactor and the A sludge reactor present equivalent overall performances, both in terms of removal rates and specific methane production. The specific methane production reaches 290 L STP CH 4 /kg VS in for the A sludge digester, and 300 L STP CH 4 /kg for the primary sludge digester.

Removal rates and specific methane production

The VS, COD removal rates and the specific methane production of the anaerobic digestion of A sludge are respectively 33%, 32% and 27% higher than for mixed sludge.

All the results are summarized in the table 2.

Overall AB process system

The mixed sludge produced by the AB process contents 85% of A sludge and 15% of B sludge. A conventional activated sludge wastewater treatment line produces a mixed sludge composed of 65% primary sludge and 35% of secondary sludge. It is also believed and to be proven on future experiments that the AB process produces more sludge than the conventional system. Adding to these facts that the AB sludge has a higher methane production than conventional mixed sludge, the whole AB process should enable to produce more methane, and thus energy per m 3 of wastewater treated, than a conventional wastewater treatment system.

CONCLUSION

-The comparison of the anaerobic digestion of Biosorption i.e. A sludge, conventional primary sludge and conventional mixed sludge was studied with semi-continuous 5-Litres digesters at mesophilic conditions (37°C). -The anaerobic digestion of A sludge is quite equivalent to the anaerobic digestion of conventional primary sludge, with equivalent or slightly better removal rates and kinetics (290 NL CH 4 /VS in ).

-Future work will aim at demonstrating that anaerobic digesters treating sludge from an AB process produce more energy than the digesters associated with a conventional activated sludge process.

Figure 1 :

 1 Figure 1: Reactor set-up.

Figure 2 :

 2 Figure 2: Cumulated methane production over VS fed for the three reactors.

Table 2 :

 2 Overall performances of the three reactors

		A Sludge	Primary Sludge	Mixed Sludge
		Reactor	Reactor	Reactor
	VS Removal	48,5%	45,3%	36,4%
	TS Removal	43,0%	40,6%	30,6%
	COD Removal	42,8%	47,2%	32,5%
	L biogas / kg VSin	470	494	380
	L CH4 / kg VS in	290	301	229
	L CH4 / kg COD in	150	156	112
	L CH4 / kg COD destroyed	345	343


Mandimby Ranaivo Rakotondravelona
email: mandimby.ranaivo@univ-reunion.fr

Bruno Baynat
email: bruno.baynat@lip6.fr

Fanilo Harivelo
email: fanilo.harivelo@univ-reunion.fr

Pascal Anelli
email: pascal.anelli@univ-reunion.fr

Table-driven geographical routing using aggregate areas

Keywords: Ad hoc networks, geographical routing, tabledriven routing, aggregation

Geographical routing protocols scale well in large ad-hoc and sensor networks but fall short with some topologies. This causes packet losses and a drop in network performance. In this paper, we propose to introduce topological information aiming to improve the quality of routing decision while keeping the protocol scalable. Every node maintains precise information for nearby nodes and aggregated information for farther nodes by means of aggregate areas. Evaluation demonstrates the scalability of the proposition and its efficiency compared to a pure geographical routing protocol.

I. INTRODUCTION

Ad-hoc and sensor networks involve nodes interacting with each other over multihop connections without central management in an infrastructure-less environment [START_REF] Dressler | Self-Organization in Ad Hoc Networks: Overview and Classification[END_REF]. One of the main issues in such self-organized networks is scalability in terms of routing. With the numerous potential applications of adhoc and sensor networks (Internet of Things, disaster recovery, military application ...) scalability is an important property. It allows the network to grow in scale without quick performance loss. For instance, a scalable routing protocol still ensures high packet delivery rate even with increased number of nodes and also without excessive cost.

Geographical routing approaches are known to scale well in large ad-hoc and sensor networks. They require nodes to maintain existing links to their one-hop neighboring nodes. But additionally, geographical information are needed to perform forwarding. The nodes have to be aware of their location by means of GPS or alternative techniques [START_REF] Zekavat | Wireless positioning systems: Operation, application, and comparison[END_REF]. Forwarding can be summed up as selecting one of the one-hop neighbors which is geographically closer to the destination. Routing decision are purely local. Because the list of existing links to be maintained by each node is small and network size-independent, these approaches are scalable. However, due to lack and/or limited use of topological information, geographical routing performs poorly with some network topologies. Nodes ignore or are not aware of the availability of some links in the network. Choosen paths are not necessarily optimal. In extreme cases, forwarding fails even if routes exist to the destination. One noticeable consequence of such failure is the appearance of dead-ends called local maxima. A local maximum is a node which has not any neighbor geographically closer to the destination, hence forwarding can make no more progress. Several techniques allows to circumvent these dead-ends as described in [START_REF] Chen | A survey of void handling techniques for geographic routing in wireless networks[END_REF], among them, flooding, geometric-based or planarization. Nonetheless, they have limitations such as the need of bidirectional links or nodes with equal radio range [START_REF] Cadger | A survey of geographical routing in wireless ad-hoc networks[END_REF].

We investigate the possibility to enrich a position-based routing with topological information. With an exclusive topology-based approach, each node has a more or less complete knowledge of every other nodes in the network [START_REF] Mauve | A survey on position-based routing in mobile ad hoc networks[END_REF]. More precisely, a list of existing links in the network is maintained at any time (proactive protocols) or determined when needed (reactive protocols). Forwarding decisions are taken based on these link information. One of the advantages of these approaches is the precision as routing is performed only on known links. But they face scaling problems. As a matter of fact, in large networks, the list of existing links can be huge and too costly to maintain.

On one hand, we aim to improve the quality of routing decisions. On the other hand, we would like to keep the scalability of the overall protocol. We provide routing with information about links, and, thus, existing nodes. But not all information is necessary for a given node. Far destinations don't require high precision routes to be reached as routing decision is taken on every hop. Let's imagine a traveler who wants to join his booked hotel room in a foreign country. First, he has to travel to the destination country. It doesn't require very detailed information so far. On arrival, he must find the town where his hotel is located by consulting, for instance, a map of the country. As he gets closer, he now needs more detailed information. And so on, arrived in the town, he will need very sharp details in order to find the quarter and the street where his hotel is located. By analogy, routing needs only more and more details as the destination comes closer. Based on this principle, we propose a scalable approach as it gives more precision to localized routing information while summarizing the global ones. For a node, the list of nearby nodes is maintained as reachable hosts in the routing table. Farther nodes are gradually aggregated into larger and larger areas as they are located away from the considered node. We aggregate nodes by their geographical coordinates.

Our proposition relies on the following components: • At some point, we need to aggregate some nodes into area. Aggregation is done on location-basis. We, then, partition the whole deployment space into geographical areas, that we call primary aggregate partitions or areas, shared and known by all nodes in the network. Primary partitions can be aggregated, in their turn, to form larger partitions to make it possible to refer to faraway nodes. They are uniquely identified and known by all nodes in the network, and so on; • Known nodes and areas are maintained in a routing table by every node; • Nodes exchange periodically routing information to keep their routing table up-to-date; • An unified forwarding strategy consists in sending a packet to the gateway of the entry (either a node or an area) closest to the destination. Mixing topology-based routing and location-based routing was already proposed in the literature [START_REF] Blažević | Self organized terminode routing[END_REF], [START_REF] Couto | Location proxies and intermediate node forwarding for practical geographic forwarding[END_REF]. We make differentiation by using only one type of routing key (the destination location), a tight integration of topology-based and locationbased components (forwarding decision relies solely on the routing table) and more flexible and user-defined settings to tune the tradeoff between precision and cost related to the maintenance of routing information thanks to aggregation. We demonstrate through simulations that our proposition makes larger use of existing paths while keeping the routing table size low.

II. TABLE-DRIVEN GEOGRAPHICAL ROUTING

A. Local area

Each node has its own local area. It is a limited geographical area surrounding a node. Its size can be defined in terms of hop counts, i.e., it is the area containing nodes at a given maximum number of hops from the considered node. The size can also be limited by geographical distances. In this case, the local area is defined by d(N, P ) ≤ l where d is a distance function, N is the considered node, P is a point in the network zone and l is the maximum allowed distance from N to be part of the local area. A node retains link information to all the nodes inside its local area. On the Fig. 1, the local area of the node N is colored in grey.

B. Aggregate areas and network zone partitioning

Link information to nodes inside the network are aggregated by means of geographical areas called aggregate areas. We use the following principles to design aggregate areas:

• We define primary aggregate areas that can be uniquely identified; • Aggregate areas can be, in turn, aggregated into larger aggregate areas that are also uniquely identified; • The aggregation level of each area is known. We define multilevel aggregate areas: order-n areas are aggregated into order-n + 1 where n is the aggregation level. The smallest aggregate areas (primary or order-1) result from a global partition of the network zone.

On Fig. 2a, the network zone has been partitioned into order-1 aggregate areas. Then on Fig. 2b, order-1 areas have been, in turn, aggregated into order-2 ones. The result is a global partition of the network into order-2 aggregate areas.

C. Table structure

From a node point of view, the routing table reflects its knowledge of, on one hand, the link information to nodes inside its local area and, on another hand, aggregated link information to nodes inside aggregate areas. Thus there are two parts inside the routing table:

• Local area part: containing route to nodes inside the local area. Each line in this part indicates how to join a specific node inside this area; • Aggregate areas part: containing information about aggregate areas known to contain at least one node. Each line of this part gives information on how to join a specific order-n aggregate area. The Table I sums up the structure of a routing table.

• Area (A) indicates the nature of the route, either leading to a node (local) or an order-n aggregate area (order-n aggr.). N is the highest order; • Destination (D) tells the ID of the destination of the route which can be a node or an aggregate area; Area ID Node ID ... to the destination according to a specific metric (e.g. hop count). It is important to note that these information are the minimal necessary to be retained in case of a static network. Additional information are required to handle mobility (which is out of the scope of this paper) such as sequence numbers [START_REF] Perkins | Highly dynamic destination-sequenced distance-vector routing (DSDV) for mobile computers[END_REF].

D. Updating tables

Because of the table-based nature of our protocol, one of the key points is the updating process of the table content. We propose an updating mechanism which is similar to what is used in DSDV [START_REF] Perkins | Highly dynamic destination-sequenced distance-vector routing (DSDV) for mobile computers[END_REF]. Every node periodically broadcasts its table content to all of its one-hop neighbors. Restraining the broadcast to the one-hop neighbors saves bandwidth yet allow the routing information to spread across the network.

Aggregation decisions. As stated earlier, a node's routing table reflects its knowledge of, one one hand, the link information to nodes inside its local area and, on another hand, aggregated link information to nodes inside aggregate areas. In other terms, information outside the local area are aggregated. In order to optimize the scalability of our solution, we pose the following principle: with great distance comes high level of aggregation. More precisely, from a node perspective, the further is the node corresponding to a link information, the higher is the order of the aggregate area containing it. The order of the aggregate area is function of the distance between the node receiving the link information and the node concerned by the latter:

order =          1 if 0 < d(N, P ) ≤ x 1 2 if x 1 < d(N, P ) ≤ x 2 ... n if x n-1 < d(N, P ) ≤ x n
where N is the center node (owner of the table), P is a point in the network zone, d is a distance function and x 1 < x 2 < ... < x n are reference distances.

The Fig. 3 shows the aggregate areas (with their respective order) considered by the node N.

Initialization. At start up, a node adds its own information to its table: the area column is filled with local, the destination with its own ID, the gateway also with its own ID and the cost with zero.

Receiving updates. Let N a node receiving a routing message for a destination D from one of its neighbors noted GW. The routing table of N will be referred as T(N). N will follow the next algorithm in order to update its routing table given this new information. In the following, C is an update of the cost C after routing update processing.

• If D is a node:

- • If D is an aggregate area:

If
-Find the smallest aggregated area known by N containing D, let A the resulting aggregate area whose order is n: * If there is not yet route to A inside T(N), add the route {Order-n aggr., A, GW, C }; * If there is already a route to A inside T(N), it will be replaced if the new route offers a lower cost. Example. On the Fig. 4 we have a representation of the local area and aggregate areas of a node N. We assume that N's table has just been initialized so contains only one entry (its own). Now, N receives the following routing messages {Destination,C} from a node that will be referred as GW:

• {A, C}: A is situated inside the local area, so N adds the line {local, A, GW, C } to its table; • {B, C}: B is outside the local area but inside area_-1. Thus, the routing table is updated with {Order-1 aggr., area_1, GW, C }; • {area_gw, C}: area_gw is an order-1 aggregate area. At this distance, N doesn't retain any order-1 area but order-2. So area_gw is aggregated into order-2 area that is area_2. Then the added entry is {Order-2 aggr., area_2, GW, C }.

E. Forwarding

Forwarding decisions rely on information contained in the routing table.

When a node receives a packet, first it checks if it is the destination. If not, it has to choose a node to which the packet will be forwarded. This process is repeated by every node until the packet reaches the final destination.

We recall the assumption that the sender node is aware of the geographical location of the destination either inferred from its ID or from a location service [START_REF] Das | Performance comparison of scalable location services for geographic ad hoc routing[END_REF]. Also, each aggregate area has a defined geographical location corresponding to that of its reference point. The reference point of an aggregate area can be defined as its geometric centroid.

We propose a greedy-like forwarding process. In pure greedy routing [START_REF] Cadger | A survey of geographical routing in wireless ad-hoc networks[END_REF], the next hop for a forwarded packet is one of the one-hop neighbors nodes that is closer to the destination. In our protocol, the next hop is chosen as the gateway to the node or the area (regarding its reference point) that is closest to the destination. We use euclidean distances for proximity comparisons.

Algorithm. Let N a node receiving a packet whose final destination is D.

• Select all the gateways offering a node or an area closer to D:

-If no gateway is found, forwarding fails; -Else, among the selected gateways, choose the one offering the node or the area closest to D. Pure greedy routing fails when no neighbor is at a smaller distance to the destination (local maximum). In our solution, forwarding fails when all the areas reachable through the gateways are at further distances from the destination than the source. But the failure rate is expected to be lower. There are two main reasons for that. First, the reachable areas are not restrained to the immediate neighborhood of a node. They can be located far from the farthest one-hop neighbors. Thus chances are high to find an area closer to the destination. Secondly, the one-hop broadcast and non geographical nature of the the updating mechanism allows the routing information to circumvent local maxima. Thus some forwarding candidates areas are made reachable even through gateways that are at greater distance from the destination than the source.

Example. On Fig. 5, the node N has to forward a packet whose final destination is D. In its routing table it has the following entries among others:

• {Order-1 aggr., area_1, A, C}: A is the gateway to reach the area area_1; • {Order-2, area_2, B, C}: B is the gateway to reach the area area_2. Because area_2 is closer to D than area_1, the packet will be forwarded to B even if B is farther from the destination than N.

F. Discussion

In this section we evaluate qualitatively our solution by discussing its performance according to some metrics.

1) Table size:

The routing table size (the number of entries inside tables) has a direct impact on the amount of memory consumed by every node in order to run the protocol properly. In classical table-based routing protocols, for a connected network, the table size is linear with the number of nodes in the network as ultimately each node will trigger an entry. Our proposition reduces the size of tables thanks to the concept of areas which aggregates node inside them. Furthermore, smaller tables leads to smaller route updates meaning reduced overhead.

2) Convergence time: Given a static network where all nodes start at the same time, the convergence time is the necessary amount of time for each node's routing table to stabilize. Compared to classical table-driven routing [START_REF] Mohseni | Comparative review study of reactive and proactive routing protocols in MANETs[END_REF], our solution needs shorter time to converge. In fact, a routing update leads to table modification in a limited range of nodes. As seen previously, when a node receives a route update, it looks for an aggregate area for the destination. The further the destination is, the higher the probability that the aggregate area is a big one. In this case, chances are high that a route with the same aggregate area has already been received, so no update occurs.

3) Packet delivery ratio (PDR): In the case of static network, our proposition is expected to be a tradeoff between a full table-driven routing protocol and a geographical one. More precisely we are improving the table-driven aspect by reducing the table size and bandwidth consumption (as seen before) with the cost of lower PDR because of routing information aggregation. But the use of routing tables can circumvent some of the local maxima making our solution more efficient than pure geographical routing protocol.

4) Mobility:

The study of mobility needs extensive simulation and is out of the scope of this paper. However we wanted to highlight a potential advantage of our solution. Thanks to the aggregate areas and routing tables having areas as entries, we tolerate some local mobility of distant nodes. In fact, for a given node, mobility of other nodes that are confined inside an aggregated area don't trigger any update. In other words, they mostly generate updates at a local level.

III. SIMULATION

A. Objectives

The main goals are to show that our proposition is more efficient compared to pure geographical routing while having good scaling capabilities in terms of routing table size and convergence time. We measured table size, convergence time and packet delivery ratio (PDR) in regard to nodes number and nodes density.

B. Local area and aggregate areas implementation

In order to partition the network zone, we chose the hierarchical area partitioning proposed in [START_REF] Ranaivo Rakotondravelona | Hierarchical area-based address autoconfiguration protocol for self-organized networks[END_REF]. Their algorithm consists in a recursive subdivision of the whole network zone into square areas: 4 order-n areas form one order-n + 1 area. The size of the primary areas (order-1) is chosen such as two nodes inside the same order-1 area can reach each other.

The local area of a node is implemented as the geographical region formed by the order-1 aggregate area containing the node and all its adjacent order-1 areas (Fig. 6a).

Order-n aggregate areas are implemented as order-n areas adjacent to the order-n area containing the node. Fig. 6b shows examples of aggregate areas with their respective orders.

C. Simulation parameters

We conducted our simulation using ns-3 (ns-3.24 release). Each node implements the 802.11b model with 2 Mbps rate and 250 m range. The nodes are randomly deployed in the network zone and remain static. They are powered on at the same time. For routing table size and convergence time measurements, we simulated 20 to 300 nodes with constant density of 1 node per 9000 m 2 . For PDR scenario, we simulated 200 nodes with density varying from 0.25 node per 9000 m 2 to 1.25 nodes per 9000 m 2 . The nodes are deployed in variable size square regions to reflect the density. 10 nodes are selected randomly to generate 2 kbps CBR flows with 64 bytes packets sent to other 10 random nodes. The flows are started after the convergence of the routing tables and last for 10 seconds. There is no loss due to interference between flows as they are not concurrent. We chose an interval of 15 s between routing periodic updates. A simulation ends 10 s after the last sent packet.

D. Results

1) Routing table size:

The Fig. 7a shows the average number of entries in the routing tables with different network sizes. As we can see, the average number of routing table entries for increasing nodes number has a logarithmic shape. This shows the scalability of our proposition as the average routing table size doesn't explode with large number of nodes in the network.

2) Convergence time: The convergence time is the necessary amount of time for all the nodes tables to stabilize after start-up for a static network. Before convergence, the tables have incomplete routing information resulting in poor forwarding decisions. Thus the convergence time has an impact on how fast a network can be deployed. It is interesting to see how this metric evolves with the size of the network as short deployment time even with high number of nodes can be considered as a criteria of scalability. The Fig. 7b shows that the convergence time grows slowly with the number of deployed nodes.

3) Packet delivery ratio (PDR): The PDR is the ratio between the number of packets successfully received by the destination and the total number of packets sent by the source. The Fig. 7c shows the average PDR with different node densities. As we can see, with a large range of node densities, our proposition brings substantial improvement to the PDR compared to pure greedy routing. It demonstrates its capacity for circumventing the local maxima, reducing forwarding failures. proaches, a comprehensive survey can be found in [START_REF] Cadger | A survey of geographical routing in wireless ad-hoc networks[END_REF]. To avoid local maximum handling, some geographical routing protocols construct landmarks [START_REF] Na | GLR: A novel geographic routing scheme for large wireless ad hoc networks[END_REF] or spanning tree-like structure [START_REF] Leong | Geographic routing without planarization[END_REF]. The idea of using areas for routing purpose has been already proposed. ZRP [START_REF] Patel | Performance analysis of zone routing protocols in mobile ad hoc networks[END_REF] defines two areas for applying respectively reactive and proactive topology-based routing protocol. Another protocol using areas not in a geographical way is DART [START_REF] Eriksson | DART: Dynamic address routing for scalable ad hoc and mesh networks[END_REF] which uses its own adressing scheme to define logical hierarchical areas. In a more geographical fashion, [START_REF] Hsiao | Geographical region summary service for geographical routing[END_REF] proposes region summaries of nodes with routing on them. Regarding variable routes details, FSR [START_REF] Iwata | Scalable routing strategies for ad hoc wireless networks[END_REF] takes into account the proximity of nodes to each other to adapt routing updates precision.

V. CONCLUSION

We propose in this paper, a geographical routing protocol provided with topological information. Every node maintains a routing table containing hosts entries for closer nodes and areas entries aggregating further nodes. To make it possible for the node to aggregate and exchange nodes information, the deployment area is splitted into partitions. Routing decisions are taken solely on routing table contents: a packet is sent to the gateway of entry nearest to the destination.

Evaluation compares the proposition performances with those of a pure geographical routing. It shows higher Packet delivery ratio validating that the proposed protocol makes a better use of link information by reducing forwarding failures. The size of the routing table and the convergence time increase slowly in logarithmic shape with the number of nodes, demonstrating a good scale factor stability.

The way our proposal introduces topological information in location-based routing, can be seen as widening the possibility for a node to choose the next hop. While on classical geographical routing, the next hop is a node, in our proposition, the next hop can be considered as an area: geographical routing is then extended to areas. Future work will study quantitatively the impact of mobility on protocol performance.

Fig. 1 .

 1 Fig. 1. Local area

2 Fig. 2 .

 22 Fig. 2. Example of aggregate areas

Fig. 3 .

 3 Fig. 3. Aggregate areas of N

  D is located inside the local area: * If there is not yet route to D inside T(N), add the route {local, D, GW, C }; * If there is already a route to D inside T(N), it will be replaced if the new route offers a lower cost. -If D is located outside the local area: * Identify the order-1 aggregate area containing D; * Restart the algorithm with the latter area as entry D.

Fig. 4 .

 4 Fig. 4. Receiving routing information

Fig. 5 .

 5 Fig. 5. Example of forwarding

Fig. 6 .

 6 Fig. 6. Local area and aggregate areas

Fig. 7 .

 7 Fig. 7. Simulation results

ACKNOWLEDGMENTS

This research was supported by Région Réunion and European Union (ERDF).


Jean-Sylvain Liénard
email: jean-sylvain.lienard@limsi.fr

Claude Barras
email: claude.barras@limsi.fr

Fine-grain voice strength estimation from vowel spectral cues

Keywords: vocal effort, vocal intensity, voice quality, discriminant analysis

This study investigates the possibility to recover the voice strength, i.e. the sound level produced by the speaker, from the signal recorded. The dataset consists of a set of isolated vowels (720 tokens) recorded in a situation where two interlocutors interacted orally at a distance comprised between 0.40 and 6 meters, in a furnished room. For each token, voice strength is measured at the intensity peak, and several sets of acoustic cues are extracted from the signal spectrum, after frequency weighting and intensity normalization. In the first phase, the tokens are grouped into increasing voice strength categories. Discriminant Analysis produces a classifier which takes into account all the signal dimensions implicitly coded in the set of cues. In the second phase, the cues of a new token are given to the classifier, which in turn produces its distances to the groups, providing the basis for estimating the unknown voice strength. The quality of the process is evaluated either in self-consistency mode or by cross-validation, i.e. by comparing the estimate with the value initially measured on the same token. The statistical margin of error is quite low, of the order of 3 dB, depending on the sets of cues used.

Introduction

It is common experience to recognize the voice strength adopted by a speaker, whatever the sound level actually received by the listener. An amplified weak voice remains perceived as a weak voice. This is true not only for the extreme modes such as shout or whisper, but also for the conversational speaking modes used everyday by any of us.

The objective of the present study is to demonstrate that, at least in the limited framework of an isolated vowel dataset, it is possible to recover the sound level emitted by the live speaker, by using acoustic cues extracted from the signal normalized in intensity.

State of the art

Our objective is a matter of several knowledge fields. The following references are not exhaustive. Relevant research work has been done on shouted speech intelligibility [START_REF] Rostolland | Acoustic features of shouted voice[END_REF], on the role of the distance between interlocutors [START_REF] Traumuller | Acoustic effects of variation in vocal effort by men, women and children[END_REF], on the effects of an excessive vocal effort [START_REF] Junqua | The Lombard reflex and its role on human listeners and automatic speech recognizers[END_REF] [START_REF] Garnier | Communiquer en environnement bruyant: de l'adaptation jusqu'au forçage vocal[END_REF], on the functioning and modeling of the vocal source [START_REF] Fant | A four parameter model of glottal flow[END_REF][6] [START_REF] Hanson | Glottal characteristics of female speakers: acoustic correlates[END_REF], on the effect of the vocal effort on the acoustical structures of the speech signal [START_REF] Liénard | Effect of vocal effort on spectral properties of vowels[END_REF] [START_REF] Huber | Formants of children, women, and men: the effects of vocal intensity variation[END_REF], on the quality of the singing voice [START_REF] Henrich | Glottal open quotient in singing: measurement and correlation with laryngeal mechanisms, vocal intensity, and fundamental frequency[END_REF], on speech prosody [START_REF]Voice source parameters and prosodic analysis[END_REF], to mention just a few.

Most authors agree on the following considerations:

• speaking louder yields a frequency shift of the so-called "glottal formant", which is not to be confused with the modes of the vocal tract. This seems to be due to the shrinking of the open phase of the glottis. Usually located in the region of the first harmonics for a weak voice, the glottal formant may shift towards the higher harmonics in loud voice.

• the magnitude difference between the 1st harmonic H1 and the second H2 is strongly related to the open quotient Oq, but it also depends on the location of the 1st formant F1, i.e. on the vowel uttered.

• speaking louder or shouting yields a relative spectral emphasis of the high frequencies named spectral tilt. This feature is often estimated as the difference between H1 and the magnitude A3 of the 3rd formant.

• the fundamental F0 increases with the vocal effort; an order of magnitude is 15% per 6 dB intensity increase.

• the frequency F1 of the first formant, correlated with the speaker mouth aperture (jaw lowering), may increase of a quantity comprised between a few Hz to 200 Hz. The other formants do not appear to be systematically altered. All of those effects are mixed with the other characteristics of the signal such as its phonemic and prosodic structures. Thus they are extremely difficult to extract from the signal without prior knowledge.

Conditions of the study

Sound intensity is the very object of our investigation. Thus it is essential to clearly distinguish the sound level produced by the speaker in live conditions, from the level of the same signal restituted by an electroacoustic device.

The concept of Vocal Effort being somewhat qualitative, we need an objective notion to represent it. We hypothesize that vocal intensity, measured in well defined conditions, is a good indicator for it. Vocal intensity should be expressed in dB SPL as in most of the above references, but in the present study we prefer the term "Voice Strength" because our data were not calibrated with a sound meter and thus cannot be taken as absolute measures, although they are reproducible throughout the database. Another reason is that we will use a spectral weighting, which greatly changes the figures representing sound intensity.

The voice strength shades occur and are recognized in everyday situations; instead of well controlled laboratory speech we will use data recorded in ordinary oral communication situations and exclude for the time being the modes of shouting and whispering which are less frequently used.

As many unknown factors may contribute to the determination of voice strength, our analysis method will not take for granted that voice strength varies linearly with any of the cues extracted from the signal.

Data and method

Data

The databases available in the speech community do not guarantee a rigorous constancy of the recording conditions, which is necessary to investigate the voice strength problem. Thus we used ours, a database named CORENC, formerly described and used in ref [START_REF] Liénard | Effect of vocal effort on spectral properties of vowels[END_REF]. It is composed of 12 French vowels (9 orals, 3 nasals) uttered in isolation by 13 speakers (6 males, 7 females, aged 19-89). Three degrees of vocal effort noted "p", "n" and ""l were elicited by varying the distance between speaker and operator (0.4 m, 1.50 m and 6 m); however the subjects did not always produce the expected degree of vocal effort, so that the actual voice intensity was only loosely related to the distance condition. Consequently the latter was not taken into account in the present study. The recording took place in a quiet furnished room in 3 sessions, the 3rd of which was done six months after the first. Each speaker participated in a number of sessions comprised between 1 and 3. The speaker was seated, his/her mouth 30 cm away from the omnidirectional microphone. The task was to repeat interactively a short sentence and the set of vowels uttered by the operator. No calibration of the absolute sound level was done; however the recorder settings were kept identical during all of the 3 sessions, allowing a consistent measurement throughout the recordings.

Signal analysis

Each of the 720 tokens is represented by two entities. The first one is the peak intensity a of the token. The 2nd entity is a set of acoustic cues measured on the spectrum after its normalization to an arbitrary value that eliminates the voice strength direct information. Then a Discriminant Analysis classifier is trained on the data, in order to give an estimate â of the voice strength from the cues alone. The performance criterion is the standard deviation of the difference <â-a>, computed on the whole database. This margin of error is expressed in dB like the voice strength. The weaker the error, the more relevant the set of cues. The study has been done with the Praat software [START_REF] Boersma | Praat: doing phonetics by computer[END_REF].

Signal intensity and frame selection

Signal intensity is computed on a 50 ms sliding Gaussian window. The frame is selected at the time where intensity reaches its peak value. This maximum intensity ax expressed in dB is considered as the voice strength of the token.

A second measure ap is performed in the same conditions, except that the signal is high pass filtered according to a function close to the A-weighting curve. The latter is used in sound level metering to take into account the reduced sensitivity of the human ear in the low frequencies, especially at low to medium level. Its origin goes back to the early work on loudness [START_REF] Fletcher | Loudness, its definition, measurement and calculation[END_REF]. The attenuation goes approximately from 0 dB at 1000 Hz to 3 dB at 500 Hz and 23 dB at 63 Hz. The intensity of the A-weighted signals is expressed in dBA.

In the following, axc and apc represent the estimates of the measured values ax and ap.

Frame analysis

The signals are analyzed by a filter bank of 18 1-Bark wide filters. The outputs are integrated on a 50 ms Gaussian window. Then the frame spectrum is normalized with respect to its maximum.

According to the indices under study one can adopt the minimal resolution (some global cues may be computed from 18 channels separated by 1 Bark) or a better resolution by interpolation (when looking for spectral details, or for a visual representation). The fundamental frequency F0 is a basic cue, generally retained among the other acoustic cues. The Bark spectral modules are treated in Praat as Ltas objects (Longterm average spectrum).

Sets of acoustic cues

Many sets of acoustic cues can be extracted from the Ltas. In the present study the following sets were selected:

• F0 and Bark channels (b1 to b18): amplitudes in dB after normalization of the maximum to an arbitrary 50 dB.

• F0 and slopes: overall slopes (p5k is the amplitude difference between the high and low parts of the spectrum, trn and trn0 both represent the slope of the Ltas trend line); slopes in the low-part of the Ltas (p1k is the amplitude difference between the two halves of the 1-9 Bark band, d0 is the difference at the origin between the trend line and the Ltas); local slopes bij between adjacent channels i and j.

• F0 and voice perturbations: phonation irregularities such as jitter (jit), shimmer (shi), mean F0 (F0m), maximum autocorrelation (cor), harmonicity (hnr), all estimated on the whole token duration.

• F0 and Centres of Gravity taken on the whole Ltas (fg0, ag0), as well as on halves and quarters (fgi, agi, order i from 1 to 6). These cues must not be confused with the physical centres of gravity defined in linear frequency and magnitude scales. For those measures an intensity threshold is set at -40 dB from the maximum.

• F0 and 10 harmonics: amplitudes in dB of the first 10 harmonics, measured from the maximum on the interpolated Ltas at frequencies multiples of F0.

Discriminant Analysis

In order to establish a potentially non-linear matching between cues and voice strength, the tokens are grouped into a set of voice strength categories, according to a 3 dB step. Each category is given a label representing the mean voice strength of its component tokens. Then labels and cues are treated by Discriminant Analysis which, in the training phase, determines a set of orthogonal axes and their eigenvectors by maximizing the ratio between intergroup and intragroup variances. The eigenvectors are arranged into a set of functions called classifier (or discriminant). In the decision phase, the cues of a new token are presented to the classifier, which computes its distances to the group centers and proposes a set of candidates in terms of probabilities. The voice strength estimate is recovered by interpolating between the closest candidates.

Self-consistency and cross-validation modes

In self-consistency mode all the tokens are used to train the classifier; the same data are used as test. This mode is equivalent to a classification. It is legitimate if the goal is to investigate the structures of a closed dataset, but it is biased if the goal is to use the classifier with new, unseen data. We check this bias by performing a cross-validation. As the corpus consists of 20 series of 36 tokens, a set of 19 series serves for training and the remaining series is used as test. The process is applied again 19 times within a circular permutation, and the error margin is computed as before on the 720 estimates.

Results

Sets of cues and A-weighting

F0 and 18 Bark channels, A-weighted signals

This set gives good results: error margins amounts to 2.94 dBA in classification, 3.25 dBA in cross-validation. This case will be used to illustrate the classification processing.

The voice strength dimension is quantified into 12 groups covering the 50 to 85 dBA interval in 3 dBA steps. The groups are labeled "05" to "16". Fig 1 (top left) displays the dispersion areas (at 1 sigma) of the groups in the plane of the 2 first classification functions. This shows the continuity and overlap of the groups, as well as the non-linearity that mainly affects the low end of the voice strength scale. The results are almost equivalent in terms of error margin: 2.86 dB in self-consistency, 3.15 in cross-validation. However the non-linearity of the dispersion areas has augmented. More importantly, the (ax,axc) distribution is much farther from the diagonal, the correlation coefficient is poorer, and the number of groups built by the Discriminant Analysis is lower: we now have 8 groups instead of 12. Besides, the voice strength dynamics is poorer (some 25 dB, instead of 35), while the error margin remains comparable, about 3dB. All this indicates that the A-weighting favors to a large extent the coding of the voice strength, especially for the normal to low voices. In the rest of the study only the Aweighted signals will be taken into consideration.

Same analysis with unweighted signals

F0 and 4 slopes: p1k, p5k, trn, d0

This set provides a poorer margin of error: 3.57 dBA, 3.75 in cross-validation. However, as only 5 cues are used, they can be considered more relevant and more general than the above 18 cues, except b1 and b2. As in the previous experiments, the weight of F0 in the main eigenvectors is low.

F0 and 6 voice perturbation cues

For this set of cues the performance is poor (4.59 dBA, 4.75 in cross-validation) and the (ap,apc) distribution diverges from the diagonal. Actually, this result is much better than expected, considering that it does not contain any spectral information. The most significant cues seem to be cor, F0m and hnr.

. F0 and 6 cues from 3 centres of gravity

This set comprises 7 cues: F0 plus the amplitudes and frequencies of 3 centres of gravity of orders 0, 1 and 2. The results are poor (3.88 dBA, 4.05 in cross-validation) and the examination of the (ap,apc) distribution shows non linearity and a lack of efficiency in its lower part. Again, F0 seems to play a secondary part. f0 b1 b2 b3 b4 b5 b6 b7 b8 b9 b10b11b12b13b14b15b16b17b18 07 0.1 0.9 0.9 1.0 1.3-1.4-1.6-0.9-1.2-1.1-0.9-1.0-0.7-0.8-0.

F0 and 10 harmonics

This set produces an acceptable performance (3.28 dBA, 3.52 in cross-validation), mainly due to F0 and h1 which play in opposite directions.

Compound sets

In this section we compose new sets by selecting some of the most significant cues found in the previous experiments.

34-cue set

This large set includes the best cues found in each of the previous sets. The results are better than with any of the initial complete sets: 2.59 dBA, 3.03 in cross-validation. The (ap,apc) distribution is closer to the diagonal and more symmetrical than in any of the previous experiments. The crescent shape observed in the main plane has no negative consequence on the performance, because its projection on the main axis, which explains 85% of the intergroup variance, is perfectly regular. It suggests that different cues contribute to the high and low parts of the voice strength interval.

Looking for the most important cues

Examining the main eigenvectors allows a ranking of the most important cues: p5k, d0, ag5, ag6, trn0, ag1. However, gathering them into a new, reduced compound set yields a margin of error that is not bad, but not much better than some of the primary sets: 3.37 dBA (3.49 in cross-validation). Thus in the present state of the study one can only propose the following qualitative observations:

• the global slope of the A-weighted Bark Ltas, taken from the medium to the high frequency bands, seems to be the main underlying factor of voice strength estimation. It is reflected in the cues that appear frequently at the top of our cues lists, such as p5k, trn0, ag5 and ag6. It is obviously related to the spectral tilt evidenced in many previous studies.

• The slope in the lower part of the Ltas also plays an important role, which appears in cues such as d0, ag1, or b1 or h1 in previous experiments. It may be related to the H1-H2 or H1-A1 parameters found in other studies.

• F0, taken as an absolute cue, seems to play a secondary role. It could be rated differently if it were taken as a difference, for a given talker, between the instantaneous and the mean F0 values.

Discussion

The main point that emerges from the present study is that in conditions close to ordinary oral communication, the speech signal, even reduced to a mere spectrum taken in the central part of a vowel, contains enough information to characterize the voice strength within a very low margin of error.

The closeness of the figures obtained in self-consistency and in cross-validation indicates that the proposed processing does not depend much on the variations due to new data. However some variability due to the talker and phonetic content may remain in the results, which leaves the door open to further improvements.

The use of a frequency weighting close to the psychoacoustically grounded A sound metering norm, is amply justified by the results, because it yields more degrees in the estimation of the voice strength.

In vocal quality studies, vocal effort in the conversational range is usually qualified with 3 timbre qualities only: modal voice, loud and low voice. This particular aspect of timbre could be replaced by voice strength which, considering the 3 dBA margin of error and the 35 dBA interval, may provide a number of degrees closer to 10 than to 3.

The link between vocal effort and F0 has been widely demonstrated experimentally in the literature. Talking louder usually yields an irrepressible raise of F0, except for especially trained talkers. The reciprocal is not true: any talker can raise F0, up to about one octave, without augmenting his/her voice strength. This capability is used to transmit intonational information. Besides, at equal voice strength, F0 is strongly related to the talker's position in the {male, female, child} vocal categorization. As we choose from the start of the study not to use any prior information on the talker and context, the secondary role found for F0 simply confirms the nonreciprocity of the relation between F0 and voice strength.

Conclusions

The approach proposed, based on discriminant analysis and on a database with controlled recording conditions, has shown that voice strength could be deduced from an uncalibrated speech signal with surprising precision. The method may prove useful in other aspects of voice and speech research, for instance for investigating the relations between phonemic systems and speaker characteristics. Some developments may be envisioned in several fields of automatic processing. In synthesis, mastering the voice strength may help to produce more natural, situation-sensitive voices. In speech recognition, knowledge of the voice strength may contribute to reduce the variability of the linguistic speech structures. In speech diarization or speech scene analysis, differences in voice strength may contribute to separating the talkers voices.

Acknowledgements

The authors would like to thank the scientific committee of LIMSI-CNRS, as well as their colleagues Drs Christophe d'Alessandro, Philippe Boula de Mareüil, Albert Rilliard from LIMSI, and Pierre Divenyi from CCRMA, Stanford Univ. f0Hz p5k trnd0 fg1 ag1 b1b3 b16 b18 ag3 ag5 ag6 trn0 trn1 trn2 jitshi brk f0m cor hnr h1h2h3h4h5h6h7h8h9 h10 070.1 -0.5 -0.5 1.0 -1.3 -0.9 0.9 1.0 -0.1 0.2 1.0 -1.2 -0.6 -0.8 -1.4 0.7 -0.6 -0.0 -0.2 0.8 1.2 -1.0 0.8 1.1 -1.4 -0.8 -1.2 -1.3 -1.2 -1.2 -1.2 -0.9 080.8 0.6 0.8 1.2 -0.9 -0.2 0.8 0.3 0.7 1.6 0.8 -0.4 0.9 0. eigensum=7.0 eigen1=0.852 eigen2=0.086

1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 2 2 2 2 2 2 2 2 2 2

Figure 1 :

 1 Figure 1: Dispersion ellipses in the main factorial plane (top left); profiles of the first two eigenvectors (bottom left); fractional eigenvalues (bottom right); voice strength estimate vs. measure distribution (top right).

Fig 1 (

 1 Fig 1 (bottom) displays the components of the two main eigenvectors in terms of the 19 cues, and the fractional eigenvalues. The main eigenvector explains 91% of the variance and the second one only 4%. The first vector is dominated by the amplitude b1 of the 1st Bark channel. The second is dominated by b2. The contribution of the other cues is less prominent, at least in the main plane.

Fig 1 (

 1 Fig 1 (top right) shows that, despite the quantification and the non-linearities, the (ap,apc) distribution is well concentrated along the diagonal. A slight asymmetry appears at the low end, meaning that this set of cues overestimates the voice strength of the weakest tokens.

Fig 2

 2 Fig 2 shows the same experiment on the unweighted signals. The results are almost equivalent in terms of error margin: 2.86 dB in self-consistency, 3.15 in cross-validation. However the non-linearity of the dispersion areas has augmented. More

Figure 2 :

 2 Figure 2: Dispersion areas and (ax,axc) distribution in the absence of the A-weighting of the signals.

Figure 3 :

 3 Figure 3: Dispersion areas and (ap,apc) distribution for the F0+6 voice cues 3.1.5. F0 and 6 cues from 3 centres of gravity

  19 difstd=2.94 corr=0.905 f0 b1 b2 b3 b4 b5 b6 b7 b8 b9 b10b11b12b13b14b15b16b17b18 10 0.1 0.3 0.8 0.0 0.7-1.3-1.6-1.1-1.3-1.3-1.2-1.3-1.1-1.3-0.9-1.0-0.6-0.7-0

Figure 4 :

 4 Figure 4: Dispersion areas and (ap,apc) distribution for the compound set of 34 cues.

  3 -1.1 1.3 0.4 1.3 -0.2 1.7 0.0 -0.2 1.3 0.3 -0.1 -1.6 -1.2 -0.8 -0.8 -0.4 0.1

difmoy=-0.05 difstd=2.59 corr=0.927


Augusto Ferrante
email: augusto@dei.unipd.it

Alexander Lanzon

Bernard Brogliato

A direct proof of the equivalence of side conditions for strictly positive real matrix transfer functions

direct way that two different side conditions, which have been used in the literature to characterize strictly positive real matrix transfer functions in the frequency domain, are equivalent.

I. INTRODUCTION

The frequency domain conditions characterizing the fact that a matrix transfer function F is strictly positive real involve a positivity constraint at infinite frequency. This constraintusually referred to as side condition -has been a source of confusion and controversy in the literature for more than a decade. As pointed out in [START_REF] Corless | On the characterization of strict positive realness for general matrix transfer functions[END_REF], the side conditions used in [START_REF] Khalil | Nonlinear Systems[END_REF], [START_REF] Tao | Strictly positive real matrices and the Lefschetz-Kalman-Yakubovich lemma[END_REF], [START_REF] Tao | Necessary and sufficient conditions for strictly positive real matrices[END_REF] were incorrect as they had some inconsistencies. To fix the problem, [START_REF] Corless | On the characterization of strict positive realness for general matrix transfer functions[END_REF] proposed a new condition at infinite frequency, i.e.

lim ω→∞ ω 2ρ det F (jω) + F (-jω) > 0 (1)
where ρ is the dimension of ker F (∞) + F (∞) . On the other hand, a different, but equally valid, condition at infinite frequency was proposed the second edition of the book by Khalil published in 1996 (see [START_REF] Khalil | Nonlinear Systems[END_REF]Lemma 10.1]); such a condition, that reads as follows, was recently used in [START_REF] Ferrante | Foundations of not necessarily rational negative imaginary systems theory: relations between classes of negative imaginary and positive real systems[END_REF] to establish a counterpart result for negative imaginary systems (see [START_REF] Ferrante | Foundations of not necessarily rational negative imaginary systems theory: relations between classes of negative imaginary and positive real systems[END_REF]Remark 1]

): ∃ δ > 0, σ 0 > 0 such that σ ω 2 F (jω) + F (-jω) ≥ σ 0 ∀|ω| ≥ δ. ( 2
)
This note is devoted to the analysis of the two sideconditions (1) and (2). We will prove that while they are in general not equivalent at infinite frequency, they are indeed equivalent under the other conditions guaranteeing that F is strictly positive real. Hence both conditions at infinite frequency are equally valid. While this could be deduced from [START_REF] Corless | On the characterization of strict positive realness for general matrix transfer functions[END_REF] and [START_REF] Khalil | Nonlinear Systems[END_REF], our results provide a direct proof of such equivalency. Notation: Let the set of real (resp. complex) numbers be denoted by R (resp. C) and the corresponding sets of matrices bernard.brogliato@inria.fr of dimension m×n be denoted by R m×n (resp. C m×n ). Given M ∈ C m×m , M denotes the complex conjugate transpose of matrix M (i.e., if M = A + jB for real matrices A and B, then M = A -jB ). A matrix M is said to be Hermitian if M = M and M > 0 denotes that the matrix M is Hermitian and positive definite. The smallest singular value of M is denoted by σ(M ). We recall that the singular values of a positive semi-definite Hermitian matrix are its nonzero eigenvalues [2, p.649].

II. MAIN RESULT

The following definition, adapted from [4, Definitions 1 and 2], is the standard definition for strictly positive real systems. It essentially states that a transfer function matrix F (s) is strictly positive real if for some > 0, the transfer function matix F (s -) is positive real and F (s) + F (-s) has full normal rank. See also [START_REF] Ferrante | Foundations of not necessarily rational negative imaginary systems theory: relations between classes of negative imaginary and positive real systems[END_REF]Lemma 2] for an equivalent re-characterization.

Definition 2.1: Let F : C -→ C m×m be a real transfer function. Then F (s) is said to be Strictly Positive Real (SPR) if there exists a real scalar > 0 such that F (s) is analytic in {s ∈ C : Re{s} > -}, F (s) + F (s) ≥ 0 for all s ∈ {s ∈ C : Re{s} > -} and F (s) + F (-s) has full normal rank. SPR matrix transfer functions can be characterized in the frequency domain by three conditions: the first two are conditions 1 and 2 in the next proposition, the third is the side condition and it has been stated in two different manners: side condition (3a) in Proposition 2.1 can be found in [START_REF] Khalil | Nonlinear Systems[END_REF], [START_REF] Ferrante | Foundations of not necessarily rational negative imaginary systems theory: relations between classes of negative imaginary and positive real systems[END_REF], while side condition (3b) can be found in [START_REF] Corless | On the characterization of strict positive realness for general matrix transfer functions[END_REF]. These side conditions can be interpreted as apparently different conditions on how F (jω) + F (-jω) approaches zero for sufficiently large |ω| in directions where it looses rank.

Proposition 2.1: Let F : C -→ C m×m be a real, rational, proper transfer function such that the following two conditions hold:

1) F (s) has no poles in {s ∈ C :

Re{s} ≥ 0}; 2) F (jω) + F (-jω) > 0 for all ω ∈ R.
Then, the following two side conditions are equivalent:

3a) ∃δ > 0, σ 0 > 0 : σ ω 2 F (jω) + F (-jω) ≥ σ 0 ∀|ω| ≥ δ (3) 3b) lim ω→∞ ω 2ρ det[F (jω) + F (-jω) ] = 0 (4)
where

ρ = dim ker(F (∞) + F (∞) ). Proof. Since F (s) is proper, let F (s) = C(sI -A) -1 B + D
for some state-space realization (A, B, C, D). Near s = ∞, the expansion

F (s) = F 0 + F 1 /s + F 2 /s 2 + . . . (5)
holds (with

F 0 = F (∞) = D), so that F (jω) + F (-jω) = Q + H ω + K(ω) ω 2 (6)
with

Q = D + D , H = (F 1 -F 1 )/j and K(ω) = ∞ i=1 Fi+1 j i+1 ω i-1 + F T i+1 (-j) i+1 ω i-1 .
Then it is easy to see that

lim ω→∞ K(ω) = K 0 ∈ R m×m , with K 0 := -F 2 -F T 2 .
If Q is positive definite, in view of the Hermitian symmetry of both H and K(ω), both conditions 3a) and 3b) are clearly satisfied and the result is obvious.

Assume than that Q is singular with rank m -ρ. Since Q is symmetric, there exists an orthogonal matrix U such that

U QU = Q 1 0 0 0 with Q 1 ∈ R (m-ρ)×(m-ρ) nonsingular.
Hence, without changing the essence of the problem, we assume that Q has the form

Q = Q 1 0 0 0 (7)
with Q 1 ∈ R (m-ρ)×(m-ρ) nonsingular, and hence Q 1 > 0.

We partition H and K(ω) accordingly, as follows:

H = H 1 H 12 H 21 H 2 , K(ω) = K 1 (ω) K 12 (ω) K 21 (ω) K 2 (ω) . (8) Thus [
F (jω) + F (-jω) ] 22 = 1 ω (H 2 + K 2 (ω) ω ). (9)
Since H = (F 1 -F 1 )/j, we have that: (i) H is Hermitian, so that H 2 is Hermitian as well, and (ii) the elements on the diagonal of H are zero, so that also the elements on the diagonal of H 2 are zero and hence H 2 is traceless. Therefore, H 2 has only real eigenvalues and the sum of the eigenvalues of H 2 is zero. Hence, either H 2 = 0 or H 2 has at least a negative eigenvalue. But for ω sufficiently large, and by continuity of the eigenvalues as functions of ω, the eigenvalues of H 2 + K2(ω) ω are arbitrarily close to those of H 2 , so that if H 2 = 0, then [F (jω)+F (-jω) ] 22 has a negative eigenvalue for a sufficiently large ω, and this is against our assumptions because [F (jω) + F (-jω) ] > 0 for all ω ∈ R so, in turn,

[F (jω) + F (-jω) ] 22 > 0 for all ω ∈ R. In conclusion, H 2 = 0. Now let Φ(ω) = Φ 1 (ω) Φ 12 (ω) Φ 21 (ω) Φ 2 (ω) = F (jω) + F (-jω) .
By continuity, as ω → ∞, m -ρ of the eigenvalues of Φ, i.e., the eigenvalues of Φ 1 (jω), tend to the eigenvalues of Q 1 (that are strictly positive) and the remaining ρ eigenvalues tend to zero. Let λ(ω) be one of the eigenvalues of Φ that tends to zero as ω → ∞. We now show that λ(ω) tends to zero at least as fast as 1/ω 2 . In fact, provided that ω is large enough so that λ(ω) is not an eigenvalue of Φ 1 (ω), then λ(ω) must be an eigenvalue of

R(ω) = Φ 2 (ω) -Φ 21 (ω)[Φ 1 (ω) -λ(ω)I] -1 Φ 12 (ω) (10)
because through Schur complements (e.g. [START_REF] Bernstein | Matrix Mathematics. Theory, Facts, and Formulas with Application to Linear Systems Theory[END_REF])

det[Φ(ω) -λ(ω)I] = det[Φ 1 (ω) -λ(ω)I] det[R(ω) -λ(ω)I].
Using the fact that

Φ 1 (ω) -λ(ω)I = Q 1 + H1 ω + K1(ω) ω 2
λ(ω)I, and using [1, Fact 9.9.43], R(ω) in (10) can be written equivalently as

R(ω) = Φ 2 (ω) -Φ 21 (ω)[Q -1 1 + ∆(ω)]Φ 12 (ω)
, where ∆(ω) → 0 as ω → ∞, since λ(ω) → 0 as ω → ∞. Finally, using (6) ( 7) and ( 8), we have

R(ω) = 1 ω 2 [K 2 (ω) -P (ω)] where P (ω) := (H 21 + K 21 (ω) ω )[Q -1 1 + ∆(ω)](H 12 + K 12 (ω) ω ).
Since for ω sufficiently large [K 2 (ω) -P (ω)] is bounded, then the eigenvalues of R(ω) tend to zero at least as fast as 1/ω 2 .

Then Φ(ω) has m strictly positive eigenvalues for all ω ∈ R, ρ of which are going down to zero at least as fast as 1/ω 2 as ω → ∞. The remaining (m -ρ) eigenvalues tend to the strictly positive eigenvalues of Q 1 as ω → ∞. Now, order the eigenvalues of Φ(ω) in non-decreasing size.

Since det(•) is the product of eigenvalues, then ω 2ρ det[Φ(ω)] is the same as the product of (ω 2 λ 1 ), . . . , (ω 2 λ ρ ), λ ρ+1 , . . . , λ m because for each i ∈ {1, . . . , ρ}, λ i tends to zero at least as fast as 1/ω 2 . Then, we have that the side condition (3b): In this case ρ = 1 and, by direct computation, we see that

lim ω→∞ ω 2ρ det[F (jω)+F (-jω) ] = lim ω→∞ 8ω 2 + 4ω 4 (1 + ω 2 ) 2 = 4 = 0.
On the other hand, we have

σ ω 2 F (jω) + F (-jω) = 2ω 2 1 + ω 2
which is clearly greater than 1 for all |ω| > 1.

The pathological case, corresponding to the situation in which some of the eigenvalues of the spectrum go to zero faster than 1 ω 2 , as ω tends to infinity, is more interesting: let (1 + ω 2 ) 3 = 0.

F (s) = s+2 (
The same conclusion is obtained with the other side condition which, in this case, reads:

σ ω 2 F (jω) + F (-jω) = 4ω 2 (1 + ω 2 ) 2 which is clearly not bounded away from zero as |ω| diverges. Finally, the trivial example F (s) = 1 0 0 -1 immediately shows that, if we do not assume the first two conditions of Proposition 2.1, the two side conditions are not necessarily equivalent.

III. CONCLUSION

This note shows that two different side conditions used in the control literature to characterize strictly positive real matrix transfer functions are equivalent.

  lim ω→∞ w 2ρ det[Φ(ω)] = 0 is equivalent to (since λ ρ+1 ,...,λ m do not tend to zero) lim ω→∞ (w 2 λ i ) = 0 for all i ∈ {1, ..., ρ}, which is equivalent to lim ω→∞ w 2 λ 1 > 0, which, finally, is equivalent to the first side condition (3a). The proposition is proved. Examples. Consider the case of F (s) = 1

.

  In this case ρ = 1 and by direct computation, we easily see thatlim ω→∞ ω 2ρ det[F (jω)+F (-jω) ] = lim ω→∞ ω 2 (16 + 8ω 2 )

This work was supported in part by the Engineering and Physical Sciences Research Council [grant number EP/R008876/1]. All research data supporting this publication are directly available within this publication.