Benjamin Elie
email: benjamin.elie@loria.fr

Yves Laprie

Role of glottal abductions in fricatives (Version by the authors) Acoustic impact of the gradual glottal abduction degree on the production of fricatives: A numerical study

Keywords: numbers: 43.70.Bk Speech production, Phonetics, Fricative, Glottal gap

The paper presents a numerical study about the acoustic impact of the gradual glottal opening on the production of fricatives. Sustained fricatives are simulated by using classic lumped circuit element methods to compute the propagation of the acoustic wave along the vocal tract. A recent glottis model is connected to the wave solver to simulate a partial abduction of the vocal folds during their self-oscillating cycles. Area functions of fricatives at the three places of articulation of French have been extracted from static MRI acquisitions. Simulations highlight the existence of three distinct regimes, named A, B, and C, depending on the degree of abduction of the glottis. They are characterized by the frication noise level: A exhibits a low frication noise level, B, which is a transitional unstable regime, is a mixed noise/voice signal, and C contains only frication noise. They have significant impacts on the first spectral moments. Simulations show that their boundaries depend on articulatory and glottal configurations. The transition regime B is shown to be unstable: it requires very specific configurations in comparison with other regimes, and acoustic features are very sensitive to small perturbations of the glottal configuration abduction in this regime.

I. INTRODUCTION

Fricatives are a class of consonants that are produced by creating a supraglottal constriction in the vocal tract so that a turbulent airflow is generated, usually downstream of the constriction. This results in the production of the so-called frication noise, which is the main characteristic of the fricative consonants. Voiceless fricatives are produced by adjusting the glottal opening area such that it is significantly greater than the area of the supraglottal constriction [START_REF] Stevens | Acoustic Phonetics[END_REF] . In that case, the glottis is completely abducted, and the generated sound contains only frication noise. In voiced fricatives, the glottis is adjusted so that both the frication noise and the voiced source contributions are mixed in the produced speech. The aeroacoustic conditions required to produce voiced fricatives are then very specific.

Studies about fricatives have focused on the spectral characteristics of the produced sound [START_REF] Jongman | Acoustic characteristics of english fricatives[END_REF] , the frication noise source 3 , the geometry of the supraglottal constriction [START_REF] Ramsay | The influence of geometry on the initiation of turbulence in the vocal tract during the production of fricatives[END_REF] , and the vocal tract configuration downstream of the supraglottal constriction [START_REF] Shadle | Revisiting the role of the sublingual cavity in the /s/-/S/ distinction[END_REF] , especially the potential obstacles (teeth, lips, etc) encountered by the airflow [START_REF] Narayanan | An articulatory study of fricative consonants using magnetic resonance imaging[END_REF] . These studies have contributed to a better understanding of the specific aeroacoustic conditions that are required to produce fricatives. It has been shown that the frication noise is generated when the airflow becomes turbulent, namely when the Reynolds number is sufficiently high [START_REF] Stevens | Acoustic Phonetics[END_REF] . High Reynolds numbers occur when the cross-sectional area of the supraglottal constriction is small and/or when the low-frequency component of the acoustic volume velocity through the constriction is large.

The latter condition implies that the glottal opening area should be sufficiently large. For voiced fricatives, in addition to these conditions for the generation of frication noise, other conditions at the vicinity of the glottis are required to produce the voicing: the vocal folds should not be completely abducted, and the transglottal pressure drop must be large enough to guarantee self-oscillations [START_REF] Ishizaka | Synthesis of voiced sounds from a two-mass model of the vocal cords[END_REF] . This implies subtle adjustments of the geometry of both the supraglottal constriction and the configuration at the glottis.

Most studies about the perceptual distinction between voiced and voiceless fricatives at the same place of articulation have focused on the duration of the voiced and the voiceless parts of the considered fricative [START_REF] Stevens | Acoustic and perceptual characteristics of voicing in fricatives and fricative clusters[END_REF][START_REF] Jesus | A parametric study of the spectral characteristics of european portuguese fricatives[END_REF][START_REF] Dmitrieva | Phonological neutralization by native and non-native speakers: The case of russian final devoicing[END_REF][START_REF] Pape | Intervocalic fricative perception in European Portuguese: An articulatory synthesis study[END_REF] . The consonant voicing is then considered as binary.

However, the amount of energy of the voicing component over that of the frication noise component is likely to vary continuously during the fricative segment. Indeed, one may assume that the motions of the articulators and the abduction movements of the vocal folds, which are relatively slow in comparison with the oscillations of the vocal folds, gradually modify both the amplitude of the frication noise source and the amplitude of the voiced source. The acoustic impact of the vocal tract configuration on frication noise sources has been widely studied [START_REF] Shadle | Revisiting the role of the sublingual cavity in the /s/-/S/ distinction[END_REF][START_REF] Narayanan | An articulatory study of fricative consonants using magnetic resonance imaging[END_REF] , but little attention has been paid to the acoustic impact of the configuration at the glottis. In a previous study [START_REF] Löfqvist | Vocal tract aerodynamics in /aCa/ utterances: Measurements[END_REF] , it has been experimentally shown that glottal pulses characteristics are subjected to large variations at vowel-consonant transitions. More interestingly, the presented results suggest that the glottal abduction movement starts before the vowel-consonant transition, and that the glottal adduction movement ends after the consonant-vowel transition. This implies a gradual abduction of the glottis, along with sustained vocal folds oscillations, by creating a posterior membranous gap, sometimes referred to the linked leak [START_REF] Cranen | Modeling a leaky glottis[END_REF] . Recently, the existence of such a glottal gap due to the partial abduction of the vocal folds has been proposed to be an important feature in the production of voiced fricatives [START_REF] Elie | A glottal chink model for the synthesis of voiced fricatives[END_REF] . Since the membranous opening directly acts on the acoustic volume velocity, it also modifies the amplitude of the frication noise source. Thus, bad coordination between the glottal opening and the geometry of the constriction may result in an uncontrolled frication noise, and is likely to produce voiced fricatives that are perceived as voiceless because of a too large amount of frication noise.

The paper develops this idea: it uses the partial abduction model introduced in our recent papers [START_REF] Elie | A glottal chink model for the synthesis of voiced fricatives[END_REF][START_REF] Elie | 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[END_REF] for a numerical study about the acoustic impact of the configuration at the glottis on the production of the voiced fricatives as a function of the geometry of the supraglottal constriction and the place of articulation. It focuses on the influence of the membranous glottal opening, which is related to the degree of abduction of the glottis, on some acoustic features of the speech signal, such as the spectral centroid, the spectral spread, and the voicing quotient. The aim is to define the boundaries of the different regimes of production of fricatives in a phonatory-articulatory space spanned by the glottal abduction degree, the position, and the geometry of the supraglottal constriction. After presenting the acoustic model in Sec. II, the paper details the configurations of the vocal tract used for the numerical simulations in Sec. III as well as the acoustic features that are investigated. Results of the simulations are presented in Sec. IV, which discusses the impact of the glottal abduction degree on the first spectral moments of the simulated fricatives, and on the amount of the generated frication noise. Finally, the boundaries of the different regimes of fricative production and their impact on the phonetic strategies to contrast voiced and voiceless fricatives are discussed and compared with data from a real speaker in Sec. V.

II. ACOUSTIC MODEL

The simulation framework used to compute the acoustic propagation inside the vocal tract is derived from the transmission line circuit analog (TLCA) approach [START_REF] Maeda | A digital simulation method of the vocal-tract system[END_REF] . It considers plane waves propagating along a spatially sampled vocal tract, modeled as a set of connected acoustic tubes, or tubelets. Unlike the other widely used approach, the reflection type line analog (RTLA) model [START_REF] Kelly | Speech synthesis[END_REF][START_REF] Story | Phrase-level speech simulation with an airway modulation model of speech production[END_REF] , it easily deals with time-varying lengths of the vocal tract, and also with uneven spatial sampling of the vocal tract. For further information, the reader may find a detailed review of existing techniques for speech synthesis in Ref. [START_REF] Kröger | Articulatory synthesis of speech and singing: State of the art and suggestions for future research[END_REF] . The framework that is used in this paper considers recent improvements of TLCA-based techniques [START_REF] Elie | 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[END_REF][START_REF] Mokhtari | Single-matrix formulation of a time domain acoustic model of the vocal tract with side branches[END_REF] , such as the possibility of connecting self-oscillating models of the vocal folds with incomplete closure due to the presence of a membranous glottal gap.

A. Acoustic propagation

Considering a general case, the vocal tract is seen as a waveguide network, where each waveguide represents a side cavity (the oral tract, the nasal tract, the piriform fossae, etc).

TLCA-based techniques have shown [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] that the wave propagation inside such networks is driven by a set of linear equations. In a matrix form, it writes

f = Zu, (1)
where

f ∈ R (N +1) is a vector containing pressure forces, Z ∈ R (N +1)×(N +1
) is a tridiagonal matrix containing impedance and loss terms associated to each tubelet, and u ∈ R N +1 is the vector containing the volume velocities inside each tubelet.

When dealing with self-oscillating models of the vocal folds, a quadratic term accounting for the pressure drop inside the glottal constriction, due to the Bernoulli resistance, should be added to the first line of the system 15 :

f = Zu Z + Qu Q , ( 2
)
where Q is a square matrix the same size as Z having only one non-zero element, that is

