Numerical simulations of head-related transfer functions (HRTFs) conventionally assume a rigid boundary condition for the pinna. The human pinna, however, is an elastic deformable body that can vibrate due to incident acoustic waves. This work investigates how sound-induced vibrations of the pinna can affect simulated HRTF magnitudes. The work will motivate the research question by measuring the sound-induced vibrational patterns of an artificial pinna with a high-speed holographic interferometric system. Then, finite element simulations are used to determine HRTFs for a tabletop model of the B&K 5128 head and torso simulator for a number of directions. Two scenarios are explored: one where the pinna is modeled as perfectly rigid, and another where the pinna is modeled as linear elastic with material properties close to that of auricular cartilage. The findings suggest that pinna vibrations have negligible effects on HRTF magnitudes up to 5 kHz. The same conclusion, albeit with less certainty, is drawn for higher frequencies. Finally, the importance of the elastic domain's material properties is emphasized and possible implications for validation studies on dummy heads 1as well as the limitations of the present work are discussed in detail.

Human sound perception is inherently spatial1 and relies on acoustical cues, which can be mathematically described/encoded as head-related transfer functions (HRTFs). Such geometric-based audio filters are highly individualized and thus expensive to acquire accurately (see, e.g., Ref. 2). The two primary approaches for obtaining reference HRTF estimations are either to measure them experimentally or to simulate them using numerical simulation techniques.3–6 

Numerical simulations within virtual acoustics, including simulations of HRTFs, generally employ wave-based models, which are believed to describe linear acoustics reasonably well.7 Most often, an impedance termination of local reaction is used on the domain boundary: this approximation was introduced about a century ago to study the sound decay and modes inside a room8,9 and presents an abundance of issues when more detailed acoustical predictions are needed.10–12 Moreover, numerical simulations of HRTFs commonly employ fully rigid boundaries:13–15 an unphysical assumption that does not even correspond to the skin impedance measured for anatomical structures way simpler than the human pinnae.16,17 For example, the effects of considering a non-rigid head material with a finite specific acoustic impedance, to account for human hair, has been previously studied.18 The findings support the hypothesis that more realistic material properties could improve the accuracy of simulated HRTF data.

Although some efforts have been made to improve the modeling approaches that use the local reaction assumption (e.g., Ref. 19), the boundary impedance model itself is, in general, an approximation, valid in limited scenarios. (For instance, local reaction is expected to be valid for highly dense solids and highly porous materials.20) It is currently unclear how such modeling assumptions (i.e., the rigid and locally reacting boundary conditions) affect the simulated HRTFs. Errors associated with these modeling approaches can contribute to the total validation error of the simulated HRTFs.21 

This work aims to improve the accuracy of HRTF simulations by modeling the pinna as an elastic, deformable body that interacts with acoustic pressure waves in air. The objective is to evaluate the validity of using a rigid boundary condition for the pinna region and consequently improve the accuracy of physics-based HRTF simulations to enhance future validation studies.

To motivate the research question, we, first, experimentally characterize the sound-induced vibrations of an average anthropometric pinna replica (KB5000; GRAS Sound & Vibration, Holte, Denmark) using a customized high-speed holographic interferometric system.22 We show that, for a genus 0 surface such as the human pinna, the sound-induced pinna vibrations are measurable and larger than 0, in contrast to the rigid boundary assumption. To evaluate the effect of such vibrations on the resulting HRTFs, we then develop physics-based elasto-acoustic models that capture the interaction of the acoustic waves with the pinna structure. The resulting HRTFs will be compared with simulated HRTFs obtained from a pure acoustic model using the conventional perfectly rigid boundary condition assumption for the pinna.

Irrespective of modeling assumptions, any numerical simulation contains numerical errors.23 To have a meaningful comparison between the numerical results of the two models, we have to distinguish between the differences caused by using different modeling assumptions (acoustic versus elasto-acoustic) versus the differences caused due to numerical errors. To make this distinction, we estimate the magnitude of the numerical errors through verification exercises. Since, to the best of our knowledge, such exercises are missing for the employed methods and models, we design and conduct the necessary verification exercises prior to comparing the two physics-based models. To this end, controlled convergence studies are conducted for both sets of models to determine a metric for the reliability of the numerical solutions. Finally, we discuss the implications of the findings on both the design and numerical modeling of head and torso simulators (HATS). To the best of our knowledge, this will be the first experimental study visualizing the sound-induced vibrational patterns of an artificial pinna and the first multiphysics modeling investigation for simulating blocked-meatus HRTFs.

In this section, we show how an acoustic source of a reasonable intensity induces measurable vibrations in an artificial anthropometric model of a human pinna (KB5000; GRAS Sound & Vibration, Holte, Denmark) with a hardness of 35 Shore 00.

Vibration measurements are done using a full-field high-speed holographic interferometry (FHHI) technique. This technique has been previously used to measure vibrations in various applications, including sound-induced motions of the tympanic membrane (TM) in the nanometer range.22,24–28 To accurately measure the evolution of nanometer-scale deformations on the complex geometry of the pinna, full-field measurements with high spatial resolution are required. Due to large low-frequency motion artifacts driven by physiological motions such as heart beat, head muscle activation, respiration, etc., when measuring on live subjects, we used the GRAS KB5000 anthropometric pinna to perform the sound-induced displacement measurements.

Figure 1(a) shows schematics of the experimental setup. The pinna is fixed at the back of its rectangular base and is mounted on a 1.5-in. Thorlabs post (Thorlabs, Newton, NJ). Our FHHI system uses a 75 mW near-infrared tunable laser (TLB6700; New Focus, San Jose, CA) with a central wavelength of 780 nm. To digitally capture the holograms with the required acquisition speed, a high-speed camera (SA-Z; Photron, Tokyo, Japan) is used. The camera has a 1 MPixel resolution (i.e., 1024 × 1024 pixels with a pixel size of 20 μm). The FHHI method relies on quantification of the optical path length variations due to the sample's motion using a phase-correlation method.22 As shown in Fig. 1(b), the FHHI setup consists of an off-axis interferometer that splits a high-coherence laser beam into a reference path and an object path. The reference path is equipped with a mirror mounted on a high-speed precision piezoelectric actuator (PI, Auburn, AL) that can vary the optical path length by a calibrated amount. The reflection from the sample is combined and interfered with the reference beam at the complementary metal-oxide semiconductor (CMOS) sensor of the high-speed camera. The excitation source (a Neumann KH120A loudspeaker; Neumann, Berlin, Germany) and the acquisition triggering of the camera are all synchronized using the master clock of a data acquisition system (DAQ USB 6366; NI, Austin, TX).

FIG. 1.

(Color online) An overview of the developed high-speed holographic system. (a) Schematic of the system and its key components. (b) Camera's field of view (FOV) showing sound-induced vibrational patterns of the artificial pinna, excited with a tonal stimuli of 500 Hz. The perimeter of the concha bowl is highlighted.

FIG. 1.

(Color online) An overview of the developed high-speed holographic system. (a) Schematic of the system and its key components. (b) Camera's field of view (FOV) showing sound-induced vibrational patterns of the artificial pinna, excited with a tonal stimuli of 500 Hz. The perimeter of the concha bowl is highlighted.

Close modal

Figure 1(b) shows a representative image of the sound-induced vibrational patterns of the artificial ear, observed from the camera's field of view.

The phase quantification approach requires spatial smoothness of optical interferometric phase (due to sample motions). The method performs correlation of speckle patterns at a 5 × 5 pixel kernel around each point before and after a series of continuously phase-shifted reference frames. A set of continuously phase-shifted reference frames are scanned within 3 ms to precisely acquire a reference frame with a π / 2 phase shift before the sample is acoustically excited. The change in the interference phase Δ ϕ at each time t R + and at each coordinate point ( x , y ) R 2 on the hologram will be determined, in radians, by
Δ ϕ ( x , y , t ) = arctan ρ ( I ref ( π / 2 ) ( x , y ) , I def ( x , y ) ) ρ ( I ref ( x , y ) , I def ( x , y ) ) ,
(1)
where ρ denotes the Pearson's correlation between two data sets that is corresponding to the cosine of the interference phase [i.e., ρ ( I ref , I def ) cos ( Δ ϕ ) ] , and I ref ( π / 2 ) ( x , y ) is a kernel centered at the desired point (x, y) having a π / 2 phase shift with respect to frame I ref ( x , y ). The I ref ( x , y ) is the last frame before the beginning of the transient sample excitation. The “ref” and “def” sub-indices represent the reference and deformed states of the sample in time, respectively. Further details of this method are well explained in previous publications.22 This change in the interference phase, Δ ϕ, is related to the displacement field, d(x, y, t), through Eq. (2),
d ( x , y , t ) = ( λ / 4 π ) × Δ ϕ ( x , y , t ) ,
(2)
where λ is the wavelength of the laser light (in meters) being used.

