Voice synthesis using power-balanced simulation of a quasi-1D model of the vocal apparatus

Abstract

The vocal apparatus is a biophysical dynamic system capable of self-oscillation, which involves fluid–structure interactions and human control. This study on the sound synthesis of voiced sounds presents a physical quasi-1D model of the vocal apparatus in the port-Hamiltonian framework and its validation through numerical experiments. The modelling ensures balanced power exchanges between fluid, tissues, and human control. Fluid is represented in the larynx and in the vocal tract using a unified 1D PDE handling transverse geometry variations. A regularisation procedure is introduced to mitigate the numerically stiff behaviour of the model observed at channel closure. Vocal folds and vocal tract walls are represented by lumped element models as well as the radiation load at the lips, which consists of a first-order high-pass filter. Spatial discretisation of the fluid model and temporal discretisation of the full system are made using structure-preserving methods to ensure energy consistency (passivity). The second part of this paper focuses on numerical experiments to progressively characterise the model and assess its validity. These experiments begin with frequency response analysis of a static vocal tract under quasi-linear conditions followed by simulations of vowel transitions (diphthongs) under forced excitation. Next, self-oscillation studies are conducted on an isolated larynx where contact parameters are adjusted. Lastly, full simulations of the self-oscillating vocal apparatus with co-articulation, representing a voice synthesizer capable of articulating vowels, are presented. The dynamics are also analysed in terms of energy transfer and passivity. Finally, these results are discussed to establish a basis for future model refinements and to identify directions for enhancing the accuracy and realism of vocal synthesis.

1 Introduction

Understanding the mechanisms of voice production has been the focus of numerous studies on the physical modelling of the vocal apparatus. Highly detailed models that capture the range of physically relevant phenomena typically involve high-dimensional representations, which can be computationally intensive. Consequently, these studies often focus on specific components or aspects of the vocal apparatus (e.g., Alipour et al., 2000; Gunter, 2003; Xue et al., 2010; Guasch et al., 2013; Zhang, 2015; Valášek, 2021). In contrast, low-order models aim to replicate the qualitative behaviour of the vocal apparatus along with certain quantitative features (such as the relationship between subglottal pressure and fundamental frequency, as demonstrated in Ruty et al. (2007)) while maintaining relatively low computational costs. This computationally efficient approach has been a key focus since the early days of sound synthesis, with foundational models like the Kelly–Lochbaum model (Kelly, 1962) or digital implementations of lumped-parameter physical models (e.g., Maeda, 1982) for the vocal tract. These models are well-suited for vocal synthesizer applications or parameter studies. For the laryngeal part, low-order models (e.g., Ishizaka and Flanagan, 1972; Story and Titze, 1994; see also the review in Erath et al., 2013) typically represent the vocal folds as assemblies of lumped mass-spring-damper elements and assume a quasi-steady behaviour of fluid flow in the glottis. However, this approach introduces a one-way coupling between the flow and the vocal folds, limiting the model’s ability to accurately represent power exchanges between the fluid and tissue, and it may introduce instabilities in numerical simulations.

Recent efforts have focused on developing power-balanced, reduced-order models of voice production within the port-Hamiltonian (pH) framework first introduced in Maschke and van der Schaft (1992). Notably, Encina et al. (2015) reformulated the body-cover vocal fold model from Story and Titze (1994), and early self-oscillating power-balanced models were introduced in Hélie and Silva (2017) and Mora et al. (2018) using lumped representations of glottal flow. Subsequent research by Mora et al. (2021b) and Wetzel (2021) proposed discrete, scalable models of compressible fluid within tubes with moving boundaries, tailored for vocal tract applications. However, this set of studies primarily focused on isolated larynx models (possibly coupled with an oversimplified acoustic load) or on isolated vocal tracts under basic configurations.

This paper presents a complete simplified power-balanced model of the vocal apparatus within the pH framework, designed for the synthesis of voiced sounds. A one-dimensional distributed fluid model is developed and described as a port-Hamiltonian system of partial differential equations, accounting for compressibility, time-varying geometry, and convective acceleration. A new regularisation procedure for handling glottal closure is introduced. This model is then spatially discretised in a structure-preserving manner, from which the discrete fluid model from Mora et al. (2021b) is recovered. The vocal fold model is based on the lumped-element body-cover model from Story and Titze (1994), while the vocal tract walls are represented by a series of simpler mass-spring damper systems. At the lips, the radiation condition is modelled as an impedance load, represented by a first-order high-pass filter, tuned according to acoustic considerations and designed to be passive and compatible with time-varying configurations. Each component is formulated within the pH framework and interconnected via ports using a two-way coupling, ensuring a global power balance. Simulations are performed using an energy-preserving time integration method across various configurations, including the vocal tract alone, the larynx alone, and the complete assembly. Notably, the ability of the complete model to reproduce the co-articulation of diphthongs in self-oscillating configurations is demonstrated—an achievement that, to the best of our knowledge, is unprecedented in the context of energy-based modelling of voice production. Energy-based modelling explicitly reveals the power exchanges and dissipated power signals, which are discussed in detail.

After positioning our contributions within the existing literature on voice production modelling in Section 1, the remainder of this paper is organised as follows. Section 2 introduces the fluid and tissue models used in the experiments. Section 3 presents a series of numerical experiments based on these models, progressively increasing in complexity and realism, to assess their validity and capabilities. Finally, Section 4 discusses the results and outlines potential directions for future work.

2 Model

