Design of a Computational Procedure for Optimizing the Acoustic Outputs of the Human Vocal Tract

1. Introduction

The subject of interest of this work is explaining the impact of geometric configuration of the human vocal tract (hereinafter referred to as the VT) on the distribution of acoustic energy on the frequency spectrum of vocal speech. For voice teachers, it is necessary to analyse the procedures facilitating the impact of geometric configuration of the VT to amplify acoustic energy values in the required frequency range [1,2]. In order to shed light on such modifications of human VT on the characteristics of vocal speech, it is necessary to assemble mathematical calculation models enabling effective examination of sensitivity on individual geometric changes in the VT [3,4]. Assembling a fully parametric model at the level of describing boundary surfaces defining the examined acoustic system with regard to the complexity of the VT is a very demanding problem [5,6,7]. Currently, the FEM is almost exclusively used for analysing the acoustic characteristics of VT, the computation complexity of which increases with the number of finite elements used, which discredits the examined volume of the acoustic system. Direct geometric reconfiguration at the level of the finite element mesh represents a very complex problem, which requires a detailed analyses of the derived model [8]. However, due to the computation demands, such models are oftentimes not suitable for subsequent optimization tasks, whose aim it is to finetune the acoustic characteristics of phonation to the required parameters.

One of the possible approaches for assembling effective computation models is to use the super-elements method [9,10], which allows for modifying the analysed acoustic system only at the point of the geometric change. To speed up the calculation, the mesh of individual elements is reduced at the level of mass, damping and stiffness matrices to a number of degrees of freedom (hereinafter referred to only as “DoF”) ensuring compliance of the frequency characteristics of a full and reduced model in the frequency range, which is important for the comprehension of voice speech.

The article also presents an alternative, yet very effective approach to assembling an isoparametric element of a higher polynomial, using typical VT geometric configuration. Such an element even allows parallel acoustic cavities to be connected. Formulating the calculation model for [a:] was subsequently used to examine the permissible geometric modifications of VT, aimed at connecting the 3rd and 4th formant and shifting it into the 2.5kHz – 3.5kHz frequency range. Bringing the formants closer together makes it possible to significantly increase the amount of acoustic energy in the field of the so-called singing formant.

2. Finite Element Model of the Vocal Tract

The VT geometry is unique for every individual. A human’s ability to change the size and shape of the VT cross-section in different places, allows us to modulate the basic signal generated by the vocal cords and create different sounds. For the purpose of VT modelling, VT geometry [11] and [12] were measured. The CT LightSpeed VCT GE-64 was used for imaging. The measurement resulted in 181 sections, 0.625mm thick. The image resolution is 512x512 pixels. It was measured in the supine position, while sounding out [a:]. The geometric model was created by analysing sagittal CT cross-sections. The surface of the geometric model is used to create the FEM of the VT model (Figure 1). The geometric model contains the parallel cavities of the piriform sinuses and the epiglottic vallecula.

The formulated calculation model for [a:] was then used to examine the permissible geometric modifications of the VT, with the aim of connecting the 3rd and 4th formants and shifting them into the 2.5kHz – 3.5kHz frequency range. Bringing these formants closer together makes it possible to significantly increase the amount of acoustic energy in the field of the so-called singing formant.

The foundation of FEM is the discretization of the solved domain with the finite number of elements, which disjointly fill the entire solved volume of VT. The mathematical description of the elements is provided by the ${\mathit{N}}_{\mathit{e}}[-]$, shape functions, which distribute the nodal values of the unknown quantities into the element’s entire volume. Multiplying by the column vector of the nodal values of the sound pressure of element ${\mathit{P}}_e\left[Pa\right]$ results in the following acoustic pressure distribution function in the volume of the element

$$p^{\prime }\left(t,x,y,z\right)={\mathit{N}}_{\mathit{e}}\left(x,y,z\right){\mathit{P}}_e\left(t\right).$$

(1)

The finite element method uses a so-called weak formulation and the Galerkin Method [13]. The differential equation of the strong formulation is a scalarly multiplied by testing ${\delta p}^{\prime }$ function and the resulting expression is integrated via the solved volume

$$\underset{V_e}{\overset{}{\int }}\left({\delta p}^{\prime }\frac{1}{c_0^2}\frac{{{\partial }^{2}p}^{\prime }}{{\partial t}^2}+{\delta p}^{\prime }\frac{r}{c_0^2{\rho }_0}\frac{\partial p^{\prime }}{\partial t}\right)dV=\underset{V_e}{\overset{}{\int }}\left({\delta p}^{\prime }∆p^{\prime }\right)dV,$$

(2)