This section provides examples of full-field vibrational measurements of a pinna induced by a loudspeaker and demonstrates how an ordinary acoustical source can induce measurable vibrational patterns in the pinna structure. To achieve this, the artificial pinna (GRAS KB5000) is acoustically excited with a loudspeaker that is placed at a distance of 1 m from the pinna, with an amplitude of nearly 80 dB sound pressure level (SPL) at the pinna location. Figure 2 shows representative results of acoustically induced pinna vibrations at different tonal stimuli and vibration phase angles with respect to the excitation signal. The top row shows the modulation intensity images determined by the absolute value of the complex optical phase. The middle row shows the corresponding wrapped mod-2π optical phase obtained by Eq. (1); and the bottom figure shows the sound-induced displacement fields corresponding to the unwrapped and scaled optical phase (Eq. (2).

FIG. 2.

(Color online) Representative sound-induced pinna vibration patterns for tonal excitation frequencies of 500, 1000, and 2000 Hz. For each of these excitation frequencies, the shown vibration map corresponds to the difference of the motion between vibration phases at 0° and 90°. The top row shows the modulation intensity data [i.e., the absolute value of the complex interferometric data (Ref. 30)] the middle row shows the mod-2π change in the optical phase induced by the sound-induced displacement as obtained by Eq. (1), and the bottom row shows the unwrapped scaled displacement data. Each of the columns depicts vibration data for a given tonal frequency of 500, 1000, and 2000 Hz. Submitted supplementary material videos for full-field pinna vibrations excited at these three tonal frequencies show the presence of both traveling and standing waves on the pinna due to incident acoustic waves.

FIG. 2.

(Color online) Representative sound-induced pinna vibration patterns for tonal excitation frequencies of 500, 1000, and 2000 Hz. For each of these excitation frequencies, the shown vibration map corresponds to the difference of the motion between vibration phases at 0° and 90°. The top row shows the modulation intensity data [i.e., the absolute value of the complex interferometric data (Ref. 30)] the middle row shows the mod-2π change in the optical phase induced by the sound-induced displacement as obtained by Eq. (1), and the bottom row shows the unwrapped scaled displacement data. Each of the columns depicts vibration data for a given tonal frequency of 500, 1000, and 2000 Hz. Submitted supplementary material videos for full-field pinna vibrations excited at these three tonal frequencies show the presence of both traveling and standing waves on the pinna due to incident acoustic waves.

Close modal

To obtain the root mean square (RMS) of the displacement field at each tonal stimulus, the sound-induced displacement versus excitation phase results, gathered at each of the image pixels, were Fourier transformed.29 The magnitudes and phases of the fundamental components were then computed (Fig. 3, left and center columns, respectively). A video of sound-induced pinna vibrations excited with a number of tonal stimuli, captured by our holographic system, is also provided in the supplementary material. The video represents the absolute cosine of the interference phase with a fixed scale of 0–1 for all frequencies. To evaluate the quality of the reconstructed magnitude and phase of the sound-induced displacement fields, we compute the Pearson product correlation coefficient between the measured phasic variations of the displacement and the sinusoidal motion reconstructed from the magnitude and angle of the Fourier-determined fundamental component of the phasic motion. As shown in the right column of Fig. 3, the square of this correlation, R2, is generally greater than 0.9, except at locations where the magnitudes of the measured displacements are less than 15 nm. This indicates that the displacement maps are accurately reconstructed.

FIG. 3.

(Color online) Representative quantitative pinna vibrations for 500, 1000, and 2000 Hz. The magnitudes (left column) and phases (center column) of the fundamental vibration components were computed using the fast Fourier transform (FFT) approach; The right column shows the Pearson product correlation coefficient between the measured phasic displacement and the sinusoidal motion reconstructed from the magnitude and angle of the Fourier-determined fundamental component of the phasic motion.

FIG. 3.

(Color online) Representative quantitative pinna vibrations for 500, 1000, and 2000 Hz. The magnitudes (left column) and phases (center column) of the fundamental vibration components were computed using the fast Fourier transform (FFT) approach; The right column shows the Pearson product correlation coefficient between the measured phasic displacement and the sinusoidal motion reconstructed from the magnitude and angle of the Fourier-determined fundamental component of the phasic motion.

Close modal

The left column of Fig. 3 clearly shows that the sound-induced displacement of the pinna is on the order of a few hundred nanometers. This indicates that the pinna vibrates under environmental acoustic stimuli. The effect of such small displacements on HRTFs is unclear; however, results in Fig. 3 question the validity of the rigid boundary condition assumption for the pinna when it comes to HRTF simulations. To evaluate the magnitude of these effects methodologically, we have developed finite element models that further explore this question.

In this section, we present the general physics-based mathematical formulation for an elastic scattering problem. The scattering problem for a rigid boundary can then be formalized as a specific case of elastic scattering. Here, we will assume small displacements together with a linear elastic material model. The models apply to isotropic and homogeneous domains. In what follows, vector and tensor-valued functions are represented by bold lowercase and uppercase letters, respectively. In this section, we use a Cartesian coordinate system in R 3. (Note that HRTFs are considered in a spherical coordinate system.) Consider an elastic body encompassing the bounded solid domain Ω s R 3 ( Γ s = Ω s R 2 is the boundary of Ωs) and impinged by a time harmonic incident wave, p inc (Fig. 4). We consider x = [ x 1 , x 2 , x 3 ] T to be the position variable and ω to represent the angular frequency. In the frequency domain, the total acoustic wave field p ( x ) is
p ( x ) = p inc ( x ) + p s ( x ) ,
(3)
where p inc is the incident wave and ps is the total scattered wave, with p inc , p , p s : R 3 R. p ( x ) satisfies the Helmholtz equation,31 
Δ p ( x ) + k 2 p ( x ) = 0 , in R 3 Ω ¯ s R ,
(4)
where Δ = · is the Laplacian operator, with = [ / x 1 , / x 2 , / x 3 ] representing the three-dimensional (3D) Cartesian gradient operator.
FIG. 4.

A two-dimensional schematic of the scattering problem. ps is the total scattered pressure field. Ωs represents a solid enclosed by an unbounded air domain. p inc and ps represent the incident and total scattered pressures, respectively.

FIG. 4.

A two-dimensional schematic of the scattering problem. ps is the total scattered pressure field. Ωs represents a solid enclosed by an unbounded air domain. p inc and ps represent the incident and total scattered pressures, respectively.

Close modal
Here, k = ω / c 0 is the angular wave number, with c0 being the speed of sound in a fluid or air, while Ω ¯ s is the closure of Ωs. For an exterior problem, the scattered field has to satisfy the Sommerfeld radiation condition at r ,
r p s ( x ) i k p s ( x ) = O ( r 1 ) , as r ,
(5)
where r = | x | and i = 1. In addition, the time-harmonic elastic wave in the solid satisfies the elastodynamics wave equation:
· S ( u ( x ) ) + ω 2 u ( x ) = 0 , in Ω s .
(6)
Here, · represents the divergence operator of a tensor-valued function. More specifically, · S may be written as S i j / x j using indicial notation, where the subscripts i and j take on the values 1, 2, 3 and the dummy index j must be summed on. u ( x ) : R 3 R 3 is the displacement vector, and S ( x ) : R 3 R 3 × 3 is the stress tensor. For an isotropic homogeneous solid, the stress tensor can be written in the following form using the generalized Hook's law for finite deformations (see, e.g., Ref. 32, p. 147):
S ( u ( x ) ) = λ tr ( E ( x ) ) I + 2 μ E ( x ) .
(7)
Here, E ( x ) = ( 1 / 2 ) ( u ( x ) + u T ( x ) ) is the infinitesimal strain tensor and μ and λ are Lamé parameters. I and tr represent the identity matrix and the trace function, respectively. Using the expression for the infinitesimal strain tensor and substituting Eq. (7) in Eq. (6) gives (see, for instance, Ref. 33, p. 275)
( λ + 2 μ ) ( · u ( x ) ) μ × × u ( x ) + ω 2 u ( x ) = 0 ,
(8)
where × represents the curl operator. In order to couple the acoustic and elastic wave equations, continuity conditions must be imposed at the interface. Continuity of the velocity dictates
n p ( x ) = ρ a ω 2 n ( x ) · u ( x ) , on Γ s .
(9)
Here, ρa is the density of air and n ( x ) R 3 R 3 is the outward normal directed from the solid to air. Continuity of traction on the surface of the elastic solid yields
p n ( x ) = S ( u ( x ) ) · n ( x ) , on Γ s .
(10)
As such, the acousto-elastic interaction boundary value problem may be formulated as follows:
{ Δ p + k 2 p = 0 , in R 3 \ Ω ¯ s , · S ( u ) + ω 2 u = 0 , in Ω s , n p = ρ a ω 2 n · u , p n = S ( u ) · n , on Γ s = Ω s , r p s i k p s = O ( r 1 ) , as r .
(11)
We assume that the above continuous interaction problem is well-posed and has a unique solution. For proof and discussion of the well-posedness of the problem, the interested reader is referred to Ref. 34.
In the case of a rigid scatterer, the Helmholtz equation is considered by itself and the boundary conditions are reduced as follows:
{ Δ p + k 2 p = 0 , in R 3 \ Ω ¯ s , n p = 0 , on Γ s = Ω s , r p s i k p s = O ( r 1 ) , as r .
(12)
HRTFs are formally defined as the complex frequency ratio between the frequency response at the center of the blocked meatus and the free-field response, which commonly corresponds to response at the center of the head with the head absent.35 In other words, for an omnidirectional point source located at ( r , θ , ϕ ) relative to the head center in a vertical polar spherical coordinate system,36 an HRTF at a frequency f = ω / 2 π is defined as
H ( r , θ , ϕ , f ) = P ( r , θ , ϕ , f ) P 0 ( r , f ) .
(13)
Here, r is the radius, θ [ 0 ° , 360 ° ) is the azimuth angle, with 0 ° in front and 90 ° to the left, and ϕ [ 90 ° , 90 ° ] is the elevation angle, with ϕ < 0 denoting directions below the horizontal plane. P represents the Fourier transform of the sound pressure at the blocked meatus, and P0 represents the Fourier transform of the sound pressure at the head center with the head absent.

We simulate far-field HRTFs for the tabletop Brüel and Kjaer high-frequency HATS type 5128 in two cases: (1) the HATS is modeled as perfectly rigid, and (2) the pinna region is modeled as linear elastic with material properties close to that of auricular cartilage, while the remainder of the head and torso is modeled as perfectly rigid. The former requires the acoustic model from Eq. (12), whereas the latter must be treated as a coupled elasto-acoustic (vibroacoustic) problem, as prescribed by Eq. (11). HRTFs are computed for a plane wave source of unit amplitude at seven spatial directions, five of which sweep the horizontal plane from 0 ° to 180 ° at 45 ° azimuth increments. The last two directions lie in the frontal plane, with an azimuth of 90 ° and a ± 45 ° elevation (see Fig. 5). We do not use reciprocity. Although the reciprocity principle holds for any linear time-invariant system in theory, and it has been shown to have negligible error in HRTF simulations with perfectly rigid boundary conditions,37 we decided against using it for the following two reasons. (1) The first reason is to enable direction-dependent error estimation. This will be more clear in Sec. IV G, where upper bounds for the error are proposed as a function of the global L2 norm of the error. (2) The second reason is to avoid any unforeseen numerical difficulties that might arise from having the source close to the surface of the elastic body in the elasto-acoustic case.

FIG. 5.

(Color online) Schematics of the comsol model showing the Cartesian coordinate system. Simulated HRTF directions are given as (azimuth, elevation) pairs.

FIG. 5.

(Color online) Schematics of the comsol model showing the Cartesian coordinate system. Simulated HRTF directions are given as (azimuth, elevation) pairs.

Close modal

We use the comsol multiphysics software (COMSOL, Inc., Burlington, MA, USA), version 6.0, to discretize the continuous models introduced in Sec. III using the finite element method (FEM). In Secs. IV A–IV E, we will describe the development of models.

A computer-aided design (CAD) drawing, provided by HBK (Chicago, IL), is used to include the geometry of the HATS in the model. The mouth and nose cavities are filled to avoid unwanted resonance frequencies that could adversely affect numerical convergence of the discrete models. In the elasto-acoustic models, separate domains for the HATS pinnae Ωs are created by intersecting a cylindrical object, tightly placed around the ears, with the HATS geometry (see Fig. 6). The HATS is aligned to a right-handed Cartesian coordinate system, as shown in Fig. 5.

FIG. 6.

(Color online) Configuration of the coupled multiphysics model and the elastic pinna domain. The golden domain shows the pinna region modeled using the solid mechanics interface. The blue surfaces show the boundaries where the fixed boundary condition is applied. The inset shows the small spherical domain Ω * (in cyan), which is used for the calculation of HRTFs.

FIG. 6.

(Color online) Configuration of the coupled multiphysics model and the elastic pinna domain. The golden domain shows the pinna region modeled using the solid mechanics interface. The blue surfaces show the boundaries where the fixed boundary condition is applied. The inset shows the small spherical domain Ω * (in cyan), which is used for the calculation of HRTFs.

Close modal

The scattering problems introduced in Eqs. (11) and (12) are defined on an unbounded domain. However, the FEM is conceptually developed for the numerical discretization of bounded domains. Therefore, its use for unbounded domains requires the introduction of an artificial boundary around the scatterer that incorporates the far-field behavior into the finite element model. On this artificial boundary, the finite element discretization may be coupled with some discrete form of the analytic solution.31 Some of the most common coupling approaches include Dirichlet-to-Neumann maps (see Refs. 38 and 39 and references therein), non-reflecting boundary conditions40,41 and perfectly matched layers (PML).42 We choose to truncate the simulation domain using the PML technique.

To truncate the simulation domain, we define a spherical domain, Ωa, around the HATS with radius Ra. Next, we introduce a PML domain by extending Ωa along its exterior with thickness dRPML. Note that the PML's physical thickness does not matter in a frequency-domain simulation. However, the thickness must be chosen such that the PML mesh is more or less regular (Ref. 43, p. 220). It is recommended that the PML's interior boundary be at a minimum distance of λ / 8 (λ represents the wavelength) from the scatterer's surface to avoid unphysical reflections.43 Therefore, for higher simulated frequencies, Ra may be reduced to save computational resources. We exploit this by considering two different model geometries, depending on whether the simulated frequencies are below or above a frequency threshold. This threshold was set to 10 kHz for the acoustic models and 5 kHz for the elasto-acoustic models (Table I).

TABLE I.

Summary of geometry and meshing specifications for the acoustic and elasto-acoustic models used in the study.

Model Ra (cm) dRPML (cm) fmax (kHz)
Acoustic       
f ≤ 10 kHz  15  35  10 
f > 10 kHz  14.5  20  20 
  
Elasto-acoustic       
f ≤ 5 kHz  18.5  35 
f > 5 kHz  15  35  10 
Model Ra (cm) dRPML (cm) fmax (kHz)
Acoustic       
f ≤ 10 kHz  15  35  10 
f > 10 kHz  14.5  20  20 
  
Elasto-acoustic       
f ≤ 5 kHz  18.5  35 
f > 5 kHz  15  35  10 

The speed of sound in air, c0, is set to 343 m/s, and the density of air, ρ0, is set to 1.225 kg/m3. In the elasto-acoustic models, the material composition of the pinna as well as its properties must be chosen. In this work, the main focus is to study the HRTF effects for a human pinna, as opposed to a silicon ear, commonly employed in mannequin fabrication.

The mature human auricle is composed of various tissues and organs, such as skin, subcutaneous tissue, (rudimentary) intrinsic muscles, elastic cartilage, cutaneous nerves, lymphatic vessels, capillaries and arterial branches, sebaceous and sweat glands, and fat (mainly in the ear lobe).44 Except for the lobule, cartilage is the dominant tissue/organ of the pinna.44 Since the mechanical response of the adipose tissues is known to be highly damped,45,46 we will focus on the cartilage tissue and consider it as the most important tissue for the elasto-acoustic simulations.

The material model and its properties for cartilage must be determined. Cartilage is mostly composed of interstitial fluid/water (70%), type II collagen (20%), and proteoglycan aggregates (6%).47 Due to the low permeability of the cartilaginous structure, the interstitial fluid cannot be squeezed out easily and its flow is restrained.47,48 As a result, for short static loading periods or for cyclic loading with moderate to high frequencies, the tissue's mechanical behavior is close to a single-phase solid: under rapidly applied loads, the fluid does not have time to flow relative to the solid cellular matrix.48 Under these circumstances, cartilage can be modeled as a homogeneous, linear elastic and nearly incompressible material.47 As such, cartilage can be characterized by its Young's modulus, Es, Poisson ratio, νs, and density, ρs.

The density of cartilage is slightly higher than that of water.49 We assume a value of ρ s = 1100 kg/m3. We use a complex modulus, E s ( 1 + i η ), common to soft tissue modeling.50 Here, Es is the storage modulus and η is a loss factor that accounts for internal damping. Measurements of the Young's modulus of cartilage have been reported in several works. However, most works focus on articular cartilage or the septum.49,51,52 A notable work by Griffin et al. characterizes the biomechanical properties of the human auricular cartilage and reports the compressional modulus for different parts of the auricle.53 Based on this last reference, a storage modulus of E s = 1.5 MPa was chosen. Measurements for the loss factor and Poisson ratio ν of auricular cartilage are scarce in the literature. Therefore, we're bound to reference values reported for other types of cartilage. The loss factor for other similar media, such as articular cartilage, TM, and skin, is reported in the range of 0.01 to 0.13.50,54 Here, we assume a loss factor of 0.1.

The selection of the Poisson ratio requires more care since it is intimately related to the compressional wave speed and to the nature of the elasto-acoustic coupling. For an isotropic linear elastic solid, ultrasound wave speed may be calculated as c p = ( λ + 2 μ ) / ρ s.49,55 This may be rewritten as c p = E s ( 1 v ) / ρ s ( 1 + v ) ( 1 2 v ). Ultrasound wave speeds in cartilage have been consistently reported in the range of 1500 to 1700 m/s in the literature.49,56,57 Measured values of cartilage's Poisson ratio range from 0.26 to 0.45.49,52,58 These values of Poisson ratio correspond to a maximum compressional wave speed of 70 m/s, which is an order of magnitude lower than the values measured and reported in the literature. The authors believe that this apparent discrepancy is due to the simplifying assumption of modeling cartilage as a homogeneous, linear elastic single-phase solid. One, then, needs to choose whether the model remains faithful to the measured values of Poisson ratio or those of the compressional wave speed. Since, as discussed before, the compressional wave speed plays a fundamental role in determining the nature of the elasto-acoustic coupling, we choose the latter and select a Poisson ratio of ν = 0.4999, leading to a compressional wave speed of 1500 m/s. This is very close to the sound speed in water and is consistent with the high water content of elastic cartilage.47 

The acoustic models are developed using the Pressure Acoustics interface. The Sound Hard Boundary feature is used to model all boundary surfaces as perfectly rigid. For the elasto-acoustic models, the Pressure Acoustics interface is used within the air domain Ωa and the Solid Mechanics interface is used within the elastic pinna domains Ωs (golden region in Fig. 6). To model a plane wave source, a Background Pressure Field is applied to the air domain. Continuity conditions [Eqs. (9) and (10)] are applied to the exterior boundary of the elastic domain using comsol's Multiphysics Acoustic-Structure Boundary node. The HATS pinna is bonded to a hard plastic core within the ear.59 Therefore, for the purposes of the current model, fixed boundary conditions are applied to the interior boundaries of the elastic domains (i.e., all surfaces not in contact with air [see Fig. 6]). This includes the circular surface at the back of the pinna structure, not captured in Fig. 6. Our studies showed that using free versus fixed boundary conditions on the interior boundaries does not make a significant difference in the HRTF response magnitude (results not shown).

A FEM solver with quadratic and cubic Lagrange elements is used for the acoustic and elasto-acoustic models, respectively. The choice for the order of shape functions was based on our findings from code verification exercises (see Sec. IV G).

As was explained in Sec. IV A, two sets of model geometries were used, depending on the simulated frequency range to reduce computational demand. For the acoustic HATS, one set of model geometry was developed for frequencies between 250 Hz and 10 000 Hz and another for frequencies above 10 000 Hz. For the elasto-acoustic HATS, a unique geometry was developed for simulated frequencies in the range of 250 Hz to 5000 Hz and another for frequencies above 5000 Hz.

The coarsest mesh for Ωa is generated such that the maximum element size, h max, is set to λ min / 5, where λ min = c 0 / f max and f max is the highest simulated frequency, except for the elasto-acoustic models with simulated frequencies above 5 kHz, where f max is set to 10 kHz instead of 20 kHz due to the limitation of computational resources (see Table I). The minimum element size is selected such that h min = h max / 2.

We rely on the bulk compressional and shear wave speeds to get a reasonable discretization for the elastic domains. The coarsest mesh is generated such that the maximum element size, h max, is set to λ min s, where λ min s = min { c p / f max , c s / f max }, with c p 1508 m/s and c s 21 m/s. Here, c s = μ / ρ s is the shear wave speed in the solid. The choice of h max, in this case, is partially driven by the availability of computational resources. Since the maximum element size for the coarsest mesh is equal to the minimum wavelength, corresponding shear waves may not be fully resolved on this mesh. However, conducting consecutive convergence studies will provide some further assurance and error estimation for the numerical solutions. Similar to the meshing strategy in the air domain, h min = h max / 2.

For the PML to be well-behaved, the mesh inside the PML domain needs to be structured. Therefore, a swept mesh (with 32 layers in the radial direction) is used inside the PML domain (see Ref. 43, p. 209).

Consecutively refined meshes are generated using a 1.1 refinement ratio, following h min i + 1 / h min i = 1.1 and h max i + 1 / h max i = 1.1, with i = 0 corresponding to the finest grid.

For the acoustic models of the HATS, a sequence of nine successively refined grids were used to conduct convergence studies below 10 kHz. Above 10 kHz, six grids were used. For the elasto-acoustic models, six and five grids were used for simulated frequencies below and above 5 kHz, respectively. The number of grids used was informed by the availability of computational resources and convergence behavior of the models.

To ensure consistency among grids, each mesh was inspected for two built-in quality metrics in comsol, namely skewness and volume versus circumradius. The results show that the refinement is consistent across grids (see  Appendix A).

The simulations are done in the frequency domain. For both acoustic and elasto-acoustic models, the simulated frequencies start at 250 Hz and are linearly spaced at 250 Hz up to 10 kHz. From 10 to 20 kHz, the spacing is increased to 500 Hz to reduce computational demand. We heavily relied on the cluster computing functionalities within comsol to parallelize jobs and speed up run times.

Unlike the usual determination of HRTFs, which involves point-wise computations of pressure at the blocked meatus (see, e.g., Ref. 60) and head center, we compute P in Eq. (13) as the volumetric average of pressure inside a small spherical domain, Ω *, located close to the center of the blocked ear canal (see Fig. 6). This is to accommodate the calculation of error bands for the solution (see the section on “Determination of solution error bands” in the supplementary material). As such, a spherical object with a radius of 4 mm (close to the average radius of the human eardrum61–63) is created at an ∼5 mm distance from the center of each blocked ear canal (Fig. 6). This spatial average will low-pass the frequency response, P. In order to recover the original HRTF, we apply an identical averaging operation on P0 prior to the division in Eq. (13). The average for P0 is calculated analytically by integrating the free-field acoustic response of a unit amplitude plane wave inside a 4 mm-radius sphere located at the head center (for more details, see the section titled “Analytic determination of average free-field pressure” in the supplementary material). The downside of doing this average is that the result will be closer to an HRTF at the center of Ω * (i.e., 5 mm away from the ear canal), which differs from a blocked-meatus HRTF at very high frequencies.

Recall that the objective is to compare numerical solutions from the acoustic and elasto-acoustic simulations to determine the effect of pinna vibrations on simulated HRTFs. Simulation results inevitably contain numerical errors. If the difference between acoustic and elasto-acoustic simulation results is comparable to numerical errors, no conclusion can be drawn about pinna vibrations' effects on HRTFs. However, if numerical errors are shown to be negligible, it suggests that modeling pinna vibrations is responsible for the observed differences. Therefore, numerical errors must be estimated—a process referred to as solution verification. To this end, we used a modified definition of the grid convergence index (GCI), introduced in Ref. 64 to propose a non-sharp upper-bound estimate for the numerical error as follows:
| E ¯ 0 * | dB < 20 ln ( 10 ) × | HRTF | GCI V * .
(14)

Here, HRTF is the simulated HRTF and V * is the volume of Ω * (see Sec. IV F). GCI is determined using convergence studies and global error norms. For more details on how the GCI is calculated, see the “Determination of solution error bands” section in the supplementary material.

In this section, the simulated HRTFs and corresponding error estimates are presented on the finest grid. Error estimates were calculated through convergence studies (see Sec. IV G and references to the supplementary material).

In this section, we will present and discuss the results from the acoustic and elasto-acoustic numerical simulations. Far-field (i.e., r ) HRTFs for the HATS will be compared for seven different directions. Figures 7 and 8 show the comparison for the directions in the horizontal plane, and the frontal plane, respectively. The top panels in each subfigure show the overlaid left/right HRTF plots for both acoustic and elasto-acoustic cases. The shaded regions show the upper bound estimates for the discretization error [Eq. (14)]. The bottom panels show the difference between HRTFs and its corresponding error band. Assuming independence of errors and normal distributions (e.g., Ref. 65, p. 144), the latter is computed as
| E ¯ 0 * | dB , Δ = | E ¯ 0 * | dB , ac · 2 + | E ¯ 0 * | dB , el · ac · 2 ,
(15)
where | E ¯ 0 * | dB , ac . and | E ¯ 0 * | dB , el . ac . are calculated according to Eq. (14) for the acoustics and elasto-acoustics problems, respectively.
FIG. 7.

(Color online) Comparison of far-field HRTFs averaged close to the blocked meatus (see Fig. 6) for the HATS in the horizontal plane using the acoustic (AC) and elasto-acoustic (El.-Ac.) models. The shaded areas show numerical error estimates [see Eqs. (14) and (15)]. The arrows in the insets show a schematic for the corresponding HRTF direction.

FIG. 7.

(Color online) Comparison of far-field HRTFs averaged close to the blocked meatus (see Fig. 6) for the HATS in the horizontal plane using the acoustic (AC) and elasto-acoustic (El.-Ac.) models. The shaded areas show numerical error estimates [see Eqs. (14) and (15)]. The arrows in the insets show a schematic for the corresponding HRTF direction.

Close modal
FIG. 8.

(Color online) Same as Fig. 7, but for directions in the frontal plane.

FIG. 8.

(Color online) Same as Fig. 7, but for directions in the frontal plane.

Close modal

A few observations about the solution uncertainty bands are made from the figures: (1) The uncertainty bands for the solutions grow as the frequency increases, (2) the bands show significant widening at the location of sharp peaks and deep notches, and (3) the error bands for the solution at the contralateral side are in general larger compared to the ones on the ipsilateral side. Such behavior of the discretization error is similar to previous estimates from other numerical methods.6,66 Figs. 7 and 8 show that the difference between HRTFs modeled using the acoustic and elasto-acoustic model is less than 1 dB across the whole frequency range for all directions considered. However, it must also be noted that the error band for this difference is large (especially for frequencies above 10 kHz). Larger error bands generally imply that the value of the mean is less reliable. This is where using a sharper error estimator than the one presented in Eq. (14) may help reduce the uncertainty in the solutions and the uncertainty of the difference as a result. For now, the value of the mean alone suggests that provided the model assumptions in this work, the vibration patterns of the pinna do not have a significant effect on HRTFs for HATS.

As mentioned in Sec. IV B, measured values of cartilage's Poisson ratio range from 0.26 to 0.45 in the literature.49,52,58 Many common pinna simulators, including the HATS 5128 ear simulator (type 4620), are made of silicon rubber.59 Silicon usually has a Poisson ratio between 0.48 to 0.495.67 To account for the difference between the value of Poisson ratio used in this work and the range listed for cartilage and silicon materials, we investigate how the selected value of the Poisson ratio will affect the results presented in Sec. V A. To this end, the elasto-acoustic models were used to run simulations using a Poisson ratio of 0.35 for the pinna: approximately the average value reported in the literature for cartilage. Simulations were only run on the finest grid. Therefore, convergence data and, as a result, error estimates were not calculated for the HRTFs obtained by these models.

Figure 9 shows the difference curves from Sec. V A overlaid by the new ones, obtained by using a Poisson ratio of 0.35, for 0 ° azimuth and elevation. For the former, the uncertainty bands have also been included in the figure. The figure shows that using a lower Poisson ratio leads to an overall larger difference (maximum of 2 dB) between the HRTF solutions. The observation holds for all the HRTF directions studied in this work (see  Appendix B).

FIG. 9.

(Color online) The effect of Poisson ratio ν on the magnitude difference between HRTFs at r obtained with the acoustic model and the elasto-acoustics model at zero azimuth and elevation. Here, ν 1 = 0.4999 refers to the incompressible case and ν 2 = 0.35. The shaded area represents the error estimate from Fig. 7 and Eqs. (14) and (15) available only for ν 1 = 0.4999.

FIG. 9.

(Color online) The effect of Poisson ratio ν on the magnitude difference between HRTFs at r obtained with the acoustic model and the elasto-acoustics model at zero azimuth and elevation. Here, ν 1 = 0.4999 refers to the incompressible case and ν 2 = 0.35. The shaded area represents the error estimate from Fig. 7 and Eqs. (14) and (15) available only for ν 1 = 0.4999.

Close modal

Near-field HRTFs (nf-HRTFs) have a significant perceptual68,69 importance for virtual reality applications, where the user is often interacting with objects in the acoustic near field.70 In this section, we show how using an elasto-acoustic model affects nf-HRTFs using a case study.

We use a monopole source located on the left side and at a 5.5 cm distance from the center of the blocked meatus. The source is located close to the pinna's “focus point”; this was done so that the spherical wave fronts hit the pinna with greater temporal coherence. The intensity of the monopole may be arbitrary for linear models, since the response is normalized with the free-field response at the head center. The original acoustic and elasto-acoustic models (Poisson ratio ν = 0.4999) are used to simulate nf-HRTFs for the tabletop HATS. We note that the models are run only for the finest grid. As a result, no convergence or uncertainty analysis is done for this case. Figure 10 shows that difference between the solutions is at the same order as the difference between far-field HRTFs.

FIG. 10.

(Color online) Comparison of nf-HRTFs within Ω * for a monopole source located at 5.5 cm away from the left ear's center of blocked ear canal (center of the circles in the inset). AC, acoustics model; El.-Ac., elasto-acoustics model. ν = 0.4999 represents computations on the finest grid.

FIG. 10.

(Color online) Comparison of nf-HRTFs within Ω * for a monopole source located at 5.5 cm away from the left ear's center of blocked ear canal (center of the circles in the inset). AC, acoustics model; El.-Ac., elasto-acoustics model. ν = 0.4999 represents computations on the finest grid.

Close modal

The primary goal of this paper is to study whether assuming a rigid boundary condition for the pinna has a significant impact on the simulated HRTFs. An experimental system was used to provide empirical evidence for sound-induced vibrations of an artificial pinna when hit by incident sound waves. Previous notable works investigating the occlusion effect have conducted vibrational measurements on the human ear tissue.71,72 However, to the best of our knowledge, full-field measurements of sound-induced vibrational patterns of pinna replicas do not exist in prior works. With this empirical evidence at hand, we set out to explore how these sound-induced vibrations affect HRTFs by using numerical simulations. Overall, our numerical results show with high confidence that modeling the vibrations of the pinna (assuming a linear elastic, homogeneous, and isotropic constitutive model for cartilage as its main constituent) has negligible effects for most practical purposes up to about 5 kHz: HRTF magnitude differences are less than 1 dB when the results of rigid vs non-rigid boundaries are compared. This agrees well with previous findings in the literature that have shown the shape of the pinna has minimal influence on HRTFs for frequencies below about 3 kHz.73,74 At higher frequencies and for the incompressible case ( ν = 0.4999), results point towards a similar conclusion, but we have less confidence in our numerical results due to a steep increase in the upper-bound estimator for the numerical errors. It is worth mentioning that HRTF measurements at high frequencies suffer from a similar deficiency. Nevertheless, given that the error bound from Eq. (14) is not sharp, it is likely that the differences are similarly small up to higher frequencies. Results hold for both the far field and the near field. Our experimental results (see Sec. II) clearly show that the acoustic field does couple mechanically with the pinna, causing it to vibrate. Such coupling, however, appears to have negligible effects on the simulated HRTF magnitude when the pinna is modeled as a linear elastic, homogeneous, and isotropic medium. In addition, the amplitude of 80 dB SPL (see Sec. II) used during measurements might have triggered nonlinearities. Such nonlinear effects have been previously reported when measuring HRTFs of mannequins with viscoelastic pinnae at increasing SPLs.75 

On the other hand, Fig. 11 shows how variations of the pinna's mechanical properties (in this case, Poisson's ratio, ν) could affect the HRTF magnitudes more significantly above 5 kHz. This suggests that at higher frequencies, the actual assumptions of the material model become more important and an impedance of local reaction might not be a good approximation of the boundary. Based on measured compressional wave speeds in cartilage,49,56,57 and due to the rich blood supply of the pinna,44 the high water content of elastic cartilage, and its low permeability for displacements at audible frequencies,47,48 we believe the pinna is better modeled as nearly incompressible ( ν = 0.4999). Previous wave-based HRTF simulations have already employed such boundary conditions (e.g., modeling the pinna with properties similar to that of water).4,76

FIG. 11.

(Color online) The effect of the Poisson ratio on difference curves for (a) k 2, (b) k 3, (c) k 4, (d) k 5, (e) k 6, and (f) k 7. Here, ν 1 = 0.4999 refers to the incompressible case and ν 2 = 0.35.

FIG. 11.

(Color online) The effect of the Poisson ratio on difference curves for (a) k 2, (b) k 3, (c) k 4, (d) k 5, (e) k 6, and (f) k 7. Here, ν 1 = 0.4999 refers to the incompressible case and ν 2 = 0.35.

Close modal

Accepting a nearly incompressible material model for the pinna has a couple of implications. First, such a model justifies why modeling the pinna as perfectly rigid is a practically valid approach to HRTF simulations, especially for living subjects. Extrapolation from the present results would indicate that a rigid boundary condition is valid at audible frequencies for any fluid-rich tissue such as human skin. A higher Poisson's ratio leads to increased compressional wave speed and higher equivalent acoustic impedance at the boundary surface, resulting in boundary conditions that more closely resemble the perfectly rigid case. Therefore, the overall differences in HRTF magnitudes are expected to be larger for materials with a lower Poisson's ratio. Poisson's ratio was found an important factor in simulated results in other linear elasto-acoustics numerical studies (e.g., it was shown to have the largest effect on numerical predictions of the occlusion effect among all other parameters considered).77 

Existing dummy heads are intentionally made flexible enough to accommodate the fitting of hearing aids,78 phones, or earphones.79 Previous HRTF (see Ref. 80, Fig. 59) and vibro-acoustical81,82 measurements suggest changes in acoustical responses with different shore hardness of the pinna.83 Additionally, our findings suggest that Poisson's ratio could affect the HRTFs. Given that the Poisson's ratio for certain types of silicon may be different than that of a real human pinna, HRTFs measured on head or ear replicas (used in, e.g., ear prosthesis84) with silicon artificial ears may not replicate the HRTFs measured in vivo on human heads. Consequently, the design of digital twins and accurate validation exercises for HRTF numerical simulations against dummy heads with silicon/rubber ears (as opposed to, e.g., mannequins with rigid pinnae) is more involved: the model validation becomes more complex since the elasto-acoustical coupling becomes harder to neglect. As such, HRTF simulations employing passive and locally reacting normal impedance measurements of silicon/rubber materials (e.g., Ref. 85) may be problematic for accurate validation studies. To have more certainty, more rigorous factorial experiments are required to quantify the influence of various material parameters.

HRTFs are most commonly employed in (dynamic) binaural synthesis, where the end goal is to employ HRTF filters to reproduce the pressure on the eardrum of an individual as closely as possible as in a real-life scenario.1,35 HRTFs acquired at the blocked meatus, as in the present study, can achieve this goal only to a limited degree due to additional approximations and assumptions that are rarely met in practice (compensation of the delivery method such as individual headphone equalization and the free-air-equivalence coupling,86 the assumption that the ear canal acts as a simple transmission line,35,87 the use of non-individual HRTFs,88 interpolation errors in dynamic scenarios,89,90 head-over-torso orientations,91 etc.). Currently, the above-mentioned approximations are generally way more difficult to account for or control (see, e.g., Refs. 86, 92, and 93) compared to error magnitudes discussed in the present work. As such, in the broader context of binaural synthesis, the present results suggest that in modeling, the detailed vibrational effects of the pinna are negligible and the error control should first focus on other aspects of the binaural synthesis chain.

The incompressible pinna material model is congruent with ultrasound compressional wave speeds in elastic cartilage (see Sec. IV B). Nevertheless, the anatomy and structure of the human pinna are way more complex than those of a simple cartilage model.44 As such, to have higher confidence in the extrapolation of the present results to in vivo HRTFs, highly controlled HRTF validation studies are needed (e.g., Ref. 6), beyond the preliminary measurements presented in this work in Sec. II. Alternatively, a more complex in silico model would also be useful in assessing the extrapolation power of the present results.

Finally, we would like to emphasize the limitations of the present study. First of all, the composition of the pinna is only a first order physics-based approximation of the highly complicated anatomy of a real in vivo pinna. (1) We consider that the entire pinna is solely composed of a single “irregular plate”44 made of elastic cartilage. (2) We model the ear as a linear elastic, homogeneous, and isotropic material, which is not entirely realistic (e.g., Ref. 48), and (3) we do not consider tissues connecting the pinna to the skull, such as muscles or ligaments: one could envision using an average ear tissue to capture the additional components. However, we believe that homogenization of a complex tissue structure, such as the ear, over a wide frequency range will require a careful study on its own. (4) We ignore potential thermoviscous losses at the surface of the skin and vellus hair. (5) We consider a frequency-independent loss factor, η, for the material. Additionally, in the absence of more reliable data, we have assumed that the loss factor value is in the range commonly referenced for the TM and skin. Dynamic mechanical analysis (DMA) on auricular cartilage samples could better inform the value of the complex moduli used in the model. (6) The effect of temperature on the mechanical behavior of cartilage as well as its other material properties (such as Poisson ratio or loss factor) has not been taken into consideration in this work. Temperature may affect the material's behavior (i.e., constitutive law) and its properties. The storage modulus used for cartilage in this work is based on measurements of the Young's modulus at 37.5 °C,53 which is close to normal body temperatures. (7) The effect of other factors, such as age or smoking, on the mechanical properties of cartilage were ignored. (8) The study did not consider vibrations of the head surface and ignored the impact of the interface between the cylindrical pinna holder and the rigid part of the HATS.

Another limitation of the study is the calculation of HRTFs as a volumetric average within Ω *, an average that correlates to the center of Ω *, 5 mm away from the more common blocked-meatus location (see Sec. IV F). Last, results are presented only for one otologically normal pinna: it is unclear if there are other ear shapes (having, e.g., larger surface areas) that are more sensitive to elasto-acoustic coupling even for ν = 0.4999. Nevertheless, the anatomy of the ear44 and the known mechanics of elastic cartilage suggest a high surface acoustic impedance of the pinna. Thus, the transmission of acoustical energy for this value of Poisson's ratio should be minimal across the ear population.

See the supplementary material for details on the solution verification exercises, including convergence studies, and error estimation for the numerical solutions.

The authors have no conflicts to disclose.

The data that support the findings of this study are available from the corresponding author upon reasonable request.

Acoustic models: Below 10 kHz, the average and standard deviation of the skewness quality metric are 0.684 × 10−3 and 5.6 × 10−3, respectively. The average and standard deviation for the volume versus circumradius quality factor are 0.733 × 10−2 and 1.21 × 10−2. A similar study was performed on the mesh sequence above 10 kHz. The average and standard deviation for the skewness (volume to circumradius) metric were 0.686 (0.735) × 10−3 and 3.1 (7.08) × 10−3, respectively. The reported metrics show that the refinement is consistent across grids (Table II).

TABLE II.

Skewness quality metric Q1 of the mesh sequences used for convergence studies on the acoustic HATS model.a

i f 10 kHz f  > 10 kHz
Ne Q ¯ 1 ( σ Q 1 ) Ne Q ¯ 1 ( σ Q 1 )
0.57  0.682 (0.110)     
0.74  0.683 (0.109)     
0.98  0.684 (0.109)     
1.29  0.684 (0.108)  3.63  0.686 (0.107) 
1.71  0.685 (0.108)  4.83  0.686 (0.106) 
2.27  0.685 (0.108)  6.43  0.687 (0.106) 
3.00  0.685 (0.107)  8.56  0.689 (0.118) 
3.99  0.686 (0.107)  11.39  0.687 (0.105) 
5.30  0.686 (0.107)  15.19  0.680 (0.126) 
i f 10 kHz f  > 10 kHz
Ne Q ¯ 1 ( σ Q 1 ) Ne Q ¯ 1 ( σ Q 1 )
0.57  0.682 (0.110)     
0.74  0.683 (0.109)     
0.98  0.684 (0.109)     
1.29  0.684 (0.108)  3.63  0.686 (0.107) 
1.71  0.685 (0.108)  4.83  0.686 (0.106) 
2.27  0.685 (0.108)  6.43  0.687 (0.106) 
3.00  0.685 (0.107)  8.56  0.689 (0.118) 
3.99  0.686 (0.107)  11.39  0.687 (0.105) 
5.30  0.686 (0.107)  15.19  0.680 (0.126) 
a

i indexes the grid, while Ne represents the number of elements (not including the PML) in millions rounded to 2 decimal places. Values are given as average (standard deviation).

Elasto-acoustic models: Below 5 kHz, the average and standard deviation of the skewness quality metric are 0.670 × 10−3 and 3.10 × 10−3. The average and standard deviation for the volume versus circumradius quality factor are 0.701 × 10−3 and 6.64 × 10−3. A similar study was performed on the mesh sequence above 5 kHz: the averages (standard deviations) for the skewness (volume to circumradius) metric were 0.676 (0.712) × 10−3 and 3.43 (7.28) × 10−3, respectively. The reported metrics show that the refinement is consistent across grids (Table III). The frequencies at which simulations are computed are the same as that specified for the acoustic problem.

TABLE III.

Same as Table II but for the elasto-acoustic model.a

i f 5 kHz f > 5 kHz
Ne Q ¯ 1 ( σ Q 1 ) Ne Q ¯ 1 ( σ Q 1 )
0.30  0.665 (0.117)     
0.36  0.668 (0.116)  0.85  0.670 (0.118) 
0.44  0.671 (0.115)  1.09  0.676 (0.116) 
0.55  0.672 (0.114)  1.42  0.678 (0.114) 
0.70  0.673 (0.114)  1.85  0.679 (0.114) 
0.89  0.670 (0.118)  2.42  0.676 (0.116) 
i f 5 kHz f > 5 kHz
Ne Q ¯ 1 ( σ Q 1 ) Ne Q ¯ 1 ( σ Q 1 )
0.30  0.665 (0.117)     
0.36  0.668 (0.116)  0.85  0.670 (0.118) 
0.44  0.671 (0.115)  1.09  0.676 (0.116) 
0.55  0.672 (0.114)  1.42  0.678 (0.114) 
0.70  0.673 (0.114)  1.85  0.679 (0.114) 
0.89  0.670 (0.118)  2.42  0.676 (0.116) 
a

i indexes the grid, while Ne represents the number of elements (not including the PML) in millions rounded to 2 decimal places. Values are given as average (standard deviation).

Figure 11 shows how the value of Poisson ratio affects the difference curves between the HRTF solutions for the remaining directions (see also Fig. 9).

1.
J.
Blauert
,
Spatial Hearing: The Psychophysics of Human Sound Localization,
revised ed. (
MIT Press
,
Cambridge, MA
,
1997
).
2.
H.
Møller
,
M. F.
Sørensen
,
D.
Hammershøi
, and
C. B.
Jensen
, “
Head-related transfer functions of human subjects
,”
J. Audio Eng. Soc.
43
,
300
321
(
1995
).
3.
M.
Morimoto
and
Y.
Ando
, “
On the simulation of sound localization
,”
J. Acoust. Soc. Jpn. (E)
1
,
167
174
(
1980
).
4.
P.
Mokhtari
,
H.
Takemoto
,
R.
Nishimura
, and
H.
Kato
, “
Comparison of simulated and measured HRTFS: FDTD simulation using MRI head data
,” in
AES 123rd Convention
, New York (
Audio Engineering Society
,
New York
,
2007
), paper no. 7240.
5.
P.
Mokhtari
,
H.
Takemoto
,
R.
Nishimura
, and
H.
Kato
, “
Computer simulation of HRTFS for personalization of 3D audio
,” in
Second International Symposium on Universal Communication
, Osaka, Japan (
IEEE
,
New York
,
2008
), pp.
435
440
.
6.
S. T.
Prepeliţă
,
J. G.
Bolaños
,
V.
Pulkki
,
L.
Savioja
, and
R.
Mehra
, “
Numerical simulations of near-field head-related transfer functions: Magnitude verification and validation with laser spark sources
,”
J. Acoust. Soc. Am.
148
,
153
166
(
2020
).
7.
M.
Vorländer
,
Auralization: Fundamentals of Acoustics, Modelling, Simulation, Algorithms and Acoustic Virtual Reality
(
Springer
,
Berlin
,
2008
).
8.
P. M.
Morse
, “
Some aspects of the theory of room acoustics
,”
J. Acoust. Soc. Am.
11
(
1
),
56
66
(
1939
).
9.
P. M.
Morse
and
R. H.
Bolt
, “
Sound waves in rooms
,”
Rev. Mod. Phys.
16
(
2
),
69
150
(
1944
).
10.
F.
Brinkmann
,
L.
Aspöck
,
D.
Ackermann
,
R.
Opdam
,
M.
Vorländer
, and
S.
Weinzierl
, “
A benchmark for room acoustical simulation. Concept and database
,”
Appl. Acoust.
176
,
107867
107871
(
2021
).
11.
C.-H.
Jeong
, “
Design, simulation, and virtual prototyping of room acoustics: Challenges and opportunities
,” in
24th International Congress on Acoustics
, Gyeongju, South Korea (
International Commission for Acoustics
,
Madrid
,
2022
), pp.
1
8
.
12.
M.
Vorländer
, “
Performance of computer simulations for architectural acoustics
,” in
Proceedings of 20th International Congress on Acoustics
, Sydney, Australia (
International Commission for Acoustics
,
Madrid
,
2010
), pp.
1
6
.
13.
B. F. G.
Katz
, “
Boundary element method calculation of individual head-related transfer function. I. rigid model calculation
,”
J. Acoust. Soc. Am.
110
,
2440
2448
(
2001
).
14.
C.
Jin
,
P.
Guillon
,
N.
Epain
,
R.
Zolfaghari
,
A.
van Schaik
,
A.
Tew
,
C.
Hetherington
, and
J.
Thorpe
, “
Creating the Sydney York Morphological and Acoustic Recordings of Ears database
,”
IEEE Trans. Multimedia
16
(
1
),
37
46
(
2014
).
15.
F.
Brinkmann
,
M.
Dinakaran
,
R.
Pelzer
,
P.
Grosche
,
D.
Voss
, and
S.
Weinzierl
, “
A cross-evaluated database of measured and simulated HRTFs including 3D head meshes, anthropometric features, and headphone impulse responses
,”
J. Audio Eng. Soc.
67
(
9
),
705
718
(
2019
).
16.
B. F. G.
Katz
, “
Acoustic absorption measurement of human hair and skin within the audible frequency range
,”
J. Acoust. Soc. Am.
108
(
5
),
2238
2242
(
2000
).
17.
P.
Plaskota
, “
Research of acoustical impedance of human skin
,”
Vib. Phys. Syst.
30
(
1
),
2019143
(
2019
).
18.
B. F. G.
Katz
, “
Boundary element method calculation of individual head-related transfer function. II. Impedance effects and comparisons to real measurements
,”
J. Acoust. Soc. Am.
110
,
2449
2455
(
2001
).
19.
F.
Pind
,
C.-H.
Jeong
,
A. P.
Engsig-Karup
,
J. S.
Hesthaven
, and
J.
Strømann-Andersen
, “
Time-domain room acoustic simulations with extended-reacting porous absorbers using the discontinuous Galerkin method
,”
J. Acoust. Soc. Am.
148
(
5
),
2851
2863
(
2020
).
20.
K. U.
Ingard
, “
Locally and nonlocally reacting flexible porous layers; A comparison of acoustical properties
,”
J. Manuf. Sci. Eng.
103
(
3
),
302
313
(
1981
).
21.
S. T.
Prepeliţă
, “
Verification and validation of wave-based simulations of head-related transfer functions
,” Ph.D. thesis,
Aalto University
,
Espoo, Finland
,
2020
.
22.
P.
Razavi
,
H.
Tang
,
J. J.
Rosowski
,
C.
Furlong
, and
J. T.
Cheng
, “
Combined high-speed holographic shape and full-field displacement measurements of tympanic membrane
,”
J. Biomed. Opt.
24
,
031008
(
2018
).
23.
W. L.
Oberkampf
and
C. J.
Roy
,
Verification and Validation in Scientific Computing
(
Cambridge University Press
,
Cambridge
,
2010
).
24.
J. T.
Cheng
,
A. A.
Aarnisalo
,
E.
Harrington
,
M.
del Socorro Hernandez-Montes
,
C.
Furlong
,
S. N.
Merchant
, and
J. J.
Rosowski
, “
Motion of the surface of the human tympanic membrane measured with stroboscopic holography
,”
Hear. Res.
263
(
1-2
),
66
77
(
2010
).
25.
J. J.
Rosowski
,
J. T.
Cheng
,
M. E.
Ravicz
,
N.
Hulli
,
M.
Hernandez-Montes
,
E.
Harrington
, and
C.
Furlong
, “
Computer-assisted time-averaged holograms of the motion of the surface of the mammalian tympanic membrane with sound stimuli of 0.4–25 kHz
,”
Hear. Res.
253
(
1-2
),
83
96
(
2009
).
26.
C.
Furlong
,
J. J.
Rosowski
,
N.
Hulli
, and
M. E.
Ravicz
, “
Preliminary analyses of tympanic-membrane motion from holographic measurements
,”
Strain
45
(
3
),
301
309
(
2009
).
27.
I.
Dobrev
,
C.
Furlong
,
J. T.
Cheng
, and
J. J.
Rosowski
, “
Full-field transient vibrometry of the human tympanic membrane by local phase correlation and high-speed holography
,”
J. Biomed. Opt.
19
(
9
),
096001
(
2014
).
28.
H.
Tang
,
P.
Razavi
,
K.
Pooladvand
,
P.
Psota
,
N.
Maftoon
,
J. J.
Rosowski
,
C.
Furlong
, and
J. T.
Cheng
, “
High-speed holographic shape and full-field displacement measurements of the tympanic membrane in normal and experimentally simulated pathological ears
,”
Appl. Sci.
9
(
14
),
2809
(
2019
).
29.
J. J.
Rosowski
,
I.
Dobrev
,
M.
Khaleghi
,
W.
Lu
,
J. T.
Cheng
,
E.
Harrington
, and
C.
Furlong
, “
Measurements of three-dimensional shape and sound-induced motion of the chinchilla tympanic membrane
,”
Hear. Res.
301
,
44
52
(
2013
).
30.
M.
Khaleghi
,
W.
Lu
,
I.
Dobrev
,
J. T.
Cheng
,
C.
Furlong
, and
J. J.
Rosowski
, “
Digital holographic measurements of shape and three-dimensional sound-induced displacements of tympanic membrane
,”
Opt. Eng.
52
(
10
),
101916
(
2013
).
31.
F.
Ihlenburg
,
Finite Element Analysis of Acoustic Scattering
(
Springer-Verlag
,
New York
,
1998
).
32.
K. D.
Hjelmstad
,
Fundamentals of Structural Mechanics
(
Springer
,
New York
,
1996
).
33.
K. F.
Graff
,
Wave Motion in Elastic Solids
(
Dover
,
Mineola, NY
,
1975
).
34.
C. J.
Luke
and
P.
Martin
, “
Fluid-solid interaction: Acoustic scattering by a smooth elastic obstacle
,”
SIAM J. Appl. Math.
55
,
904
922
(
1995
).
35.
H.
Møller
, “
Fundamentals of binaural technology
,”
Appl. Acoust.
36
(
3
),
171
218
(
1992
).
36.
V.
Algazi
,
R.
Duda
,
D.
Thompson
, and
C.
Avendano
, “
The CIPIC HRTF database
,” in
Proceedings of the2001 IEEE Workshop on the Applications of Signal Processing to Audio and Acoustics
, New Platz, NY (
IEEE
,
New York
), pp.
99
102
.
37.
Y.
Kahana
,
P. A.
Nelson
,
M.
Petyt
, and
S.
Choi
, “
Numerical modelling of the transfer functions of a dummy-head and of the external ear
,” in
Audio Engineering Society 16th International Conference: Spatial Sound Reproduction
, Rovaniemi, Finland (
Audio Engineering Society
,
New York
,
1999
).
38.
K.
Gerdes
and
F.
Ihlenburg
, “
On the pollution effect in FE solutions of the 3D-Helmholtz equation
,”
Comput. Methods Appl. Mech. Eng.
170
(
1-2
),
155
172
(
1999
).
39.
I.
Harari
,
I.
Patlashenko
, and
D.
Givoli
, “
Dirichlet-to-Neumann maps for unbounded wave guides
,”
J. Comput. Phys.
143
,
200
223
(
1998
).
40.
D.
Givoli
and
B.
Neta
, “
High-order non-reflecting boundary scheme for time-dependent waves
,”
J. Comput. Phys.
186
,
24
46
(
2003
).
41.
A.
Bayliss
,
M.
Gunzburger
, and
E.
Turkel
, “
Boundary conditions for the numerical solution of elliptic equations in exterior regions
,”
SIAM J. Appl. Math.
42
,
430
451
(
1982
).
42.
B.
Kaltenbacher
,
M.
Kaltenbacher
, and
I.
Sim
, “
A modified and stable version of a perfectly matched layer technique for the 3-D second order wave equation in time domain with an application to aeroacoustics
,”
J. Comput. Phys.
235
,
407
422
(
2013
).
43.
comsol, 6.0 Acoustics Module User's Guide (
COMSOL AB
,
Burlington, MA
,
2021
).
44.
F. E.
Lucente
, “
Anatomy, histology, and physiology
,” in
Diseases of the External Ear
(
W. B. Saunders
,
Philadelphia
,
1995
), pp.
1
18
.
45.
C.
Blondé-Weinmann
,
P.
Hamery
,
V.
Zimpfer
,
T.
Joubaud
, and
S.
Roth
, “
Effects of impedance mismatch on impulse wave's propagation at the outer ear: A preliminary numerical investigation
,”
Comput. Methods Biomech. Biomed. Engin.
24
,
S147
S148
(
2021
).
46.
C.
Blondé-Weinmann
,
T.
Joubaud
,
V.
Zimpfer
,
P.
Hamery
, and
S.
Roth
, “
Numerical and experimental investigation of the sound transmission delay from a skin vibration to the occluded ear canal
,”
J. Sound Vib.
542
,
117345
(
2023
).
47.
D. R.
Carter
and
M.
Wong
, “
Modelling cartilage mechanobiology
,”
Philos. Trans. R. Soc. London, Ser. B
358
(
1437
),
1461
1471
(
2003
).
48.
J. M.
Mansour
, “
Biomechanics of cartilage
,” in
Kinesiology: The Mechanics and Pathomechanics of Human Movement
, edited by
C. A.
Oatis
(
Lippincott Williams & Wilkins
,
Baltimore, MD
(
2003
), Vol. 2, pp.
66
79
.
49.
H.
Nieminen
, “
Acoustic properties of articular cartilage: Effect of composition, structure and mechanical loading
,” Ph.D. thesis,
University of Kuopio
,
Kuopio, Finland
,
2007
.
50.
D. D.
Greef
,
J.
Aernouts
,
J.
Aerts
,
J. T.
Cheng
,
R.
Horwitz
,
J. J.
Rosowski
, and
J. J.
Dirckx
, “
Viscoelastic properties of the human tympanic membrane studied with stroboscopic holography and finite element modeling
,”
Hear. Res.
312
,
69
80
(
2014
).
51.
M. K.
Brummund
,
F.
Sgard
,
Y.
Petit
, and
F.
Laville
, “
Three-dimensional finite element modeling of the human external ear: Simulation study of the bone conduction occlusion effect
,”
J. Acoust. Soc. Am.
135
,
1433
1444
(
2014
).
52.
W.
Grellmann
,
A.
Berghaus
,
E.-J.
Haberland
,
Y.
Jamali
,
K.
Holweg
,
K.
Reincke
, and
C.
Bierögel1
, “
Determination of strength and deformation behavior of human cartilage for the definition of significant parameters
,”
J. Biomed. Mater. Res. A
78
,
168
174
(
2006
).
53.
M. F.
Griffin
,
Y.
Premakumar
,
A. M.
Seifalian
,
M.
Szarko
, and
P. E.
Butler
, “
Biomechanical characterisation of the human auricular cartilages; implications for tissue engineering
,”
Ann. Biomed. Eng.
44
,
3460
3467
(
2016
).
54.
W.
Kabir
,
C. D.
Bella
,
P. F.
Choong
, and
C. D.
O'Connell
, “
Assessment of native human articular cartilage: A biomechanical protocol
,”
Cartilage
13
(
Suppl. 2
),
427S
437S
(
2021
).
55.
A. P.
Sarvazyan
,
M. W.
Urban
, and
J. F.
Greenleaf
, “
Acoustic waves in medical imaging and diagnostics
,”
Ultrasound Med. Biol.
39
,
1133
1146
(
2013
).
56.
J.
Töyräs
,
M. S.
Laasanen
,
S.
Saarakkala
,
M. J.
Lammi
,
J.
Rieppo
,
J.
Kurkijärvi
,
R.
Lappalainen
, and
J. S.
Jurvelin
, “
Speed of sound in normal and degenerated bovine articular cartilage
,”
Ultrasound Med. Biol.
29
(
3
),
447
454
(
2003
).
57.
S.
Patil
,
Y.
Zheng
,
J.
Wu
, and
J.
Shi
, “
Measurement of depth-dependence and anisotropy of ultrasound speed of bovine articular cartilage in vitro
,”
Ultrasound Med. Biol.
30
(
7
),
953
963
(
2004
).
58.
Y.
Chang
,
N.
Kim
, and
S.
Stenfelt
, “
The development of a whole-head human finite-element model for simulation of the transmission of bone-conducted sound
,”
J. Acoust. Soc. Am.
140
(
3
),
1635
1651
(
2016
).
59.
High-frequency head and torso simulator type 5128, Brüel & Kjaer
,” available at https://www.bksv.com/en/transducers/simulators/head-and-torso/hats-type-5128 (Last viewed March 13, 2024).
60.
H.
Ziegelwanger
,
P.
Majdak
, and
W.
Kreuzer
, “
Numerical calculation of listener-specific head-related transfer functions and sound localization: Microphone model and mesh discretization
,”
J. Acoust. Soc. Am.
138
(
1
),
208
222
(
2015
).
61.
M.
Khaleghi
,
J. T.
Cheng
,
C.
Furlong
, and
J.
Rosowski
, “
In-plane and out-of-plane motions of the human tympanic membrane
,”
J. Acoust. Soc. Am.
139
,
104
117
(
2016
).
62.
J. T.
Cheng
,
M.
Hamade
,
S. N.
Merchant
,
J. J.
Rosowski
,
E.
Harrington
, and
C.
Furlong
, “
Wave motion on the surface of the human tympanic membrane: Holographic measurement and modeling analysis
,”
J. Acoust. Soc. Am.
133
,
918
937
(
2013
).
63.
M.
Hiipakka
,
T.
Korkeakoulu
, and
D.
Tiivistelmä
, “
Measurement apparatus and modelling techniques of ear canal acoustics
,” Master's thesis,
Helsinki University of Technology
,
Helsinki, Finland
,
2008
.
64.
P. J.
Roache
, “
Perspective: A method for uniform reporting of grid refinement studies
,”
J. Fluids Eng.
116
(
3
),
405
413
(
1994
).
65.
J. R.
Taylor
,
An Introduction to Error Analysis: The Study of Uncertainties in Physical Measurements
, 2nd ed. (
University Science Books
,
Sausalito, CA
,
1996
).
66.
J.
Meyer
,
S.
Prepeliţă
,
A.
Khajeh-Saeed
,
M.
Smirnov
, and
P.
Hoffmann
, “
Verification on head-related transfer functions of a snowman model simulated using the finite-difference time-domain method
,”
IEEE/ACM Trans. Audio. Speech Lang. Process.
31
,
2579
2591
(
2023
).
67.
K.
Larson
, “
Can you estimate modulus from durometer hardness for silicone? yes, but only roughly… and you must choose your modulus carefully!
,” [White Paper], available at https://www.dow.com/content/dam/dcc/documents/en-us/tech-art/11/11-37/11-3716-01-durometer-hardness-for-silicones.pdf (Last viewed August 19, 2023).
68.
D. S.
Brungart
and
W. M.
Rabinowitz
, “
Auditory localization of nearby sources. Head-related transfer functions
,”
J. Acoust. Soc. Am.
106
(
3
),
1465
1479
(
1999
).
69.
W. L.
Martens
, “
Psychophysical calibration for controlling the range of a virtual sound source: Multidimensional complexity in spatial auditory display
,” in
7th International Conference on Auditory Display (ICAD)
, Espoo, Finland (
International Community for Auditory Display
,
Arlington, VA
,
2001
).
70.
M.
Marschall
,
J. G.
Bolaños
,
S. T.
Prepeliţă
, and
V.
Pulkki
, “
A database of near-field head-related transfer functions based on measurements with a laser spark source
,”
Appl. Acoust.
203
,
109173
(
2023
).
71.
C.
Blondé-Weinmann
,
T.
Joubaud
,
P.
Hamery
,
V.
Zimpfer
, and
S.
Roth
, “
Involvement of the outer ear's cartilage and soft-tissues in hearing protection limitation during high-level impulse noise exposition
,” in
27th International Congress on Sound and Vibration
, Prague, Czech Republic (
International Institute of Acoustics and Vibration
,
Auburn, AL
,
2021
).
72.
C.
Blonde-Weinmann
,
T.
Joubaud
,
V.
Zimpfer
,
P.
Hamery
,
S.
De Mezzo
, and
S.
Roth
, “
Experimental evaluation of earplug behavior in front of high-level impulse noise using laser doppler vibrometer
,”
J. Acoust. Soc. Am.
154
,
792
800
(
2023
).
73.
V. R.
Algazi
,
R. O.
Duda
,
R.
Duraiswami
,
N. A.
Gumerov
, and
Z.
Tang
, “
Approximating the head-related transfer function using simple geometric models of the head and torso
,”
J. Acoust. Soc. Am.
112
(
5
),
2053
2064
(
2002
).
74.
K.
Iida
,
M.
Yairi
, and
M.
Morimoto
, “
Role of pinna cavities in median plane localization
,” in
Proceedings of the 16th International Congress on Acoustics
, Seattle, WA (
International Commission for Acoustics
,
Madrid
,
1998
), pp.
845
846
.
75.
A.
Jost
and
D. R.
Begault
, “
Observed effects of HRTF measurement signal level
,” in
International Conference on Architectural Acoustics & Sound Reinforcement
, St. Petersburg, Russia (
Audio Engineering Society
,
New York
,
2002
), pp.
1
5
.
76.
P.
Mokhtari
,
Y.
Hirota
, and
D.
Morikawa
, “
PRTF notch enhancement by selective modification of pinna surface reflection coefficient
,” in
Proceedings of the 24th International Congress on Acoustics, Ica 2022
, Gyeongju, Korea (
International Commission for Acoustics
,
Madrid
,
2022
), pp.
1
8
.
77.
M. K.
Brummund
,
F.
Sgard
,
Y.
Petit
, and
F.
Laville
, “
On the influence of the material properties of the external ear on occlusion effect simulations
,”
Can. Acoust.
40
(
3
),
110
111
(
2012
).
78.
IEC 60318-7:2022
, Head and Torso Simulator for the Measurement of Sound Sources Close to the Ear, Standard IEC 60318-7 (
International Electrotechnical Commission
,
Geneva
,
2022
).
79.
P.57:Artificial ears
, Artificial Ears, Standard (
International Telecommunication Union—Telecommunication
,
Geneva
,
2021
).
80.
GRAS
,
KEMAR 45BB & 45BC Manual
[Manual Revision 19] (
GRAS Sound and Vibration
,
Holte, Denmark
,
2020
).
81.
R.
Shimokura
,
T.
Nishimura
, and
H.
Hosoi
, “
Vibrational and acoustical characteristics of ear pinna simulators that differ in hardness
,”
Audiol. Res.
11
(
3
),
327
334
(
2021
).
82.
R.
Shimokura
,
T.
Nishimura
, and
H.
Hosoi
, “
Manipulating the hardness of HATS-mounted ear pinna simulators to reproduce cartilage sound conduction
,”
Appl. Sci.
12
(
24
),
12532
(
2022
).
83.
C.
Blondé-Weinmann
,
P.
Hamery
,
T.
Joubaud
, and
V.
Zimpfer
, “
Effects of the hardness of acoustic test fixtures' ears on the evaluation of earplug's direct transmissions facing high-level impulse noises
,” in
Proceedings of Acoustics Week in Canada 2023 Conference
, Montreal, Canada (
Canadian Acoustical Association
,
Concord, Canada
,
2017
).
84.
G.
Wersényi
,
V.
Scheper
,
S.
Spagnol
,
T.
Eixelberger
, and
T.
Wittenberg
, “
Cost-effective 3D scanning and printing technologies for outer ear reconstruction: Current status
,”
Head Face Med.
19
(
1
),
46
(
2023
).
85.
F.
Di Giusto
,
D.
Sinev
,
K.
Pollack
,
S.
van Ophem
,
E.
Deckers
, and
F.
Make
, “
Analysis of impedance effects on head-related transfer functions of 3D printed pinna and ear canal replicas
,” in
Forum Acusticum 2023
, Turin, Italy (
European Acoustics Association
,
Madrid
,
2023
), pp.
1
8
.
86.
D.
Hammershøi
and
H.
Møller
, “
Binaural technique—basic methods for recording, synthesis, and reproduction
,” in
Communication Acoustics
, edited by
J.
Blauert
(
Springer
,
Berlin
,
2005
), pp.
223
254
.
87.
D.
Hammershøi
and
H.
Møller
, “
Sound transmission to and within the human ear canal
,”
J. Acoust. Soc. Am.
100
(
1
),
408
427
(
1996
).
88.
L.
Picinali
and
B. F.
Katz
, “
System-to-user and user-to-system adaptations in binaural audio
,” in
Sonic Interactions in Virtual Environments: Human-Computer Interaction Series
, edited by
M.
Geronazz
and
S.
Serafin
(
Springer International Publishing
,
Cham, Switzerland
,
2023
), pp.
115
143
.
89.
K.
Hartung
,
J.
Braasch
, and
S. J.
Sterbing
, “
Comparison of different methods for the interpolation of head-related transfer functions
,” in
16th International Conference: Spatial Sound Reproduction
, Rovaniemi, Finland (
Audio Engineering Society
,
New York
,
1999
), pp.
319
329
.
90.
C.
Pörschmann
,
J. M.
Arend
,
D.
Bau
, and
T.
Lübeck
, “
Comparison of spherical harmonics and nearest-neighbor based interpolation of head-related transfer functions
,” in
2020 AES International Conference on Audio for Virtual and Augmented Reality
, online (
Audio Engineering Society
,
New York
,
2020
), pp.
1
10
.
91.
F.
Brinkmann
,
R.
Roden
,
A.
Lindau
, and
S.
Weinzierl
, “
Audibility and interpolation of head-above-torso orientation in binaural technology
,”
IEEE J. Sel. Top. Signal Process.
9
(
5
),
931
942
(
2015
).
92.
D.
Pralong
and
S.
Carlile
, “
The role of individualized headphone calibration for the generation of high fidelity virtual auditory space
,”
J. Acoust. Soc. Am.
100
(
6
),
3785
3793
(
1996
).
93.
F.
Wightman
and
D.
Kistler
, “
Measurement and validation of human HRTFs for use in hearing research
,”
Acta Acust. united Acust.
91
(
3
),
429
439
(
2005
).

Supplementary Material