Q (1,1) = R b , and u Q ∈ R (N +1) = [U 2 1 , U 2 2 , . . . , U 2 N +1
] T is the vector containing the square power of the volume velocities. The term R b is the Bernoulli resistance [START_REF] Elie | 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[END_REF] . Eq. 2 is computed at each time step, and the volume velocities u are the solution of the equation. The acoustic pressure P Out is computed as the first time derivative of the volume velocities at the lip termination, i.e. U N +1 , where N is the number of tubelets that models the vocal tract.

The domain of validity of the plane waves assumption depends on the radius of the maximal cross-sectional area of the vocal tract. Indeed, the cut-off frequency under which the assumption is considered as valid [START_REF] Eriksson | Higher order mode effects in circular ducts and expansion chambers[END_REF] is

f c = 0.5861c s 2r max , (3)
where c s is the speed of sound and r max is the radius of the largest vocal tract section. In this study, f c differs according to the considered place of articulation: it is around 10 kHz, 8 kHz, and 7 kHz for the palato-alveolar, the alveolar, and the labiodental fricatives, respectively.

B. Frication noise generation model

The frication noise is generated by a turbulent air flow that appears downstream of the supraglottal constriction 1 . Interactions with obstacles, such as the teeth, the lips, and the walls of the vocal tract, usually occur and impact the acoustic nature of the turbulent noise source. The aeroacoustic mechanisms that are involved in the production of the frication noise are not totally comprehended yet. Consequently, their complicated nature makes them challenging to integrate into simplified acoustic models. Indeed, small variations of geometric or biomechanic parameters, such as the jet angle, the nature of the obstacle, or wall discontinuities downstream of the supraglottal constriction, may significantly modify the spectral characteristics of the produced turbulent noise [START_REF] Shadle | Articulatory-acoustic relationships in fricative consonants[END_REF][START_REF] Narayanan | An articulatory study of fricative consonants using magnetic resonance imaging[END_REF] . Simplified models for frication noise sources commonly consider them as acoustic poles, dipoles, or quadripoles, depending on the cause of the turbulence [START_REF] Shadle | Articulatory-acoustic relationships in fricative consonants[END_REF]22 .

In this paper, the frication noise source is modeled as an acoustic dipole that is activated when the Reynolds number is above a certain critical value, arbitrarily chosen at Re c = 1700, referring to the previous study by Sondhi and Schroeter [START_REF] Sondhi | A hybrid time-frequency domain articulatory speech synthesizer[END_REF] . The noise source is a bandpass filtered Gaussian white noise [START_REF] Maeda | Phoneme as concatenable units: VCV synthesis using a vocal tract synthesizer[END_REF] . According, to Stevens 1 (p.388, 2nd edition), the amplitude of the noise source is proportional to U 3 DC /a 2.5 c , where a c is the cross-section area of the supraglottal constriction. Consequently, considering the i th section of the spatially sampled vocal tract, the amplitude of the noise source P n i (t) at section i and instant t is

P n i (t) = max 0, ξw c (t) Re 2 (t) -Re 2 c U 3 DC (t) a 5/2 i-1 (t) , ( 4
)
where ξ is an arbitrarily adjustable real constant used to control the noise level, and w c (t) is a colored noise function, U DC is the air flow volume velocity inside the vocal tract, and

a i-
+ 8π µd c a 2 c U DC -P Sub = 0, (5)
where κ is a scaling factor, which is 1.42 for a rectangular duct 7 , a, l, and d denote the area, the length, and the width of the glottis (index gl) and the supraglottal constriction (index c), and P sub is the subglottal pressure.

The frication noise source function w c (t) is computed from a Gaussian white noise function w g (t), shaped by the impulse response g(t) of a finite impulse response filter in order to simulate a colored noise, hence w c (t) = w g (t) * g(t).

Following previous studies 22, [START_REF] Birkholz | Enhanced area functions for noise modeling in the vocal tract[END_REF] , g(t) is the impulse response of a low-pass filter having a cut-off frequency of f c = 0.15

√ π U DC a 3/2 c
. The gain ξ has been empirically chosen. It has been tuned so that, when vowel-fricative-vowel (VFV) pseudowords are simulated, the relative level of energy of the fricative in comparison with the surrounding vowels is similar to that observed in natural speech signals. As expected, the amplitude ξ varies according to the place of articulation, and is higher for sibilant fricatives than for non-sibilant fricatives (ξ = 4.3×10 -7 for the palato-alveolar fricatives /Z,S/, ξ = 0.78×10 -7 for the alveolar fricatives /z,s/, and ξ = 0.61×10 -7 for the labiodental fricatives /v,f/).

C. Glottis model with a membranous glottal gap

The model used to simulate the glottis is similar to the one introduced in recent extensions of the Single-Matrix Formulation paradigm [START_REF] Elie | 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[END_REF] and used to simulate voiced fricatives [START_REF] Elie | A glottal chink model for the synthesis of voiced fricatives[END_REF] . It considers two distinct portions of the glottis, represented in Fig. 1: an oscillating part that is computed using a classic two-mass model of the vocal folds [START_REF] Lous | A twomass vocal-fold coupled to vocal tract and trachea, with application to prothesis design[END_REF] , and a partially abducted part, the so-called membranous glottal gap, that allows an incomplete closure along the length of the vocal folds during the oscillation cycles. The two-mass model of the vocal folds considers recent improvements to take into account smooth contours, a mobile separation point [START_REF] Lous | A twomass vocal-fold coupled to vocal tract and trachea, with application to prothesis design[END_REF] , the viscous losses and the unsteady flow effects 27,[START_REF] Bailly | Influence of a constriction in the near field of the vocal folds: Physical modeling and experimental validation[END_REF] . Note that other models of vocal folds with consideration of the membranous glottal gap have been proposed in the past [START_REF] Kröger | Simulation of vocal fold oscillation behaviour by a self-oscillating glottis model[END_REF][START_REF] Birkholz | Articulatory synthesis of words voice qualities using a modified two-mass model of the folds[END_REF][START_REF] Zañartu | Modeling the effects of a glottal opening vocal fold dynamics with implications for vocal hyperfunction[END_REF] .

We chose our model, which has been specifically designed to be used with the Single-Matrix Formulation paradigm. It should be noted that lumped models, such as the one used in this study, are simplifications of the realistic vocal folds dynamics, and that finite element models might provide a more realistic description of the acoustic impact of the abduction/adduction movement. Yet, lumped models are still useful for speech research as they allow qualitative investigations of the global behavior of the vocal folds to be made using a few parameters and with a very low computation cost, in comparison with finite element models. A review of existing self-oscillating models for the vocal folds may be found in Ref. [START_REF] Erath | A review of lumped-element models of voiced speech[END_REF] . Similarly to Birkholz et al. [START_REF] Birkholz | Articulatory synthesis of words voice qualities using a modified two-mass model of the folds[END_REF] , our model assumes the angle made by the vocal folds at the membranous gap tip to be constant, hence the length l ch is controlled by the abduction h ab :

l ch = l t h ab h max , (6)
where h max is the maximal abduction length, i.e. the length h ab for the full glottal abduction.

The glottal opening area is then A g = h ab l ch . In order to present generic results, for the rest of the paper, the abduction is controlled by a parameter called degree of abduction, denoted by D ab and expressed in percent, which represents the ratio between the area of the glottal gap and that of the fully abducted glottis. That is, D ab = 0% when there is no gap (full adduction), and D ab = 100% when the glottis is fully abducted. Note that this ratio is the square power of the ratio h ab /h max . We chose to use an area ratio rather than a length ratio because of the possibility to compare with the experimental measurements of the glottal opening area presented in Sec. V C.

The glottal gap has been shown to have significant acoustic impact on speech production, whether on the spectral tilt [START_REF] Cranen | Modeling a leaky glottis[END_REF] , or on self-oscillating movements of the vocal folds [START_REF] Zañartu | Modeling the effects of a glottal opening vocal fold dynamics with implications for vocal hyperfunction[END_REF] . Except for a few studies [START_REF] Löfqvist | Vocal tract aerodynamics in /aCa/ utterances: Measurements[END_REF][START_REF] Elie | A glottal chink model for the synthesis of voiced fricatives[END_REF] , the acoustic impact of the glottis configuration on the fricative production has been given little interest. However, its role during fricatives is important, since the anticipating abduction movement during the vowel preceeding the consonant modifies the voice quality, which becomes breathy at the vowel offset [START_REF] Löfqvist | Vocal tract aerodynamics in /aCa/ utterances: Measurements[END_REF] . It is thus likely that the breathy voice at the consonant onset favors the frication noise appearance during the fricative. Besides, external lighting and sensing photo-glottography (ePGG) measurements [START_REF] Honda | Glottal-opening and airflow pattern during production of voiceless fricatives: a new non-invasive instrumentation[END_REF] (cf. Fig. 2)

show a glottal opening waveform which is the superimposition of two components, a low frequency component that reaches the maximum around the middle of the fricative segment, and a higher frequency component that corresponds to the oscillations of the vocal folds.

The ePGG signal displayed in Fig. 2 is an example of time evolution of the glottal opening chosen among a set of VFV pseudowords recorded for several speakers. The coexistence of these two components confirms the idea of a partial abduction of the glottis, which results in a glottal configuration similar to the glottis configuration shown in Fig. 1.

