Benjamin Elie
email: benjamin.elie@inria.fr

Yves Laprie

Extension of the single-matrix formulation of the vocal tract: consideration of bilateral channels and connection of self-oscillating models of the vocal folds with a glottal chink

Keywords: Speech synthesis, Vocal folds, Glottal chink, Lateral consonants

The paper presents extensions of the single-matrix formulation (Mokhtari et al., 2008, Speech Comm. 50(3) 179 -190) that enable self-oscillation models of vocal folds, including glottal chink, to be connected to the vocal tract.

They also integrate the case of a local division of the main air path into two lateral channels, as it may occur during the production of lateral con-

sonants. Provided extensions are detailed by a reformulation of the acoustic conditions at the glottis, and at the upstream and downstream connections of bilateral channels. The simulation framework is validated through numerical simulations. The introduction of an antiresonance in the transfer function due to the presence of asymmetric bilateral channels is confirmed by the simulations. The frequency of the antiresonance agrees with the theoretical predictions. Simulations of static vowels reveal that the behavior of the vocal folds is qualitatively similar whether they are connected to the

Introduction

Time-domain continuous speech synthesizers are commonly based on simplified physical models to compute the acoustic propagation along the vocal tract [START_REF] Kelly | Speech synthesis[END_REF][START_REF] Maeda | A digital simulation method of the vocal-tract system[END_REF][START_REF] Mokhtari | Single-matrix formulation of a time domain acoustic model of the vocal tract with side branches[END_REF][START_REF] Story | Phrase-level speech simulation with an airway modulation model of speech production[END_REF][START_REF] Birkholz | Enhanced area functions for noise source modeling in the vocal tract[END_REF] and/or the self-sustaining oscillations of the vocal folds [START_REF] Ishizaka | Synthesis of voiced sounds from a two-mass model of the vocal cords[END_REF][START_REF] Pelorson | Theoretical and experimental study of quasisteady-flow separation within the glottis during phonation. Application to a modified twomass model[END_REF][START_REF] Erath | A theoretical model of the pressure field arising from asymmetric intraglottal flows applied to a two-mass model of the vocal folds[END_REF][START_REF] Moisik | Modeling the biomechanical influence of epilaryngeal stricture on the vocal folds: A low-dimensional model of vocalventricular fold coupling[END_REF].

In comparison with finite-element-based methods [START_REF] Fleischer | Formant frequencies and bandwidths of the vocal tract transfer function are affected by the mechanical impedance of the vocal tract wall[END_REF], which require a huge amount of time, their low computation time makes them interesting for continuous speech synthesis.

Simplified acoustic models use the plane wave assumption to compute the acoustic propagation along a set of acoustic tubes. The dimensions of the elementary tubes (or tubelets) approximate the geometry of the vocal tract.

In regards to the typical dimensions of the human vocal tract, these models are valid up to frequencies around 5 kHz [START_REF] Story | Vowel acoustics for speaking and singing[END_REF].

Articulatory synthesis bridges the gap between the articulatory and acoustic domains of speech. This is thus an invaluable tool to apprehend the acoustic impact of the speech articulator's gestures and their temporal coordination, that of the anatomic characteristics of the human vocal tract, and of the interactions between the vocal folds and the vocal tract. In order to enable speech production to be studied via articulatory synthesis, several aspects should be covered by the numerical simulations of the speech aerodynamic/acoustic phenomena. First, the complexity of the vocal tract should be accurately modeled so that the various cavities (nasal tract, paranasal sinuses, sublingual cavities, etc) can be taken into account during the simulation. Then, the simulation framework should be able to deal with time-varying geometries of the vocal tract in order to simulate word-level or phrase-level utterances. This constrains the time trajectory of each articulator to be accurately modeled. Finally, the acoustic coupling between the glottal source, i.e. the vocal folds, and the vocal tract needs to be realistically modeled.

So far, there is no known time-domain continuous speech synthesizer that can deal with all these constraints. The scientific literature about speech synthesis based on simplified physical models identifies two main techniques: the reflection-type line analog method [START_REF] Kelly | Speech synthesis[END_REF], which is called RTLA in this paper, and the transmission line circuit analog [START_REF] Maeda | A digital simulation method of the vocal-tract system[END_REF] method, called TLCA.

The reader may find a detailed review of existing techniques for speech synthesis in [START_REF] Kröger | Articulatory synthesis of speech and singing: State of the art and suggestions for future research[END_REF]. Basically, RTLA has the advantage of accurately accounting for the frequency dependence of acoustic losses, but suffers from the constraints on the tubelet dimensions in regards to the chosen simulation frequency. As a consequence, the total length of the vocal tract cannot be modified during the simulation. This is an important issue for continuous speech synthesis since the length of the vocal tract varies during natural speech production. Its use is usually limited to studies about the self-sustained motion of the vocal folds [START_REF] Story | Voice simulation with a body-cover model of the vocal folds[END_REF][START_REF] Lous | A symetrical two-mass vocal-fold model coupled to vocal tract and trachea, with application to prothesis design[END_REF][START_REF] Bailly | Influence of a constriction in the near field of the vocal folds: Physical modeling and experimental validation[END_REF] coupled with simplified acoustic resonators. Using RTLA to simulate running speech constrains the vocal tract to unrealistic simplified geometries [START_REF] Story | Phrase-level speech simulation with an airway modulation model of speech production[END_REF]. On the other hand, many continuous speech synthesizers use TLCA [START_REF] Maeda | A digital simulation method of the vocal-tract system[END_REF][START_REF] Ishizaka | Synthesis of voiced sounds from a two-mass model of the vocal cords[END_REF][START_REF] Birkholz | Influence of temporal discretization schemes on formant frequencies and bandwidths in the time-domain simulation of the vocal tract system[END_REF][START_REF] Laprie | Comment faire parler les images aux rayons X du conduit vocal (How to make X-ray images speak)[END_REF]. It is based on the electric-acoustic analogy: the vocal tract acoustics is seen as a lumped electric circuit. The main advantage of TLCA is its flexibility of use with time-varying geometries of the vocal tract, including length variation and uneven spatial sampling of the vocal tract. However, this analogy does not allow the frequency dependence of the acoustic losses and the acoustic radiation to be accurately taken into account. Recently, the Single-Matrix Formulation (SMF) [START_REF] Mokhtari | Single-matrix formulation of a time domain acoustic model of the vocal tract with side branches[END_REF] has been a major contribution to TLCA models: it is now possible to compute the acoustic propagation along a vocal tract modeled as a waveguide network, where each waveguide represents a side cavity. Consequently, using SMF is a useful tool to study the acoustic effects of the numerous side cavities in the context of continuous speech synthesis.

Another important challenge to tackle when dealing with articulatory synthesis is the glottal source model. Indeed, in order to thoroughly study the phonatory mechanisms, one should include a glottis model that is able to realistically account for the coupling between the vocal folds and the vocal tract. Many efforts to simulate the production of the glottal source have been made, and self-oscillating models of vocal folds, based on lumped massspring systems, have been of particular interest. The reader may find detailed reviews of these models in [START_REF] Birkholz | A survey of self-oscillating lumpedelement models of the vocal folds[END_REF][START_REF] Erath | A review of lumped-element models of voiced speech[END_REF]. Curiously, most of these studies have neglected the possibility of including a posterior glottal opening to simulate air leakage. This glottal opening, also called glottal chink, allows a DC component to appear in the glottal flow waveform, as commonly observed in vivo [START_REF] Klatt | Analysis, synthesis, and perception of voice quality variations among female and male talkers[END_REF]. A partial closure of the glottis during the oscillating cycle of the vocal folds may be useful to simulate voiced fricatives [21] and breathiness, which is an important acoustic cue, especially for gender identification [START_REF] Klatt | Analysis, synthesis, and perception of voice quality variations among female and male talkers[END_REF].

First attempts at modeling a glottal chink have used parametric models [START_REF] Cranen | Modeling a leaky glottis[END_REF][START_REF] Cranen | Physiologically motivated modelling of the voice source in articulatory analysis/synthesis[END_REF]. In these papers, two models of glottal leakage are proposed: firstly, the glottal chink is caused by a partial abduction of the vocal folds, i.e. only a portion of the vocal folds vibrates, the other part is abducted and forms a triangular glottal chink, and secondly, the glottal chink is formed in the inter-arytenoid portion of the glottis. In the second case, the vocal folds vibrate along their whole length. Later on, Wilhelms-Tricarico [START_REF] Wilhelms-Tricarico | A modified two-mass model of the vocal folds with a chink and gradual closure, speech Communication Group[END_REF] proposed a modification of the classic two-mass model by Ishizaka and Flanagan [START_REF] Ishizaka | Synthesis of voiced sounds from a two-mass model of the vocal cords[END_REF] to include the glottal chink by connecting an electric branch in parallel to the vocal fold model. The glottal system is then connected to the first resonance of the vocal tract. Recently, Zañartu et al. [START_REF] Zañartu | Modeling the effects of a posterior glottal opening on vocal fold dynamics with implications for vocal hyperfunction[END_REF] have studied the effect of the glottal chink on self-oscillating movements of the vocal folds. Although this study is a significant advance in glottis modeling, its connection with RTLA models of acoustic propagation does not make it suitable for continuous speech synthesis with realistic time-varying geometries of the vocal tract.