where ${\rho }_0\left[kg{m}^{-3}\right]$ is the density of the environment, $c_0\left[m{s}^{-1}\right]$ is the velocity of sound, $r\left[kg{m}^{-3}{s}^{-1}\right]$ is the volume damping coefficient, $V_e\left[m^3\right]$ is the element volume and $∆\left[m^{-2}\right]$ is the Laplace operator. The acoustic pressures and their variations expressed using nodal values and their values are put into the equation. The equation is converted into the frequency domain and Green’s first identity is applied. The pressure gradient in the resulting surface integral can be expressed using acoustic velocity

$$\delta {\mathit{P}}_e^T\underset{V_e}{\overset{}{\int }}\left(-{\omega }^2\frac{1}{c_0^2}{\mathit{N}}_e^T{\mathit{N}}_{\mathit{e}}+i\omega {\mathit{N}}_e^T\frac{r}{c_0^2{\rho }_0}{\mathit{N}}_{\mathit{e}}+{\left(\nabla {\mathit{N}}_{\mathit{e}}\right)}^T\nabla {\mathit{N}}_{\mathit{e}}\right)dV{\mathit{P}}_e=-\delta {\mathit{P}}_e^T\underset{{\partial V}_e}{\overset{}{\int }}\left({\mathit{N}}_e^T{\mathit{n}}^T{\rho }_{0}i\omega {\mathit{v}}^{\prime }\right)dS,$$

(3)

where $\omega \left[s^{-1}\right]$ is the angular frequency, $\nabla \left[m^{-1}\right]$ is the nabla operator, $\mathit{n}[-]$ is the unit normal to the boundary ${\partial V}_e\left[m^2\right]$ and ${\mathit{v}}^{\prime }\left[m{s}^{-1}\right]$ is the acoustic velocity. The expressions can be simplified as follows

$$\delta {\mathit{P}}_e^T\left(\left(-{\omega }^2{\mathit{M}}_e+i\omega {\mathit{B}}_e+{\mathit{K}}_e\right){\mathit{P}}_e-{\mathit{V}}_e\right)=0,$$

(4)

where the ${\mathit{M}}_e\left[s^{2}m\right]$ is the element’s mass matrix

$${\mathit{M}}_e=\frac{1}{c_0^2}\underset{V_e}{\overset{}{\int }}\left({\mathit{N}}_e^T{\mathit{N}}_{\mathit{e}}\right)dV,$$

(5)

the ${\mathit{B}}_e\left[sm\right]$ is the element’s damping matrix (the relationship between the specific acoustic impedance $Z_s$ and acoustic velocity ${\mathit{v}}^{\prime }$ was applied here)

$${\mathit{B}}_e=\frac{1}{c_0^2{\rho }_0}\underset{V_e}{\overset{}{\int }}\left(r{\mathit{N}}_e^T{\mathit{N}}_{\mathit{e}}\right)dV+{\rho }_0\underset{{\partial V}_e}{\overset{}{\int }}\left({\mathit{N}}_e^T\frac{1}{Z_s}{\mathit{N}}_{\mathit{e}}\right)dS,$$

(6)

the ${\mathit{K}}_e\left[m\right]$ is the element’s stiffness matrix

$${\mathit{K}}_e=\underset{V_e}{\overset{}{\int }}\left({\left(\nabla {\mathit{N}}_{\mathit{e}}\right)}^T\left(\nabla {\mathit{N}}_{\mathit{e}}\right)\right)dV,$$

(7)

and ${\mathit{V}}_e\left[kg{s}^2\right]$ is the vector of the right sides

$${\mathit{V}}_e=\underset{{\partial V}_e}{\overset{}{\int }}\left({\mathit{N}}_e^T{\mathit{n}}^T{\rho }_{0}i\omega {\mathit{v}}^{\prime }\right)dS.$$

(8)

Relationship (4) mathematically defines the acoustic properties of one element. After discretization of the space of finite elements, global matrices describing the whole system solved are then assembled.

Figure 1 illustrates the FEM mesh model of the VT for vowel [a:]. This model has a total of 36553 nodes. The mesh model is not structed and the types of elements used are quadratic tetrahedrons with ten nodes. There are 186570 elements in the model. The model is created using commercial software Ansys.

Air flow is not stationary due to vocal cord oscillations. The air flow volume is periodic and can be approximated using an analytical model. The Liljencrants-Fant model (L-F model) [14] is one of the basic models of air pulses. This model has four parameters. The first derivative of this pulse is in Figure 2. A velocity-type boundary condition is modelled at the level of the vocal cords.

Acoustic energy is emitted by the mouth into the surrounding space. The external environment model is replaced by the emitted acoustic impedance. This is derived under the assumption of a circular oral opening, constant acoustic velocity throughout the oral opening and an infinite external environment [15]. In such a case, the specific acoustic impedance at the output of the VT can be defined by the following relationship.

$$Z_s\left(\omega \right)=c_0{\rho }_0\left(1-\frac{J_1\left(\frac{2\omega R}{c_0}\right)}{\frac{\omega R}{c_0}}+i\frac{H_1\left(\frac{2\omega R}{c_0}\right)}{1+\frac{\omega R}{c_0}}\right),$$