We have shown that the glottal gap due to the incomplete closure of the vocal folds may be connected to the Single-Matrix Formulation as an acoustic waveguide connected in parallel to the oscillating part of the glottis [START_REF] Elie | 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[END_REF] . One may refer to previous 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 details about the numerical computation, and to our previous paper [START_REF] Elie | 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[END_REF] for details about the integration of the glottal leakage side branch.

The ability of our model to reproduce the actual time evolution of the glottal opening area has been tested by running several simulations. For each simulation, the time varying parameter D ab is set as the low frequency component of the ePGG data. The presented Then, the contours of the vocal tract have been extracted by hand on the midsagittal slice to compute the midline. This line, which should be perpendicular to the propagation of a plane wave in the vocal tract, is used to decompose it into tubelets. The method used to compute the midline is based on dynamic programming to select the best path of segments connecting the larynx to the lips [START_REF] Laprie | Articulatory copy synthesis from cine X-ray films[END_REF] . The determination of the midline is applied either on the whole vocal tract when there is no occlusion, which is the case with the vocal tract shapes of fricatives exploited in the paper, or on all the open sections delimited by the vocal tract extremities or occlusions.

The last step consists in dividing the vocal tract into tubelets perpendicularly to the midline. Attention is paid to the fact that two consecutive tubelets cannot cross in high curvature regions of the vocal tract. Recovering the third dimension from the 2D information, namely the midsagittal distance and the length of each tubelet, in order to estimate the area function has given rise to a number of works [START_REF] Mcgowan | Analyses of vocal tract crossdistance to area mapping: An investigation of a set of vowel images[END_REF][START_REF] Soquet | Mid-sagittal cut to area function transformations: Direct measurements of mid-sagittal distance and area with MRI[END_REF][START_REF] Perrier | Vocal tract area function estimation from midsagittal dimensions with two sets of coefficients[END_REF] . However, the improvement with Role of glottal abductions in fricatives (Version by the authors) respect to methods derived from that proposed by Heinz and Stevens [START_REF] Heinz | On the relations between lateral cineradiographs, area functions, and acoustic spectra of speech[END_REF] is not very marked.

We thus accepted the α β parameters proposed by Soquet et al. [START_REF] Soquet | Mid-sagittal cut to area function transformations: Direct measurements of mid-sagittal distance and area with MRI[END_REF] . The choice of the α -β transformation over 3D segmentation of the vocal tract is motivated by its generic nature that fits our study. Our main interest is to determine a global shape of the vocal tract, with realistic position, length, and area of the constriction. Fig. 3 shows the extracted area functions. For the sake of clarity, only the cardinal vowel contexts, /i,a,u/, are displayed for each place of articulation.

Trachea

The incomplete closure of the glottis during the oscillation cycle of the vocal folds leads to a constant coupling between the vocal tract and the trachea. It is therefore important to account for this coupling by connecting a subglottal waveguide upstream of the glottis.

Such connection with the single-matrix formulation has been previously proposed by Ho et al. [START_REF] Ho | An anatomically based, time-domain acoustic model of the subglottal system for speech production[END_REF] , but it uses a very complex geometry of the subglottal system, modeled as a tree-like structure that accounts for all the branching patterns of the bronchial airways. Achieving such a fine degree of modeling is not in the scope of the paper, since it is hard to get such information for real speakers. Consequently, the subglottal system is modeled here as a single waveguide connected to a pressure source located at the input. The area function that is used is borrowed from a previous study by Story [START_REF] Story | Phrase-level speech simulation with an airway modulation model of speech production[END_REF] .

Glottal parameters

The input parameters used for the glottis model are the same than the ones used in our previous paper [START_REF] Elie | 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[END_REF] . Values of the mass and stiffness of the vocal folds model have been chosen so that the fundamental frequency of the simulated voiced signal is approximately 150 Hz.

This corresponds roughly to the resonance frequency of the lumped mass-spring system.

The maximal abduction length h max = 1 mm, and the total length of the vocal folds l t is 2.2 cm, so that the maximal glottal opening A max = l t h max = 0.22 cm 2 . The length of the vocal folds is set to a rather high value, in comparison with the typical value for males ( 1.6 cm), in order to fit with our subject, who has a low pitched voice and a long vocal tract ( 18 cm, see Fig. 3).

B. Simulated configurations

The numerical study investigates the acoustic impact of various phonatory and articula-

Subglottal pressure

In this paper, the subglottal pressure refers to the constant value assigned to the pressure source connected to the trachea input. In the simulations presented in Sec. IV A, it varies from 500 to 1500 Pa, with an increment step of 100 Pa. These values correspond to subglottal pressures encountered in normal [START_REF] Giovanni | Estimated subglottic pressure in normal and dysphonic subjects[END_REF] and loud or singing [START_REF] Sundberg | How constant is subglottal pressure in singing[END_REF] voices. For the simulations presented in Sec. IV B, it is set to a nominal value P Sub = 1000 Pa.

Supraglottal constriction

The acoustic impact of the geometry of the supraglottal constriction is studied by modifying its cross-sectional area, denoted by a c , in the area functions extracted from MR images.

For the simulations presented in Sec. IV B, a c is set to different values ranging from 0.1 cm 2 to 0.5 cm 2 . For those presented in Sec. IV A, a c is taken as the minimal cross-section area extracted from MR images.

C. Investigated features

Regimes of fricatives: voicing quotient and minimal abduction of the vocal folds

In a previous study [START_REF] Elie | A glottal chink model for the synthesis of voiced fricatives[END_REF] , it has been shown that, for a given articulatory condition corresponding to a fricative, the simulated speech signal could exhibit three regions with distinct acoustic properties according to the glottal opening: i) an almost purely voiced signal when the glottis is almost entirely adducted (little frication noise is generated, similar to an approximant consonant), ii) a mixed voiced/noisy signal when the level of both voiced component and the frication noise share a similar order of magnitude, and iii) a purely noisy signal, similar to the voiceless fricative, when the voiced component is negligible in comparison with the frication noise. These three regions, respectively denoted A, B, and C in the rest of the paper, are studied as a function of the glottal abduction degree D ab . This is done by using a voicing index, named the voicing quotient [START_REF] Elie | Robust tonal and noise separation in presence of colored noise, and application to voiced fricatives[END_REF] . The voicing quotient (denoted V Q in this paper) is defined as the proportion, expressed in percent, of the energy of the periodic component in the speech signal, hence

V Q(%) = 100 × ||s p || 2 2 ||s p + s n || 2 2 , (7)
where s p is the periodic (or voiced) component of the signal, and s n is the frication noise signal. Both periodic and aperiodic components of the simulated signals are computed thanks to a specifically designed periodic/aperiodic decomposition technique of the speech signal, that has been proven to be robust in the case of colored noise such as in voiced fricatives [START_REF] Elie | Robust tonal and noise separation in presence of colored noise, and application to voiced fricatives[END_REF] . A value of 100 % indicates a speech signal containing only periodic components, and a value of 0 % indicates a speech signal containing only noisy components. Note that the voicing quotient is directly related to the Harmonics-to-Noise Ratio (HNR) that quantifies the harmonicity of the signal, but the definition of the voicing quotient is more adapted to our study, since it directly quantifies the amount of voicing in fricatives.

For a given area function, varying the glottal abduction degree leads to a decreasing voicing quotient as the glottis opens up. A typical curve, shown in Fig. 4, highlights the presence of these regions. In A, the voicing quotient is constantly high, around 90%, since just a little frication noise is generated so that the voiced components dominate in the simulated speech signal. In B, the frication noise is generated along with a voiced source, the latter being produced by the oscillating part of the vocal folds. The amount of frication noise that is generated is very sensitive to small perturbations of the abduction degree, hence large variations of the voicing quotient. Then, in C, the voicing quotient vanishes. In this regime, the noise source is predominant over the voiced source, leading to a voiceless fricative. This can be seen in the spectrum in Fig. 4, where the harmonics at f = nF

Spectral characteristics

The studied spectral characteristics investigated in this paper are the first two spectral moments, namely the spectral centroid, and the spectral spread.

The spectral centroid is a measure of the balance between the low and high frequency contributions in a spectral distribution, it writes

S 1 = N n=1 f n A n N n=1 A n , (8)
where f n , with n = 1, 2, . . . , N , indicates the positive frequency bins of the 2N -order FFT, and A n are the corresponding magnitude.

The spectral spread is a measure of the square root of the variance of the spectral distribution. A small spectral spread indicates a spectrum that concentrates its energy in a small

N n=1 A n (f n -S 1 ) 2 N n=1 A n . (9)
Both S 1 and S 2 are expressed in Hz and are computed in the frequency domain ranging from 50 Hz to 10 kHz.

IV. ACOUSTIC FEATURES

A. Effect of the subglottal pressure 1. Voicing quotient From left to right: palato-alveolar fricatives /S,Z/, alveolar fricatives /s,z/, and labiodental fricatives /f,v/. Each row of figures corresponds to a vowel context, /i,a,u/, from top to bottom.