Energy-based models for fluid dynamics and vocal folds are proposed within the port-Hamiltonian framework. These models are subsequently interconnected, resulting in power-preserving fluid–structure interactions. Symmetry around the laryngeal midplane is assumed, allowing reduction of the problem to the case of a single vocal fold and simplifying the presentation.

2.1 Fluid flow model

The modelling of the fluid flow is proposed within the framework of port-Hamiltonian systems. It relies on the strict separation of the intrinsic behaviour of the fluid (described by means of the Hamiltonian, summing up kinetic energy and internal energy, or more specifically by the fluid behaviour law relating the pressure to the mass density) from the evolution laws (provided as mass and momentum conservation laws). As an original contribution of this paper, the latter are formulated as 1D partial differential equations under a detailed set of assumptions. Spatial discretisation is applied to obtain a discrete model, similar to that used in Mora et al. (2021b), with the change of state variables proposed in Risse et al. (2024). Here, the model with this revised set of state variables is presented. The model is constructed as a conservative system (the energy conservation naturally derives from the port-Hamiltonian systems), and lumped dissipation laws are added after spatial discretisation.

2.1.1 Formulation as a continuous 1D PDE system

The time-varying fluid domain is described in 3D cartesian space as follows:

Figure 1 presents the domain as well as its boundaries. Boundaries and represent the inlet and outlet of the tube through which fluid can flow. In contrast, denotes a rigid, sealed, and impermeable boundary, while corresponds to a mobile impervious fluid–structure interface.

FIGURE 1

The following assumptions are made for the fluid in :

2.1.2 Spatial discretisation: finite differences on staggered grids

The 1D continuous equation is semi-discretised using a power-preserving finite difference approximation on a staggered grid (Trenchant et al., 2018). For the sake of brevity, only the configuration (C1) where mass flows inputs are considered at both ends of the domain is detailed here, but other configurations can be obtained using the same procedure. The 1D mesh is first presented, and then field approximations are made on that mesh. The resulting dynamics obtained from injection of these approximations in Equation 8 and integration on mesh volumes are presented using a set of chosen discrete state variables.

2.1.2.1 Meshing of the 1D domain

A total of edges of lengths are naturally obtained for the dual grid (Figure 2).

FIGURE 2

2.1.2.2 Approximation of the fields

2.1.2.3 Choice of state variables and discrete Hamiltonian

The state is chosen to bewith the degrees of freedom of velocity on the primal grid and of masses of fluid in the volumes of the dual grid. This choice is motivated by the simplicity of the port-Hamiltonian formulation resulting from it and has no impact on the dynamics of the system.

The energy of the discrete fluid system is obtained by injecting field approximations in Equation 6, yieldingmaking use of functions of the state variables to reconstruct the volumetric mass densities in volumes of the dual grid and primal grid and the masses of fluid in volumes of the primal grid . We also introduce the function of the state reconstructing the mean squared velocity in volumes of the dual grid . Effort are naturally obtained as the gradient of the Hamiltonian:and provide the mass flow rates averaged over the volumes of the primal grid, the enthalpy averaged over the volumes of the dual grid, and the net reaction force of the fluid at the fluid-structure interfaces of the primal grid.

2.1.2.4 Semi-discrete dynamical system

Injecting the field approximations in the 1D conservation equations and integrating Equation 3b on the edges of the primal grid and Equation 3a on the edges of the dual grid yields the semi-discrete dynamical system. one obtains the port-Hamiltonian formulation:where is the input vector composed of the ingoing mass flow at , outgoing mass flow at , and fluid structure interface velocities . The finite difference matrix and the input matrix are written as

The power-conjugated outputs of the system are obtained by skew-symmetry as , such that .

2.1.3 Lumped dissipation laws

The dynamics of the fluid system are modified to account for these phenomena through the introduction of the dissipation matrix where is the dissipation vector. Identification with Equation 17 yieldsfor friction forces, to be completed by an additional term when for the kinetic energy loss term at the exit of the glottis (at the -th term of ). Finally, the acoustic radiation model presented in Supplementary Appendix SA4.3 is connected in a power-balanced manner to the fluid model at the exit of the vocal tract, writing and .

2.1.4 Equivalent circuit representation of a section of the semi-discretised model

An equivalent circuit representation of a portion of the fluid domain (more precisely to an edge of the primal grid) is depicted in Figure 3, evidencing two sub-circuits. The lower sub-circuit accounts for axial inertia and mass accumulation as in classical transmission line structures, while the upper one handles the reaction of the fluid to a change of geometry (connected to a model of tissues to be described in the next sections). Even if the two sub-circuits seem to be disconnected, they are in reality coupled through the multivariate non-separable Hamiltonian that intrinsically links the various energy-storing components.

FIGURE 3

2.2 Vocal fold model

The classical body-cover model of the vocal fods presented in Story and Titze (1994) is used. It consists of a body mass coupled with two cover masses. The port-Hamiltonian formulation has been developed in Encina et al. (2015).

2.2.1 Presentation

The assembly, presented in the left section of Figure 4, comprises the three masses coupled together by four springs and three dampers. Derivation of the port-Hamiltonian formulation of this assembly is recalled here. The body mass describes the displacement of the in-depth vocalis muscle while the lower and upper cover masses and correspond to movements of upstream and downstream parts of the superficial layers (vocal ligament and epithelium). The state of this system, identified by the subscript for “tissues”, can be chosen aswhere , , and are the displacement of the three masses around the equilibrium positions, and , , and are the momentum of the three masses. From the state, potential energy functions and total potential energy are denoted:

FIGURE 4

The resultant force of the springs on each mass are directly recovered by differentiation with respect to masses displacements:

The total kinetic energy of the system is computed assuch that the Hamiltonian is given by

The port-Hamiltonian formulation of the dynamics including linear dissipation effects and external interactions iswithand where dissipation coefficients , , and are associated with linear dampers. In the numerical experiments, these coefficients are tuned from a single damping ratio as , , and . The inputs are the external forces applied on the lower and upper cover masses, with the associated input matrix:

The power-conjugated output vector contains the velocities and of the cover masses.

2.2.2 Setting the characteristics of the springs

The above formulations leave us with the choice of potential energy functions . State-of-the-art lumped models of the vocal folds commonly use nonlinear springs with an added third order term (e.g., Ishizaka and Flanagan, 1972; Story and Titze, 1994,), such that the potential energies, restoring forces, and effective stiffness of individual springs arewhere is a reference elongation for the nonlinearity, such that the effective stiffness is multiplied by a factor of four for .

2.3 Vocal tract tissue

A better representation of low-frequency resonances of the vocal tract requires the introduction of a model for tissues of the vocal tract. We here consider a simplified local mechanical wall surface impedance , where , , and are, respectively, the equivalent unit surface resistance, mass, and stiffness. Their values are adapted from those given in Birkholz et al. (2022). Parameter values for lumped components depicted in the right section of Figure 4 are naturally obtained by integration of their unit surface equivalent on the fluid-structure interface surfaces (, , ). The outer wall velocity is an input of the wall dynamical model, such that the overall geometry can still be controlled. The port-Hamiltonian formulation of this model, ready for connection to the fluid model, is presented in Supplementary Appendix SA4.4.

2.4 Assembly

2.4.1 Fluid–tissue connection

In the context of the numerical experiments, the models presented are tested in various configurations. The vocal tract will consist of the fluid model alone, the geometry of which will be controlled by inputs defined in order to produce articulations of different vowels. In the larynx, the fluid model will be connected to the vocal fold model by assuming the continuity of normal velocities and net forces at the fluid–structure interface. Note that several fluid cells may be connected to a single mass of the vocal folds. In this case, connected fluid cells share the same geometrical velocity, and the force received by the vocal fold mass is the result of pressure forces of all fluid cells connected to it. The resulting assembled system can be written as a new port-Hamiltonian system, the state of which is composed of the concatenation of the states of the fluid and the vocal folds.

2.4.2 Glottal closure

In many cases, self-oscillations of the vocal folds include phases of glottal closure and thus contact between the vocal folds. This phenomenon is handled in the modelling by introducing direct coupling between the vocal folds model that is only activated when the vocal folds are detected to be in contact—when the heights of the fluid elements connected to the vocal fold masses are negative. The implementation is straightforward by adding the compression energy of these springs to the energy of the assembled system without any change in the matrices. However, the fluid modelling described in Section 2.1 prevents the full closure of the channel, as the reaction force diverges to infinity when goes to 0, making the system stiff—which makes time integration difficult.

Mora et al. (2021b) proposed a switching port-Hamiltonian model that disconnects the fluid from the vocal folds when , effectively taking the stiffness of the fluid layer away from the dynamics of the fold during contact. Considering the anatomical 3D geometry of the glottal channel and the imperfect contact between vocal folds in the human larynx (possibly with glottal leakage), we propose an alternative non-idealised approach to model contact without the need for logical switches checks during the simulation. Considering a given distance threshold , an effective height of the glottal flow channel is introduced, such that , with specifying the smooth transition width. The Hamiltonian for the fluid domain (Equation 13) is then modified to account for the effective glottal heights: .

Moreover, the discrepancy can be interpreted as a measure of compression of the vocal fold tissues and used within a modified Hamiltonian for the vocal fold model , with to account for the additional stiffness term corresponding to fold interpenetration. Additional linear dissipation is also added to the folds when they are interpenetrating, taking the form of an augmentation of the damping ratio to its critical value when . The effect of parameters , , , and will be discussed in the numerical experiments.

3 Numerical experiments

The energy-consistent model established in the previous section as a port-Hamiltonian system is now tested and discussed through a sequence of numerical experiments of increasing complexity. The final experiment simulates the full vocal apparatus exhibiting self-oscillations, where vowel articulation is achieved through geometrical control of the vocal tract.

3.1 General numerical setting and experimental configurations

For time discretisation, a continuous Galerkin time-stepping method is used which preserves a discrete power balance and is unconditionally stable. The method is presented in Chapter 5 of Müller (2021). Second-order polynomials are used for the projection. All simulations are run at a standard audio sampling rate .

3.2 Frequency response

In the first experiment, the frequency response of the fluid model in static configurations is studied. Firstly, a straight tube of total length m, homogeneous height m, and with rigid walls is considered. The model is excited by a mass flow impulse of amplitude on the left end side. On the right end side, a zero total enthalpy is prescribed (no radiation) so that the waveguide should theoretically exhibit resonances at frequencies under acoustic approximations (small-amplitude fluctuations of velocity and enthalpy). As the proposed model is nonlinear, the results of the time-domain simulation are analysed and the frequency response is estimated from the ratio of the spectrum of the right mass flow to that of the prescribed impulse of left mass flow , evidencing as many resonances as the number of edges. The convergence of resonance frequencies with respect to the spatial discretisation step is studied and shown in Figure 5, with a convergence of order 2.

FIGURE 5