(9)

where $J_1[-]$ and $H_1[-]$ are the first order Bessel and Struve functions and $R\left[m\right]$ is the radius of the considered oral opening. The disadvantage of this procedure is the non-linear frequency dependance of the acoustic impedance.

A disadvantage of the described FEM model is the large number of DoF. The dimensions of the model’s mass, damping and stiffness matrices equal to the number of the DoF model. Although the matrices are symmetrical and sparse, their dimensions are generally so large (in this case: 36 553) that the model can be described as very computationally demanding. The computation demands depend on the type of analysis, the computer’s performance and the performance optimization of the solver.

3. Aim of This Work

Calculating the sensitivity of acoustic properties on the geometric modifications of the VT or calculating the optimal geometry VT for achieving the required properties are important tasks for clinical practice. Such tasks are generally more demanding and are oftentimes solved using optimization approaches, which require ample calculations given various modification models. The number of analyses is then considerable, requiring more time for calculations.

The aim of this work is to optimize the model’s computation demands, allowing complex optimization calculations to be done. The first procedure uses concept of super-elements obtained by reducing an extensive FEM model divided into smaller subsystems. The second method uses a specially designed element for FEM. This element is designed in such a way so that it enables the discretization of the VT geometry. The models assembled using this element have a very low number of DoF. They only contain the necessary nodes on the surface of the model for defining the geometry and nodes on the curve passing through the centre of the VT cross-section. Both procedures generate models with a similar number of DoF and node spacing.

4. IRS Reduction

The commonly used method for reducing the computing demands of large models are reductions at the level of the M and K matrices. The principle consists in decreasing the number of DoF. The resulting reduced model can be understood as one super-element. The basic reduction method is the Guyan Reduction Method [16], which is generalized by the IRS (improved reduced system) reduction.

IRS reduction is based on a model containing stiffness and mass matrices. The model’s DoF is divided into those that are preserved, ${\mathit{P}}_m$ and those that are reduced, ${\mathit{P}}_s$ [17]

$$\mathit{P}=\left[\begin{array}{c}{\mathit{P}}_m\\ {\mathit{P}}_s\end{array}\right].$$

(10)

The mathematical model transforms into the following shape

$$\left(-{\omega }^2\left[\begin{array}{cc}{\mathit{M}}_{m,m}& {\mathit{M}}_{m,s}\\ {\mathit{M}}_{s,m}& {\mathit{M}}_{s,s}\end{array}\right]+\left[\begin{array}{cc}{\mathit{K}}_{m,m}& {\mathit{K}}_{m,s}\\ {\mathit{K}}_{s,m}& {\mathit{K}}_{s,s}\end{array}\right]\right)\left[\begin{array}{c}{\mathit{P}}_m\\ {\mathit{P}}_s\end{array}\right]=\left[\begin{array}{c}{\widehat{\mathit{V}}}_m\\ 0\end{array}\right].$$

(11)

For the reduction, it is necessary to find map $\mathit{P}$ in the ${\mathit{P}}_m$ set. This is ensured by transformation matrix ${\mathit{T}}_{IRS}$ by the following relationship

$$\mathit{P}=\left[\begin{array}{c}{\mathit{P}}_m\\ {\mathit{P}}_s\end{array}\right]={\mathit{T}}_{IRS}{\mathit{P}}_m..$$

(12)

Transformation matrix ${\mathit{T}}_{IRS}$ is used to transform the whole mathematical description into the reduced DoF

$${\mathit{K}}_{red}={\mathit{T}}_{IRS}^T\mathit{K}{\mathit{T}}_{IRS},$$

(13)

$${\mathit{M}}_{red}={\mathit{T}}_{IRS}^T\mathit{M}{\mathit{T}}_{IRS}$$

(14)

is defined by the following relationship [18]

$${\mathit{T}}_{IRS}={\mathit{T}}_s+\mathit{S}\mathit{M}{\mathit{T}}_s{\mathit{M}}_{red}^{-1}{\mathit{K}}_{red},$$

(15)

where ${\mathit{T}}_s$ is the transformation matrix of the static Guyan reduction

$${\mathit{T}}_s=\left[\begin{array}{c}\mathit{E}\\ -{\mathit{K}}_{s,s}^{-1}{\mathit{K}}_{s,m}\end{array}\right],$$

(16)

and matrix $\mathit{S}$ is defined by the following relationship

$$\mathit{S}=\left[\begin{array}{cc}0& 0\\ 0& {\mathit{K}}_{s,s}^{-1}\end{array}\right].$$

(17)

The formula for ${\mathit{T}}_{IRS}$ is implicit because it depends on the reduced system. Therefore, it has to be calculated iteratively. The first iteration step is solved using the Guyan reduction ${\mathit{T}}_{IRS,1}={\mathit{T}}_s$. After obtaining the first reduced mass and stiffness matrices, iteration can continue using the IRS method