Fig. 5 shows the voicing quotient as a function of the glottal abduction degree and the subglottal pressure. For the sake of clarity, it does not show results of all of the 21 configurations (7 contexts for each of the 3 places of articulation), but only those corresponding to the context of cardinal vowels, namely /i,a,u/. The same simplification applies in the rest of the paper. The general behavior of the voicing quotient is similar to the typical curve displayed in Fig. 4: it is constantly high for weak abductions, then it suddenly plunges at a certain point and vanishes at a second critical point. For palato-alveolar fricatives, these critical abduction degrees tend to be smaller when the subglottal pressure rises. Indeed, when the subglottal pressure increases, this raises the Reynolds number due to an higher DC component of the airflow volume velocity inside the vocal tract. Consequently, the Reynolds critical number Re c above which the frication noise is generated is reached at smaller glottal openings. In a D ab -P Sub plane, as shown in Fig. 5, it results in the left part exhibiting high values of voicing quotient, and the right part exhibiting low values, both parts being separated by a small transition area. Note that there is no significant variations with the subglottal pressure for the alveolar and labiodental fricatives.

Also, for palato-alveolar fricatives, the glottal abduction degree for which the voicing quotient starts to significantly decrease is smaller than for other places of articulation. The lowering of this critical abduction seems to be more important for palato-alveolar fricatives.

The vowel context seems to have little influence on the voicing quotient. Indeed, for each place of articulation, the shapes of the voicing quotient values for the three vowel contexts are very similar.

Spectral centroid

The spectral centroid as a function of phonatory conditions shows distinct areas, as evidenced in Fig. 6. At the bottom left, where the abduction degree and the subglottal pressure are small, the spectral centroid is constantly between 1000 and 2000 Hz. At the top right, when both the abduction degree and the subglottal pressure are large, the spectral centroid is much higher, and lies around 4500 Hz for palato-alveolar fricatives, and up to 6000 Hz for alveolar fricatives. these areas, the spectral centroid is relatively stable in regards with variations of the phonatory configurations. The values in this region depend on the place of articulation considered: alveolar fricatives have the highest values of spectral centroid, around 6 kHz, then the labiodental fricatives have a spectral centroid around 5.6 kHz, and finally, the palato-alveolar fricatives have a spectral centroid around 4.5 kHz. These values are in agreement with data observed in natural speech by Jongman et al. [START_REF] Jongman | Acoustic characteristics of english fricatives[END_REF] for English, and Lonchamp for French [START_REF] Lonchamp | Description acoustique (acoustic description)[END_REF] . The alveolar fricatives of English and French speakers exhibit the higher spectral centroid than other fricatives, whereas palato-alveolar fricatives have the lowest spectral centroid. Same data also showed that, like in the presented simulations, voiced fricatives exhibit lower spectral centroids than their voiceless counterpart. A small transition area may be seen. In this range, small perturbations of the glottal opening or of the subglottal pressure lead to large variations of the spectral centroid.

Similarly to the voicing quotient, the vowel context seems to have no significant influence on the spectral centroid of the simulated fricative. However alveolar and palato-alveolar fricatives simulated in the /u/ context exhibit slightly different patterns of spectral centroid values than for /i/ and /a/ context. For palato-alveolar fricatives, the high values (above 4 kHz) are reached for weaker abductions and smaller subglottal pressures than for other places of articulation. Alveolar fricatives in the /u/ context show smaller spectral centroids (around 4.8 kHz, while it is around 6 kHz for /i,a/ contexts). This may be due to the strong lip protusion and the small lip opening of the speaker in that case, as shown by the area functions in Fig. 3, in comparison with the other contexts. The main effect is then to significantly lower the formant frequencies, and consequently the spectral centroid.

Variation patterns of the spectral centroid as a function of the glottal abduction may be explained by the fact that when the subglottal pressure and the glottal opening are small, the aerodynamics conditions are not fulfilled to generate a frication noise, so that the simulated speech signal contains only periodic, or voiced components. In that case, the energy is mainly concentrated in the low frequency range, hence a low spectral centroid.

Then, when the glottal opening and the subglottal pressure are sufficiently large to generate a frication noise (central region), noisy components arise in the simulated speech signals.

This has for effect to enhance the high frequency range of the spectrum. Hence the rise of the spectral centroid. A slight increase of either the subglottal pressure or the glottal opening in this regime leads to a decrease of the voiced component level, and an increase of the noise level, hence the rise of the spectral centroid. When the voiced components have disappeared, at high values of D ab and P Sub , and when the simulated signals contain only noise, the spectral centroid reaches its maximal value, since the low frequency domain of the spectrum is significantly weakened by the absence of the voiced contributions. This is confirmed by the similarities between the pattern of voicing quotient values in Fig. 5 and that of spectral centroid values in Fig. 6.

Spectral spread

The plot of the spectral spread as a function of P Sub and D ab , shown in Fig. 7, also highlights the presence of areas that are related to both the voicing quotient and the spectral centroid. Again, the palato-alveolar fricative exhibits a different behavior than other places of articulation. For alveolar and labiodental fricatives, the spectral spread values at the bottom left corner is low, and then suddenly increases to reach a maximum around 3 kHz.

At the top right corner, the spectral spread is also small, with values at the same order of magnitude as at the bottom left corner. They are also relatively stable. These three areas are approximately located at the same place as the three areas of the spectral centroid. Palatoalveolar fricatives present smaller values of the spectral spread, and they are relatively constant independently of the position in the D ab -P Sub plane. Similarly to the spectral centroid, the existence of three different regions, as well as their characteristics may be explained by the degree of contribution of the frication noise. In weak abduction and low subglottal pressure conditions, the absence of frication noise yields to an almost purely harmonic spectrum with a marked spectral slope. Hence a small spectral spread. As long as no frication noise is generated, the spectral shape remains the same independently of the values of the glottal opening and the subglottal pressure. As soon as the frication noise is generated, the spectrum of the simulated fricative is enhanced in the high frequency domain. At the same time, the voiced component is still predominant in the low frequency range. As a consequence, the spectral spread is large, since the signal contains energy both in the low frequency domain (voiced component) and in the high frequency domain (noisy component). The spectral spread value is very sensitive to the balance between the voiced component energy level and the noise level, hence quick changes in the spectral spread values in the central area. Finally, when values of both P Sub and D ab are large enough so that there is no voiced in the simulated fricative any longer, the spectrum in the low frequency domain is significantly weakened so that the energy is concentrated in the mid-and high-frequency domains, hence smaller values of the spectral spread.

In the transition regime, where the spectral spread is unstable, the alveolar fricative in the /u/ context (around 2150 Hz) shows smaller spectral spreads than in other contexts (around 2600 Hz). However, in the stable region corresponding to the maximal glottal abduction, the spectral spread is then larger in the /u/ context (around 1800 Hz) than in other contexts (around 1300 HZ). There is no significant differences among vowel contexts for the labiodental fricative. The differences observed for the /u/ context are certainly due to the fact that the lip protrusion is stronger, and that the lips are closer each other than in other contexts. However, labiodental fricatives are characterized by small lip opening whatever the vowel context, hence little acoustic differences among contexts.

B. Influence of the cross-sectional area of the supraglottal constriction

Voicing quotient

The voicing quotient as a function of the glottal abduction degree D ab and the supraglottal constriction area a c is shown in Fig. 8. As expected, the voicing quotient is high for weak abductions, then dramatically decreases at a critical abduction degree, and vanishes for large values of D ab . The influence of the supraglottal constriction area is limited when a c is larger than 0.3 cm 2 . Under this value, a c has a more marked influence: the value of D ab above which the voicing quotient vanishes exhibits a local minimum for a c around 0.1 cm 2 .

Spectral centroid

Plots of the spectral centroid of the simulated signals as a function of D ab and a c , shown in Fig. 9, exhibit complicated patterns. A strong glottal abduction leads to higher values of the spectral centroid. The supraglottal constriction area has also a significant influence on the spectral centroid. Basically, a small a c leads to a high value of the spectral centroid. It results in low values of the spectral centroid at the top left corner, and high values at the bottom right corner of the D ab -a c plane. There is no significant differences according to the vowel context, except for /su/, which exhibits smaller spectral centroids, as already observed in Fig. 6. It is also worth noting than for a given a c , the spectral centroid suddenly rises for D ab values that roughly correspond to quick variations of the voicing quotient in Fig. 8. The local minimum at a c = 0.1 cm 2 is then clearly visible for alveolar and labiodental fricatives.

It is less marked for palato-alveolar fricatives, as they exhibit more complicated patterns. This suggests that the cross-sectional area of the supraglottal constriction in postalveolar fricatives seems to have a more significant impact on the spectral properties.

Spectral spread

Fig. 10 shows the spectral spread as a function of D ab and the constriction area a c . For a given constriction area, one can observe the same behavior as that in Fig. 7: the spectral spread is the smallest at both extremities of the glottal abduction degree, and reaches a maximum in a central area. The spectral spread is also larger for narrow constrictions. Similarly to what has been observed in Fig. 7, palato-alveolar fricatives have smaller spectral spreads than other places of articulation. This is certainly due to the fact that postalveolar fricatives present a predominant peak 2 , located between 2 and 3 kHz, which concentrates the energy of the spectrum in this frequency domain. Like for the spectral centroid, there are no significant differences in the spectral spread patterns among the different vowel contexts.

V. MINIMAL ABDUCTION DEGREES

