1 Introduction
The shape of the human upper vocal tract is quite intricate and the study of its acoustics proves a challenging research topic. The cavity that connects the larynx to the mouth opening forms a curved tube with varying irregular crosssectional shape, in turn connected to several side branches, e.g., the nasal tract and the piriform fossae. Since every detail of this complex structure has an effect on voice production, the most typical approach to the analysis of vocal tract acoustics is the employment of threedimensional (3D) simulations [1, 2, 3]. Such simulations capitalize on the computation of pressure wave propagation within precise 3D geometrical reconstructions of real vocal tracts, and are capable of producing highly accurate results. However, the computational power required to solve wave equations in three dimensions is remarkably high, yielding to simulation times that range from several minutes to more than a day for just a few milliseconds of audio [1, 4].
Several research groups explored twodimensional (2D) simulations as an alternative to a precise but costly 3D approach [5, 6, 7, 4]. The reason for this interest is twofold. Compared to the 3D case, 2D wave equation solvers are characterized by much lower computational requirements and can potentially achieve realtime or quasi realtime performance [8, 9]. Furthermore, this class of simulations appears to be capable of preserving most of the geometrical details of the vocal tract that is being analyzed. In 2D a vocal tract is represented by a flat contour, typically extracted from the midsagittal slice of a 3D tube built from its area function [7, 4]; it can include the vocal tract’s original curvature as well as 2D equivalents of its sidebranches [6]. In other words, the only missing information regards the varying crosssectional shape of the 3D vocal tract, which is replaced by a set of circumferences with the same area progression.
However, in a practical scenario the advantages of 2D simulations prove quite limited. As thoroughly discussed by Arnela and Guasch for the case of vowel synthesis [4], a direct midsagittal representation of the tube built from the vocal tract’s area function produces erroneous formant locations and bandwidths, i.e., not matching the results obtained from the 3D acoustic simulation of the same 3D shape. To fix these discrepancies, the extracted 2D contours have to undergo a nonlinear deformation process that leads to a significant downsizing of the vocal tract’s constrictions, in some cases up to an order of magnitude [9]. A direct consequence of this modification regards simulation times. To model such narrow constrictions, 2D simulations need to run with extremely high spatiotemporal resolutions, in turn increasing the computational load of the solver and producing waiting times far from realtime performance. Moreover, the methodology proposed by Arnela and Guasch is applicable only to the case of straight tubes, and it is still unclear how to obtain 2D contours that acoustically match curved and/or branching 3D geometries.
In this paper, we present a novel approach to vocal tract modeling, that stems from the 2D rationale [9] and improves upon it. It capitalizes on an extended 2D FiniteDifference TimeDomain solver (2.5D FDTD), capable of simulating how pressure propagates in 3D geometries that are symmetric at least in one dimension, typically along the axis (with and being the dimensions of the starting 2D scheme). At the core of this approach there is the inclusion within the model of extra impedance terms, that derive from the tube’s depth (i.e., its continuous extension along ) and that are sampled in every point of the scheme; the result is a depth map, that combined with the 2D midsagittal contour of the original 3D geometry allows for a fast simulation of its acoustics, with computational requirements comparable to the case of standard 2D numerical schemes. Depth maps can be retrieved from area functions as well as from full 3D models of real vocal tracts, like Magnetic Resonance Imaging scans; furthermore, the overall 2.5D representation of the analyzed geometry leaves its original proportions intact, discarding the need for the extremely high spatiotemporal resolutions that characterizes 2D simulations. However, as discussed in the next section, the development of the 2.5D model is still in progress and at its current stage of development it can be employed only for the acoustic simulation of straight tubes with circular crosssectional shape.
2 Methods
2.1 Acoustics Equation
Let us consider a sound wave propagating in air, and let us call
the displacement vector of the medium caused by the motion of the wave and
the volume of a generic element of the medium (particle). Now suppose that the medium is contained between two rigid surfaces, that extend towards infinity along both the and the axis without ever intersecting and that are one the reflected image of the other (being the plane of symmetry). In the described scenario, the wave moves in a space that is 3D but constrained on one axis (i.e., ); such a space can be described by the scalar field , which contains the euclidean distance between the points of the two surfaces at the same coordinates. In other words, describes the space in terms of its depth along .If we choose to focus on a particle of volume , with and relatively small, will only have and components due to the presence of the enclosing surfaces. Consequently, during the passage of the sound wave, the particle’s volume change can be expressed as:
(1) 
if is quite small (reasonable assumption in a realistic scenario) and is slowly varying in space. The two products in parentheses represent the first order Taylor Series approximations of the increment of and , on and respectively. Hence, by choosing , the fractional volume change can be written as:
(2) 
If we describe the motion of the particle via Newton’s equation and we combine it with the equation of state (where is the bulk modulus of air), it is possible to prove the following relationship between the displacement and the fractional volume change just obtained:
(3) 
with defined as the speed of sound in air. Finally, by following Fletcher and Rossing [10] and considering independent of time, we can obtain from Equation 3 the following acoustic wave equation:
(4) 
where is sound pressure. The equation describes air wave motion in a 3D space constrained in one dimension; given the similarity with the standard 2D acoustic wave equation, this 3D space can be considered as an extended 2D (2.5D) space. By symmetrically constraining the 2.5D space along the dimension too (hence turning into a circular area), and by assuming planewave propagation only (due to large distance from the source), Equation 4 reduces to Webster’s horn equation [11].
2.2 2.5D FDTD Wave Solver
The procedure followed for the design of the 2.5D wave solver shares many similarities with the case of the 2D vocal tract model described in [9]. We decomposed the 2.5D acoustic wave equation into extended versions of the 2D equation of continuity and 2D equation of motion:
(5)  
(6) 
where is the 2D acoustic particle velocity. Equation 6 does not include any explicit dependence on the depth of the considered space (via ). This is due to the fact that the equation of motion derives from Newton’s second law, whose 2D and 2.5D forms coincide. However, similarly to what proposed by Allen and Raghuvanshi [12], we augmented the obtained standard 2D equation with the scalar field . By varying between 1 and 0, this term allows for the transition between the momentum equation and the enforcement of a prescribed velocity . Through this mechanism it is possible to enforce boundary conditions (more details in the next paragraph) and to simulate dynamic geometries [9].
We applied a standard 2D Yee scheme, where each grid point consists of a squared 2D cell [12, 13]. Per each time step , pressure values are sampled across the whole domain at the center of the cells, while the and components of the velocity vectors are sampled on the edges shared with the right and top neighbor cell respectively. Hence, in a generic cell at discrete coordinates , will be sampled at , while at and at . This leads to the following discrete update rules of Equations 5 and 6, where we denote with the standard discrete spatial derivatives as performed in 2D FDTD:
(7)  
(8) 
with:
(9) 
The terms , and are the components of the depth map of the 2.5D space and they effectively act as extra impedance terms in two dimensions. Their relation with the field will be explained in detail in the next subsection.
In line with the work of Takemoto and Mokthari [1], the vocal tract’s walls are simulated by adapting the local reactive boundary approach originally proposed by Yokota et al. [14]. This is done by means of setting and , with equal to the pressure value sampled from the cell in front of the wall, the unit vector normal to it and directed towards the wall itself and the boundary admittance. Similarly, arbitrary excitation can be injected into the domain via boundary cells (i.e, ), by setting equal to the output velocity of a glottal model; this can be coupled with the vocal tract by feeding back the pressure value of the cell in front of the excitation velocity stream. Alike its 2D version [9], the 2.5D vocal tract model allows for the coupling with both 1D [15, 16] and 2D glottal models [17].
2.3 Depth Map
The depth map is a discretized version of the continuous field , processed according to the described 2D scheme. Since every individual FDTD cell holds three acoustic parameters sampled in three different locations (center, right and top), the discretization of will likewise produce three depth values per each cell, namely , and . As a result, in the FDTD grid the pressure values are aligned with values, while the velocity vector components are aligned with and values.
In a practical scenario, the field
can be replaced by the 3D model of the geometry whose acoustics is being analyzed. Even if typically composed of discrete elements (e.g., vertices), the model can be treated as a continuous geometry via interpolation. However,
always represents a symmetric domain; this is not necessarily true for a 3D model, hence additional processing is required to derive the depth map from it. In the case of a generic 3D tube, the depth map is computed as follows:
a 2D contour is extracted as the intersection between the tube and its midsagittal plane (i.e., 2D vocal tract geometry). The contour itself represents the vocal tract’s walls, while the midsagittal plane becomes the plane of the 2.5D space;