$${\mathit{T}}_{IRS,i+1}={\mathit{T}}_s+\mathit{S}\mathit{M}{\mathit{T}}_s{\mathit{M}}_{red,i}^{-1}{\mathit{K}}_{red,i}.$$

(18)

The convergence of this iteration scheme is accelerated if ${\mathit{T}}_s$ is replaced by ${\mathit{T}}_{IRS,i}$ in the second part. The model’s final reduction matrix is determined after the convergence. The shape of the model is then

$$\left(-{\omega }^2{\mathit{M}}_{red}+{\mathit{K}}_{red}\right){\mathit{P}}_m={\mathit{V}}_m.$$

(19)

5. Properties of IRS Reduction

The derivation of IRS reduction [18] utilizes simplified assumptions, which cause errors in the calculated model. This error depends on the amount and location of the reduced and preserved nodes and on the number of iterations of the IRS method. Therefore, the IRS reduction effect test is conducted on the model’s properties. The motivation to use reduction methods to modify models was to reduce computation demands, which allow optimal calculations to be carried out in reasonable times. The expression for the IRS method transformation matrix is implicit and is solved iteratively, thus making the calculation more complex. Therefore, the convergence test of the IRS iteration scheme was conducted. This test is important for optimally selecting the number of iterations, which guarantee compromise between computation demands and the generated error.

The dominant frequency formants are used to study the importance of VT for forming vowels. From the acoustics viewpoint, the formant is a local frequency extreme in the voice spectrum, which is created by amplifying the acoustic pressure by transfer function of the VT between the vocal cords and the mouth. These local extremes are always in the area around eigenfrequency of the acoustic oscillations of the solved geometry. Therefore, the task then turns into finding eigenfrequencies.

For the purposes of this work, the geometric model of the measured VT had to be divided into subsequent volumes. The cross-sections are always planar and the resulting geometry is illustrated in Figure 3. The main VT cavity is made up of 16 volumes. Furthermore, the model contains another two parts forming the piriform sinuses and 4 parts that close the model. The vallecula are one part of the main cavity. All parts of the divided VT were reduced by the IRS method. The aim of the reduction was to only preserve the nodes at the intersection of the dividing planes with the VT surface. This reduced the total number of DoF from 36553 to 1154. Almost 97% of the model size is reduced. Figure 4 illustrates the size of the relative error of eigenfrequency models in relation to the number of reduced DoF. These graphs are used to evaluate the smallest eigenfrequencies for every VT part and every analysed degree of reduction. To reduce the amount of data in the image, the eigenfrequency with the worst result is always displayed. It can be seen that there is practically no error up to a reduction of 70%-80% DoF. In the case of only preserving the DoF on the perimeter of the dividing surface, the errors are in tenths of percent. Different degrees of reduction were achieved at different parts of the VT. This is given by the fact that the ratio between the number of preserved DoF and the total number of DoF is different. In some parts, only 86% of DoF were reduced and maximum reduction reached 97%.

The computation times increase with the number of iterations. The optimal number of iterations was determined by a numerical experiment. The relative error of eigenfrequency was monitored depending on the number of iterations (Figure 5). The eigenfrequencies were re-evaluated and the worst result was always used for illustration purposes. The error is up to 68% after the first iteration and the average error is approximately 20%. The error decreases steeply with the number of iterations. The decrease decelerates after the fifth iteration. The error does not exceed 1.7% after the 10th iteration and the maximum error after the 20th iteration is 1.08%.

The test conducted using the IRS method only evaluated eigenfrequency, showing a low distortion of results. After the IRS reduction, the model has a lower number of DoF and thereby even a lower number of eigenfrequencies and oscillation shapes. The IRS preserves the lowest frequency. The model preserves the lowest eigenfrequency even when it does not contain the nodes necessary to describe the respective oscillation shape.

Super-elements for all parts of the VT were created. According to Figure 3, this represented a total of 22 super-elements, based on which the entire VT calculation model was assembled.

Such a created model has 1154 nodes (Figure 6), which represents about 3% of the original model. The model’s computation demands increase with its size faster than linearly. The processing time savings are high. In this case, no emitted acoustic impedance was used at the level of the mouth. This boundary condition was modelled as the null sound pressure. This can be interpreted as emission without any losses. The quality of the results is evaluated by comparing the values of eigenfrequencies. The frequencies of the reduced models are as follows

$$f=\left[600,1721,1904,2704,3071,3193,3689,4257\right]Hz.$$

(20)