In this section, two non-dimensional quantities are defined according to the method explained in Sec. III C: D 1 is the minimal abduction degree from which frication noise is generated, namely the boundary between A and B, and D 2 is the minimal abduction degree from which the noise component is predominant over the voiced component, namely the boundary between B and C. A third quantity is ∆D = D 2 -D 1 , which is the width of the A. Effect of the subglottal pressure Fig. 11 shows the median value of D 1 , D 2 , and ∆D as a function of P Sub for the three places of articulation, as well as the median absolute deviation. Increasing the subglottal pressure has for main effect to decrease D 1 , while there is no significant impact on D 2 , which remains constant. For D 1 , this is due to the fact that when P Sub is high, the lowfrequency component of the acoustic volume velocity inside the vocal tract increases as the pressure drop between the subglottal region and the mouth increases. This results in the rise of the Reynolds number, hence smaller minimal abduction degrees. The evolution of D 1 and D 2 as a function of P Sub are qualitatively similar for all places of articulation. The median value among the group of palato-alveolar fricatives is the lowest for both D 1 and D 2 . This suggests that small variations of the vocal tract geometry, and especially at the supraglottal constriction, may significantly modify the acoustic features of the produced fricative. Interestingly, ∆D slightly increases with P Sub : it is 10% for P Sub = 500 Pa, and goes up to approximately 20% for P Sub = 1500 Pa. It is also interesting to note that the values are globally very similar for each place of articulation and for each realization. This result suggests that for any articulatory configuration, the abduction range required to produce voiced fricatives depends mostly on the subglottal pressure and is not dependent on the phonological context.

B. Effect of the cross-section area of the supraglottal constriction

C. Experimental observations

This section reports experimental observations of the relationship between acoustic features of speech signal and glottal opening. For that purpose, simultaneous acquisitions of audio speech signals and glottal opening measurements, via the external lighting and sensing photo-glottography (ePGG) [START_REF] Honda | Glottal-opening and airflow pattern during production of voiceless fricatives: a new non-invasive instrumentation[END_REF] , have been performed on a series of vowel-fricative-vowel (VFV) pseudowords uttered by a French native, male speaker.

For the experiments, both ePGG and audio signals are sampled at 20 kHz. The ePGG device [START_REF] Honda | Glottal-opening and airflow pattern during production of voiceless fricatives: a new non-invasive instrumentation[END_REF] is a non-invasive technique that consists in converting the light going through the glottis, from a light source located on the surface of the side neck, into electric current, thanks to a photosensor unit located on the speaker's front neck. The glottal aperture is then deduced from the electric current generated by the photosensor, given the reasonable assumption that it is in its linear domain. Consequently, by normalizing the glottal opening by its maximal value, the glottal opening A g provided by the ePGG data can be directly related to the glottal abduction degree D ab used in the numerical study. The ePGG device has been previously successfully used to investigate glottal opening gestures (e.g. in Ref. 44 ).

The speaker is asked to utter VFV pseudowords where the vowels are the same, chosen among the three cardinal vowels /i,a,u/, and the fricatives are the 6 French fricatives /v,f,s,z,S,Z/. Fig. 13 shows an example of recorded data, with ePGG measurements, the voicing quotient V Q, and the first spectral moments S 1 and S 2 computed from the audio recordings.

It corresponds to 6 pseudowords representing one example of each fricative: /uSu/, /aZa/, /usu/, /izi/, /ufu/, /ava/. Interestingly, the acoustic features computed from the audio recording behave accordingly to the observations from the numerical study: the voicing quotient plunges at the beginning of the fricative, and then suddenly increases again at its end. During this very fast evolution of the voicing quotient, the glottal opening does not exhibit such inflections. From the V Q curves, it is possible to estimate the moment where the speaker is in the transition regime B: it starts when the voicing quotient starts to suddenly drop, and finishes at the moment when V Q vanishes, namely when the speaker enters into regime C. Symmetrically, the speaker quits C to produce B when V Q starts to rise again.

Experimental minimal areas A 1 and A 2 can then be estimated by reporting the values of the glottal opening at these instants. They have been carefully estimated by hand and are Role of glottal abductions in fricatives (Version by the authors) shown in Fig. 13 with circle (A 1 ), and square (A 2 ) marks. Again, with the sole exception of the /f/ offset, the experimental values A 1 and A 2 are similar to those from the numerical study: A 1 lies between 49 and 65%, and A 2 is between 71 and 85%. The difference ∆A = A 2 -A 1 is constantly between 14 and 20%. The critical areas A 1 and A 2 are systematically lower at the fricative offset than at the fricative onset, which suggests a hysteresis effect, not investigated by our numerical study. This is particularly true for /f/, where V Q rises at the fricative offset almost at the end of the glottal adduction movement. It is also worth noting that the glottal opening during the /ava/ utterance has not reached the critical glottal opening areas, hence a constantly high V Q, and that it denotes an alternative way to pronounce voiced fricatives, in comparison with the realization of the other fricatives shown in Fig. 13.

The behavior of the spectral moments is also interesting in regards with our numerical simulations. The spectral centroid and the spectral spread are low and constant during the vowel, as in regime A. Once the glottal opening reaches the critical opening A 1 , both spectral moments suddenly rise: the spectral centroid constantly rises up to a maximal value when the glottal opening reaches the critical value A 2 , while the spectral spread reaches a local maximum, then decreases until a constant value when A g = A 2 . This corresponds to the behavior of the spectral centroid and the spectral spread in regime B, as shown in Figs. 6 and7. This phenomenon occurs symmetrically at the fricative offset. For voiceless fricatives, when A g is larger than A 2 , i.e. in regime C, both spectral moments are stable.

D. Discussion in regards to articulatory strategies

Our simulations give a new point of view about strategies available to implement the voice/unvoiced contrasts of fricatives. Here we exploit the results of our simulations together with the observations of French speakers and German learners of French 45 . In this previous work dedicated to language learning we recorded sentences with voiced fricatives in a word final position uttered by French and German speakers, and the same words (in an isolated condition) by French speakers to investigate voicing without the influence of the initial vowel of the next word. The final devoicing of voiced fricatives of French by German learners of French allows the differences in the realization of voiced fricatives to be highlighted.

In our simulations, the computed values of D 1 and D 2 , whether as a function of the Role of glottal abductions in fricatives (Version by the authors) subglottal pressure or the cross-section area of the supraglottal constriction, lead to very small values of ∆D, namely under 20 %. In this range, corresponding to the unstable transition regime B, small perturbations of the vocal tract and the glottal configurations lead to significant modifications of the acoustic characteristics. Consequently, sustaining B to produce voiced fricatives may be very difficult for the speaker, and falling into regime C is very likely. On the contrary, voiceless fricatives is much easier because it requires to open the glottis sufficiently to be in regime C, which is very stable. However, since the abduction/adduction movement of the vocal folds is relatively slow, the speaker goes from A to C by transiting through B during a short moment, leading to voiced frames in the fricative segment. Then, in order to contrast voiced and voiceless fricatives, the speaker may have several strategies. A first strategy would consist in producing a weak frication noise, i.e. maintaining regime A for the whole duration of the voiced fricative, as for the /ava/ utterance in Fig. 13. Secondly, the speaker may stay in the transition regime B.

Since the boundary between the voiced and unvoiced regimes is almost vowel independent for each of the places of articulation (see Fig. 11), that probably enables the control of this regime to be mastered by a speaker. However, this requires to reduce the duration of the fricative, since B is very unstable. Our acoustic measurements 45 show that some French speakers (approximately one third for postalveolar fricatives) sustain voicing for the whole fricative even in a final position. This corresponds to one of these two strategies. However, there is a substantial part of French speakers who fail to produce voicing for the whole fricative even if perception tests confirmed that the corresponding fricatives are perceived as voiced. The third production strategy resorts on a less precise control of the glottal opening which results in a brief excursion in regime C, hence some unvoiced frames. The abduction/adduction movement corresponding to the devoicing and then revoicing of speech is all the easier since there is a vowel after the fricative as exhibited by Fig. 2. This probably explains why French speakers realize a release vocal schwa (/@/) after a voiced fricative in final position. This also explains that short voiceless segments of fricatives have been shown to be an acoustic clue for voicing perception [START_REF] Stevens | Acoustic and perceptual characteristics of voicing in fricatives and fricative clusters[END_REF] .

The first and second strategies have been observed in some French speakers, among a corpus of 45 speakers, for producing voiced fricatives in final position 45 . Interestingly, these strategies are used by 27 % of the speakers for /Z/, by 62 % of the speakers for /z/, and by 75% of the speakers for /v/. In regards with the results presented in this paper, this makes sense since palato-alveolar is the place of articulation that systematically presents the shorter region A ∪ B (smaller values of D 1 and D 2 ), followed by alveolars, and by labiodentals. Hence difficulties to maintain these regimes.

VI. CONCLUSION

The paper has presented a numerical study about the influence of the existence of a membranous glottal gap due to the gradual abduction/adduction movement of the vocal folds on several acoustic features of produced fricatives. Simulations used a recent glottis model that is connected with the classic 1D wave solver based on a transmission line circuit analog framework. It is, to the best of our knowledge, the first study about the acoustic impact of fine glottal configurations, such as the partial abduction and the incomplete closure of the glottis, on the production of voiced fricatives.