in every wall cell, , and are set to zero;

in every cell outside of the the tube’s walls, , and are assigned a fixed open space depth;

for every cell with discrete coordinates inside the tube’s walls, the two lines perpendicular to the plane and passing by and are intersected with the model (Figure 1). The length of the two resulting segments is assigned to and respectively;

for every cell with discrete coordinates inside the tube’s walls, and are interpolated as and respectively;

for every cell with discrete coordinates inside the tube’s walls, the depth is obtained as ;

a minimum depth threshold is defined and applied in every cell (depth values are clamped). A typical threshold is at least an order of magnitude lower than the smallest nonzero depth value found across the domain.
The result of this process is a singleaxis symmetric equivalent of the original 3D model. The crosssections of such a 2.5D shape have same areas and widths (i.e., extension along ) as the 3D model’s sections; moreover, their shapes can be irregular, but always symmetric along , as showed in Figure 1.
2.4 Model Applicability
The development of the 2.5D vocal tract model is still in progress, yet some of its features can be analyzed already.
The use of the depth map imposes to slightly change the crosssectional shapes of the original geometry by means of singleaxis symmetry. As a result, the acoustics of the 3D and of the 2.5D vocal tracts will be somewhat different in the portion of the spectrum that is beyond 5 kHz [3]. However, the 2.5D representation is capable of preserving curvature on the plane, area progression as well as several details of the varying irregular crosssectional shapes of real vocal tract geometries. These are features beyond the capabilities of 1D and 2D models, and whose combined acoustics effects are yet to be studied.
At the current stage of development, using the 2.5D model to analyze the acoustics of realistic vocal tracts is impractical though. The model capitalizes on the solution of a lossless system (Equations 5, 6). The solver adds wall impedance by means of the 2D walls boundary condition (Section 2.2), but the losses happening on the surface described by the depth map (i.e., 2.5D boundaries) are largely ignored, as the extra depthrelated impedance terms , and model only the effects of the surface’s spatial derivative^{1}^{1}1This scenario is analogous to the case of Webster’s horn equation, where the absolute size of the modeled tube does not affect propagation.. In other words, there is no consistent way of matching the nonplanar modes of a generic 3D shape.
This issue can be easily overcome for the modeling of 3D tubes with circular crosssections, which can also include bending. By extending the calculation of modes in plates proposed by Fletcher and Rossing [10], it is possible to prove that nonplanar modes in each section of a lossless 2.5D tube are identical to the ones found in its 2D midsagittal crosssection. As a result, we can apply the methodology proposed by Arnela and Guasch [4] to match the first nonplanar mode of a circular crosssection, by scaling the radii of every section of the tube by the constant factor and by modifying the depth map accordingly.
3 Experiments and Results
3.1 Model Validation
A comparative study was carried out to validate the model’s precision, by using as reference a highaccuracy 3D Finite Element Method (FEM) [18]. To this end, we computed the transfer functions of a set of vocal tractlike geometries, we extracted the corresponding formants’ positions and compared them with the results obtained with the 3D FEM. The study was carried out for the following static vowels: /a/, /i/ and /u/. We chose to construct the 3D tubes from Story’s area functions [19]. This dataset defines a standard in vocal tracts’ acoustic analysis and is ideal to work with circular crosssections (see Section 2.4).
As further validation, we compared the model’s time performance with two different versions of the same 2D FDTD: a standard serial implementation (2DS), and a highly optimized parallel implementation (2DP). Both the models were timed while computing the transfer function of Story’s /u/. The target of such comparisons was to estimate the extra computational cost introduced by the depth terms, as well as the time boost made possible by a parallel pipeline.
3.2 Experimental Setup
The 2.5D FDTD acoustic wave solver was implemented in MATLAB. The resolution of the FDTD grid was set to mm; this is the largest value that preserves the geometrical details of the tubes and was determined empirically. Any resolution above this threshold produces less accurate results. Figure 2 illustrates the resulting 2.5D tube geometry for vowel /a/. We set the size of the grid to cells, to obtain a domain that could fit any of the three tubes. The temporal resolution of the simulation is restricted by the CourantFriedrichsLewy condition in 2D, , where is the speed of sound ( m/s within the vocal tract). We therefore set s, equal to a simulation rate of 661,500 Hz.
Air density was set to kg/m^{2}s and all wall cells were assigned a standard boundary admittance . Since the 2.5D FDTD does not yet include a way to precisely model the radiation losses, we implemented an openend termination at the mouth opening by imposing a Dirichlet boundary condition in the same fashion of what done in [4]. Following the procedure described in [9], the impulse responses were obtained by exciting the model with a bandpassed velocity pulse injected from the glottal end; the pressure variation was recorded via a microphone placed 3 mm inside the mouth opening, keeping track of computational times. Per each vocal tract, the MATLAB solver serially iterated across the full grid to simulate a total of 50 ms of audio. The application ran on a workstation equipped with an Intel Core i78700K processor.
The same parameters were maintained also for the 2DS and the 2DP, but all the depthrelated impedance terms were removed from their solvers^{2}^{2}2These settings produce incorrect 2D formants. Yet, the 2D models were used for time comparisons only and their output was ignored.. Furthermore, while the 2DS was implemented in MATLAB, the 2DP was implemented as a shader as described in [9] and ran on a Nvidia GTX 1080 graphics card.
3.3 Results
We applied the Fast Fourier Transform to the impulse responses to obtain the vocal tracts’ transfer functions for the three vowels. The positions of the first 8 formants were extracted and compared with the results of the 3D FEM approach. The exact formants’ positions for the 3D FEM can be found in Arnela’s PhD dissertation
[20]. Table 1 shows the formants’ positional differences and the percentage error between the 2D and the 3D simulation results for /a/, /i/ and /u/ respectively. For each vowel, the 2.5D runtime for producing 50 ms of audio was 15 minutes and 21 seconds; for the 2DS it was 14 minutes and 41 seconds, while for the 2DP 1.4 seconds.Formants  /a/  /i/  /u/  

