# A comparative study of two-dimensional vocal tract acoustic modeling based on Finite-Difference Time-Domain methods

Debasish Ray Mohapatra, Victor Zappi, Sidney Fels

► To cite this version:

Debasish Ray Mohapatra, Victor Zappi, Sidney Fels. A comparative study of two-dimensional vocal tract acoustic modeling based on Finite-Difference Time-Domain methods. International Seminar on Speech Production, Dec 2020, Virtual conference, United States. hal-05174815

HAL Id: hal-05174815

https://hal.science/hal-05174815v1

Submitted on 21 Jul 2025

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
# A comparative study of two-dimensional vocal tract acoustic modeling based on Finite-Difference Time-Domain methods

Debasish Ray Mohapatra¹, Victor Zappi², Sidney Fels¹

¹Electrical and Computer Engineering Department, University of British Columbia, Canada

²College of Arts, Media and Design, Northeastern University, USA

debasishray@ece.ubc.ca, v.zappi@northeastern.edu, ssfels@ece.ubc.ca

## Abstract

The two-dimensional (2D) numerical approaches for vocal tract (VT) modelling can afford a better balance between the low computational cost and accurate rendering of acoustic wave propagation. However, they require a high spatio-temporal resolution in the numerical scheme for a precise estimation of acoustic formants at the simulation run-time expense. We have recently proposed a new VT acoustic modelling technique, known as the 2.5D Finite-Difference Time-Domain (2.5D FDTD), which extends the existing 2D FDTD approach by adding tube depth to its acoustic wave solver. In this work, first, the simulated acoustic outputs of our new model are shown to be comparable with the 2D FDTD and a realistic 3D FEM VT model at a low spatio-temporal resolution. Next, a radiation model is developed by including a circular baffle around the VT as head geometry. The transfer functions of the radiation model are analyzed using five different vocal tract shapes for vowel sounds /a/, /e/, /i/, /o/ and /u/.

**Keywords:** computational acoustics, vocal tract, FDTD, articulatory speech synthesis

## 1. Introduction

Various articulatory models have been developed to approximate the sophisticated upper vocal tract geometry and capture its acoustic features using a physics-based acoustic wave solver. The 3D acoustic analysis (Takemoto, Mokhtari, and Kitamura 2010; Vampola et al. 2015) can precisely characterize the vocal tract's spectral features while offering better geometrical flexibility. However, such models are not pragmatic for designing real-time speech synthesizers due to their high computational cost. In contrast to 3D, the 1D vocal tract models are computationally lightweight. Nevertheless, their over-simplified representation of a tube structure substantially eliminates its geometrical details. Hence, the acoustic analysis of 1D models is not accurate at higher formants.

As an alternative, the 2D acoustic wave solvers (Arnela and Guasch 2014; Speed, Murphy, and Howard 2009) represent the complex vocal tract geometry as a 2D contour (a mid-sagittal cut of a 3D VT geometry) and capture the wave interaction only along the mid-sagittal plane. Due to the dimensionality reduction, the 2D models are lightweight compare to 3D ones. The most well-known 2D acoustic analysis methods for vocal tract modelling are the finite element method (FEM) and FDTD. Since 2D models do not lump off-plane wave interaction, direct use of cross-sectional areas from the vocal tract MRI images to represent 2D contours inside a simulation grid yields erroneous acoustic output. Alternatively, a nonlinear area func-

tion transformation method needs to be implemented to tune the acoustic features of 2D models that can match the more complex 3D ones. However, this procedure kills the performance boost of 2D models. Consequently, we have proposed a new vocal tract acoustic model (2.5D FDTD) that improves the 2D model by capturing the wave interaction across the mid-sagittal plane. The next section briefly discusses the numerical implementation of a 3D tube using 2D and 2.5D FDTD schemes.

## 2. Methodology

### 2.1. Acoustic Wave Solver