Simulations have highlighted the existence of three distinct regimes of fricative production, depending on the amount of frication noise that is generated. The first regime, labeled A in this paper, corresponds to an almost purely voiced signal, corresponding to an approximant consonant, where the DC component of the volume velocity in the vocal tract is too low to generate a frication noise with a significantly high acoustic level. In this regime the spectral centroid and the spectral spread are relatively low, and perturbations of the speaker configurations, such as the glottal abduction degree, or the geometry of the supraglottal constriction, do not significantly modify the spectral features. The second one, labeled B, is an unstable transition regime that corresponds to the situation where the voiced and the frication noise components of the speech signal have similar energy. In the transition regime B, the spectral centroid and the spectral spread are higher than in the first regime, because of the presence of the noise component that enforces the high frequency domain of the uttered voice. Unlike regime A, small perturbations of the speaker configurations significantly modify the acoustic features. Finally, the third regime, labeled C, corresponds to voiceless fricatives, i.e. when the frication noise component is predominant over the voiced component. In regime C, like in regime A, the spectral features are very stable in regards with variation of the speaker configurations. However, unlike in regime A, the spectral centroid is high and the voicing quotient is almost null. The existence of these regimes has been evidenced for each of the three places of articulation of French fricatives. The presence and the characteristics of these regimes have also been observed and evidenced experimentally thanks to simultaneous audio and glottal opening recordings on a real speaker.

In the articulatory-phonatory space, the transition regime B is the one with the smallest extent, confirming the fact that voiced fricatives are a difficult-to-produce class of consonants, because of the very specific aeroacoustic conditions required. Simulations have shown that the range of glottal abduction for producing voiced fricatives do not vary significantly with the place of articulation, nor with the phonological context, and not significantly with the geometry of the supraglottal constriction. On the contrary, this range varies with the subglottal pressure: the greater the pressure the longer the range of the glottal abduction degree.

In terms of articulatory strategy, B is not suitable to sustain fricatives since it is an unstable transition regime, hence several alternatives to contrast voiced fricative with voiceless fricatives: i) staying in A by favoring the voicing over the frication noise, or ii) reducing the length of the fricative segment to maximize the proportion of the amount of time in B in relation to the whole fricative segment. The latter strategy, although used by many speakers [START_REF] Stevens | Acoustic and perceptual characteristics of voicing in fricatives and fricative clusters[END_REF] , may lead to inappropriate coordination between the articulatory configurations and the configurations at the glottis. This may explain the presence of final devoicing in fricatives, which is frequently observed in many languages [START_REF] Dmitrieva | Phonological neutralization by native and non-native speakers: The case of russian final devoicing[END_REF][START_REF] Charles-Luce | Word-final devoicing in german and the effects of phonetic and sentential contexts[END_REF] , or the partial devoicing of voiced fricatives 45 .

The presented study also shows that the consideration of the partial glottis abduction in articulatory synthesis is important to simulate natural running speech. However, there is still a lack of experimental measurements of the precise time evolution of the glottal opening, which limits the use of realistic time scenarios of the coordination between the vocal tract and the glottis. This paper provides a few examples of External lighting and sensing photoglottography (ePGG) data that confirm the observations made from the numerical study.

Additional data, in addition to articulatory gestures, should be acquired in the next future in order to thoroughly study the tract-glottis coordination.

Finally, the proportion of direct and fine control by the speaker on the glottal abduction degree, due to active muscle control, over that of the fluid-structure interactions resulting from the aeroacoustic conditions is not fully comprehended yet. Possible future improvements of glottis models that account for the potential uncontrolled glottal abduction, due to fluid-structure interactions at the vicinity of the glottis, would be a step forward the full

Figure 1 .

 1 Figure 1. a) View of the glottis model, adapted from Cranen and Schroeter 13 . In this model, the incomplete closure is due to a partial abduction of the vocal folds. l g is the length of the vibrating part of the vocal fold, l ch is the length of the membranous gap and l t = l g + l ch is the total length of the vocal folds. The abduction of the vocal folds 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, the energy loss, and the air inertance inside the membranous gap respectively.

Figure 2

 2 Figure 2. a) Top: example of external lighting and sensing photo-glottography (ePGG) curve of a VFV sequence. Here is the minimal pair /aZa/-/aSa/. The phonetic segmentation is represented by vertical dashed lines. Middle: wide band spectrogram of the acoustic signal of the utterance /aZa/. Bottom: wide band spectrogram of the acoustic signal of the utterance /aSa/. b) Zoom-in view (300 ms window) of the time evolution of the glottal opening area simulated from ePGG data measured for the three places of articulation. The reference time t = 0 is the time at the maximal glottal opening.

Figure 3 .

 3 Figure 3. Area functions of the three places of articulation (from top to bottom: /S,s,f/) extracted from MRI data, for 3 cardinal vowel contexts, /i,a,u/.

  tory configurations on the simulated signals, namely the subglottal pressure, the position, and the geometry of the supraglottal constriction, each of them as a function of the glottal abduction degree D ab . It consists in simulating voice signals with static vocal tract area functions. For each simulated configuration, D ab varies from 0 to 100% with an increment step of 2%, which results, for each of the configurations, in 51 simulated signals with various glottal openings. Each simulated signal is 200 ms long, with a simulation frequency of 60 kHz, and with an abduction degree D ab that remains constant. The purpose of setting a simulation frequency to 60 kHz, which is far above the frequency domain of interest, is to avoid strong frequency warping 16 .

  0 , with F 0 = 150 Hz, are visible in the low frequency domain, and vanish in regime C. This observation suggests the existence of two stable regimes (A and C) and an unstable transition regime, i.e. B. From the variation of V Q as a function of the abduction degree, it is possible to define two quantities, called the minimal abduction, denoted D 1 and D 2 , and corresponding to the boundaries between regime A and B, and between B and C, respectively. They are computed as the two values D ab = {D 1 , D 2 } so that it leads to the best linear fitting in the three regions of the function V Q = f (D ab ), as shown in Fig. 4.

Figure 4 .

 4 Figure 4. Top: example of curve representing the voicing quotient V Q as a function of the glottal abduction degree D ab . The area function is that of the alveolar fricative in the /i/ context. The subglottal pressure value is P sub = 1000 Pa, and a c = 0.15 cm 2 . The optimized linear fitting used to estimate the minimal abductions D 1 and D 2 is represented by the dashed line. Bottom: Spectrogram representing the evolution of the simulated signal spectrum as a function of D ab .

Figure 5 .

 5 Figure 5. Voicing quotient of simulated voice signals as a function of the subglottal pressure P Sub and the glottal abduction degree D ab . Each column of figures corresponds to a place of articulation.

Figure 6 .

 6 Figure 6. Spectral centroid of simulated voice signals as a function of the subglottal pressure P Sub and the glottal abduction degree D ab . Each column of figures corresponds to a place of articulation. From left to right: palato-alveolar fricatives /S,Z/, alveolar fricatives /s,z/, and labiodental fricatives /f,v/. Each row of figures corresponds to a vowel context, /i,a,u/, from top to bottom. Contour lines at a 1000 Hz-interval are represented by solid lines.

Figure 7 .

 7 Figure 7. Spectral spread of simulated voice signals as a function of the subglottal pressure P Sub and the glottal abduction degree D ab . Each column of figures corresponds to a place of articulation.From left to palato-alveolar fricatives /S,Z/, alveolar fricatives /s,z/, and labiodental fricatives /f,v/. Each row of figures corresponds to a vowel context, /i,a,u/, from top to bottom. Contour lines at a 500 Hz-interval are represented by solid lines.

Figure 8 .

 8 Figure 8. Voicing quotient of simulated voice signals as a function of the supraglottal constriction area a c and the glottal abduction degree D ab . Each column of figures corresponds to a place of articulation. From left to right: palato-alveolar fricatives /S,Z/, alveolar fricatives /s,z/, and labiodental fricatives /f,v/. Each row of figures corresponds to a vowel context, /i,a,u/, from top to bottom.

Figure 9 .

 9 Figure 9. Spectral centroid of simulated voice signals as a function of the supraglottal constriction area a c and the glottal abduction degree D ab . Each column of figures corresponds to a place of articulation. From left to right: palato-alveolar fricatives /S,Z/, alveolar fricatives /s,z/, and labiodental fricatives /f,v/. Each row of figures corresponds to a vowel context, /i,a,u/, from top to bottom. Contour lines at a 1000 Hz-interval are represented by solid lines.

Figure 10 .

 10 Figure 10. Spectral spread of simulated voice signals as a function of the supraglottal constriction area a c and the glottal abduction degree D ab . Each column of figures corresponds to a place of articulation. From left to right: palato-alveolar fricatives /S,Z/, alveolar fricatives /s,z/, and labiodental fricatives /f,v/. Each row of figures corresponds to a vowel context, /i,a,u/, from top to bottom. Contour lines at a 500 Hz-interval are represented by solid lines.

Figure 11 .

 11 Figure 11. Top: median values of D 1 , D 2 , and ∆D as a function of P Sub for the three fricatives. Bottom: values of D 1 , D 2 , and ∆D at P Sub = 1000 Pa for each extracted area function. For each group of 7 area functions, they correspond to the following vowels /i,E,a,o,u,y,ø/, respectively.