Using super-elements and reduction methods confirmed the applicability for decreasing computation demands. Comparing the size of the original model with the reduced model is not sufficient for evaluating the effectiveness of these procedures. It is also necessary to evaluate the computation demands of the reduction itself and the number of runs of the created model. If the number of runs of the reduced calculation model is high, it can be assumed that this procedure is worthwhile. However, it is necessary to add up the reduction demands for a single run and any subsequent calculations. Before starting the iteration, it is necessary to calculate ${\mathit{K}}_{s,s}^{-1}$ and it is necessary to calculate the reduced mass and stiffness matrices (simply by multiplying the matrices) and invert the reduced mass matrix for every iteration. Most of the computing time during the IRS is spent on these inversions. The ${\mathit{K}}_{s,s}$ matrix is square-shaped with a dimension corresponding to the number of reduced DoF. On the contrary, the dimension of the reduced mass matrix is the same as the number of preserved DoF. If the number of reduced DoF is relatively low, calculating ${\mathit{K}}_{s,s}^{-1}$ is not demanding but repeatedly calculating ${\mathit{M}}_{red}^{-1}$ is very demanding. On the contrary, if the DoF is primarily reduced, then calculating ${\mathit{K}}_{s,s}^{-1}$ is demanding and calculating ${\mathit{M}}_{red}^{-1}$ is not demanding. Generally speaking, a reduction in which less DoF are preserved is even more effective from the viewpoint of the demands of the reduction itself. In this case, approximately 97% of the DoF are reduced, which means that the dimensions of the ${\mathit{K}}_{s,s}$ matrix is almost the entire stiffness matrix. If the majority of the DoF are reduced, then the demands correspond to the number of inversions of the stiffness matrix. This calculation is comparable to the demanding calculations of general eigenvalues or eigenfrequencies. The IRS method of creating super-elements is only practical when the resulting model is to be used repeatedly. In addition, implementing IRS reduction makes most modifications of the model impossible.

6. The Designed Isoparametric Element

The human vocal cord has a characteristic geometry, which can be described as a curved canal, with a dominant longitudinal dimension and less dominant transverse dimensions. The cross-section is different from the circular one. The requirements for a new FEM element are that it is suitable for VT discretization given the model’s low number of DoF. With regard to the defined requirements, the designed element has a triangular base in the reference state, which extends to the third dimension (Figure 7). The triangular base is based on the shape functions of a higher degree polynomial but most of the nodes are analytically eliminated. The base only has nodes at the vertices and in one wall.

The element’s topology is chosen in such a way, so that the side, where more nodes are located, forms the surface of the VT and the nodes that are outside of it lie in the centre of the cross-section. The VT is longitudinally divided, with two cross-sections along every part, which is discretized by several elements. Since the outer sides of the surface contain many nodes, it is possible to take into consideration the complex shape of the VT’s cross-section. The base functions of the 6th degree and 6 elements were used in every cross-section. Overall, every cross-section is described using 36 nodes on the surface and one central node. The elements in the longitudinal direction are linear. Therefore, it is necessary to cut a sufficient number of cross-sections in the VT geometry.

The derivation of the base functions is based on a 2D triangular 6th degree polynomial element (Figure 8 on the left). The derivation of such a 2D element is described in [19]. Area coordinates are used here (Figure 8 on the right), which are defined as a proportion of the defined area to the area of the element [20]. The first step in the algorithm is assigning the $L_1,L_2,L_3$ coordinates to all nodes. Parameters $\alpha ,\beta ,\gamma$ for each node are calculated from the node coordinates

$$\alpha =n{L}_1,\beta =n{L}_2,\gamma =n{L}_3,$$

(21)

where $n$ is the degree of the element. It is now possible to express an explicit relationship for the node’s shape function

$$N_i\left(L_1,L_2,L_3\right)={\eta }_{{\alpha }_i}\left(L_1\right){\eta }_{{\beta }_i}\left(L_2\right){\eta }_{{\gamma }_i}\left(L_3\right),$$

(22)

where

$$\begin{gathered} {\eta }_{{\alpha }_i}\left(L_1\right)=\prod _{i=1}^{{\alpha }_i}\left(\frac{n{L}_1-i+1}{i}\right)pro{\alpha }_i\ge 1, \\ {\eta }_{{\alpha }_i}\left(L_1\right)=1pro{\alpha }_i=0. \end{gathered}$$

(23)

Functions ${\eta }_{{\beta }_i}$ and ${\eta }_{{\gamma }_i}$ are defined analogically to function ${\eta }_{{\alpha }_i}$ [19]. Parameters $\alpha ,\beta ,\gamma$ determine for which node the given shape function is valid. This procedure creates a full-fledged element. The desired element does not contain any nodes. Therefore, the redundant nodes are eliminated. Redundant nodes are eliminated by dividing its shape function into proportions, which are added, in the correct ratio, to the preserved nodes.

Once the redundant nodes are eliminated, the redundant surface coordinates are converted into independent coordinates $\eta ,\xi$ [20]

$$L_1=1-\eta -\xi ,L_2=\eta ,L_3=\xi .$$

(24)

The resulting system of the triangular base shape functions of the derived element is as follows