We implemented the FDTD numerical analysis technique to discretize the 2D acoustic wave equation on a staggered Yee grid. The 2D Yee scheme consists of a rectangular grid, define acoustic parameters (pressure ($p^{(n)}$) and velocity ($v_x^{(n)}, v_y^{(n)}$)) at each grid point, a squared 2D cell (Allen and Raghuvanshi 2015; Yee 1966), as shown in Figure 1. However, the employed technique for acoustic simulation does not facilitate dynamic boundary conditions. Hence, a new scalar field $0 \le \beta(x, y, t) \le 1$ was introduced to the solver, which could transit smoothly between $\beta = 1$ (air) and $\beta = 0$ (boundary). At $\beta = 0$, a prescribed velocity $v = v_b$ was enforced to handle the boundary condition. Specifically, per each time step ($n$), the wave solver updates pressure component ($p^{(n)}$) and velocity components ($v_x^{(n)}, v_y^{(n)}$) at each grid point, by solving the below linear wave equations in the time domain (Zappi et al. 2016):

$$ \frac{\partial p}{\partial t} + (1 - \beta)p = -\rho c^2 \left( \frac{\partial v_x}{\partial x} + \frac{\partial v_y}{\partial y} \right) \quad (1) $$

$$ \beta \frac{\partial v}{\partial t} + (1 - \beta) v = -\beta^2 \frac{\nabla p}{\rho} + (1 - \beta) v_b \quad (2) $$

The 2.5D FDTD follows the above rationale but improves upon it by incorporating new impedance terms in its 2D acoustic wave solver, known as tube depth. The tube's depth ($D$), i.e., its continuous extension along the $z$ axis (with $x$ and $y$ axis being the dimension of starting 2D scheme), is derived from the cross-sectional area of the tube and sampled at each point of the scheme as shown in Figure 1. Then the resulting depths are mapped to their respective acoustic parameters ($p, v_x, v_y$) as demonstrated by Mohapatra, Zappi, and Fels (2019). The inclusion of tube depth allows the wave solver to apprehend wave interactions across the mid-sagittal plane, in turn eliminating the nonlinear cross-sectional area functions transformation requirement. Therefore, it is expected that the time complexity of a 2.5D FDTD tube model will be better than the existing 2D models. Basically, with the help of following discrete update rules, we adopt a time marching algorithm (see Algorithm 1) to sample the acoustic parameters at every grid point:
<figure><img src="image_2.png" /><figcaption>Figure 1: Representation of acoustic parameters ($p$, $v_x$, $v_y$) and tube depth ($\bar{D}$, $D_{(x)}$ and $D_{(y)}$) inside the 2D Yee scheme</figcaption></figure>

$$p^{(n+1)} = \frac{\bar{D}p^{(n)} - \rho c^2 \Delta t \tilde{\nabla} \cdot \mathbf{V}^{(n)}}{\bar{D}} \quad (3)$$

$$\mathbf{v}^{(n+1)} = \frac{\beta \mathbf{v}^{(n)} - \beta^2 \Delta t \tilde{\nabla} p^{(n+1)} / \rho + \Delta t (1 - \beta) \mathbf{v}_b}{\beta + \Delta t (1 - \beta)} \quad (4)$$

with:

$$\mathbf{V} = (D_{(x)}v_x, D_{(y)}v_y) \quad (5)$$

## Algorithm 1 FDTD Time-Marching algorithm

Input: VT area function $a(x)$, audio time $t$

1: Initialize physical constants: air density $\rho$, sound speed $c$, boundary admittance coefficient $\mu$.

2: Initialize the simulation sampling rate $R$.

3: Set the temporal resolution ($\Delta t$) and grid resolution ($\Delta s$) with $R$ and CFL condition.

4: Create an optimized grid having size ($M \times N$).

5: Set the total number of time steps $T$ of the simulator using $t$ and $\Delta t$.

6: Define boundary cells with $a(x)$ and normal acoustic impedance $Z$.

7: Set depth values ($\bar{D}$, $D_{(x)}$ and $D_{(y)}$) for each grid cell from $a(x)$.

8: Define source excitation cells

9: Initialize source excitation velocity $v_e$.

10: Initialize acoustic components ($p$, $v_x$ and $v_y$) for each grid cell.

11: **for** $n = 1...T$ **do**

12:     **for** $i = 1...M$ **do**

13:         **for** $j = 1...N$ **do**

14:             Update $p^{n+1}(i, j)$ with $v_x^n(i, j)$, $v_y^n(i, j)$ and $D(i, j)$ (Eq. 3 and Eq.5)

15:                 **if** $(i, j)$ = excitation cell **then**

16:                 $v_x^n \leftarrow v_x^n + v_e^n$ and $v_y^n \leftarrow v_y^n + v_e^n$

17:             **end if**

18:             **if** $(i, j)$ = boundary cell **then**