Fig. 12

 12 Fig. 12 represents the minimal abduction degrees D 1 , D 2 , and ∆D = D 2 -D 1 , as a function of the cross-sectional area of the supraglottal constriction a c . The values of the minimal abduction degree are still systematically lower for the palato-alveolar fricatives than for the other fricatives. Although very close, they are also smaller for alveolar fricatives than for labiodental fricatives. D 1 admits a local minimum at a c = 0.1 cm 2 . D 2 does not show

Figure 12 .

 12 Figure 12. From left to right: D 1 , D 2 , and ∆D, as a function of the supraglottal constriction area a c , and for the three places of articulation.

Figure 13 .

 13 Figure13. From top to bottom: spectrogram, acoustic pressure waveform P out , voicing quotient V Q, relative glottal area A g , and the spectral moments S 1 and S 2 , computed from the audio recordings of 6 pseudowords, /uSu/, /aZa/, /usu/, /izi/, /ufu/, /ava/. Circle marks denote the boundaries between A and B, and square marks denote the boundaries between B and C.

  1 is the area of the upstream tubelet. Re = 2ρ

						µ	U DC π √ ac is the Reynolds number of the air
	flow inside the vocal tract, where µ and ρ are the shear viscosity and the mass density of
	the air, respectively.
	The computation of U DC follows the method proposed by Maeda 24 : it is the positive
	solution at each time instant t of the quadratic equation of the low frequency model of the
	vocal tract, namely
	κρ	1 a 2 gl	+	1 a 2 c	U 2 DC + 12	l 2 gl d gl a 3 gl

ACKNOWLEDGEMENTS