Static vocal tract configurations are now considered. The geometry is set to correspond to a vowel articulation; area functions from Story et al. (1996) are interpolated to produce a set of equidistant cross-section areas from which the values of used for the simulation are determined. The total length is set to the vocal tract length corresponding to the simulated vowel. The soft wall model is now used, as well as the radiation condition at the lips. Figure 6 presents the discretised area function for the vowel /ɑ/ and the associated equivalent impedance and mass flow transfer function. Table IV of Story et al. (1996) presents the first three formant frequencies for different vowels, extracted from audio recordings of the individual whose vocal tract geometry was measured. Figure 7 compares these measured values to those obtained from simulations of the model for different vowels. Values obtained from simulation show significant differences with those measured, indicating the need for a better loss model. In particular, and as pointed out in studies such as Fleischer et al. (2015), we observed that the vocal tract wall model shifts up the first formant frequency. The choices made for the fluid domain description in this paper mean that only a part of the fluid domain boundary (corresponding to in Figure 1) corresponds to a fluid–structure interface, reducing the effective area of mobile vocal tract walls. An alternate formulation of the fluid model for arbitrary shape of cross-sections, which would remove this issue, is briefly presented in Supplementary Appendix SA4.5.

FIGURE 6

FIGURE 7

Story et al. (1996).

3.3 Diphthong articulation

This experiment presents the use of the model as a vocal tract with forced geometry to reproduce co-articulation of diphthongs. The general configuration is similar to the previous one, with a total length m divided into edges. A target sequence of vowels is chosen. From this sequence, a target time trajectory for is derived by linear interpolation between individual vowel geometries and smoothed using an averaging filter. The corresponding input velocity to control the model is obtained using a finite difference operator, ensuring that the effective trajectory of obtained by simulation is the same as the target trajectory. An example target sequence for diphthong /ɑo/ is given at the top of Figure 8, presenting the wideband spectrogram of the radiated pressure for a simulation corresponding to this trajectory and using a synthetic glottal flow model (Klatt & Klatt model, e.g., Doval et al., 2006). The spectrogram exhibits the time evolution of formant frequencies as the geometry is modified.

FIGURE 8

3.4 Isolated larynx

One of the main objectives of this model is to be able to reproduce self-oscillations of the assembly composed of the laryngeal fluid flow and the vocal folds. We first study a configuration corresponding to an isolated larynx, with no vocal tract model or load. Note that this type of configuration is similar to that presented in Mora et al. (2021b) and can be interpreted as an excised larynx experiment. We use a total of segments for the discretisation of . Segments are connected to the lower mass of the vocal folds and segments to the upper mass. Segments 0 and 7 are fixed to a height of . The geometry of the fluid system at rest is defined as

The subglottal total enthalpy is prescribed to be a smoothed step function (as shown in the top axis of Figure 10). An ideal boundary condition is set at the supra-glottal end, enforcing a zero total enthalpy.

3.4.1 Oscillating configuration

Physical parameters from Table 1 are used for simulation, corresponding to case C of Story and Titze (1994). Using their value of for the damping ratio, the model is unable to oscillate for realistic sub-glottal pressure values. This can be expected as no vocal tract loading is considered in this experiment. However, Story and Titze (1994) mention that estimated values of the damping ratio range from 0.1 to 0.4. Figure 9 presents a cartography of simulation points corresponding to varying values of the sub-glottal pressure step, the approximated value of which is computed from the sub-glottal enthalpy as and of the damping ratio . In the studied range of sub-glottal pressures, oscillations exist only for .

FIGURE 9

3.4.2 Signal analysis

Figure 10 shows the distances of the simulated masses to mid-plane and mass flow under the cover masses, obtained for the configuration indicated by a diamond marker in Figure 9. During the first s, the oscillation builds up to reach a periodic regime, lasting until the input enthalpy is decreased. Glottal closing happens both under the lower and upper superficial masses when their positions are negative. As expected from the high numerical value of , the body mass oscillates with a smaller amplitude. The bottom axis of the figure presents a zoom on one period of the oscillation. A clear phase difference is visible between displacements of the lower and upper masses. The model accounts for compressibility and unsteady effects in the fluid. Mass flow signals are similar under both masses, which is expected as the fluid flow in the larynx is, in fact, almost incompressible and stationary. Nonetheless, nonstationary effects related to geometry variations do show a visible effect, such as on the portion of the cycle where the channel is closed under the lower mass but not under the upper one. The mass flow then vanishes under the lower mass but not under the upper mass as its displacement still pushes fluid out of the larynx.

FIGURE 10

FIGURE 11

The time-averaged powers are plotted on the right column of Figure 11. After establishing the periodic regime, it is clear that mean variations of the Hamiltonian are null, as expected. Kinetic energy dissipation is the most prominent source of energy loss whereas viscous losses and dissipation in the folds are of the same order of magnitude. From the energy budget of the larynx, the power received from the lower airways is mostly used to install a glottal flow able to produce flow-induced oscillations at the expensive cost of strong kinetic energy losses at the glottal exit. Once the flow is installed, a small extra effort is required to provide the power needed to compensate for viscous losses within the glottis and damping in the tissues.

3.4.3 Numerical experiment on the effect and choice of height regularisation parameters

The parameters and specifying the regularisation of the constriction height (Section 2.4.2) do not have a direct physical interpretation, making their values difficult to set. In the following, a choice is made that reduces the numerical stiffness and preserves some important characteristics of the signals. Figure 12 presents snapshots of three simulations corresponding to decreasing values of , all other parameters being fixed to the same values as for the previous simulation.