19:                 $v_b^n \leftarrow p^n(i, j)/Z$

20:             **end if**

21:         Update $v_x^{n+1}(i, j)$ with $p^{n+1}(i, j)$ and $v_b^n$ (Eq. 4)

22:         Update $v_y^{n+1}(i, j)$ with $p^{n+1}(i, j)$ and $v_b^n$ (Eq. 4)

23:     **end for**

24:     **end for**

25: **end for**

## 2.2. Vocal Tract Wall Losses

To simulate vocal tract wall reflection, we adopted the local reactive boundary approach as proposed here (Yokota, Sakamoto, and Tachibana 2002; Takemoto, Mokhtari, and Kitamura 2010). The method estimates wall losses (locally reacting soft walls) and enforces boundary condition by using $v_w$, the normal particle velocity going into or coming from a wall.

$$v_w = \frac{p_w}{Z_n} \hat{n} \quad (6)$$

where $p_w$ is the current pressure value in an air cell, located in front of the wall cell. And $Z_n$, the normal acoustic impedance can be computed by using the normal sound absorption coefficient $\alpha_n$ as follows,

<figure><img src="image_3.png" /><figcaption>Figure 2: (a) Simplified head geometry representation with circular mouth aperture (b) Snapshot of acoustic wave propagation inside the vocal tract tube and emanating from mouth aperture at time instant $t = 6.1ms$ for vowel /a/.</figcaption></figure>

$$Z_n = \rho c \frac{1 + \sqrt{1 - \alpha_n}}{1 - \sqrt{1 - \alpha_n}} \quad (7)$$

The boundary absorption coefficient $\alpha_n$ can be derived from the boundary admittance coefficient $\mu \in (0, 1)$. During the simulation, the solver combines equation (6) with equation (4) by setting $v_b = v_w$ for all the velocity vectors which are positioned in between an air cell and wall cell. However, this method does not comprehend frequency dependent wall losses.

## 2.3. Excitation Model

For acoustic analysis and speech production, articulatory models require a source excitation function just above the glottal end to input acoustic energy into the tube (Mohapatra and Fels 2018). In this work, we are only characterizing the acoustic behaviour of a vocal tract tube through its transfer function analysis. Hence, a band-passed velocity pulse having frequency range 2Hz-20kHz was injected near the glottis for various vocal tract shapes, and the corresponding transfer functions (impulse response) were obtained at the mouth-end.

## 2.4. Radiation Model

The outward wave propagation at the mouth-end (lips) contributes to acoustic energy losses in the vocal tract. It is one of the essential dissipation mechanisms in voice production. Hence, a nonzero radiation model of a vocal tract influences the synthesized audio output by lowering the formant frequencies and increasing their bandwidths. It might be challenging to build a model that can precisely attain the radiation effects, considering a realistic human head geometry. However, it has already been demonstrated that the inclusion of a spherical baffle (a simplified approximation of the actual head geometry) around the 3D vocal tract model offers promising results for vowel production (Arnela, Guasch, and Alias 2013). Since the 2D model uses a mid-sagittal cut of a 3D vocal tract, the spherical baffle can be replaced by a circular baffle for the 2.5D radiation model as shown in Figure 2(a). Moreover, to illustrate acoustic waves’ propagation emanating from the mouth toward infinity (Figure 2(b)), the computational domain can be extended out of the vocal tract by using perfectly matched layers (PMLs). PMLs help in eliminating spurious wave reflection at the domain boundaries.

We create a fixed-sized circular baffle with a diameter of 0.20m around the 2.5D vocal tract contour inside the simulation domain. The circular baffle center has been adjusted to fit in the entire vocal tract contour while connecting its exit cross-section. Inside the 2D scheme, the region enclosed by the vocal tract and baffle boundary is set to air cells. In order to absorb
<figure><img src="image_4.png" /><figcaption>Figure 3: 2.5D vocal tract tube representation for vowel /a/.</figcaption></figure>

outgoing waves and approximate infinite space, 6 layers of PML are added at all sides of the domain boundary.

### 3. Experimental Setup