$${\mathit{N}}^T\left(\eta ,\xi \right)=\left[\begin{array}{c}\begin{array}{c}1-\xi -\eta \\ -\frac{\eta }{20}\left(648{\eta }^4\xi -972{\eta }^3\xi +558{\eta }^2\xi -117\eta \xi +648{\xi }^5-1620{\xi }^4+1530{\xi }^3-675{\xi }^2+157\xi -20\right)\end{array}\\ \frac{9}{5}\eta \xi \left(108{\eta }^4-180{\eta }^3+105{\eta }^2-25\eta +108{\xi }^4-252{\xi }^3+213{\xi }^2-77\xi +12\right)\\ -\frac{9}{4}\eta \xi \left(216{\eta }^4-396{\eta }^3+246{\eta }^2-61\eta +216{\xi }^4-468{\xi }^3+354{\xi }^2-107\xi +15\right)\\ 2\eta \xi \left(324{\eta }^4-648{\eta }^3+441{\eta }^2-117\eta +324{\xi }^4-648{\xi }^3+441{\xi }^2-117\xi +20\right)\\ -\frac{9}{4}\eta \xi \left(216{\eta }^4-468{\eta }^3+354{\eta }^2-107\eta +216{\xi }^4-396{\xi }^3+246{\xi }^2-61\xi +15\right)\\ \frac{9}{5}\eta \xi \left(108{\eta }^4-252{\eta }^3+213{\eta }^2-77\eta +108{\xi }^4-180{\xi }^3+105{\xi }^2-25\xi +12\right)\\ -\frac{\xi }{20}\left(648{\eta }^5-1620{\eta }^4+1530{\eta }^3-675{\eta }^2+648\eta {\xi }^4-972\eta {\xi }^3+558\eta {\xi }^2-117\eta \xi +157\eta -20\right)\end{array}\right].$$

(25)

The 2D triangular element is extracted into the spatial dimension by a linear shape function. The planar element exists in the $\eta ,\xi$ plane and its extraction is on the $\zeta \in 〈\mathrm{0,1}〉$ axis. Two new shape functions are created from the shape function of every node. This results in a final system of 16 shape functions of the element from Figure 7

$${\mathit{N}}^T\left(\eta ,\xi ,\zeta \right)=\left[\begin{array}{c}\left(1-\zeta \right){\mathit{N}}^T\left(\eta ,\xi \right)\\ \zeta {\mathit{N}}^T\left(\eta ,\xi \right)\end{array}\right].$$

(26)

The derived shape functions are valid in the element’s volume. The element’s volume is defined by the set, $V$ in the space of local coordinates as

$$V=\left\{\begin{array}{c}\eta \in \mathbb{R}|0\le \eta \le 1\\ \xi \in \mathbb{R}|0\le \xi \le 1-\eta \\ \zeta \in \mathbb{R}|0\le \zeta \le 1\end{array}\right\}.$$

(27)

The resulting element is a polynomial of the 6th degree in a planar base and the 1st degree in the third dimension. This allows complex shapes to be modelled by only using 6 elements in every cross-section. In order for the derived element to be correct for FEM, its conditionality has to be checked in all geometric configurations, in which it is used. This was verified by numeric tests, which verified its conditionality on the tested cross-sectional shapes. Conditionality was only tested in 2D. Figure 9 illustrates the conditionality of the elements with a circular cross-section shape. The conditionality ranges at intervals of approximately (1.7, 2.6). Ideally, the conditionality is 1 and infinite conditionality means singularity. When using numbers with a floating decimal point and 8 bytes per number, the accuracy is approximately 15 decimal places. The conditionality is proportional to the increase in the range of the matrix values during inversion and thus decreasing calculation accuracy. The conditionality values in tens are not an obstacle for the calculation.

A human’s vocal tract is geometrically more complex than a rotating cylinder (circular cross-section). Therefore, a second conditionality test was conducted. Figure 10 illustrates the test for significantly non-convex geometry. Here, it depends on how the elements are distributed in the geometry. Therefore, 2 variants are shown. The worst conditionality found is illustrated in 10.3.

The element mass matrices and stiffness matrices for $c_0=1m{s}^{-1}$ are full. The structure of the mass and stiffness matrices are illustrated in Figure 11. This shows that the larger values of the stiffness matrix are located in the node positions, which lies in the element’s wall, creating the VT shell. On the contrary, larger mass matrix values are found in the node positions that do not lie in this wall.

7. Properties of the Derived Isoparametric Element

Combining six derived elements results in a prismatic volume, with a node in the centre of the bases and 36 nodes along the perimeter of the bases, so 74 nodes in total. This is similar to the super-element amount created by the IRS method. Each of these groups of six elements can fill the VT volume surrounded by two cross-sections.