FIGURE 12

It is preferable to allow folds to interpenetrate as it allows for better control of their behaviour during contact by tuning and , effectively reducing the stiffness of the numerical problem. When the channel is closed, fluid flow should be stopped. Here, the value of m is used for the simulations as it yields qualitatively good behaviour at channel closing by letting folds interpenetrate and stopping fluid flow. Further studies are necessary to assess the impact of this choice on phonation pressure thresholds and to find an optimal value from comparison with measurements or higher order models.

3.5 Larynx connected to vocal tract

This last experiment focuses on the simulation of the larynx coupled to the vocal tract. This effectively represents a simplified vocal apparatus capable of self-oscillations and articulations. A total of fluid segments are used, with the first seven segments corresponding to the larynx and the remaining 20 representing the vocal tract (Figure 4). Segments 1 to 3 are connected to the lower mass of the vocal folds, and segments 4 to 6 are connected to the upper mass. The lengths of the segments are defined as follows:

This configuration reproduces the concatenation of the larynx and vocal tract fluid models used in the previous experiments. The soft wall model is used in the vocal tract, and the boundary condition at the lips is set using the high-pass radiation model. Like the isolated larynx case, the sub-glottal total enthalpy serves as the control input and is set to a smoothed step function. The vocal tract geometry follows the trajectory defined in Figure 8 to produce the diphthong /ɑo/. Initial laryngeal heights are retained from Section 3.4 (see Equation 29), the damping ratio for the vocal folds set to , and other physical parameters are provided in Table 1.

Time-domain results are presented in Figure 13. The top axis of the figure displays the mass flow and cover mass displacement signals over the full simulation duration. Oscillations build up at the start of the simulation, followed by a stable periodic regime, which is then modified by the transition of the vocal tract shape from /ɑ/ to /o/ between 0.4s and 0.6s. After 0.6s, the vocal tract geometry remains fixed, leading to another stable periodic regime. The output mass flow amplitude, shown in grey in the middle row, is significantly altered after the change in vocal tract shape. This change is likely due to the substantial modification in the vocal tract area at the lips, between the wide opening for /ɑ/ and the narrow opening for /o/.

FIGURE 13

The middle axis of the figure depicts zoomed-in views of two stable periodic regimes corresponding to vowels /ɑ/ and /o/. The mass flow signals in the glottis are strongly influenced by the vocal tract, as shown by comparisons of the respective waveforms. In addition, the fundamental frequency is noticeably lower for /o/ (99 Hz) than for /ɑ/ (119 Hz). However, the displacement of the masses is less perturbed by the acoustic coupling with the vocal tract than the mass flow signals.

Averaged power signals are shown at the bottom of the figure, with the averaging filter length set to ten periods of oscillation for configuration /ɑ/. Unfiltered signals, containing high-frequency content (amplified by the vocal tract), are not displayed; however, the global power balance is verified up to machine precision at each time step. The distinct physical origins of dissipated power—jet dissipation at the larynx output, viscous fluid losses, dissipation in the folds, radiated power at the lips, and dissipation within the vocal tract walls—have varying contributions. In particular, viscous fluid losses and fold dissipation are more important here than in cases without a vocal tract (Figure 11). Moreover, the relative contributions of each dissipation mechanism shift substantially after the vocal tract shape modification—notably the increase of the dissipated power related to viscous effects in the glottis and the vocal tract. Further investigation is needed to explain this effect.

4 Discussion and perspectives

This study demonstrates the effectiveness of a quasi-1D port-Hamiltonian model in replicating key characteristics of voice production, ensuring well-posed energy transfer dynamics between fluid, tissues, and human control. The model exhibits self-oscillatory laryngeal behaviours and incorporates a variable-geometry vocal tract. It thus serves as a simple voice synthesizer capable of articulating vowels. This positions the approach as a basis for developing a powerful tool for voice synthesis. This statement is discussed in two parts: firstly, through an analysis of the proposed model, and second, by considering motivated refinements for future model development.

4.1 Simulations of the proposed model

In Section 3, we presented several simulations of the proposed model corresponding to different configurations.

Vocal tract simulations show the evolution of formant frequencies with variations of geometry of the conduit. The incorporation of dissipative effects, such as radiation losses and wall vibrations, proved crucial for obtaining realistic formant bandwidths and frequencies—in line with Fleischer et al. (2015) and Birkholz et al. (2022). Story et al. (1996) provides a set of MRI measurements of vocal tract area functions corresponding to vowel sounds, as well as the corresponding first three formant frequencies extracted from acoustic recordings. Significant error is observed between these measured formant frequencies and our model predictions, which seems related to our choice of vocal tract wall and radiation parameters. The discrepancies observed between these measured formant frequencies and the model predictions—likely stemming from parameter choices related to vocal tract wall damping and radiation losses—highlight the need for further research to refine the modelling of dissipation mechanisms or to optimize these parameters for enhanced realism.

Simulations of the larynx alone and of the complete system show self-oscillations of the vocal folds with constant sub-glottal total enthalpy. We proposed a regularisation procedure to eliminate the unrealistic behaviour of the discrete fluid model at glottal closure. Qualitative analysis of simulation results for different values of the regularisation parameter emphasises the effects of this procedure. Consequently, the proposed model is suitable for the sound synthesis of voiced sounds with coupling between glottis and vocal tract. It enables vowel articulation during self-oscillation and, theoretically, also consonant articulation, given the model’s capability to handle local closures of the vocal tract. The study of this last point will require the careful design of an appropriate control of the geometry.