A comparative study was carried out for two different scenarios to characterize the acoustic features of the 2.5D vocal tract model: (1) numerical simulation with computational domain ended at the mouth aperture by imposing Dirichlet boundary conditions (open-end boundary condition) and (2) numerical simulation with free-field radiation, which includes head geometry and PML layers (radiation condition). These conditions were examined using the transfer functions analysis approach for five different vocal tract shapes(vowel sounds: /a/, /e/, /i/, /o/ and /u/). The 2.5D vocal tract shapes were derived using Story's area function dataset (Story 2008), as it is the standard one for the vocal tract acoustic analysis. We have only considered the vocal tract geometry as a straight tube with circular cross-sections for this study as shown in Figure 3.

For open-end boundary conditions, we use a precise 3D FEM vocal tract model as the reference for the transfer function analysis. A virtual microphone is placed 3mm inside the mouth opening to record pressure variations. In the case of free-field radiation, we compare the transfer functions of the 2.5D model with Story's calculated formants, obtained from the frequency-domain analysis of vocal tracts (Sondhi and Schroeter 1987). And the microphone is positioned 3mm outside of the mouth to capture the radiating acoustic waves. Courant-Friedrichs-Lewy (CFL) stability condition in two dimensions: $\Delta t \le \Delta s/\sqrt{2c}$ is imposed, where c is the speed of sound.

The simulation generates 50ms of synthesized audio for every vowel sounds. During the simulation, we set the following physical parameters fixed: air density $\rho = 1.14$ kg/m$^2$s, boundary admittance $\mu = 0.005$ and sound speed c = 350 m/s. We implement both the 2D and 2.5D vocal tract models in the MATLAB environment. The simulation runs on a workstation equipped with an Intel Core i7-8700K processor. The custom code for the 2.5D vocal tract model is publicly available here<fnref n="1" />.

### 4. Model Validation

We follow the below steps to analyze the acoustic behaviour of the 2.5D model:

(A) Transfer function analysis with open-end boundary condition for the 2D FDTD and 2.5D FDTD vocal tract models (vowel sound: /a/, /i/ and /u/).

(B) Model performance evaluation (simulation run-rime).

(C) Transfer function analysis for the 2.5D vocal tract model with free-field radiation (vowel sound: /a/, /e/, /i/. /o/ and /u/).

For Step (A) and (B), we simulate each vowel sound under three different spatial grid resolution for both 2D and 2.5D FDTD methods: low ($\Delta s = 0.74$mm), mid ($\Delta s = 0.28$mm) and high ($\Delta s = 0.18$mm). The high spatial resolution should produce precise acoustic output as it yields vocal tract geometry with greater details. However, it requires a larger simulation domain size. The grid resolution values are empirically derived. Step (A) and (B) examine the improvement in a 2.5D vocal tract model over its 2D simulation due to the addition of tube depth. Step (C) investigates the change in transfer function with the inclusion of a circular baffle for the radiation model.

## 5. Result

We applied Fast Fourier Transformation (FFT) to the recorded pressure waves to obtain their transfer function. Formants were extracted by considering local maxima in transfer function curves. Here we are only analyzing the first three formants since they are mainly responsible for distinguishing vowel sounds (Vampola et al. 2015). The formants for the 3D FEM model can be found here (Arnela Coll 2015).

### 5.1. Step A: Transfer Function Analysis

Transfer function analysis for the 2D and 2.5D simulation with open-end boundary condition<fnref n="2" />,<fnref n="3" />.

<table><caption>Table 1: Positional errors of the first three formants computed for vowel sound /a/, with respect to 3D FEM model.</caption><thead><tr><th>Resolution (Δs)</th><th>F1</th><th>F2</th><th>F3</th></tr></thead><tbody><tr><td>low (% Error)</td><td>660(-5.16)<br>700(0.57)</td><td>1040(-2.61)<br>1040(-2.61)</td><td>2960(-2.33)<br>3020(-0.36)</td></tr><tr><td>mid (% Error)</td><td>660(-5.16)<br>680(-2.29)</td><td>1040(-2.61)<br>1060(-0.74)</td><td>3000(-1.02)<br>3060(0.95)</td></tr><tr><td>high(% Error)</td><td>700(0.57)<br>680(-2.29)</td><td>1060(-0.74)<br>1040(-2.62)</td><td>3020(-0.36)<br>3040(0.29)</td></tr></tbody></table>