There is a significant difference between both types of elements. The IRS approach explicitly does not contain internal nodes, but their effect is included in the resulting super-element, while the derived isoparametric element does not contain these nodes. This results in it being able to only affect such oscillation shapes, which can be described using the existing nodes. This means that when only one group of six elements is used, it is not possible to find a higher longitudinal oscillation shape. The same problem arises with the transverse oscillation shapes. When using a sufficient number of layers assembled of a group of six derived elements, the qualitative results are the same as for the full-scale FEM models reduced by the IRS method. An advantage of the super-elements from the IRS method is the ability to find eigenfrequency of the oscillation shape, which cannot be described by the nodes used. In contrast, the model from the derived elements is incomparably numerically more effective. Assembling one group of six derived elements is direct. The shape functions are derived in general and transforming them into the necessary geometry is very easy. Therefore, creating their stiffness and mass matrices is also easy. This is done only 6x for one layer, while creating a super-element represents assembling a full FEM model from ordinary elements contains thousands of DoF. When reducing nodes, it is necessary to invert the matrix by a dimension equal to the number of reduced DoF and then invert the resulting mass matrix in every iteration. If tens of nodes are to exist after the reduction, it is then necessary to invert the matrix by thousands of rows. This is computationally demanding. Overall, assembling a model directly from derived elements is 250x more effective than simply reducing a full model, where the time it takes to assemble a full model is not taken into consideration. Moreover, if the time to assemble a new model were taken into consideration, this ratio would be even greater. This result was calculated from real computing times on a regular computer. Since it is a ratio, the same ratio can be expected on different hardware and on differently optimal codes.

The VT model (Figure 12) for vowel [a:], divided into parts according to Figure 3, can be discretized by the derived isoparametric elements, using sixteen layers of elements. The VT is modelled to include parallel cavities. The piriform sinuses and the vallecula are modelled by six additional elements that are attached to the main cavity from the side. Overall, the model contains 949 nodes. There are 17 planes with nodes in the main VT cavity, which is sufficient for modelling lower longitudinal oscillation shapes, which are the most important for voice creation. The velocity of sound is chosen at 350m/s, the VT walls are assumed to be stiff. The VT is excited by speed pulses of the L-F model. It is not necessary to know the shape of these pulses to find eigenfrequencies because the boundary condition used will be of homogeneous velocity. The boundary condition at the VT output in the mouth is defined as the emitted acoustic impedance.

The model’s eigenfrequencies are illustrated in Table 1, along with a description of the respective oscillation shape. The first 4 oscillation shapes are longitudinal. The first transverse oscillation shape is in the widest part of the VT with a frequency of 4261.9Hz. For the human voice, the most significant frequency ranges between 300Hz to 5kHz. There are 5 longitudinal oscillations and one transverse oscillation in this area. The first 4 longitudinal shapes are illustrated in Figure 13 and the first transverse shape is illustrated in Figure 14.

$$F_1=530.4Hz,F_2=1146.6Hz,F_3=3061.2Hz,F_4=3876.1Hz$$

.

The character of the vowel is given by the formants. These are the most represented frequencies in the voice spectrum. When the VT is excited with speed pulses, the spectrum of the pulses is modulated by the VT transmission characteristic. The frequency peaks represent formants if the excitation is broadband. The transmission characteristic is in Figure 15, which shows the significance of the character of the oscillation shape. The first 6 resonance peaks correspond to the longitudinal oscillation frequencies. There are also 2 transverse oscillations among these two shapes, which are not visible on the transmission characteristics. This is given by the position of their nodes at the excitation and emission points. Therefore, the longitudinal oscillation shapes are important for voice creation.

8. Energy Maximalization in the Frequency Band

From the viewpoint of the vocal abilities of professional singers, it’s advantageous to be able to amplify certain frequencies of the spectrum when singing. The human ear is most sensitive around 4kHz. If the singing dominantly represents frequencies around 3kHz, then this singing is much clearer and when sung with an orchestra, the orchestra does not drown out the singing. This is called a singing formant and it is primarily important for opera singers.

The singing formant is the subject of research [21,22]. It was found that it can be achieved by an appropriate ratio of the VT cross-section in the epiglottis and the piriform sinuses cross-section [23]. Therefore, the assembled model was used to verify the hypothesis. The FEM model was parametrized to allow for changes in the cross-section of the main VT cavity in all cross-sections and the piriform sinuses cross-section. The ratio of these cross-sections was changed and the effect on the frequency of the formants was monitored (Figure 16).

The analysis showed that given a suitable cross-section ratio, the third and fourth formants are brought closer together, to a distance of 224Hz from the original 815Hz. Bringing the two resonance peaks closer together results in a higher amplification in the transmission characteristics in Figure 17. This can be interpreted as the singing formant. This double peak is located in the area around 3kHz, which is consistent with observations done on opera singers.

