Applied Science, Faculty of
Electrical and Computer Engineering, Department of
DSpace
UBCV
Mohapatra, Debasish Ray
Master of Applied Science - MASc
University of British Columbia
The human voice is a complex but unique physiological process. It involves the neuromuscular control of articulators to form an intricate upper vocal tract geometry, which yields different speech sounds. The existing computational vocal tract models have significant limitations concerning acoustic precision and simulation performance. The high-dimensional vocal tract models can compute precise acoustic wave propagation at the expense of simulation run-time. This thesis aims to fill these lacunae through two major contributions. Firstly, we introduce a novel vocal tract that extends the existing two-dimensional (2D) vocal tract modelling approach while having three-dimensional behaviour. The proposed model (2.5D FDTD) employs the Finite-Difference Time-Domain numerical scheme to discretize and compute acoustic components on a staggered grid. The simulated acoustic outputs of our new model are shown to match with the 2D FDTD vocal tract model at a low spatiotemporal resolution for open static vocal tract shapes. Contrary to 2D FDTD, the model adds tube depth as an additional impedance parameter to the acoustic wave solver by lumping off-plane waves. This technique offers an excellent balance between computational cost and acoustic precision while promising better geometrical flexibility for vocal tract modelling. Secondly, we tested the model's basic capabilities through the acoustic simulation of cardinal English vowel sounds. For realistic modelling of vowel sounds, we built a vocal tract radiation model. We also couple the 2.5D vocal tract with a self-oscillatory lumped-element vocal fold model to illustrate a fully connected articulatory speech synthesizer. This study offers a speech synthesis tool that can generate static vowel sounds and set up a new pathway for lightweight vocal tract modelling and other computational acoustic research.
Talking Tube: A Novel Approach for Vocal Tract AcousticModelling using The Finite-Difference Time-DomainMethodbyDebasish Ray MohapatraB.Tech, Electronics and Communication Engineering, S‘O’A University, 2013A THESIS SUBMITTED IN PARTIAL FULFILLMENTOF THE REQUIREMENTS FOR THE DEGREE OFMaster of Applied ScienceinTHE FACULTY OF GRADUATE AND POSTDOCTORALSTUDIES(Electrical and Computer Engineering)The University of British Columbia(Vancouver)April 2021c© Debasish Ray Mohapatra, 2021The following individuals certify that they have read, and recommend to the Fac-ulty of Graduate and Postdoctoral Studies for acceptance, the thesis entitled:Talking Tube: A Novel Approach for Vocal Tract Acoustic Modellingusing The Finite-Difference Time-Domain Methodsubmitted by Debasish Ray Mohapatra in partial fulfillment of the requirementsfor the degree of Master of Applied Science in Electrical and Computer Engi-neering.Examining Committee:Sidney Fels, Electrical and Computer Engineering, UBCSupervisorBryan Gick, Department of Linguistics, UBCSupervisory Committee MemberiiAbstractThe human voice is a complex but unique physiological process. It involves theneuromuscular control of articulators to form an intricate upper vocal tract ge-ometry, which yields different speech sounds. The existing computational vocaltract models have significant limitations concerning acoustic precision and simu-lation performance. The high-dimensional vocal tract models can compute preciseacoustic wave propagation at the expense of simulation run-time. This thesis aimsto fill these lacunae through two major contributions.Firstly, we introduce a novel vocal tract that extends the existing two-dimensional(2D) vocal tract modelling approach while having three-dimensional behaviour.The proposed model (2.5D FDTD) employs the Finite-Difference Time-Domainnumerical scheme to discretize and compute acoustic components on a staggeredgrid. The simulated acoustic outputs of our new model are shown to match with the2D FDTD vocal tract model at a low spatiotemporal resolution for open static vo-cal tract shapes. Contrary to 2D FDTD, the model adds tube depth as an additionalimpedance parameter to the acoustic wave solver by lumping off-plane waves. Thistechnique offers an excellent balance between computational cost and acoustic pre-cision while promising better geometrical flexibility for vocal tract modelling.Secondly, we tested the model’s basic capabilities through the acoustic simu-lation of cardinal English vowel sounds. For realistic modelling of vowel sounds,we built a vocal tract radiation model. We also couple the 2.5D vocal tract witha self-oscillatory lumped-element vocal fold model to illustrate a fully connectedarticulatory speech synthesizer. This study offers a speech synthesis tool that cangenerate static vowel sounds and set up a new pathway for lightweight vocal tractmodelling and other computational acoustic research.iiiLay SummaryNaturalness and intelligibility are the two primary characteristics associated withevery speech synthesizer. Articulatory speech synthesizers use the realistic geo-metrical description of the human speech anatomy to produce synthesized speechutterances. Theoretically, they can produce high-quality speech output that soundsnatural and intelligible as they directly follow the human speech production mech-anism. The vocal tract, a critical component of these synthesizers, has an irregulartube structure that deforms during speech production to modulate the source acous-tic energy. The existing vocal tract models are either computationally expensive oroversimplify while approximating the actual vocal tract geometry. In this thesis,we propose a new model (a.k.a, 2.5D FDTD vocal tract), which has better geo-metrical adaptability than a regular two-dimensional vocal tract model. The modelshows the capability of precise simulation of acoustic wave propagation for openstatic vocal tract shapes at a low computational overhead. We couple the pro-posed model with a self-oscillatory vocal fold to produce cardinal vowel soundsand demonstrate an articulatory speech synthesis tool.ivPrefaceThis thesis is an original intellectual product of the author, Debasish R. Mohapatra.And, it is part of the Brain-to-Speech (Brain2Speech) project.Following are the publications resulting from the work described in this thesis.Conference Publications[C1] P. Saha, D. R. Mohapatra, V. P. Srungarapu and S. Fels, “Sound-Stream II:Towards Real-Time Gesture-Controlled Articulatory Sound Synthesis,” CanadianAcoustics 46 (4), pp58-59, 2018.[C2] D. R. Mohapatra and S. Fels, “Limitations Of Source-Filter Coupling InPhonation,” Canadian Acoustics 46 (4), pp60-61, 2018.[C3] D. R. Mohapatra, V. Zappi and S. Fels, “An Extended Two-DimensionalVocal Tract Model for Fast Acoustic Simulation of Single-Axis Symmetric Three-Dimensional Tubes,” Proceedings of Interspeech, pp3760-65, 2019.[C4] D. R. Mohapatra, V. Zappi and S. Fels, “A comparative study of two-dimensionalvocal tract acoustic modeling based on Finite-Difference Time-Domain methods,”ISSP 2020 (Accepted)[C5] P. Saha, D. R. Mohapatra and S. Fels, “Speak With Your Hands - UsingContinuous Hand Gestures to control Articulatory Speech Synthesizer,” ISSP 2020(Accepted)vAuthor’s ContributionIn [C1], all the authors contributed equally towards the conceptualization of ideasand manuscript writing. My main contribution was the development of a hardwareinterface using Arduino to control tongue muscle excitations. Srungarapu helpedextract an inventory of cross-sectional areas from the Artisynth tongue model andcoupled it to the JASS audio synthesis engine.In [C2], I reviewed existing methods to couple vocal fold and vocal tract modelsand their limitations. I wrote the original manuscript, and Fels acted in a supervi-sory role. A modified version of this work is included in Chapter 2.In [C3] and [C4], Zappi and I equally contributed as authors. I developed the2.5D FDTD vocal tract model with Zappi’s guidance. Fels helped in reviewingthe manuscripts and supervising the project. The [C4] characterizes the proposedvocal tract model [C3] with more vocal tract area functions dataset. The work ispartly included in Chapter 4 and 5.In [C5], my main contribution was in interfacing the hand glove to the Pink Trom-bone and data collection. Shah developed the methodology to control the PinkTrombone web interface parameters using hand gestures. I partially contributedtowards writing the manuscript and revising it.Project CodeWe used the MATLAB environment to implement the 2.5D vocal tract model andthe self-oscillatory two-mass vocal fold model. The custom code for this project isavailable here:https://github.com/Debasishray19/Talking-TubeviTable of ContentsAbstract . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . iiiLay Summary . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . ivPreface . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . vTable of Contents . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . viiList of Tables . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . xiList of Figures . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . xivList of Abbreviations . . . . . . . . . . . . . . . . . . . . . . . . . . . . xxiiAcknowledgments . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . xxiiiDedication . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . xxiv1 Introduction . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 11.1 Speech Synthesis Technology . . . . . . . . . . . . . . . . . . . . 11.2 Motivation . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 31.3 Articulatory Speech Synthesis Model . . . . . . . . . . . . . . . 41.4 Vocal Tract Modelling Road Map . . . . . . . . . . . . . . . . . . 81.5 Contribution . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 101.6 Thesis Organisation . . . . . . . . . . . . . . . . . . . . . . . . . 11vii2 Background and Previous Work . . . . . . . . . . . . . . . . . . . . 122.1 Mechanics of Voice Production . . . . . . . . . . . . . . . . . . . 142.2 Articulatory Voice Models . . . . . . . . . . . . . . . . . . . . . 162.2.1 Glottal Source Models . . . . . . . . . . . . . . . . . . . 162.2.2 Vocal Tract Structural Model . . . . . . . . . . . . . . . . 202.2.3 Vocal Tract Acoustic Model . . . . . . . . . . . . . . . . 292.3 Discretizing Two-Dimensional Acoustic Wave Equation Using TheFDTD Method . . . . . . . . . . . . . . . . . . . . . . . . . . . 332.4 Data Acquisition For Vocal Tract Structural Model . . . . . . . . 392.5 Summary . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 393 Limitations of 2D FDTD Vocal Tract Models . . . . . . . . . . . . . 403.1 2D FDTD Modelling Methodology . . . . . . . . . . . . . . . . . 433.1.1 Problem Statement . . . . . . . . . . . . . . . . . . . . . 433.1.2 Step 1: Nonlinear Scaling of 3D Vocal Tract Geometry . . 443.1.3 Step 2: Tuning Grid Resolution . . . . . . . . . . . . . . 473.2 2D FDTD Vocal Tract Transfer Function . . . . . . . . . . . . . . 483.2.1 Transfer Function Analysis . . . . . . . . . . . . . . . . . 513.3 2.5D VT Model - An Extended 2D Vocal Tract Model . . . . . . . 523.4 Summary . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 534 2.5D FDTD Vocal Tract Acoustic Model For Static Vowels . . . . . 544.1 2.5D Vocal Tract Model . . . . . . . . . . . . . . . . . . . . . . . 564.2 2.5D Vocal Tract Acoustic . . . . . . . . . . . . . . . . . . . . . 574.2.1 Acoustic Equation . . . . . . . . . . . . . . . . . . . . . 574.2.2 2.5D FDTD Wave Solver . . . . . . . . . . . . . . . . . . 594.3 Depth Map . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 604.4 Vocal Tract Wall Losses . . . . . . . . . . . . . . . . . . . . . . . 624.5 Model Applicability . . . . . . . . . . . . . . . . . . . . . . . . . 644.6 2.5D Vocal Tract Modelling Methodology . . . . . . . . . . . . . 654.7 Summary . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 68viii5 An Articulatory Model Using The 2.5D FDTD Vocal Tract Model . 695.1 Articulatory Model Formulation . . . . . . . . . . . . . . . . . . 695.2 Model Characterization . . . . . . . . . . . . . . . . . . . . . . . 725.2.1 Characterization Procedure . . . . . . . . . . . . . . . . . 735.3 Experimental Setup . . . . . . . . . . . . . . . . . . . . . . . . . 765.3.1 Simulation Parameters . . . . . . . . . . . . . . . . . . . 765.3.2 Excitation Function . . . . . . . . . . . . . . . . . . . . . 775.3.3 Mouth-End Conditions . . . . . . . . . . . . . . . . . . . 795.3.4 Vocal Tract Transfer Functions . . . . . . . . . . . . . . . 805.3.5 System Specifications . . . . . . . . . . . . . . . . . . . 825.4 Results . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 835.4.1 Acoustic analysis for the 1D one-end closed lossless cylin-drical tube . . . . . . . . . . . . . . . . . . . . . . . . . . 835.4.2 Acoustic Analysis For Open Static Vocal Tract Shapes . . 855.4.3 Acoustic Analysis For Vocal Tract Radiation Models . . . 925.4.4 Time Complexity Analysis . . . . . . . . . . . . . . . . . 945.5 Discussion . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 955.5.1 SA-Analytical: One-end Closed Cylindrical Tube Model . 965.5.2 SA-Vowel: Open Static Vocal Tract Shapes For Vowel Sounds 965.5.3 SA-Radiation: 2.5D Vocal Tract Radiation Model For openstatic vocal tract shapes . . . . . . . . . . . . . . . . . . . 1025.5.4 Formant Bandwidth Analysis . . . . . . . . . . . . . . . 1035.6 Vocal Fold Model . . . . . . . . . . . . . . . . . . . . . . . . . . 1055.6.1 Two-Mass Vocal Fold Model . . . . . . . . . . . . . . . . 1055.7 2.5D FDTD Vocal Tract’s Feasibility For An Articulatory Model . 1075.7.1 Coupled vs Uncoupled Vocal Folds . . . . . . . . . . . . 1075.7.2 Synthesized Vowel Sound Generation . . . . . . . . . . . 1085.7.3 Model Analysis . . . . . . . . . . . . . . . . . . . . . . . 1085.8 Summary . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 1126 Conclusions . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 1136.1 Summary of Contributions . . . . . . . . . . . . . . . . . . . . . 1136.2 Limitations . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 115ix6.2.1 2.5D Vocal Tract Modelling Limitations . . . . . . . . . . 1166.2.2 Experimental Study Limitations . . . . . . . . . . . . . . 1186.3 Future Directions . . . . . . . . . . . . . . . . . . . . . . . . . . 1186.4 Concluding Remarks . . . . . . . . . . . . . . . . . . . . . . . . 121Bibliography . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 122A General Mathematics . . . . . . . . . . . . . . . . . . . . . . . . . . 137A.1 1D Acoustic Wave Equation . . . . . . . . . . . . . . . . . . . . 137A.2 Approximation of first-order derivatives . . . . . . . . . . . . . . 139B 2D Vocal Tract Structural Model . . . . . . . . . . . . . . . . . . . . 141C Acoustic Analysis of 2.5D FDTD Model . . . . . . . . . . . . . . . . 145C.1 1D one-end closed lossless cylindrical tube . . . . . . . . . . . . 145C.2 2.5D simulation of vowel sounds with different grid resolutions . . 147D 2.5D FDTD Time Marching Algorithm . . . . . . . . . . . . . . . . 149xList of TablesTable 2.1 Summary of FRANK model components: component name andcomponent type (FE models (FE), rigid body (Rigid)), adoptedfrom [119] . . . . . . . . . . . . . . . . . . . . . . . . . . . . 23Table 3.1 Formant locations Fi (Hz) and formant deviations ∆Fi/Fi (%)with respect to the 3D FEM model in the existing 2Dlow (with-out scaling); and 2Dlow, 2Dmid and 2Dhigh (with scaling) forvowel /A/. . . . . . . . . . . . . . . . . . . . . . . . . . . . . 50Table 5.1 Tabular representation of a stepwise procedure to analyze theacoustic features of the 2.5D FDTD vocal tract model. Thischapter produces results only for the highlighted section. . . . 72Table 5.2 Various vocal tract geometrical shape types for the 2.5D model. 75Table 5.3 Simulation parameters for 2.5D FDTD vocal tract acoustic model. 77Table 5.4 Positions of the first four formants, Fi (Hz) for vowel /A/, ex-tracted from its frequency response of band-pass and Gaussianvelocity pulses. These input signals are applied separately as anexcitation function to the 2.5D vocal tract model. . . . . . . . . 77Table 5.5 System specifications for 2.5D FDTD numerical simulations . 82Table 5.6 Formant locations Fi (Hz) of the 1D and 2.5Dlow vocal tractmodel for vowel sounds - /A/, /i/, /u/ and /æ/. ∆Fi/Fi (%) presentsthe percentage of formant deviations of Kee’s 1D vocal tractmodel with respect to the 2.5Dlow model. . . . . . . . . . . . . 87xiTable 5.7 Formant locations Fi (Hz) and formant deviations ∆Fi/Fi (%) of2.5Dlow, 2Dhigh and 1D simulations of vowel /A/ with respect tothe 3D FEM vocal tract model. . . . . . . . . . . . . . . . . . 88Table 5.8 Formant locations Fi (Hz) and formant deviations ∆Fi/Fi (%) of2.5Dlow, 2Dhigh and 1D simulations of vowel /i/ with respect tothe 3D FEM vocal tract model. . . . . . . . . . . . . . . . . . 89Table 5.9 Formant locations Fi (Hz) and formant deviations ∆Fi/Fi (%) of2.5Dlow, 2Dhigh and 1D simulations of vowel /u/ with respect tothe 3D FEM vocal tract model. . . . . . . . . . . . . . . . . . 90Table 5.10 Formant locations Fi (Hz) and formant deviations ∆Fi/Fi (%)of 2.5Dlow and 1D simulations of vowel /æ/ with respect to the2Dhigh FDTD vocal tract model. . . . . . . . . . . . . . . . . . 91Table 5.11 Formant locations Fi (Hz) of the 3D and 2.5Dlow vocal tractradiation model for the vowel sound - /u/. ∆Fi/Fi (%) presentsthe percentage of formant deviations of the 2.5Dlow vocal tractmodel with respect to the 3D FEM model. . . . . . . . . . . . 93Table 5.12 Formants derived Fi (Hz) using the 2.5D model simulations.Averaged formants Fi(±s.d) across several seconds of recordedspeech utterances and s.d. is the standard deviation. . . . . . . 93Table 5.13 Simulation run-time (in hh:mm:ss) for the 2D FDTD GPU wavesolver with different domain sizes (grid resolutions) to generate50 ms of audio (open-end boundary condition). Reproducedwith permission from Victor Zappi [146]. . . . . . . . . . . . . 94Table 5.14 Simulation run-time (in hh:mm:ss) for the 2.5D FDTD wavesolver with different domain sizes (grid resolutions) to generate50 ms of audio (open-end boundary condition). . . . . . . . . . 95Table 5.15 Simulation run-time (in hh:mm:ss) for the 2.5D FDTD wavesolver with different domain sizes (grid resolutions) to generate50 ms of audio (free field radiation boundary condition). . . . . 95Table C.1 Formant locations Fi (Hz) and formant deviations ∆Fi/Fi (%)of the 2.5Dlow simulations of one-end closed cylindrical tubemodels (L1 = 17.07cm) with respect to their analytical solution. 145xiiTable C.2 Formant locations Fi (Hz) and formant deviations ∆Fi/Fi (%)of the 2.5Dlow simulations of one-end closed cylindrical tubemodels (L2 = 50.16cm) with respect to their analytical solution. 146Table C.3 Formant locations Fi (Hz) and formant deviations ∆Fi/Fi (%) ofthe 1D and 2.5Dlow simulation of a one-end closed cylindricaltube model (L1 = 17.07cm, A1 = 0.91cm2) with respect to itsanalytical solution. . . . . . . . . . . . . . . . . . . . . . . . . 146Table C.4 Formant locations Fi (Hz) and formant deviations ∆Fi/Fi (%)of the 2.5Dlow, 2.5Dmid and 2.5Dhigh simulations for vowel /A/with respect to the 3D FEM vocal tract model. . . . . . . . . . 147Table C.5 Formant locations Fi (Hz) and formant deviations ∆Fi/Fi (%) ofthe 2.5Dlow, 2.5Dmid and 2.5Dhigh simulations for vowel /i/ withrespect to the 3D FEM vocal tract model. . . . . . . . . . . . . 147Table C.6 Formant locations Fi (Hz) and formant deviations ∆Fi/Fi (%)of the 2.5Dlow, 2.5Dmid and 2.5Dhigh simulations for vowel /u/with respect to the 3D FEM vocal tract model. . . . . . . . . . 148Table C.7 Formant locations Fi (Hz) and formant deviations ∆Fi/Fi (%)of the 2.5Dlow, 2.5Dmid and 2.5Dhigh simulations for vowel /æ/with respect to the 2Dhigh FDTD vocal tract model. . . . . . . 148xiiiList of FiguresFigure 1.1 A graphical illustration of vocal tract modelling roadmap togenerate and evaluate synthesized speech utterances [11, 106].The different steps in the roadmap nomenclature have beencategorized using the following abbreviations- SS: Stationarystructure related, DS: Dynamic structure related, SA: Station-ary acoustic related, DA: Dynamic acoustic related, FT: Vocalfold and vocal tract, RS: Recorded and synthesized speech. Foreach step in our roadmap, we compare the proposed 2.5D vo-cal tract model with the existing models. We analyze geometryto study structural models, whereas, we do spectrum analysisto compare the acoustic models. . . . . . . . . . . . . . . . . 6Figure 1.2 General outline of an articulatory speech synthesizer. . . . . . 7Figure 2.1 Basic Anatomy of the Human Speech Production: larynx, sub-glottal component, supraglottal component, n Anatomy &Physiology, OpenStax CNX, adopted from [88]. . . . . . . . . 13Figure 2.2 Inferior view of the vocal folds from the laryngopharynx, nAnatomy & Physiology, OpenStax CNX, adopted from [88]. . 14Figure 2.3 Variants of one-mass vocal fold models, reproduced with per-mission from Peter Birkholtz [20]. . . . . . . . . . . . . . . . 18Figure 2.4 Coronal view of the two-mass (left) and body-cover (right) vo-cal fold models, reproduced with permission from Peter Birkholz[20]. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 18xivFigure 2.5 Finite element mesh of vocal fold model shown in coronalplane, adopted from [132] . . . . . . . . . . . . . . . . . . . 19Figure 2.6 Illustrations of variations in vocal tract shape for different staticvowel sounds with mid-sagittal vocal tract MRI images : (left)Vowel [A], (middle) Vowel [i], (right) Vowel [u], reproducedwith permission from Peter Birkholz [106] . . . . . . . . . . 21Figure 2.7 (left) An illustration of articulators from the upper vocal tract,adopted from image by Patrick J. Lynch, medical illustratorn Wikimedia Commons. (right) A midsagittal view of thebiomechanical model providing an overview of the upper air-way, adopted from [5]. . . . . . . . . . . . . . . . . . . . . . 22Figure 2.8 Biomechanical tongue models in ArtiSynth for different pho-netic segments: (left) Vowel [A], (middle) Vowel [i], (right)Vowel [u] . . . . . . . . . . . . . . . . . . . . . . . . . . . . 23Figure 2.9 Vocal tract structural model with a number of articulatory pa-rameters, reproduced with permission from Peter Birkholz [22] 23Figure 2.10 Straight vocal tract tube structure for vowel [A]: (top) 3D tubewith various cross-sectional shapes (realistic, elliptical and cir-cular), reproduced with permission from Marc Arnela [10],(bottom) 2D tube with circular cross-section . . . . . . . . . . 25Figure 2.11 (a)A 3D mesh representation of the upper airway and cross-sectional area with red circles. (b) Volumetric vocal tract ge-ometry of vowel [A], intersected by grid planes (blue), repro-duced with permission from Marc Arnela [11]. . . . . . . . . 25Figure 2.12 An illustration of 2D midsagittal vocal tract contour with cen-terline (white) and semi-polar grid axis (red) to extract the areafunction. . . . . . . . . . . . . . . . . . . . . . . . . . . . . 27Figure 2.13 1D Vocal tract acoustic tube models. (a) ACBD matrix basedmodel, (b) Digital Waveguide Model (c) circuit analogy model,reproduced with permission from Marc Arnela [12] . . . . . . 31xvFigure 2.14 Representing pressure and velocity components in a staggeredYee-grid. The acoustic pressure p(x, y) is sampled at the cen-ter of the cell (red circle), vx is sampled at the right edge of thecell (orange) and vy is sampled at the top edge of the cell (blue) 33Figure 2.15 An illustration of the FDTD time-marching algorithm. . . . . 36Figure 2.16 An illustration of acoustic pressure wave propagation insidea rectangular grid with (500× 500) domain size and havingPML layers at the domain border. A sinusoidal function isused as the acoustic source inside the 2D domain. . . . . . . . 37Figure 2.17 An illustration of acoustic pressure wave propagation inside arectangular grid with (500× 500) domain size and having noPML layers at the domain border. A sinusoidal function is usedas the acoustic source inside the 2D domain. . . . . . . . . . . 38Figure 3.1 (left) An illustration of 3D vocal tract geometry for vowel [A]with a straight cylindrical tube having circular cross-sections.(right) A mid-sagittal cut of the 3D vocal tract geometry toretrieve the 2D vocal tract contour. . . . . . . . . . . . . . . . 42Figure 3.2 Due to acoustic wave reflections at the tube junction, the pres-sure pi(t) at any tube segment (i) is the resultant sum of thefi(t) and bi(t), the forward propagating and backward prop-agating pressure wave components. These phenomena (theeffects of wave scattering at the tube junction) are modelledusing the reflection coefficient r′, which depends (top) on theadjacent tube cross-sectional areas ratio for a 3D tube (bottom)on the adjacent tube diameter/tube width ratio for a 2D tube. . 42Figure 3.3 Determining a scaling function r(x) to retrieve the area func-tion for the 2D vocal tract contour, a2D(x) using the original3D vocal tract cross-sections, a3D(x). x represents the lengthof the vocal tract from the glottal-end. This scaling functionensures the spectral features of a 2D model match a 3D vocaltract model. . . . . . . . . . . . . . . . . . . . . . . . . . . . 43xviFigure 3.4 2D vocal tract area function (red line) for vowels /a/, /i/, /u/and /æ/ after nonlinear scaling of the original 3D area function. 46Figure 3.5 Illustration of changes in tube geometry for a 2D vocal tractcontour (vowel [A]) due to area functions scaling and grid res-olution tuning. The x and y axis of the above plots representthe number of grid cells. . . . . . . . . . . . . . . . . . . . . 48Figure 4.1 Highlighted steps (yellow) in the modified roadmap which areaddressed in this chapter. . . . . . . . . . . . . . . . . . . . 55Figure 4.2 Graphical representation of wave propagation in xy-plane boundedby two rigid surfaces (A and B). The euclidean distance be-tween the two surfaces is |a1b1|= D(x,y). . . . . . . . . . . . 57Figure 4.3 on the left, depth map extraction (D(x) solid line, D(y) dashedline) from a slice C of an irregular tube, at coordinate x = xˆ(for simplicity we assume the same shape at xˆ+1/2). On theright, equivalent mirror-plane symmetric cross-section C′, asrepresented in the 2.5D FDTD. . . . . . . . . . . . . . . . . . 61Figure 4.4 Vocal tract boundary losses due to the (a) 2D contour on themid-sagittal plane (b) tube surface across the mid-sagittal plane. 62Figure 4.5 At a low grid resolution (∆s = 0.74mm) the vocal tract geom-etry for vowel /a/ (a) before scaling (b) after scaling. . . . . . 64Figure 4.6 2.5D tube discretization (mid-sagittal contour and depth map)for vowel /a/. . . . . . . . . . . . . . . . . . . . . . . . . . . 65Figure 4.7 Representing acoustic parameters (p, vx & vy) and depth com-ponents (D¯, D(x) & D(y)) for each grid cell. . . . . . . . . . . 66Figure 4.8 The colormap shows the mid-sagittal view of acoustic pres-sure wave propagation inside the 2.5D vocal tract tube at threedifferent time instants for vowel /a/. A sinusoidal excitationfunction is used at the glottal end. The x and y axis of theabove plots represent the number of grid cells. . . . . . . . . . 67Figure 5.1 General outline of an articulatory speech synthesizer. . . . . . 70xviiFigure 5.2 Highlighted steps (yellow) in the modified roadmap which areaddressed in this chapter. . . . . . . . . . . . . . . . . . . . 71Figure 5.3 (a) Vowel quadrilateral representation and (b) Tabular view ofEnglish diphthongs as per the IPA (International Phonetic Al-phabet) chart. . . . . . . . . . . . . . . . . . . . . . . . . . . 73Figure 5.4 (left) Input band-pass velocity pulse (right) and its frequencycharacteristics. In the frequency domain the spectral amplituderemains constant from 2Hz-22kHz. . . . . . . . . . . . . . . 78Figure 5.5 Snap shot of acoustic wave propagation at t = 15ms for twodifferent mouth-end boundary conditions (a) open-end bound-ary condition (b) free field radiation condition. . . . . . . . . 78Figure 5.6 Vocal tract shape for vowel /A/ with open-end boundary con-dition. The black dot at the right hand side represents the po-sition of virtual microphone/listener (3mm inside the mouthaperture). . . . . . . . . . . . . . . . . . . . . . . . . . . . . 81Figure 5.7 Vocal tract shape for vowel /A/ with radiation condition. Theright hand side close-up image shows the detail of 2D schemeat the mouth aperture. The black dot in it represents the po-sition of virtual microphone/listener (3mm outside the mouthaperture). . . . . . . . . . . . . . . . . . . . . . . . . . . . . 81Figure 5.8 (a) Illustration of a one-end closed cylindrical tube with fun-damental mode (b) Tabular representation of cylindrical tubeswith different lengths of Li (in cm) and cross-sectional areas ofAi(in cm2). . . . . . . . . . . . . . . . . . . . . . . . . . . . 82Figure 5.9 Formant deviations ∆Fi/Fi (%) of the 2.5Dlow simulations ofone-end closed cylindrical tube models (L1 = 17.07cm) withrespect to their analytical solution. The percentage of devia-tions have been derived using the formant positions from Ta-ble C.1. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 84xviiiFigure 5.10 Formant deviations ∆Fi/Fi (%) of the 2.5Dlow simulations ofone-end closed cylindrical tube models (L1 = 50.16cm) withrespect to their analytical solution. The percentage of devia-tions have been derived using the formant positions from Ta-ble C.2. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 84Figure 5.11 Formant deviations ∆Fi/Fi (%) of the 1D and 2.5Dlow simula-tion of a one-end closed cylindrical tube model (L1 = 17.07cm,A1 = 0.91cm2) with respect to its analytical solution. The per-centage of deviations have been derived using the formant po-sitions from Table C.3. . . . . . . . . . . . . . . . . . . . . . 87Figure 5.12 Formant deviations ∆Fi/Fi (%) of the 2.5D simulations withdifferent grid resolutions for vowel /A/. The percentage of de-viations have been derived using the formant positions fromTable C.4. . . . . . . . . . . . . . . . . . . . . . . . . . . . . 88Figure 5.13 Formant deviations ∆Fi/Fi (%) of the 2.5D simulations withdifferent grid resolutions for vowel /i/. The percentage of de-viations have been derived using the formant positions fromTable C.5. . . . . . . . . . . . . . . . . . . . . . . . . . . . . 89Figure 5.14 Formant deviations ∆Fi/Fi (%) of the 2.5D simulations withdifferent grid resolutions for vowel /u/. The percentage of de-viations have been derived using the formant positions fromTable C.6. . . . . . . . . . . . . . . . . . . . . . . . . . . . . 90Figure 5.15 Formant deviations ∆Fi/Fi (%) of the 2.5D simulations withdifferent grid resolutions for vowel /æ/. The percentage of de-viations have been derived using the formant positions fromTable C.7. . . . . . . . . . . . . . . . . . . . . . . . . . . . . 91Figure 5.16 Vocal tract transfer functions (VTTFs) for vowel /A/ generatedfrom the computational simulation of (red dotted line) Kee’s1D model (blue line) 2.5D FDTD model. . . . . . . . . . . . 98xixFigure 5.17 Transformed vocal tract area function (vowel sound /A/) for the2D and 2.5D vocal tract model with reference to the originalarea function (3D). The highlighted area inside the (dotted)blue box shows area functions of different vocal tract modelsat a narrow constriction. . . . . . . . . . . . . . . . . . . . . 98Figure 5.18 Vocal tract transfer functions for vowel sounds (/A/, /i/, /u/ and/æ/) generated from the computational simulation of (1) 2Dhigh(red dotted lines) (2) 2.5Dlow (blue lines). . . . . . . . . . . . 100Figure 5.19 The transfer function of vowel sound /A/ for the (a) 2.5D non-radiation [blue line] and (b) 2.5 radiation vocal tract model[red dotted line]. We run both of these simulations at a fixedboundary admittance coefficient of µ = 0.010. . . . . . . . . 102Figure 5.20 Bandwidth variations (BW1) of the first formant F1 for the vowel/A/ by the 2.5D FDTD simulation using the open-end boundarycondition. . . . . . . . . . . . . . . . . . . . . . . . . . . . . 104Figure 5.21 Coronal (left) and superior (right) view of the two mass ap-proximation of the vocal cords [61]. . . . . . . . . . . . . . . 105Figure 5.22 Normalized output pressure for the 2.5D FDTD simulation ofvowel sound /A/ with a (left) uncoupled VF model (right) cou-pled VF model. . . . . . . . . . . . . . . . . . . . . . . . . . 109Figure 5.23 Glottal waveform (normalized volume velocity - ug) of a two-mass vocal fold model for vowel simulation /A/ with (left) un-coupled system (right) coupled system. . . . . . . . . . . . . 109Figure 5.24 Mid-sagittal view of acoustic pressure wave propagation insidean uncoupled 2.5D vocal tract tube at three different time in-stants for vowel /a/ simulation. The two-mass vocal fold modelis positioned at the glottal-end. The x and y axis of the aboveplots represent the number of grid cells. . . . . . . . . . . . . 110Figure 5.25 Mid-sagittal view of acoustic pressure wave propagation insidea coupled 2.5D vocal tract tube at three different time instantsfor vowel /a/ simulation. The two-mass vocal fold model ispositioned at the glottal-end. The x and y axis of the aboveplots represent the number of grid cells. . . . . . . . . . . . . 111xxFigure A.1 As acoustic wave displaced a distance of ξ , the air particles onplane AB move to A’B’ and on plane CD move to C’D’ . . . . 137Figure B.1 An illustration of vocal tract contour in a 2D rectangular grid.Number of cells between the upper and lower contour is alwaysodd. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 141Figure B.2 Mid-sagittal cross-sectional area function against vocal tractlength and the corresponding vocal tract structural models forvowels: [A, i, I, U, u, o] . . . . . . . . . . . . . . . . . . . . . 144xxiList of Abbreviations1D: One-dimensional2D: Two-dimensional3D: Three-dimensionalBC/BCs: Boundary condition/ Boundary conditionsCFL: Courant-Friedrichs-Lewy conditionCSA: Cross-sectional areaCT: Computed tomographyMRI: Magnetic Resonance ImagingrtMRI: real-time Magnetic Resonance ImagingEMA: Electromagnetic articulographyFEM: Finite Element Model/ Finite Element MethodFDTD: Finite-Difference Time-Domain methodFFT: Fast Fourier TransformPa: Pascal (unit of pressure)Pr: Pressure componentUA: Upper airwayPDE: Partial differential equationDOF: Degree of freedomSPG: Semi-polar gridGPUs: Graphics Processing UnitsVTTFs: Vocal Tract Transfer FunctionsxxiiAcknowledgmentsI express my sincere gratitude to my supervisor, Dr. Sidney Fels; this thesis wouldhave remained incomplete without his continuous support, guidance and mentor-ship. From the outset, his belief in my ability and providing opportunity with allthe academic freedom are aspects of my Master’s that I will always treasure. Thelessons that I learned about cherishing small moments in life and having a sense ofhumour in stressful times, I will take forward in my career.A heartfelt thanks to Dr. Victor Zappi for his endless guidance and persever-ance. Victor’s extensive expertise in physical modelling and computational acous-tic; and the generosity with which he helped me were significant factors in mymotivation and accomplishments. I am deeply indebted to you for always takingout your valuable time to answer my naive questions and detailed discussions overnumerous phone calls and emails. I hope someday I will meet you in person toexpress my gratitude.I would also like to thank all my colleagues in the Brain2Speech project andHCT lab group. Many of you have inspired and enlightened me throughout myresearch. I will never forget the good moments that we all have shared duringour weekly meetings and outings. A special thanks to the x427 lab members forhelping me in my difficult times and always being kind to me.My appreciation goes to my dear friends who are staying miles away but neverlet me feel the distance. I am also thankful to all my roommates who helped meduring my stay in Canada. Finally, I would like to acknowledge that I have beenfortunate to have my parents and my sister supporting and encouraging me through-out the journey. I would not have come so far without your patience and prayersfor me. Thank you for everything!xxiiiDedicationTo my parents, for your perseverance and endless supportTo my sister, for being there with mexxivChapter 1IntroductionThe human voice is one of the most unique and natural ways of communication.From speaking to whispering, grunting to gabbing or crying to laughing, we exer-cise voice to express our thoughts, emotions, and individuality. It is a quintessentialmedium through which we share knowledge and pass information, making humanexistence possible. Through the years, significant efforts have been made to under-stand the underlying mechanism behind how humans perceive and produce speech.Though both of these processes are interrelated, speech perception has led us tobuild computer systems that can recognize speech, whereas studying speech pro-duction has motivated us to build computational models that can generate artificialvoice. This thesis exclusively focuses on the latter part and introduces a uniquephysics-based computational model to design a real-time articulatory speech syn-thesis.1.1 Speech Synthesis TechnologyNumerous speech synthesizers exist to generate artificial voices. The ultimate goalof all these models is to improve the naturalness and intelligibility of syntheticspeech, which can resemble the human voice. Depending on their fundamentalprinciples, they can be categorized as follows:• Formant synthesis• Time-domain concatenative synthesis1• Articulatory synthesisFormant synthesis, also known as rule-based synthesis, generates the artificialspeech waveform from an acoustic model by controlling various adjustable voiceparameters [30]. Typically for a vocal fold model, these parameters are the funda-mental glottal frequency (F0), degrees of voicing and noise levels. For a vocal tractmodel, they are formant (resonance) frequencies and their respective bandwidths.The formant synthesis can produce reliably intelligible speech due to its flexibilityand the high degree of controls. A wide variety of emotional speech, prosodiesand intonations are generated [87, 95] using this technique. However, the speechoutputs from a formant synthesizer are not natural, rather robotic-sounding speech.The introduction of PSOLA (Pitch-Synchronous Overlap-Add) in 1985 moti-vated the research and development of the time-domain concatenative synthesissystems. These synthesizers depend upon the concatenation of segments of prere-corded speech at run-time. These segments such as phonemes, di-phones, sylla-bles, words or sentences are also called speech units. The time-domain concate-native synthesizers are the most natural-sounding as they output artificial voicesgenerated from the actual human speech samples. The various subtypes of con-catenative synthesis [60, 92, 134] are unit section synthesis, diphone synthesis anddomain-specific synthesis. However, there is a large discrepancy between naturalvariations in our voice while speaking and the automated techniques for segment-ing and concatenating speech units. This technique might lead to audible glitchesin the synthesized audio outputs. We can overcome this limitation by increasingthe speech unit length at the expense of computer storage memory. Moreover, thequality of the generated speech output from the time-domain concatenative synthe-sizers relies upon a database that can accommodate a broad set of speech units thatrequire huge inventory space.Compared with other speech synthesis methods, the articulatory speech syn-thesizer is the most sophisticated and computationally expensive approach. Intheory, it could yield high-quality synthetic speech utterances. The fundamentalprinciple of an articulatory synthesis is to build a computational model that closelyfollows the speech anatomy structure and emulates the human speech productionmechanism. Hence, the potential of such methods [99, 101] rely on the accurate2representation of articulatory and acoustic components involved in the speech pro-duction process. Unlike formant synthesis, we can not produce any sound usingarticulatory synthesis due to the limitation of human speech physiology, and thatbrings naturalness to the synthetic sound. Similarly, instead of prerecorded voicedata, the output of an articulatory synthesizer depends only on the model simu-lation, which increases the voice intelligibility. However, one of the significantchallenges with this approach is to generate speech-like sound in real-time with-out compromising the acoustic features of the model-generated speech waveforms.Over the years, many low-dimensional articulatory models are proposed to reducethe computational cost; still, most of them are not flexible enough to adopt variousphysiological features associated with the human speech production mechanism.1.2 MotivationThe explosion of computing capabilities in the past decade and substantial researchin numerical mathematics have created an opportunity to design high-dimensionalspeech production models. Moreover, the emerging high-quality technologies inthe medical imaging field [40] have opened the door to capture the real-time vocalfold oscillation, the intricate geometry of a vocal tract and its dynamic bound-ary variation while speaking. These data have the utmost significance and madethe path easy to study several physical phenomena such as the collision and self-oscillatory characteristics of vocal folds, the nonlinear coupling between vocalfolds and vocal tract, generation of turbulent airflows, the transformation of aerody-namic energy into acoustic energy and vocal tract articulations to produce speechsounds. Therefore, articulatory speech synthesis is once again become an activearea of research.Besides offering deep insights into various guiding principles of speech pro-duction and the characteristic differences in the human voice, such physics-basedspeech synthesis models have wide-ranging applications such as graphical talking-head model as a virtual assistant, speaker-specific speech synthesizer, real-timesinging synthesizer, etc. In a clinical environment [93, 122], they can also be usedto examine the impact of inborn defects in speech anatomy and evaluate outcomesof post-surgical interventions, which will open new possibilities for the diagnosis3and amelioration of various speech pathologies. Finally, better models that imi-tate the physiology of the human speech system can contribute to language peda-gogy [16, 42] (e.g., visualizing the articulatory movements and how they influenceacoustic wave propagation and speech output).1.3 Articulatory Speech Synthesis ModelAn articulatory model consists of several components representing different speechanatomies. Though each component has its significance on the output speech wave-forms, the major specialized areas of study in the field of articulatory speech syn-thesis are vocal fold modelling, vocal tract modelling and their coupling mecha-nism. Several vocal fold models (e.g., lumped-element and continuum models)already exist in the literature [61, 132], which can precisely emulate the structuraland flow characteristics of the human vocal fold. Nevertheless, these models stillsuffer from various computational limitations as well as being challenging to deriveparameters from the actual vocal fold data. On the other hand, the existing low-dimensional vocal tract models are now tractable on GPUs for real-time or quasi-real-time performance. However, they also face obstacles while directly using theactual human vocal tract geometry as input parameters for the computational sim-ulation. Therefore, in this thesis, we primarily focus on a novel low-dimensionalvocal tract modelling approach that directly uses the actual vocal tract geometryand illustrate an end-to-end fully connected articulatory speech synthesizer.The vocal tract model requires much focus as it filters out the glottal acousticenergy to produce speech sounds. Figure 1.1 presents a flowchart representationof the modelling roadmap to build a computational vocal tract model using therealistic vocal tract geometry extracted from MRIs and evaluate the synthesizedspeech outputs. Primarily, our roadmap highlights several features of an actualvocal tract (i.e., in terms of modelling) in a stepwise fashion that need to be ad-dressed while designing a precise structural and acoustic computational vocal tractmodel. Throughout this thesis, we also use this roadmap to discuss many existingchallenges in the vocal tract modelling research domain. It mainly consists of twosegments known as the vocal tract structural and vocal tract acoustic models. How-ever, we should not look at these segments independently, as the structural varia-4tions (e.g., side branches, curvature details, dynamic walls, etc.) in a vocal tractregulate the internal acoustic wave propagation, consequently influence acousticproperties in the synthesized speech outputs. As the human vocal tract geometry isintricate, our roadmap provides a detailed internal nomenclature of these segmentsfor two different situations- the stationary (static wall) and dynamic (moving wall)boundary conditions. A vocal tract model should adopt these features to synthesizerealistic speech utterances.Primarily, this study explores the mapping of an actual vocal tract shape andits geometrical characteristics to acoustic (i.e., synthesized speech sounds) throughnumerical analysis. Multiple approaches exist in the literature to address this re-search problem. They can be categorized based on the dimension of the numer-ical domain and acoustic wave propagation. For example, sophisticated high-dimensional (3D) vocal tract models have high geometrical fidelity and matchthe realistic vocal tract geometry while producing precise acoustic simulations ofsynthesized speech. However, such models are not pragmatic for designing real-time speech synthesizers due to their high computational cost. Similarly, the low-dimensional models (1D) have fast acoustic simulation, but they do not adopt thehigh-quality spectral features of a 3D model due to their oversimplified represen-tation of the tube geometry. Moreover, it is harder to model boundary losses in1D vocal tracts [131], as we have in realistic 3D models. Some techniques exist toenhance the acoustic features of 1D models by augmenting the characteristics of acylindrical tube and the effects of additional geometric features like curvature andbranching (e.g., piriform fossae, subglottal tract) [57, 82]. Still, they are ineffec-tive in capturing the complete wave interactions of a 3D tube, especially nonplanarwaves. As a result, the acoustic formants of 1D models are generally precise upto the lower end of the speech spectrum, i.e. 5 kHz, while missing the higher-order modes and antiresonances. These missing features significantly influence thenaturalness of synthesized speech sounds [12]. Therefore, there is no direct path-way from the actual human vocal tract geometry to the 1D model that gives goodperceptual agreement of synthesized speech utterances even for simple open staticvocal tract configurations.Alternatively, 2D models provide a more excellent balance of vocal tract geo-metrical fidelity, and the 2D acoustic wave solver can be accelerated using modern5RecordedSpeechGeometry Extraction2.5D Model1D Model2D Model3D ModelVocal Tract Structural ModelSS-Depth Contour & depthSS-Curvature Tube CurvatureSS-Branches Side BranchesDS-Diphthongs Dynamic VTSS-Vowel Static VT Shape1D Model2D Model3D ModelVocal Tract Acoustic ModelSynthesized Speech2.5D ModelVocal Fold ModelFT-Coupling CouplingRS-Perception Perceptual EvaluationExisting ModelsExisting ModelsStationary BoundaryDynamic BoundarySA-Loss VT Boundary LossesSA-Vowel Vowel SoundsSA-Radiation Mouth RadiationSA-Analytical Analytical TubeStationary BoundaryDynamic BoundaryDA-Loss VT Boundary LossesDA-Diphthongs DiphthongsDS-Depth Contour & depthDS-Curvature Tube CurvatureDS-Branches Side Branches DA-Radiation Mouth RadiationRecordedAudio RS-Spectrum Spectrum AnalysisSynthesized AudioFigure 1.1: A graphical illustration of vocal tract modelling roadmap to gen-erate and evaluate synthesized speech utterances [11, 106]. The differ-ent steps in the roadmap nomenclature have been categorized using thefollowing abbreviations- SS: Stationary structure related, DS: Dynamicstructure related, SA: Stationary acoustic related, DA: Dynamic acous-tic related, FT: Vocal fold and vocal tract, RS: Recorded and synthesizedspeech. For each step in our roadmap, we compare the proposed 2.5Dvocal tract model with the existing models. We analyze geometry tostudy structural models, whereas, we do spectrum analysis to comparethe acoustic models.6Trachea Tube Area FunctionTrachea PressurePropagationGeometryTrachea ModelVocal Fold Structural ModelVocal Fold Flow ModelGeometryGeometryPressurePPVocal TractAreaFunction2.5D FDTDWave SolverGeometryUPVocal Fold Model Vocal Tract ModelSynthesizedSpeechLungPressureFigure 1.2: General outline of an articulatory speech synthesizer.GPUs. They allow the onset of transverse modes while simulating the acousticwave propagation inside a 2D tube. Theoretically, compared to the 1D, the 2D vo-cal tract model can offer better geometrical flexibility though its full potentials areyet to be studied. Therefore, in this study, we choose the 2D vocal tract model forour articulatory speech synthesizer. However, the existing 2D models go througha nonlinear geometrical transformation [8, 146] to achieve the acoustic behaviourof a higher dimensional vocal tract model, which exterminates their performanceboost. They also ignore the acoustic wave interactions across the 2D sagittal plane.Hence, we propose a new vocal tract modelling technique (2.5D FDTD) that im-proves the existing 2D models while directly using the actual vocal tract geomet-rical parameters from MRIs, and the model also embeds the effect of nonplanarwaves in its wave solver. The proposed model upholds acoustic features of a real-istic 3D vocal tract model without compromising the simulation performance.In this thesis, we then illustrate an articulatory synthesizer that uses our 2.5Dvocal tract model to generate speech utterances for open static vocal tract shapes.Figure 1.2 outlines the architecture of a generic articulatory model which con-tains the vocal tract modelling roadmap (see Figure 1.1) as one of its principalcomponents. The articulatory model consists of the trachea, vocal fold and vo-cal tract model. In this work, we ignore the influence of the trachea model anduse a lumped-element vocal fold model for our speech synthesizer. The following7section provides a walk-through of our modelling roadmap.1.4 Vocal Tract Modelling Road MapIn this section, we refer to Figure 1.1 to outline a step-by-step approach that wehave planned to follow for introducing a novel vocal tract modelling technique(2.5D FDTD Vocal tract) that can generate and evaluate synthesized speech out-puts. As there is an overlap of multiple steps in the roadmap between the station-ary and dynamic boundary conditions, we have not outlined them independently.Instead, we present the steps connected to the stationary vocal tract boundary andonly include the varying steps of the dynamic boundary condition. However, thedynamic boundary condition requires attention to other additional details such asvocal tract geometry extraction from rt-MRIs, a smooth transition of vocal tractwalls in a computational domain, etc. Here, we do not highlight RS-Spectrum step(i.e., transfer function analysis) from our roadmap separately, as it is a standard ap-proach for characterizing the acoustic features of a vocal tract model with variousconfigurations (e.g., different shapes, side branches, lossy and lossless tubes). Es-sentially, it applies to all the below steps having the stationary vocal tract boundarycondition.We also refer to the roadmap to track our proposed model’s current progress,characterize its acoustic output, and set future goals. It is to be noted that thoughevery step in our roadmap is indispensable, their orders are not necessarily fixed.SS-Depth Use mid-sagittal contour and tube depth to capture the effect of non-planar acoustic wave interactions across the mid-sagittal plane inside a 2D compu-tational domain.SA-Analytical Characterize acoustic output of a lossless analytical tube model(e.g., one-end closed cylindrical tube) using the proposed numerical scheme.SA-Loss Incorporate vocal tract boundary losses for the 2D mid-sagittal contourand the tube surface across the sagittal plane (tube depth).8SS/SA-Vowels Characterize acoustic output of the complex open static vocal tractshapes (e.g., vowel sounds) using the 2.5D FDTD vocal tract model.SA-Radiation Characterize acoustic features of the proposed model for the mouthradiation condition.FT-Coupling Implement a nonlinear coupling mechanism to connect a lumped-element vocal fold model to the vocal tract model.SS-Curvature Analyze the geometrical flexibility of the proposed vocal tract model(e.g., tube structure: straight and bent tubes; tube curvature: circular, elliptical andirregular cross-sections etc.).SS-Branches Study the influence of vocal tract side branches with the 2.5 model.RS-Perception Perceptually evaluate the synthesized speech output generated fromthe proposed model.DS/DA-Diphthongs Implement dynamic vocal tract boundary condition to gen-erate diphthong sounds.Our proposed 2.5D vocal tract model is still in progress. Therefore, this the-sis only presents a selective number of steps from our roadmap to demonstrate themodel’s designing procedures and illustrates its basic capabilities by characteriz-ing acoustic features for the stationary boundary condition. In this study, we onlymodel open static vocal tract shapes due to their simplified geometry, and theircross-sectional datasets are already available in the literature. The existing low-dimensional models produce a good approximation of acoustic spectra for thesevocal tract shapes. Thus, our 2.5D FDTD model needs to be characterized forthese simple vocal tract shapes at first. Similarly, we only characterize the 2.5Dmodel using the transfer function analysis method (RS-Spectrum) instead of theperceptual evaluation (RS-Perception) of synthesized audios. Though correct posi-tioning of formants in the transfer functions of the synthesized speech outputs does9not translate to perceptual accuracy, it is still an essential part of the real speechwaveforms. Therefore, these preliminary steps are critical for the future develop-ment of the 2.5D vocal tract model.Hence, the scopes of this thesis are threefold. First, we improve the exist-ing 2D FDTD acoustic wave solver by introducing a new vocal tract model (2.5DFDTD model), which allows capturing the effect of off-plane wave interactions(SS-Depth). However, at the current stage, the model only imposes boundarylosses due to the 2D mid-sagittal contours (partially SA-Loss). Second, we char-acterize the model’s acoustic properties with reference to an analytical tube (SA-Analytical) and open static vocal tract shapes for cardinal vowels (SS/SA-Vowels).We also propose a 2.5D vocal tract radiation model by including a simplified headgeometry inside the computational domain (SA-Radiation). Third, we illustrate anend-to-end articulatory speech synthesizer by coupling our vocal tract model to alumped-element vocal fold model (FT-Coupling).1.5 ContributionMajor Contributions:• Developed the 2.5D FDTD vocal tract acoustic model. We derived theacoustic wave solver equation for the 2.5D vocal tract model. The modelextends the existing 2D wave solver by adding depth (i.e., tube height acrossthe mid-sagittal plane) to the numerical domain and facilitating the nonpla-nar acoustic wave interactions. We demonstrated a strategy for extracting thevocal tract depth parameters from its cross-sectional area function and thenshowcased how to map them to their corresponding acoustic components ina 2D domain.• Characterized acoustic features of the synthesized speech sounds gener-ated from the 2.5D vocal tract model. In the context of acoustic simula-tion, we first presented the acoustic features of our proposed model matchwith the analytical solution of a 1D one-end closed lossless cylindrical tube.We then carried out transfer function analysis for the 2.5D model for fourEnglish vowel sounds. With this approach, we showed that the 2.5D FDTD10simulation-generated acoustic spectra are comparable to a high-quality 3Dvocal tract model and the existing 2D and 1D models.Secondary Contribution:• Designed the 2.5D vocal tract radiation model and examined its acousticcharacteristics. We designed a circular baffle (a simplified head geometry)around the 2.5D vocal tract model and then added PMLs at the computationaldomain boundaries to illustrate the vocal tract radiation model. The PMLsemulate an infinite space for the outward propagating pressure waves at themouth aperture. We compared the acoustic formants of our radiation modelwith the 3D FEM radiation model and the recorded speech signals that areavailable in the literature.1.6 Thesis OrganisationThis thesis entirely focuses on the 2.5D FDTD vocal tract modelling technique andthe computational simulation of open static vowels using the proposed model asthey are the essential step for the model characterization and future works. Follow-ing are the thesis chapters:• Chapter 2 presents a brief introduction to the physiology of human phona-tion with a detailed literature review of the articulatory voice models.• Chapter 3 explores the limitations of the existing 2D vocal tract modelsusing the FDTD numerical scheme.• Chapter 4 formulates the acoustic wave solver PDEs and introduces a newmethodology for the 2.5D FDTD vocal tract acoustic modelling.• Chapter 5 characterizes acoustic features of the 2.5D model by compar-ing its acoustic spectra to other vocal tract models in the literature. It alsodemonstrates the coupling of the 2.5D vocal tract model to a lumped-elementvocal fold model to design an articulatory speech synthesizer.• Chapter 6 summarizes our work, describes potential future directions andprovides concluding remarks.11Chapter 2Background and Previous WorkThe detailed study of speech science began in the eighteenth century when Kratzen-stein and Von Kempelen demonstrated the voicing of vowels and consonants inspeech by using an artificial speaking machine1. By the mid-nineteenth century,Müller introduced the source-filter theory [86], which still guides modern phonet-ics research in many cases. The source-filter theory considers speech productionas a two-stage process, which involves the periodic oscillation of vocal folds asthe source of acoustic energy having its own spectral shape. As an independentcomponent, the resonant property of a vocal tract (filter) amplifies or attenuatesthe source spectrum to form an output speech waveform. The complex vocal tractgeometry and its constrictions [43, 73] govern the resonant characteristics of thefilter. The shape of these constrictions can be determined by muscular contractionsdriven by motor neurons [50]. Later, with the advent of new digital technologies,quantitative studies such as Stevens and House (1955) and Fant (1970) refined thetraditional phonetic theories and provided a new comprehensive framework for thestudy of modern speech science. This chapter offers a brief overview of a few se-lected speech organs and their physiology to describe the contemporary phonetictheories. It also analyzes the existing computational speech synthesis models.Section 2.1 briefly introduces a particular set of speech organs that play a vitalrole in human voicing. It reviews their functionalities and interactions regardingthe physiology of human speech production. Though all of these speech organs1http://www.haskins.yale.edu/featured/heads/SIMULACRA/kempelen.html12Figure 2.1: Basic Anatomy of the Human Speech Production: larynx, sub-glottal component, supraglottal component,nAnatomy & Physiology,OpenStax CNX, adopted from [88].(shown in Figure 2.1) have indispensable contributions toward the speech produc-tion, it will be appropriate to categorize them into three components for furtherdiscussion: the larynx, the subglottal component (below the larynx, e.g. lungs,trachea, and associated respiratory musculature) and the supraglottal component(above the larynx, e.g. upper vocal tract, nasal tract). This approach will serve thepurpose of explaining the speech production mechanism on a high level withoutgoing into great detail about the function of each speech organ.The larynx contains the vocal folds and works as a reference point for the three-way splits. The vocal folds rapidly open and close, converting steady airflow fromthe lungs to a series of quasi-periodic air puffs. These periodic air puffs flow intothe upper vocal tract in the supraglottal region. An articulatory voice model followsthis mechanism by building independent computational models for the vocal folds,upper vocal tract and connecting them through a nonlinear coupling approach. Sec-tion 2.2 reviews the existing state-of-the-art articulatory models with a particularemphasis on vocal tract modelling, which is the essence of this thesis. Section 2.313Figure 2.2: Inferior view of the vocal folds from the laryngopharynx, nAnatomy & Physiology, OpenStax CNX, adopted from [88].demonstrates the implementation of the finite-difference time-domain numericalscheme to the generic 2D acoustic wave equation and the process of discretizingthem for further numerical analysis. Section 2.4 shows the construction of vocaltract structural models from cross-sectional area data sets of various vowel sounds.2.1 Mechanics of Voice ProductionThe human vocal system starts with the lungs and lower airway (trachea), whichwork like a pair of bellows. The lungs supply air pressure and continuous air-flow through the trachea. This functionality can be described using two principles:Boyle’s law, which states the air pressure and volume are inversely proportional ata constant temperature, and the equalization of pressure, which governs the flowof air due to the pressure difference. During respiration, the chest cavity or thoraxincreases or decreases the pressure on the lungs, causing the depletion or expan-sion of the lungs’ volume. There is continuous airflow from the lungs due to thesevolumetric changes.The larynx, a cartilaginous structure inferior to the laryngopharynx, is locatedat the end of the trachea. It connects the trachea to the pharynx. The vocal folds(a.k.a., ‘vocal cords’), located in the larynx (shown in Figure 2.2), form a constric-14tion to the airway as it stretches along the posterior and anterior direction across thelarynx. This movement regulates the continuous airflow from the lungs. The vocalfolds are connected by muscle to the thyroid and arytenoid cartilages of the larynxand are often represented as two-layer body-cover structures. The body layer con-sists of the muscular and deep layers of the lamina propria, whereas the cover layerincludes the intermediate and superficial lamina propria and the epithelium layer.The space between the two vocal cords is known as glottis.The phonation starts as the periodic oscillation of vocal folds modulates thecontinuous airflow from the lungs to produce acoustic energy. Van den Bergfirst proposed the classic myoelastic-aerodynamic theory [130] to explain the self-oscillatory behaviour of the vocal folds during phonation. His theory was moti-vated by two concepts: the myoelastic nature of vocal folds and the airflow aero-dynamics through the glottis. According to the myoelastic-aerodynamic theoryof voice production, the complete adduction of the vocal folds closes the glottis.However, due to the lungs’ continuous airflow, air pressure (sub-glottal pressure)builds up at the bottom surface of the vocal folds. When the sub-glottal pressuresurpasses a particular threshold value, known as the onset pressure, it force opensthe glottal walls and allows the air to escape through the glottis. This phenomenoncan be explained using Bernoulli’s principle of fluid dynamics and the laws ofconservation of mass and energy through a computational vocal fold model. Themodel demonstrates the drop in air pressure at the intraglottal and supragolttal re-gions as the air rapidly flows through vocal folds. Due to this negative pressure,the elastic recoil force pulls back the vocal folds and closes the glottis. The cyclethen repeats to maintain the oscillation of vocal folds.Though the traditional myoelastic-aerodynamic theory demonstrates the vo-cal fold oscillation, it does not explain how energy is transferred from the airflowinto the vocal folds to sustain this oscillation. Later, Ishizaka and Matsudaira [61]showed the vertical phase difference in the vocal fold surface could help to achievepressure asymmetry and the energy transfer between the airflow and vocal folds tocreate a sustained self-oscillatory motion. This phenomenon is generally consid-ered the primary mechanism for voice production and is widely adopted in manycomputational models [120, 149] for the physics-based speech synthesizers.152.2 Articulatory Voice ModelsAn articulatory voice model comprises three different computational components,each addressing separate issues: a glottal model as the acoustic source, a vocal tractstructural model to define the intricate boundaries for various vocal tract shapes anda vocal tract acoustic model to emulate the pressure wave propagation.• Glottal Source Model: The main objective here is to match the self-oscillatorybehaviour of vocal folds while simulating their flow characteristics. Glottalsource models can be classified into parametric [44], kinematic [121] andself-oscillating biomechanical models. However, we only discuss the exist-ing biomechanical vocal fold models as they can simulate the flow-inducedoscillation.• Vocal Tract Structural Model: This model creates a sophisticated geo-metrical boundary corresponding to a vocal tract shape that changes at eachinstant of time depending upon the steadiness of voiced sounds.• Vocal Tract Acoustic Model: Once the source model induces the acousticenergy into the upper vocal tract, the acoustic model approximates the acous-tic pressure wave propagation inside a predefined vocal tract geometry. Themodel measures and stores the variations in pressure wave amplitudes in anaudio file at the mouth end to produce the synthetic speech sound generatedfrom the model simulations.2.2.1 Glottal Source ModelsThe biomechanical glottal models are a combination of structural and flow models.The structural model illustrates the glottal opening and closing and the change inthe glottal area during vocalization. The flow model describes airflow through theglottis and the induction of energy to the vocal fold structure to maintain its self-oscillatory behaviour. Two main types of glottal models exist in literature—first,the lumped-element models, which consider vocal folds as point masses. Second,the continuum models, which use numerical meshes to describe vocal fold surfaces.16A. Lumped-Element Vocal Fold ModelsThe lumped-element vocal fold models [20, 81] are the most simplistic modelswith a low computational cost. They approximate the vocal fold tissues to lumpedmasses connected to a rigid body through spring and damping elements. The springand damping elements intend to account for the tissue elasticity and energy losses.Despite their simplicity, many of these models are capable enough to reproduceoriginal human vocal folds characteristics [70]. The earliest lumped-element modelwas introduced in 1968 by Flanagan and Landgraf [45]. This model, also knownas the one-mass vocal fold model (shown in Figure 2.3a), represents vocal foldsas a single linear spring-mass system. With single DOF, a one-mass model onlysimulates lateral displacement of vocal folds but fails to introduce the phase differ-ence between upper and lower vocal fold edges during oscillation [107]. Hence, toachieve a sustaining self-oscillation motion, the model needs to be coupled with avocal tract while the supraglottal pressure and vocal fold tissue velocity need to bein phase [145]. Later, Adachi and Yu [1] proposed a new one-mass glottal modelhaving two DOF (shown in Figure 2.3b), which can oscillate both in parallel andperpendicular to the airflow. Unlike Flanagan’s model, this vocal fold model cansustain self-oscillation under capacitive acoustic load. In literature, there are othervariants [14, 97] of the one-mass vocal fold model exist. Nevertheless, they are nolonger in use due to their over-simplified representation of vocal folds and manycomputational limitations.The multi-mass glottal models use more than one mass element to representeach side of the vocal fold. Ishizaka and Flanagan [62, 96] first proposed the clas-sic two-mass vocal fold model (shown in Figure 2.4a), which was able to duplicateprincipal features of human vocal cords. The model successfully mimics the phasedifference along the vocal cord edges and produces perceptually better speech syn-thesized sounds than the one-mass model. As a low-dimensional vocal fold model,the two-mass model has high relevance for articulatory voice modelling. Tokudaet al. [124] added one more mass on each side of the vocal fold to simulate themodal voice and falsetto. With an additional vertical DOF, Koizumi et al. [70] pro-posed a new two-mass vocal fold model, which was more natural-sounding. Themulti-mass models successfully discretize each side of the vocal fold by adding17(a) with single DOF (b) with two DOFFigure 2.3: Variants of one-mass vocal fold models, reproduced with permis-sion from Peter Birkholtz [20].(a) two-mass model (b) body-cover modelFigure 2.4: Coronal view of the two-mass (left) and body-cover (right) vocalfold models, reproduced with permission from Peter Birkholz [20].18more mass elements. These models capture the oscillatory motion across vocalfold edges in detail. However, they do not distinguish between the vocal fold bodyand cover tissues. Therefore, Story and Titze [110] introduced a three-mass model(shown in Figure 2.4b), which follows the layered body-cover structure of the vo-cal folds. Story’s model uses two spring-mass elements and a single spring mass-element to represent the cover and body tissues, respectively. It is also known asthe body-cover vocal fold model. Zhang [148] demonstrated the phonation onsetand sound production efficiency by controlling the body stiffness in a body-coverglottal model. These lumped-element models use the classic Bernoulli-based flowmodel or 2D Navier-Stroke solvers to characterize the airflow in the vocal fold anddemonstrate its self-oscillatory motion.B. Continuum Vocal Fold ModelsWith an increase in the number of mass elements and DOF, the lumped-elementmodels may provide better results. However, there is a severe challenge in adopt-ing the detailed geometry of the vocal folds and their viscoelastic properties. Thereis a major shift in the vocal folds modelling approach from the last decades, whichincludes continuum models due to a boost in computing capabilities. The highlycomplex continuum vocal folds models take advantage of the FEM approach (shownFigure 2.5: Finite element mesh of vocal fold model shown in coronal plane,adopted from [132]19in Figure 2.5) to showcase the dynamics of complex boundaries and the flow-induced vibration of the structural model. These models [2, 38] are based on thefundamental laws of continuum mechanics. They generally integrate the aeroa-coustic flow and the vocal folds’ structural model by solving their PDE equationsnumerically either in the time or frequency domain.Zheng et al. [150] employed an immersed boundary method flow solver cou-pled with a 2D FEM vocal fold structural model to demonstrate a robust solutionfor the flow-structure interaction (FSI) in the laryngeal. The biomechanical modelinvestigates several flow characteristics like phonation onset, glottal jet dynamics,etc. Nonetheless, in many cases, the 2D models do not match with the naturalphysiological characteristics of phonation. Hence, a 3D model is more beneficialfor laryngeal study. Allipour first introduced a 3D continuum glottal model [2],which uses Bernoulli’s equation as the flow server. Rosa et al. [37] first proposeda continuum model that fully couples the 3D structural dynamics of the vocal foldwith a flow model, described by the steady Navier-Stroke equation. Xue extendedZheng’s work and proposed a 3D vocal fold structural model [140] which imple-ments the geometrical and structural details of a realistic human larynx. Mittalet al. [80] used the FEM laryngeal model to design a tool that can guide duringsurgical intervention. Though computationally expensive, these continuum models[63, 64, 136] are considered as the critical steps towards the study of phonationbehaviour and the motor-controlled vocal sound in humans and other mammals.2.2.2 Vocal Tract Structural ModelThe speech-production task requires coordinating movements of articulators (e.g.,jaw, lips, tongue blade, nasal cavity, palate, etc.), laryngeal, and pharyngeal com-ponents in a unified manner. These positional changes shape the upper vocal tractgeometry (airway), influencing its resonance property and speech sound quality.The vocal tract geometry remains fix for vowels or steady voice (shown in Fig-ure 2.6), whereas it changes while producing diphthongs. Consequently, an artic-ulatory speech synthesizer necessitates a structural model of the vocal tract, whichcan approximate the shape of the upper airway. This section discusses the existingvocal tract structural models in detail.20 Figure 2.6: Illustrations of variations in vocal tract shape for different staticvowel sounds with mid-sagittal vocal tract MRI images : (left) Vowel[A], (middle) Vowel [i], (right) Vowel [u], reproduced with permissionfrom Peter Birkholz [106]A. Biomechanical ModelsThe biomechanical or articulatory physiological model characterizes the phoneticsegments in terms of articulatory targets. So, a particular phonetic segment canbe defined as a set of precise articulatory configurations that articulators try toreach through temporal deformations. The model requires a relational knowledgebetween muscle activation and articulatory movements for simulation. Generally,the biomechanical model uses a 2D or 3D numerical scheme (typically FEM) todesign (shown in Figure 2.7) soft tissues (e.g., tongue muscles and soft palate) andrigid bodies (e.g., jaw, hyoid bone, etc.). It then applies muscle activation to movethe articulators and extract the corresponding upper airway geometry.Perkell [89] first constructed a 2D biomechanical tongue model along the sagit-tal plane. Later, Kakita et al. [66], and Hashimoto et al. [54] implemented a 3Dtongue model using the FEM technique. Still, the biomechanical models createdin these studies only incorporate a part of the upper vocal tract, such as a singletongue model [49] or a single facial model [77, 98]. Hirari et al. [56] developeda 2D physiological model composed of both tongue soft tissues and rigid organssuch as the jaw and hyoid bone. It supported the mechanical principle of the dy-namic balance between forces and moments to interface soft and rigid structures.However, the model simulation was computationally expensive as the structureswere modelled independently. Dang and Honda [36] proposed a new FEM-based21Figure 2.7: (left) An illustration of articulators from the upper vocal tract,adopted from image by Patrick J. Lynch, medical illustratornWikime-dia Commons. (right) A midsagittal view of the biomechanical modelproviding an overview of the upper airway, adopted from [5].partial-3D tongue model in the sagittal plane for a faster simulation and to examinemuscle activation functions from the model dynamics.Recently, a new framework (Functional Reference ANatomical Knowledge- FRANK) has been developed to create a generic 3D biomechanical model ofthe human head and neck, which expedites the study of human motor controlbehaviours in speech production. The framework incorporates several previousanatomical models [4, 29, 31, 103] using a comprehensive biomechanical toolkit,ArtiSynth [76]. ArtiSynth is an open-source toolkit that offers an internal physics-based engine to create biomechanical models and simulate their dynamics. TheFRANK model is comprised of FE models, rigid bodies, tongue muscles, bonesand cartilages. Table 2.1 provides a summary of the model components, includingtheir name and modelling type (e.g., FE model, rigid body).Using the FRANK model, a tongue configuration for a specific phonetic seg-ment can be achieved by activating different tongue muscles (shown in Figure 2.8).Although all these previously discussed physiological models provide upper vocaltract geometry with pinpoint accuracy, they are highly sophisticated and compu-tationally expensive. Besides, these models have many degrees of freedom whichmakes them difficult to control. Therefore, deriving upper airway geometry usingthe biomechanical model is not a pragmatic approach at this point.22Component Name TypeFace FETongue FEJaw, Hyoid, Maxilla RigidSoft-Palate FEPharynx FELarynx FELarynx Cartilages RigidTable 2.1: Summary of FRANK model components: component name andcomponent type (FE models (FE), rigid body (Rigid)), adopted from[119]Figure 2.8: Biomechanical tongue models in ArtiSynth for different phoneticsegments: (left) Vowel [A], (middle) Vowel [i], (right) Vowel [u]Figure 2.9: Vocal tract structural model with a number of articulatory param-eters, reproduced with permission from Peter Birkholz [22]23B. Articulatory ModelsArticulatory modelling [22, 23, 79] is the most flexible and straightforward ap-proach to replicate any speaker’s upper vocal tract shape irrespective of their ageand sex. These models heavily depend on a corpus of data such as volumetricMRI images or EMA to match the model anatomy and articulatory information toa speaker. Unlike biomechanical models, articulatory models control articulatormovements through a specific set of chosen parameters (e.g., tongue height, lip po-sition, jaw movement, etc.), as shown in Figure 2.9. Once parameters’ values areset, the corresponding airway geometry can be extracted. However, these valuesare not predetermined and need to be evaluated empirically. Birkholz et al. [23]used MRI images to adjust the controlling parameters. Engwall [41] employed acombinatorial search to estimate them. Nonetheless, like a biomechanical model,an articulatory model may require several controlling parameters to create an artic-ulatory configuration that can match the real upper vocal tract geometry.C. Geometrical ModelsThe biomechanical and articulatory models consider each articulator in the up-per vocal tract as a separate entity. Whereas the geometrical model approximatesthe profile of a complete vocal tract as a single entity, i.e. upper airway cavity(shown in Figure 2.11(a)). To do so, the geometrical model juxtaposes multiplecylindrical tubes of the same length with varying cross-sections. The length andcross-sections of these tubes depend upon the vocal tract shape. Fant first proposed[43] an oversimplified four-tube model for vowels and a three-tube model for con-sonants sounds. However, for a better approximation of the upper airway cavity,the geometrical model should consist of several tube segments. Usually, these ge-ometrical models use a 1D area function data set extracted from the actual vocaltract to determine the tubes’ cross-sections. These tube cross-sections can be mod-elled in 1D, 2D and 3D domains with circular, elliptical or more realistic shapes (asshown in Figure 2.10) for acoustic analysis. Thus, in geometrical models, acousticcharacteristics are only governed by vocal tract area functions but not the individ-ual articulatory movements. This approach allows better control of an upper vocaltract’s acoustic behaviour. It is to be noted that the vocal tract acoustic model can24Figure 2.10: Straight vocal tract tube structure for vowel [A]: (top) 3D tubewith various cross-sectional shapes (realistic, elliptical and circular),reproduced with permission from Marc Arnela [10], (bottom) 2D tubewith circular cross-section(b)(a)Figure 2.11: (a)A 3D mesh representation of the upper airway and cross-sectional area with red circles. (b) Volumetric vocal tract geometryof vowel [A], intersected by grid planes (blue), reproduced with per-mission from Marc Arnela [11].25not directly use the biomechanical and articulatory models for acoustic simulation.It requires the upper airway geometry to approximate the acoustic wave propaga-tions. Therefore in this thesis, we adopt the actual vocal tract area functions as thegeometrical parameters and investigate an efficient way of mapping them to theircorresponding acoustic properties using the numerical analysis method. This fun-damental step needs to be explored first before implementing the articulatory andbiomechanical models for speech synthesis.In general, the 1D area function [15, 111] describes variations in the cross-sectional area of a vocal tract against its length. It can be represented as (x,y) pairs,where x as the vocal tract length (Ln) from the glottis and y as its correspondingcross-sectional area (An). The measurement of area function starts just above theglottis (k = 0) and end at the lips (k = n).(x,y) = (Ln =n∑k=0Lk,An) (2.1)Several methods exist in the literature to compute the area function. Although themost followed are:• Extracting area function from the vocal tract volumetric imageModern medical imaging techniques [15, 83, 112] (typically MRI is pre-ferred over CT scan to avoid radiation risk [28]) are used to capture the up-per vocal tract geometry while voicing. Next, the method applies an imagesegmentation algorithm to separate the upper airway from the surroundingtissue (e.g., sinus piriformis, vallecula and lips) and constructs a volumetricimage out of image data points. For cross-sectional area function analysis,the model introduces a set of 2D grid planes that intersect the 3D volumet-ric vocal tract image, as shown in Figure 2.11(b). Arnela [11] proposed asemi-polar grid (SPG) planes technique to approximate the vocal tract cross-sectional shapes. The SPG consists of three types of planes - horizontalplanes in the pharynx, polar planes in the velar region and vertical planes inthe oral cavity, that do not change their orientation over time. These planesstart from the glottis and are placed in equal intervals till the mouth opening.However, such planes do not coincide with the planar wavefronts across the26Figure 2.12: An illustration of 2D midsagittal vocal tract contour with center-line (white) and semi-polar grid axis (red) to extract the area function.vocal tube and may not provide better acoustic output. Hence, an adaptivegrid strategy [10, 33] was suggested to extract the cross-sectional tube area.In particular, the adaptive grid planes are perpendicular to the vocal tractmidline while extracting the vocal tract cross-sections. An initial midlinecould be estimated using the SPG and a smoothing and resampling techniqueas proposed here [11]. Story presented an iterative bisection method [111]to determine the vocal tract midline and then implemented an algorithm tomeasure its cross-sectional area [58], that was originally used to determineupper airway geometry and volume with regard to sleep disorders.• Extracting area function from the vocal tract mid-sagittal sectionThe mid-sagittal section of a vocal tract provides important information onarticulatory movements and the shape of the upper airway. It has been atthe root of multiple vocal tract models. The mid-sagittal slice of a 3D vo-cal tract tube represents its anterior-inferior and posterior-superior outlines(vocal tract contours) on a 2D plane (shown in Figure 2.12). Once the 2Dvocal tract contour is determined, the method can extract the area functionby estimating a contour centerline. Otherwise, a mathematical relationshipcan be formulated empirically to compute the cross-sectional area at differ-ent vocal tract coordinates against their distance from the glottal end. Heinz,27and Stevens [55] first hypothesized the cross-sectional area (A) as a functionof (d1), the mid-sagittal distance from the glottis, through "αβ model":A = αdβ1 (2.2)where α and β coefficients have different values depending on the regionof the vocal tract. But there were not adequate answers to determine theirvalues. Denis [17] implemented an optimization algorithm based on gradientdescent technique and a back-propagation network to calculate the valuesof α and β . Still, the value of coefficients may change, depending on thesubject. Sundberg [113] analyzed the geometry of pharynx in four horizontalplanes and proposed a linear relation between the area function and squareof the mid-sagittal distance:A = c1d21 + c2 (2.3)where the coefficients, c1 and c2 depend on the vocal tract region and the sexof a subject. Perrier [90] put forward a new nonlinear relationship to deter-mine the area function from the mid-sagittal profile for the french cardinalvowels.On the other hand, the vocal tract outlines in a 2D mid-sagittal plane candirectly extract the area function. Birkholz [21] superimposed a semi-polargrid and applied a smoothing technique to find the centreline of a vocal tract.Then at every equally-spaced point along the centreline, a perpendicularplane is introduced, intersecting the posterior-superior and anterior-inferioroutlines at different coordinates. The distance between the two intersectionpoints for each plane provides the diameter (d) of a single cross-section,which can be used to compute the area-function (A):A = pi(d/2)2 (2.4)Goldstein [51] proposed a medial-axis transform or Blum transform algo-rithm that represents the vocal tract contour with a set of circles of maximum28radius, such that the circles are tangent to the vocal tract outlines. Then amidline can be estimated by joining the locus of all these circles. Laprie[72] designed a loss energy function, which needs to be minimized for betterestimation of the vocal tract centerline.The use of area function simplifies the acoustic analysis since it assumes planewave propagation through the structural model. However, this assumption limitsthe vocal tract model’s acoustic precision, depending on the maximum diameter ofthe tube.2.2.3 Vocal Tract Acoustic ModelAfter having an acoustic source (glottal model) and a vocal tract structural model,the physics-based speech synthesizer requires an acoustic model to emulate soundwaves propagating across the vocal tract. In general, sound waves propagate throughthe air particles as a series of compression and rarefactions. The vibration in spacedisplaces air particles from their original positions, creating pressure gradients dueto the change in the air particle density. Subsequently, the elastic property of airparticles forces them to restore their original position by adding a velocity compo-nent in the resulting wave. This phenomenon can be described using the governinglinear wave equations as follows,Equation of continuity:∂ p∂ t=−ρc2(∇.v) (2.5)Equation of motion:∂v∂ t=− 1ρ(∇p) (2.6)or, can be expanded to represent in terms of the axial components∂ p∂ t=−ρc2(∂vx∂x+∂vy∂y+∂vz∂ z) (2.7)∂vx∂ t=− 1ρ(∂ p∂x) (2.8)∂vy∂ t=− 1ρ(∂ p∂y) (2.9)29∂vz∂ t=− 1ρ(∂ p∂ z) (2.10)where, the parameters p, v, c and ρ represent the scalar pressure field, particlevelocity, speed of sound and medium density respectively. Using Equation 2.5and Equation 2.6, a generic acoustic wave equation can be derived. The completederivation of 1D acoustic wave equation has been provided in Appendix A.1.The physics-based vocal tract acoustic model employs a wave solver that uses anumerical analysis method to solve the above acoustic wave equations either in thetime or frequency domain to generate vowels and diphthongs sounds. The solvercomputes the time-domain acoustic wave equation either in the irreducible form(acoustic pressure) or in the mixed form (acoustic pressure and particle velocity),or its frequency counterpart, the Helmholtz equation. Matsuzaki [78] analyzed theformants and radiation impedance of Japanese vowels in 3D elliptical tube modelsusing the FEM technique in the frequency domain (FD-FEM). Fleischer et al. [47]study the impact of vocal tract wall impedance on formants’ position and band-width using the FD-FEM model. Kako and Touda [67] used the FEM analysisto solve the Helmholtz equation and find the corresponding complex eigenvaluesthat control the peaks of the frequency response function known as formants. Themodel sets a functional relationship that can determine the variation in complexeigenvalues with respect to vocal tract boundary changes. The other frequencydomain FEM implementations can be found here [53, 84, 147]. Kagawa [65] cre-ated a bend vocal tract model and the head geometry using the boundary elementmethod in the frequency domain (FD-BEM).In computational acoustic, the time-domain analysis is always preferred as itallows the visualization of wave propagation that evolves with time. Moreover,various speech phenomena like changes in vocal tract geometry and vocal foldvibration occur over time as we speak. So, time-domain analysis fits well forstudying the voiced sound mechanism. Svancara [114, 115] studied the effect oftonsillectomy on speech production using the time domain finite element analy-sis (TD-FEM). Vampola [126, 128, 129] used the time-domain FE modelling tocharacterize the vocal tract acoustic for Czech vowels. Arnela et al. [7, 8] havedemonstrated the acoustic properties of the 3D and 2D vocal tract models for staticEnglish vowel sounds /A/, /i/ and /u/ using the FEM method. The acoustic proper-301. Introduction 91.2.4 Vocal tract acoustic modelOnce having the vocal tract geometry and the glottal source model (vocal tract inflow),we can focus on simulating the acoustics of the vocal tract, which will finally yield thegenerated voice. Vocal tract acoustic models can be divided in those that rely on one-dimensional strategies and those that follow a three-dimensional approach.One-dimensional modelsFor historical reasons and thanks to their simplicity, one-dimensional models (aka tubemodels) have become widely known and used. These models approximate the vocal tractgeometry as a finite set of concatenative tubes, each one having constant cross-section(i.e., vocal tract area functions, e.g. Story et al., 1996). Then, assuming plane wavepropagation within the vocal tract, which holds up to about 5 kHz, the acoustic pressureat the exit (voice) is computed. They can be subdivided into the ABCD matrix basedmodels (see Fig. 1.4a), the Digital Waveguides Models (DWM, see Fig. 1.4b), circuitanalogy models (see Fig. 1.4c) and computational models:• Tube models of the ABCD matrix type (e.g., Sondhi and Schroeter, 1987) computethe acoustic transfer function of the vocal tract as the products of the individualtransfer functions for each tube. Then, the voice signal is computed in the frequencydomain using the vocal tract transfer function, and the inflow provided by thesource model. In contrast to time domain models, ACBD matrix based modelsܣ ܤܥ ܦ ܣ ܤܥ ܦ ܲ݋ݑݐܷ݋ݑݐ = ܣ ܤܥ ܦ ܷܲ݅݊݅݊ ܲ ܷ = Acoustic pressure = Volume velocity (a)f = forward signalb = backward signalʏ = delay timef1(t) f1(t-ʏͿf2(t) f2(t- ʏͿL Lb1(t) b1(ƚнʏ)b2(t) b2(t+ ʏͿ(b)= Resistor= Inductor= Capacitor(c)Fig. 1.4: Vocal tract acoustic models based on tubes (1D). (a) ACBD matrix based model,(b) Digital Waveguide Model or KL model and (c) circuit analogy model.Figure 2.13: 1D Vocal tract acoustic tube models. (a) ACBD matrix basedmodel, (b) Digital Waveguide Model (c) circuit analogy model, repro-duced with permission from Marc Arnela [12]ties of these models heavily depend on the accurate representation of the intricategeometry of a vocal tract. The 3D FEM tube models accurately represent the cir-cular, elliptical, and irregular tube cross-sections with straight and bent structures.However, as the number of elements in the FEM tube model increases to definecross-sectional surfaces precisely, the simulation run-time also rises. In most cases,the FEM simulation takes hours to several days to simulate a few seconds of audiooutput. Hence, for the real-time articulatory speech synthesizer, the FEM-basednumerical approach is not favoured.Similarly, in contrast to 3D, 1D acoustic tube models are much simpler andcomputationally less expensive due to dimensionality reduction. Numerous 1Dvocal tract acoustic models have been proposed in the literature. Sondhi andSchroeter [101] proposed a hybrid time-frequency domain model that computesthe vocal tract model’s transfer function by multiplying the transfer function ofeach tube. The model uses an ABCD matrix-type tube (shown in Figure 2.13a),where the matrix elements depend on the tube parameters. The voice signal iscomputed by using the vocal tract transfer function and a two-mass source model.The digital waveguide models, also known as KL models [68] proposed by Kellyand Lochbaum (shown in Figure 2.13b), approximate the vocal tract as a set oftransmission lines. It computes the acoustic wave propagation inside a vocal tracttube by solving the one-dimensional acoustic wave equation. But the 1D waveg-31uide models do not allow a smooth change of the vocal tract geometry, which isrequired during the generation of diphthong sounds. Birkholz [24] proposed atube model based on circuit analogies (shown in Figure 2.13c), where the voltageis equivalent to pressure, and the current is equivalent to volume velocity in theacoustic wave equation. The tube parameters like the inertance of the mass of airin the tube section, air compressibility and vocal tract boundary losses are repre-sented as Li, Ci and Ri, respectively. The model computes voice signals in the timedomain. Compared to KL models, circuit models can approximate the acousticwave propagation more precisely as the vocal tract geometry changes [71].However, 1D acoustic models oversimplify geometrical features of an actualvocal tract by assuming plane wave propagation in the vocal tract model. The 2Dacoustic model preserves much of the geometrical features as it allows capturingboth fundamental and transverse modes, improving the acoustic output of the vo-cal tract model. Additionally, 2D acoustic models computationally less expensivethan 3D models. With regard to numerical schemes, the finite-difference time-domain (FDTD) method computes the acoustic wave propagation faster than theFEM and BEM numerical approach [27, 138]. The FDTD technique uses an or-thogonal Yee grid system that can directly load voxel data as an analysis model,whereas the FEM and BEM models need to build the mesh from the voxel datato define elements. Previously, Yokota et al. [143] used the FDTD techniqueto visualize the acoustic propagation in a closed room. Webb and Bilbao [139]demonstrated the simulation of realistic room acoustic by implementing 3D FDTDnumerical scheme. The model uses GPU architecture with parallel computation forbetter computing performance. Yuan and Berggren [144] carried out the acousticfield analysis using FDTD simulation. Recently, the FDTD method is extensivelyused in measuring the acoustic properties of musical instruments and building theircomputational models [3, 18, 19]. However, the regular orthogonal grid in FDTDmay not conform to the irregular geometry neatly all the time. That may impactthe boundary conditions and the corresponding acoustic output. Moreover, it ischallenging to include frequency-dependent losses using FDTD. Despite its limita-tions, the FDTD is the most simplistic approach for acoustic analysis. Thus, in thisthesis, we implement the FDTD scheme to solve the acoustic wave solver equationinside a 2D numerical domain for open static vocal tract shapes. The following32p(x, y) vx(x+1/2, y)vy(x, y +1/2)Figure 2.14: Representing pressure and velocity components in a staggeredYee-grid. The acoustic pressure p(x, y) is sampled at the center of thecell (red circle), vx is sampled at the right edge of the cell (orange) andvy is sampled at the top edge of the cell (blue)section demonstrates the implementation of a generic FDTD acoustic wave solver.2.3 Discretizing Two-Dimensional Acoustic WaveEquation Using The FDTD MethodThe FDTD algorithm is widely used in the electromagnetics domain to solve thedifferential time-domain Maxwell’s equations [142], where both the electric andmagnetic fields are vector quantities. In acoustic, the pressure is a scalar field, andvolume velocities are the vector quantities. The spatial rendering of the acous-tic signal using FDTD can be obtained by discretizing the pressure and velocitycomponents of the acoustic equations both in space and time. Hence, the acousticpressure and velocity components are represented on a rectangular staggered grid,as in the Yee algorithm [141]. The spatial positions of these discretized acousticcomponents have a specific arrangement on the staggered Yee grid. As shown inFigure 2.14, each grid cell has three acoustic components, having pressure values psampled at the center of the cells, while the vx and vy of the velocity vector sampledon the edges shared with the right and top neighbour cell respectively. Hence, in ageneric cell at a discrete coordinate (x, y), p will be sampled at (x, y), while vx at33(x+1/2, y) and vy at (x, y+1/2).The FDTD algorithm employs the second-order central difference approxima-tion (as shown in Appendix A.2) to solve the acoustic wave equations (2.5 and2.6). The wave solver primarily estimates the pressure and velocity componentsat each grid cell inside a 2D Yee grid composed of (Nx×Ny) cells 2. In a com-putational domain3, the grid contains uniform cells with a spatial dimension ∆x inthe x-direction and ∆y in the y-direction. In general, the spatial resolution needsto be uniform along both the axial direction (i.e., ∆x=∆y=δ )4, and the algorithmcomputes the acoustic components at discrete time instants. Since we use the ex-plicit time integration, the time step (∆t) and the cell width (δ ) follow the relatedCourant–Friedrichs–Lewy stability condition (CFL) in two dimensions as follows,∆t ≤ δ√c2(2.11)where, c = speed of sound.By considering the simulation time step = ∆t and grid resolution (∆x, ∆y) = δ , thediscretization of acoustic equation in two dimensions yields:Discretizing Equation 2.7,p(x,y)t+∆t − p(x,y)t∆t=−ρc2(v(x+δ/2,y)t − v(x−δ/2,y)tδ+v(x,y+δ/2)t − v(x,y−δ/2)tδ)By rearranging the above equation yields,p(x,y)t+∆t = p(x,y)t − ∆tρc2δ(cxvx + cyvy) (2.12)where,cxvx = v(x+δ/2,y)t − v(x−δ/2,y)t2In this thesis, we also refer grid cells as grid points and grid pixels.3Here we consider 2D domain only, but the discretization of acoustic wave equations can beextended to 3D in a similar manner.4In the following chapters, we represent spatial resolution δ as ∆s34cyvy = v(x,y+δ/2)t − v(x,y−δ/2)tDiscretizing Equation 2.8,v(x+δ/2,y)t+∆t − v(x+δ/2,y)t∆t= − 1ρ(p(x+δ ,y)t+∆t − p(x,y)t+∆tδ)By rearranging the above equation yields,v(x+δ/2,y)t+∆t = v(x+δ/2,y)t − ∆tρδ cx p (2.13)where,cx p = p(x+δ ,y)t+∆t − p(x,y)t+∆tSimilarly, discretizing Equation 2.9 and rearranging it obtains,v(x,y+δ/2)t+∆t = v(x,y+δ/2)t − ∆tρδ cy p (2.14)where,cy p = p(x,y+δ )t+∆t − p(x,y)t+∆tThe FDTD wave solver (Equations 2.12, 2.13 and 2.14) follows a time-marchingalgorithm as demonstrated in Figure 2.15 to update the acoustic pressure and ve-locity components at each time step. The wave solver first sets up the grid size andother physical parameters of the modelling environment, such as air density, soundspeed, and spatiotemporal resolution. It also initializes the acoustic components(i.e., p, vx and vy) to zero across the grid. Then, for a faster acoustic simulation,the coefficient terms appearing in FDTD update equations (e.g., ∆tρc2δ and∆tρδ ) arecomputed and stored in arrays before the simulation begins. The 2D wave solvermainly solves the pressure and velocity components for all the grid points in anumerical domain during the simulation. At each time step (n), the pressure com-ponent is updated for the next time instant (n+1)∆t using Equation 2.12; and the35StartDefine modelling parametersCompute field coefficientsStart time step (n=0)Update p using 𝒗𝒙𝒏 and 𝒗𝒚𝒏 for the time instant (n+1) ∆tUpdate vx and vy using 𝒑𝒏ା𝟏 for the time instant (n+1) ∆tVerify the boundary conditionInject source velocity to vx and vyIncrement time step n=n+1n = End of time step?StopNoYesFigure 2.15: An illustration of the FDTD time-marching algorithm.velocity components are updated for the next time instant (n+1)∆t using Equations2.13 and 2.14. This process will be repeated as defined in the time-marching al-gorithm. During the simulation, the boundary condition can be enforced inside thecomputational domain, and source acoustic energy can be induced to visualize thewave propagation. The algorithm implements perfectly matched layers (PMLs),an artificial absorbing boundary condition at the domain boundary, to simulate thefree-field radiation. Figure 2.16 and 2.17 below illustrate the acoustic renderingof a sinusoidal signal in free space using the FDTD numerical scheme with andwithout PML layers, respectively.36Figure 2.16: An illustration of acoustic pressure wave propagation inside arectangular grid with (500×500) domain size and having PML layersat the domain border. A sinusoidal function is used as the acousticsource inside the 2D domain.37Figure 2.17: An illustration of acoustic pressure wave propagation inside arectangular grid with (500×500) domain size and having no PML lay-ers at the domain border. A sinusoidal function is used as the acousticsource inside the 2D domain.382.4 Data Acquisition For Vocal Tract Structural ModelFor this thesis work, we build the area function-based vocal tract structural models.We use Story’s area function [109] dataset to construct the 2D vocal tract contourson a mid-sagittal plane for the following static vowel sounds: /A/, /i/, /u/ and /æ/.The primary reason behind using this dataset is that it is at the root of many articula-tory models, especially the 3D FEM vocal tract model [12]. Since we compare theproposed model’s acoustic characteristics (2.5D FDTD) with 3D FEM, 2D FDTDand 1D vocal tract models, we use Story’s vocal tract area functions for all of them.Story used the MRI imaging technique to capture the American English vowelsounds. The procedures for acquiring image sets have been discussed here [108]in detail. Image processing algorithms like segmentation and shape-based inter-polation were used to separate tissues and construct the 3D airway, respectively.As discussed in Section 2.2.2, the cross-sectional area was extracted at a regularinterval just above the glottis to lips.As we implement the finite-difference time-domain numerical scheme to com-pute the pressure wave propagation, the model uses a rectangular 2D Yee grid tobuild the vocal tract geometry. Appendix B shows the area function and the cor-responding vocal tract structural models inside the numerical domain for differentstatic vowel sounds.2.5 SummaryThis chapter has reviewed the existing state-of-the-art articulatory speech synthe-sis models, emphasizing the physics-based vocal tract modelling methods. It hasalso highlighted the advantages and current limitations of a vocal tract’s structuraland acoustic modelling with great details. To alleviate some of these challenges,we have demonstrated the acoustic analysis using the FDTD in 2D, which is com-putationally less expensive than the FEM-based approach. The following chaptersexplore the constraints in the existing 2D vocal tract modelling methods and pro-pose a new technique that strikes a balance between the computational load andunerring acoustic feature representation in a fixed size domain.39Chapter 3Limitations of 2D FDTD VocalTract ModelsIrrespective of different numerical approaches, the 2D vocal tract models provide abetter trade-off between computational overhead and precise estimation of acous-tic features for the open static vocal tract shapes (e.g., cardinal vowel sounds in theEnglish language) compared to the 3D and 1D models [85, 146], respectively. Thisenhanced effect derives from its ability to simulate transverse modes of acousticwaves in two dimensions only. Hence, in contrast to 1D, the 2D articulatory speechsynthesizer promises better geometrical flexibility for vocal tract modelling. The-oretically, the geometrical versatility and the onset of transverse modes in the 2Dapproach help estimate acoustic formants correctly beyond the planar wave regime.Therefore, the 2D models are always preferred for acoustic simulation of openstatic vocal tract shapes having complex contours. Similarly, the 2D models arecomputationally cheaper than the 3D vocal tract models due to the lower dimen-sion. It should be noted that the 3D models facilitate the design of realistic vocaltract geometry inside a computational domain, which can not be achieved usingthe current 1D and 2D vocal tract models. However, the existing 2D models donot use the original vocal tract cross-sectional area functions extracted from MRIsfor the vocal tract structural model. They are also not suitable for designing thereal-time articulatory speech synthesizer. In this chapter, we explore the modellingtechnique of a 2D vocal tract and discuss its current limitations.40Speed [102] and Wang [137] demonstrated the synthesis of vowel sounds us-ing the lightweight 2D FDTD method, but far from real-time. Nevertheless, thesemethods face challenges when the vocal tract boundary is dynamically modifiedto produce diphthong utterances. Zappi et al. [146] proposed a novel 2D FDTDacoustic wave solver that takes advantage of massive core parallelism available on acommodity graphics card (GPUs) for a faster acoustic simulation of 2D vocal tractswhile maintaining the model stability for the dynamic boundary condition. How-ever, Zappi’s proposed 2D model does not capture the nonplanar acoustic wave in-teractions across the sagittal plane, and a direct representation of 3D tube geometryinside a 2D computational domain does not yield precise acoustic output. Alter-natively, a high spatial grid resolution can be enforced along with an area functionscaling method to tune the model’s acoustic output at the expense of GPU’s per-formance [8]. This area function transformation helps adjust the wave reflectioncoefficients of a 2D vocal tract model to match a 3D model and, consequently,improves its acoustic features.This chapter considers Zappi’s 2D FDTD vocal tract as a reference to analyzethe shortcomings in the current 2D modelling technique and set a framework fromwhich our proposed vocal tract model (2.5D FDTD) stems. It is structured as fol-lows. Section 3.1 presents the development of a tuned 2D vocal tract model. In sec-tion 3.2, we compare the transfer function of a vocal tract shape (vowel sound /A/),generated using the 2D FDTD scheme to a complex 3D FEM vocal tract model.This analysis exposes the underlying limitations of the 2D vocal tract model anddiscusses the feasibility of the present solution to overcome them. Section 3.3 pro-poses a novel approach that improves upon the existing 2D model by addressingits limitations through its modelling technique. Section 3.4 summarizes the signifi-cance of the proposed solution relative to the 2D FDTD scheme in terms of overallvocal tract acoustic modelling.41Figure 3.1: (left) An illustration of 3D vocal tract geometry for vowel [A]with a straight cylindrical tube having circular cross-sections. (right) Amid-sagittal cut of the 3D vocal tract geometry to retrieve the 2D vocaltract contour.llDi Di+1fi(t)fi+1(t)bi(t)bi+1(t)llAi Ai+1fi(t)fi+1(t)bi(t)bi+1(t)𝑟′3𝐷 =𝐴𝑖 − 𝐴𝑖+1𝐴𝑖 + 𝐴𝑖+1𝑟′2𝐷 =𝐷𝑖 − 𝐷𝑖+1𝐷𝑖 + 𝐷𝑖+1Figure 3.2: Due to acoustic wave reflections at the tube junction, the pres-sure pi(t) at any tube segment (i) is the resultant sum of the fi(t) andbi(t), the forward propagating and backward propagating pressure wavecomponents. These phenomena (the effects of wave scattering at thetube junction) are modelled using the reflection coefficient r′, which de-pends (top) on the adjacent tube cross-sectional areas ratio for a 3D tube(bottom) on the adjacent tube diameter/tube width ratio for a 2D tube.42?a3D(x) a2D(x) –Figure 3.3: Determining a scaling function r(x) to retrieve the area functionfor the 2D vocal tract contour, a2D(x) using the original 3D vocal tractcross-sections, a3D(x). x represents the length of the vocal tract fromthe glottal-end. This scaling function ensures the spectral features of a2D model match a 3D vocal tract model.3.1 2D FDTD Modelling Methodology3.1.1 Problem StatementA 3D vocal tract model (juxtaposed cylindrical tubes/ducts) extends the 1D acous-tic wave solver to capture wave interactions across all three planes (axial, coronaland sagittal) inside the vocal tract boundary. In contrast, a 2D model uses the mid-sagittal slice of a 3D geometry (see Figure 3.1) to define the 2D vocal tract contoursinside a computational domain. Due to the dimensionality reduction in a 2D model(absence of geometrical features across the sagittal plane), the acoustic wave solverestimates wave propagations along the sagittal plane only. Moreover, the dissimi-larities in the tube geometry of a 3D and 2D model influence their wave reflectioncoefficients between two consecutive tube segments as shown in Figure 3.2. Dueto their disparities in the acoustic wave reflection coefficients, which help in mod-elling the wave scattering at the tube boundary, the 2D model’s acoustic features,generated using the original cross-sections of a tube are not comparable with a re-alistic 3D model. This problem could be addressed by changing the r′2D to matchthe r′3D, consequently modifying tube segments’ widths of a 2D vocal tract. Hencefor a 2D model, we need to find a new area function configuration by using a scal-ing function (Figure 3.3) that matches the acoustic behaviour of a 3D model. Forthis purpose, the next section demonstrates Arnela’s proposed algorithm to derivethe scaling function and the modified 2D vocal tract area functions [8].40r(x)Figure 3.3: Determining a scaling function r(x) to retrieve the area functionfor the 2D vocal tract contour, a2D(x) using the original 3D vocal tractcros -sections, a3D(x). x represents the length of the vocal tract fromthe glottal-end. This scaling function ensures the spectral features of a2D model match a 3D vocal tract model.3.1 2D FDTD Model ing Methodology3.1.1 Problem StatementA 3D vocal tract model (juxtaposed cylindrical tubes/ducts) extends the 1D acous-tic wave solver to capture wave interactions acros all thre planes (axial, coronaland sagittal) inside the vocal tract boundary. In contrast, a 2D model uses the mid-sagittal slice of (see Figure 3.1) to define the vocal tract contoursinside a co putational i . to the dimensionality reduction ( bsence ofgeom trical featur s across he gittal plane), he 2D acoustic wave solver esti-ma es wave propagations along the sagitt l pl ne o ly. Moreov r, the dissimilar-ities n the tube g ometry of a 3D and 2D model influence their wave reflectioncoefficients betwe n two consecutive tube segments as shown in Figure 3.2. Dueto their disparities in the ti i help modelthe wave scattering at the tube junction, the 2D model’s acoustic features gener-ated using the original cross-sections of an actual vocal tract are not comparablwith the 3D model. However, this problem could be solved by changing r′2D (i.e.,wave reflection coefficient o a 2D tube) to match r′3D (i.e., wave reflec ion coeffi-cient of a 3D tube), consequently modifying t be segme ts’ widths of the originvocal tra t. Hence for the 2D voc l tract, we need to find a new area function con-figuration by using a scaling function, named r(x) (see Figure 3.3) that matchesthe acoustic behaviour of a 3D vocal tract model. For this purp se, the following43section demonstrates Arnela’s proposed algorithm to derive the scaling function,which will then be used to extract the modified 2D vocal tract area functions fromthe actual vocal tract geometry [8].3.1.2 Step 1: Nonlinear Scaling of 3D Vocal Tract GeometryConsidering we already have the 3D cross-sectional details of a cylindrical tube(Section 2.4), Arnela proposed a nonlinear scaling method [8] to obtain the 2Darea functions. This approach ensures the wave reflection coefficients of a 2Dtube match well with its 3D description [25]. As a result, the acoustic features ofsynthesized audio outputs generated using the 2D vocal tract model improve. Thearea functions transformation method is stated as follows:Consider 3D and 2D vocal tracts (straight tube) of length L and circular cross-sections. We can divide each of these vocal tracts into N ducts, i.e., x1,x2, ...,xi,xi+1, ...,xN , where x1 and xN represent the glottal-end and mouth-end cross-sections, re-spectively. These discretized duct segments can be visualized as expansion cham-bers having the same length but different cross-sections. Considering, the ductarea functions for a 3D vocal tract are a13D,a23D, ...,aN3D and for a 2D vocal tract area12D,a22D, ...,aN2D respectively, Boij and Nilsson [25] have proposed that the 2D nu-merical simulation matches well to 3D experimental results whenever the 2D areaexpansion ratio of two consecutive duct segments is equivalent to the 3D one forthe plane wave propagation. Imposing this condition implies,ai+12D =[ai+13Dai3D]ai2D (3.1)or,ri+12D =[ri+13Dri3D]2ri2D (3.2)where, ri3D and ri2D denote cross-sectional radius for a 3D and 2D vocal tract re-spectively. The Equation 3.2 allows to construct the 2D vocal tract geometry oncethe radius rd2D at an arbitrary duct segment (xd), with 0≤ xd ≤ L is determined. Fora known rd2D, the Equation 3.2 can be rearranged as follows,44ri2D =ri+12D mi, 1≤ i < drd2D, i = dri−12D mi, d < i≤ N(3.3)where,mi =(ri3D/ri+13D)2, i < d(ri3D/ri−13D)2, i > dFinally, we need to determine a specific duct segment (xd) and its cross-sectionalradius (rd2D) to implement the nonlinear scaling function (Equation 3.3). This canbe achieved by considering the onset position of the first nonplanar mode in 2D and3D vocal tracts and their corresponding frequency values [94]. In a vocal tract, thefirst nonplanar mode is induced at a duct segment having the largest cross-sectionalarea. Hence, rd2D can be precisely computed by taking into account the frequencyvalue of the first non-planar mode in a 3D and a 2D duct as follows [48]:f 3Dc =1.84pic02rd3D, f 2Dc = 0.5c02rd2D(3.4)where, c0 is the speed of sound. We can enforce f 2Dc = f3Dc to obtain rd2D, whichneeds to be used in Equation 3.3.rd2D = rd3D0.5pi1.84(3.5)Here, xd is the distance of the duct segment having the largest cross-section fromthe glottis (x1) in the 3D tube. Figure 3.4 shows the new area function configura-tions of static vowels (/a/, /i/, /u/ and /æ/) for the 2D vocal tract model.450 0.02 0.04 0.06 0.08 0.1 0.12 0.14 0.16 0.18Distance from the glottis [m]00.10.20.30.40.50.60.70.80.91Cross-sectional area [m2 ]10-33D Area Function2D Area Function(a) vowel /A/0 0.02 0.04 0.06 0.08 0.1 0.12 0.14 0.16 0.18Distance from the glottis [m]0123456Cross-sectional area [m2 ]10-43D Area Function2D Area Function(b) vowel /i/0 0.02 0.04 0.06 0.08 0.1 0.12 0.14 0.16 0.18Distance from the glottis [m]0123456Cross-sectional area [m2 ]10-43D Area Function2D Area Function(c) vowel /u/0 0.02 0.04 0.06 0.08 0.1 0.12 0.14 0.16 0.18Distance from the glottis [m]00.511.522.533.544.5Cross-sectional area [m2 ]10-43D Area Function2D Area Function(d) vowel /æ/Figure 3.4: 2D vocal tract area function (red line) for vowels /a/, /i/, /u/ and/æ/ after nonlinear scaling of the original 3D area function.463.1.3 Step 2: Tuning Grid ResolutionThe nonlinear scaling function (Step 1) adjusts the formants of the 2D vocal tractto match with a 3D model. However, this approach tapers off the contour widthof a 2D tube as it scales down the original area function. Therefore, the 2D vocaltract model loses its geometrical resolution and boundary complexity, which are theessence of its acoustic properties. Alternatively, the spatial grid resolution of thecomputational domain can be tuned by decreasing grid cell size for a precise repre-sentation of vocal tract boundary (shown in Figure 3.5). As the acoustic output inan articulatory speech synthesis model heavily depends upon the correct approx-imation of vocal tract boundary, we can conclude that higher spatial resolution isalways preferred in the 2D modelling environment.20 40 60 80 100 120 140 160 180 200 22010203040Before scaling area function [grid cell size of ∆s = 0.74 mm, domain size = 232×49]20 40 60 80 100 120 140 160 180 200 220510152025303540After scaling area function [grid cell size of ∆s = 0.74 mm, domain size = 232×43]47100 200 300 400 500 60020406080100After scaling area function [grid cell size of ∆s = 0.28 mm, domain size = 612×109]100 200 300 400 500 600 700 800 90020406080100120140160After scaling area function [grid cell size of ∆s = 0.18 mm, domain size = 917×161]Figure 3.5: Illustration of changes in tube geometry for a 2D vocal tract con-tour (vowel [A]) due to area functions scaling and grid resolution tuning.The x and y axis of the above plots represent the number of grid cells.3.2 2D FDTD Vocal Tract Transfer FunctionWe employ the 2D FDTD acoustic wave solver to derive the vocal tract transferfunction for the vowel sound /A/. During the numerical simulation, we set the com-putational domain and system parameters as proposed in Section 5.3. An open-endboundary condition is enforced at the mouth-end, and a band-pass velocity pulseis positioned at the glottal end as a source excitation model. To account for the 2Dvocal tract boundary losses, we follow a method as proposed in Section 4.4. Theextracted formants from the vocal tract transfer function are then compared with a48highly accurate 3D FEM model [12]. Basically, the 2D FDTD scheme employs thefollowing discrete update rules at each grid point to sample the acoustic pressureand velocity components at each time step1,2.p(n+1) =p(n)−ρc2∆t ∇˜ ·V (n)1+(1−β )∆t (3.6)(3.7)v(n+1) =βv(n)−β 2∆t ∇˜p(n+1)/ρ + (1−β )∆tvbβ +(1−β )∆t (3.8)with:V = (vx, vy ) (3.9)To analyze the impact of the nonlinear area function scaling method on the 2Dvocal tract model’s acoustic outputs, we build vocal tract contours of vowel sound/A/ inside the 2D numerical grid in the following manners: (1) by directly using 3Dcross-sectional area functions (2) by using the resultant transformed cross-sectionsobtained from Step 1 (Section 3.1.2). From now on, they will be referred to asarea functions without scaling (WOS) and with scaling (WS), respectively. Thisexperiment primarily examines how the area function scaling method adjusts theformant frequencies in the synthesized speech waveforms generated using the 2DFDTD scheme. Next, we choose three different grid resolutions for the vowelsimulation as different test conditions to study the influence of precise vocal tractboundary representation in a 2D numerical domain: ∆s = 0.74 mm (661.5 kHz),∆s= 0.28 mm (1764 kHz) and ∆s= 0.18 mm (2646 kHz). Henceforth, these valueswill be referred to as 2Dlow, 2Dmid and 2Dhigh, respectively. We have chosen thesegrid resolutions empirically for our experimental study of the 2D FDTD model.However, setting a grid resolution higher than 0.18 mm is not a pragmatic approachdue to the computational simulation run time. Furthermore, the 2D vocal tractmodel loses much of its geometrical complexity at a grid resolution lower than 0.74mm. The simulation results for the above conditions are provided in Table 3.1.1∇˜ the standard discrete spatial derivatives2The definition of β parameter has been provided in Section 4.2.249Formants 3DFEM2Dlow(WOS)2Dlow(WS)2Dmid(WS)2Dhigh(WS)F1 696 Hz660.02 Hz(-5.16%)660.02 Hz(-5.16%)660.00 Hz(-5.17%)700.00 Hz(0.57%)F2 1068 Hz1220.03 Hz(14.23%)1040.03 Hz(-2.61%)1040.01 Hz(-2.62%)1060.00 Hz(-0.74%)F3 3031 Hz2680.08 Hz(-11.57%)2960.08 Hz(-2.33%)3000.03 Hz(-1.02%)3020.02 Hz(-0.36%)F4 4064 Hz3740.11 Hz(-7.97%)3680.11 Hz(-9.44%)4000.04 Hz(-1.57%)3960.02 Hz(-2.55%)F5 5039 Hz4660.14 Hz(-7.52%)4940.14 Hz(-1.96%)4820.05 Hz(-4.34%)4900.03 Hz(-2.75%)F6 5703 Hz5560.16 Hz(-2.50%)5580.16 Hz(-2.15%)5560.06 Hz(-2.50%)5600.04 Hz(-1.80%)F7 7065 Hz6720.20 Hz(-4.88%)6360.19 Hz(-9.97%)6840.07 Hz(-3.18%)6720.05 Hz(-4.88%)F8 7595 Hz7520.22 Hz(-0.98%)7480.22 Hz(-1.51%)7320.08 Hz(-3.61%)7420.05 Hz(-2.30%)Table 3.1: Formant locations Fi (Hz) and formant deviations ∆Fi/Fi (%) withrespect to the 3D FEM model in the existing 2Dlow (without scaling); and2Dlow, 2Dmid and 2Dhigh (with scaling) for vowel /A/.503.2.1 Transfer Function AnalysisThe following section analyzes acoustic formants of the 2D FDTD wave solver forthe vocal tract shape /A/ (results from Table 3.1) to expose its limitations in termsof vocal tract modelling technique.• The percentage of formant deviations between the 2D FDTD and 3D FEMmodels were high with the direct use of 3D vocal tract cross-sections insidethe 2D computational domain (2Dlow-WOS) for acoustic simulation of vowelsound /A/. These deviations can be explained by considering the mismatchof wave reflection coefficients between the 3D and 2D vocal tract models.• The 2D wave solver exhibited a remarkable improvement in its acoustic for-mants’ positions after using the scaled area functions as the actual vocaltract geometry inside the 2D numerical domain. Though this technique didnot capture the nonplanar acoustic waves, it adjusted the acoustic formantsof the 2D model, which are comparable to the 3D vocal tract model. Hence,we conclude the existing 2D models always require a nonlinear area func-tion transformation to compensate for the missing geometrical information.However, the proposed solution adversely influences the boundary complex-ity, especially the vocal tract shapes with narrow constrictions.• It is to be noted that the area function scaling approach for the dynamicboundary condition is much more complex. In this case, during the simu-lation, the 2D vocal tract contours change temporally in a step-wise fashioninside the computational domain. Therefore, we need to determine a newarea function for the 2D vocal tract model at each timestep as the vocal tractboundary walls move. Additionally, the area function scaling method doesnot work for the bent tube structure as the tube segments (ducts) are notjuxtaposed anymore.• Though the proposed solution by Arnela works for a straight tube with cir-cular cross-sections, it demands a higher grid resolution to produce preciseacoustic formants for a vocal tract shape. For example, the formants’ de-viations for the 2Dhigh (WS) simulation were much lower compared to the2Dlow (WS) vocal tract model.51• Contrarily the higher grid resolution requires a longer simulation durationas the computational overload is notably high, which is not favourable fordesigning a real-time articulatory speech synthesizer.3.3 2.5D VT Model - An Extended 2D Vocal Tract ModelTo improve the existing 2D model, we propose a new vocal tract model (2.5D VTmodel), incorporating the below features to the existing 2D acoustic wave solver.A detailed modelling methodology has been provided in the next chapter.• The new model lumps the effect of off-plane acoustic waves by capturingtube geometry along z-axis (i.e., tube height/tube depth). In other words,the 2D wave solver seizes the acoustic wave propagation across the sagittalplane virtually.• This additional geometrical information (i.e., tube depth) helps character-ize the wave reflection coefficient of adjacent tube segments as the same forboth the 2.5D and 3D models. Therefore, unlike the 2D model, the 2.5D vo-cal tract does not suffer from the same constraints, and we can directly usearea functions that closely represent the actual vocal tract geometry. This ap-proach eliminates the requirement of scaling the original 3D cross-sectionalareas for a 2D vocal tract model.• The acoustic wave solver of the proposed model also adopts a mechanism toaccount for the boundary losses due to the tube surface that extends acrossthe sagittal plane. That allows a detailed study of vocal tract tubes with dif-ferent cross-sectional shapes and their impact beyond the plane wave propa-gation regime.523.4 SummaryThis chapter thoroughly discussed the 2D FDTD vocal tract modelling approachfor an articulatory speech synthesizer. We started by exploring the primary limita-tion in a 2D structural model since it does not capture the acoustic interaction ofoff-plane waves. Hence, a direct mid-sagittal representation of a tube built fromthe 3D vocal tract’s area function produces erroneous formant locations, i.e., notmatching the results obtained from the sophisticated 3D acoustic simulation of thesame vocal tract shape. Thus, the extracted 2D contours have to undergo a non-linear transformation to fix these discrepancies. However, this approach leads to asignificant downsizing of the vocal tract’s constrictions. The 2D simulations needto run with high spatiotemporal resolutions to model such narrow constrictions,thus increasing the computational load of the solver and producing waiting timesfar from real-time performance. Moreover, the methodology proposed by Arnelaand Guasch applies only to the case of straight tubes, and it is still unclear how toobtain 2D contours that acoustically match curved and branching 3D geometries.To improve upon the existing 2D vocal tract model, we have proposed a newvocal tract modelling technique (2.5D FDTD), which can avoid the current imped-iments of the 2D FDTD wave solver. The modelling methodology and wave solverequations for the new model are discussed in detail in the next chapter.53Chapter 42.5D FDTD Vocal Tract AcousticModel For Static VowelsIn this chapter, we present a novel approach to vocal tract modelling that stemsfrom the 2D rationale [146] and improves upon it. It capitalizes on an extended 2DFinite-Difference Time-Domain wave solver (2.5D FDTD), capable of simulatinghow pressure wave propagates in 3D geometries that are symmetric at least acrossa plane, typically across the xy-plane (with x and y being the dimensions of thestarting 2D scheme). At the core of this approach, there is the inclusion within themodel of an additional impedance term which is derived from the tube’s depth (i.e.,its continuous extension along z axis). The tube depth is sampled in every pointof the scheme; the result is a depth map, that combined with the 2D mid-sagittalcontour of the original 3D geometry, allows for fast simulation of its acoustics,with computational requirements comparable to the case of standard 2D numericalschemes.Depth maps can be retrieved from area functions as well as from full 3D mod-els of actual vocal tracts, like Magnetic Resonance Imaging scans. Furthermore,the overall 2.5D representation of the analyzed geometry leaves its original pro-portions intact, discarding the need for nonlinear transformation of the vocal tractarea functions. Consequently, this approach removes the requirement for extremelyhigh spatiotemporal resolutions that characterize 2D simulations and amelioratessimulation run-time.54Geometry Extraction2.5D Model1D Model2D Model3D ModelVocal Tract Structural ModelSS-Depth Contour & depthSS-Curvature Tube CurvatureSS-Branches Side BranchesDS-Diphthongs Dynamic VTSS-Vowel Static VT Shape1D Model2D Model3D ModelVocal Tract Acoustic ModelSynthesized Speech2.5D ModelExisting ModelsExisting ModelsStationary BoundaryDynamic BoundarySA-Loss VT Boundary LossesSA-Vowel Vowel SoundsSA-Radiation Mouth RadiationSA-Analytical Analytical TubeStationary BoundaryDynamic BoundaryDA-Loss VT Boundary LossesDA-Diphthongs DiphthongsDS-Depth Contour & depthDS-Curvature Tube CurvatureDS-Branches Side Branches DA-Radiation Mouth RadiationFigure 4.1: Highlighted steps (yellow) in the modified roadmap which areaddressed in this chapter.554.1 2.5D Vocal Tract ModelAs discussed in Chapter 3, the 2D simulation of a vocal tract tube does not capturethe interaction of nonplanar waves across the mid-sagittal plane. The modellingtechnique itself presumes acoustic waves propagate along the xy-plane only (mid-sagittal plane) under a restrictive boundary condition (2D tube contours), but theyare unbounded along the z axis. However, the 2.5D simulation of a vocal tract tubeallows incorporating wave propagation effects for off-plane waves by adding tubedepth (i.e., the extension of tube height along the z axis) to the generic 2D FDTDacoustic wave solver equations. This approach has already been demonstrated forcomputational modelling of brass instrument [18], but using 1D acoustic wavesolver only. Like a 3D vocal tract, the 2.5D simulation also considers vocal tractwall losses due to the nonplanar wave interactions. Hence, theoretically, it is ex-pected that a 2.5D vocal tract tube’s acoustic characteristics should be comparableto its corresponding 3D representation. Moreover, contrary to a 2D model, the2.5D vocal tract does not require nonlinear cross-sectional area function transfor-mation as it captures the complete geometrical details of the 3D vocal tract tube,which are symmetric across the xy-plane.Since the development of the 2.5D FDTD acoustic wave solver is still under-way, in this chapter, we will go through some of the crucial steps of our roadmap(see Figure 4.1), which guide the underlying principles of the 2.5D vocal tractmodel. First, SS-Depth: Deriving tube depths for a vocal tract model from its MRI-based cross-sectional area functions and mapping them correctly to the 2D acousticwave solver. This technique is a major step toward improving the traditional 2Dwave solvers [8, 146] as it allows capturing the effect of nonplanar acoustic wavesvirtually. Second, SA-Loss: Modelling vocal tract boundary losses to match theacoustic characteristics of an actual vocal tract. However, the current implemen-tation of the 2.5D model only accounts for wall losses due to the 2D mid-sagittalcontours and ignores losses for tube surfaces across the xy-plane. Consequently,we herein have partially addressed the requirements of SA-Loss. Though bound-ary losses influence the tube’s acoustic properties, it has been left for the futuredevelopment of the 2.5D model. Finally, we demonstrate the structural modellingof open static vocal tract shapes using the 2.5D FDTD scheme (SS-Vowel).56x-axisy-axisz-axisa1b1Surface ASurface BFigure 4.2: Graphical representation of wave propagation in xy-planebounded by two rigid surfaces (A and B). The euclidean distance be-tween the two surfaces is |a1b1|= D(x,y).4.2 2.5D Vocal Tract AcousticThis section illustrates the implementation of the first step in our roadmap: SS-Depth. We first introduce the acoustic wave equation for the proposed 2.5D modeland then outline a detailed procedure to combine it with the existing 2D FDTDscheme.4.2.1 Acoustic EquationLet us consider a sound wave propagating in air, and let us call ξ the displacementvector of the medium caused by the motion of the wave and V the volume of ageneric element of the medium (particle). Now suppose that the medium is con-tained between two rigid surfaces that extend towards infinity along both the x andthe y axis without ever intersecting and that are one the reflected image of the other(being x/y the plane of symmetry). In the described scenario, the wave moves ina space that is 3D but constrained on one axis (i.e., z); such a space can be de-scribed by the scalar field D(x,y), which contains the euclidean distance betweenthe points of the two surfaces at the same x,y coordinates (shown in Figure 4.2). Inother words, D describes the space in terms of its depth along z.57If we choose to focus on a particle of volume V = dxdyD, with dx and dyrelatively small, ξ will only have x and y components due to the presence of theenclosing surfaces. Consequently, during the passage of the sound wave, the parti-cle’s volume change can be expressed as:dV = (∂Dξx∂xdx) dy+(∂Dξy∂ydy) dx (4.1)if ξ is quite small (reasonable assumption in a realistic scenario) and D isslowly varying in space. The two products in parentheses represent the first orderTaylor Series approximations of the increment of ξ and D, on x and y respectively.Hence, by choosing dx = dy = ds, the fractional volume change can be written as:dVV=1D(∂Dξx∂x+∂Dξy∂y) (4.2)If we describe the motion of the particle via Newton’s equation and we combineit with the equation of state p = −K dVV (where K is the bulk modulus of air), itis possible to prove the following relationship between the displacement and thefractional volume change as follows:∂ 2ξ∂ t2=c2D∇(∂Dξx∂x+∂Dξy∂y) (4.3)with c defined as the speed of sound in air. Equation 4.3 can be derived by extend-ing the acoustic wave equation (see Appendix A.1) to 2D and replacing the frac-tional volume change as given in Equation 4.2. Finally, by following Fletcher andRossing [94] and considering D independent of time, we can obtain from Equation4.3 the following acoustic wave equation:∂ 2 p∂ t2=c2D(∂∂x(D∂ p∂x) +∂∂y(D∂ p∂y) ) (4.4)where p is sound pressure. The equation describes airwave motion in a 3D space,58constrained in 1D; given the similarity with the standard 2D acoustic wave equa-tion, this 3D space can be considered as an extended 2D (2.5D) space. By sym-metrically constraining the 2.5D space along the y dimension too (hence turning Dinto a circular area), and by assuming plane-wave propagation only (due to consid-erable distance from the source), Equation 4.4 reduces to Webster’s horn equation[125].4.2.2 2.5D FDTD Wave SolverThe procedure followed for the design of the 2.5D wave solver shares many sim-ilarities with the case of the 2D vocal tract model described in [146]. We decom-posed the 2.5D acoustic wave equation into extended versions of the 2D equationof continuity and 2D equation of motion:∂ p∂ t+ σ ′p = −ρc2D(∂Dvx∂x+∂Dvy∂y) (4.5)β∂v∂ t+ σ ′v = −β 2∇pρ+ σ ′vb (4.6)where v is the 2D acoustic particle velocity. And, we define the effective absorp-tivity σ ′ = 1−β +σ , where σ is the PML absorptivity used at domain boundariesas described in Section 5.3.3. However, similarly to what proposed by Allen andRaghuvanshi [3], we augmented the obtained standard 2D equation with the scalarparameter field 0≤ β (x, t)≤ 1. By varying between 1 (air) and 0 (boundary), thisterm allows for the transition between the momentum equation and the enforce-ment of a prescribed velocity vb . Through this mechanism, it is possible to enforceboundary conditions (more details in the Section 4.4) and to simulate dynamicgeometries to generate utterances, like diphthongs [146]. Equation 4.6 does notinclude any explicit dependence on the depth of the considered space (via D). Thisis due to the fact that the equation of motion derives from Newton’s second law,whose 2D and 2.5D forms coincide.We applied a standard 2D Yee scheme, where each grid point consists of asquared 2D cell [3, 141]. Per each time step n, pressure values p(n) are sampledacross the whole domain at the center of the cells, while the v(n)x and v(n)y compo-59nents of the velocity vectors are sampled on the edges shared with the right and topneighbor cell respectively. Hence, in a generic cell at discrete coordinates (xˆ, yˆ), pwill be sampled at (xˆ, yˆ), while vx at (xˆ+1/2, yˆ) and vy at (xˆ, yˆ+1/2). This leadsto the following discrete update rules of Equations 4.5 and 4.6, where we denotewith ∇˜ the standard discrete spatial derivatives as performed in 2D FDTD:p(n+1) =D¯p(n)−ρc2∆t ∇˜ ·V (n)D¯(1+σ ′∆t)(4.7)v(n+1) =βv(n)−β 2∆t ∇˜p(n+1)/ρ + σ ′∆tvbβ +σ ′∆t(4.8)with:V = ( D(x)vx, D(y)vy ) (4.9)The terms D¯, D(x) and D(y) are the components of the depth map of the 2.5D space,and they effectively act as additional impedance terms in two dimensions. Theirrelation with the field D will be explained in detail in the following subsection.4.3 Depth MapThe depth map is a discretized version of the continuous field D, processed ac-cording to the described 2D scheme. Since every individual FDTD cell holds threeacoustic parameters sampled in three different locations (center, right and top), thediscretization of D will likewise produce three depth values per each cell, namelyD¯, D(x) and D(y). As a result, in the FDTD grid, the pressure values are alignedwith D¯ values, while the velocity vectors are aligned with D(x) and D(y) values.In a practical scenario, the field D can be replaced by the 3D model of the ge-ometry whose acoustics is being analyzed. Even if typically composed of discreteelements (e.g., vertices), the model can be treated as a continuous geometry viainterpolation. However, D always represents a symmetric domain; this is not nec-essarily true for a 3D model, hence additional processing is required to derive thedepth map from it. In the case of a generic 3D tube, the depth map is computed asfollows:60Figure 4.3: on the left, depth map extraction (D(x) solid line, D(y) dashed line)from a slice C of an irregular tube, at coordinate x = xˆ (for simplicity weassume the same shape at xˆ+1/2). On the right, equivalent mirror-planesymmetric cross-section C′, as represented in the 2.5D FDTD.1. a 2D contour is extracted as the intersection between the tube and its mid-sagittal plane (i.e., 2D vocal tract geometry). The contour itself representsthe vocal tract’s walls, while the mid-sagittal plane becomes the x/y plane ofthe 2.5D space;2. in every wall cell, D¯, D(x) and D(y) are set to zero;3. in every cell outside of the the tube’s walls, D¯, D(x) and D(y) are assigneda fixed open space depth of 0.50m. We select a larger open space depthvalue compare to the vocal tract cross-sectional width to emulate the effectof acoustic wave propagation along the z axis in an open space;4. for every cell with discrete coordinates (xˆ, yˆ) inside the tube’s walls, thetwo lines perpendicular to the x/y plane and passing by (xˆ + 1/2, yˆ) and(xˆ, yˆ+ 1/2) are intersected with the model (Figure 4.3). The length of thetwo resulting segments is assigned to D(x)(xˆ, yˆ) and D(y)(xˆ, yˆ) respectively;5. for every cell with discrete coordinates (xˆ, yˆ) inside the tube’s walls, D(x)(xˆ, yˆ)and D(y)(xˆ, yˆ) are interpolated as ( D(x)(xˆ, yˆ)+Dx(xˆ+1, yˆ) )/2 and ( D(y)(xˆ, yˆ)+D(y)(xˆ, yˆ+1) )/2 respectively;61(a) (b)Figure 4.4: Vocal tract boundary losses due to the (a) 2D contour on the mid-sagittal plane (b) tube surface across the mid-sagittal plane.6. for every cell with discrete coordinates (xˆ, yˆ) inside the tube’s walls, thedepth D¯(xˆ, yˆ) is obtained as ( D(x)(xˆ, yˆ) + D(x)(xˆ− 1, yˆ) + D(y)(xˆ, yˆ) +D(y)(xˆ, yˆ−1) )/4;7. a minimum depth threshold of 0.001m is defined and applied in every cell(depth values are clamped) to maintain the stability of the system. A typicalthreshold is at least an order of magnitude lower than the smallest non-zerodepth value found across the domain.The result of this process is a single-axis symmetric equivalent of the original3D model. The cross-sections of such a 2.5D shape have the same areas and widths(i.e., extension along y) as the 3D model’s sections; moreover, their shapes can beirregular but always symmetric along z, as shown in Figure 4.3. The interpolationsat steps 5 and 6 are introduced to satisfy the assumption that D is slowly varyingin space, required to obtain Equation 4.1. This enhances the overall stability ofthe solver too. We can use this lossless 2.5D FDTD wave solver to characterize itsacoustic features for an analytical tube model. This study will allow us to explorethe fundamental capabilities of the 2.5D model for a tube with simplistic geometry.A detailed discussion of this has been provided in the next chapter.4.4 Vocal Tract Wall LossesThis section looks over the partial implementation of SA-Loss to account for thevocal tract boundary losses. Unlike 2D, the 2.5D numerical scheme should takecare of boundary losses due to (1) the mid-sagittal 2D vocal tract contours (Fig-62ure 4.4a) and (2) tube surface across the mid-sagittal plane (Figure 4.4b). As de-scribed in Section 4.1, at this stage, the 2.5D FDTD wave solver accounts for theacoustic energy loss due to the wave interactions along the mid-sagittal plane only.In a bounded domain (cavity or fluid-filled duct/tube), like a human vocal tract,the fully reflective boundary condition does not achieve realistic acoustic results.The standard acoustic wave equations assume laminar flow inside the tube model,and they completely ignore the thermal and vortical diffusion at the proximity ofboundary walls [26]. Previously, for many vocal tract models, the complex acous-tic impedance boundary condition has been followed to address the visco-thermalwall losses at the mouth opening and the boundary of the acoustical domain insidethe vocal tract [10, 47]. Since the time-domain acoustic solver governs the novelvocal tract model that we have introduced here, it cannot implement the frequency-dependent losses on the wall. Hence, we can not employ the complex impedanceboundary method using the 2.5D FDTD scheme. Therefore, we adapted the localreactive boundary approach, as proposed by Takemoto and Mokhtari [116], whichpartially considers the wall losses and simulates the wall reflection.The boundary conditions at the vocal tract wall are enforced using vw , thenormal particle velocity going to or coming from a wall:vw =pwZnnˆ (4.10)In equation (4.10), pw represents the acoustic pressure measured on the air cellin front of the wall, while nˆ is the unit vector normal to the wall. These normalparticle velocities on the air cells near the wall are combined with equation (4.8) bysetting vb = vw for every velocity vector vx and vy. The normal acoustic impedance,Zn in equation (4.10) can be computed by using the normal sound absorption coef-ficient αn as follows,Zn = ρc1+√1−αn1−√1−αn(4.11)αn =10.5+0.25(µ+µ−1)(4.12)In equation (4.11), ρ and c are the air density and speed of sound respectively.63(a) (b)Figure 4.5: At a low grid resolution (∆s = 0.74mm) the vocal tract geometryfor vowel /a/ (a) before scaling (b) after scaling.The dimensionless boundary absorption coefficient, αn is related to the boundaryadmittance coefficient µ ∈ (0,1) as represented in equation (4.12).4.5 Model ApplicabilityThe use of the depth map imposes to slightly change the cross-sectional shapes ofthe original geometry by means of mirror-plane symmetry. As a result, the acous-tics of the 3D and of the 2.5D vocal tracts will be somewhat different in the portionof the spectrum that is beyond 5 kHz [10]. However, the 2.5D representation is ca-pable of preserving curvature on the x/y plane, area progression as well as severaldetails of the varying irregular cross-sectional shapes of actual vocal tract geome-tries. These are features beyond the capabilities of 1D and 2D models, and whosecombined acoustics effects are yet to be studied.At the current stage of development, using the 2.5D model to analyze theacoustics of realistic vocal tracts is impractical. The model capitalizes on the so-lution of a lossless system (Equations 4.5, 4.6). The solver adds wall impedanceby means of the 2D walls boundary condition (Section 4.4), but the losses happen-ing on the surface described by the depth map (i.e., 2.5D boundaries) are largelyignored, as the extra depth-related impedance terms D¯, D(x) and D(y) model onlythe effects of the surface’s spatial derivative1. In other words, there is no consistent1This scenario is analogous to the case of Webster’s horn equation, where the absolute size of themodelled tube does not affect propagation.64Figure 4.6: 2.5D tube discretization (mid-sagittal contour and depth map) forvowel /a/.way of matching the nonplanar modes of a generic 3D shape.This issue can be easily overcome for the modelling of 3D tubes with circularcross-sections, which can also include bending. By extending the calculation ofmodes in plates proposed by Fletcher and Rossing [48], it is possible to prove thatnonplanar modes in each section of a lossless 2.5D tube are identical to the onesfound in its 2D mid-sagittal cross-section. As a result, we can apply the method-ology proposed by Arnela and Guasch [8] to match the first nonplanar mode of acircular cross-section by scaling the radii of every section of the tube by the con-stant factor 0.5pi1.84 and by modifying the depth map accordingly (see Equation 3.5).However, this scaling approach does not impact the overall vocal tract boundarycomplexity inside a computational domain. As shown in Figure 4.5, the vocal tractgeometry approximately remains the same after the scaling even at a low resolu-tion.4.6 2.5D Vocal Tract Modelling MethodologyWe follow the below procedures to design and discretize the 2.5D vocal tract tubemodel for open static vocal tract shapes (SS-Vowel) using the 2.5D FDTD numer-ical scheme.• Derive the 2D vocal tract contours from its original 3D cross-sectional areafunctions, obtained through MRI imaging of the actual vocal tract (see Sec-65vxpഥ𝐷D(x)D(y) vyvxpvyDepth MappingFigure 4.7: Representing acoustic parameters (p, vx & vy) and depth compo-nents (D¯, D(x) & D(y)) for each grid cell.tion 2.4) and scale the contour widths by a constant factor of 0.5pi1.84 .• Define the 2D vocal tract contours using the 2D Yee-grid scheme. Thenstagger the scalar pressure field and velocity components for each grid pixel.During the simulation, the number of grid pixels determines the computa-tional domain size.• The grid resolution (spatial resolution) of the acoustic simulation can be de-rived using the CFL stability condition (Equation 2.11). For 2D acousticwave modelling, the maximum Courant–Friedrichs–Lewy (CFL) number ofstable modelling is about 0.707 [74, 75].• Derive the tube depth at each grid point inside the computational domain (asshown in Figure 4.6).• Compute and map the tube depth components (D¯, D(x) and D(y)) for theircorresponding acoustic parameters (p, vx and vy) as shown in Figure 4.7.• Apply boundary impedance condition at the mid-sagittal 2D vocal tract walls(Section 4.4).• Place a source excitation function at the glottal end and a virtual microphonenear the mouth-end to record the outward pressure wave propagation insidethe vocal tract tube. And then switch on the 2.5D FDTD acoustic wave solverto activate the model simulation as shown in Figure 4.8.6650 100 150 2001020304050t=10ms50 100 150 2001020304050t=30ms50 100 150 2001020304050t=50msFigure 4.8: The colormap shows the mid-sagittal view of acoustic pressurewave propagation inside the 2.5D vocal tract tube at three different timeinstants for vowel /a/. A sinusoidal excitation function is used at theglottal end. The x and y axis of the above plots represent the number ofgrid cells.674.7 SummaryThis chapter addresses the shortcomings of existing 2D vocal tract models by in-troducing a novel approach that extends the notion of 2D FDTD wave solver whileadding a new impedance parameter called tube depth to the original 2D acousticwave equations. Section 4.1 conceptualized fundamental characteristics of a 2.5Dvocal tract tube and its potential to estimate precise acoustic wave propagation ef-fects without compromising the simulation performance. Section 4.2 introducedthe continuous and discretized version of acoustic wave equations derived fromWebster’s Horn equation. These equations (equation of continuity and equation ofmotion) are used in designing the 2.5D FDTD wave solver for acoustic simulationof tubes. Moreover, the discretized version encompasses the tube depth parameters(i.e., tube heights at each grid pixel) in the original 2D scheme. Section 4.3 dis-cussed procedures to compute the tube depth components and map them to theircorresponding acoustic parameters in each grid cell inside the simulation domain.However, as discussed in Section 4.5, the development of the 2.5D model isstill in progress, and the current lossless acoustic wave equations do not accountfor boundary losses due to off-plane wave interactions with the tube curvature.Moreover, we have not dealt with the non-symmetric tube depths in this study.These limitations force only to test the basic features of the proposed model. As astepping stone to fully realize the potential of the 2.5D FDTD vocal tract model,acoustic features of the 2.5D FDTD wave solver needs to be validated against theexisting 2D and the more realistic 3D vocal tract models. This comparative studywill provide a clear overview of the impacts of an extra geometrical feature (tubedepth) in wave solver equations on simulated acoustic output and its performance.We have discussed the acoustic characterization of the 2.5D vocal tract model indetail in the next chapter.68Chapter 5An Articulatory Model Using The2.5D FDTD Vocal Tract ModelThe primary goal of our new vocal tract model is to design an articulatory speechsynthesizer that can simulate speech utterances. The human vocal tract geometry isintricate; consequently, it influences the propagating acoustic flow characteristicsand speech outputs. Therefore, a speech synthesizer needs to be examined forvarious speech sounds - e.g., vowel sounds, diphthongs, fricative consonants, stopconsonants, voiceless sounds, etc. This chapter illustrates an articulatory model tosynthesize cardinal English vowel sounds for open static vocal tract shapes. In thisthesis, we only study cardinal vowels as they have the uncomplicated open vocaltract geometries that are conducive to making progress on our roadmap first, andthey are the crucial steps for further development of the 2.5D vocal tract model.5.1 Articulatory Model FormulationAn articulatory speech synthesizer consists of multiple components (biomechani-cal models) with different functionalities (see Figure 5.1). For example, the lungsproduce a continuous airflow and channel them into the larynx through the sub-glottal trachea model. The vocal fold model in the larynx oscillates to convert thecontinuous airflow into a series of periodic air puffs and generates acoustic energy.Next, the vocal tract model filters out the acoustic energy to produce synthesized69Trachea Tube Area FunctionTrachea PressurePropagationGeometryTrachea ModelVocal Fold Structural ModelVocal Fold Flow ModelGeometryGeometryPressurePPVocal TractAreaFunction2.5D FDTDWave SolverGeometryUPVocal Fold Model Vocal Tract ModelSynthesizedSpeechLungPressureFigure 5.1: General outline of an articulatory speech synthesizer.audio outputs. Though each of these components can be modelled independently,their acoustic characteristics depend upon each other. Hence a proper couplingmechanism needs to be followed to assemble these biomechanical models. It hasalready been demonstrated that the subglottal trachea model and lungs have a mi-nor influence on the acoustic output of an articulatory model [46, 135]. Similarly,we do not require a nasal tract model as we are only examining cardinal vowelsin this study. Therefore, we only illustrate an articulatory speech synthesizer bycoupling a vocal fold biomechanical model to the 2.5D vocal tract model.To build an articulatory model, we begin by examining the 2.5D FDTD scheme’scapability to precisely estimate acoustic wave propagation inside a vocal tract liketube shape. As a preliminary step toward characterizing the 2.5D wave solver’sacoustic features, we simulate vocal tract shapes for English vowel sounds andconduct a comparative study against the existing vocal tract models. We also in-vestigate the 2.5D vocal tract radiation model. Figure 5.2 highlights the steps fromthe roadmap, which are addressed in this chapter. Section 5.2 discusses the proce-dures that we have followed to characterize the 2.5D vocal tract model’s spectralfeatures. Section 5.3 specifies all the experimental requirements and system spec-ifications that we have adopted throughout our experiment. Section 5.4 largelyproduces a comparative study of acoustic spectra as formants’ positions generatedfrom the 2.5D FDTD numerical simulations and other vocal tract models (e.g., 1D70Geometry Extraction2.5D Model1D Model2D Model3D ModelVocal Tract Structural ModelSS-Depth Contour & depthSS-Curvature Tube CurvatureSS-Branches Side BranchesDS-Diphthongs Dynamic VTSS-Vowel Static VT Shape1D Model2D Model3D ModelVocal Tract Acoustic Model2.5D ModelVocal Fold ModelFT-Coupling CouplingExisting ModelsExisting ModelsStationary BoundaryDynamic BoundarySA-Loss VT Boundary LossesSA-Vowel Vowel SoundsSA-Radiation Mouth RadiationSA-Analytical Analytical TubeStationary BoundaryDynamic BoundaryDA-Loss VT Boundary LossesDA-Diphthongs DiphthongsDS-Depth Contour & depthDS-Curvature Tube CurvatureDS-Branches Side Branches DA-Radiation Mouth RadiationRecordedAudio RS-Spectrum Spectrum AnalysisSynthesized AudioFigure 5.2: Highlighted steps (yellow) in the modified roadmap which areaddressed in this chapter.[131], 2D [146] and 3D [12]) for different vocal tract geometries. It also estimatesthe time complexity of the proposed model. Section 5.5 discusses the significanceof the results in the context of the overall field of vocal tract modelling. Section 5.6provides a brief discussion of the two-mass vocal fold biomechanical model as thesource excitation function for our articulatory speech synthesizer. The vocal foldmodel has been adopted from the literature. In the end, we illustrate a completeend-to-end articulatory synthesizer by coupling the two-mass vocal fold model tothe 2.5D FDTD vocal tract model.715.2 Model CharacterizationThis section offers a comprehensive guideline to characterize the acoustic fea-tures of our vocal tract model at an introductory level. For that, Table 5.1 high-lights SA-Analytical to SA-Radiation and SS-Curvature to DA-Diphthongs fromour roadmap. Since the 2.5D vocal tract acoustic model is still in progress, thetable below only presents a subset of contexts to examine the 2.5D model as perour roadmap. However, Table 5.1 does not provide an extensive list of contexts forstudying all the steps in our roadmap using the 2.5D model. Therefore, the resultsin this chapter can not be generalized for all kinds of vocal tract geometries andspeech utterances.Step ContextsSA-AnalyticalComparative study for an analytical tube model.with open-end boundary condition.(e.g., 1D one-end closed lossless cylindrical tube)SA-Vowel Comparative study for open static VT geometry withopen-end boundary condition. (e.g., cardinal vowels)SA-RadiationCharacterize acoustic features of the open VTradiation model with free-field boundary conditionand head geometry.(e.g., English cardinal vowels)FT-CouplingCharacterize the structural and flow interactions betweenthe vocal fold and vocal tract models.(e.g., vowel sounds and fricative consonants)SS-Curvature Acoustic Analysis for various open VT geometries.(e.g., straight and bent tube with different curvature details)SS-BranchesStudying effects of vocal tract side brancheson the model’s acoustic properties.(e.g., nasal tract and piriform fossa)RS-PerceptionPerceptual evaluation of synthesized speech output.(e.g., synthesized sound from 2.5D vs 3D vs 1D;recorded speech vs synthesized speech)DA-Diphthongs Comparative study for dynamic VT geometry.(e.g., English diphthongs)Table 5.1: Tabular representation of a stepwise procedure to analyze theacoustic features of the 2.5D FDTD vocal tract model. This chapter pro-duces results only for the highlighted section.72ɑai uəBackFrontCloseOpenDiphthong Example/aʊ/ Town/aɪ/ Light/eɪ/ Play/eə/ Pair/ɪə/ Deer/oʊ/ Slow/ɔɪ/ Toy/ʊə/ Sure(a) (b)æFigure 5.3: (a) Vowel quadrilateral representation and (b) Tabular view of En-glish diphthongs as per the IPA (International Phonetic Alphabet) chart.5.2.1 Characterization ProcedureHere we outline a detailed procedure that needs to be followed for each step inTable 5.1 to ensure a complete characterization of the 2.5D vocal tract’s acousticfeatures. The suggested methods are generic; hence they can be applied to othervocal tract models.• SA-Analytical: In a cylindrical tube with uniform cross-sections, the travers-ing longitudinal pressure waves reflect from either its closed or open ends toset up standing wave patterns, consequently producing pressure nodes (i.e.,formants). Due to their simplified geometry, we can easily map their acous-tic formants using analytical approaches. They have been widely used tostudy wind instruments and come with various types such as both-end open,both-end closed and one-end closed tubes. However, the human vocal tractcomes relatively close to being described as one-end closed cylindrical tubes(open-closed tubes). Hence, as a preliminary study, we first analyze for-mants’ locations of a computational one-end closed cylindrical tube modelusing the 2.5D FDTD scheme. It is to be noted that the analytical solutionof a 1D one-end closed cylindrical tube does not take into account the tubeboundary losses. Hence, we build a lossless 2.5D cylindrical tube model in-side the numerical domain to compare the acoustic formants. A lossless tubemodel should not significantly influence the acoustic formants’ positions inthe planar wave propagation regime.73• SA-Vowel: In contrast to a regular cylindrical tube, the acoustic wave in-teractions are quite intricate inside a human vocal tract due to its complexboundary geometry. Hence, we need to study the 2.5D FDTD scheme’sacoustic features using the actual vocal tract shapes. In this work, we onlysynthesize speech sounds that are associated with open static vocal tract ge-ometries. The complicated nature of the acoustic wave interactions in theclosed vocal tract geometries is primarily ignored while modelling the cur-rent 2.5D acoustic wave solver.As open static vocal tract shapes essentially produce vowel sounds,we simulate English vowels, for which vocal tract cross-sectional area func-tion datasets are already available [109]. We select cardinal vowels at thefour corners of the vowel quadrilateral space: /a/, /A/, /u/ and /i/ and at thecentre (/@/) as they represent the extreme points of articulation (see Fig-ure 5.3(a)). However, we do not have the cross-sectional area datasets for/a/ and /@/. Hence, for numerical simulations, we only choose the followingvowel sounds: /A/, /u/, /i/ and /æ/ and present their acoustic outputs in thisthesis. We select the vowel sound /æ/ because it is close to /a/ in the vowelquadrilateral space. We also compare the model’s acoustic spectra with theexisting 1D, 2D and 3D models for complete characterization.• SA-Radiation: The vocal tract radiation mechanism is an essential factor inthe speech production process. It influences the output speech waveform bychanging the formant frequencies and bandwidth. Therefore, we investigatethe 2.5D vocal tract radiation model with a similar approach proposed byArnela et al. [9]. It records the change in the model’s acoustic featuresas the pressure waves emanating from the mouth aperture with a simplifiedhead geometry (e.g., semi-realistic head, circular baffle with circular mouthaperture, and circular baffle with elliptical mouth aperture). However, weonly simulate the radiation model with a circular baffle for the open staticvocal tract shapes in this work.• RT-Coupling: The nonlinear coupling mechanism is an essential factor indetermining the flow-structure interactions between the vocal fold and vocaltract for various speech sounds. It also determines the naturalness of the74Tube Structure Straight, BentTube Curvature Circular, Elliptical,Miror-axial plane symmetryTable 5.2: Various vocal tract geometrical shape types for the 2.5D model.synthesized audio. Thus, we plan to study the acoustic characteristics of the2.5D vocal tract coupled with a vocal fold model. Though we illustrate acoupling mechanism in this study, we do not study its impact due to variousspeech utterances.• SS-Curvature: The human vocal tract’s geometry is quite sophisticated. Dueto its bent and irregular cross-sections, we generally observe antiresonancesand the onset of higher-order modes in the recorded speech waveforms.Hence, every vocal tract model must be verified for its geometrical flexi-bilities and their influence on the model’s acoustic properties. For a specificspeech sound, the 2.5D vocal tract can be modelled in several ways con-sidering its geometrical structure and cross-section, as shown in Table 5.2.Therefore, we conduct a comparative study starting with simulating speechutterances for all the possible vocal tract shapes and then matching their out-puts with a more realistic 3D vocal tract model. However, we have onlyadopted straight tubes with circular cross-sections for all the numerical anal-yses in this thesis work.• SS-Branches: The voice quality and acoustic individuality of a person aregreatly influenced by physiological details such as the geometry of the hy-popharynx, which includes the piriform fossae [59, 69]. Moreover, depend-ing on the type of sound being produced, it is likely that the other articula-tory features (e.g., sublingual cavity, spaces between the teeth, nasal sinuses)may also influence the acoustic output. Therefore, for a realistic speech syn-thesizer, the model should implement as many of these side branches. Theacoustic effects of the nasal tract and piriform fossae have been widely stud-ied using existing vocal tract models [35, 39, 131]. Hence, we plan to incor-porate them into our modelling environment using the 2.5D FDTD schemeand characterize their spectral properties.75• RS-Perception: We propose a psycho-acoustic experiment [133] for the per-ceptual evaluation of the synthesized acoustic outputs. For this, we askhuman subjects to distinguish between synthesized speech outputs gener-ated from different numerical approaches (1D, 2D, 2.5D and 3D vocal tractmodels). Furthermore, we also estimate the similarity between the recordedspeech and synthesized outputs through a user study.• DA-Diphthongs: For an acoustic analysis of the dynamic vocal tract geome-try, we choose 8 English diphthongs (a combination of two adjacent vowelsin the quadrilateral vowel space) as shown in Figure 5.3(b). This strategyoffers an extensive study of the acoustic behaviour of our proposed modeldue to the articulatory movements. The model’s acoustic outputs for dy-namic boundary conditions also need to be compared with the other vocaltract models.To examine the above contexts (except FT-Coupling) for our roadmap, we fol-low the transfer function analysis approach (RS-Spectrum), a standard practice toanalyze the acoustic features of area-based articulatory models over a range of fre-quencies. The acoustic analysis of the 2.5D vocal tract model has been presentedin Section 5.4.5.3 Experimental SetupThe following section outlines the simulation parameters, modelling requirementsand conditions, and system specifications.5.3.1 Simulation ParametersFollowing are the parameters that remain fixed throughout the 2.5D vocal tractmodel simulations. During the spectra analysis for formants’ positions, we setthe boundary admittance coefficient to µ=0.005 at the vocal tract walls, whichcorresponds to the wall impedance of the vocal tract tissue Zn = 83666 kg/m2s[115]. However, while demonstrating the tuning of formants’ bandwidths, we ad-just µ ∈ (0,1) to achieve the desired range.76Symbol Description Valuec sound speed 350m/sρ air density 1.140kg/m3η Courant’s Number (CFL) 0.707max(σ∆t) maximum attenuation coefficient at PML layers 0.5µ boundary admittance coefficient 0.005− Number of PML layers 6Table 5.3: Simulation parameters for 2.5D FDTD vocal tract acoustic model.5.3.2 Excitation FunctionFor acoustic analysis and speech production, articulatory models require a sourceexcitation function just above the glottal end to input acoustic energy into the vo-cal tract tube [81]. We are only characterizing the acoustic behaviour of vocal tracttubes through transfer function analysis method for SA-Analytical, SA-Vowel andSA-Radiation. Hence, in our study, we inject a band-passed velocity pulse (shownin Figure 5.4) having a frequency range 2Hz-22kHz near the glottis. The band-pass excitation function is created using the windowed-sinc filter as it can be usedto separate a band of frequency from others and remains very stable [100]. Al-ternatively, a Gaussian pulse can also be implemented as shown by Takemoto etal. [116] to estimate the frequency response of a tube. Table 5.4 shows the fre-quency response of a vocal tract model (vowel /A/) remains the same irrespectiveof the band-pass and Gaussian pulses as the excitation functions. However, we usea vocal fold biomechanical model at the glottal end as the excitation function toillustrate an articulatory speech synthesizer (FT-Coupling).Excitation Function F1 F2 F3 F4Band-pass 700 1040 3020 4120Gaussian 700 1040 3020 4120Table 5.4: Positions of the first four formants, Fi (Hz) for vowel /A/, extractedfrom its frequency response of band-pass and Gaussian velocity pulses.These input signals are applied separately as an excitation function to the2.5D vocal tract model.770 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9Time [second]020406080100120140160180200Amplitude [m/s]0 5 10 15 20Frequency [kHz]102030405060708090100Amplitude [dB]Figure 5.4: (left) Input band-pass velocity pulse (right) and its frequencycharacteristics. In the frequency domain the spectral amplitude remainsconstant from 2Hz-22kHz.(a) (b)Figure 5.5: Snap shot of acoustic wave propagation at t = 15ms for two dif-ferent mouth-end boundary conditions (a) open-end boundary condition(b) free field radiation condition.785.3.3 Mouth-End ConditionsWhile speaking, the acoustic pressure waves terminate lips and propagate outward.The outward wave propagation (radiation) at the mouth-end contributes to acousticenergy losses in the vocal tract. It is one of the essential dissipation mechanismsin voice production. Hence, a vocal tract radiation model influences the transferfunction by changing its formant frequencies and bandwidths. In our study, weconsider two different situations at the vocal tract mouth-exit cross-section.In the first case, the computational domain ends at the mouth exit. For that,we close the vocal tract exit cross-section by a baffle and enforce the Dirichletboundary condition (open-end condition - Figure 5.5(a)) on it (p = 0) to emulateopen-end termination. In the second case, we consider free-field radiation by al-lowing acoustic waves to emanate from the mouth aperture (radiation condition -Figure 5.5(b)). We adopt the open-end occurrence to analyze the acoustic charac-teristics of SA-Analytical and SA-Vowel. We also use the open-end terminationcondition while demonstrating an end-to-end articulatory speech synthesizer (FT-Coupling). We only impose the free field radiation condition while investigatingSA-Radiation.It might be challenging to build a vocal tract radiation model that can pre-cisely capture all the radiation effects, considering a realistic human head geome-try. However, it has already been demonstrated that the inclusion of a spherical baf-fle (a simplified approximation of the actual head geometry) around the 3D vocaltract model offers promising results for the vowel production [9]. As the 2D modeluses a mid-sagittal slice of a 3D vocal tract, the spherical baffle can be replaced bya circular baffle for the 2.5D radiation model as shown in Figure 5.5(b). Moreover,to illustrate acoustic wave propagation emanating from the mouth toward infinity,the computational domain can be extended out of the vocal tract using perfectlymatched layers (PMLs). PMLs help in absorbing outward acoustic pressure wavesand eliminating spurious wave reflection at the domain boundaries.We create a fixed-sized circular baffle with a diameter of 0.20m around the2.5D vocal tract contour inside the simulation domain. The baffle walls have beenassumed as hard walls without any boundary impedance (µ = 0). The circularbaffle center has been adjusted to fit the entire vocal tract contour while connect-79ing its exit cross-section. We close the circular baffle near the glottal end abruptlyto reduce the domain size, representing the backside of the head geometry. Suchsemi-circular baffle geometry has negligible influence on the acoustic output of theradiation model [8]. However, a detailed analysis of the complete head geometry’srole in the speech production mechanism needs to be done. Nevertheless, we caneasily obtain a full circle as baffle shape by extending the computational domainsize at the expense of simulation run-time. Inside the 2.5D scheme, the region en-closed by the vocal tract and baffle boundary is set to air cells. To absorb outgoingwaves and approximate free field space, we add 6 PMLs at all sides of the domainboundary. The 6-layers of PML are sufficient to remove any audible acoustic wavereflection, with PML absorptivity σ decreasing linearly from 0.5/∆t at the domainedges to 0 inside. Outside the PML absorbing layers, σ = 0 (no PML absorption).5.3.4 Vocal Tract Transfer FunctionsWe place a virtual microphone at a grid point to sample the propagating pressurewaves (p0(n)) at each time step. For the open-end condition, we place the micro-phone 3mm inside the mouth-exit (Figure 5.6), whereas for the radiation condition,the microphone is positioned 3mm outside the mouth-exit (Figure 5.7). For eachexperiment, we run the computation simulation to generate an acoustic output of50 ms. Then we apply the Fast Fourier transform (FFT) to the recorded pressuresamples (p0(n)) and excitation velocity samples (ve(n)) at glottal cross-section toobtain the vocal tract transfer function (VTTF = H( f )) and formants (local max-ima in the transfer function or impulse response curve) in the following manner:H( f ) =p0( f )ve( f )(5.1)where p0( f ) and ve( f ) stands for the FFT of p0(n) and ve(n) respectively. Weonly analyze the first eight formants in VTTFs (up to 10 kHz approximately). Be-yond that limit, the acoustic behaviour of the vocal tract becomes intricate due tothe onset of higher-order nonplanar acoustic modes. To attend those modes, the2.5D model requires the inclusion of vocal tract side branches and complex cross-sectional geometry.8020 40 60 80 100 120 140 160 180 200 220510152025303540Figure 5.6: Vocal tract shape for vowel /A/ with open-end boundary condi-tion. The black dot at the right hand side represents the position ofvirtual microphone/listener (3mm inside the mouth aperture).Figure 5.7: Vocal tract shape for vowel /A/ with radiation condition. The righthand side close-up image shows the detail of 2D scheme at the mouthaperture. The black dot in it represents the position of virtual micro-phone/listener (3mm outside the mouth aperture).815.3.5 System SpecificationsWe use a serial machine to implement the 2.5D FDTD model and run the acousticsimulations for various tube shapes. Table 5.5 lists out the system specificationsand the 2.5D FDTD modelling environment.Technical Specification DescriptionOperating System Windows 10 HomeProcessor Intel(R) Core(TM) i7-8750H CPU @ 2.20GHzSystem Type 64-bitInstalled Memory(RAM) 15.9 GB (usable)Modelling Environment MATLAB 2019aTable 5.5: System specifications for 2.5D FDTD numerical simulationsClosedendOpen endLAA1 = 0.91 cm2 A2 = 1.31 cm2L1 = 17.07 cm Table C.1 Table C.1L2 = 50.16 cm Table C.2 Table C.2(a) (b)Figure 5.8: (a) Illustration of a one-end closed cylindrical tube with funda-mental mode (b) Tabular representation of cylindrical tubes with differ-ent lengths of Li (in cm) and cross-sectional areas of Ai(in cm2).825.4 Results5.4.1 Acoustic analysis for the 1D one-end closed lossless cylindricaltubeDue to the external turbulence, the air column inside a cylindrical tube vibrates andforms standing wave patterns. Consequently, this phenomenon generates pressurenodes or acoustic formants. In a 1D one-end closed cylindrical tube (shown inFigure 5.8 (a)), these formant frequencies can be extracted analytically using thefollowing rule:fn =nc4L(5.2)where c is the speed of sound (c = 350m/s) and L is the tube length. n = 1,3,5...(2n+1), as a one-end closed cylindrical tube resonate at a fundamental frequency andodd harmonics. It must be noted that Equation 5.2 does not take into account tubecross-sectional areas. Hence, the derived formants’ positions using the analyticalmethod should be considered as the close approximation of actual ones.We study four one-end closed lossless cylindrical tubes (µ = 0) having differ-ent lengths (L) and cross-sectional areas (A) inside the 2.5D computational domain.Each cylindrical tube model consists of 44 ducts, and all the ducts in a tube have thesame length (L/44) and cross-sections (A). We have chosen the tube lengths andtheir cross-sections empirically. However, we avoid simulating long cylindricaltubes having very narrow constrictions as they are considerably more susceptibleto numerical noises. The sampling rate for the computational simulation has beenfixed to 661.5kHz (2.5Dlow, i.e., 2.5D simulation at a low grid resolution). Weapply the open-end boundary condition at the exit cross-section of the tube. Fig-ure 5.8 (b) presents all the test scenarios that we have followed for our experimentsand maps their corresponding tabular results.83F1 F2 F3 F4 F5 F6 F7 F8Formants-1-0.500.511.52 Formant Deviations: F i/Fi (%)A1=0.91 cm2A2=1.31 cm2Figure 5.9: Formant deviations ∆Fi/Fi (%) of the 2.5Dlow simulations of one-end closed cylindrical tube models (L1 = 17.07cm) with respect to theiranalytical solution. The percentage of deviations have been derived us-ing the formant positions from Table C.1.F1 F2 F3 F4 F5 F6 F7 F8Formants-1-0.500.511.522.53 Formant Deviations: F i/Fi (%)A1=0.91 cm2A2=1.31 cm2Figure 5.10: Formant deviations ∆Fi/Fi (%) of the 2.5Dlow simulations ofone-end closed cylindrical tube models (L1 = 50.16cm) with respectto their analytical solution. The percentage of deviations have beenderived using the formant positions from Table C.2.845.4.2 Acoustic Analysis For Open Static Vocal Tract ShapesWe have analyzed the simulation outputs of four open static vocal tract shapes(vowel sounds: /a/, /i/, /u/ and /æ/) to characterize the acoustic features of the2.5D FDTD scheme for complex tube geometries. We build the vocal tract modelinside the computational domain as discussed in Section 4.6. A band-pass excita-tion function is positioned near the glottis, and the open-end boundary conditionis enforced at the exit cross-section of our vocal tract model. Since the numericalcomputation ends at the mouth aperture, no PML layers are added to the domainboundaries. We have simulated the 2.5D vocal tract model with three different sam-pling rates ( fs = 1/∆t) and their corresponding grid resolutions: 2.5Dlow = 661.5kHz (∆s = 0.74mm), 2.5Dmid = 1764 kHz (∆s = 0.28mm), 2.5Dhigh = 2646 kHz(∆s = 0.18mm). The grid resolutions are set to fulfill the CFL stability condition.2.5D with reference to 1D vocal tract modelFirst, we compare the 2.5D simulation results to a 1D vocal tract model. 1D sys-tems can only simulate straight tubes with circular cross-sections, and they gener-ally produce a good match of the lower end of the speech spectrum. In this study,we use Kee’s time-domain-based 1D vocal tract model as our reference system[131]. The model has been implemented in the JASS simulation environment (AJava Audio Synthesis System), and it numerically (finite volume method) solvesthe linearized Navier-Stokes acoustic wave equation (closely related to Webster’shorn equation) inside a 1D domain. It linearly interpolates the original area func-tion to generate precise speech spectra. We set the sampling rate ( fs) of the 1Dsystem to 88.2 kHz during the numerical simulation. It has already been demon-strated that with this sampling rate, the 1D model produces good results for vowelsound simulations in the planar wave propagation regime (below 5 kHz) [131]. Weapply the open-end boundary condition at the mouth exit for the 1D model. Forboth the 2.5D and 1D systems, we use Story’s area function dataset to generatethe vocal tract geometry. However, these models follow different approaches forestimating vocal tract boundary losses.For a comparative study, Table 5.6 presents acoustic formants of the 1D and2.5Dlow vocal tract model for vowel sounds: /A/, /i/, /u/ and /æ/. We only analyze85the first three formants of transfer functions generated from the 1D model simula-tions as they are unreliable at higher formants (≥ 5kHz). We expect differences inthe formants’ positions between these two models due to their different modellingmethodology. Moreover, in contrast to the 2.5D model, the 1D vocal tract allowspressure wave propagation along the axial direction only.Hence, we first examine the 1D and 2.5Dlow systems’ acoustic spectra for ananalytical cylindrical tube model and measure their deviation from the analyticsolution (see Figure 5.11). Due to the regular geometry of a cylindrical tube, dis-cretization noises in the simulation environment remains minimum for both sys-tems. Thus, this analysis approximates the existing deviation between the 2.5Dand 1D systems due to the differences in the modelling approach.2.5D with reference to 2D and 3D vocal tract modelsNext, the extracted formant frequencies from the 2.5D FDTD simulation are ana-lyzed with reference to the 3D FEM and the existing 2D FDTD vocal tract models.For the 3D model, we directly utilize the simulation results, produced by Arnela forthe following vowels only: /A/, /i/ and /u/ [12]. The 3D model uses a sampling rateof fs = 2000 kHz. It also imposes a boundary losses condition at the vocal tractwalls, similar to the 2D/2.5D vocal tract model. The 3D computational domainhas been meshed using tetrahedral elements, and the mesh size changes across thedomain by having high resolution inside the vocal tract boundary (h = 0.001 m)and low resolutions (h≥ 0.001 m) outside of it.For the vowel sound /æ/, we compare acoustic outputs between the 2.5D and2D models only. For the 2D model, we only consider the higher grid resolution(2Dhigh) since it is more reliable, as demonstrated earlier. The 3D FEM and 2DFDTD models consider vocal tract geometry as a straight tube with circular cross-sections. Table 5.7 to Table 5.10 produce formants for the 3D, 2.5D, 2D and 1Dvocal tract models. We consider Marc’s 3D FEM vocal tract as the reference modelfor the comparative analysis of acoustic formants. In general, if the percentage offormant deviations lie between 4%–5% between these models, we could considerthem perceptually not relevant [9, 46]. Therefore, we use this deviation range todefine models’ accuracy during the spectral analysis of the synthesized speech.86F1 F2 F3Formants-4-202468 Formant Deviations: F i/Fi (%)1D2.5DFigure 5.11: Formant deviations ∆Fi/Fi (%) of the 1D and 2.5Dlow simu-lation of a one-end closed cylindrical tube model (L1 = 17.07cm,A1 = 0.91cm2) with respect to its analytical solution. The percent-age of deviations have been derived using the formant positions fromTable C.3.Vowels⇒ /A/ /i/ /u/ /æ/Formants ⇓ 1D 2.5D 1D 2.5D 1D 2.5D 1D 2.5DF17375.2%70029714.2%2602714.2%2608392.3%820F2991−4.7% 10402089-3.2%2160573-20.4%7201712-0.4%1720F32898−4.0% 30202908−4.9% 306023994.3%23002383-3.1%2460Table 5.6: Formant locations Fi (Hz) of the 1D and 2.5Dlow vocal tract modelfor vowel sounds - /A/, /i/, /u/ and /æ/. ∆Fi/Fi (%) presents the percentageof formant deviations of Kee’s 1D vocal tract model with respect to the2.5Dlow model.87Formants 3D 2.5Dlow 2Dhigh 1DF1 696 700 (0.57%) 720 (3.44%) 737 (5.89%)F2 1068 1040 (-2.62%) 1040 (-2.62%) 991 (-7.20%)F3 3031 3020 (-0.36%) 3020 (-0.36%) 2898 (-4.3%)F4 4064 4120 (1.37%) 4000 (-1.57%) NAF5 5039 4980 (-1.17%) 4940 (-1.96%) NAF6 5703 5680 (-0.40%) 5600 (-1.80%) NAF7 7065 7060 (-0.07%) 6860 (-2.90%) NAF8 7595 7620 (0.32%) 7440 (-2.04%) NATable 5.7: Formant locations Fi (Hz) and formant deviations ∆Fi/Fi (%) of2.5Dlow, 2Dhigh and 1D simulations of vowel /A/ with respect to the 3DFEM vocal tract model.F1 F2 F3 F4 F5 F6 F7 F8Formants-3-2.5-2-1.5-1-0.500.511.52 Formant Deviations: F i/Fi (%)2.5Dlow2.5Dmid2.5DhighFigure 5.12: Formant deviations ∆Fi/Fi (%) of the 2.5D simulations with dif-ferent grid resolutions for vowel /A/. The percentage of deviations havebeen derived using the formant positions from Table C.4.88Formants 3D 2.5Dlow 2Dhigh 1DF1 263 260 (-1.14%) 260 (-1.14%) 297 (12.92%)F2 2111 2160 (2.32%) 2120 (0.42%) 2089 (-1.04%)F3 3010 3060 (1.66%) 3000 (-0.33%) 2908 (-3.38%)F4 4138 4220 (1.98%) 4120 (-0.43%) NAF5 5019 5020 (0.01%) 4920 (-1.97%) NAF6 5753 5740 (-0.22%) 5560 (-3.35%) NAF7 6575 6480 (-1.44%) 6440 (-2.05%) NAF8 7645 7680 (0.45%) 7600 (-0.58%) NATable 5.8: Formant locations Fi (Hz) and formant deviations ∆Fi/Fi (%) of2.5Dlow, 2Dhigh and 1D simulations of vowel /i/ with respect to the 3DFEM vocal tract model.F1 F2 F3 F4 F5 F6 F7 F8Formants-3-2-10123 Formant Deviations: F i/Fi (%)2.5Dlow2.5Dmid2.5DhighFigure 5.13: Formant deviations ∆Fi/Fi (%) of the 2.5D simulations with dif-ferent grid resolutions for vowel /i/. The percentage of deviations havebeen derived using the formant positions from Table C.5.89Formants 3D 2.5Dlow 2Dhigh 1DF1 259 260 (0.38%) 260 (0.38%) 271 (4.63%)F2 757 720 (-4.88%) 700 (-7.52%) 573 (-24.30%)F3 2264 2300 (1.59%) 2260 (-0.17%) 2399 (5.9%)F4 3603 3620 (0.47%) 3600 (-0.08%) NAF5 4173 4280 (2.56%) 4120 (-1.27%) NAF6 5062 4960 (-2.01%) 4960 (-2.01%) NAF7 6126 6080 (-0.75%) 6100 (-0.42%) NAF8 6594 6620 (0.39%) 6520 (-1.12%) NATable 5.9: Formant locations Fi (Hz) and formant deviations ∆Fi/Fi (%) of2.5Dlow, 2Dhigh and 1D simulations of vowel /u/ with respect to the 3DFEM vocal tract model.F1 F2 F3 F4 F5 F6 F7 F8Formants-3-2-10123 Formant Deviations: F i/Fi (%)2.5Dlow2.5Dmid2.5DhighFigure 5.14: Formant deviations ∆Fi/Fi (%) of the 2.5D simulations with dif-ferent grid resolutions for vowel /u/. The percentage of deviations havebeen derived using the formant positions from Table C.6.90Formants 2Dhigh 2.5Dlow 1DF1 820 820 (0.0%) 839 (2.31%)F2 1700 1720 (1.17%) 1712 (0.70%)F3 2460 2460 (0.0%) 2383 (-3.13%)F4 4320 4380 (1.38%) NAF5 4760 4800 (0.84%) NAF6 5820 5880 (1.03%) NAF7 7080 7200 (1.69%) NAF8 8300 8440 (1.68%) NATable 5.10: Formant locations Fi (Hz) and formant deviations ∆Fi/Fi (%) of2.5Dlow and 1D simulations of vowel /æ/ with respect to the 2DhighFDTD vocal tract model.F1 F2 F3 F4 F5 F6 F7 F8Formants-3-2-10123 Formant Deviations: F i/Fi (%)2.5Dlow2.5Dmid2.5DhighFigure 5.15: Formant deviations ∆Fi/Fi (%) of the 2.5D simulations with dif-ferent grid resolutions for vowel /æ/. The percentage of deviations havebeen derived using the formant positions from Table C.7.915.4.3 Acoustic Analysis For Vocal Tract Radiation ModelsThe free-field radiation condition at the mouth aperture is crucial in emulating therealistic voice production mechanism. There are multiple strategies to emulatethe radiation effect. Two of the widely followed approaches are the impedanceload models and the free-field radiation condition. The impedance load model[126, 127, 151] uses the frequency-dependent radiation impedance in its acousticwave solver to approximate the radiation effect and terminates the computationaldomain at the vocal tract exit cross-section. Alternatively, the free-field radiationcondition [9, 114] extends the computational domain outside of the vocal tract andallows the pressure wave to propagate outward to imitate the vocal tract radiationeffects more naturally. Thus, we use the free-field condition for the proposed 2.5Dmodel to emulate the radiation mechanism.This section demonstrates the acoustic characteristics of the 2.5D vocal tractradiation model by including a simplified head geometry (circular baffle) and PMLlayers at the domain boundaries, as discussed in Section 5.3.3. However, the inclu-sion of a circular baffle and PML layers increases the computational domain size.Consequently, the simulation run-time for the vocal tract radiation model is muchhigher compare to a non-radiation model (open-end boundary condition). Hence,we only study the radiation characteristics of the 2.5D model at a low grid resolu-tion (∆s = 0.74 mm). Nonetheless, to provide a complete overview of the influenceof the free-field radiation condition on the computational cost, Table 5.15 approxi-mates the simulation run-time for all the grid resolutions (low, mid and high).Similar to the open-end boundary condition, we approximate the vocal tractgeometry as a straight tube with circular cross-sections for the radiation condi-tion. We derive the acoustic formants using the proposed 2.5D radiation modelfor vowel sounds /A/, /i/, /u/ and /æ/ and then compare them with reference tothe 3D FEM vocal tract radiation model [9] and the formants estimated from therecorded speech signals1 [109]. For the 3D FEM model, we have the simulationresults for the vowel sound /u/ only. The 3D radiation model uses the same mod-elling parameters as discussed for the open-end boundary condition. However, itincorporates a spherical baffle around the bent vocal tract geometry with circu-1For the radiation model, we only match the first three formants as per the data availability.92lar cross-sections. The spectral data from the recorded speech has already beenextracted using PRAAT’s formant analysis module and provided here [109]. Thetables below present the simulation results (formants’ positions) for the 2.5D vocaltract radiation model2. Table 5.11 presents a comparative study of acoustic for-mants for the 2.5D FDTD radiation model with reference to the 3D FEM model.Table 5.12 shows the acoustic formants from the 2.5D models and the recordedspeech signals.Vowel⇒ /u/Formants ⇓ 3D 2.5DF1 256 260 (1.56%)F2 729 720 (-1.23%)F3 2275 2180 (-4.17%)Table 5.11: Formant locations Fi (Hz) of the 3D and 2.5Dlow vocal tract radia-tion model for the vowel sound - /u/. ∆Fi/Fi (%) presents the percentageof formant deviations of the 2.5Dlow vocal tract model with respect tothe 3D FEM model.Vowel F1 F2 F3/A/715 (±9)6401180 (±31)10202691 (±21)3000/i/295 (±17)2601923 (±28)20802797 (±60)2960/u/430 (±5)260917(±14)7202087(±14)2180/æ/765 (±15)8001724 (±22)16002501 (±59)2420Table 5.12: Formants derived Fi (Hz) using the 2.5D model simulations. Av-eraged formants Fi(±s.d) across several seconds of recorded speech ut-terances and s.d. is the standard deviation.2Colormap: Acoustic formants from the recorded speech utterances, Acoustic formants from the2.5D vocal tract radiation model simulations.935.4.4 Time Complexity AnalysisThe computational complexity is one of the most significant factors in determin-ing an articulatory model’s feasibility as a real-time speech synthesizer. It can beanalyzed by approximating the simulation run-time or time-complexity of the nu-merical algorithm. The FDTD method follows a time marching algorithm (see Ap-pendix D) to update grid cells’ pressure and velocity components at each timestep.Hence, the time complexity of the FDTD algorithm is dependent upon the numberof grid points or the computational domain size.The time complexity of the FDTD computational simulation can be defined asO(n2), considering a 2D domain of size (n× n). As both the 2D and 2.5D wavesolvers use a 2D numerical domain, their time complexity remains the same. Dueto the quadratic nature of the algorithm, the computational costs of these acousticwave solvers are high. However, during the simulation, the wave solver updatesacoustic components at each grid point using its neighbouring cells’ information.Therefore, the algorithm can be mapped to parallel architectures, like GPUs, andaccelerate the computation. Zappi et al. [146] have already demonstrated a fastparallel 2D FDTD GPU wave solver with quasi-real-time simulation performance(see Table 5.13). Our proposed 2.5D FDTD scheme extends Zappi’s approachwithout compromising the parallel nature of the algorithm. Hence, it is safe toassume that the optimized parallel implementation of the 2.5D running on a moderngraphics card will attain close to the 2D simulation run-time for a fixed domainsize.Domain Size Simulation Run-time (Approx.)(341×43) low 00 : 00 : 03(690×97) mid 00 : 00 : 12(1049×146) high 00 : 00 : 29Table 5.13: Simulation run-time (in hh:mm:ss) for the 2D FDTD GPU wavesolver with different domain sizes (grid resolutions) to generate 50 msof audio (open-end boundary condition). Reproduced with permissionfrom Victor Zappi [146].The tables below present the approximate simulation run-times for the 2.5D94FDTD scheme with different domain sizes. We run each simulation to produce50 ms of audio using a serial machine (specifications are provided in Table 5.5)for both the open-end and free field radiation boundary conditions. As both the2D and 2.5D FDTD schemes use a 2D numerical domain for acoustic simulations,their time complexity will remain approximately identical.Domain Size Simulation Run-time (Approx.)(232×43) low 00 : 05 : 38(612×109) mid 01 : 37 : 00(917×161) high 06 : 01 : 40Table 5.14: Simulation run-time (in hh:mm:ss) for the 2.5D FDTD wavesolver with different domain sizes (grid resolutions) to generate 50 msof audio (open-end boundary condition).Domain Size Simulation Run-time (Approx.)(253×283) low 00 : 40 : 20(647×729) mid 12 : 58 : 20(962×1085) high 43 : 03 : 00Table 5.15: Simulation run-time (in hh:mm:ss) for the 2.5D FDTD wavesolver with different domain sizes (grid resolutions) to generate 50 msof audio (free field radiation boundary condition).5.5 DiscussionThis section looks at the 2.5D FDTD simulation results in the context of observa-tions noticed while conducting experiments for each step of our road map. Imagingvocal tract geometrical details and mapping them to their corresponding acousticfeatures are complicated. Hence, the available datasets of vocal tract area functionsand their spectral analysis with various speech utterances are very limited. Thus,most of the comparisons presented here to the other vocal tract models and experi-mental data in the literature are purely qualitative, and they only explore the 2.5Dvocal tract model’s basic capabilities as possible pathways for future developments.955.5.1 SA-Analytical: One-end Closed Cylindrical Tube Model• In the case of a 1D one-end closed cylindrical tube, the analytical solutionassumes an ideal tube of zero diameters for the planar wave propagationcondition. Therefore the computational tube model only matches the ana-lytic solution as long as it allows the plane wave propagation inside the tubegeometry. This condition can be reached by setting the duct cross-sectionssufficiently small compared to the wavelength at the frequencies of interest.Thus, we chose a small cross-section (A) for all of our computational tubemodels.• Figure 5.9 and 5.10 show for two different tube lengths, the derived formantsfrom the 2.5Dlow simulations are identical to their analytic solutions. In mostcases, the absolute error percentages (∆Fi/Fi) remain below 2%.• The existing formants’ deviations between the numerical and analytic solu-tion are most likely a consequence of Equation 5.2, which disregard the tubediameter and the numerical noises entirely.• As per Equation 5.2, the simulation results show the formants’ locationsfor the computational tube models having the same length (L) but differentcross-sections (A) are almost identical to each other with few exceptions.From the above analysis, we can conclude that our 2.5D FDTD acoustic modeldoes an excellent job replicating the pressure wave propagation inside a 1D one-end closed cylindrical tube model and matching its formants’ frequencies to theanalytic solution at the plane wave propagation regime.5.5.2 SA-Vowel: Open Static Vocal Tract Shapes For Vowel SoundsWe compare the 2.5D simulation results of cardinal vowels with the existing 1D,2D and 3D vocal tract models. The current 2.5D model is a simplified representa-tion (straight tube with circular cross-sections) of the actual vocal tract geometry.Therefore, it completely ignores the significance of the intricate vocal tract shape,side branches, mouth cross-sections, the source model, etc., on acoustic outputs.Thus, at this point, we avoid attempting a comparative study of our proposed modelwith the original spectral features derived from the recorded speech signals.962.5D vs 1D:• First, we carried out a numerical assessment of Kee’s 1D model by compar-ing its simulation results to the analytical tube model. Figure 5.11 shows, fora cylindrical tube of L1 = 17.07cm and A1 = 0.91cm2, the 1D system reason-ably produces comparable results with few exceptions. The transfer functionhas a higher error percentage (∆Fi/Fi) for F1 compare to F2 and F3. Whereas,in the case of 2.5Dlow simulation, the error percentages for all the three for-mants are below 2% in relation to their analytic solutions (see Figure 5.11).Besides, due to the dimensionality reduction, the percentage of formant de-viations between the 1D and 2.5D simulations remain within 4%−7% (canbe derived using Table C.3) for the cylindrical tube model. Thus, we alsoexpect approximately similar discrepancies in acoustic spectra of the 2.5Dand 1D model for vowel simulations.• Table 5.6 shows that both the 2.5D and 1D models are in good agreement inthe case of vowel sound simulations with deviations that tend to stay between4%−5% for the lower formants. Vowels /A/ and /æ/ display the best matches.The F1 and F2 of /i/ and /u/ are characterized by larger shifts of the formants’positions, respectively. However, these differences are apparent consideringthat the 1D model simplifies the complex 3D vocal tract geometry to a 1Darea function and allows wave propagation only along the axial direction.• At a higher acoustic spectrum (beyond 5 kHz), the wavelength of acousticwaves becomes of the same order of size as the transverse dimension of thevocal tract (tube aperture). Consequently, the transverse acoustic modes areintroduced, and the requirement for the planar wave propagation conditiondoes not hold up at the high-frequency range for the 1D system. In thiscase, the 1D model displays a larger deviation as we are outside of the planewave propagation regime. Table 5.7-5.9 show Kee’s 1D system produces agood match of the first three formants in the transfer function to the proposed2.5D and 3D models with few exceptions. However, beyond the plane waveregime (i.e., 5 kHz), the transfer function of the 1D vocal tract does notmatch at all to the proposed 2.5D model (see Figure 5.16).970 1000 2000 3000 4000 5000 6000 7000 8000Frequency [Hz]00.20.40.60.811.2Normalized Amplitude [dB]2.5D Model1D ModelFigure 5.16: Vocal tract transfer functions (VTTFs) for vowel /A/ generatedfrom the computational simulation of (red dotted line) Kee’s 1D model(blue line) 2.5D FDTD model.Figure 5.17: Transformed vocal tract area function (vowel sound /A/) for the2D and 2.5D vocal tract model with reference to the original area func-tion (3D). The highlighted area inside the (dotted) blue box shows areafunctions of different vocal tract models at a narrow constriction.982.5D vs 2D vs 3D:• Due to the area function transformation, the 2D vocal tract model loses itsboundary complexity at narrow constrictions. However, the 2.5D modellumps the effects of nonplanar wave propagation, which means the adjacenttube segments for both the 2.5D and 3D model can be characterized with thesame wave reflection coefficients. Hence, we do not require such area func-tion transformation for the 2.5D model. Figure 5.17 shows that in contrastto the 2D model, the 2.5D vocal tract’s cross-sectional areas remain approx-imately the same as the original 3D vocal tract for narrow constrictions (i.e.,area of the tube cross-sections at 4−8 cm from the glottal end).• We observed a higher percentage of formants’ deviations (> 9%) betweenthe 2Dlow (WS) and 3D model for the vowel sound /A/ simulation (Table 3.1).However, even at a low grid resolution, the 2.5Dlow simulation results (Ta-ble 5.7-5.9) are relatively comparable to the high-quality 3D FEM vocal tractmodel for vowel sounds /A/, /i/ and /u/ with a percentage of deviations below3%. We only noticed a slightly higher shift of the second formant (F2 > 4%)for the vowel /u/ in the case of the 2.5Dlow simulation.• As discussed in Section 3.1.3, a higher grid resolution is always preferredfor the 2D vocal tract simulation to retrieve the boundary complexity at theexpense of simulation run-time. Therefore, the acoustic spectrum of 2Dhighfor vowel /A/ matches well with the 3D FEM model as shown in Table 3.1.Whereas in the case of 2.5Dlow simulation, the model retains its geometryintact even at a lower grid resolution due to the vocal tract depth mapping,consequently eliminating the grid resolution tuning requirement. Thus, theformants’ locations (peaks) in the acoustic spectra of the 2.5Dlow and 2Dhighare nearly identical as shown in Figure 5.18 for all four vowel sounds.• A higher grid resolution does not produce any substantial changes in thevocal tract geometry for the 2.5D model. Therefore, we did not observe anysignificant differences in the formants’ locations for different spatial gridresolutions (2.5Dlow, 2.5Dmid and 2.5Dhigh) during the all four vowel soundsimulations (see Figure 5.12-5.15). Consequently, we can safely conclude991 2 3 4 5 6 7 8 9 10Frequency [kHz]050100Amplitude [dB]vowel /a/1 2 3 4 5 6 7 8 9 10Frequency [kHz]050100Amplitude [dB]vowel /i/1 2 3 4 5 6 7 8 9 10Frequency [kHz]050100Amplitude [dB]vowel /u/1 2 3 4 5 6 7 8 9 10Frequency [kHz]050100Amplitude [dB]vowel /æ/Figure 5.18: Vocal tract transfer functions for vowel sounds (/A/, /i/, /u/ and/æ/) generated from the computational simulation of (1) 2Dhigh (reddotted lines) (2) 2.5Dlow (blue lines).100that the 2.5Dlow simulations of open static vocal tract shapes produce preciseestimation of acoustic wave propagation.• We observed the formant deviations (∆Fi/Fi) of the 2.5Dlow model at higherformants (F6, F7 and F8) are comparatively less than the 2Dhigh model forvowel sounds (/A/, /i/ and /u/ ). A possible explanation for these differencescould be the inclusion of vocal tract depth in the 2D acoustic wave solver.This additional geometrical feature allows capturing nonplanar acoustic waveinteractions, consequently excites the higher-order acoustic modes.The above analysis shows the acoustic characteristics (F1-F3 in VTTFs) ofthe 2.5D and 1D vocal tract models are reasonably similar with few discrepanciesunder the planar wave cut-off frequency. Moreover, the 2.5Dlow also produces rel-atively identical formants’ positions with reference to the 2Dhigh and 3D vocal tractmodels for the four English vowel sounds even at a low grid resolution. Therefore,the 2.5D vocal tract model is computationally cheaper for precise acoustic simu-lation. It is to be noted that we can not endorse the accuracy of the 2.5D schemeat this point as we only simulate four English vowel sounds. An extensive studycan be carried out to characterize the 2.5D model for static speech sounds in otherlanguages.1015.5.3 SA-Radiation: 2.5D Vocal Tract Radiation Model For openstatic vocal tract shapesThe vocal tract radiation model is arguably one of the most critical acoustic en-ergy dissipation mechanism that influences the spectral features of the synthesizedacoustic outputs. Figure 5.19 shows the differences in the VTTFs of the 2.5Dradiation and non-radiation model simulations for vowel sound /A/. The transferfunctions clearly show the effects of the mouth radiation boundary condition as theformant position and bandwidths change, especially beyond 5 kHz. Thus, a vocaltract model that precisely emulates the radiation loss mechanism is an importantfactor for the articulatory speech synthesis models.500 1000 1500 2000 2500 3000 3500 4000 4500 5000 5500Frequency [Hz]20406080Amplitude [dB]Figure 5.19: The transfer function of vowel sound /A/ for the (a) 2.5D non-radiation [blue line] and (b) 2.5 radiation vocal tract model [red dottedline]. We run both of these simulations at a fixed boundary admittancecoefficient of µ = 0.010.• Kee’s 1D vocal tract model has not been characterized for the mouth radia-tion condition using Story’s area functions previously. Hence, in this study,we have not compared the acoustic features of the 2.5D radiation model withthe 1D vocal tract model. As demonstrated in Table 3.1, the 2D FDTD modelrequires higher spatial resolution to produce accurate VTTFs for vowel sim-ulations. Consequently, the 2Dhigh vocal tract model is computationally veryexpensive with the mouth radiation boundary condition due to the inclusionof the head geometry and PMLs. Thus, the 2D radiation model is not a real-istic approach for synthesizing speech outputs. Therefore, in this study, weonly examine the acoustic outputs of the 2.5D radiation model with referenceto the recorded speech signals and Arnela’s 3D FEM radiation model.102• Table 5.11 shows a comparative study of formants’ positions of the 3D and2.5D radiation model for vowel sound /u/. Vowel /u/ displayed a good matchof the first two formants, but we observed a larger shift in the third formant.Since the vocal tract geometry plays a vital role in determining the acousticcharacteristics of a radiation model [9], these differences are apparent. The3D model approximates the vocal tract geometry as a bent tube inside thecomputational domain, whereas the 2.5D model considers a straight tube.• Table 5.12 shows the acoustic formants of vowel sounds (/A/, /i/, /u/ and/æ/) generated from the 2.5D radiation model and the recorded speech sig-nals. We observed noticeable differences in acoustic formants between therecorded speech and 2.5D simulation-generated outputs. However, these ex-isting deviations in VTTFs are obvious since our proposed 2.5D model over-simplifies the realistic vocal tract cross-section and head geometry. It rep-resents the vocal tract as a straight tube with circular cross-sections insidea circular baffle. Arnela et al. [9] have already demonstrated that the 3DFEM radiation model’s spectral characteristics match well with a realisticvocal tract (realistic head with irregular VT cross-sections) after consideringa bent vocal tract tube with elliptical cross-sections coupled with lips. Asimilar approach can also be implemented for the 2.5D vocal tract radiationmodel to mitigate the current formants’ positional differences.5.5.4 Formant Bandwidth AnalysisThe vocal tract geometry and wall impedance, glottal input, physical properties ofthe air and the mouth-end boundary conditions influence the position, bandwidthand amplitude of the formant frequencies [47, 52]. Therefore, computational vocaltract models demand accurate modelling of these features to generate better acous-tic outputs. Previously, we presented the formant positions for the 2.5D simulationof open static vocal tract shapes. However, their bandwidths also play an essentialrole in determining the naturalness of vowel sounds. As the current 2.5D vocaltract model partly incorporates the wall impedance, the simulation-generated for-mants’ bandwidths do not match the actual bandwidth of natural speech, which isabout 80-120Hz [137]. However, with the inclusion of mouth radiation condition,103we observed some improvements in formants’ bandwidth, as shown in Figure 5.19.The formants’ bandwidths in transfer function can be tuned relatively easily bychanging the boundary admittance coefficient (µ) as shown in Figure 5.20; conse-quently, the wall impedance (see Equation 4.11 and 4.12). However, that compro-mises the relationship of speech being driven by measurable parameters from theactual vocal tract. Therefore, such methods should always be avoided by preciselymodelling the vocal tract energy loss mechanisms and other geometrical features.550 600 650 700 750 800 850Frequency [Hz]50556065707580Amplitude [dB]=0.005=0.015=0.025Figure 5.20: Bandwidth variations (BW1) of the first formant F1 for the vowel/A/ by the 2.5D FDTD simulation using the open-end boundary condi-tion.1042.2. Lumped-Element Models of the Vocal FoldsFigure 2.7: Coronal (left) and superior (right) view of the vocal folds withquantities highlighted. Ishizaka et al, [64]2.2.2 Multi-Mass ModelMulti-Mass models added extra mass elements to each vocal fold to over-come some of the drawbacks of the single-mass model. The classic two-massmodel introduced by Ishizaka and Flanagan [64] added a second mass su-perior to the one-mass model (shown in fig.2.6), while restricting motion tomedial-lateral translation. While the airflow fluid loading occurs only on thelower mass, this configuration allows for the experimentally observed phasedifference between the inferior and superior edges of the vocal folds. Dueto the importance of this model in terms of understanding the field of vocalfold modeling better, let us briefly look at the mechanics behind this model.Two-Mass ModelFigure 2.7, shows the different quantities involved in the two-mass vocalfold structure in the coronal view. The displacements of the two vocalfolds define the vocal fold shape, and as a consequence, the flow at everytime instant. The lateral displacements x1 and x2 are the solution to thedifferential equation:m1x¨1 + r1x˙1 + s1 + kc(x1 − x2) = f1 (2.1)m2x¨2 + r2x˙2 + s2 + kc(x2 − x1) = f2 (2.2)12Figure 5.21: Coronal (l r i ht) iew of the two mas approx-imation of the vocal cords [61].5.6 Vocal Fold ModelTo analyze our vocal tract model’s feasibility for an articulatory speech synthe-sizer, we create an implementation of the two-mass vocal fold model [62] in theliteratur . It has already been demonstrated that the two-mass vocal fold modelcan sustain the self-oscillatory characteristics of human vocal folds and produceglottal waveform to a reasonable xtent. Moreover, Sodhi et al. [101], and Kee etal. [131] have demonstrated the coupling of the 1D vocal tracts to the two-massmodel. Likewise, Zappi et al. [146] have illustrated the coupling of the 2D FDTDvocal tract to the two-mass model to produce synthesized vowel sounds: /A/, /i/ and/u/. Hence, the choice of this model is logical in the context of the 2.5D vocal tractsimulation for an articul tory model.5.6.1 Two-Mass Vocal Fold ModelThe two-mass vocal fold model is based on lumped elements. It approximates eachside of the vocal cord as a structure of two mass elements (i.e., superior and inferiormass) connected by damper-spring systems as shown in Figure 5.21. The modelallows the displacement of point masses in the medial-lateral direction only andintroduces a phase difference between the superior and inferior edges of the vocalfolds. The two-mass model assumes the phase differences in the vocal cord edgesallow energy transfer from the airflow to the vocal cord, whereas the lung pressure105drives the vocal cord vibration. The displacement of the two vocal folds definesthe vocal fold shape and, consequently, the flow at every time instant. Hence, themathematical model of the system can be represented with the following second-order differential equations, which basically solve the lateral displacements x1 andx2 of mass elements in a discretized manner at each timestep.m1x¨1 + r1x˙1 + s1 + kc(x1− x2) = f1 (5.3)m2x¨2 + r2x˙2 + s2 + kc(x2− x1) = f2 (5.4)where:m1 and m2 the superior and inferior masses of the vocal cord,r1 and r2 the damping resistances,s1 and s2 the restoring spring forces (k1x1, k2x2),x1 and x2 the lateral displacement from the rest position, x0 > 0kc the stiffness of the coupling spring,f1 and f2 the aerodynamic forces on the vocal fold masses (m1,m2).The two-mass vocal fold model then uses the well-known Bernoulli equation todescribe the one-dimensional glottal flow, which helps obtain the pressure distribu-tion (pi) along the glottal flow region and the volume velocity near the glottal end(ug). We implement the two-mass vocal fold model in the MATLAB environmentand set the vocal fold structural parameters as described here [61]. Our articulatorymodel assumes a standard lung impedance pressure of Ps = 686.46 Pa (7 cm H2O)as the subglottal pressure during the voiced speech utterances. For the vocal foldmodel, we set the air density of ρ = 1.14 kg/m3 and air viscosity coefficient ofµ = 1.86× 10−5 kg/(m.s) respectively. The sampling rate of the simulations hasbeen fixed to 661.5 kHz.1065.7 2.5D FDTD Vocal Tract’s Feasibility For AnArticulatory ModelThe central purpose of this section is to demonstrate an articulatory speech syn-thesizer using the 2.5D vocal tract model as one of the components (FT-Coupling).Basically, we couple the two-mass vocal fold to our vocal tract model as an acousticsource that injects a periodic volume velocity pulse into the vocal tract to gener-ate synthesized speech utterances. However, we do not carry out any numericalassessment of the synthesized audio outputs in this study. We have also not ex-amined the flow and structural interactions between the vocal fold and vocal tractmodels. Instead, the primary objective here is to illustrate the 2.5D vocal tractmodel’s feasibility for designing an articulatory speech synthesizer as an applica-tion. Hence, we only examine the model’s stability throughout the simulation byproducing synthesized audio outputs of English vowels which are studied in Sec-tion 5.4.5.7.1 Coupled vs Uncoupled Vocal FoldsWe consider two different conditions at the glottal end to drive the proposed 2.5Dvocal tract model. In the first case, we couple the two-mass vocal fold with thevocal tract model via the back-propagating supraglottal pressure (p1), which comesfrom the 2.5D simulation domain. The p1 is the time-varying pressure samples atthe laryngeal cavity (the first cross-section of the vocal tract tube). Therefore,we update the two-mass vocal fold parameters at each timestep and compute therespective volume velocity (ug) using p1. The volume velocity obtained from thevocal fold model is set at the glottal end as the excitation velocity for the 2.5Dvocal tract model.We carried out the same experiment with the uncoupled version of the two-mass vocal fold model in the second case, considering both the vocal fold and vocaltract models as independent components. This approach would imply consideringthe supraglottal pressure p1 as a fixed parameter and simulating the two-mass vocalfold model separately to obtain the glottal volume velocity (ug). Subsequently, thederived glottal volume velocity was used as the source excitation function for the2.5D vocal tract model to synthesize vowel sounds.1075.7.2 Synthesized Vowel Sound GenerationWe run simulations to produce vowel sounds (/A/, /i/, /u/ and /æ/) of 50 ms. A vir-tual microphone is placed 3mm inside the mouth-exit as a listener, and an open-endboundary condition is applied at the mouth-exit cross-section. The virtual micro-phone collects outward pressure samples at the listener end in an array. For vowelsounds, we lowpass the normalized output-pressure samples above 22kHz using asecond-order Butterworth filter. The feed-forward and feedback coefficients of thefilter are 0.0754×{1,2,1} and {−1.0875,0.3890} respectively [3]. To producea smooth audio output, we downsample the filtered output to match the computerspeaker’s sampling rate (44100 kHz). The synthesized vowel sounds for both thecoupled and uncoupled system can be found here3.5.7.3 Model AnalysisBoth the coupled and uncoupled articulatory models remain stable throughout thesimulation. Figure 5.24 and 5.25 show the snapshots of pressure wave propaga-tion inside the vocal tract geometry (vowel sound: /A/) at different time instants forboth of these models. We observed a significant amount of nonlinearity (see Fig-ure 5.22) in the output pressure waveform for the coupled 2.5D vocal tract modelcompared to the uncoupled one. Moreover, we also noticed higher negative peaksin the normalized pressure waveform for the coupled model. Similar character-istics were also observed while simulating the coupled 2D FDTD vocal tract andtwo-mass vocal fold model [132]. To further understand if this nonlinearity is afunction of the coupling or the simulation artifacts, we captured the glottal volumevelocity, the input acoustic energy to the vocal tract model. Figure 5.23 shows thedifferences in the volume velocity waveforms for both of these models in termsof nonlinearity and negative dips. Thus, we conclude the coupled model producesa certain degree of nonlinearity due to flow interactions between the 2.5D vocaltract and vocal fold model. These results confirm the feasibility of the 2.5D vocaltract model for an articulatory speech synthesizer and its ability in emulating thenonlinear characteristics in the coupling mechanism, as highlighted by Titze [123].3https://github.com/Debasishray19/Talking-Tube/tree/master/audioFiles1080 5 10 15 20 25 30 35 40 45 50Time [ms]-0.6-0.4-0.200.20.40.60.81P norm [Pa]5 10 15 20 25 30 35 40 45 50Time [ms]-1-0.500.51P norm [Pa]Figure 5.22: Normalized output pressure for the 2.5D FDTD simulation ofvowel sound /A/ with a (left) uncoupled VF model (right) coupled VFmodel.5 10 15 20 25 30 35 40 45 50Time [ms]00.20.40.60.81Normalized Volume Velocity [m3 /s]0 5 10 15 20 25 30 35 40 45Time [ms]-0.200.20.40.60.81Normalized Volume Velocity [m3 /s]Figure 5.23: Glottal waveform (normalized volume velocity - ug) of a two-mass vocal fold model for vowel simulation /A/ with (left) uncoupledsystem (right) coupled system.10920 40 60 80 100 120 140 160 180 200 220510152025303540t=10ms20 40 60 80 100 120 140 160 180 200 220510152025303540t=30ms20 40 60 80 100 120 140 160 180 200 220510152025303540t=50msFigure 5.24: Mid-sagittal view of acoustic pressure wave propagation insidean uncoupled 2.5D vocal tract tube at three different time instants forvowel /a/ simulation. The two-mass vocal fold model is positionedat the glottal-end. The x and y axis of the above plots represent thenumber of grid cells.11020 40 60 80 100 120 140 160 180 200 220510152025303540t=10ms20 40 60 80 100 120 140 160 180 200 220510152025303540t=30ms20 40 60 80 100 120 140 160 180 200 220510152025303540t=50msFigure 5.25: Mid-sagittal view of acoustic pressure wave propagation insidea coupled 2.5D vocal tract tube at three different time instants forvowel /a/ simulation. The two-mass vocal fold model is positionedat the glottal-end. The x and y axis of the above plots represent thenumber of grid cells.1115.8 SummaryThis chapter characterized the acoustic features of the 2.5D FDTD vocal tractmodel for different vocal tract shapes in the quadrilateral vowel space and dis-cussed the importance of these results for the preliminary study of the 2.5D vocaltract model. Throughout the discussion, we also explored the current limitationsin our proposed vocal tract model. We began by introducing a number of con-texts to examine the following steps in our roadmap: SA-Analytical, SA-Voweland SA-Radiation.Our spectral analysis results (RS-Spectrum) showed that the 2.5D simulation-generated acoustic formants for a 1D one-end closed cylindrical tube model matchperfectly with the analytical solutions. The 2.5D model also outperformed the 1Dvocal tract model for cardinal vowels simulations in the planar wave regime. Wethen illustrated the impact of tube depth mapping by comparing the 2.5D and 2Dsimulation results to the precise 3D FEM model. This comparative study high-lighted the ability of the 2.5D model for precise estimation of acoustic wave prop-agation in a vocal tract tube without scaling the original area function nonlinearly.Thus, we concluded that the 2.5D model could produce synthesized acoustic out-put at a low spatial resolution comparable to a high-dimensional 3D FEM model.Consequently, the 2.5D model is computationally less expensive than the existing2D and 3D models. However, these results can not be generalized for all speechutterances due to the current limitations in our experimental study and the 2.5Dmodelling environment. We have discussed these limitations in detail in the nextchapter.In the end, we demonstrated an articulatory speech synthesizer toolkit by cou-pling the 2.5D model with the two-mass vocal fold model and produced synthe-sized vowel sounds (FT-Coupling). The next chapter summarizes the contributionof our work and suggests some of the future works for developing the 2.5D vocaltract model further.112Chapter 6ConclusionsIn this thesis, we introduced a new vocal tract acoustic model using the finite-difference time-domain technique. The proposed 2.5D vocal tract model has beendesigned in order to precisely simulate acoustic wave propagation by lumping non-planar wave interactions inside a 2D computational domain without compromisingthe time-complexity of the 2D wave solver. The model also maintains the abil-ity of the 2D wave solver to be accelerated using GPU architecture. The primaryobjective of our research was to balance the competing priorities between acous-tic precision and computational cost while directly using geometrical features ofan actual vocal tract. Since the vocal tract model works as a filter to attenuate thesource acoustic energy and produce synthesized audio outputs, our 2.5D vocal tractwill enable speech researchers to design real-time (or quasi-real-time) articulatoryspeech synthesizers.This chapter revisits and summarizes the contribution of our study. Here, wealso discuss some of the experimental and modelling limitations in our currentresearch work. We conclude this chapter with suggestions of potential applicationsand promising future directions, along with final remarks.6.1 Summary of ContributionsThis thesis has primarily focused on vocal tract acoustic modelling and improvingthe existing 2D models using the finite-difference numerical analysis method. In113Chapter 3, we first discussed current limitations in the 2D models by using Zappi’s2D vocal tract acoustic model implementation [146] as a reference. We concludedthat the existing 2D models do not approximate the wave propagation across themid-sagittal plane. Besides, they require a nonlinear scaling of the original vocaltract cross-sections and a high spatiotemporal resolution during the computationalsimulations to match the complex 3D vocal tract model’s acoustic spectral features[12]. In Chapter 4, we thus introduced a new vocal tract model (2.5D FDTD) byadding tube depth as an impedance parameter to the 2D acoustic wave solver. InChapter 5, we characterized the proposed model’s acoustic outputs using a subsetof the vocal tract shapes for the English vowel sounds and illustrated an articulatoryspeech synthesizer using the 2.5D model. A detailed summary of contributions inthis work is as follows.• Developed the 2.5D FDTD vocal tract acoustic model. We began by defin-ing the vocal tract tube depth, i.e., the tube’s height across the sagittal plane,as an additional geometrical parameter for the 2D acoustic wave solver. Wedemonstrated the procedure for measuring depths (D¯, D(x) and D(y)) fromthe vocal tract cross-sectional area functions and mapped them to their re-spective acoustic components (p, vx and vy) at each grid point inside a 2Dcomputational domain. We then derived the acoustic wave solver equationfor the 2.5D vocal tract model, that is, the 2D representation of Webster’shorn equation. The major advantage of our model is it allows capturing theeffect of off-plane wave interactions across the sagittal plane like a 3D vo-cal tract model. Therefore, unlike the 2D models, we do not require thearea function transformation and spatial grid resolution tuning with the 2.5Dvocal tract model.• Characterized acoustic features of the synthesized speech sounds gen-erated from the 2.5D vocal tract model. We illustrated the 2.5D vocaltract model’s basic capabilities by characterizing its acoustic properties foran analytical tube model; and open static vocal tract shapes (English vowelsounds), representing the extreme point of articulations in the quadrilateralvowel space. We built the lossless one-end closed cylindrical tubes using the2.5D FDTD scheme and examined the simulation-generated transfer func-114tions for the analytical tube model. The result showed the 2.5D scheme pro-duced an excellent match of acoustic formants in relation to their analyticsolutions.To study the acoustic characteristics of open static vocal tract shapeswith the 2.5D FDTD scheme, we used Story’s 1D area functions [109] toconstruct the vocal tract geometry. A comparative study of VTTFs was car-ried out between the 2.5D, 1D, 2D and 3D vocal tract models. In the planarwave regime (below 5 kHz), we noticed the acoustic formants of the 2.5Dand 1D vocal tract model were in good agreement, with a few exceptions.Overall, the result showed that the 2.5D model produced acoustic outputs(first 8 formants) comparable to the existing 2D FDTD and high-quality 3DFEM models even at a low spatial grid resolution.• Designed the 2.5D vocal tract radiation model and examined its acousticcharacteristics. We built the 2.5D vocal tract radiation model by construct-ing a circular baffle around the vocal tract as a simplified head geometryand imposed the free-field radiation condition at the mouth aperture to allowthe outward pressure wave propagation. Our vocal tract radiation model’sacoustic characteristics (formants’ positions) were consistent with the 3DFEM model, with few exceptions. However, they showed noticeable differ-ences compared to the formants retrieved from the recorded speech signalsdue to the oversimplified representation of the head geometry and mouthaperture in the 2.5D vocal tract radiation model.6.2 LimitationsIn this work, we made the simplifying assumptions of the vocal tract geometryinside the computational domain. We then characterized acoustic features of itssynthesized audio outputs only for a few preliminary steps from the roadmap (seeTable 5.1). Though these steps are crucial to verify the basic capabilities of the2.5D model, their acoustic outputs can not be generalized for all speech utterances.This section presents some of the current limitations to highlight the potential av-enues for the future development of the 2.5D model. We divide these details into115two categories: limitations of the 2.5D vocal tract model and those of the experi-mental study.6.2.1 2.5D Vocal Tract Modelling Limitations• Semi-realistic vocal tract geometry: The current 2.5D model assumes thevocal tract geometry as a straight tube with circular cross-sections. In theplanar wave regime, the influences of the vocal tract shapes having the samelength and cross-sectional areas are negligible on acoustic formants’ posi-tion. Thus, for a specific vowel sound, the acoustic outputs of our model willbe approximately similar to a realistic vocal tract geometry below 5 kHz.However, at a higher frequency, the appearance of higher-order modes in theVTTFs are eliminated due to the radial symmetry of the 2.5D vocal tracthaving circular cross-sections [6]. Similarly, a straight tube without any vo-cal tract side branches does not allow the appearance of antiresonances inthe VTTFs. Therefore, the current 2.5D model only produces resonances inits VTTFs, and its acoustic characteristics do not match the output of a realvocal tract geometry. As discussed in Section 4.5, it is to be noted that fullyirregular vocal tract cross-sections can not be represented using the 2.5Dmodel.• Boundary losses due to tube depths: The vocal tract wall impedance isone of the influential factors in determining the formant bandwidths in theVTTFs. As discussed previously, the current implementation of the 2.5Dmodel completely ignores the boundary losses due to the tube surface thatlies across the sagittal plane (tube depth). Thus, we have to manually tunethe formant bandwidths by changing the boundary admittance coefficient.However, a more realistic approach needs to be taken by including the walllosses for the tube depth. This loss mechanism could be incorporated byconsidering elastic vocal tract walls and allowing the tube vibration alongthe z axis as demonstrated for the 1D model [131].• Vocal tract mouth aperture: The acoustic characteristics of the vocal tractradiation model depend on the head geometry, vocal tract cross-sections andmouth aperture. The current 2.5D model approximates the geometry of the116mouth aperture using circular cross-sections only. Whereas, Arnela has al-ready demonstrated the positive impact of semi-realistic vocal tract cross-sections (elliptical) with lips on the acoustic formants of the 3D FEM radia-tion model [9]. Hence, we need to incorporate the elliptical mouth apertureand lips into the proposed 2.5D model to precisely model the vocal tractradiation mechanism.• Turbulent noises for fricative consonants: In this work, we have only in-vestigated the acoustic characteristics of the open static vocal tract shapesfor voiced sounds. For this case, articulators stay apart from each other at afixed position for a smooth flow of acoustic pressure waves. However, forsibilant sounds (fricative consonants: /s/, /z/), the articulators come close toproduce a wide range of turbulent flows. For such cases, the current linearacoustic wave equation may not address the aerodynamic factors as it gov-erns the mere propagation of sound waves. Instead, we need a flow equationthat can describe the sound generation by turbulent flows, such as the incom-pressible Navier-Stroke equation [91]. Therefore, the current 2.5D acousticwave solver needs to be augmented to generate the synthesized fricative con-sonants.• Closed vocal tract for stop consonants: The voiced and voiceless stopconsonants (e.g., /p/, /t/, /b/ and /d/) are produced by closing a segment ofthe airway for a short period so that no air can be passed through and thenreleasing that closure again [104]. A particular set of articulators (e.g., lipsand tongue blade) generally forms the closure. In such cases, the vocal tractcross-sectional area for the closed segments reduces to zero. As the vocaltract area function ultimately governs the current 2.5D model, the acousticsimulation of vocal tract shapes with closed constrictions may not producethe correct acoustic output. Therefore, our model needs to be investigatedfurther for the closed vocal tract shapes.1176.2.2 Experimental Study LimitationsIn this study, we have characterized acoustic outputs of the 2.5D vocal tract modelfor three English cardinal vowel sounds (/A/, /i/ , /u/) and /æ/. However, we havenot examined acoustic features of /a/ and /@/ vowel sounds due to the absence oftheir 1D cross-sectional area function datasets in literature. To utter /a/, the mouthaperture needs to be completely opened; simultaneously, the tongue should be po-sitioned at the extreme front, and for /@/, the tongue must be positioned halfwaybetween the front and back vowel. Vowel /a/ and /@/ are a few of the major casesin the vowel space that need to be verified, as they are the open-front unroundedand mid-central vowels, respectively. Hence, the current result can not be general-ized for all vowel sounds in the vowel quadrilateral space. Similarly, we are yet tocharacterize the 2.5D vocal tract model for vowel sounds in other languages thatforce the articulators to move toward extreme positions but do not come under thepre-defined vowel quadrilateral space of the English language. For example, theEnglish and mandarin vowels’ acoustic features are different due to the differentarticulatory positions [32].6.3 Future DirectionsIn Section 1.4, we outlined a complete roadmap with several milestones to developour vocal tract model. We used it throughout the previous chapters to discuss ourcurrent progress in terms of the model implementation. We are yet to finish all thesteps from our roadmap to present a fully developed 2.5D vocal tract model. Thus,the roadmap’s unfinished steps should be considered future directions for furtherdeveloping the 2.5D model. In this section, instead of going through each of thosesteps individually, we provide a general overview of future works concerning the2.5D vocal tract structural and acoustic modelling; and the possible future applica-tions.• Development of complex vocal tract geometries such as elliptical andmirror-plane symmetric cross-sections with a bent structure should be car-ried out to characterize acoustic features of the 2.5D model beyond the planewave cut-off frequency. Unlike circular cross-sections, the semi-irregulargeometry eliminates radial symmetry of a tube, and it allows the onset of118higher-order acoustic modes responsible for other vowel features such astimber (vowel naturalness) [6]. The inclusion of side branches affects VT-TFs by producing dips (antiresonances in both low and high frequencies) inthe simulation-generated transfer functions due to the interaction betweenthe right and left piriform fossa [117]. Similarly, incorporating a nasal tracthelps in synthesizing nasal vowel sounds. Therefore, these structural fea-tures need to be included in the current 2.5D vocal tract model to characterizeits acoustic outputs more precisely.• Designing dynamic vocal tract boundary walls still needs to be consideredfor the 2.5D vocal tract model to produce diphthong sounds (vowel-vowelutterances). The dynamic boundary condition requires the change in vocaltract cross-sections at runtime. Our current model adds a time-varying ge-ometry field β (x, t) to the acoustic wave solver equations (see Section 4.2.2),which can smoothly transition between fully solid (wall) and fully open (air)states to modify the vocal tract cross-sections during the simulation. How-ever, in this study, β has been remained fixed at each grid point to simu-late static vocal tract shapes. Another challenge in synthesizing diphthongsounds is the extraction of vocal tract cross-sectional areas during the kine-matic movements of articulators between the two static vocal tract shapesfor adjacent vowel sounds. The rt-MRI scan can produce multiple framesof still vocal tract images at different time samples. However, due to theirlow temporal resolution, they are not preferred for area function extraction.Alternatively, different interpolation methods can be used (e.g., semi-polargrid strategy and adaptive grid strategy) [11] to approximate the change vo-cal tract cross-sections. We can also use a machine learning algorithm forthis purpose [13, 118].• The nonlinear acoustic interactions between the vocal folds and the vo-cal tract are yet to be emulated using our proposed vocal tract model. In thisthesis, we demonstrated an articulatory speech synthesizer by coupling thetwo-mass vocal folds with the 2.5D vocal tract model to generate the synthe-sized vowel sounds. However, we have not examined the impact of couplingon the flow characteristics for both of these biomechanical models and the119output pressure waves at the listener-end for various speech utterances. Al-ternatively, to study the fluid-structure interactions between the vocal foldsand vocal tract model, the compressible Navier-Stokes fluid solvers can alsobe considered. However, this method requires solving more complex equa-tions in the computational domain.• The parallel implementation of the 2.5D FDTD acoustic wave solver us-ing modern GPUs needs to be done to synthesize speech utterances in real-time or quasi-real-time. As discussed earlier, Zappi et al. [146] has alreadydemonstrated the parallelization of the 2D FDTD acoustic wave solver. Asimilar approach can also be followed since the 2.5D model inherits the par-allel nature of the 2D FDTD wave solver algorithm. This implementationwill allow a faster acoustic simulation of vocal tract shapes than the currentone, which uses a serial machine to update the pressure and velocity compo-nents at each timestep in a sequential manner.• Unification of a biomechanical upper airway model to an acoustic modelallows a realistic control of vocal tract geometry through articulatory move-ments and producing the corresponding synthesized speech sounds. Dab-baghchian et al. have already demonstrated a biomechanical-acoustic modelby linking an adapted version of the Artisynth FRANK biomechanical modelwith the 3D FEM-based acoustic model [34]. A similar approach can also befollowed for the 2.5D acoustic model. This technique will help in the fasteracoustic simulation of synthesized speech outputs using the biomechanicalmodel. Moreover, it will create an opportunity for further exploration indesigning the neural-to-speech interfaces (Brain-Computer Interfaces-BCI).1206.4 Concluding RemarksWith the advancement of computing capabilities and high-resolution imaging tech-nology, the 3D articulatory models have come a long way for an accurate repre-sentation of upper airway geometry and acoustic simulation of various vocal tractshapes to synthesize speech sounds. However, such methods are still not feasiblefor clinical research and speech synthesis interfaces as they are computationallyexpensive.In this work, we presented a new vocal tract model for articulatory speechsynthesis. The model produces acoustic formants comparable to a high-quality 3Dvocal tract model while maintaining the time complexity of a 2D model for openstatic vocal tract shapes. An articulatory speech synthesis toolkit was presentedby coupling the model with a lumped-element vocal fold model for computationalsimulations of synthesized vowel sounds. This work provides a starting point forfurther developing an accurate and efficient 2.5D vocal tract model, which willeventually help design better articulatory speech synthesizers.121Bibliography[1] S. Adachi and J. Yu. Two-dimensional model of vocal fold vibration forsound synthesis of voice and soprano singing. The Journal of theAcoustical Society of America, 117(5):3213–3224, 2005. → page 17[2] F. Alipour, D. A. Berry, and I. R. Titze. A finite-element model ofvocal-fold vibration. The Journal of the Acoustical Society of America, 108(6):3003–3012, 2000. → page 20[3] A. Allen and N. Raghuvanshi. Aerophones in flatland: Interactive wavesimulation of wind instruments. ACM Transactions on Graphics (TOG), 34(4):1–11, 2015. → pages 32, 59, 108[4] P. Anderson, N. M. Harandi, S. Moisik, I. Stavness, and S. Fels. Acomprehensive 3d biomechanically-driven vocal tract model includinginverse dynamics for speech research. In Sixteenth Annual Conference ofthe International Speech Communication Association, 2015. → page 22[5] P. Anderson, S. Fels, I. Stavness, W. G. Pearson Jr, and B. Gick. Intravelarand extravelar portions of soft palate muscles in velic constrictions: Athree-dimensional modeling study. Journal of Speech, Language, andHearing Research, 62(4):802–814, 2019. → pages xv, 22[6] M. Arnela and O. Guasch. Finite element computation of elliptical vocaltract impedances using the two-microphone transfer function method. TheJournal of the Acoustical Society of America, 133(6):4197–4209, 2013. →pages 116, 119[7] M. Arnela and O. Guasch. Three-dimensional behavior in the numericalgeneration of vowels using tuned two-dimensional vocal tracts. InProceeding of 7th Forum Acousticum, 2014. → page 30[8] M. Arnela and O. Guasch. Two-dimensional vocal tracts withthree-dimensional behavior in the numerical generation of vowels. The122Journal of the Acoustical Society of America, 135(1):369–379, 2014. →pages 7, 30, 41, 44, 56, 65, 80[9] M. Arnela, O. Guasch, and F. Alías. Effects of head geometrysimplifications on acoustic radiation of vowel sounds based on time-domainfinite-element simulations. The Journal of the Acoustical Society ofAmerica, 134(4):2946–2954, 2013. → pages 74, 79, 86, 92, 103, 117[10] M. Arnela, S. Dabbaghchian, R. Blandin, O. Guasch, O. Engwall,A. Van Hirtum, and X. Pelorson. Influence of vocal tract geometrysimplifications on the numerical simulation of vowel sounds. The Journalof the Acoustical Society of America, 140(3):1707–1718, 2016. → pagesxv, 25, 27, 63, 64[11] M. Arnela, S. Dabbaghchian, O. Guasch, and O. Engwall. A semi-polargrid strategy for the three-dimensional finite element simulation ofvowel-vowel sequences. In 18th Annual Conference of the InternationalSpeech Communication Association, INTERSPEECH 2017, Stockholm,Sweden, 20 August 2017 through 24 August 2017, volume 2017, pages3477–3481. The International Speech Communication Association (ISCA),2017. → pages xiv, xv, 6, 25, 26, 27, 119[12] M. Arnela Coll et al. Numerical production of vowels and diphthongs usingfinite element methods. PhD thesis, Universitat Ramon Llull, 2015. →pages xv, 5, 31, 39, 49, 71, 86, 114[13] S. Asadiabadi and E. Erzin. Vocal tract contour tracking in rtmri using deeptemporal regression network. IEEE/ACM Transactions on Audio, Speech,and Language Processing, 28:3053–3064, 2020. → page 119[14] J. Awrejcewicz. Numerical analysis of the oscillations of human vocalcords. Nonlinear Dynamics, 2(1):35–52, 1991. → page 17[15] T. Baer, J. C. Gore, L. C. Gracco, and P. W. Nye. Analysis of vocal tractshape and dimensions using magnetic resonance imaging: Vowels. TheJournal of the Acoustical Society of America, 90(2):799–828, 1991. →page 26[16] O. Bälter, O. Engwall, A.-M. Öster, and H. Kjellström. Wizard-of-oz testof artur: a computer-based speech training system with articulationcorrection. In Proceedings of the 7th international ACM SIGACCESSconference on Computers and accessibility, pages 36–43, 2005. → page 4123[17] D. Beautemps, P. Badin, and R. Laboissière. Deriving vocal-tract areafunctions from midsagittal profiles and formant frequencies: A new modelfor vowels and fricative consonants based on experimental data. SpeechCommunication, 16(1):27–47, 1995. → page 28[18] S. Bilbao and J. Chick. Finite difference time domain simulation for thebrass instrument bore. The Journal of the Acoustical Society of America,134(5):3860–3871, 2013. → pages 32, 56[19] S. Bilbao and C. J. Webb. Physical modeling of timpani drums in 3d ongpgpus. Journal of the Audio Engineering Society, 61(10):737–748, 2013.→ page 32[20] P. Birkholz. A survey of self-oscillating lumped-element models of thevocal folds. Studientexte zur Sprachkommunikation: ElektronischeSprachsignalverarbeitung 2011, pages 47–58, 2011. → pages xiv, 17, 18[21] P. Birkholz. Enhanced area functions for noise source modeling in thevocal tract. In Proc. 10th Int. Seminar Speech Prod, pages 32–40, 2014. →page 28[22] P. Birkholz and B. J. Kröger. Vocal tract model adaptation using magneticresonance imaging. In Proceedings of the 7th International Seminar onSpeech Production, pages 493–500, 2006. → pages xv, 23, 24[23] P. Birkholz, D. Jackèl, and B. J. Kroger. Construction and control of athree-dimensional vocal tract model. In 2006 IEEE InternationalConference on Acoustics Speech and Signal Processing Proceedings,volume 1, pages I–I. IEEE, 2006. → page 24[24] P. Birkholz, D. Jackèl, and B. J. Kroger. Simulation of losses due toturbulence in the time-varying vocal system. IEEE Transactions on Audio,Speech, and Language Processing, 15(4):1218–1226, 2007. → page 32[25] S. Boij and B. Nilsson. Reflection of sound at area expansions in a flowduct. Journal of Sound and Vibration, 260(3):477–498, 2003. → page 44[26] R. Bossart, N. Joly, and M. Bruneau. Hybrid numerical and analyticalsolutions for acoustic boundary problems in thermo-viscous fluids. Journalof Sound and Vibration, 263(1):69–84, 2003. → page 63[27] J. Botts and L. Savioja. Integrating finite difference schemes for scalar andvector wave equations. In 2013 IEEE International conference on124acoustics, speech and signal processing, pages 171–175. IEEE, 2013. →page 32[28] A. S. Brody, D. P. Frush, W. Huda, R. L. Brent, et al. Radiation risk tochildren from computed tomography. Pediatrics, 120(3):677–682, 2007.→ page 26[29] S. Buchaillard, P. Perrier, and Y. Payan. A biomechanical model of cardinalvowel production: Muscle activations and the impact of gravity on tonguepositioning. The Journal of the Acoustical Society of America, 126(4):2033–2051, 2009. → page 22[30] J. E. Cahn. Generating expression in synthesized speech. PhD thesis,Massachusetts Institute of Technology, 1989. → page 2[31] H. Chen, S. Fels, T. Pang, L. Tsou, F. R. de Almeida, and A. A. Lowe.Three-dimensional reconstruction of soft palate modeling fromsubject-specific magnetic resonance imaging data. Sleep and Breathing, 16(4):1113–1119, 2012. → page 22[32] H. C. Chen and M. J. Wang. An acoustic analysis of chinese and englishvowels. 2011. → page 118[33] S. Dabbaghchian, M. Arnela, and O. Engwall. Simplification of vocal tractshapes with different levels of detail. In 18th International Congress ofPhonetic Science, Glasgow, UK, 10-14 Aug, 2015, pages 1–5. University ofGlasgow, 2015. → page 27[34] S. Dabbaghchian, M. Arnela, O. Engwall, and O. Guasch. Simulation ofvowel-vowel utterances using a 3d biomechanical-acoustic model.International Journal for Numerical Methods in Biomedical Engineering,37(1):e3407, 2021. → page 120[35] J. Dang and K. Honda. Acoustic characteristics of the piriform fossa inmodels and humans. The Journal of the Acoustical Society of America, 101(1):456–465, 1997. → page 75[36] J. Dang and K. Honda. Construction and control of a physiologicalarticulatory model. The Journal of the Acoustical Society of America, 115(2):853–870, 2004. → page 21[37] M. de Oliveira Rosa, J. C. Pereira, M. Grellet, and A. Alwan. Acontribution to simulating a three-dimensional larynx model using the finite125element method. The Journal of the Acoustical Society of America, 114(5):2893–2905, 2003. → page 20[38] M. De Vries, H. Schutte, and G. Verkerke. Determination of parameters forlumped parameter models of the vocal folds using a finite-element methodapproach. The Journal of the Acoustical Society of America, 106(6):3620–3628, 1999. → page 20[39] B. Delvaux and D. Howard. A new method to explore the spectral impactof the piriform fossae on the singing voice: benchmarking using mri-based3d-printed vocal tracts. PLoS One, 9(7):e102680, 2014. → page 75[40] D. Demolin, T. Metens, and A. Soquet. Three-dimensional measurement ofthe vocal tract by mri. In Proceeding of Fourth International Conferenceon Spoken Language Processing. ICSLP’96, volume 1, pages 272–275.IEEE, 1996. → page 3[41] O. Engwall. Combining mri, ema and epg measurements in athree-dimensional tongue model. Speech Communication, 41(2-3):303–329, 2003. → page 24[42] O. Engwall, O. Bälter, A.-M. Öster, and H. Kjellström. Designing the userinterface of the computer-based speech training system artur based on earlyuser tests. Behaviour & Information Technology, 25(4):353–365, 2006. →page 4[43] G. Fant. Acoustic theory of speech production. Number 2. Walter deGruyter, 1970. → pages 12, 24[44] G. Fant, J. Liljencrants, and Q.-g. Lin. A four-parameter model of glottalflow. STL-QPSR, 4(1985):1–13, 1985. → page 16[45] J. Flanagan and L. Landgraf. Self-oscillating source for vocal-tractsynthesizers. IEEE Transactions on Audio and Electroacoustics, 16(1):57–64, 1968. → page 17[46] J. L. Flanagan. Speech analysis synthesis and perception, volume 3.Springer Science & Business Media, 2013. → pages 70, 86[47] M. Fleischer, S. Pinkert, W. Mattheus, A. Mainka, and D. Mürbe. Formantfrequencies and bandwidths of the vocal tract transfer function are affectedby the mechanical impedance of the vocal tract wall. Biomechanics andmodeling in mechanobiology, 14(4):719–733, 2015. → pages 30, 63, 103126[48] N. H. Fletcher and T. D. Rossing. The physics of musical instruments.Springer Science & Business Media, 2012. → pages 45, 65[49] J.-M. Gérard, R. Wilhelms-Tricarico, P. Perrier, and Y. Payan. A 3ddynamical biomechanical tongue model to study speech motor control.arXiv preprint physics/0606148, 2006. → page 21[50] B. Gick, I. Wilson, and D. Derrick. Articulatory phonetics. John Wiley &Sons, 2012. → page 12[51] U. G. Goldstein. An articulatory model for the vocal tracts of growingchildren. PhD thesis, Massachusetts Institute of Technology, 1980. → page28[52] N. Hanna, J. Smith, and J. Wolfe. Frequencies, bandwidths and magnitudesof vocal tract and surrounding tissue resonances, measured through the lipsduring phonation. The Journal of the Acoustical Society of America, 139(5):2924–2936, 2016. → page 103[53] A. Hannukainen, T. Lukkari, J. Malinen, and P. Palo. Vowel formants fromthe wave equation. The Journal of the Acoustical Society of America, 122(1):EL1–EL7, 2007. → page 30[54] K. Hashimoto and S. Suga. Estimation of the muscular tensions of thehuman tongue by using a three-dimensional model of the tongue. Journalof the Acoustical Society of Japan (E), 7(1):39–46, 1986. → page 21[55] J. M. Heinz and K. N. Stevens. On the derivation of area functions andacoustic spectra from cineradiographic films of speech. The Journal of theAcoustical Society of America, 36(5):1037–1038, 1964. → page 28[56] H. Hirai. A physiological model of speech organs incorporatingtongue-laynx interaction. J Acoust. Soc. Jpn.(J), 52:918, 1995. → page 21[57] J. C. Ho, M. Zañartu, and G. R. Wodicka. An anatomically based,time-domain acoustic model of the subglottal system for speechproduction. The Journal of the Acoustical Society of America, 129(3):1531–1547, 2011. → page 5[58] E. A. Hoffman, D. Gnanaprakasam, K. B. Gupta, J. D. Hoford, S. D.Kugelmass, and R. S. Kulawiec. Vida: An environment formultidimensional image display and analysis. In Biomedical ImageProcessing and Three-Dimensional Microscopy, volume 1660, pages694–711. International Society for Optics and Photonics, 1992. → page 27127[59] K. Honda, H. Takemoto, T. Kitamura, S. Fujita, and S. Takano. Exploringhuman speech production mechanisms by mri. IEICE TRANSACTIONS onInformation and Systems, 87(5):1050–1058, 2004. → page 75[60] A. Iida, N. Campbell, S. Iga, F. Higuchi, and M. Yasumura. A speechsynthesis system with emotion for assisting communication. In ISCAtutorial and research workshop (ITRW) on speech and emotion, 2000. →page 2[61] K. Ishizaka. Fluid mechanical considerations of vocal cord vibration.SCRL Monogr., 8:28–72, 1972. → pages xx, 4, 15, 105, 106[62] K. Ishizaka and J. L. Flanagan. Synthesis of voiced sounds from atwo-mass model of the vocal cords. Bell system technical journal, 51(6):1233–1268, 1972. → pages 17, 105[63] W. Jiang, X. Zheng, and Q. Xue. Computational modeling offluid–structure–acoustics interaction during voice production. Frontiers inbioengineering and biotechnology, 5:7, 2017. → page 20[64] W. Jiang, J. H. Rasmussen, Q. Xue, M. Ding, X. Zheng, and C. P. Elemans.High-fidelity continuum modeling predicts avian voiced sound production.Proceedings of the National Academy of Sciences, 117(9):4718–4723,2020. → page 20[65] Y. Kagawa, R. Shimoyama, T. Yamabuchi, T. Murai, and K. Takarada.Boundary element models of the vocal tract and radiation field and theirresponse characteristics. Journal of sound and vibration, 157(3):385–403,1992. → page 30[66] Y. Kakita and O. Fujimura. Computation of mapping from muscularcontraction patterns to formant patterns in vowel space. The Journal of theAcoustical Society of America, 74(S1):S117–S117, 1983. → page 21[67] T. Kako and K. Touda. Numerical method for voice generation problembased on finite element method. Journal of Computational Acoustics, 14(01):45–56, 2006. → page 30[68] J. Kelly. Speech synthesis. Proc. 4th Int. Congr. Acoustics, 1962, pages1–4, 1962. → page 31[69] T. Kitamura, K. Honda, and H. Takemoto. Individual variation of thehypopharyngeal cavities and its acoustic effects. Acoustical science andtechnology, 26(1):16–26, 2005. → page 75128[70] T. Koizumi, S. Taniguchi, and S. Hiromitsu. Two-mass models of the vocalcords for natural sounding voice synthesis. The Journal of the AcousticalSociety of America, 82(4):1179–1192, 1987. → page 17[71] B. J. Kröger and P. Birkholz. Articulatory synthesis of speech and singing:State of the art and suggestions for future research. In Multimodal Signals:Cognitive and Algorithmic Issues, pages 306–319. Springer, 2009. → page32[72] Y. Laprie, M. Loosvelt, S. Maeda, R. Sock, and F. Hirsch. Articulatorycopy synthesis from cine x-ray films. 2013. → page 29[73] P. Lieberman and S. E. Blumstein. Speech physiology, speech perception,and acoustic phonetics. Cambridge University Press, 1988. → page 12[74] Y. Liu. Finite-difference numerical modeling of any even order accuracy.Oil Geophysical Prospecting, 33:1–10, 1998. → page 66[75] Y. Liu. Maximizing the cfl number of stable time–space domain explicitfinite-difference modeling. Journal of Computational Physics, 416:109501,2020. → page 66[76] J. E. Lloyd, I. Stavness, and S. Fels. Artisynth: A fast interactivebiomechanical modeling toolkit combining multibody and finite elementsimulation. In Soft tissue biomechanical modeling for computer assistedsurgery, pages 355–394. Springer, 2012. → page 22[77] J. C. Lucero and K. G. Munhall. A model of facial biomechanics forspeech production. The Journal of the Acoustical Society of America, 106(5):2834–2842, 1999. → page 21[78] H. Matsuzaki, N. Miki, and Y. Ogawa. 3d finite element analysis ofjapanese vowels in elliptic sound tube model. Electronics andCommunications in Japan (Part III: Fundamental Electronic Science), 83(4):43–51, 2000. → page 30[79] P. Mermelstein. Articulatory model for the study of speech production. TheJournal of the Acoustical Society of America, 53(4):1070–1082, 1973. →page 24[80] R. Mittal, X. Zheng, R. Bhardwaj, J. H. Seo, Q. Xue, and S. Bielamowicz.Toward a simulation-based tool for the treatment of vocal fold paralysis.Frontiers in physiology, 2:19, 2011. → page 20129[81] D. R. Mohapatra and S. Fels. Limitations of source-filter coupling inphonation. arXiv preprint arXiv:1811.07435, 2018. → pages 17, 77[82] P. Mokhtari, H. Takemoto, and T. Kitamura. Single-matrix formulation of atime domain acoustic model of the vocal tract with side branches. SpeechCommunication, 50(3):179–190, 2008. → page 5[83] C. A. Moore. The correspondence of vocal tract resonance with volumesobtained from magnetic resonance images. Journal of Speech, Language,and Hearing Research, 35(5):1009–1023, 1992. → page 26[84] K. Motoki. Three-dimensional acoustic field in vocal-tract. AcousticalScience and Technology, 23(4):207–212, 2002. → page 30[85] J. Mullen, D. M. Howard, and D. T. Murphy. Waveguide physical modelingof vocal tract acoustics: flexible formant bandwidth control from increasedmodel dimensionality. IEEE Transactions on Audio, Speech, and LanguageProcessing, 14(3):964–971, 2006. → page 40[86] J. Müller and W. Baly. The physiology of the senses, voice, and muscularmotion, with the mental faculties. Taylor, Walton & Maberly, 1848. →page 12[87] I. R. Murray and J. L. Arnott. Toward the simulation of emotion insynthetic speech: A review of the literature on human vocal emotion. TheJournal of the Acoustical Society of America, 93(2):1097–1108, 1993. →page 2[88] OpenStax. Anatomy & physiology.http://cnx.org/contents/14fb4ad7-39a1-4eee-ab6e-3ef2482e3e22@8.24.,Feb 26, 2016. OpenStax CNX. → pages xiv, 13, 14[89] J. S. Perkell. A physiologically-oriented model of tongue activity in speechproduction. 1974. → page 21[90] P. Perrier, L.-J. Boë, and R. Sock. Vocal tract area function estimation frommidsagittal dimensions with ct scans and a vocal tract cast: modeling thetransition with two sets of coefficients. Journal of Speech, Language, andHearing Research, 35(1):53–67, 1992. → page 28[91] A. Pont, O. Guasch, and M. Arnela. Finite element generation ofsibilants/s/and/z/using random distributions of kirchhoff vortices.International journal for numerical methods in biomedical engineering, 36(2):e3302, 2020. → page 117130[92] E. Rank and H. Pirker. Generating emotional speech with a concatenativesynthesizer. In Fifth International Conference on Spoken LanguageProcessing, 1998. → page 2[93] C. A. Rosen. Phonosurgical vocal fold injection: procedures and materials.Otolaryngologic Clinics of North America, 33(5):1087–1096, 2000. →page 3[94] T. D. Rossing and N. Fletcher. The physics of musical instruments.Springer-Verlag, 1991. → pages 45, 58[95] M. Schröder. Emotional speech synthesis: A review. In Seventh EuropeanConference on Speech Communication and Technology, 2001. → page 2[96] D. Sciamarella and C. d’Alessandro. A study of the two-mass model interms of acoustic parameters. In Seventh International Conference onSpoken Language Processing, 2002. → page 17[97] P. Šidlof and J. HORÁ CEK. Vocal fold motion and voice production:mathematical modelling and experiment. In Proceedings of ForumAcusticum, pages 2737–2742. → page 17[98] E. Sifakis, A. Selle, A. Robinson-Mosher, and R. Fedkiw. Simulatingspeech with a physics-based facial muscle model. In Proceedings of the2006 ACM SIGGRAPH/Eurographics symposium on Computer animation,pages 261–270. Citeseer, 2006. → page 21[99] J. O. Smith. Physical modeling using digital waveguides. Computer musicjournal, 16(4):74–91, 1992. → page 2[100] S. W. Smith et al. The scientist and engineer’s guide to digital signalprocessing, volume 14. California Technical Pub. San Diego, 1997. →page 77[101] M. Sondhi and J. Schroeter. A hybrid time-frequency domain articulatoryspeech synthesizer. IEEE Transactions on Acoustics, Speech, and SignalProcessing, 35(7):955–967, 1987. → pages 2, 31, 105[102] M. Speed, D. Murphy, and D. M. Howard. Characteristics oftwo-dimensional finite difference techniques for vocal tract analysis andvoice synthesis. In Tenth Annual Conference of the International SpeechCommunication Association, 2009. → page 41131[103] I. Stavness, J. E. Lloyd, Y. Payan, and S. Fels. Coupled hard–soft tissuesimulation with contact and constraints applied to jaw–tongue–hyoiddynamics. International Journal for Numerical Methods in BiomedicalEngineering, 27(3):367–390, 2011. → page 22[104] K. N. Stevens. Models for the production and acoustics of stop consonants.Speech communication, 13(3-4):367–375, 1993. → page 117[105] K. N. Stevens and A. S. House. Development of a quantitative descriptionof vowel articulation. The Journal of the Acoustical Society of America, 27(3):484–493, 1955. → page 12[106] S. Stone, M. Marxen, and P. Birkholz. Construction and evaluation of aparametric one-dimensional vocal tract model. IEEE/ACM Transactions onAudio, Speech, and Language Processing, 26(8):1381–1392, 2018. →pages xiv, xv, 6, 21[107] B. H. Story. An overview of the physiology, physics and modeling of thesound source for vowels. Acoustical Science and Technology, 23(4):195–206, 2002. → page 17[108] B. H. Story. Synergistic modes of vocal tract articulation for americanenglish vowels. The Journal of the Acoustical Society of America, 118(6):3834–3859, 2005. → page 39[109] B. H. Story. Comparison of magnetic resonance imaging-based vocal tractarea functions obtained from the same speaker in 1994 and 2002. TheJournal of the Acoustical Society of America, 123(1):327–335, 2008. →pages 39, 74, 92, 93, 115, 141[110] B. H. Story and I. R. Titze. Voice simulation with a body-cover model ofthe vocal folds. The Journal of the Acoustical Society of America, 97(2):1249–1260, 1995. → page 19[111] B. H. Story, I. R. Titze, and E. A. Hoffman. Vocal tract area functions frommagnetic resonance imaging. The Journal of the Acoustical Society ofAmerica, 100(1):537–554, 1996. → pages 26, 27[112] A. M. Sulter, D. G. Miller, R. F. Wolf, H. K. Schutte, H. P. Wit, and E. L.Mooyaart. On the relation between the dimensions and resonancecharacteristics of the vocal tract: a study with mri. Magnetic ResonanceImaging, 10(3):365–373, 1992. → page 26132[113] J. Sundberg, C. Johansson, H. Wilbrand, and C. Ytterbergh. From sagittaldistance to area. Phonetica, 44(2):76–90, 1987. → page 28[114] P. Švancara and J. Horácˇek. Numerical modelling of effect of tonsillectomyon production of czech vowels. Acta Acustica united with Acustica, 92(5):681–688, 2006. → pages 30, 92[115] P. Švancara, J. Horácek, J. Vokrál, and L. Cerny`. Computational modellingof effect of tonsillectomy on voice production. Logopedics PhoniatricsVocology, 31(3):117–125, 2006. → pages 30, 76[116] H. Takemoto, P. Mokhtari, and T. Kitamura. Acoustic analysis of the vocaltract during vowel production by finite-difference time-domain method. theJournal of the acoustical society of America, 128(6):3724–3738, 2010. →pages 63, 77[117] H. Takemoto, S. Adachi, P. Mokhtari, and T. Kitamura. Acousticinteraction between the right and left piriform fossae in generating spectraldips. The Journal of the Acoustical Society of America, 134(4):2955–2964,2013. → page 119[118] H. Takemoto, T. Goto, Y. Hagihara, S. Hamanaka, T. Kitamura, Y. Nota,and K. Maekawa. Speech organ contour extraction using real-time mri andmachine learning method. In Interspeech, pages 904–908, 2019. → page119[119] K. Tang. A feasibility study of template-based subject-specific modellingand simulation of upper-airway complex. PhD thesis, University of BritishColumbia, 2017. → pages xi, 23[120] I. R. Titze. The physics of small-amplitude oscillation of the vocal folds.The Journal of the Acoustical Society of America, 83(4):1536–1552, 1988.→ page 15[121] I. R. Titze. A four-parameter model of the glottis and vocal fold contactarea. Speech communication, 8(3):191–201, 1989. → page 16[122] I. R. Titze. Mechanical stress in phonation. Journal of Voice, 8(2):99–105,1994. → page 3[123] I. R. Titze. Nonlinear source–filter coupling in phonation: Theory. TheJournal of the Acoustical Society of America, 123(4):1902–1915, 2008. →page 108133[124] I. T. Tokuda, J. Horácˇek, J. G. Švec, and H. Herzel. Comparison ofbiomechanical modeling of register transitions and voice instabilities withexcised larynx experiments. The Journal of the Acoustical Society ofAmerica, 122(1):519–531, 2007. → page 17[125] G. Ungeheuer. Elemente einer akustischen Theorie der Vokalartikulation.Springer-Verlag, 2013. → page 59[126] T. Vampola, J. Horácˇek, and J. G. Švec. Fe modeling of human vocal tractacoustics. part i: Production of czech vowels. Acta acustica united withAcustica, 94(3):433–447, 2008. → pages 30, 92[127] T. Vampola, J. Horácˇek, J. Vokárˇál, and L. Cˇerny`. Fe modeling of humanvocal tract acoustics. part ii: Influence of velopharyngeal insufficiency onphonation of vowels. Acta Acustica United with Acustica, 94(3):448–460,2008. → page 92[128] T. Vampola, A.-M. Laukkanen, J. Horácˇek, and J. G. Švec. Vocal tractchanges caused by phonation into a tube: a case study using computertomography and finite-element modeling. The Journal of the AcousticalSociety of America, 129(1):310–315, 2011. → page 30[129] T. Vampola, J. Horácˇek, A.-M. Laukkanen, and J. G. Švec. Human vocaltract resonances and the corresponding mode shapes investigated bythree-dimensional finite-element modelling based on ct measurement.Logopedics Phoniatrics Vocology, 40(1):14–23, 2015. → page 30[130] J. Van den Berg. Myoelastic-aerodynamic theory of voice production.Journal of speech and hearing research, 1(3):227–244, 1958. → page 15[131] K. van den Doel and U. M. Ascher. Real-time numerical solution ofwebster’s equation on a nonuniform grid. IEEE transactions on audio,speech, and language processing, 16(6):1163–1172, 2008. → pages5, 71, 75, 85, 105, 116[132] A. Vasudevan. The voice box: a fast coupled vocal fold model forarticulatory speech synthesis. PhD thesis, University of British Columbia,2017. → pages xv, 4, 19, 108[133] R. Veldhuis. A computationally efficient alternative for theliljencrants–fant model and its perceptual evaluation. The Journal of theAcoustical Society of America, 103(1):566–571, 1998. → page 76134[134] J. Vroomen, R. Collier, and S. Mozziconacci. Duration and intonation inemotional speech. In Third European Conference on SpeechCommunication and Technology, 1993. → page 2[135] H. Wakita and G. Fant. Toward a better vocal tract model. STL-QPSR, 19(1):9–29, 1978. → page 70[136] X. Wang, W. Jiang, X. Zheng, and Q. Xue. A computational study of theeffects of vocal fold stiffness parameters on voice production. Journal ofVoice, 2019. → page 20[137] Y. Wang, H. Wang, J. Wei, and J. Dang. Mandarin vowel synthesis basedon 2d and 3d vocal tract model by finite-difference time-domain method.In Proceedings of The 2012 Asia Pacific Signal and Information ProcessingAssociation Annual Summit and Conference, pages 1–4. IEEE, 2012. →pages 41, 103[138] C. J. Webb. Parallel computation techniques for virtual acoustics andphysical modelling synthesis. 2014. → page 32[139] C. J. Webb and S. Bilbao. Computing room acoustics with cuda-3d fdtdschemes with boundary losses and viscosity. In 2011 IEEE InternationalConference on Acoustics, Speech and Signal Processing (ICASSP), pages317–320. IEEE, 2011. → page 32[140] Q. Xue, R. Mittal, X. Zheng, and S. Bielamowicz. Computationalmodeling of phonatory dynamics in a tubular three-dimensional model ofthe human larynx. The Journal of the Acoustical Society of America, 132(3):1602–1613, 2012. → page 20[141] K. Yee. Numerical solution of initial boundary value problems involvingmaxwell’s equations in isotropic media. IEEE Transactions on antennasand propagation, 14(3):302–307, 1966. → pages 33, 59[142] K. S. Yee and J. S. Chen. The finite-difference time-domain (fdtd) and thefinite-volume time-domain (fvtd) methods in solving maxwell’s equations.IEEE Transactions on Antennas and Propagation, 45(3):354–363, 1997.→ page 33[143] T. Yokota, S. Sakamoto, and H. Tachibana. Visualization of soundpropagation and scattering in rooms. Acoustical science and technology, 23(1):40–46, 2002. → page 32135[144] X. Yuan, D. Borup, J. W. Wiskin, M. Berggren, R. Eidens, and S. A.Johnson. Formulation and validation of berenger’s pml absorbing boundaryfor the fdtd simulation of acoustic scattering. IEEE Transactions onUltrasonics, Ferroelectrics, and Frequency Control, 44(4):816–822, 1997.→ page 32[145] M. Zanartu, L. Mongeau, and G. R. Wodicka. Influence of acoustic loadingon an effective single mass model of the vocal folds. The Journal of theAcoustical Society of America, 121(2):1119–1129, 2007. → page 17[146] V. Zappi, A. Vasuvedan, A. Allen, N. Raghuvanshi, and S. Fels. Towardsreal-time two-dimensional wave propagation for articulatory speechsynthesis. In Proceedings of Meetings on Acoustics 171ASA, volume 26,page 045005. Acoustical Society of America, 2016. → pagesxii, 7, 40, 41, 54, 56, 59, 71, 94, 105, 114, 120[147] C. Zhang, W. Zhao, S. H. Frankel, and L. Mongeau. Computationalaeroacoustics of phonation, part ii: Effects of flow parameters andventricular folds. The Journal of the Acoustical Society of America, 112(5):2147–2154, 2002. → page 30[148] Z. Zhang. Characteristics of phonation onset in a two-layer vocal foldmodel. The Journal of the Acoustical Society of America, 125(2):1091–1102, 2009. → page 19[149] Z. Zhang. Mechanics of human voice production and control. The journalof the acoustical society of america, 140(4):2614–2635, 2016. → page 15[150] X. Zheng, Q. Xue, R. Mittal, and S. Beilamowicz. A coupledsharp-interface immersed boundary-finite-element method forflow-structure interaction with application to human phonation. Journal ofbiomechanical engineering, 132(11), 2010. → page 20[151] X. Zhou, C. Y. Espy-Wilson, S. Boyce, M. Tiede, C. Holland, and A. Choe.A magnetic resonance imaging-based articulatory and acoustic study of“retroflex” and “bunched” american english/r. The Journal of theAcoustical Society of America, 123(6):4466–4481, 2008. → page 92136Appendix AGeneral MathematicsA.1 1D Acoustic Wave EquationArea (S)ABCDA’B’C’D’ξdx𝝃 +𝝏𝝃𝝏𝒙𝒅𝒙𝟏 +𝝏𝝃𝝏𝒙𝒅𝒙wave propagationFigure A.1: As acoustic wave displaced a distance of ξ , the air particles onplane AB move to A’B’ and on plane CD move to C’D’Considering acoustic plane waves propagate in a tube (as shown in Figure A.1)along the axial direction (perpendicular to the tube width), ξ represents the dis-placement of air due to the acoustic wave along the x axis. Due to the air displace-ment, the element ABCD of thickness dx moves to A’B’C’D’. The displacement ofair particles at CD can be determined using the Taylor series. Assuming the cross-137sectional area of the tube as S, we can define the initial volume of the elementas V = Sdx. As the air displaces, the resulting change in volume of the elementbecomesdV = sdx(∂ξ∂x)(A.1)And, the total volume of this element can be represented asV +dV = sdx(1+∂ξ∂x)(A.2)Now considering the change in air pressure as d pa due to the displacement of theelement, then the bulk modulus (K) can be defined by the following relation,d pa =−K dVV (A.3)We can describe the changing air pressure d pa simply as acoustic pressure of p.Comparing Equation A.1 and A.3 givesp =−K ∂ξ∂x(A.4)Finally, the motion of the element ABCD can be expressed using Newton’s equa-tion. This can done by setting the pressure-gradient force along the x direction asmass times acceleration as follows−S(∂ p∂xdx)= ρSdx∂ 2ξ∂ t2(A.5)or−∂ p∂x= ρ∂ 2ξ∂ t2(A.6)where ρ is the density of the element. Differentiating Equation A.4 with respect to138x and equating it with Equation A.6 yields∂ 2ξ∂ t2=Kρ∂ 2ξ∂x2(A.7)where Kρ = c2 and c is the speed of sound. Equation A.7 can also be stated usingthe acoustic pressure as follows∂ 2 p∂ t2= c2∂ 2 p∂x2(A.8)Equations A.7 and A.8 are two different versions of the 1D acoustic wave equationswhile referring to the acoustic displacement ξ and acoustic pressure p, respectively.A.2 Approximation of first-order derivativesThis section approximates the derivative of a function at a point x0, which fallsbetween two sampling points. The sampling points are uniformly spaced, havingspacial distance δ . The function can be defined at the sampling points (i.e., x0 + δ2and x0− δ2 ) using the Taylor series expansion as follows,f (x0 +δ2) = f (x0)+δ2f ′(x0)+12!(δ2)2f ′′(x0)+13!(δ2)3f ′′′(x0)+ ... (A.9)f (x0− δ2 ) = f (x0)−δ2f ′(x0)+12!(δ2)2f ′′(x0)− 13!(δ2)3f ′′′(x0)+ ... (A.10)Subtracting (A.9) from (A.10) yields,f (x0 +δ2)− f (x0− δ2 ) = δ f′(x0)+23!(δ2)3f ′′′(x0)+ ... (A.11)Rearranging (A.11) yields,139f ′(x0) =(d f (x)dx)x=x0=f (x0 + δ2 )− f (x0− δ2 )δ+O(δ 2) (A.12)The "Big-Oh" term represents all the terms which has not been shown here explic-itly. As the spatial width (δ ) between two sampling points are sufficiently small, itshigher power terms (≥ 2) can be neglected for a reasonable approximation. Hencethe second-order accurate differential equation can be represented as,f ′(x0) =f (x0 + δ2 )− f (x0− δ2 )δ(central difference) (A.13)Equation A.13 can also be written using the forward difference as follows,f ′(x0) =f (x0 +δ )− f (x0)δ(forward difference) (A.14)140Appendix B2D Vocal Tract Structural ModelHere we illustrate the construction of the mid-sagittal vocal tract geometry us-ing the 2D Yee-Grid scheme for vowel sounds. The geometry uses Story’s cross-sectional area function [109] to create straight tubes with circular cross-sections.Inside the grid, the diameter of each tube segment is approximated by odd numberof cells (shown in Figure B.1). The CFL condition is followed to determine thegrid spatial resolution.Grid Edges Vocal Tract ContourFigure B.1: An illustration of vocal tract contour in a 2D rectangular grid.Number of cells between the upper and lower contour is always odd.141vowel [A]vowel [i]142vowel [I]vowel [ U]143vowel [u]vowel [o]Figure B.2: Mid-sagittal cross-sectional area function against vocal tractlength and the corresponding vocal tract structural models for vowels:[A, i, I, U, u, o]144Appendix CAcoustic Analysis of 2.5D FDTDModelHere we present the acoustic formants’ positions for various simulations.C.1 1D one-end closed lossless cylindrical tubeFormants AnalyticalsolutionNumericalsolution(A1 = 0.91cm2)Numericalsolution(A2 = 1.31cm2)F1 512 520 (1.56%) 520 (1.56%)F2 1537 1560 (1.49%) 1540 (0.19%)F3 2562 2580 (0.70%) 2580 (0.70%)F4 3588 3620 (0.89%) 3620 (0.89%)F5 4613 4660 (1.01%) 4640 (0.58%)F6 5638 5780 (0.74%) 5680 (0.74%)F7 6663 6720 (0.85%) 6720 (0.85%)F8 7688 7760 (0.93%) 7740 (0.67%)Table C.1: Formant locations Fi (Hz) and formant deviations ∆Fi/Fi (%) ofthe 2.5Dlow simulations of one-end closed cylindrical tube models (L1 =17.07cm) with respect to their analytical solution.145Formants AnalyticalsolutionNumericalsolution(A1 = 0.91cm2)Numericalsolution(A2 = 1.31cm2)F1 175 180 (2.85%) 180 (2.85%)F2 523 520 (-0.57%) 520 (-0.57%)F3 872 880 (0.91%) 880 (0.91%)F4 1221 1240 (1.55%) 1240 (1.55%)F5 1569 1580 (0.70%) 1580 (0.70%)F6 1918 1940 (1.14%) 1940 (1.14%)F7 2267 2300 (1.45%) 2280 (0.57%)F8 2616 2640 (0.91%) 2640 (0.91%)Table C.2: Formant locations Fi (Hz) and formant deviations ∆Fi/Fi (%) ofthe 2.5Dlow simulations of one-end closed cylindrical tube models (L2 =50.16cm) with respect to their analytical solution.Formants Analyticalsolution1D Numericalmodel2.5D NumericalmodelF1 512 553 (8.00%) 520 (1.44%)F2 1537 1504 (-2.19%) 1560 (1.44%)F3 2562 2474 (-3.43%) 2580 (0.66%)Table C.3: Formant locations Fi (Hz) and formant deviations ∆Fi/Fi (%) ofthe 1D and 2.5Dlow simulation of a one-end closed cylindrical tubemodel (L1 = 17.07cm, A1 = 0.91cm2) with respect to its analytical solu-tion.146C.2 2.5D simulation of vowel sounds with different gridresolutionsFormants 3D 2.5Dlow 2.5Dmid 2.5DhighF1 696 700 (0.57%) 680 (-2.29%) 680 (-2.29%)F2 1068 1040 (-2.62%) 1060 (-0.74%) 1040 (-2.62%)F3 3031 3020 (-0.36%) 3060 (0.95%) 3040 (0.29%)F4 4064 4120 (1.37%) 4080 (0.39%) 4100 (0.88%)F5 5039 4980 (-1.17%) 4960 (-1.56%) 4960 (-1.56%)F6 5703 5680 (-0.40%) 5680 (-0.40%) 5640 (-1.10%)F7 7065 7060 (-0.07%) 6960 (-1.50%) 6940 (-1.76%)F8 7595 7620 (0.32%) 7520 (-0.99%) 7520 (-0.99%)Table C.4: Formant locations Fi (Hz) and formant deviations ∆Fi/Fi (%) ofthe 2.5Dlow, 2.5Dmid and 2.5Dhigh simulations for vowel /A/ with respectto the 3D FEM vocal tract model.Formants 3D 2.5Dlow 2.5Dmid 2.5DhighF1 263 260 (-1.14%) 260 (-1.14%) 260 (-1.14%)F2 2111 2160 (2.32%) 2140 (1.37%) 2160 (2.32%)F3 3010 3060 (1.66%) 3040 (0.99%) 3040 (0.99%)F4 4138 4220 (1.98%) 4100 (-0.91%) 4120 (-0.43%)F5 5019 5020 (0.01%) 4920 (-1.97%) 4940 (-1.57%)F6 5753 5740 (-0.22%) 5640 (-1.96%) 5620 (-2.31%)F7 6575 6480 (-1.44%) 6560 (-0.22%) 6540 (-0.53%)F8 7645 7680 (0.45%) 7660 (0.19%) 7680 (0.45%)Table C.5: Formant locations Fi (Hz) and formant deviations ∆Fi/Fi (%) ofthe 2.5Dlow, 2.5Dmid and 2.5Dhigh simulations for vowel /i/ with respectto the 3D FEM vocal tract model.147Formants 3D 2.5Dlow 2.5Dmid 2.5DhighF1 259 260 (0.38%) 260 (0.38%) 260 (0.38%)F2 757 720 (-4.88%) 720 (-4.88%) 720 (-4.88%)F3 2264 2300 (1.59%) 2300 (1.59%) 2300 (1.59%)F4 3603 3620 (0.47%) 3600 (-0.08%) 3620 (0.47%)F5 4173 4280 (2.56%) 4260 (2.08%) 4260 (2.08%)F6 5062 4960 (-2.01%) 4960 (-2.01%) 5000 (-1.22%)F7 6126 6080 (-0.75%) 6180 (0.88%) 6180 (0.88%)F8 6594 6620 (0.39%) 6700 (1.60%) 6700 (1.60%)Table C.6: Formant locations Fi (Hz) and formant deviations ∆Fi/Fi (%) ofthe 2.5Dlow, 2.5Dmid and 2.5Dhigh simulations for vowel /u/ with respectto the 3D FEM vocal tract model.Formants 2Dhigh 2.5Dlow 2.5Dmid 2.5DhighF1 820 820 (0.0%) 820 (0.0%) 820 (0.0%)F2 1700 1720 (1.17%) 1720 (1.17%) 1720 (1.17%)F3 2460 2460 (0.0%) 2480 (0.81%) 2500 (1.62%)F4 4320 4380 (1.38%) 4340 (0.46%) 4300 (-0.46%)F5 4760 4800 (0.84%) 4800 (0.84%) 4720 (-0.84%)F6 5820 5880 (1.03%) 5840 (0.34%) 5860 (0.68%)F7 7080 7200 (1.69%) 7180 (1.41%) 7180 (1.41%)F8 8300 8440 (1.68%) 8340 (0.48%) 8320 (0.24%)Table C.7: Formant locations Fi (Hz) and formant deviations ∆Fi/Fi (%) ofthe 2.5Dlow, 2.5Dmid and 2.5Dhigh simulations for vowel /æ/ with respectto the 2Dhigh FDTD vocal tract model.148Appendix D2.5D FDTD Time MarchingAlgorithmAlgorithm 1 FDTD Time-Marching algorithmInput: VT area function a(x), audio time t1: Initialize physical constants: air density ρ , sound speed c, boundary admittance coefficient µ .2: Initialize the simulation sampling rate R.3: Set the temporal resolution (∆t) and grid resolution (∆s) with R and CFL condition.4: Create a optimized grid having size (M×N).5: Set the total number of time steps T of the simulator using t and ∆t.6: Define boundary cells with a(x) and normal acoustic impedance Z.7: Set depth values (D¯, D(x) and D(y)) for each grid cell from a(x).8: Define source excitation cells9: Initialize source excitation velocity ve.10: Initialize acoustic components (p, vx and vy) for each grid cell.11: for n = 1...T do12: for i = 1...M do13: for j = 1...N do14: Update pn+1(i, j) with vnx(i, j), vny(i, j) and D(i, j) (Eq. 3.6 and Eq.3.9)15: if (i, j) = excitation cell then16: vnx ← vnx + vne and vny ← vny + vne17: end if18: if neighbouring cell = boundary cell then19: vnb← pn(i, j)/Z20: end if21: Update vn+1x (i, j) with pn+1(i, j), vnx and vnb (Eq. 3.8)22: Update vn+1y (i, j) with pn+1(i, j),vny and vnb (Eq. 3.8)23: end for24: end for25: end for149
Thesis/Dissertation
eng
Electrical and Computer Engineering
Vancouver : University of British Columbia Library
Attribution-NonCommercial-NoDerivatives 4.0 International
Graduate
Talking tube : a novel approach for vocal tract acoustic modelling using the finite-difference time-domain method
Text