The reduced-dimensional approach adopted in this work also proves valuable for exploring the parameter space, particularly in determining oscillation pressure thresholds and pressure-fundamental frequency relationships (Alipour and Scherer, 2007; Ruty et al., 2007). This offers a natural extension of the work presented, allowing comparison with in vivo measurements (in the context of the ANR project AVATARS) or other models. This ability to perform parameter studies is an advantage of the low order fluid models such as that proposed here. Models with two- or three-dimensional flows are more precise but generally too expensive (days of computation for a few oscillation cycles) to perform parameter studies. Voice synthesis based on machine learning is efficient and produces very convincing results. However, any access to the physical variables involved in voice production is lost.

To facilitate future numerical experiments, the development of a fast numerical integration method is of great interest. Indeed, the current in-house iterative solver based on the method presented in Müller (2021), chapter 3 is fully implicit and results in long simulation times (in the order of 10 min of computation time for 1 s of sound). Ongoing research is focused on the development of a fast (explicit or semi-explicit) time-stepping method that satisfies a discrete power balance. Energy quadratisation techniques involving auxiliary variables, as introduced in Shen et al. (2018), have been successfully applied or extended to various problems in musical acoustics (Ducceschi et al., 2023; Bilbao et al., 2023; Castera and Chabassier, 2023; Russo et al., 2024) and show promise for enabling real-time synthesis, which remains one of our primary future objectives.

4.2 Possible refinements for future models

The use of the port-Hamiltonian (pH) framework allows the procedural connection of multiple sub-systems in a power-balanced way. In this study, subsystems for fluid dynamics, tissues in the larynx and vocal tract, and acoustic radiation were considered. Each subsystem can be refined independently without breaking the overall energy balance and passivity, offering significant flexibility for future model improvements. We propose a list of perspectives in this regard.

4.2.1 Fluid model refinements

Use of the PDE presented in 2.1.1 is, to the best of our knowledge, new in the context of vocal apparatus modelling. Unlike Bernoulli-type glottal flow solvers commonly used in the literature, this model accounts for unsteady effects and compressibility, making it suitable for representing both wave propagation in the vocal tract and glottal flow dynamics. The spatial discretisation of this equation using finite differences results in the model presented in Mora et al. (2021b). However, alternative structure-preserving discretisation methods based on a variational formulation, such as partitioned finite elements (Cardoso-Ribeiro et al., 2018; Brugnoli et al., 2020), could also be employed. Specifically, the cross-section area of the fluid domain may be discretised using methods other than piecewise constant approximations. Our intuition is that a smooth laryngeal profile is expected to improve the behaviour of the discrete fluid model during channel closure and may alleviate the need for the regularisation procedure introduced in Section 2.4.2.

A first necessary step to apply these more involved discretisation methods is to include distributed losses directly in the continuous fluid model. In the larynx, viscous friction terms can be derived from a non-homogeneous axial velocity profile on cross-sections, modelled using a boundary layer approach. This is incompatible with the homogeneity assumption (H8) on , requiring a relaxation of this hypothesis in future work. A basic possible solution could be to approximate in the Hamiltonian (or to consider both the average and the standard deviation of as more detailed degrees of freedom). Additionally, kinetic energy losses due to partial pressure recovery at abrupt changes in cross-section area arise from complex turbulent phenomena. These phenomena cannot be captured well by a 3D-to-1D reduction but may be partially represented by localising an ideal point where airflow separates from the glottal walls. In this study, the separation point is fixed at the exit of the glottis. Future work could incorporate a dynamic separation point, as proposed in Pelorson et al. (1994).

4.2.2 Vocal fold model refinements

The vocal fold model we use is the body-cover model from Story and Titze (1994). Numerous lumped-element vocal fold models, often based on mass-spring assemblies, have been developed over the years (see Erath et al., 2013 for a review). Most of these models are straightforward to implement within the port-Hamiltonian framework and could be used with the fluid model presented here without major difficulty. The proposed framework also allows for the consideration of asymmetric vocal folds, which could be useful for studying pathological configurations (Xue et al., 2010).

Following the earlier discussion around the fluid equations, the use of “aerodynamically smooth” vocal fold models seems particularly promising. However, one of the main challenges of these lumped-element models is the identification of accurate physical parameters. Modal representations derived from the reduction of complex finite element models, as presented in Yokota et al. (2021), may help improve this.

4.2.3 Control issues

Controlling vocal fold parameters based on muscle activations is crucial for achieving realistic pitch modulation and natural vibrato mechanisms in sound synthesis. Integrating control rules, such as those proposed in Titze and Story (2002), into the current framework while maintaining the preservation of energetic properties also represents a promising avenue for future development.

References

• 1

  AlipourF.BerryD. A.TitzeI. R. (2000). A finite-element model of vocal-fold vibration. J. Acoust. Soc. Am.108, 3003–3012. 10.1121/1.1324678

  • CrossRef
  • Google Scholar
  • View reference in article

• 2

  AlipourF.SchererR. C. (2007). On pressure-frequency relations in the excised larynx. J. Acoust. Soc. Am.122, 2296–2305. 10.1121/1.2772230

  • CrossRef
  • Google Scholar
  • View reference in article

• 3

  BilbaoS.WebbC.WangZ.DucceschiM. (2023). “Real-time gong synthesis,” in Proceedings of the 26th Conference of Digital audio effects (DAFx-23) (Copenhagen, Denmark).

  • Google Scholar
  • View reference in article

