Under the given assumptions the system of voice production can be considered to be equivalent to the model in Fig. 43.
The solution of the system in Fig. 43 describes the volume velocities Ug,U1..Un for each segments in terms of differential equations:
(39)
These equations must be solved to determine Ug, which in turn is necessary to solve the motion equation for both masses of the vocal folds. The force acting on the surfaces of both parts of the vocal fold depends on the pressure differences along the exposed surface of the fold. The exposed areas are lgd1 and lgd2 respectively for m1 and m2 and the pressure values follow from eq. (33) (Ishizaka & Flanagan, 1972:1244):
(40)
However, the pressure values depend also on the actual displacements x1, x2 of the folds, which can be presented by the following table (Ishizaka & Flanagan, 1972:1244)
(41)
where x1min=-(Ag01/2lg) and x2min=-(Ag02/2lg) are the initial displacements of the folds.
These considerations lead to the formulation of the motion equations for both masses (Ishizaka & Flanagan, 1972:1244):
(42)
The forces F1 and F2 are defined in eq. (40) and (41).
The equations of motion (42) are linearized, in the sense that the coupling between the masses depends only on the difference between the lateral displacements of both parts of the vocal fold.
The equations (36)-(39) form the set of differential equations needed to describe the motion of both masses as well as the volume velocity of the air Ug. Additionally, the model provides the sizes of the glottal areas at the level of the first and second mass: Ag1 and Ag2. The equations have the form of ordinary differential equations, however, the motion equations involve the second order derivative of the folds' displacements. The numerical solution of the equations can be achieved both by an approximation using difference equations and by the integration of the equations using numerical methods (Press et al., 1986:547-577).
The difference equation can be used to approximate the relations (36)-(39) in the following manner:
(43)
where T denotes the sampling interval. When this discrete approximation is applied to the equations, it leads to the set of discrete difference equations, which can be solved iteratively for a given time resolution depending on the sampling interval T. The difference equations are presented by Ishizaka and Flanagan (1972:196-197).
T should be chosen carefully. Its limitations are due to the duration of the sound propagation through the vocal tract segments (T should be considerably shorter than the time sound takes to travel through the shortest tube segment, which constitutes T's lower boundary) and the stability of the numerical solution (upper boundary). Ishizaka and Flanagan (1972:1248) used sampling rates between 10 and 30 kHz.
The modelling of the vocal folds motion is computed iteratively. The current positions of the masses, velocities, forces and pressures are calculated based on the values found in the precedings steps. The initial position is at rest, i.e. the glottal areas are Ag01, Ag02 (x1=x2=0) and no flow occurs. The first loop is used to calculate the first samples of the glottal airflow, the forces and the displacements x1 and x2. The displacements determine the new values of the glottal areas Ag1, Ag2, which in turn are used in the computation of new samples of the glottal airflow. There is continued iteration until the desired number of solution samples is obtained.
The integration solution of the system (36)-(39) demands their rearrangement in order to form the set of the first order differential equation only. This can be achieved using the following substitution:
(44) 
(45) 
where z is a new variable. After this substitution the equations should be modified to have the following general form:
The method for solving the first order differential equations depends mainly on the boundary conditions, which are usually divided into two categories (Press et al., 1986:548):
The two-mass model of vocal
fold vibration involves the first of the two categories. There
are numerous methods of solving this type of equation systems. They are based
on Euler's method (
with
tn+1=tn+h), but are numerically more stable than Euler's
original method. Today, the Runge-Kutta method is usually used (Press et
al., 1986:550-554). It is not the fastest, but one of the safest methods
(the procedure converges almost always to a solution). In the Runge-Kutta
method the order of the error estimation is used as an additional parameter.
Usually the Runge-Kutta method of the fourth or fifth order is implemented
in computer programs. The procedure is available in public mathematical libraries
(LINPACK, EISPACK) or in general purpose programs (MATLAB, Mathematica).
The main advantage of the integration solution compared to the direct discretization is the adaptive control of step size in the integration method. This allows a quicker computation of the equation, however not for all samples. In the present study, both methods were used and the results were almost equal, however it was sometimes more convenient to use direct discretization to obtain the normal sampled signal. The integration solution was implemented within the SIMULINK Toolbox of Matlab environment, the discrete approximation was coded in Pascal.
|