F1 




F2 




F3 




F4 




F5 




F6 




F7 




F8 



4 Discussion and Conclusion
The acoustic analysis of vowels’ area functions is a standard validation methodology that allows to test the fundamental properties of a wave propagation model. In the case of the 2.5 FDTD, all the three vowels are characterized by positional errors that tend to stay below 2%, and in several cases do not go above 1%. To achieve comparable results, a standard 2D FDTD requires a spatial resolution that is almost three times higher than the value used in this test [9]. Such a level of accuracy suggests that 2.5D equations are a valid means to model how waves travel in a constrained space, and that the proposed solver is correct. Furthermore, the affordable spatial resolution remarkably decreases the computational load of the simulation compared to 2D models’ massive grids [4, 7, 9].
The current serial MATLAB implementation of the 2.5D FDTD proves to be extremely lightweight, even when compared with the performance of the equallysized 2DS. The time difference between the two systems is 20 s (2%); this value is quite small and derives from the fact that the instruction sets of the 2.5D and 2DS solvers differ by four multiplications only (see [9]). Given the similar computational load of the two systems, it is safe to assume that an optimized parallel implementation of the 2.5D running on a modern graphics card will reach quasirealtime performance, close to what achieved with the 2DP.
These results gain particular interest in the context of the future development of the 2.5D vocal tract model. We are currently working on the inclusion of freeradiation effects (via 2.5D perfect matching layers) and on the validation of the modeling of the vocal tract’s curvature. The inclusion of these features is relatively straightforward and will make the system capable of simulating geometrical and acoustic details still unavailable in onedimensional and 2D models. The following step will consist of the introduction in the acoustic equations of extra loss terms, to model 2.5D boundaries. By doing so, it will be possible to simulate the effects of varying irregular crosssections, thus fully exploring the potential of the 2.5D model.
5 Acknowledgements
This work is supported by the Natural Sciences and Engineering Research Council (NSERC) of Canada and Canadian Institutes for Health Research (CIHR).
References
 [1] H. Takemoto, P. Mokhtari, and T. Kitamura, “Acoustic analysis of the vocal tract during vowel production by finitedifference timedomain method,” The Journal of the Acoustical Society of America, vol. 128, no. 6, pp. 3724–3738, 2010.
 [2] T. Vampola, J. Horáček, A.M. Laukkanen, and J. G. Švec, “Human vocal tract resonances and the corresponding mode shapes investigated by threedimensional finiteelement modelling based on ct measurement,” Logopedics Phoniatrics Vocology, vol. 40, no. 1, pp. 14–23, 2015.
 [3] M. Arnela, S. Dabbaghchian, R. Blandin, O. Guasch, O. Engwall, A. Van Hirtum, and X. Pelorson, “Influence of vocal tract geometry simplifications on the numerical simulation of vowel sounds,” The Journal of the Acoustical Society of America, vol. 140, no. 3, pp. 1707–1718, 2016.
 [4] M. Arnela and O. Guasch, “Twodimensional vocal tracts with threedimensional behavior in the numerical generation of vowels,” The Journal of the Acoustical Society of America, vol. 135, no. 1, pp. 369–379, 2014.
 [5] J. Mullen, D. M. Howard, and D. T. Murphy, “Waveguide physical modeling of vocal tract acoustics: flexible formant bandwidth control from increased model dimensionality,” IEEE Transactions on Audio, Speech, and Language Processing, vol. 14, no. 3, pp. 964–971, 2006.
 [6] M. Speed, D. Murphy, and D. M. Howard, “Characteristics of twodimensional finite difference techniques for vocal tract analysis and voice synthesis,” in Tenth Annual Conference of the International Speech Communication Association, 2009.
 [7] Y. Wang, H. Wang, J. Wei, and J. Dang, “Mandarin vowel synthesis based on 2d and 3d vocal tract model by finitedifference timedomain method,” in Proceedings of The 2012 Asia Pacific Signal and Information Processing Association Annual Summit and Conference. IEEE, 2012, pp. 1–4.
 [8] J. Mullen, D. M. Howard, and D. T. Murphy, “Realtime dynamic articulations in the 2d waveguide mesh vocal tract model,” IEEE Transactions on Audio, Speech, and Language Processing, vol. 15, no. 2, pp. 577–585, 2007.
 [9] V. Zappi, A. Vasuvedan, A. Allen, N. Raghuvanshi, and S. Fels, “Towards realtime twodimensional wave propagation for articulatory speech synthesis,” in Proceedings of Meetings on Acoustics 171ASA, vol. 26, no. 1. ASA, 2016, p. 045005.
 [10] N. H. Fletcher and T. D. Rossing, The physics of musical instruments. Springer Science & Business Media, 2012.
 [11] G. Ungeheuer, Elemente einer akustischen Theorie der Vokalartikulation. SpringerVerlag, 1962.
 [12] A. Allen and N. Raghuvanshi, “Aerophones in flatland: Interactive wave simulation of wind instruments,” ACM Transactions on Graphics (TOG), vol. 34, no. 4, p. 134, 2015.
 [13] K. Yee, “Numerical solution of initial boundary value problems involving maxwell’s equations in isotropic media,” IEEE Transactions on antennas and propagation, vol. 14, no. 3, pp. 302–307, 1966.
 [14] T. Yokota, S. Sakamoto, and H. Tachibana, “Visualization of sound propagation and scattering in rooms,” Acoustical science and technology, vol. 23, no. 1, pp. 40–46, 2002.
 [15] K. Ishizaka and J. L. Flanagan, “Synthesis of voiced sounds from a twomass model of the vocal cords,” Bell system technical journal, vol. 51, no. 6, pp. 1233–1268, 1972.
 [16] B. H. Story and I. R. Titze, “Voice simulation with a bodycover model of the vocal folds,” The Journal of the Acoustical Society of America, vol. 97, no. 2, pp. 1249–1260, 1995.
 [17] A. Vasudevan, V. Zappi, P. Anderson, and S. Fels, “A fast robust 1d flow model for a selfoscillating coupled 2d fem vocal fold simulation.” in INTERSPEECH, 2017, pp. 3482–3486.
 [18] M. Arnela and O. Guasch, “Finite element computation of elliptical vocal tract impedances using the twomicrophone transfer function method,” The Journal of the Acoustical Society of America, vol. 133, no. 6, pp. 4197–4209, 2013.
 [19] B. H. Story, “Comparison of magnetic resonance imagingbased vocal tract area functions obtained from the same speaker in 1994 and 2002,” The Journal of the Acoustical Society of America, vol. 123, no. 1, pp. 327–335, 2008.
 [20] M. Arnela Coll, “Numerical production of vowels and diphthongs using finite element methods,” Ph.D. dissertation, Universitat Ramon Llull, 2015.
Comments
There are no comments yet.