• 4

  BirkholzP.HasnerP.KurbisS. (2022). “Acoustic comparison of physical vocal tract models with hard and soft walls,” in IEEE international conference on acoustics, speech and signal processing (ICASSP 2022) (IEEE), 8242–8246. 10.1109/ICASSP43922.2022.9746611

  • CrossRef
  • Google Scholar
  • View reference in article

• 5

  BrugnoliA.Cardoso-RibeiroF. L.HaineG.KotyczkaP. (2020). Partitioned finite element method for structured discretization with mixed boundary conditions. IFAC-PapersOnLine53, 7557–7562. 10.1016/j.ifacol.2020.12.1351

  • CrossRef
  • Google Scholar
  • View reference in article

• 6

  Cardoso-RibeiroF. L.MatignonD.LefèvreL. (2018). A structure-preserving partitioned finite element method for the 2D wave equation. IFAC-PapersOnLine51, 119–124. 10.1016/j.ifacol.2018.06.033

  • CrossRef
  • Google Scholar
  • View reference in article

• 7

  CasteraG.ChabassierJ. (2023). Numerical analysis of quadratized schemes. Application to the simulation of the nonlinear piano string. Tech. Rep.RR-9516. Inria.

  • Google Scholar
  • View reference in article

• 8

  DovalB.d’AlessandroC.Henrich BernardoniN. (2006). The spectrum of glottal flow models. Acta Acustica united Acustica92, 1026–1046.

  • Google Scholar
  • View reference in article

• 9

  DucceschiM.HamiltonM.RussoR. (2023). “Simulation of the snare-membrane collision in modal form using the scalar auxiliary variable (SAV) method,” in Proceedings of the forum acusticum 2023. 10.61782/fa.2023.0995

  • CrossRef
  • Google Scholar
  • View reference in article

• 10

  EncinaM.YuzJ.ZanartuM.GalindoG. (2015). “Vocal fold modeling through the port-Hamiltonian systems approach,” in 2015 IEEE conference on control applications (CCA) (Sydney, Australia: IEEE), 1558–1563. 10.1109/CCA.2015.7320832

  • CrossRef
  • Google Scholar
  • View reference in article

• 11

  ErathB. D.ZañartuM.StewartK. C.PlesniakM. W.SommerD. E.PetersonS. D. (2013). A review of lumped-element models of voiced speech. Speech Commun.55, 667–690. 10.1016/j.specom.2013.02.002

  • CrossRef
  • Google Scholar
  • View reference in article

• 12

  FlanaganJ. L. (1965). Speech analysis synthesis and perception. Berlin, Heidelberg: Springer Berlin Heidelberg. 10.1007/978-3-662-00849-2

  • CrossRef
  • Google Scholar
  • View reference in article

• 13

  FleischerM.PinkertS.MattheusW.MainkaA.MürbeD. (2015). Formant frequencies and bandwidths of the vocal tract transfer function are affected by the mechanical impedance of the vocal tract wall. Biomechanics Model. Mechanobiol.14, 719–733. 10.1007/s10237-014-0632-2

  • CrossRef
  • Google Scholar
  • View reference in article

• 14

  GuaschO.TernstromS.ArnelaM.AlıasF. (2013). Unified numerical simulation of the physics of voice:The EUNISON project. The International Society for Computers and Their Applications (ISCA).

  • Google Scholar
  • View reference in article

• 15

  GunterH. E. (2003). A mechanical model of vocal-fold collision with high spatial and temporal resolution. J. Acoust. Soc. Am.113, 994–1000. 10.1121/1.1534100

  • CrossRef
  • Google Scholar
  • View reference in article

• 16

  HélieT.SilvaF. (2017). “Self-oscillations of a vocal apparatus: a port-Hamiltonian formulation,” in Proceedings of geometric science of information: third international conference (GSI 2017). Editors Nielsen,F.BarbarescoF. (Paris, France: Springer International Publishing), 375–383.

  • Google Scholar
  • View reference in article

• 17

  IshizakaK.FlanaganJ. L. (1972). Synthesis of voiced sounds from a two-mass model of the vocal cords. Bell Syst. Tech. J.51, 1233–1268. 10.1002/j.1538-7305.1972.tb02651.x

  • CrossRef
  • Google Scholar
  • View reference in article