Figure 17 compares the frequency characteristics of the model modified for creating the singing formant with the original geometry model. Here, it can be seen that the geometric modification neither affected the area below 2 kHz nor the area above 7 kHz. Since it did not affect the frequency of the first 2 formants, the phonated sound did not change. The maximum amplification in the area of the singing formant (Figure 18) increased from the original 60 dB to 79 dB, which corresponds to a ten-fold amplification. The average amplification in the 2750-3250 Hz band increased from 43 dB by 17 dB to 60 dB, which also corresponds to about a ten-fold increase in acoustic pressure.

To evaluate the amount of acoustic energy in the monitored frequency band, it is necessary to excite the transmission function. Setting up the excitation of the L-F model was selected by basic frequency pulses of 100Hz, 𝑡_𝑒 = 0.006𝑠, 𝑡_𝑎 = 0.00023𝑠, 𝑡_𝑝 = 0.0045𝑠 and 𝐸_𝑒 = 0.4 𝑚³𝑠⁻². Acoustic energy is defined by the following relationship [15]

$$w=\frac{1}{2}{\rho }_0{\left|\mathit{v}\prime \left(t\right)\right|}^2+\frac{1}{2}\frac{{p\prime \left(t\right)}^2}{{\rho }_0{c}_0^2}.$$

(28)

Acoustic energy is made up of a potential (pressure) and a kinetic part, which changes over time. Only the pressure part of the energy is used to evaluate acoustic energy. After converting it to the frequency domain, acoustic energy is defined as follows

$$w\left(\omega \right)=\frac{1}{2}\frac{p\prime \left(\omega \right)\tilde{p}\prime \left(\omega \right)}{{\rho }_0{c}_0^2}.$$

(29)

Exciting VT pulses generated by the vocal cords are periodic, they contain a fundamental frequency and higher harmonic components. Their spectrum is not continuous. Acoustic energy is then simply a sum of the frequency components. A comparison on the singing formant shows a 111x higher amount of energy in the monitored frequency band.

9. Conclusion

The new isoparametric element was derived for the direct assembly of FEM VT models. In parallel, a spatial FEM VT model was assembled from ordinary elements. Super-elements were created from it using the IRS method. The created VT models were compared. The aim of these procedures is to reduce the computation demands of the model. The VT geometry for vowel [a:] was obtained from CT measurements and it is very detailed. All models even include piriform sinuses and epiglottic vallecula. All of the models used achieve satisfactory compliance in terms of the results of the formant position of [a:]. However, the calculation demands vary significantly.

The FEM model created from ordinary elements has 36553 DoF and 186570 number of elements. Assembling the model requires calculating the stiffness and mass matrices for each element. This is done by numerically integrating the shape functions over a selected number of integration points. The solution then lies in analysing the resulting matrices.

The analyses of the resulting matrices can be more effective by reducing the matrices using reduction methods. In such a case, it is necessary to first assemble a full FEM model and reduce it. If the degree of reduction is high, then the largest number of demands used on the inversion of the stiffness sub-matrix takes place during the IRS procedure. In terms of size, it is close to a full matrix. Even though the IRS method reduces the DoF from 36553 to 1154, thereby decreasing the calculation demands of the formants to negligible, this saving is negatively compensated by calculating the inversion stiffness sub-matrix. As a result, this procedure is not suitable for a one-time analysis because the savings correspond to the demands for the reduction itself.

The last option used is applying a new element for FEM. The obvious advantage of this procedure is that it does not require assembling the model from ordinary elements and a calculation-efficient model is created directly, not by simplifying an extensive model. This model contains 132 elements and a 7th accuracy order integration scheme is used (4 layers with 13 integration points in a triangular plane). The total number of integration points is 6864, which is more than 100x less compared to the assembled ordinary model. The total number of DoF is 949, which represents significant savings even in assembling the global matrix since there are lower demands on the computer’s memory. When comparing the demands at the computation time level, the model assembled from the new isoparametric elements is approximately 250x faster than the simple IRS reduction. If the time needed to assemble the full FEM model were added, the ratio would be even greater.

This work compares two approaches for creating effective FEM VT models. It turns out that the generally known IRS method is less beneficial than using the new element described here. There is a considerable difference in efficiency. The model with the new element is an order of magnitude faster. One disadvantage of the new element is its design for the VT geometry. In the acoustic model, which does not geometrically resemble the VT, this element does not offer any advantages. In contrast the IRS method is applicable everywhere. The designed isoparametric element was tested by practical calculation of the singing formant. The calculation confirmed generally valid assumptions for its creation.

The model using the derived isoparametric element was used for a sensitivity study aimed at finding the optimal geometric configuration for creating a singing formant. It has been confirmed that a suitable ratio between piriform sinuses cross-section and the laryngeal part of the VT will bring the 3rd and 4th formant closer together in the area of 3kHz. This singing formant provides a 17dB higher amplification than the original model with the geometry obtained from CT measurement. The amount of acoustic energy is 111x greater than the original model when the LF model is excited by the model.