<table><caption>Table 2: Positional errors of the first three formants computed for vowel sound /i/, with respect to 3D FEM model.</caption><thead><tr><th>Resolution (Δs)</th><th>F1</th><th>F2</th><th>F3</th></tr></thead><tbody><tr><td>low (% Error)</td><td>240(-8.74)<br>260(-1.14)</td><td>2120(0.42)<br>2160(2.32)</td><td>3020(0.33)<br>3060(1.66)</td></tr><tr><td>mid (% Error)</td><td>260(-1.14)<br>260(-1.14)</td><td>2140(1.37)<br>2140(1.37)</td><td>3000(-0.33)<br>3040(0.99)</td></tr><tr><td>high(% Error)</td><td>260(-1.14)<br>260(-1.40)</td><td>2140(1.37)<br>2160(2.32)</td><td>3020(0.33)<br>3040(0.99)</td></tr></tbody></table>

<table><caption>Table 3: Positional errors of the first three formants computed for vowel sound /u/, with respect to 3D FEM model.</caption><thead><tr><th>Resolution (Δs)</th><th>F1</th><th>F2</th><th>F3</th></tr></thead><tbody><tr><td>low (% Error)</td><td>260(0.38)<br>260(0.38)</td><td>700(-7.52)<br>720(-4.88)</td><td>2260(-0.17)<br>2300(1.59)</td></tr><tr><td>mid (% Error)</td><td>220(-15.0)<br>260(0.38)</td><td>600(-20.7)<br>700(-7.52)</td><td>2260(-0.17)<br>2300(1.59)</td></tr><tr><td>high(% Error)</td><td>260(0.38)<br>260(0.38)</td><td>700(-7.52)<br>720(-4.88)</td><td>2280(0.70)<br>2300(1.59)</td></tr></tbody></table>

<footnote n="1">¹https://github.com/Debasishray19/Talking-Tube</footnote>

<footnote n="2">²Colormap: For 2D simulation, For 2.5D simulation.</footnote>

<footnote n="3">³All the formant frequencies are in Hz.</footnote>
## 5.2. Step B: Model Performance

The computational cost of FDTD vocal tract model depends upon the spatial grid resolution and simulation domain size. We implement an algorithm that enforces an optimized domain size for FDTD models. The table below presents the simulation domain size and time duration corresponding to each grid resolution for vowel /u/.

<table><caption>Table 4: Simulation run-time at different grid resolutions.</caption><thead><tr><th>Resolution (Δs)</th><th>Domain Size (width×height)</th><th>Duration (Approx.) (in seconds)</th></tr></thead><tbody><tr><td>low (0.74mm)</td><td>266 × 35</td><td>280</td></tr><tr><td>mid (0.28mm)</td><td>702 × 85</td><td>4.98e+3</td></tr><tr><td>high (0.18mm)</td><td>1051 × 127</td><td>1.85e+4</td></tr></tbody></table>

## 5.3. Step C: Transfer Functions Analysis

The table below presents transfer function analysis for the 2.5D vocal tract models with the free-field radiation. We have only considered the low grid resolution (Δs = 0.74mm) for simulation of vocal tract radiation models.

<table><caption>Table 5: Formant frequencies and the percentage positional errors of the first three formants for the vocal tract radiation model.</caption><thead><tr><th>Vowel</th><th>F1</th><th>F2</th><th>F3</th></tr></thead><tbody><tr><td>/a/</td><td>640<br/>-7.78%</td><td>1020<br/>8.28%</td><td>3000<br/>1.76%</td></tr><tr><td>/e/</td><td>340<br/>-29.75%</td><td>1920<br/>1.74%</td><td>2320<br/>-3.49%</td></tr><tr><td>/i/</td><td>260<br/>-20.00%</td><td>2080<br/>-2.75%</td><td>2960<br/>-0.53%</td></tr><tr><td>/o/</td><td>480<br/>-3.80%</td><td>760<br/>-1.04%</td><td>2340<br/>-3.34%</td></tr><tr><td>/u/</td><td>260<br/>-17.19%</td><td>720<br/>2.56%</td><td>2180<br/>-5.13%</td></tr></tbody></table>

# 6. Discussion & Conclusion

The transfer function analysis is a standard practice to validate acoustic features of area-based articulatory models. In the first step of our evaluation, the 2.5D FDTD simulation produces positional errors below 2% at a low grid resolution. Moreover, this characteristic is consistent across all three vowels. However, the positional error is high for the 2D FDTD simulation at a low grid resolution in most cases. **Table 1, 2 and 3** show that the acoustic outputs of the 2D FDTD model can be tuned by increasing its computational grid resolution whereas this approach doesn't make much difference for the 2.5D FDTD simulation. **Table 4** shows for a real-time speech synthesizer the low grid resolution is always desirable.