Starting from the single-matrix formulation presented in [START_REF] Mokhtari | Single-matrix formulation of a time domain acoustic model of the vocal tract with side branches[END_REF], this paper details the theory and the methodology for extending it by overcoming its limitations. The limitations of SMF are the following: it does not offer the possibility to connect a self-oscillating model of the vocal folds, and the configuration of anastomozing waveguides, i.e. the local division of the main oral tract into two lateral channels, as it may occur during the production of lateral consonants, is not discussed. The aim is then to propose a complete simulation framework for speech synthesis that can account for the complexity of the vocal tract geometry and its numerous cavities taken simultaneously, that can deal with a time-varying realistic model of the vocal tract, including length variation, and that realistically models the acoustic coupling between the glottal source and the vocal tract, including a glottal leakage.

The paper is organized as follows. The transmission line circuit analog and the original single matrix formulation are detailed in Sec. 2. It also includes the required acoustic conditions at the glottis for integrating self-oscillating models of the vocal folds. The main aspects of the simulation framework, called Extended Single-Matrix Formulation (ESMF), are detailed in Sec. 3. They consist in the mathematical formulations for introducing the case of anastomozing waveguides into the single-matrix formulation, as well as the mathematical formulations to connect self-oscillating models of the vocal folds and a glottal chink to the single-matrix formulation. Finally, Sec. 4 presents numerical simulations that illustrate the accuracy of the extended single-matrix formulation to deal with the new features.

Theoretical background

This section summarizes the single-matrix formulation of the vocal tract by Mokhtari [START_REF] Mokhtari | Single-matrix formulation of a time domain acoustic model of the vocal tract with side branches[END_REF], which is itself derived from the transmission line circuit analog model by Maeda [START_REF] Maeda | A digital simulation method of the vocal-tract system[END_REF]. The present paper provides modifications in the formulation, taking into account the internal resistance of a noise source pressure to simulate frication noise. This reformulation is motivated by the fact that many quantities introduced in this section are used to demonstrate the contributions detailed in the next section. Yet, for the sake of brevity, not all computation details have been provided, and one may refer to the original papers [START_REF] Maeda | A digital simulation method of the vocal-tract system[END_REF][START_REF] Mokhtari | Single-matrix formulation of a time domain acoustic model of the vocal tract with side branches[END_REF] for more details.

Transmission line circuit analog model

The vocal tract is modeled as a concatenation of cylindrical tubes (or tubelets) along which plane waves propagate. The length and the crosssectional areas of the tubelets are such that they approximate the vocal tract geometry. TLCA represents each tubelet as lumped circuit elements. Fig. 1 shows the lumped circuit elements of a single tube section and Tab. 1 details the acoustic-electric analogy. The terms W R , W C , and W L are constant terms denoting respectively the resistance, the stiffness, and the mass of the vocal tract walls per area unit.

Chosen values for this study are those provided in [START_REF] Birkholz | Influence of temporal discretization schemes on formant frequencies and bandwidths in the time-domain simulation of the vocal tract system[END_REF], namely W R = 8000 kg.m -2 .s -1 , W C = 8.45 × 10 6 kg.m -2 .s -2 , and W L = 21 kg.m -2 . By convention,

(R i = 4πµl i a i ) C i Air compliance (C i = a i l i (ρcs) 2 ) L i Air inertance (L i = ρl i 2a i ) R w i Wall resistance (R w i = W R 2l i √ πa i ) C w i Wall compliance (C w i = 2l i √ πa i W C ) L w i Wall inertance (L w i = W L 2l i √ πa i ) U d i Flow source (-∂ ∂t l i a i ) P n i
Frication noise source (P n i = max{0, ξw

U 3 DC a 3/2 i-1 (Re 2 -Re 2 c )}) R n i Internal resistance of noise source (R n i = κρ U DC a 2 i-1 + 8πµ l i-1 a 2 i-1
) indices follow the air flow direction. For instance, considering the vocal tract, index 1 denotes the glottis connection, and index N denotes the lip termination.

Note that unlike in [START_REF] Maeda | A digital simulation method of the vocal-tract system[END_REF][START_REF] Mokhtari | Single-matrix formulation of a time domain acoustic model of the vocal tract with side branches[END_REF], lumped circuit elements include a frication noise source. It is made up of a pressure source P n i , with an internal resistance R n i [START_REF] Birkholz | Enhanced area functions for noise source modeling in the vocal tract[END_REF][START_REF] Maeda | Phoneme as concatenable units: VCV synthesis using a vocal tract synthesizer[END_REF][START_REF] Sondhi | A hybrid time-frequency domain articulatory speech synthesizer[END_REF], that becomes active when the air flow is considered as turbulent, namely when the Reynolds number Re is above a certain threshold Re c . Mathematical expressions for computing P n i and R n i are in Tab. 1, where ξ is an arbitrarily adjustable real constant used to control the noise level, and w is a Gaussian white noise to which first-order lowpass and third-order highpass filters have been applied [START_REF] Maeda | Phoneme as concatenable units: VCV synthesis using a vocal tract synthesizer[END_REF]. The chosen Re c threshold for the study is 2700, and ξ = 10 -6 , as suggested by Sondhi and Schroeter [START_REF] Sondhi | A hybrid time-frequency domain articulatory speech synthesizer[END_REF].

The frication noise source is usually located at the next point downstream of the supraglottal constriction [START_REF] Birkholz | Enhanced area functions for noise source modeling in the vocal tract[END_REF][START_REF] Maeda | Phoneme as concatenable units: VCV synthesis using a vocal tract synthesizer[END_REF][START_REF] Sondhi | A hybrid time-frequency domain articulatory speech synthesizer[END_REF], hence the spatial lag i -1. The term κ denotes a scaling coefficient set to 1.42, as used in [START_REF] Ishizaka | Synthesis of voiced sounds from a two-mass model of the vocal cords[END_REF][START_REF] Maeda | Phoneme as concatenable units: VCV synthesis using a vocal tract synthesizer[END_REF], and U DC is the low-frequency component of the air flow [START_REF] Maeda | Phoneme as concatenable units: VCV synthesis using a vocal tract synthesizer[END_REF]. Considering the lumped circuit elements displayed in Fig. 1, the continuous time equations relating the pressure and volume velocities inside all of the tubelets are

P i-1 -P i = ∂ ∂t [(L i-1 + L i ) U i ] + R i-1 + R i + R n i-1 U i + P n i-1 U i -U i+1 = u 1 + u 2 + u 3 , (1)
The paper follows the discrete representation of Eq. ( 1) introduced by Maeda [START_REF] Maeda | A digital simulation method of the vocal-tract system[END_REF], which yields to the following set of linear equations at any simulation time step:

             F 1 = Z 1 U 1 + b 1 U 2 F i = b i-1 U i-1 + Z i U i + b i U i+1 for 2 ≤ i ≤ N F N +1 = b N U N + Z N +1 U N +1 , (2)
where F i , b i , and Z i are pressure forces, loss terms, and impedance terms, associated with tubelet i, and N is the number of tubelets modeling the considered tract. One may refer to [START_REF] Maeda | A digital simulation method of the vocal-tract system[END_REF] for detailed steps leading to Eq. ( 2).

The acoustic propagation is computed by estimating the values of the volume velocity U i inside tubelets i = 1, . . . , N + 1 at each simulation time step. The previous system of equations forms a well-determined system, the N + 1 unknown volume velocities are governed by a set of N + 1 linear equations. It can be rewritten into the following matrix form

f = Zu, (3)
where

f ∈ R (N +1) = [F 1 , . . . , F N +1 ] T u ∈ R N +1 = [U 1 , U 2 , . . . , U N +1 ] T , Z ∈ R (N +1)×(N +1) =         Z 1 b 1 0 b 1 Z 2 b 2 0 0 . . . . . . . . . 0 b N Z N +1         .

Single-matrix formulation of the vocal tract

The SMF, as defined by Mokhtari et al. [START_REF] Mokhtari | Single-matrix formulation of a time domain acoustic model of the vocal tract with side branches[END_REF], is a reformulation of the TLCA model to be used with a vocal tract seen as a waveguide network.

The equations driving the acoustic propagation along every waveguide are merged into a single matrix. Each waveguide of the network represents a side cavity, represented as a parallel side branch in the analog lumped circuit.

In the waveguide network, the oral subcavity, going from the glottis to the lips, is called the root node. Quantities derived from the root node are denoted by the symbol (1) as exponent. Each waveguide connected to the oral tract is one of its children, and children of the root nodes may also have children themselves. The whole vocal tract can then be seen as a tree structure, where the root node is the oral subcavity. According to Mokhtari et al. [START_REF] Mokhtari | Single-matrix formulation of a time domain acoustic model of the vocal tract with side branches[END_REF], the system of equations driving the acoustic propagation inside N waveguides of the network is

        f (1)