• 18

  KellyJ. L. (1962). “Speech synthesis,” in Proc. Speech communication seminar. Stockholm (Sep. 1962.

  • Google Scholar
  • View reference in article

• 19

  MaedaS. (1982). A digital simulation method of the vocal-tract system. Speech Commun.1, 199–229. 10.1016/0167-6393(82)90017-6

  • CrossRef
  • Google Scholar
  • View reference in article

• 20

  MaschkeB.van der SchaftA. (1992). Port-controlled Hamiltonian systems: modelling origins and systemtheoretic properties. IFAC Proc. Vol.25, 359–365. 10.1016/S1474-6670(17)52308-3

  • CrossRef
  • Google Scholar
  • View reference in article

• 21

  MoraL. A.Le GorrecY.MatignonD.RamirezH.YuzJ. I. (2021a). On port-Hamiltonian formulations of 3-dimensional compressible Newtonian fluids. Phys. Fluids33, 117117. 10.1063/5.0067784

  • CrossRef
  • Google Scholar
  • View reference in article

• 22

  MoraL. A.RamirezH.YuzJ. I.Le GorecY.ZañartuM. (2021b). Energy-based fluid–structure model of the vocal folds. IMA J. Math. Control Inf.38, 466–492. 10.1093/imamci/dnaa031

  • CrossRef
  • Google Scholar
  • View reference in article

• 23

  MoraL. A.YuzJ. I.RamirezH.GorrecY. L. (2018). A port-Hamiltonian fluid-structure interaction model for the vocal folds this work was supported by CONICYT-PFCHA/2017-21170472, and AC3E CONICYT-basal project FB-0008. IFAC-PapersOnLine51, 62–67. 10.1016/j.ifacol.2018.06.016

  • CrossRef
  • Google Scholar
  • View reference in article

• 24

  MüllerR. (2021). Time-continuous power-balanced simulation of nonlinear audio circuits: realtime processing framework and aliasing rejection. Sorbonne Université: Ph.D. thesis.

  • Google Scholar
  • View reference in article

• 25

  OlverP. J. (1986). Applications of lie groups to differential equations, 107. New York, NY: Graduate Texts in Mathematics. 10.1007/978-1-4684-0274-2

  • CrossRef
  • Google Scholar
  • View reference in article

• 26

  PelorsonX.HirschbergA.Van HasselR. R.WijnandsA. P. J.AureganY. (1994). Theoretical and experimental study of quasisteady-flow separation within the glottis during phonation. Application to a modified two-mass model. J. Acoust. Soc. Am.96, 3416–3431. 10.1121/1.411449

  • CrossRef
  • Google Scholar
  • View reference in article

• 27

  RisseT.HélieT.SilvaF.FalaizeA. (2024). Minimal port-Hamiltonian modeling of voice production: choices of fluid flow hypotheses, resulting structure and comparison. IFAC-PapersOnLine58, 238–243. 10.1016/j.ifacol.2024.08.287

  • CrossRef
  • Google Scholar
  • View reference in article

• 28

  RussoR.BilbaoS.DucceschiM. (2024). Scalar auxiliary variable techniques for nonlinear transverse string vibration. IFAC-PapersOnLine58, 160–165. 10.1016/j.ifacol.2024.08.274

  • CrossRef
  • Google Scholar
  • View reference in article

• 29

  RutyN.PelorsonX.Van HirtumA.Lopez-ArteagaI.HirschbergA. (2007). An in vitro setup to test the relevance and the accuracy of low-order vocal folds models. J. Acoust. Soc. Am.121, 479–490. 10.1121/1.2384846

  • CrossRef
  • Google Scholar
  • View reference in article

• 30

  ShenJ.XuJ.YangJ. (2018). The scalar auxiliary variable (SAV) approach for gradient flows. J. Comput. Phys.353, 407–416. 10.1016/j.jcp.2017.10.021

  • CrossRef
  • Google Scholar
  • View reference in article

• 31

  StevensK. N. (1971). Airflow and turbulence noise for fricative and stop consonants: static considerations. J. Acoust. Soc. Am.50, 1180–1192. 10.1121/1.1912751

  • CrossRef
  • Google Scholar
  • View reference in article

• 32

  StoryB. H.TitzeI. R. (1994). Voice simulation with a body-cover model of the vocal folds. J. Acoust. Soc. Am.97, 1249–1260. 10.1121/1.412234

  • CrossRef
  • Google Scholar
  • View reference in article

• 33

  StoryB. H.TitzeI. R.HoffmanE. A. (1996). Vocal tract area functions from magnetic resonance imaging. J. Acoust. Soc. Am.100, 537–554. 10.1121/1.415960

  • CrossRef
  • Google Scholar
  • View reference in article

• 34

  TitzeI. R.StoryB. H. (2002). Rules for controlling low-dimensional vocal fold models with muscle activation. J. Acoust. Soc. Am.112, 1064–1076. 10.1121/1.1496080

  • CrossRef
  • Google Scholar
  • View reference in article

• 35

  TrenchantV.RamirezH.Le GorrecY.KotyczkaP. (2018). Finite differences on staggered grids preserving the port-Hamiltonian structure with application to an acoustic duct. J. Comput. Phys.373, 673–697. 10.1016/j.jcp.2018.06.051

  • CrossRef
  • Google Scholar
  • View reference in article

• 36

  ValášekJ. (2021). Numerical simulation of fluid-structure-acoustic interaction in human phonation. Czech Technical University. Ph.D. thesis.

  • Google Scholar
  • View reference in article

• 37

  WetzelV. (2021). Lumped power-balanced modelling and simulation of the vocal apparatus: a fluid-structure interaction approach. Sorbonne Université: Ph.D. thesis.

  • Google Scholar
  • View reference in article

• 38

  XueQ.MittalR.ZhengX.BielamowiczS. (2010). A computational study of the effect of vocal-fold asymmetry on phonation. J. Acoust. Soc. Am.128, 818–827. 10.1121/1.3458839

  • CrossRef
  • Google Scholar
  • View reference in article

• 39

  YokotaK.IshikawaS.TakezakiK.KobaY.KijimotoS. (2021). Numerical analysis and physical consideration of vocal fold vibration by modal analysis. J. Sound Vib.514, 116442. 10.1016/j.jsv.2021.116442

  • CrossRef
  • Google Scholar
  • View reference in article

• 40

  ZhangZ. (2015). Regulation of glottal closure and airflow in a three-dimensional phonation model: implications for vocal intensity control. J. Acoust. Soc. Am.137, 898–910. 10.1121/1.4906272

  • CrossRef
  • Google Scholar
  • View reference in article