In the case of free-field of radiation, **Table 5** presents the first three formants and their positional errors. The positional error of the first formant is high across all the vowels. And the second and third formants produce relatively lesser error values. However, these positional errors are determined relative to Story's 1D vocal tract model, an oversimplified representation of the actual vocal tract. Hence, as future work, the 2.5D radiation model's acoustic characteristics need to be compared with a 3D vocal tract radiation model.

Additionally, our current 2.5D model should only be seen as simple approximations of the voice radiation mechanism. Nev-

ertheless, the influence of few head details like the nose and cross-section of the mouth aperture is worth studying. As well as, our current 2.5D FDTD model does not include boundary losses for the tube depth. All these current limitations need to be incorporated for the future development of the 2.5D vocal tract acoustic model.

# 7. Acknowledgements

This work is supported by the Natural Sciences and Engineering Research Council (NSERC).

# 8. References

Allen, Andrew and Nikunj Raghuvanshi (2015). "Aerophones in flat-land: Interactive wave simulation of wind instruments". In: *ACM Transactions on Graphics* (TOG) 34.4, pp. 1–11.

Arnela Coll, Marc (2015). "Numerical production of vowels and diphthongs using finite element methods". PhD thesis. Universitat Ramon Llull.

Arnela, Marc and Oriol Guasch (2014). "Two-dimensional vocal tracts with three-dimensional behavior in the numerical generation of vowels". In: *The Journal of the Acoustical Society of America* 135.1, pp. 369–379.

Arnela, Marc, Oriol Guasch, and Francesc Alias (2013). "Effects of head geometry simplifications on acoustic radiation of vowel sounds based on time-domain finite-element simulations". In: *The Journal of the Acoustical Society of America* 134.4, pp. 2946–2954.

Mohapatra, Debasish Ray and Sidney Fels (2018). "Limitations Of Source-Filter Coupling In Phonation". In: *arXiv preprint arXiv:1811.07435*.

Mohapatra, Debasish Ray, Victor Zappi, and Sidney Fels (2019). "An extended two-dimensional vocal tract model for fast acoustic simulation of single-axis symmetric three-dimensional tubes". In: *arXiv preprint arXiv:1909.09585*.

Sondhi, Man and Juergen Schroeter (1987). "A hybrid time-frequency domain articulatory speech synthesizer". In: *IEEE Transactions on Acoustics, Speech, and Signal Processing* 35.7, pp. 955–967.

Speed, Matt, Damian Murphy, and David M Howard (2009). "Characteristics of two-dimensional finite difference techniques for vocal tract analysis and voice synthesis". In: *Tenth Annual Conference of the International Speech Communication Association*.

Story, Brad H (2008). "Comparison of magnetic resonance imaging-based vocal tract area functions obtained from the same speaker in 1994 and 2002". In: *The Journal of the Acoustical Society of America* 123.1, pp. 327–335.

Takemoto, Hironori, Parham Mokhtari, and Tatsuya Kitamura (2010). "Acoustic analysis of the vocal tract during vowel production by finite-difference time-domain method". In: *The Journal of the Acoustical Society of America* 128.6, pp. 3724–3738.

Vampola, Tomáš, Jaromír Horáček, Anne-Maria Laukkanen, and Jan G Švec (2015). "Human vocal tract resonances and the corresponding mode shapes investigated by three-dimensional finite-element modelling based on CT measurement". In: *Logopedics Phoniatrics Vocology* 40.1, pp. 14–23.

Yee, Kane (1966). "Numerical solution of initial boundary value problems involving Maxwell's equations in isotropic media". In: *IEEE Transactions on antennas and propagation* 14.3, pp. 302–307.

Yokota, Takatoshi, Shinichi Sakamoto, and Hideki Tachibana (2002). "Visualization of sound propagation and scattering in rooms". In: *Acoustical science and technology* 23.1, pp. 40–46.

Zappi, Victor, Arvind Vasuvedan, Andrew Allen, Nikunj Raghuvanshi, and Sidney Fels (2016). "Towards real-time two-dimensional wave propagation for articulatory speech synthesis". In: *Proceedings of Meetings on Acoustics* 171ASA. Vol. 26. 1. ASA, p. 045005.