f (2) . . .

f (N )         =         Z (1) C (2)T 1 . . . C (N )T 1 C (2) 1 Z (2) C (N )T 2 . . . . . . . . . C (N ) 1 C (N ) 2 Z (N )         .         u (1)
u (2) . . .

u (N )         , (4)
where f (m) , Z (m) , and u (m) contain the elements of Eq. ( 3) driving the acoustic propagation along the waveguide (m), independently of the rest of the network. Note that the number of waveguides N included in the network should not be confused with N , which denotes the number of tubelets included in a single waveguide. The number N may differ according to the waveguide it is related to.

The single-matrix formulation is then a concatenation of the linear systems driving the acoustic propagation inside each individual network, where sparse coupling matrices C

(n) m are added to account for the junction between two waveguides (m) and (n). In [START_REF] Mokhtari | Single-matrix formulation of a time domain acoustic model of the vocal tract with side branches[END_REF], the formulation provided by the authors distinguishes three cases: i) the waveguides are not directly connected, hence

C (n) m = 0, ii) (n) is a child of (m), connected at the K th tubelet of (m). Then C (n) m is a sparse matrix whose non-zero elements are c 1,K = b (m) K and c 1,K+1 = Z (m,n) C
, and iii) (m) and (n) are twins, both connected to the waveguide (p). This is the case for the piriform fossae, for instance: they are usually both connected to the sampled oral tract at the same tubelet. Then

C (n) m is a sparse matrix whose sole non-zero element is c 1,1 = Z (p,n) C .

Details of the mathematical expressions for computing C

(n) m may be found in [START_REF] Mokhtari | Single-matrix formulation of a time domain acoustic model of the vocal tract with side branches[END_REF]. The term Z (m,n) C accounts for the coupling between the pair of waveg-uides {m, n} at the junction K of the waveguide (m), and

Z (m,n) C = -b (m) K + R (m) K + R (m) n K + 2 T L (m) K .

Two-mass model and aeroacoustic considerations at the glottis

The method presented in this paper is intended to support any kind of self-oscillating model of the vocal folds. For the sake of brevity, only one will be used for the simulations as an example. We arbitrarily chose the 2 × 2-mass model with smooth contours [START_REF] Pelorson | Theoretical and experimental study of quasisteady-flow separation within the glottis during phonation. Application to a modified twomass model[END_REF][START_REF] Lous | A symetrical two-mass vocal-fold model coupled to vocal tract and trachea, with application to prothesis design[END_REF], but other models could be considered, such as the one-mass model [START_REF] Zañartu | Influence of acoustic loading on an effective single mass model of the vocal folds[END_REF], the three-mass model [START_REF] Tokuda | Comparison of biomechanical modeling of register transitions and voice instabilities with excised larynx experiments[END_REF],

or even models dealing with more masses [START_REF] Schwarz | Spatio-temporal quantification of vocal fold vibrations using high-speed videoendoscopy and a biomechanical model[END_REF]. It is based on two springmass systems, representing the rear and front ends of the vocal folds. It stems from the basic two-mass model by Ishizaka and Flanagan [START_REF] Ishizaka | Synthesis of voiced sounds from a two-mass model of the vocal cords[END_REF]. The model presented in this section considers recent improvements: contours are smooth, allowing a mobile separation point [START_REF] Pelorson | Theoretical and experimental study of quasisteady-flow separation within the glottis during phonation. Application to a modified twomass model[END_REF][START_REF] Lous | A symetrical two-mass vocal-fold model coupled to vocal tract and trachea, with application to prothesis design[END_REF], and it adds corrective terms to take into account the viscous losses and the unsteady flow effects [START_REF] Bailly | Influence of a constriction in the near field of the vocal folds: Physical modeling and experimental validation[END_REF][START_REF] Vilain | Experimental validation of a quasi-steady theory for the flow through the glottis[END_REF]. The geometry of the glottal constriction used in this paper is illustrated in Fig. 2.

The Bernoulli equation for unsteady flow, with an additive Poiseuille corrective term, gives the pressure P (x, t) along the glottal constriction:

P (x, t) = P sub + Be(x, t) + P o(x, t) + In(x, t) x < x s P (x, t) = P sup x > x s , (5)
where x s is the flow separation point, and Be(x, t), P o(x, t), and In(x, t)

are respectively the steady term of the Bernoulli equation, the Poiseuille corrective term and the unsteady term of the Bernoulli equation. They are defined as:

Be(x, t) = - ρU 2 g (t) 2l 2 g 1 h 2 (x, t) - 1 h 2 (x 0 , t) , P o(x, t) = - 12µU g (t) l g x x 0 dx h 3 (x, t) , In(x, t) = - ρ l g ∂ ∂t U g (t) x x 0 dx h(x, t) , (6)
where h(x) is the glottal opening along the x-coordinate, l g is the length of the vocal folds, ρ and µ are respectively the mass density and the shear viscosity of the air. The position of the flow separation point x s varies according to the glottal constriction geometry:

x s is such that h s = h(x s ) = 1.2h 1 if 1.2h 1 < h 2 , and x s = x 2 otherwise.
The value 1.2 is an ad-hoc criterion, as used in [START_REF] Pelorson | Theoretical and experimental study of quasisteady-flow separation within the glottis during phonation. Application to a modified twomass model[END_REF][START_REF] Lous | A symetrical two-mass vocal-fold model coupled to vocal tract and trachea, with application to prothesis design[END_REF]. At the flow separation point x = x s , the pressure drop between the upstream and the downstream parts of the glottis is given by Eq. ( 7)

P sub (t) -P sup (t) = R b (t)U 2 g (t) + R v (t)U g (t) + ∂ ∂t [L g (t)U g (t)] , (7)
where

R b (t) = ρ 2l 2 g 1 h 2 (x s , t) - 1 h 2 (x 0 , t) , R v (t) = 12µ l g xs x 0 dx h 3 (x, t) , L g (t) = ρ l g xs x 0 dx h(x, t) , (8)
From the determination of the glottal flow U g , Eq. ( 5) gives the pressure distribution P (x) along the glottal constriction. The pressure forces are then used to derive the mass positions at each simulation step following the classic system of differential equations

Mÿ + R ẏ + Ky = F, (9)
where

M ∈ R 4×4 + = diag (m 1u , m 2u , m 1d , m 2d ), R ∈ R 4×4 + = diag (r 1u , r 2u , r 1d , r 2d ),
and

F ∈ R 4×4 = diag (F 1u , F 2u , F 1d , F 2d
) are diagonal matrices containing the values of respectively the mass, the damping and the pressure forces applied to each mass, y ∈ R 4 = [y 1u , y 2u , y 1d , y 2d ] T is the vector containing the displacement of each mass from its rest position, and K ∈ R 4×4 is a matrix containing stiffness coefficients. Due to the presence of a coupling spring k c ,

K =         k 1u + k cu -k cu 0 0 -k cu k 2u + k cu 0 0 0 0 k 1d + k cd -k cd 0 0 -k cd k 2d + k cd        
Pressure forces F i,j (t), with j = {u, d}, are derived from the pressure applied to the mass at instant t:

F 1,j = l g x 1 x 0 x -x 0 x 1 -x 0 P (x)dx + l g x 2 x 1 x -x 2 x 1 -x 2 P (x)dx F 2,j = l g x 2 x 1 x -x 1 x 2 -x 1 P (x)dx + l g x 3 x 2 x -x 3 x 2 -x 3 P (x)dx

Extended Single-Matrix Formulation of the Vocal Tract

General remark:

In this paper, the subglottal pressure is imposed as an input parameter.

Although its temporal evolution should be considered for time-domain continuous speech synthesis, especially in order to simulate natural prosody, it has not been taken into account in this paper since the compatibility with the single-matrix formulation has already been proven, even for complex geometries, by Ho et al. [START_REF] Ho | An anatomically based, timedomain acoustic model of the subglottal system for speech production[END_REF].

Anastomozing waveguides: bilateral consonants

In some cases, the air path inside the vocal tract may be locally divided into two lateral channels. For instance, this is observed in the case of lateral consonants [START_REF] Narayanan | Toward articulatory-acoustic models for liquid approximants based on mri and epg data. part i. the laterals[END_REF][START_REF] Zhang | Acoustic modeling of American English lateral approximants[END_REF][START_REF] Zhang | A vocal-tract model of american english /l[END_REF]. The acoustic effect is not fully apprehended, mainly because of the lack of appropriate acoustic models, and also because of insufficient relevant articulatory data. Indeed, lateral consonants are usually modeled by a single lateral channel, and additionally a supralingual cavity [START_REF] Stevens | Acoustic Phonetics[END_REF][START_REF] Narayanan | Geometry, kinematics, and acoustics of tamil liquid consonants[END_REF]. Zhang et al. [START_REF] Zhang | Acoustic modeling of American English lateral approximants[END_REF][START_REF] Zhang | A vocal-tract model of american english /l[END_REF] studied the effect of bilateralization with a frequency-based method to compute the transfer function of bilateral vocal tracts. However, to the best of our knowledge, there is still no existing model of bilateral vocal tracts to be used in the context of time-domain continuous speech synthesis. This section details the theoretical and mathematical developments for integrating this condition into the single-matrix formulation.