The work belongs to the project ArtSpeech, with the financial support of the French ANR (Agence Nationale de la Recherche). The authors are grateful to the laboratory IADI (Imagerie Adaptative Diagnostique et Interventionnelle) and the CIC-IT (Centre d'Intervention Clinique -Innovations Technologiques) of Nancy for the MRI acquisitions, and to Dr. Shinji Maeda, Dr. Angélique Amelot, and Pr. Didier Demolin for the ePGG measurements.


Kristin Mcleod

Maxime Sermesant

Xavier Pennec

Improving Understanding of Long-Term Cardiac Functional Remodelling via Cross-Sectional Analysis of Polyaffine Motion Parameters

come   L'archive ouverte pluridisciplinaire

Introduction

Cardiac function can be altered under disease states, reducing the long-term ability of the heart to pump efficiently. While there has been effort dedicated towards understanding typical functional abnormalities occurring for different diseases, the functional remodelling processes that occur over time are poorly understood. A good understanding of functional remodelling is key in optimising treatment to, for example, predict long-term functional remodelling for decision making of cardiac implants, or to determine the best timing to intervene for valve replacement surgery.

State-of-the-art in Cardiac Motion Tracking

Cardiac motion tracking via non-rigid registration has been used in recent years to quantify 3D motion of the heart. Quantitative measures of motion can be used to detect abnormal regions and to provide more detailed measures beyond simple abstractions of function such as volumes or ejection fraction. An overview of recent methods is given in [START_REF] Tobon-Gomez | Benchmarking framework for myocardial tracking and deformation algorithms: An open access database[END_REF]. Population-based analysis of motion has been used to quantify changes over time in cardiac function by performing groupwise analysis of strain derived from tagged magnetic resonance imaging [START_REF] Qian | Identifying regional cardiac abnormalities from myocardial strains using nontracking-based strain estimation and spatio-temporal tensor analysis[END_REF] or by studying spatio-temporally aligned displacements [START_REF] Rougon | Building and using a statistical 3D motion atlas for analyzing myocardial contraction in MRI[END_REF]. While informative, these methods only give a snapshot into the function of the heart and do not provide insight into how the function remodels over time.

State-of-the-art in Cardiac Remodelling

Functional changes have been studied by, for example, comparing strain values before and shortly after myocardial infarction [START_REF] Gerber | Microvascular obstruction and left ventricular remodeling early after acute myocardial infarction[END_REF], and pre-and post-surgery for aortic stenosis patients [START_REF] Bauer | Acute improvement in global and regional left ventricular systolic function after percutaneous heart valve implantation in patients with symptomatic aortic stenosis[END_REF]. However, few studies have investigated long-term functional changes, even for 1D measures of strain [START_REF] Kuznetsova | Longitudinal changes in left ventricular diastolic function in a general population[END_REF]. Analysis of changes in the heart over time is an ongoing area of research, which has traditionally been addressed using longitudinal patient data. In light of the fact that longitudinal studies require long-term follow-up of the same patients, potentially over years, the prospect of using the results of such studies to guide therapy planning remains low. In contrast, cross-sectional analysis of patients within the same population has the ability to describe long-term changes without the need to wait years for a patient to age. Long-term cardiac structural remodelling was analysed by Mansi et al. [START_REF] Mansi | A statistical model for quantification and prediction of cardiac remodelling: Application to tetralogy of fallot[END_REF] to study the long-term shape remodelling that occurs in the right ventricle of Tetralogy of Fallot patients in response to increased load caused by regurgitated blood from the lungs. This study was performed only on the end diastolic phase of the cardiac cycle, and is therefore not able to capture functional remodelling that (most likely) occurs over time in these patients.

Aim and Paper Organisation

In the present work we propose to analyse long-term remodelling of cardiac function using a method that is able to represent the full 3D dynamics of the heart. The proposed method makes use of a recent method to track cardiac motion by describing motion regionally as affine transformations, and fusing these to a global smooth transformation [START_REF] Mcleod | Spatio-temporal tensor decomposition of a polyaffine motion model for a better analysis of pathological left ventricular dynamics[END_REF]. In contrast to previous work, where the motion was analysed only at a single time point, the present work builds upon recent effort in cardiac motion tracking and population-based analysis to analyse changes in motion over time. Here, we are interested in studying motion over two time scales: the 3D dynamics during the heart beat (motion parameters), and the remodelling over years (changes in these motion parameters). The contributions of the present work are:

1. Derivation of a long-term functional model using image registration, polyaffine projection, and cross-sectional statistical analysis 2. Quantitative analysis of the accuracy of predicting age given motion parameters using leave-one-out numerical validation 3. Qualitative analysis of the accuracy of predicting motion given age by analysing motion sequences and image-derived parameters

The remainder of the paper is organised as follows: Sec 2 summarises the techniques used to track the motion of the left ventricle and project this onto a polyaffine space, in which the polyaffine motion parameters can be analysed. Sec. 3 describes the proposed method to perform cross-sectional analysis of the polyaffine motion parameters. Sec. 4 describes experiments performed to quantitatively and qualitatively evaluate the proposed method, and in Sec. 5 the methods are evaluated, discussed, and future perspectives are given.

Background: Polyaffine Motion Tracking

A global (affine) description of motion is not sufficient to characterise the motion of the heart. On the other hand, voxel-wise descriptors most likely over-represent the motion when considering that the heart consists of connected tissue which moves in an elastic manner. Describing the motion at a regional level can provide a suitable compromise between these two, to sufficiently describe the motion while representing the motion with as few parameters as possible. Polyaffine projection allows motion to be described at a regional level while maintaining a globally smooth transformation by means of weight functions that smooth the transformation between regions. Moreover, for the purpose of estimating the remodelling over time, it is important to have a robust yet informative description of the motion. We believe that the Polyaffine transformations provide a good trade-off.

Polyaffine projection from a dense deformation field was proposed by Seiler et al. [START_REF] Seiler | Capturing the multiscale anatomical shape variability with polyaffine transformation trees[END_REF] via log-affine matrices M i described in each region i that are fused to a global deformation field by v(x) = N i ω i (x)M i x. ω(x) are weight parameters in the N regions, describing the weight of the transformation of each region of a given voxel x written in homogeneous coordinates. As described in [START_REF] Seiler | Capturing the multiscale anatomical shape variability with polyaffine transformation trees[END_REF], the polyaffine projection can be formulated as a linear least squares problem by solving ∇C M = 0, where C is described by the following formula:

C(M ) = Π i ω i (x) • M i x -v(x) 2 dx.
The polyaffine projection technique used in this work was implemented following the method described in [START_REF] Mcleod | Spatio-temporal tensor decomposition of a polyaffine motion model for a better analysis of pathological left ventricular dynamics[END_REF], which was developed specifically for motion tracking of the left ventricle of the heart.

Methods: Regression of Polyaffine Motion Parameters

The compact parameterisation of motion parameters using the polyaffine formulation described in the previous section enables cross-sectional analysis of motion. Common methods for cross-sectional analysis include principal component analysis (PCA) and partial least squares (PLS). PLS has the advantage over PCA of computing directly the components most related to the output variable. In both cases, the regression model is performed by computing the motion parameters from age since the converse requires prediction of many output parameters: motion described by 5916 parameters per subject (12 affine parameters ×17 regions ×29 time frames), given few input parameters: age described by a single value (in years). Following the method for computing shape features from age proposed by Mansi et al. [START_REF] Mansi | A statistical model for quantification and prediction of cardiac remodelling: Application to tetralogy of fallot[END_REF], canonical correlation analysis (CCA) can be used to reverse the relationship in order to obtain a model of motion given age. The proposed pipeline is summarised in Fig. 1.

Projec'on to reduced-order space Cross-sec'onal analysis of 3D mo'on Voxel-wise nonrigid mo'on tracking

Fig. 1: Proposed pipeline from image to motion model. Non-rigid motion tracking is performed to describe the voxel-wise motion in the heart. This is then projected to a polyaffine space in order to represent the motion with few parameters, consistently defined from one subject to another. Finally, cross-sectional analysis is performed on the motion parameters to formulate a generative model of functional remodelling.

Experiments: Polyaffine Motion Parameter Regression

Two types of experiments were performed to validate the proposed method. The first experiment was a quantitative validation of the accuracy of predicting age given the motion parameters. While this is not a clinically meaningful validation (since the age of patients is generally known), it serves as a means to validate the regression since it is straightforward to compute the accuracy when predicting a single output value. The second experiment was a qualitative assessment of the prediction of motion given age, which is difficult to validate quantitatively given the large number of predicted parameters, but can be analysed qualitatively by visualising the output to determine if it matches physiological motion dynamics.

Testing Datasets

The motion parameters used in the experiments were taken directly from the previous work of [START_REF] Mcleod | Spatio-temporal tensor decomposition of a polyaffine motion model for a better analysis of pathological left ventricular dynamics[END_REF], which consists of two datasets; a healthy control dataset, and a pathological group of Tetralogy of Fallot patients. The healthy control subjects from the openly available STACOM 2011 cardiac motion tracking challenge dataset [START_REF] Tobon-Gomez | Benchmarking framework for myocardial tracking and deformation algorithms: An open access database[END_REF] (n = 15, 3 female, mean age ± SD = 28 ± 5) were used since these have already been used to validate the polyaffine cardiac motion tracking. The pathological subjects were a dataset of Tetralogy of Fallot patients (n = 10, 5 female, mean age ± SD = 21 ± 7), which were already used to analyse the polyaffine motion parameters in [START_REF] Mcleod | Descriptive and intuitive population-based cardiac motion analysis via sparsity constrained tensor decomposition[END_REF]. The Tetralogy of Fallot cohort forms a suitable testing set for this problem since there are known functional changes over time in these patients in response to poor pulmonary valve function or in the absence of a pulmonary valve. Furthermore, modelling functional remodelling in this population complements the previous analysis of shape remodelling performed by Mansi et al. [START_REF] Mansi | A statistical model for quantification and prediction of cardiac remodelling: Application to tetralogy of fallot[END_REF]. In both cohorts, only left ventricular motion was tracked due to the absence of open testing sets to validate bi-ventricle motion tracking methods. The motion tracking results from these experiments were consistent with those presented in [START_REF] Mcleod | Spatio-temporal tensor decomposition of a polyaffine motion model for a better analysis of pathological left ventricular dynamics[END_REF].

Statistical Analysis

Following the techniques described in [START_REF] Mansi | A statistical model for quantification and prediction of cardiac remodelling: Application to tetralogy of fallot[END_REF] to predict shape remodelling, partial least squares (PLS) was applied to each vectorised motion array to reduce the parameter space to a vector of PLS component loadings, from which CCA between the loading vector (representing motion) and age was computed. CCA was performed iteratively on the number of components, to find the minimal number of components required to obtain a statistically significant combination of PLS components with respect to the dependent variable. Statistical significance was computed using the Bartlett-Lawley hypothesis test. Since we were interested in modelling motion (and not age), the important variable to explain is Y (vector of Polyaffine motion parameters from M ), rather than X (age). PLS captured more than 90% of the variability in Y with only 2 components for the ToF group, and 7 components for the control group. The canonical correlations are plotted for each component from PLS of each group in Fig. 2. The canonical correlations for PLS are ordered in magnitude by design since PLS is maximising the variance in Y and covariance in X and Y .

Quantitative Validation of the Prediction of Age Given Motion

The testing error of predicting age given motion was computed by performing leaveone-out validation (training on n -1 subjects, and predicting the age of the n th subject. In all error computations, the number of CCA components for the full training set were used. The training error of the Tetralogy of Fallot group was 0.30 years and the testing error was 5.85 years. For the control group, the respective training and testing errors were 0.47 years and 6.81 years. For the ToF group, the yielded errors were below the population standard deviation of age (6.49). For the control group, however, the testing errors were above the population standard deviation of 5.0 years. needed for each group to obtain statistically significant combinations of the modes. Fewer components were required in the more heterogeneous pathological group (i.e. the ToF group) than relatively homogeneous control group.

Qualitative Analysis of Prediction of Motion Given Age Quantitative validation of the motion parameters computed from age would require longitudinal data, which is a challenge to obtain given the long time intervals required between time points. Furthermore, quantitative validation voxel-wise is challenging due to the fact that the heart remodels structurally over time in addition to the functional remodelling. Since there is currently no method to identify how/where voxels move over time as a result of structural remodelling compared to those that move as a result of functional remodelling, quantitative validation of the motion estimation is not possible. Therefore, only qualitative analysis of motion via visualisation of the motion sequences was performed. Fig. 3 shows the predictive (healthy or pathological) deformation models applied to a single magnetic resonance image using the motion visualisation method described in [START_REF] Mcleod | Spatio-temporal tensor decomposition of a polyaffine motion model for a better analysis of pathological left ventricular dynamics[END_REF], shown here at two phases of the cardiac cycle to compare the remodelling of motion in the Tetralogy of Fallot group to the control group. As expected, the control group shows little changes in motion through ageing, since the population is healthy the motion should remain stable over time. In contrast, the Tetralogy of Fallot group show large differences in motion over time, particularly in the free wall.

As a preliminary means of quantitative validation, the predicted motion sequences for the Tetralogy of Fallot group were segmented to obtain the left ventricular volumes, plotted in Fig. 4. As the curves show, the minimum of the curves increase with age, which ultimately decreases ejection fraction; a commonly used measure of function by cardiologists. Ejection fraction is known to decrease with age after initial repair in these patients [START_REF] Waien | Serial follow-up of adults with repaired tetralogy of fallot[END_REF].

Discussion and Perspectives

The experiments suggest that the proposed methods are promising for computing the long-term functional changes in the left ventricle. The computation of age from motion parameters was more accurate in the Tetralogy of Fallot group than the healthy control group, which is as expected since the control group consists of a range of ages in which the motion parameters should be reasonably stable (thus making regression challenging). The qualitative validation suggests that the methods are able to predict feasible motion dynamics from age alone, despite the small sizes of the datasets. These findings should be, however, confirmed in a larger population in which longitudinal data is available to quantitatively validate the prediction of motion from age. Additionally, with larger datasets the models could be sub-grouped to ensure that the change due to the evolution is much larger than the cross-sectional variability.

An important limitation of the present work is that no longitudinal data was used for the regression. As mentioned by Eng et al. [START_REF] Eng | Adverse left ventricular remodeling and age assessed with cardiac mr imaging: the multi-ethnic study of atherosclerosis[END_REF], cross-sectional analysis can underestimate size changes over time compared to longitudinal analysis. Since longitudinal data was not available at the time of this study they were not used, therefore future studies will focus on a combined longitudinal and cross-sectional study design.

This work presents only functional changes without consideration of structure, despite the fact that these are inherently coupled and could be more informative when used in combination. A very interesting future work will address coupling shape (for example using the methods in [START_REF] Mansi | A statistical model for quantification and prediction of cardiac remodelling: Application to tetralogy of fallot[END_REF]) with function.

In the present work, regression was performed on matricised tensors. However, the regression could be performed directly on the tensors. Performing re- gression directly on the tensor may be an advantage for retaining the structure of the tensor in the regression, and to reduce the number of parameters required to describe the generative motion models.

Conclusions

In this work we propose a method to analyse motion changes in a population over time using a non-rigid registration algorithm in conjunction with polyaffine projection and statistical analysis. The proposed method was validated quantitatively and qualitatively on a control dataset and on a dataset of Tetralogy of Fallot patients and complements the work of [START_REF] Mansi | A statistical model for quantification and prediction of cardiac remodelling: Application to tetralogy of fallot[END_REF] to allow future potential work on joint shape and motion evolution modelling. This is, to our knowledge, the first method that is able to predict long-term functional remodelling at a level that retains the full complexity of motion dynamics.

Fig. 2 :

 2 Fig.2: Canonical correlations for PLS for the Tetralogy of Fallot group (left) and the control group (right). The red line indicates the minimum number of comneeded for each group to obtain statistically significant combinations of the modes. Fewer components were required in the more heterogeneous pathological group (i.e. the ToF group) than relatively homogeneous control group.

Fig. 3 :

 3 Fig. 3: Snapshots of the sequences shown at two phases of the cardiac cycle: frame 10 (top row) and frame 20 (bottom row) of a total of 30 frames for the Tetralogy of Fallot group (left) and healthy control group (right). Three ages are shown for each group at -2 standard deviations of the age, the mean age, and +2 standard deviations of age. Note that the motion in the control group remains stable over the ages, since very little motion changes occur over normal ageing. In contrast, noticeable changes in the motion of the Tetralogy of Fallot group are visible, as highlighted by the yellow arrows.

Fig. 4 :

 4 Fig. 4: Segmented volumes computed from sequences predicted from the Tetralogy of Fallot group model. Ejection fraction (used as a measure of cardiac function) decreased over time, consistent with clinical findings in this patient group.

Acknowledgements: This project was carried out as a part of the Centre for Cardiological Innovation (CCI), Norway, funded by the Research Council of Norway. The authors would like to thank Tommaso Mansi for providing the code for computing CCA.