Authors: Wayne Zhao, Rita Singh
Categories: Article, dynamical systems, parameter estimation, phonation models, vocal fold oscillation, voice profiling
Source: Entropy
Doi: 10.3390/e25071039
During phonation, the vocal folds exhibit a self-sustained oscillatory motion, which is influenced by the physical properties of the speaker’s vocal folds and driven by the balance of bio-mechanical and aerodynamic forces across the glottis. Subtle changes in the speaker’s physical state can affect voice production and alter these oscillatory patterns. Measuring these can be valuable in developing computational tools that analyze voice to infer the speaker’s state. Traditionally, vocal fold oscillations (VFOs) are measured directly using physical devices in clinical settings. In this paper, we propose a novel analysis-by-synthesis approach that allows us to infer the VFOs directly from recorded speech signals on an individualized, speaker-by-speaker basis. The approach, called the ADLES-VFT algorithm, is proposed in the context of a joint model that combines a phonation model (with a glottal flow waveform as the output) and a vocal tract acoustic wave propagation model such that the output of the joint model is an estimated waveform. The ADLES-VFT algorithm is a forward-backward algorithm which minimizes the error between the recorded waveform and the output of this joint model to estimate its parameters. Once estimated, these parameter values are used in conjunction with a phonation model to obtain its solutions. Since the parameters correlate with the physical properties of the vocal folds of the speaker, model solutions obtained using them represent the individualized VFOs for each speaker. The approach is flexible and can be applied to various phonation models. In addition to presenting the methodology, we show how the VFOs can be quantified from a dynamical systems perspective for classification purposes. Mathematical derivations are provided in an appendix for better readability.
Keywords: vocal fold oscillation, phonation models, dynamical systems, parameter estimation, voice profiling
Phonation is a complex bio-mechanical process wherein the glottal airflow, mediated by the muscles in the larynx and driven by an intricate balance of aerodynamic and mechanical forces across the glottis, maintains the vocal folds in a state of self-sustained vibration [1,2]. During this process, depending on the physical state of the vocal folds, their eigenmodes of vibration synchronize, or strive to do so. This self-sustained motion of the vocal folds during phonation is highly sensitive to perturbations caused by many possible influencing factors, which may affect the speaker during speech production. In recent years, there has been a surge of interest in building voice-based diagnostic aids—computational models based on artificial intelligence and machine learning that can infer the speaker’s state (and thereby the factors that are affecting the speaker) from voice. Such applications can benefit greatly from being able to deduce the fine-level nuances in the motion of the vocal folds, and being able to measure the response to various perturbing factors to infer their nature. However, doing this using traditional methods on an individual basis for each speaker is very difficult. Traditional methods to observe and record vocal fold oscillations (VFOs) are based on actual physical measurements taken using various instruments in clinical settings, e.g., [3,4,5].
The primary focus of this paper is to derive the VFOs for a speaker directly from recorded voice signals, alleviating the need for taking physical measurements. The solution we propose is an analysis-by-synthesis approach, based on physical models of phonation. We propose a methodology to deduce the parameters of a chosen phonation model from a speaker’s voice recording, which can then be substituted into a VFO model to obtain speaker-specific solutions, which represent the vocal fold oscillations of the speaker.
In the paragraphs below, we first review some relevant facts about the process of phonation. In the sections that follow, we present our proposed methodology in two in the first, we show how we can infer the parameters of a physics-based model of phonation from measurements of glottal excitation. Since the parameters of such models represent the physical properties of the speaker’s vocal folds, using them to obtain solutions to the models gives us an estimate of the speaker’s VFOs. We subsequently show how to extend the model to include the physics of mucosal wave propagation in the vocal tract, and propose a forward-backward algorithm to estimate the parameters of the joint model. These can then be used in the corresponding phonation model to obtain its solutions.
Later on, we also explain how the solutions of these models (which are the deduced VFOs) can be characterized from a dynamical systems perspective to derive discriminative information in the form of features that are useful for classification tasks.
In this context, it is important to note at the outset that the models we use to demonstrate our methodology are simple and well-established models in the literature. It is not the goal of this paper to propose new models of phonation, but rather to propose a novel methodology to derive their parameters from speech signals, so that they can be solved to yield VFOs on a speaker-by-speaker basis. The main contribution of this paper lies in the derivation of the parameters of these models, and their individualized VFO solutions. The viability of these solutions is demonstrated experimentally using classification experiments.
By the myoelastic-aerodynamic theory of phonation, the forces in the laryngeal region that initiate and maintain phonation relate to (a) pressure balances and airflow dynamics within the supra-glottal and sub-glottal regions and (b) muscular control within the glottis and the larynx. The balance of forces necessary to cause self-sustained vibrations during phonation is created by physical phenomena such as the Bernoulli effect and the Coandǎ effect [6,7,8]. Figure 1 illustrates the interaction between these effects that is largely thought to drive the oscillations of the vocal folds.
Figure 1 Schematic of the balance of some key forces through one cycle of the self-sustained vibrations of the vocal folds. Sequential “snapshots” of the cycle are numbered 1–5. Arrows depict the direction and type of forces in the glottal area. The color codes for the arrows depict net forces due to the Pink–muscular; Green–Bernoulli effect; Yellow–Coandǎ effect; Blue–vocal fold elasticity and other factors; Black and Red–air pressure. Lighter shades of each color depict weaker forces. Figure from [9] with permission.
The process of phonation begins with the closing of the glottis. This closure is voluntary and facilitated by the laryngeal muscles. Once closed, the muscles do not actively play a role in sustaining the vibrations. Glottal closure is followed by a contraction of the lungs which pushes out air and causes an increase in pressure just below the glottis. When this subglottal pressure crosses a threshold, the vocal folds are pushed apart, and air rushes out of the narrow glottal opening into the much wider supra-glottal region, creating negative intra-glottal pressure (with reference to atmospheric air pressure) [9].
The exact physics of the airflow through the glottis during phonation is well studied, e.g., [2,10,11,12,13,14]. The current understanding from these, from the airflow perspective, is that the glottis forms a flow separation plane. The air expansion in this region and the low pressure created in the vicinity of the glottis through the Coandǎ effect-induced entrainment cause a lowering of pressure close to the glottis and a net downward force on the glottis. At the same time, lowered pressure in the glottal region due to the Bernoulli effect that ensues from the high-velocity air volume flow through the glottis exerts a negative force on the glottis. The negative Bernoulli pressure causes elastic recoil, causing it to begin to close again. The closing reduces the volume flow through the glottis, diminishing the downward forces acting on it. Increased pressure build-up in the sub-glottal region causes the glottis to open again. This chain of oscillations continues in a self-sustained fashion throughout phonation until voluntary muscle control intervenes to alter or stop it or as the respiratory volume of air in the lungs is exhausted.
Physical models of phonation, e.g., [9,11,15,16,17,18,19,20], attempt to explain this complex physical process using relations derived from actual physics, especially aerodynamics and the physics of mechanical structures.
For modeling purposes, we note that phonation is not the only source of excitation of the the vocal tract in producing speechsounds, which comprise both voiced and unvoiced sounds. However, phonation is indeed the primary source of excitation of the vocal tract in the production of voiced sounds, wherein the oscillation of the vocal folds modulates the pressure of the airflow to produce a (quasi-) periodic glottal flow wave at a fundamental frequency (the pitch), which in turn results in the occurrence of higher order harmonics. The resultant glottal flow further excites the vocal tract, which comprises the laryngeal cavity, the pharynx, and the oral and nasal cavities, to produce individual sounds. The vocal tract serves as a resonance chamber that produces formants. The identities of the different sounds produced within it are derived from these resonances, which in turn are largely dependent on the configurations of the vocal tract specified by their time-varying cross-sectional area.
From this perspective, phonation modeling has typically involved the modeling of two sub-processes: the self-sustained vibration of the vocal folds, and the propagation of the resultant pressure wave through the vocal tract [21]. Each sub-process model has associated parameters that determine the model output, given an input.
Following the above division of the process, we identify the two following model types, each modeling one of the sub-processes: (i) Vocal fold oscillation (VFO) models (also called vocal folds models, or oscillation models), and (ii) vocal tract (VT) models.
The VFO models describe the vibration of vocal folds and their aerodynamic interaction with airflow. Such models are of four broad one-mass models, e.g., [2,16,22,23,24], two-mass models, e.g., [11,15], multi-mass models [18], and finite element models, e.g., [17]. Each of these has proven to be useful in different contexts. The VT models describe the interaction of the glottal pressure wave with the vocal tract, which is turn has been described in the literature by varied models, such as statistical models, e.g., [25], geometric models, e.g., [26] and biomechanical models, e.g., [27].
In addition, different models are also applied to describe the aero-acoustic interaction of the glottal airflow and the vocal tract. Some examples of these are reflection-type line analog models and transmission line circuit analog models, e.g., [28], hybrid time-frequency domain models, e.g., [29] and finite-element models, e.g., [30].
Each of the models mentioned in Section 1.2 includes a set of parameters that determine its state, and output, given an input. For instance, given the parameters for a VFO model, the glottal flow waveform can be determined; given the glottal flow waveform as an input, and the parameters for a VFO model, the acoustic signal can be determined.
The problem we tackle in this paper is the inverse given a model output, we must estimate the model parameters from it.
This can be of great practical use. For example, with speaker-specific parameter setting, the output of these models can be used as a proxy for the actual vocal fold motion of the corresponding speaker. To obtain the parameters of such models, the traditional solution has been to take actual physical measurements of the vocal fold oscillations, or of the glottal flow using techniques such as high-speed videostroboscopy, as mentioned in Section 1. This is not always feasible.
On the other hand, the inverse problem of parameter estimation of phonation models is quite difficult to solve through purely computational means. For example, in order to estimate the parameters of a vocal tract (VT) model, one must take into account the vocal tract coupling, the effect of the lossy medium that comprises the walls of the vocal tract, lip radiation, etc. Without serious approximations, the inverse problem in this case becomes eventually intractable. Some approaches simplify the solution by discretizing the vocal tract as a sequence of consecutive tubes of varying cross-sectional area, or with a mesh-grid. However, these approximations invariably increase the estimation error.
This paper proposes a methodology for solving the inverse problem of phonation models through purely computational means. As mentioned earlier, the methodology follows an analysis-by-synthesis approach. We explain this by first reviewing our previously proposed Adjoint Least Squares Estimation (ADLES) algorithm [31] that estimates the parameters of a VFO model by minimizing the error between a reference glottal flow waveform and the signal generated by the physical VFO model. We then describe our proposed ADLES-VFT algorithm to estimate the parameters of a joint VFO and VT model (also called a body-cover model). Instead of comparing the model-generated excitation at the glottis to a reference glottal flow, the flow is propagated through the vocal tract to generate a signal at the lips, which is compared to a recorded voice signal which is used as a reference. The algorithm proposed iteratively re-estimates the model parameters by minimizing the error between the reference voice sample and this generated signal.
Once estimated, these parameters are used with the VFO model to generate the speaker’s VFOs.
In this section we describe the VFO and VT models that we use as examples in this paper. We also explain the formulation of the joint model that we ultimately use to solve the inverse problem of parameter estimation.
A schematic illustration of a general mass-spring oscillator model for the vocal folds is shown in Figure 2. This is used to model the phonation process as described below.
Figure 2 Approximating the vocal folds with mass-spring oscillators in the phonation process. Airflow from the lungs, driven by the subglottal pressure Ps, passes through the glottis, and vocal folds are set into a state of self-sustained vibration, producing the glottal flow ug which is a quasi-periodic pressure wave. The vibration of vocal folds is analogous to a pair of mass-spring-damper oscillators. Further, the glottal flow resonates in the speaker’s vocal tract and nasal tract and produces voiced sound.
One-mass models describe the vibration of the vocal folds as that of a single mass-damper-spring oscillator.
where x is lateral displacement of a mass M, B and K are damping and stiffness coefficients, respectively, f is the driving force, and t is time [2]. The driving force is velocity-dependent and can be estimated by Bernoulli’s energy
where Pg is the mean glottal pressure, Ps is sub-glottal pressure, ρ is air density, and v is the air particle velocity. The kinetic pressure in the supra-glottal region is neglected [2].
Other models, namely two-mass, multi-mass and finite element models can also be used as the basis for the VFO model, and are described briefly in the Appendix A for reference.
For our paper, we adopt the version of the one-mass model of the vocal folds proposed in [24], illustrated in Figure 3. This is an asymmetric body-cover model which models the left and right vocal folds individually as one-mass components of a coupled dynamical system. It incorporates an asymmetry parameter, which can emulate the asymmetry in the vibratory motions of left and right vocal folds, and hence is also ideally suited to modeling pathological or atypical phonation [32].
Figure 3 Diagram of the one-mass body-cover model for vocal folds. The lateral displacements at the midpoint of the left and right vocal folds are denoted as ξl and ξr, and ξ0 represents the half glottal width at rest.
The key assumptions made in formulating this model
This assumption allows us to represent the mucosal wave as a one-dimensional surface wave on the mucosal surface (the cover) and treat the remainder of the vocal folds (the body) as a single mass or safely neglect it. As a result, the oscillation model can be linearized, and the oscillatory conditions are much simplified while maintaining the model’s accuracy.
In the one-mass asymmetric model proposed in [24], with reference to Figure 3, the center-line of the glottis is denoted as the z-axis. At the midpoint (z=0) of the thickness of the vocal folds, the left and right vocal folds oscillate with lateral displacement ξl and ξr, resulting in a pair of coupled Van der Pol
where β is the coefficient incorporating mass, spring and damping coefficients, α is the glottal pressure coupling coefficient, and Δ is the asymmetry coefficient. For a male adult with normal voice, the reference values for the model parameters (from clinical measurements) are usually approximately set to α=0.5, β=0.32 and Δ=0.
The literature describes a number of different approaches to modeling the vocal-tract, including bio-mechanical models, statistical models, and geometric models. For reference, they are described briefly in the Appendix B.
In our work we use an acoustic wave propagation model described by PDEs for the vocal tract. The vocal tract itself is represented as tube of length L, beginning at the glottis and ending at the lips. Representing the distance along the central axis of the vocal tract as x, x∈(0,L) (where x=0 at the glottis and x=L at the lips) and the time-varying volume velocity of air at any position x along the vocal tract as u(x,t), it can be shown that the PDE that describes u(x,t) is given by
where f(x,t) represents the vocal tract profile, which models the characteristics of the vocal tract, including the effect of the nonuniform yielding wall on the acoustic flow dynamics, the effect of vocal tract coupling, lip radiation, etc., and must also be estimated by our algorithm. The derivation of Equation (4) is given in Appendix B. Note that if the vocal tract is assumed to be a rigid tube f(x,t)=0 and Equation (4) reduces to the well-known Webster-Horn equation [33]. In deriving Equation (4) we have assumed a static vocal tract, i.e., that the cross-sectional area A(x) of the vocal tract at any position x is constant. This assumption is valid during phonation, in particular during the steady state of sustained phonation; however our solution can also be extended to consider time-varying vocal tracts A(x,t), although we have not done so in this paper.
The oscillation of the vocal folds results in the movement of air with a time-varying volume velocity u0(t)=u(0,t) at the glottis. The vocal tract modulates this to result in the volume velocity uL(t)=u(L,t) at the lips and the corresponding pressure wave pL(t), where [34]
where A(L) is the opening area at the lip, c is the speed of sound, and ρ is the ambient air density.
Our objective is to derive vocal fold oscillations during phonation directly from the speech signal. In order to do so, we will utilize the VFO and VT models.
The VFO model represents the actual dynamics of the vocal folds. Given the model parameters, which are α, β and Δ for the coupled one-mass model of Equation (3), it can be used to compute the vocal fold oscillations and the volume velocity of air at the glottis. The VT model represents the dynamics of the vocal tract. Given an excitation (at the glottis), it can be used to compute the pressure wave at the lips, which manifests as sound.
Ours is the inverse given the the pressure wave at the output of the vocal tract (i.e., the recorded speech signal) we must estimate the VFO parameters that could generate it. This is an analysis-by-synthesis in order to analyze the given voice signal, we must identify the model parameters that synthesize the closest (in a metric sense) facsimile to it.
We present the solution in a two-step manner. In the first, given the actual excitation of the vocal tract to produce the voice signal, the parameters of the VFO model are estimated to minimize the error between its output and the reference excitation signal. We refer to this as the backward approach (and the corresponding estimation as the “backward” problem), since the reference excitation signal itself must first be derived by passing the voice signal through an inverse model of the vocal tract, i.e., “backward” through the vocal tract. We have previously described our solution to the backward problem in [31], and restate it in Section 4.1 for completeness.
The second, more complete solution considers the joint model, i.e., both the motions of the vocal folds and the propagation of the resulting air flow (the excitation) through the vocal tract. The model parameters are estimated by comparing the signal produced at the lips by the joint model to the recorded voice signal. We refer to this as the “forward-backward” approach since this requires forward propagating the output of the VFO through the VT model, prior to applying the backward approach. The solution to this problem is the primary contribution of this paper.
The two approaches are illustrated in Figure 4.
Figure 4 The VFO model models the generation of glottal signals by the movements of the vocal folds. The VT model models the transformation of the glottal signal generated by the vocal folds to the final voice signal. The joint VFO-VT model combines the two, using the output of the VFO model as the input to the VT model. ADLES compares the glottal signal u0(t) generated by the VFO model to a reference glottal signal ug(t) to estimate VFO parameters. ADLES-VFT compares the output of the joint model, uL(t), to a reference signal um(t) obtained from an actual voice recording, to estimate both VFO and VT parameters. The output of the VFO model is the desired vocal-fold oscillation.
Once estimated, the VFO model parameters can be used to compute the oscillations of the vocal folds and their phase-space trajectories.
As the first step, we describe the backward how to derive the VFO model parameters that best explain the glottal excitation for a given phonated signal. We use the approach proposed in [31]—we estimate the VFO model parameters to minimize the error between the volume velocity of air predicted by the model and a reference signal representing the actual glottal excitation to the vocal tract. If the VFO model were to be considered in isolation, this reference could be obtained through actual physical measurements, e.g., through photography [5], physical or numerical simulations [7,17], or by inverse filtering the speech signal using a technique such as [35] (the approach used in [31]).
For the purpose of our discussion in this section, however, we do not specify where this reference excitation is obtained from, since the estimation of VFO model parameters from a given glottal excitation is only a waypoint towards estimation from the joint model that includes both the VFO and VT components. As we will see in Section 4.2, this does not in fact require explicit knowledge of the reference signal at the glottis.
Let ug(t) be the reference signal representing the true air-volume velocity at the glottis that excites the vocal tract. The volume velocity of air u0(t) (we remind the reader that u0(t)=u(0,t)) at the glottis can also be computed from the vocal fold opening at the glottis as
where ξ0 is the half glottal width at rest, 2ξ0+ξl(t)+ξr(t) is the complete glottal opening, d is the length of the vocal fold, and c˜ is the air particle velocity at the midpoint of the vocal folds.
We assume that the movement of the vocal folds follows the VFO model of Section 3.1. Correspondingly, ξl(t) and ξr(t) must obey Equation (3), subject to boundary conditions. The model parameters α, β and Δ can hence be computed to minimize the difference between the air volume velocity u0(t) predicted by the model and the reference u0g(t).
We define the instantaneous residual R as the error between u0(t) and ug(t):
The overall ℓ2 error between u0(t) and ug(t) is given by the integral
where ϑ=[α,β,Δ] represents the parameters of the VFO model, and T represents the complete length of the reference signal.
The actual estimation can now be stated as
where Cr and Cl are constants representing the quiescent positions of the vocal folds, and the folds are assumed to be at rest prior to the onset of phonation. For the computation we set ξ0 to a typical value of 0.1 cm. The length of the vocal folds d may be set to 17.5 mm (which is within the range of normal lengths for both male and female subjects), and the air particle velocity c˜ to 5000 cm/s [2].
Note that given α,β and Δ the differential equations of model Equations (11)–(15) (the constraints) can be solved by any ODE solver to obtain ξl(t) and ξr(t). So, in principle, we could solve the constrained optimization problem of Equation (9)–(15) by a grid search over (α,β,Δ) to identify the specific values that minimize the squared error of Equation (8). This would, however, be highly inefficient.
Instead we propose the ADLES (“ADjoint LEast Squares”) algorithm, which restates the constraints (11)–(15) as Lagrangians on the objective, and derives a gradient descent solution. The detailed derivation of ADLES is given in Appendix C. We summarize the key steps below.
Incorporating constraints (11)–(15) into the objective, we define the Lagrangian:
where λ, η, μl,μr, νl and νr are Lagrangian multipliers. Note that λ and η are also functions of time (we have not explicitly shown the “(t)” above for brevity of notation).
We obtain the Lagrangian parameters λ and η as the solution to the following For 0<t<T:
with initial conditions at t=T:
Note that given α, β, Δ, ξr and ξl, Equations (17)–(24) represent a differential-algebraic system of equations and can be solved by any DAE solver to obtain λ and η.
Given ξr, ξl, λ and η, the derivatives of F(ξl,ξr,ϑ) w.r.t. α, β and Δ can now be obtained
The derivatives from Equation (25) are plugged into a gradient descent update rule for the model
where τ· is the step-size.
The overall ADLES algorithm is summarized in Algorithm
The backward approach, solved by the ADLES algorithm in Section 4.1, derives the VFO parameters by minimizing the error between the output of the VFO model u0(t) and the glottal excitation ug(t). However, in general, ug(t) is not available, and this error cannot actually be computed.
Instead, in the forward-backward approach, we further propagate the generated excitation u0(t) through the vocal tract, represented by the VT model of Equation (4) to obtain a signal uL(t)=u(L,t) at the lips. This is the output of the joint VFO and VT models. We compute the error between the generated signal uL(t) and the air velocity measurement derived from the recorded voice signal, which is available, and propagate this error backward through the vocal tract, to obtain the error at the glottis. The VT and VFO model parameters are estimated to minimize this error. Thus, the algorithm itself proceeds through iterations of two a forward step in which the VFO-generated excitation is propagated through the VT model to generate an output signal, and a backward step in which the error between the generated signal and the recorded speech is propagated backward through the VT to adjust the model parameters. We explain the entire procedure below.
The recorded voice signal is, in fact, a pressure wave and records the pressure wave emitted at the lips. Let pm(t) be the measured acoustic pressure at the lip. The corresponding volume velocity is given by [34]
um(t) is now our reference signal at the lips to which uL(t) must be matched, in order to estimate model parameters.
The propagation of u0(t)=u(0,t) through the vocal tract is assumed to follow the dynamics of Equation (4). Let Hf be the nonlinear operator representing acoustic wave propagation through the vocal tract from the glottis to the lip. The subscript f in Hf represents the vocal-tract profile f(x,t) in Equation (4) and indicates the dependence of Hf on f(x,t). Thus the vocal-tract output uL(t) is given by
Our objective is to minimize the difference between the measured volume velocity um(t) and the predicted volume velocity uL(t) near the lip subject to constraint that the dynamics of the vocal folds must follow the VFO model of Equation (3). Note that the parameters of the joint model include the VFO model parameters α, β and Δ, and the vocal tract profile f(x,t) required by the VT model. Although we only require the VFO model parameters to determine vocal fold oscillation, the minimization must be performed against all of these. Thus the estimation problem becomes
subject to
where, as before, (30) and (31) represent the asymmetric vocal folds displacement model (3), I.C. stands for initial condition, and Cs are constants. The minimization is performed against the complete set of parameters of the joint VFO-VT model, i.e., α, β, Δ and f(x,t).
Unlike in Equations (11)–(15), this cannot be solved, even in principle, by simply scanning for the optimal α, β and Δ, since Hf is characterized by f(x,t) which is also unknown and must be determined.
To solve the optimization problem of (29)–(35), we derive an efficient gradient-descent solution which we term the ADLES-VFT algorithm. The essential steps of the solution are given below. The details of the derivation are in Appendix D.
First, note that, as before, the constraint Equations (30)–(35) are ordinary differential equations with initial conditions that, given α, β and Δ, can be solved by any ODE solver. The solution will give us the VFO model generated glottal excitation u0(t).
Next, we propagate the generated excitation u0(t) through the VT model. For this, we must solve
where B.C. stands for boundary condition, and I.C. stands for initial condition. The vocal tract is assumed to be circular at the glottis and the lips. Here, nΓ is the outward unit normal to the vocal tract boundary Γ, at the glottis.
Equations (36)–(40) represent a set of partial differential equations. The boundary conditions relate to the air volume velocity at the glottal end of the vocal tract. The initial conditions relate to air volume velocity at the initial time, t=0, when the generated glottal signal u0(t) enters the vocal tract. Given u0(t) and f(x,t) (36)–(40) can be solved using a PDE solver. In our work we use the finite-element method described in Appendix E.
Solving (36)–(40) gives us u(x,t) at all positions x∈(0,L) and time t∈(0,T). In the process, it also gives us Hf(u0(t))=uL(t)=u(L,t).
The backward pass updates all model parameters including the VT term f(x,t), and VFO parameters based on the error at the output.
We denote the estimation residual
We must propagate this residual backward through the VT model. To do so, we use a time reversal technique [36] and backpropagate the difference (41) into the vocal tract, which
where z is the time reversal of u. Note that the boundary conditions and initial conditions in (42)–(45) are now defined at the lip, and the equation itself propagates backward through the vocal tract.
As before, Equations (42)–(45) can be solved by the finite-element method of Appendix E to give us z(x,t). The gradient update rule for f(x,t) is then obtained as
where ι is a learning rate parameter (see Appendix D).
As in the case of the backward approach of Section 4.1 we define a residual
Note that unlike in Section 4.1 the residual in Equation (47) is defined at the lips, rather than at the glottis. As before, we can define the total squared residual error as
where ϑ=[α,β,Δ] are the parameters of the vocal folds model (3). F must be minimized with respect to ϑ, subject to the constraints imposed by the VFO model.
Once again, as in Section 4.1 we fold in the constraints into the objective through Lagrangian multipliers as
where λ, η, μ and ν are multipliers, and λ and η are themselves functions of time. Optimization requires minimization of L(ϑ) with respect to α, β and Δ.
This leads us (Appendix D) to the following set of equations for λ and η:
with initial conditions (at t=T, i.e., at the lips):
Given R(t), u0(t) and uL(t) Equations (50)–(57) represent a set of differential-algebraic equations and can be solved with a DAE solver.
We finally obtain the derivative of F w.r.t. α, β and Δ (represented below as Fα, Fβ and FΔ, respectively) as
The gradient descent update rules for the VFO model parameters are finally obtained as
The overall ADLES-VFT algorithm for solving the parameter estimation problem (29)–(35) is summarized in Algorithm 2.
In this solution, we have adopted the simple gradient descent method. However, other gradient-based optimization approaches, such as the conjugate gradient method, can also be used.
We note here that Algorithm 2 requires several terms to be specified. In our implementation, the quiescent positions of the vocal folds, Cr and Cl were set to 0. We initialize [α,β,Δ]=[0.8,0.32,1.0]—these values were empirically found to work best. f(x,t) is initialized to 0. This effectively initializes Equation (4) to the Webster Horn equation. The step sizes τα,τβ and τΔ are all adaptively set to 0.01/max(Fα,Fβ,FΔ), and ι is set to 1. The actual objective minimize, Equation (29), requires scaling pm(t) by A(L)/ρc prior to comparison to uL(t). In practice, since pm(t) is derived from a voice signal recorded at a distance from the lips, the unknown transmission loss between the lips and the microphone must also be considered. To deal with this, we simply normalize both sequences to 0 mean and unit variance, and do not apply any additional scaling.
In the next section, we demonstrate the usefulness of the ADLES and ADLES-VFT algorithms experimentally.
Unfortunately, until the time of writing this paper, we could not obtain actual electroglottographic measurements or video data of vocal fold motion to compare our derived VFOs to. However, from a computational perspective the algorithms proposed can still be validated in different ways. We explain these below.
Our first validation approach is to use the proxy of showing that the solutions obtained are indeed discriminative of fine-level changes in glottal flow dynamics of the phonation process.
Having recovered the model parameters by our backward or forward approach, we can solve the models to obtain the time-series corresponding to the oscillations of each vocal fold, as estimated from recorded speech samples. We note that the models we have discussed in this paper are essentially dynamical systems represented by coupled nonlinear equations that may not have closed-form solutions, but can be numerically solved.
To interpret these in discriminative settings, we can utilize some well-established methods for characterizing dynamical systems, borrowing them from chaos theory and other areas of applied mathematics (e.g., flow, orbit, attractor, stability, Poincaré map, bifurcation, Lyapunov exponents, etc.). These are described in Appendix F.
Appendix F describes the concepts and tools used to study the behaviors (e.g., flow, orbit, attractor, stability, Poincaré map, bifurcation) of nonlinear dynamical systems such as Equation (3). The phase space of the system in Equation (3) (representing vocal fold motion) is four-dimensional and includes states (ξr,ξ˙r,ξl,ξ˙l).
For this nonlinear system, it is expected that attractors such as limit cycles or toruses will appear in the phase space. Such phenomena are consequences of specific parameter settings. Specifically, the parameter β determines the periodicity of oscillations; the parameter α and Δ quantify the asymmetry of the displacement of left and right vocal folds and the degree to which one of the vocal folds is out of phase with the other [24,37]. We can visualize them by plotting the left and right displacements and the phase space portrait.
The coupling of right and left oscillators is described by their entrainment; they are in n:m entrainment if their phases θr, θl satisfy |nθr−mθl|<C where n,m are integers, and C is a constant [24]. Such entrainment can be revealed by the Poincaré map, where the number of trajectory crossings of the right or left oscillators within the Poincaré section indicates the periodicity of its limit cycles. Therefore, their ratio represents the entrainment. We can use the bifurcation diagram to visualize how the entrainment changes with parameters. An example of such a bifurcation diagram is shown in Figure 5 [15,37].
Figure 5 (a) A 3D Bifurcation diagram of the asymmetric vocal fold model. The third dimension is perpendicular to the parameter plane shown, and depicts the entrainment ratio n:m (encoded in different shades of gray) as a function of model parameters α and Δ, where n and m are the number of intersections of the orbits of right and left oscillators across the Poincaré section ξ˙r,l=0 at stable status. This is consistent with the theoretical results in [24]); (b) Phase-space trajectories (or phase portraits) corresponding to the points A (left panel), B (center panel) and C (right panel). The horizontal axis is displacement of a vocal fold, and the vertical axis is its velocity.
As we will see later (and as indicated in Figure 5), model parameters can characterize voice pathologies, and these can be visually evident in in phase portraits and bifurcation plots.
We use the backward algorithm to estimate the asymmetric model parameters for clinically acquired pathological speech data. The data comprise speech samples collected from subjects suffering from three different vocal pathologies. Our goal is to demonstrate that the individualized phase space trajectories of the asymmetric vocal fold model are discriminative of these disorders.
The data used in our experiments is the FEMH database [38]. It comprises 200 recordings of the sustained vowel //. The data were obtained from a voice clinic in a tertiary teaching hospital, and the complete database includes 50 normal voice samples (control set) and 150 samples that represent common voice pathologies. Specifically, the set contains 40/60/50 samples for glottis neoplasm, phonotrauma (including vocal nodules, polyps, and cysts), and unilateral vocal paralysis, respectively.
Figure 6 shows some phase portraits showing the coupling of the right and left vocal folds obtained using the ADLES solution. We observe that the attractor behaviors are typical and even visually differentiable for different types of pathologies.
Figure 6 Phase portraits showing the coupling of the left and right oscillators (ADLES-based estimation) for (a) normal 1 limit cycle, (b) 1 limit cycle, (c) 2 limit cycles, (d) vocal limit torus. The convergence trajectory is also shown, and the limit cycles can be observed as the emergent geometries in these plots.
Table 1 shows the results of deducing voice pathologies by simple thresholding of parameter ranges. Specifically, the ranges of model parameters in each row of Table 1 correspond to regions in the bifurcation diagram in Figure 5. Each region has distinctive attractors and phase entrainment, representing distinct vocal fold behaviors and thereby indicating different voice pathologies. By extracting the phase trajectories for the speech signal and, thereby, the underlying system parameters, the ADLES algorithm can place the vocal-fold oscillations in voice production on the bifurcation diagram and thus deduce the pathology.
Our second validation approach is to compare the excitation signal obtained through inverse filtering with the glottal flow signal (VFO) obtained through the backward or forward-backward algorithm. The rationale behind this is that within reasonable bounds of error, the glottal flow signal obtained through our model is expected to conform to the oscillation patterns seen in the excitation signal for each speaker.
Figure 7 shows the glottal flow obtained by inverse filtering and those obtained by the asymmetric model with the parameters estimated by our ADLES method. We observe consistent matches, showing that the ADLES algorithm does achieve its objectives in individualizing the asymmetric model to each speaker instance.
Figure 7 Glottal flows from inverse filtering and ADLES estimation for (a) normal speech (control), (b) neoplasm, (c) phonotrauma, and (d) vocal palsy.
Our third validation approach is to compare the estimation precision of the backward approach and the forward-backward approach. Table 2 shows the mean absolute error (MAE) of calculating glottal flows and parameters for four voice types (normal, neoplasm, phonotrauma, vocal palsy) obtained by backward (ADLES) and forward-backward (ADLES-VFT) algorithms. The glottal flows obtained by inverse filtering the speech signals are treated as ground truths. Since there is no ground truth for model parameters, we treat the parameters obtained by ADLES as ground truth. These results suggest that our forward-backward algorithm can effectively recover the vocal tract profile, glottal flow, and model parameters.
Our fourth validation comes indirectly from prior studies. Information from a dynamical systems perspective can give insights about the underlying mechanisms and principles that govern the vocal fold dynamics. Examples of features in this category are recurrence analysis features, Lyapunov exponents, Hurst exponents, etc. These are mentioned in Appendix F.
Some of these features have been used in real-world applications and proven to be effective. For example, in [39], the authors hypothesize that since COVID-19 impairs the respiratory system, effects on the phonation process could be expected, and signatures of COVID-19 could manifest in the vibration patterns of the vocal folds. In this paper, features have been derived from a signal processing perspective.
This study used the ADLES method to estimate the asymmetric vocal folds model parameters. It further used the parameters and estimation residuals as features to other binary classifiers such as logistic regression, support vector machine, decision tree, and random forest, achieving around 0.8 ROC-AUC (area under the ROC curve) in discriminating positive COVID-19 cases from negative instances, on clinically collected and curated data. The data used contained recordings of extended vowel sounds from affected speakers and control subjects. The authors also discovered that COVID-19 positive individuals display different phase space behaviors from negative the phase space trajectories for negative individuals were found to be more regular and symmetric across the two vocal folds, while the trajectories for positive patients were more chaotic, implying a lack of synchronization and a higher degree of asymmetry in the vibrations of the left and right vocal folds.
In a companion study, the authors in [40] used the ADLES-estimated glottal flows as features to CNN-based two-step attention neural networks. The neural model detects differences in the estimated and actual glottal flows and predicts two classes corresponding to COVID-19 positive and negative cases. This achieved 0.9 ROC-AUC (normalized) on clinically collected vowel sounds. Yet another study used higher order statistics derived from parameters, and Lyapunov and Hurst exponents derived from the phase space trajectories of the individualized asymmetric models, to detect Amyotrophic Lateral Sclerosis (ALS) from voice with high accuracy (normalized ROC-AUC of 0.82 to 0.99) [41].
In this paper we have presented a dynamical system perspective for physical process modeling and phase space characterization of phonation, and proposed a framework wherein these can be derived for individual speakers from recorded speech samples. The oscillatory dynamics of vocal folds provide a tool to analyze different phonation phenomena in many real-world task settings. We have proposed a backward approach for modeling vocal fold dynamics, and an efficient algorithm (the ADLES algorithm) to solve the inverse problem of estimating model parameters from speech observations. Further, we have integrated the vocal tract and vocal folds models, and have presented a forward-backward paradigm (the ADLES-VFT algorithm) for effectively solving the inverse problem for the coupled vocal fold-tract model. Extensions of these approaches can use other physical models of voice production, and other physical processes including phonation.
We have shown that the parameters estimated by these algorithms allow the models to closely emulate the vocal fold motion of individual speakers. Features and statistics derived from the model dynamics are (at least) discriminative enough for use in regular machine-learning based classification algorithms to accurately identify various voice pathologies from recorded speech samples. In future, these approaches are expected to be helpful in deducing many other underlying influences on the speaker’s vocal production mechanism.
The phase space characterization presented in this paper is based on phase space trajectories (a topological perspective)—the left and right vocal fold oscillations, velocities or accelerations. Measurements can also be performed from statistical, signal processing, information-theoretic and other perspectives. Another direction of suggested research is characterizing the phase space from algebraic perspectives. We can recast the study of the topological structures of the phase space to the study of its algebraic constructs, such as homotopy groups and homology/cohomology groups, which are easier to classify. For example, algebraic invariants can characterize the homeomorphisms between phase spaces (e.g., evolution maps, Poincaré maps) and reveal large-scale structures and global properties (e.g., existence and structure of orbits). We can also build upon the deep connection between dynamical systems and deep neural models. We can study deep learning approaches for solving and analyzing dynamical systems, and explore the integration of dynamical systems with deep neural models to analyze and interpret the behaviors of the vocal folds. We delegate these explorations to future work.
The following is a brief description of various VFO
Two-mass models describe vocal fold motion as two coupled mass-damper-spring
where xi, Mi, Bi are the i-th oscillator’s displacement, mass, and viscous damping coefficient, K is the coupling stiffness between the two masses, Fi is the driving force, and Ri is the elastic restoring force [15]. This model assumes (1) small air inertia and quasi-steady glottal flow, (2) negligible supra-glottal pressure, and (3) that the nonlinearity induced by vocal fold collision is small. These assumptions lead to small-amplitude oscillations and model simplification.
Multi-mass models:
Finite element models:
The following is a brief description of various VT
Bio-mechanical models:
Statistical models:
Geometric models:
The vocal tract can be viewed as a compact, orientable, differentiable manifold M embedded in R3. Its boundary ∂M includes the wall of the vocal tract. Consider the tangent bundle TM. Denote the set of all vector fields on TM as Γ(TM), which is a C∞(M)-module [43]. A vector field is a smooth section on TM,Γ(TM)∋X:M→TM. It associates each point p∈M with a tangent vector v¯(p):=Xp:C∞(M)→∼R[43]. Let γ(t):R⊇I→M be a maximal integral curve [43] through p at t0, which is a solution
The curve γ(t) is a one-parameter group. When acting on the Lie group M, it gives the flow Φ:R×M→M.Φt(p)=γ(t). The particle velocity at p is given by v(p,t):=γ′(t)=v¯(p)∘γ(t)
The planar motion of the pressure wave in the vocal tract is governed by the equations [36]:
where p^(p,t) is the acoustic pressure, div is the divergence operator, grad is the gradient operator, ρ is the ambient air density, and c is the speed of sound. Equation (A11) describes the conservation of mass, and (A12) describes the conservation of momentum [36].
For notational convenience, we use cylindrical coordinates p=(r,θ,x), where the x direction aligns with the central axis of vocal tract. Denote the inner surface of vocal tract as Σ, and the shape function of inner surface as r=R(θ,x). Then the cross-sectional area of the vocal tract
The average acoustic pressure
and the volume velocity
where vx is the x component of v. Integrating (A11) over the volume of vocal tract bounded by cross sections at x0 and x
where we substitute Equations (A17) and (A18) into (A14) and (A15), and apply Stokes’ theorem [28,36]; nv is the component of v normal and outward to the inner surface Σ.
The element of area dσ is given by [28,36]:
where Sdθdx is a top 2-form on Σ [43]. Substituting (A19) into (A18) and differentiating w.r.t. x
Following similar steps, integrating the x component of (A12) over the cross section at x
where pw is the pressure acting on the wall of the vocal tract.
To simplify our problem, we combine the wave Equations (A20) and (A21) into a single vocal tract model. Differentiating (A20) w.r.t x and (A21) w.r.t. t, and canceling out the pressure term
where the vocal tract profile is absorbed into a single term f(x,t). This represents the characteristics of the vocal tract, i.e., the effect of the nonuniform yielding wall on the acoustic flow dynamics, which needs to be estimated by our algorithm.
In this appendix, we derive the ADLES algorithm to estimate VFO model parameters from a given glottal excitation signal.
Let ug(t) be the actual air volume velocity at the glottis. We can can also derive u0(t) from the displacement of the vocal folds as
where ξ0 is the half glottal width at rest, d is the length of vocal fold, and c˜ is the air particle velocity at the midpoint of the vocal fold.
We assume the movements of the vocal folds to follow the dynamics specified by Equation (3) in Section 3.1. Our objective is then to estimate the parameters of the model to minimize the
such that the movements of the vocal folds follow the VFO
where Cr and Cl are constants representing the quiescent positions of the vocal folds.
To solve the functional least squares in (A26), we require the gradients of (A26) w.r.t. the model parameters α, β and Δ. Subsequently, we can adopt any gradient-based (local or global) method to obtain the solution.
Considering the residual
We have
and
where ϑ=[α,β,Δ].
Incorporating the constraints into the objective using Lagrangian multipliers, we define the Lagrangian:
where λ, η, μ and ν are Lagrangian multipliers. Taking the derivative of the Lagrangian w.r.t. the model parameter α
Integrating the term λ∂αξ¨r twice by parts
Applying the same to η∂αξ¨l, substituting into (A37) and simplifying the final expression we
Since the partial derivative of the model output ξ w.r.t. the model parameter α is difficult to compute, we eliminate the related terms by setting
For 0<t<T:
with initial conditions at t=T:
As a result, we obtain the derivative of F w.r.t. α
The derivatives of F w.r.t. β and Δ are similarly obtained
Having calculated the gradients of F w.r.t. the model parameters, we can now apply the gradient-descent rule to optimize our objective (A26):
where τ· is the step-size.
The overall algorithm is summarized as
In this appendix we formulate the problem of estimating the parameters of the combined vocal fold-tract model from speech measurements. Let Ω×T be the domain of volume velocity u, where Ω is the spatial domain, and T is the time domain. In the one-dimensional case, Ω=[0,L] where L is the length of vocal tract, and T=[0,T]. Given measured acoustic pressure pm(t) at the lip, the corresponding volume velocity is given by [34]:
where A(L) is the opening area at the lip, c is the speed of sound, and ρ is the ambient air density. We denote u0(t)u(0,t), uL(t)u(L,t). The glottal flow u0(t) can be derived from the vocal folds displacement model (3)
where ξ0 is the half glottal width at rest, d is the length of the vocal fold, and c˜ is the air particle velocity at the midpoint of the vocal fold (see Figure 3).
Let H be the nonlinear operator representing acoustic wave propagation from the glottis to the lip. We have the forward propagation
where f is the vocal tract profile in (A22), and Γ=∂Ω is the boundary.
We can split Γ into two parts Γ=Γ0⋃Γ1, Γ0⋂ΓL=∅ corresponding to x=0 and x=L. However, we neglect the difference to ease our derivation. Note that in the one-dimensional case u(t) and uL(t) are only functions of t. However, more generally, they are functions of both x on the boundary Γ and t. We now define two nonlinear operators
and
Note that both Hf and F are bounded. Our objective is to (estimate model parameters α, β and Δ such that they) minimize the difference between the measured volume velocity um and the predicted volume velocity uL near the lip subject to
where (30) and (A60) represent the asymmetric vocal folds displacement model (3), I.C. stands for initial condition, and Cs are constants.
We solve the inverse problem for the combined vocal folds-tract model using the Forward-Backward method proposed below. We call this algorithm the Adjoint Least Squares—Vocal Folds-Tract (ADLES-VFT) parameter estimation algorithm.
In order to solve the parameter estimation problem represented by (A60)–(A64), first we need to estimate the vocal tract profile f in H and (A23). Specifically, we need to
where B.C. stands for boundary condition, ug is the volume velocity at the glottis and nΓ is the outward unit normal to the boundary Γ. We now derive the solution to (A65)–(A69).
In order to estimate f∈L2(Ω×T), we take an iterative approach, i.e.,
where δfk∈L2(Ω×T) is a small variation, and τ is a step size. A Taylor expansion of F (A56) at fk
where F′ is the Fréchet derivative [44]. Omitting higher order terms, we
where F′(f) is a nonlinear operator
Correspondingly, the adjoint operator [44,45,46]
We would like F(fk+δfk)=uLk+δuLk→k→∞um. This is equivalent to
It is simple to obtain the solution to (A75):
where F′(fk)* is the adjoint operator. It is difficult to compute F′(fk)F′(fk)*. We use its property of positive-definiteness to approximate it by γI where I is the identity matrix.
We denote the estimation residual
We have
Now consider the wave Equation (A65). Let u+δu be a solution with variation f+δf. Substitution into (A65)
Subtracting (A65)
Next, we use a time reversal technique [36] and backpropagate the difference (A76) into the vocal tract, which
where z is the time reversal of u. It follows [47]
where we substitute from (A87)–(A89) into (A79) and (A83); from (A89)–(A92) we apply initial conditions (A81), (A82), (A85) and (A86); from (A92)–(A94) we integrate by parts; from (A94)–(A95) we apply boundary condition (A80); from (A95)–(A96) we use boundary condition (A84); from (A96)–(A97) we use definition (A73); from (A97)–(A98) we use definition (A74) and the duality property
from (A98)–(A99), we assume that the second-order variation is small, i.e.,
(or δ(fu)=δ(f)u+fδ(u)≈0). By the arbitrariness of δf, it follows
and hence
Substitution into (A77) and (A70)
where ι=τ/γ.
Hence, we obtain an iterative forward-backward approach for solving the vocal tract profile f.
Now, we derive solution to the parameter estimation problem (A60)–(A64) using the adjoint least squares method of Appendix C as follows.
Denote the estimation error
and
where ϑ=[α,β,Δ] are the parameters of the vocal folds model (3). We would like to obtain update rules for the model parameters α, β, and Δ, i.e.,
where the the partial derivatives F·∂·F≡∂F∂· and τ· is the step size. We now define the Lagrangian:
where λ, η, μ and ν are multipliers. Taking the derivative of the Lagrangian w.r.t. the model parameter α yields
Integrating the term λ∂αξ¨r by parts twice
Defining the estimation residual R :=Hf(u0)−A(L)ρcpm(t), applying this to η∂αξ¨l, and substituting into (A112), after simplification
where the term Hf′|u0≈uL/u0 by linearization. Since the partial derivatives of the displacement ξ w.r.t. the model parameter α are difficult to compute, we cancel out the related terms by
For 0<t<T :
with initial
At t=T :
Assuming the boundary conditions on ξl and ξr to hold—an assumption we empirically find to be valid for Cr=Cl=0, we obtain the derivative of F w.r.t. α:
Similarly, we obtain the derivatives of F w.r.t. β and Δ
The algorithm for solving the parameter estimation problem (A60)–(A64) is outlined below.
In this appendix we describe a finite-element solution for time-dependent systems of PDEs examplified by the wave propagation
and the reverse propagation
where the solution must be defined for x∈(0,L) and t∈(0,T).
First, for the time-dependent system of PDEs, we discretize it along time t with the backward Euler method [48], yielding a sequence of differential equations. We split the time domain T into N uniform length intervals Δt. For time step n, 0≤n≤N−1, applying the backward Euler method to the left side of (A65)
where DtDt− is a finite difference operator w.r.t. time, and the superscript n indicates that the differencing is performed at time step n [48,49]. Substitution into (36)
Next, we define the residual at time step n
Applying Galerkin’s method [48,50]
where v∈Wk,2 (Wk,2 is the Sobolev space of functions with bounded L2 norm and k-th order weak derivatives) is a qualified test function. Galerkin’s method orthogonally projects the residual to the function space Wk,2. Expanding (A139)
Integration by parts for the second-order term in (A140)
where nΓ is the outward normal unit vector of the boundary Γ, and ds is the 1-form [43] on Γ. For problem (A126)–(A130), applying the boundary condition (A128) and substituting (A141) back into (A140) yields the variational
For the reverse problem represented by (A131)–(A134), applying the boundary condition (A132) and substituting (A141) back into (A140) yields a similar variational
We can split the variational problem (A142) into two
where we have interchanged the unknown un with u. (A144) is the bilinear form, and (A145) is the linear form [48].
Our original problem (A126)–(A130) and (A131)–(A134) then reduce to
for each time step. By the Lax-Milgram Lemma [49], solving (A146) is equivalent to solving the functional minimization
By the calculus of variations, and taking the variation of the functional gives (A146), hence the name variational form [48,49].
For each time step, we solve (A146) with finite element method. We discretize the domain Ω with a mesh of uniformly spaced triangular cells. We take the P2 elements as the basis function space, which contains piece-wise, second-order Lagrange polynomials defined over a cell. Each basis function has a degree-of-freedom (DoF) of 6 over a two-dimensional cell [48,51].
Each element is associated with a coordinate map that transforms local coordinates to global coordinates and a DoF map that maps local DoF to global DoF [48,51]. Each cell is essentially a simplex and can be continuously transformed into the physical domain. We use the following two previously established points in this context, and proceed as explained thereafter.
Assume a solution u=B+cjψj (using Einstein summation convention) with basis ψj∈P2 and coefficients cj. Here the ψj are piece-wise second-order Lagrange polynomials [48]. Also, the mesh of cells, P2 elements of basis functions, coordinate map, DoF map, and the assembly of the linear systems, are implemented and solved with the FEniCS software [52].
The function B(x) incorporates the boundary condition and, as an example, can take the
We also project B(x) over the basis functions P2 and express it as B(x)=bjψj. As a result, we obtain an unified expression u=Ujψj with Uj incorporating bj and cj. Similarly, we have fn=Fnjψj, un−1=Un−1jψj, un−2=Un−2jψj. Set the test function as v=ψ^i (where ψ^i is the ith Lagrange polynomial basis in a test cell). Substitution into (A144) and (A145)
where ψi′ is the derivative of ψi w.r.t. x.
Setting Mi,j=∫Ωψ^iψjdx, Ki,j=∫Ωψi′ψj′dx and collecting (A149) and (A150) into matrix-vector form, we
where A=M+Δt2c2K, and b=Δt2MFn+2MUn−1−MUn−2. Hence, we have reduced the problem of (A146) into solving the linear system (A151), with the solution described above. Furthermore, the matrices M (known as the mass matrix) and K (known as the stiffness matrix) can be pre-calculated for efficiency.
This appendix provides some essential definitions used for the characterization of dynamical system models for analysis.
A phonation model outputs a phase space trajectory of state variables that describes the movements of the vocal folds. The trajectories tend to fall into orbits with regular or irregular behaviors that explain observed behavior patterns of the vocal folds. The possible types and distributions of these orbits depend on the system parameters.
A dynamical system can be instantiated with ordinary or partial differential equations with initial conditions, and the evolution function Φ is the solution to the ODE or PDE.
Write Φx(t)≡Φt(x)≡Φ(t,x). The map Φt:M→M is a diffeomorphism (i.e., differentiable, invertible, bijection map between manifolds).
The behaviors of flows can be described by their attractor/attraction sets.
Namely, the orbit γx is “trapped” in the interior of A. A dynamical system can have more than one attractor set depending on the choice of initial points (or the choice of parameters, as we will see later). Locally we can talk about a basin of attractionB(A), which is a neighborhood of A satisfying for any initial point x∈B(A), and its orbit is eventually trapped in A.
There are different types of attractor sets. Some are shown in Figure A1. The simplest one is a fixed point or an equilibrium point, to which a trajectory in phase space converges regardless of initial settings of the variables (or their starting point). To study vocal fold behaviors, we are particularly interested in those attractors revealing the periodic motion of the flow in phase space. Such attractors include the limit cycle or the limit torus, which are isolated periodic or toroidal orbits, respectively. Some attractor sets have a fractal structure resultant from a chaotic state of the dynamical system [20,54], and are called strange attractors.
Figure A1 Illustration of different attractors in a dynamical system.
Chaos is a characteristic state of a nonlinear dynamical system. There are different definitions of chaos. A simple one is as
Thus, by definition, chaos is a state characterized by extreme sensitivity to initial conditions (trajectories starting from any two arbitrarily close initial conditions diverge exponentially). The Lypunov exponent is used to quantify this divergence. It also measures the measures the sensitivity of the evolution of the dynamical system to initial conditions.
Attractor sets (of all types) also characterize the stability of dynamical systems.
To study the trajectory (or orbit) structure of dynamical systems, we use the Poincaré map or Poincaré section.
Since the flow of a dynamical system in its phase space is a function of its parameters, the topological structure of the trajectories (including attractor sets) in phase space changes as the parameters change. To see how the topological structure changes with system parameters, we study the bifurcation map of the system.
At a bifurcation point, the system stability may change as the topological structure splits or merges, such as the periodic doubling or halving of a limit cycle.
Methodology, W.Z. and R.S.; Validation, W.Z.; Formal analysis, W.Z. and R.S.; Investigation, W.Z.; Writing—original draft, W.Z. and R.S.; Supervision, R.S. All authors have read and agreed to the published version of the manuscript.
Not applicable.
This study used the Far Eastern Memorial Hospital (FEMH) database, which is available to all participants through the annual IEEE FEMH Voice Data Challenges. A similar dataset called “VOICED Database” is openly available through PhysioNet https://physionet.org/content/voiced/1.0.0/ (accessed on 1 June 2023).
The authors declare no conflict of interest.
This research received no external funding.
This study used the Far Eastern Memorial Hospital (FEMH) database, which is available to all participants through the annual IEEE FEMH Voice Data Challenges. A similar dataset called “VOICED Database” is openly available through PhysioNet https://physionet.org/content/voiced/1.0.0/ (accessed on 1 June 2023).