The aim is then to define the mathematical expression of the coupling matrix

C (n)
m of Eq. ( 4) when (m) and (n) are two parallel waveguides, connected to each other at two different points. Following the medical and the hydrography terminology, this relation is called anastomosis in this paper. This relationship distinguishes a main anastomozing waveguide and a secondary one, called anabranch. For instance, the main anastomozing waveguide in Fig. 3 is waveguide (1), while (n) is its anabranch.

It is clear that the conditions at the upstream connection are similar to a parent-child relationship, hence similar expressions in the coupling matrix

C (n)
m . In that case, three equations of the system, corresponding to the three tubelets around the junction, are modified:

F (1) Ku = b (1) Ku-1 U (1) Ku-1 + Z (1) Ku U (1) Ku + b (1) Ku U (1) Ku+1 + b (1) Ku U (n) 1 , F (1)
Ku+1 = b (1) Ku U (1) Ku + Z (1) Ku+1 U (1) Ku+1 + b (1) Ku+1 U (1) Ku+2 + Z (1,n) C U (n) 1 , F (n) 1 = Z (n) 1 U (n) 1 + b (n) 1 U (n) 2 + b (1) Ku U (1) Ku + Z (1,n) C U (1)
Ku+1 . [START_REF] Fleischer | Formant frequencies and bandwidths of the vocal tract transfer function are affected by the mechanical impedance of the vocal tract wall[END_REF] Similarly, at the downstream connection K d , the equations corresponding to the tubelets around the junction are modified:

F (1) K d +1 = b (1) K d U (1)
K d + Z (1) K d +1 U (1) K d +1 + b (1) K d +1 U (1)
K d +2 + Z (1,n) C+1 U (n) N +1 , F (1)
K d +2 = b (1) K d +1 U (1) K d +1 + Z (1) K d +2 U (1) K d +2 + b (1) K d +2 U (1) K d +3 + b (1) K d +1 U (n) N +1 , F (n) N +1 = b (n) N U N (n) + Z (n) N +1 U (n) N +1 + b (1) K d +1 U (1) K d +2 + Z (1,n) C+1 U (1) K d +1 , (11)
where

Z (1,n) C+1 = -b (1) K+1 + R (1) K+1 + 2 T L (1) K+1 , Z (n)
N +1 = -b (n) N -b (1)
K+1 -

2 T L (n) N + L (1) K+1 -R (n) N -R (1)
K+1 -R (n) n N . (12)
Note that in that case, the expression of Z

N +1 differs from other configurations because of different boundary conditions at the waveguide termination.

The matrix formulation is then

  f (1) f (n)   =   Z (1) C (n)T 1 C (n) 1 Z (n)   .   u (1) u (n)   , (13)
where C

(n)

1 is the matrix accounting for the coupling between the anastomozing waveguides. It is a sparse matrix with 4 non-zero elements c 1,Ku = b

(1) Ku , c 1,Ku+1 = Z (1,n) C , c N +1,K d +2 = b (1) K d +1 , and c N +1,K d +1 = Z (1,n)
C+1 , where N is the generic number of tubelets that model the anabranch (n), and K u and K d are the upstream and downstream locations on the main oral tract where (n) is connected.

Integration of the self-oscillating model of vocal folds

In the original papers [START_REF] Maeda | A digital simulation method of the vocal-tract system[END_REF][START_REF] Mokhtari | Single-matrix formulation of a time domain acoustic model of the vocal tract with side branches[END_REF], the glottal source is generated by an imposed oscillating glottal input area. This may be sufficient to simulate utterances with a good quality, but it cannot be used to study the acoustic coupling between the vocal folds and the vocal tract. In [START_REF] Ho | An anatomically based, timedomain acoustic model of the subglottal system for speech production[END_REF], the authors connect the original two-mass model by Ishazaka and Flanagan [START_REF] Ishizaka | Synthesis of voiced sounds from a two-mass model of the vocal cords[END_REF] to the single-matrix formulation, in order to simulate the self-oscillating motion of the vocal folds.

However, the quadratic term in Eq. ( 7) is not taken into account. The present section details the mathematical considerations to account for a more realistic aeroacoustic model at the glottis when the single-matrix formulation is used in the context of time-domain speech synthesis.

To connect the self-oscillating model of the vocal folds, the pressure distribution along the glottal constriction should be known at each simulation step. Thus, Eq. ( 7) should be integrated into the single-matrix formulation.

Once U g is known, pressure forces are derived from Eq. ( 5), and the mass positions are computed following Eq. ( 9). Introducing Eq. ( 7) into the single-matrix formulation of Eq. ( 22) requires the first line of the system to be modified into a quadratic equation

F 1 = Z 1 U 1 + b 1 U 2 + R b U 2 1 . (14)
The matrix form is then

f = Zu Z + Qu Q , ( 15
)
where Q is a square matrix the same size as Z having only one non-zero element, i.e. Q (1,1) = R b , and

u Q ∈ R (N +1) = [U 2 1 , U 2 2 , . . . , U 2 N ]
T is the vector containing the square power of the volume velocities. Eq. ( 15) is also valid in the case of a waveguide network, since it does not directly modify the coupling equations between the different side cavities modeling the VT.

The system is almost entirely linear: only the first line is a quadratic equation. To solve the system, one should first solve the quadratic equation separately. It could be straightforward if the first line of Z contained only one non-zero element. Unfortunately, it is not the case since Z (1,1) = 0 and Z (1,2) = 0. A practical solution consists in finding an equivalent system in which the matrix of linear coefficients is a diagonal matrix. Left-multiplying both sides of Eq. ( 15) by Z -1 yields

Z -1 f = Iu Z + Z -1 Qu Q , ( 16
)
where I is the identity matrix. In this formulation, the first line of the system depends only on U 1 . The quadratic equation can then be solved, and the value U 1 that is accepted is the largest positive solution. If both roots are negative, U 1 is set to 0. Once the glottal volume velocity U 1 is known, Eq. ( 15) can be rearranged

f = Zũ Z , (17)
where

f ∈ R N = [f 2 -b 1 U 1 , f 3 , . . . , f n , . . . , f N +1 ] T , Z ∈ R N ×N is the matrix
Z of which the first line and first column have been withdrawn, and ũZ ∈

R N = [U 2 , U 3 , . . . , U n , . . . , U N +1 ] T .
Eq. ( 17) is then a well-determined tridiagonal linear system. Any classical method to solve such systems gives the solutions U i with i = 2, . . . , N + 1.

A glottal chink model

The self-sustaining model of the vocal folds presented in the previous section considers them to vibrate uniformly along their length. Consequently, at this stage, it is not possible to account for a partial closure of the glottis.

This may be an issue for the synthesis of several types of phonation, including breathy voice and voiced fricatives. In these cases, the glottis is never completely closed, and an offset appears in the glottal flow waveform, so that it is never null.

The model presented in this paper connects an electric branch in parallel to the vocal fold model. The electric analogy and the parallel branch make this model interesting for the extended single-matrix formulation. Hence the equivalent electric circuit represented in Fig. 4. Unlike in [START_REF] Wilhelms-Tricarico | A modified two-mass model of the vocal folds with a chink and gradual closure, speech Communication Group[END_REF], U g and U ch separate upstream of the glottal contraction, and merge downstream of the glottal expansion. These assumptions follow the model detailed in [START_REF] Cranen | Physiologically motivated modelling of the voice source in articulatory analysis/synthesis[END_REF].

The formulation assumes that the partial closure of the glottis is due to a partial abduction of the vocal folds and that the chink is linked to the membranous part of the glottis, as shown in Fig. 4. Since the model considers one-dimensionnal inviscid flow, the chink and the oscillating part of the vocal folds may be represented as independent channels.

The glottal chink area is a ch = l ch h ab , where l ch is the length of the glottal chink and h ab is the abduction of the vocal folds. Note that in the original model [START_REF] Cranen | Physiologically motivated modelling of the voice source in articulatory analysis/synthesis[END_REF], the glottal chink may be extended with a constant opening area, due to the inter-arytenoid portion of the glottis. The presented model does not account for it, but the implementation is straightforward: the glottal chink area is then a ch = l ch h ab + a ia , where a ia is the opening area of the inter-arytenoid portion of the glottis. When the glottal chink is considered, U 1 is no longer connected to U g . Indeed, as seen in Fig. 4 b), U 1 is the merged flow downstream the glottal expansion, and is therefore the sum of the volume velocities through the chink and through the vibrating part of the vocal folds (U 1 = U g + U ch ).

The presence of the glottal chink modifies the first line of the system defined in Eq. ( 16). Indeed, applying Eq. ( 7) to both glottal branches and reorganizing equations yields to the following boundary conditions at the glottis

