## 2 Methodology
### Vocal Tract Impulse Response  
The impulse response of each vocal tract was obtained from the time-domain FEM simulations [11]. The propagation of acoustic waves within the 3D vocal tracts was provided by the FEM solution to the acoustic wave equation:  
\[\partial_{tt}^{2}p-c_{0}^{2}\
abla^{2}p=0, \tag{1}\]  
where \(p(x,t)\) stands for the acoustic pressure, \(\partial_{tt}^{2}\) for the second order time derivative, and \(c_{0}\) for the speed of sound. \(c_{0}\) was set to the usual value of 350 m/s. An exterior domain was included to let waves emanate from the mouth and account in this way for radiation losses. A perfectly matched layer  
Figure 1: Synthesis of vowels [a], [i], and [u] with realistic vocal tract geometries (above) and their simplified counterparts of circular cross-sections set in a straightened midline (below). The output pressure signal \(p(t)\) is computed as the convolution of the glottal source \(u_{g}(t)\) with the vocal tract impulse response \(h(t)\) obtained from a 3D FEM (finite element method) simulation. Three phonation type examples are represented in the figure: Tense (dashed red line), modal (solid black line), and lax (dotted green line).  
(PML) was imposed on the computational domain boundaries to prevent wave reflections. Wall losses were considered by prescribing a boundary admittance coefficient \(\mu\) on the vocal tract walls which was set to \(\mu=0.005\). Sound waves were generated within the vocal tract imposing a volume velocity \(u_{g}(t)\) on the glottal cross-sectional area. Specifically, the following Gaussian pulse was used:  
\[u_{g}(t)=e^{-\left[(t-T_{gp})/0.29T_{gp}\right]^{2}}[\text{m}^{3}/\text{s}], \tag{2}\]  
with \(T_{gp}=0.646/f_{c}\) and \(f_{c}=10\) kHz.  
Numerical simulations were performed with a sampling frequency of \(f_{s}=8000\) kHz. This unusually large value was selected to ensure the stability of the explicit discrete time scheme used to solve the wave Equation (1). Time events of 20 ms were simulated, tracking the acoustic pressure, \(p_{0}(t)\) at a mesh node located 4 cm away from the mouth exit. The vocal tract transfer function \(H(f)\) was then obtained as:  
\[H(f)=\frac{P_{o}(f)}{U_{g}(f)}, \tag{3}\]  
with \(P_{o}(f)\) and \(U_{g}(f)\) respectively being the Fourier transforms of \(p_{o}(t)\) and \(u_{g}(t)\). Note that this compensates for the slight spectral decay introduced by the Gaussian pulse. \(H(f)\) was computed up to 12 kHz in order to generate speech at 24 kHz. This sampling frequency allowed us to cover the whole 8 kHz octave band, in which the HFE levels would be computed.  
Figure 2 shows the computed vocal tract transfer functions \(H(f)\) of [a], [i], and [u] for the realistic and simplified geometries. Observe that below 5 kHz the two representations behaved very similarly, whereas above that frequency strong differences emerged. This is consistent with the results presented in [5], where it was observed that plane wave propagation, which dominates below 5 kHz, is barely affected by the cross-sectional shape and vocal tract bending. Beyond that limit, however, higher order modes also propagate and play a significant role in the realistic configurations. In contrast, radial symmetry prevents the onset of most higher order modes in the simplified vocal tracts for the examined frequency range [4; 5]. As observed in Figure 2, the depicted vocal tract transfer functions show an almost flat global trend in contrast to the spectral characteristics of speech because they do not include the effect of the glottal source.  
The vocal tract impulse responses \(h(t)\) were finally obtained from the inverse Fourier transform of the vocal tract transfer functions \(H(f)\) (see Figure 1).