F 1 = R b U 2 g + Z 1 U g + b 1 U 2 + Z (1,ch) C U ch , F 2 = b 1 U g + Z 2 U 2 + b 2 U 3 + b 1 U ch , F ch = Z ch U ch + Z (1,ch) C U g + b 1 U 2 , (18)
where

Z (1,ch) C = b 1 + R 1 + 2 T L 1 (19)
Z ch = b 1 + R ch + R 1 + 2 T (L ch + L 1 ). (20)
Introducing Eqs. [START_REF] Birkholz | A survey of self-oscillating lumpedelement models of the vocal folds[END_REF] in the single matrix formulation yields

  f (1)
F ch   =   Z (1) C (ch)T 1 C (ch) 1 Z ch   .   u (1)
U ch   + R b U 2 g , (21)
where C

(ch) 1

is the matrix accounting for the coupling between the glottal chink and the vocal tract. It is a sparse row vector with two non-zero elements:

C (ch) 1 = [Z (1,ch) , b 1 , 0, . . . , 0].

General form of the extended single-matrix formulation of the vocal tract

Finally, the general form of the extended single-matrix formulation (ESMF) of the vocal tract is

           f (1)
f (2) . . .

f (N ) F ch            =            Z (1) C (2)T 1 . . . C (N )T 1 C (ch)T 1 C (2) 1 Z (2) C (N )T 2 0 . . . . . . . . . C (N ) 1 C (N ) 2 Z (N ) 0 C (ch) 1 0 . . . 0 Z ch            .            u (1)
u (2) . . .

u (N ) U ch            + Qu Q (22)
where N is the number of waveguides modeling the whole vocal tract. In addition to the single-matrix formulation defined in Eq. ( 4), two new configurations detailed in this section are considered: i) (n) is an anabranch of (m),

i.e. it is connected to (m) at two different points, the upstream connection K u and the downstream connection

K d . Then C (n) m is a sparse matrix whose non-zero elements are c 1,Ku = b (m) Ku , c 1,Ku+1 = Z (m,n) C , c N +1,K d +2 = b (m) K d +1 , and c N +1,K d +1 = Z (m,n) C+1
, and where N is the generic number of tubelets that model (n), and ii) (n) is the glottal chink. In that case, it is connected to

Numerical simulations

In order to validate the simulation framework presented in the previous section, this section provides examples of synthesized vowels in static configurations. Area functions are derived from X-ray films comprising short French sentences [START_REF] Sock | DOC-VACIM an X-ray database and tools for the study of coarticulation, inversion and evaluation of physical models[END_REF]. They are obtained by dividing the vocal tract shape into tubelets perpendicular to the vocal tract centerline, determined via a specified algorithm [START_REF] Laprie | Articulatory copy synthesis from cine X-ray films[END_REF], and then applying α β transformations to recover the area [START_REF] Soquet | Mid-sagittal cut to area function transformations: Direct measurements of mid-sagittal distance and area with MRI[END_REF]. Chosen parameters of the vocal folds model (see Tab. 2) lie in the typical range of values found in the literature [START_REF] Pelorson | Theoretical and experimental study of quasisteady-flow separation within the glottis during phonation. Application to a modified twomass model[END_REF][START_REF] Sciamarella | On the acoustic sensitivity of a symmetrical two-mass model of the vocal folds to the variation of control parameters[END_REF][START_REF] Bailly | Vocal fold and ventricular fold vibration in period-doubling phonation: Physiological description and aerodynamic modelinga)[END_REF]. In the following numerical simulations, the vocal folds are supposed symmetric, namely the upper vocal fold has the same mechanical parameters as the lower one.

Note that when a waveguide is closed at its end, e.g. the piriform fossae, a termination area set to 0 may cause the numerical computation to break because of infinite or NaN values. Setting the termination area to a very small value, for instance the order of magnitude of the unit roundoff, is sufficient to efficiently approximate a closed termination. A similar technique can be used when sudden modifications of the global structure of the network occur, such as a sudden occurrence of nasalization or bilateralization: setting the input area function of a certain waveguide to very small values makes it temporally shunted from the network. Doing so prevents discontinuities that may cause undesired artifacts.

Anastomozing waveguides

The acoustics of lateral consonants has been previously studied mainly by using single lateral channel models [START_REF] Stevens | Acoustic Phonetics[END_REF][START_REF] Narayanan | Geometry, kinematics, and acoustics of tamil liquid consonants[END_REF]. Consequently, only a few previous works focused on the acoustic effect of bilateralization. To validate the bilateral model presented in this paper, it is compared to results provided by Zhang et al. [START_REF] Zhang | Acoustic modeling of American English lateral approximants[END_REF][START_REF] Zhang | A vocal-tract model of american english /l[END_REF], using the same area functions. Fig. 5 shows the singletube model used for the simulation in [START_REF] Zhang | A vocal-tract model of american english /l[END_REF]. The area functions of the different parts are derived from MRI [START_REF] Zhang | Acoustic modeling of American English lateral approximants[END_REF], and are averaged over the length of each section of the single-tube model (see Fig. 5). Dimensions can be found in the original paper [START_REF] Zhang | Acoustic modeling of American English lateral approximants[END_REF]. When the supralingual cavity is taken into account, the latter is connected to the main oral tract at the same location as the lateral channel. Consequently, the supralingual cavity is a twin waveguide of the lateral channel.

From the area functions provided by the original paper, transfer functions of the vocal tract are computed thanks to the method used in [START_REF] Maeda | A digital simulation method of the vocal-tract system[END_REF][START_REF] Mokhtari | Single-matrix formulation of a time domain acoustic model of the vocal tract with side branches[END_REF], which consists in simulating a sudden glottal closure and computing the Fourier transform of the acoustic response to the step-down glottal excitation. The vocal tract acoustic response functions are computed for different configura- tions of bilateralization. They are defined by an asymmetry factor

γ = l c 2 l c 1 ,
where l c 2 and l c 1 denote the length of both lateral channels. During the simulations, l c 1 is set to 4.2 cm, and l c 2 is increased to raise the asymmetric factor γ. Simulations include two main configurations: with and without the consideration of the supralingual cavity.

Fig. 6 shows transfer functions for several configurations obtained with ESMF. The asymmetric factor γ varies from 1 (bottom curve), corresponding to the symmetric configuration, to 1.37 (top curve). The increment step of γ between two successive curves is 0.01. The thick line at the bottom is the transfer function of the vocal tract where both lateral channels are merged into a single channel, and without the consideration of the supralingual cavity. The computed transfer functions agree with those obtained in [START_REF] Zhang | A vocal-tract model of american english /l[END_REF] with an independent frequency-based method. For instance, in both configurations (with and without the suparalingual cavity), the bilateralization introduces a pole/zero around 4 kHz. The frequency of the pole/zero pair drops as the length of the anabranch increases. This is in agreement with the theory, which relates the frequency of the pole/zero introduced by the bilateralization with the resonance frequency of an equivalent tube having a length equal to the sum of the lateral channels [START_REF] Prahler | Analysis and synthesis of the American English lateral consonant[END_REF], hence a drop of the frequency as the length increases. The theoretical zero frequencies are represented in Fig. 6 by circle marks: they match the zero frequencies of the simulated vocal tract acoustic response functions.

The length asymmetry of the bilateral channels slightly impacts the formant frequencies. F 5 is the formant for which the effect is the most predominant. It seems to be due to its proximity to the introduced pole/zero pair.

Finally, the introduction of the supralingual cavity in the model gives rise to another pole/zero pair slightly above 2.7 kHz. As observed in [START_REF] Zhang | A vocal-tract model of american english /l[END_REF], this corresponds to the first quarter-wavelength resonance of a 3.1 cm-long tube, around 2700 Hz. The supralingual cavity lowers the formant frequencies, and especially F 3 and F 4 , which are close to the zero introduced by the supralingual cavity.

The agreement between the global effects of bilateralization, and the supralingual cavity on the transfer function obtained with ESMF and those obtained with the independent frequency based technique used in [START_REF] Zhang | A vocal-tract model of american english /l[END_REF], proves the efficiency of the method to deal with bilateralization. Since it is suitable for time-domain continuous speech synthesis, the method could be useful to thoroughly investigate the acoustic effect of the tongue movement during the production of bilateral consonants.

Effect of the acoustic model on the motion of the vocal folds

In this section, numerical simulations are used to study the validity of ESMF when connected to a self-sustaining model of the vocal folds. The method is compared with the concurrent approach using reflection-type line analog (RTLA) models [START_REF] Kelly | Speech synthesis[END_REF][START_REF] Meyer | A quasiarticulatory speech synthesizer for german language running in real time[END_REF]. This approach is widely used and validated when connected with self-oscillating models of the vocal folds, including the 2 × 2-mass model with smooth contours used in this paper [START_REF] Lous | A symetrical two-mass vocal-fold model coupled to vocal tract and trachea, with application to prothesis design[END_REF][START_REF] Bailly | Influence of a constriction in the near field of the vocal folds: Physical modeling and experimental validation[END_REF][START_REF] Vilain | Experimental validation of a quasi-steady theory for the flow through the glottis[END_REF]. Thus, simulations consist in computing the motion of the vocal folds via the model described in Sec. /i/, /e/, /A/, /ø/, /o/, and /u/. Fig. 7 shows the spectrum and the spectral envelope of the synthesized output pressure (P Out ) radiated at the lips, as well as the glottal flow and the glottal opening at the location of mass 1 and mass 2, both by RTLA and ESMF, for 2 French vowels (/a/ and /e/). For each vowel, the peaks of the spectral envelopes, i.e. the formants, obtained from both methods are similar. One can also observe similarities in the global shape of the other waveforms, namely P Out , U g , and the motion of the vocal folds.

However, the comparison also highlights the differences between both approaches. For instance, the obtained fundamental frequency is lower with ESMF than with RTLA. The amplitude of the vocal folds' motion, and the phase shift between both masses are also lower with ESMF (see Fig. 8).

This yields to a lower amplitude of the glottal flow. Such differences are somewhat expected since the coupling with the acoustic propagation may be 29 The right column displays the same quantities computed thanks to ESMF. The spectral envelope is computed using a 10-order LPC (Linear Predictive Coding) [START_REF] Mccandless | An algorithm for automatic formant extraction using linear prediction spectra[END_REF] modified, due to the different models of acoustic losses. It is accepted that RTLA has the advantage of better accounting for the frequency dependence of the acoustic losses and the acoustic radiation. This may also explain the differences between formant bandwidths obtained with RTLA and ESMF.

It is worth noting that slight modifications of the input parameters of

Conclusions

This paper has presented the theoretical aspects for extending the singlematrix formulation [START_REF] Mokhtari | Single-matrix formulation of a time domain acoustic model of the vocal tract with side branches[END_REF] of the vocal tract to a more general and complete tool. The presented framework allows the connection of self-oscillating models of the vocal folds, or the connection of a glottal system made up of self-oscillating vocal folds and a glottal chink. It can also account for anastomozing waveguides to study the acoustic effect of bilateralization.

The accuracy of the simulation framework to account for bilateral channels has been studied by computing transfer functions of vocal tract configurations derived from a previously published paper [START_REF] Zhang | A vocal-tract model of american english /l[END_REF]. The effects of the bilateralization on the simulated transfer functions quantitatively agrees with the theoretical predictions [START_REF] Zhang | A vocal-tract model of american english /l[END_REF][START_REF] Prahler | Analysis and synthesis of the American English lateral consonant[END_REF]. Acoustic features of lateral consonants, such as the introduction of an antiresonance due to asymmetric bilateral channels, may thus be reproduced using ESMF.

The the developed approach is compatible with the subglottal network proposed in [START_REF] Ho | An anatomically based, timedomain acoustic model of the subglottal system for speech production[END_REF]. Combining the subglottal and supraglottal networks would be a step forward in the implementation of a complete and realistic speech synthesizer.

This would create a useful tool to study the phenomena involved in speech production. This paper may be considered as a basis to make such investigations. More specifically, it may be used to relate acoustic cues of the produced speech signal to their articulatory or phonatory origins, thanks to analysis-by-synthesis techniques. This could be a great benefit for phonetic sciences, and/or language training.

Figure 1 :

 1 Figure 1: Lumped circuit elements of an acoustic tube. The acoustic-electric analogy is detailed in Tab. 1.

Figure 2 :

 2 Figure 2: Geometry of the glottal constriction, as introduced by Lous et al. [14]. Indices u and d stand for up and down respectively.

Fig. 3

 3 Fig.3shows the waveguide connections in the case of bilateralization. It includes a secondary waveguide that is connected to the oral branch at two points.

Figure 3 :

 3 Figure 3: Equivalent diagram of anastomozing waveguides.

Figure 4 :

 4 Figure 4: a) View of the partially closed glottis, extracted from Cranen and Schroeter [22]. In this model, the partial closure is due to a partial abduction of the vocal folds. l g is the length of the vibrating part of the vocal folds, l ch is the length of the glottal chink and l t = l g + l ch is the total length of the vocal folds. The abduction of the vocal folds is assumed to be constant and is denoted by h ab . b) Electric-circuit analogy of the partially closed glottis. U ch , R ch , and L ch are the volume velocity through the glottal chink, the energy loss, and the air inertance inside the glottal chink, respectively.

( 1 )

 1 and C (ch) 1 = [Z (1,ch) , b 1 , 0, . . . , 0].

Figure 5 :

 5 Figure 5: Single-tube model used for the simulation of bilaterals. Figure adapted from [35].

Figure 6 :

 6 Figure 6: Transfer functions of the vocal tract in several configurations using the extended single-matrix formulation, with (right) and without (left) consideration of a supralingual cavity. The bottom curve (thick line) is the transfer function of the vocal tract where both lateral channels are merged into a single channel and without the supralingual cavity. The transfer functions of the vocal tract with bilateralization are plotted from bottom to top, where the asymmetric length factor varies from 1 (bottom) to 1.37 (top). The increment step of the length asymmetry factor between two successive curves is 0.01. Theoretical positions of zeros introduced by the bilateralization are denoted by circle marks.

  2.3 connected to several configurations of the vocal tract, using two models of acoustic propagation: ESMF and RTLA. The configurations of the vocal tract are area functions corresponding to 6 French vowels:

(

  FigureThe results of the simulations for 2 French vowels, /a/ and /e/. From top to bottom, the left column shows the spectrum and spectral envelope of the synthesized vowel via RTLA, as well as the output pressure (P Out ) radiated at the lips, the glottal flow U g and the glottal opening at the location of mass 1 (solid line), and mass 2 (dashed line).

  Fig. 8 shows that, both methods simulate qualitatively similar oscillations of the vocal folds and glottal flow, they are quantitatively different. The main differences lie in the phase shift and in the open quotient: when connected to RTLA, the phase shift between the rear and front parts of the vocal folds is larger than when they are connected to ESMF. There is no significant when the values of the mass and stiffness are modified. The open quotient is also significantly different when the vocal folds are connected to RTLA. However, in this case, the modification of the input parameters modifies the open quotient. Also, for both techniques, the ampli-

Figure 8 :

 8 Figure 8: Comparison of quantities derived from the motion of the vocal folds obtained with the different methods. ESMF* corresponds to the quantities obtained with modified mass and stiffness of the vocal folds model.

Figure 9 :

 9 Figure 9: Formant frequency of the first 4 formants obtained with the different methods for 6 French vowels. Results are compared with values obtained with the chain-matrix paradigm (CMP)[START_REF] Sondhi | A hybrid time-frequency domain articulatory speech synthesizer[END_REF] for the corresponding area functions.

Fig. 9

 9 Fig.9shows the effect of the methods on the formant frequencies. The values obtained with the different methods are compared with values obtained from the transfer function of the corresponding vocal tracts using the chain-matrix paradigm (CMP)[START_REF] Sondhi | A hybrid time-frequency domain articulatory speech synthesizer[END_REF]. This frequency-based method is used as

Figure 10 :

 10 Figure10: The results of the simulation for two French vowels /a/ (left column), and /i/ (right column). From top to bottom, the left column represents the narrow-band spectrogram, the output pressure, radiated at the lips, the volume velocity inside the glottis (U t = U g + U ch ), the glottal opening at the location of mass 1 (solid line) and mass 2 (dashed line), and the length l ch of the glottal chink, computed for /a/. The right column represents these quantities computed for /i/.

Figure 11 :

 11 Figure 11: Details of phonatory cycles of the simulated sustained /i/ at two time instants: (left) closed glottal chink, and (right) open glottal chink, l ch = 0.5 cm. From top to bottom: P Out , opening height of the vocal folds, and U t .

  connection of the self-oscillating model of the vocal folds has been tested with a 2 × 2-mass model with smooth contours. In comparison with the widely accepted reflection-type line analog model of acoustic propagation, the connection of the vocal folds to the single-matrix formulation yields to qualitatively similar behaviors of the vocal folds, though it may quantitatively modify their oscillations. However, at this stage, it is not possible to conclude whether these differences are important for quantitative studies. Simulations also revealed that the formant frequencies of the vowels obtained with ESMF are closer to the resonance frequencies of the vocal tract than those of the vowels simulated via the reflection-type line analog model. The paper also provides simulations with various lengths of the glottal chink. They clearly show the acoustic effects, including the appearance of a DC component in the glottal flow waveform. As the glottal chink becomes larger, this generates a large DC acoustic flow inside the vocal tract, which may lead to the generation of a frication noise in addition to the voiced glottal source if there is a supraglottal constriction in the vocal tract. This feature is important to simulate voiced fricatives.The simulation framework is intended to be used for the synthesis of natural continuous speech, i.e. natural phrase-level utterances. For instance, it has been used to validate a two-dimensional model of the velum for the copy synthesis of utterances containing nasal phonemes[START_REF] Laprie | 2D articulatory velum modeling applied to copy synthesis of sentences containing nasal phonemes[END_REF]. The transmission line circuit analog model enables the dynamic variations of the vocal tract geometry and its complexity to be accurately taken into account. Consequently, it can easily support realistic geometries and dynamic deformations of the vocal tract, which constitutes its main advantage. Thanks to the contributions of the paper, it is now possible to account for a realistic coupling between the vocal folds and the vocal tract, as well as the possibility of including glottal leakage. This last point allows the different types of phonation to be realistically simulated. Although it has only been briefly discussed in the paper,

Table 1 :

 1 Acoustic-electric analogy

	Electric	Acoustic
	Current	Volume velocity u
	Voltage	Acoustic pressure p
	R i	Energy loss

Table 2 :

 2 Input parameters for the vocal folds model

	Parameter	Unit	Value
	Subglottal pressure P sub	Pa	800
	Position of mass 1 x 1	mm	0.2
	Position of mass 2 x 2	mm	3.2
	Vocal fold thickness d g	mm	3
	Vocal fold length l g	mm	10
	Opening at point 0 h 0	mm	40
	Vocal folds abduction h ab	mm	2.5
	Nominal mass m 1	g	0.1
	Nominal stiffness k 1	N/m	80
	Nominal mass m 2	g	0.125
	Nominal stiffness k 2	N/m	80
	Nominal damping coefficient r Coupling spring k c	N/m	k/2

i kg.rad.s -1 0.2 k i m i /2

Acknowledgements

This study is supported by the ANR (Agence Nationale de la Recherche) ArtSpeech project. The authors would like to sincerely thank Dr. Shinji Maeda and Dr. Parham Mokhtari, for their useful advice and fruitful discussions.

Effect of the glottal chink

This section presents a short study about the effect of the glottal chink on the acoustic parameters. It consists in computing the glottal flow and the motion of the vocal folds coupled with a static configuration of the vocal tract and a linearly increasing length of the glottal chink. We chose to represent the effect of the glottal chink on /a/ and /i/. To modify the size of the glottal chink, we modify the length l ch via a ratio α ∈ [0,1] such that l ch = αl t (see Fig. 4). The quantities l t and h ab are kept constant, and by deduction, l g = (1 -α)l t . Fig. 10 shows the simulation results. It displays the narrow-band spectrogram, the output pressure radiated at the lips P Out , the total volume velocity through the glottis U t , the motion of the vocal folds, and the length of the glottal chink l ch . The effect of the glottal chink on these physical quantities is clearly seen in both configurations (/a/ for the left column, and /i/ for the right column). For instance, as expected, the presence of the glottal chink stops U t from being null, since the glottis is never totally closed. The waveform of U t is then the superimposition of the glottal flow U g inside the vibrating part of the vocal folds and an offset corresponding to U ch . The U g component behaves similarly to the case of a totally closed chink, and U ch increases as l ch increases, i.e. when the glottal chink opens up. The effect is more visible in Fig. 11, which shows details of a few phonatory cycles at two different time instants, corresponding to a glottal configuration without chink, and with a chink of length 0.


N J Cira

A Benusiglio

& M Prakash

Vapour mediated sensing and motility in two-component droplets

 and chemotaxing cells 17 .

Text:

We observed that pure water and pure propylene glycol (PG) spread completely when placed on corona treated clean glass slides (Supplementary Information section 1). This is expected on such a high energy surface for which the spreading parameter, defined as , is larger than zero, where represents the surface energy of the solid/vapor, liquid/vapor, and solid/liquid interfaces [START_REF] De Gennes | Capillarity and wetting phenomena : drops, bubbles, pearls, waves[END_REF] . Surprisingly, mixtures of PG and water formed droplets with apparent contact angles, , even though . The trend in went from zero to a maximum value and back to zero as PG was added to water (Fig. 2a) which cannot be simply explained by the monotonically decreasing liquid/vapour surface tension (Extended Data Fig. 1) [START_REF] Hoke | Surface Tensions of Propylene-Glycol Plus Water[END_REF] . Breathing on a droplet noticeably modified the contact angle. To quantify this observation, we deposited droplets in controlled humidity chambers and found that apparent contact angle decreased with relative humidity, and droplets spread under saturated humidity (Fig. 2b), suggesting vapour played a role in droplet stabilization.

Using tracer beads (1 μm diameter) we visualized an internal flow from centre to edge along the bottom of the droplet, similar to the flow in the 'coffee ring' effect [START_REF] Deegan | Capillary flow as the cause of ring stains from dried liquid drops[END_REF] . We also observed a flow from the edge to the centre along the top of the droplet, at higher velocity than the outward flow (Fig. 2d and2e). This less commonly seen 'counter flow' has been observed with surfactant or thermal gradients only in pinned droplets [START_REF] Hu | Analysis of the effects of Marangoni stresses on the microflow in an evaporating sessile droplet[END_REF][START_REF] Truskett | Influence of surfactants on an evaporating drop: Fluorescence images and particle deposition patterns[END_REF] . It collects tracer beads at the liquid-vapour interface into a prominent ring (Fig. 2e). Microscopic observation of the droplets revealed a thin film extending tens of micrometres from the edge of the bulk droplet into which the 1 μm tracer beads did not enter (Fig. 2d, Supplementary Video 4). For the same droplets on a lower energy surface the counter flow was confined to the border of the droplet (Fig. 2d). No tracer bead ring appeared (Supplementary Video 4). There was no thin film around the droplets, and the droplets were less mobile and did not interact.

From these observations we can understand the mechanism which prevents complete spreading. The high energy surface favours spreading of the droplet, as seen for pure liquids [START_REF] Tanner | Spreading of Silicone Oil Drops on Horizontal Surfaces[END_REF] . For a twocomponent droplet, the more volatile compound (water here) evaporates more quickly than the less volatile compound (PG). Evaporation is faster at the border of the droplet than the bulk [START_REF] Deegan | Capillary flow as the cause of ring stains from dried liquid drops[END_REF] , and the border of the droplet has a higher surface area to volume ratio. Therefore PG, with a lower than water, is left in higher concentration at the border than the bulk. The resulting gradient of surface tension pulls liquid towards the centre along the top of the droplet, an effect shown to slow down or stop spreading [START_REF] Bascom | Dynamic surface phenomena in the spontaneous spreading of oils on solids[END_REF][START_REF] Bernett | Prevention of Liquid Spreading or Creeping. Contact Angle, Wettability and Adhesion[END_REF][START_REF] Pesach | Marangoni Effects in the Spreading of Liquid-Mixtures on a Solid[END_REF] . Here the spreading is stopped resulting in a droplet with a stable apparent contact angle (Extended Data Fig. 2) surrounded by a thin film(Fig. 2d,3a). [START_REF] Bascom | Dynamic surface phenomena in the spontaneous spreading of oils on solids[END_REF] Next we build a simple model to test this mechanism of droplet stabilization. We assume a sharp transition of surface tension between the bulk droplet (

) and the surrounding thin film (

). We introduce a quasi-static horizontal force balance at the intersection of the thin film and the bulk droplet, . To calculate , we model the water loss from the thin film due to evaporation, estimating the water fraction and surface tension of the film as a function of external relative humidity ( ) and water fraction of the droplet (Supplementary Information section 2.3, Extended Data Fig. 3 and4). Using this model we fit a single parameter for 40% and observe that the prediction globally captures the nonmonotonic contact angle curve and accounts for variation in this curve as a function of relative humidity (Fig. 2a). Our current model only accounts for water evaporation, and is therefore less accurate at high PG concentration and high humidity.

Based on this model, for any two miscible chemicals on a high-energy substrate, droplets should form if and only if one of the chemicals in the mixture has both a higher surface tension and higher vapour pressure (quadrant I and III, Fig. 2c). To test this law, we placed various twocomponent mixtures on corona treated slides. In about 200 unique combinations (Extended Data Table 1), droplet formation vs spreading was well predicted, excluding reactive pairs (Fig. 2c), and these droplets had similar attributes to the PG/water system, such as high mobility and interactions. We also deposited PG/water droplets on other high-energy substrates -piranha treated glass, clean silicon wafers, and freshly scraped steel -and found similar behaviour. These two-component droplets have characteristics of both wetting and non-wetting liquids: they maintain a defined contact angle, but sit on a thin fluid film. As long as , the droplets should not 'feel' the solid surface, and chemical inhomogeneities and roughness should not cause pinning. The droplet contact angle is also independent of the substrate ( ) and surface roughness. Without pinning, the droplets display high mobility and hence move under the influence of minute forces. We do not observe high mobility on low energy surfaces with a three phase contact line, where high hysteresis inhibits droplet motion.

When two droplets were deposited at distances of up to several radii apart, they moved toward each other from over a wide range of concentrations, even when both droplets had the same concentration (Fig. 1b,1c,3c). Droplets increased speed as they approached each other (Fig. 3b). These long-range interactions were preserved even across a break in the glass slide (Supplementary Video 5). PG/water droplets followed a pipette tip containing water placed near but not touching the droplet or the glass slide (Supplementary Video 6). These observations and our measurements of vs (Fig. 2b) led us to the surprising conclusion that long-range interactions were vapour mediated.

Based on the observations above, we propose a mechanism for vapour-mediated interactions different than in previous systems [START_REF] Bangham | The behaviour of liquid drops and adsorbed films at cleavage surfaces of mica[END_REF][START_REF] Carles | Spreading Involving the Marangoni Effect -Some Preliminary-Results[END_REF] . Evaporation from a sessile droplet is known to produce a vapour gradient [START_REF] Hu | Evaporation of a sessile droplet on a substrate[END_REF] . Since the vapour pressure of water is 100 times larger than the vapour pressure of PG, the dominant vapour is water. Two neighbouring droplets each lie in a gradient of water vapour produced by the other (Fig. 3a). This gradient causes a local increase in relative humidity and thus decreased evaporation of the thin film on the adjacent portions of the droplets, breaking symmetry. The decreased evaporation leads to an increased water fraction in the thin film, hence increasing locally. Asymmetric around the droplet causes a net force that drives the droplets towards each other.

To test this mechanism, we next propose a mathematical model to calculate the expected distance between two identical droplets as a function of time (Fig. 3a, Supplementary Information section 2.4). We start with the diffusion equation to estimate the relative humidity profile around a droplet. By utilizing our prior measurements of of a static droplet as a function of uniform external humidity, we estimate the local around each droplet as a , where is the slope of the vs plot (Fig. 2b), is the radius of the droplet, is the distance between the droplet centres, and is the parameter of integration. This net force causes droplet motion and is balanced by a viscous drag force ( ). Here we neglect inertia since the Reynolds number ( ) is smaller than 1 (for typical droplet velocity mm/s and droplet radius 1 mm, ).

We calibrated by measuring droplet speed on ramps of known angle, observing that it scaled linearly with the velocity ( ) as (Extended Data Fig. 6). The drag coefficient was a linear function of the droplet perimeter, consistent with existing theory based on viscous dissipation at three-phase contact lines 29 (Supplementary Information section 2.1 and 2.2, Extended Data Fig. 5 and6 ). Equating with , we obtain and integrate the instantaneous velocity to arrive at the distance between the two droplets, . Plotting as a function of with as the time of droplet contact, we observe a good agreement between model and data, with no adjustable parameters (Fig. 3b). In a log-log plot behaves as a scaling law of exponent at long distance, which is also captured by the model (Fig. 3b inset).

In Fig. 3c, we present a phase diagram of long-range interactions between one pinned droplet and one mobile droplet, as a function of concentration of both droplets. Over a large concentration range the mobile droplet was attracted to the pinned droplet. However, when [PG] pinned >>[PG] mobile , the mobile droplet fled, indicating a repulsive force. We hypothesize that at high PG concentration, the gradient of PG vapour begins to play a role, decreasing and driving the mobile droplet away.

At short range, two droplets of like concentrations coalesce upon contact. Droplets of sufficiently different concentrations can undergo a prolonged 'chasing phase' [START_REF] Bangham | The behaviour of liquid drops and adsorbed films at cleavage surfaces of mica[END_REF] as explained by Riegler and Lazar [START_REF] Riegler | Delayed coalescence Behavior of droplets with completely miscible liquids[END_REF] (Fig. 1b). Fluid is directly exchanged between the droplets, as visualized by a fluorescent dye (Supplementary Video 7). This exchange of fluid leads to a surface tension gradient across both the droplets, where the droplet of lower surface tension 'chases' the droplet of higher surface tension, which in turn 'flees' away [START_REF] Riegler | Delayed coalescence Behavior of droplets with completely miscible liquids[END_REF] . Additional subtleties of short-range interactions can be obtained by adjusting concentrations and volumes (Supplementary Information, Extended Data Fig. 7).

Using the fundamental understanding we developed for this system, we built several self-fueled surface tension driven fluidic machines out of the everyday materials such as food colouring, glass slides, and Sharpie™ marker (Supplementary Information section 1.5). First, we used the long-range interactions to create a droplet self-aligner, which aligns randomly placed droplets of identical concentrations in different 'lanes' into a single straight line (Fig. 4a, Supplementary Video 8). Second, we utilized the short-range interactions to create sustained droplet chasing, during which droplets circled around a single loop for several minutes (Fig. 4b, Supplementary Video 9). We note that since the droplets are unaffected by prior trajectories, the droplets are able to repeatedly cross over their own paths. Third, we created a completely vertical droplet oscillator by placing a large low surface tension droplet beneath a higher surface tension droplet bounded in a lane on a glass slide (Fig. 4c, Supplementary Video10). Finally we created a concentration based autonomous sorter that segregated small droplets into reservoirs based on their surface tension. In this device, we relied on gravity to bring droplets down a ramp, where they sampled wells from low to high surface tension, merging only when they reached a like concentration (Fig. 4d, Supplementary Video 11). These examples illustrate the wide variety of autonomous sensing and motility based devices that can be created using this system Online Content Methods, along with any additional Extended Data display items and Source Data, are available in the online version of the paper; references unique to these sections appear only in the online paper

Figures:Figure 1 |

 1 Figures:

Figure 2 |

 2 Figure 2 | Individual droplet characteristics. a, Isolated droplets (0.5 μL) on a clean glass surface display a non-monotonic apparent contact angle as a function of PG %. Crosses and triangles indicate data taken at 75% RH and 40% RH respectively. Dashed lines indicate the model's fit to the data. b, The cosine of the apparent contact angle varies linearly (line of best fit shown) with external humidity for 0.5 μL 10% PG droplets. Error bars are the range of three measurements at 75% RH. c, Behaviour of two-component mixtures of all nonreactive combinations of 21 miscible fluids (see supplementary Table S1 for chemical list) on clean corona treated glass. For each liquid pair, difference in surface tension is plotted against difference in vapour pressure. Red dots indicate droplet formation while black crosses indicate complete wetting. d, Important differences between two-component droplets deposited on high and low energy solid substrates. From top to bottom: accumulation of beads at the liquid/vapour

Figure 3 |

 3 Figure 3 | Long-range droplet interactions. a, Schematic of vapour gradients and evaporation from two droplets distance apart. Increased vapour concentration between the droplets leads to less evaporation. b, Distance between droplets as a function of time before contact for two freely moving 0.5 μL 10% PG droplets. The error bars represent the standard deviation of 12 experiments, and the dashed line is the model prediction, inset: log-log scale, solid line is the power law fit. is the time of droplet contact. c, Phase diagram of interactions between a single pinned and a single free 0.5 μL droplet (the axes are volume % of PG). Each dot represents an experiment; each colour indicates the direction of motion of the free droplet.

Figure 4 |

 4 Figure 4 | Droplet based devices. We created four devices by drawing permanent marker (Sharpie™, black) lines, which are hydrophobic enough that droplets do not cross them (Supplementary Videos 8-11). a, Spontaneous droplet aligner. Upper left inset: 0.5 μL green droplets of 10% PG are dispensed at random initial positions separated by 5 mm spaced sharpie lines. Upper right inset: droplets have automatically aligned into final positions. The graph shows the y position of each droplet as a function of time. The colour code represents the x position in the aligner. b, Short-range chasing between a 1% PG droplet (red) and a 25% PG droplet (blue) in a 2.1 cm mean diameter circle. Inset: a three image time lapse (10 s spacing, arrows representing direction of motion). The graph shows the travelled distance as a function of time. c, We deposited a 25% PG droplet (blue) above a 1% PG droplet (red) bounded in a 4 mm lane on a vertical glass slide. The top droplet oscillates up and down. Top: one oscillation, images are separated by 1 s. Bottom: vertical position of the top droplet as a function of time. d, Surface tension sorter. Schematic: wells of various concentrations of PG (colours) are confined by sharpie lines (black). Concentrations from top to bottom are 30%, 25%, 20%, 15%, 10%, 5% PG. Each image shows the time lapse trajectory of a droplet as it is deposited at the top and moves down due to gravity, sampling each well, only merging with a well of like concentration. Sorting happens by purely passive means.

Acknowledgments:

We thank all members of the Prakash Lab for discussions. We thank Justin Williams for early support of this work, Bertrand Buisson for discussions, and Graham Dick for discussions and reagents. N.J.C. is supported by an NSF GRFP fellowship. A.B. is supported by Pew Foundation. M.P. is supported by the Pew Program in Biomedical Sciences, Terman Fellowship and Gordon and Betty Moore Foundation. Author Contribution: N.J.C made the original observation. All authors designed the research. N.J.C and A.B. conducted experiments, and all authors interpreted the data; N.J.C. and A.B. developed the models. N.J.C and A.B. wrote the manuscript, and all authors commented on it.

Supplementary Information is available in the online version of the paper.

Author

Information Reprints and permissions information is available at www.nature.com/reprints. The authors declare no competing financial interests. Readers are welcome to comment on the online version of the paper. Correspondence and requests for materials should be addressed to N.J.C. (ncira@stanford.edu) or M.P. (manup@stanford.edu).