Numerical simulations of near-field head-related transfer functions

: Despite possessing an increased perceptual signiﬁcance, near-ﬁeld head-related transfer functions (nf-HRTFs) are more difﬁcult to acquire compared to far-ﬁeld head-related transfer functions. If properly validated, numerical simulations could be employed to estimate nf-HRTFs: the present study aims to validate the usage of wave-based simulations in the near-ﬁeld. A thorough validation study is designed where various sources of error are investigated and controlled. The present work proposes the usage of a highly-omnidirectional laser-induced breakdown (LIB) of air as an acoustic point source in nf-HRTF measurements. Despite observed departures from the linear regime of the LIB pressure pulse, the validation results show that asymptotically-estimated solutions to a lossless model (wave-equation and rigid boundaries) agree in magnitude with the LIB-measured nf-HRTF of a rigid head replica approximately within 1–2 dB up to about 17 kHz. Except a decreased reliability in notch estimation, no signiﬁcant shortcoming of the continuous model is found relative to the measurements below 17 kHz. The study also shows the difﬁculty in obtaining accurate surface boundary impedance values for accurate validation studies.


I. INTRODUCTION
A head-related transfer function (HRTF) pair represents the Fourier representation of two impulse responses which characterize the sound transmission from a point source to the ears of the listeners under stationary conditions in the free field.HRTFs represent an encoding of the auditory cues resulted from the interaction of the acoustic field with the morphology of the listener.
Near-field HRTFs (nf-HRTFs) are the set of HRTFs for which the source is placed in the vicinity of the listener and which possess a number of distinctive features [1][2][3] compared to the far-field HRTFs.The far-field is considered to start at distances where HRTFs do not change significantly with distance; such radius is usually taken to be 1 m based on spherical head studies 4 that do not fully extrapolate to a real head. 5he particular nf-HRTFs features offer extra auditory cues which are used in, e.g., distance perception 6 or the perception of personal space which is expected to carry special meaning 7 to the listeners.Given a point source around a fixed head, such particular nf-HRTFs features/cues are generated by acoustical mechanisms specific to the near-field: the increased source proximity emphasizes the headshadowing and the (inverse-square law) attenuation effects which in turn translate into considerable increased interaural level differences, 1,2 while the increased changes in the orientation of the ears relative to the source as the source approaches the head 1,3 causes the so-called acoustical parallax effect, 8 which results in an increased variation in the high-frequency nf-HRTF features with source location.The added HRTF variability in the near-field is expected to cause an additional increase in the discretization error of the simulated HRTFs.Despite such added perceptual value, nf-HRTFs and near-field auditory perception are not as studied as the far-field mainly due to additional experimental difficulties. 1,9,10redicting nf-HRTFs through wave-based simulations represents a convenient solution to studying the auditory localization in the near-field.Although investigations of wavebased nf-HRTFs simulations exist, 3,[11][12][13][14] the validity of such simulations has not yet been established except, to some very limited degree, for a very simple snowman model. 15he present study aims to address this gap and validate wave-based simulations of nf-HRTFs by employing the finite-difference time-domain (FDTD) method and reliable verification and validation (V&V) procedures. 16To reduce the presentation size, only the HRTF magnitude will be considered in the current work.

II. BACKGROUND A. HRTF
For an ideal omnidirectional point source at r 1 relative to the center of the head, the HRTFs are formally defined as the following free-field corrected transfer function: 17 HRTFðr where x represents the angular frequency, P ear ðxÞ represents the Fourier transform of a captured pressure signal at a location of interest in the external ear, while P ref ðxÞ is the Fourier transform of a pressure captured at a reference location without the scattering structure of interest.

B. Near-field acoustical source
Compared to far-field HRTF measurements, the deviations of the acoustical source from an ideal point source could significantly bias the resulting nf-HRTFs.To begin with, the physical size of the source can either acoustically couple with the anatomical structures of the listener 10 or can superimpose multiple formal HRTFs (e.g., a volumetric source such as a balloon pop).Second, the directivity of the source will cause deviations from the expected spherical wavefront that in turn will bias the obtained nf-HRTFs.Third, as with any measurement, the acoustical excitation signal must have a good degree of repeatability in amplitude, phase, and directivity.Moreover, the source must be efficient enough in the audible frequency band to provide the minimal required signal-to-noise ratio (SNR) for the nf-HRTF problem.Finally, HRTFs seem to only address the linear acoustical effects around the listener; as such, the sound source must behave acceptably linear.
Previous nf-HRTFs measurements employed sound sources that only partially satisfied such requirements.Brungart and Rabinowitz 1 designed a special sound source based on a horn driver fitted with a tube.Such an inverse horn approach is common in many acoustical measurements requiring a point-source. 18,19Loudspeakers are also a common choice (either single element 20 or a polyhedron radiator, usually, a dodecahedron 21,22 ).Finally, electric sparks have also been used in various nf-HRTF studies. 23,24Other acoustical point sources were designed for various reasons [25][26][27][28][29] but are not found in the scientific literature related to nf-HRTF measurements.
Considering the present validation study, at least two options are possible regarding the acoustical source: (i) using a source with well-characterized acoustical properties and replicating it in the simulation domain; and (ii) finding a source which is as close as possible to an ideal point source which is straightforward to implement in a simulation.The former option would entail either a separate source validation study (a much more difficult multiphysics problem than the present validation study; see, e.g., the work by Karjalainen et al. 30 ) or a source calibration procedure which would significantly lower the reliability of a nf-HRTF validation study.Moreover, employing a non-point source does not adhere to the formal HRTF definition and will narrow the application domain of the validation results.Thus, the latter alternative is presently more attractive.][33][34] One type of sound source which seems to satisfy most of the requirements for nf-HRTFs measurements is the laser-induced breakdown (LIB) of air. 29,35LIBs are formed in a very narrow time window, 36 quickly become massless, 29 offer higher control with more concentrated plasma volumes compared to electric sparks, 33,37,38 could be less invasive to the acoustic field, 39 and seem to be highly omnidirectional in the audible range. 29Consequently, the LIB sound source was chosen in the present validation study.Nevertheless, nonlinearities, extra scattering, and its coupling with the scatterer are of concern for a rigorous nf-HRTF validation study and need to be addressed.
At each boundary point r b , a resistive boundary condition (BC) of local reaction is employed,

D. Discrete models
The same discrete models as in Prepelit ¸ a et al. 41 are used.Briefly, the standard rectilinear (SRL) FDTD update scheme is used to discretize Eqs. ( 2) and (3).The scheme is run at maximal stable Courant number on a uniform 3D Cartesian grid which tessellates the space into voxels.The boundaries are discretized in a stair-stepped fashion on the same grid with a conservative voxelization algorithm.
The used FDTD solver is implemented on graphics processing units 42 (GPUs) and parallelized further using the Message Passing Interface.This implementation was previously verified. 41he continuous asymptotic solution and its precision at a ¼ 0:05 significance level are estimated based on a 1storder weighted regression model and a convergence study. 41ere, with the risk of decreased precision in the asymptotic predictions, no interpolation of the source/receivers is done.

E. Dissimilarity metric
To compare two magnitudes of the transfer functions H 1 ; H 2 2 C at a certain angular frequency x on the decibel (dB) scale, the ratio-scale measure of dissimilarity is employed 43 where E represents the expected value operator, and U dB;95% ½x represents the total uncertainty on the dB scale from both H 1 and H 2 at a significance level of a ¼ 0:05.For example, if the absolute uncertainties of jH 1 j and jH 2 j are known, then U dB;95% ½x can be calculated based on uncertainty propagation laws (Ref.44, pp.19-52), where U jH 1 j;95% ½x and U jH 2 j;95% ½x represents the 95% absolute uncertainty intervals for the magnitude of H 1 and H 2 at x, respectively.Here, U jHj;95% are given by confidence intervals (CIs) at a ¼ 0:05.If, e.g., H 1 is an HRTF as in Eq. ( 1), then the relative error propagation law can be re-applied to obtain U jH 1 j;95% based on U jH 1;ear j;95% and U jH 1;ref j;95% III. METHODS

A. Laser set-up
The used pulsed laser (CFR 400, Quantel laser, Les Ulis, France) is a Q-switched solid-state (Nd:YAG gain medium) laser, pumped by a flashlamp.The wavelength of the 7 mm-diameter laser beam is 1064 lm, while the electromagnetic (EM) pulse duration is 8 ns with a total energy of 400 mJ.No wavelength separation optical parts were attached to the laser.The laser was controlled remotely via the serial link and worked on internal triggering mode (i.e., both flashlamp and Q-switch signals were triggered internally).To reduce acoustical scattering and increase LIB repeatability, 45 the laser was fitted with a Galilean beam expander (all parts manufactured by Thorlabs, Newton, NJ) having a diameter of 50.8 mm.The expander was fitted with a 90 broadband dielectric elliptical mirror and ended with an adjustable lens tube containing a plano-convex lens of 30 cm focal length.Figure 1 shows a 3D rendering of the lasing set-up.

LIB acoustical characterization
While the general free-field acoustical properties of the generated pressure pulse were previously investigated, 29,37,45 the exact acoustical characteristics are highly dependent on the laser set-up (e.g., laser type and energy, 46,47 used optics 45 ), environment (e.g., gas composition, 48 presence of impurities, 49 thermodynamic properties 49,50 ) and other lasing parameters (e.g., EM pulse duration, 51 strength, laser wavelength).Moreover, a nf-HRTF validation study requires a more thorough acoustical characterization of the source.Finally, the linear acoustics model in Eqs. ( 2)-(3) cannot account for the inherent nonlinearities 52 in the acoustical pulse; the modeling error in the present validation study might be unreasonable.As such, it is important to assess the acoustical characteristics of the resulting pressure pulse, given the set-up at hand.
The detailed results of the acoustical analysis of the resulting LIB pressure wave can be found in the supplemental material. 52The main results are summarized here: (1) A window length of t LIB ¼ 91 ls was found sufficient to include the pressure pulse with 8 ls before the main peak and 83 ls after the main peak.(2) A 180 ls Q-switch time and 3 Hz EM pulse repetition rate were found to yield the smallest acoustic-pulse variability without any LIB misses.Measurements indicate an average pressure-magnitude standard deviation of the pulse around r ¼ 0:4 dB.(3) Repeatability of the pulse improved after some usage, likely due to thermal effects in the lasing system.To account for this, the laser was "pre-heated" by firing a number of sparks (note the exact minimal usage could not be easily determined and some subsequent results might be affected).(4) The resulting pressure pulse is highly omnidirectional with a standard deviation of its magnitude in the audible range r % 0:4 dB.Thus, the average magnitude of the generated LIB pulse has a 4r direction-dependent variability within the audible range of 60.9 dB.(5) The resulting pressure pulse is also highly stable in phase for each direction; measurements show a 2r variability in group delay below about 1 ls.(6) The nonlinearity of the acoustic pulse was assessed based on four markers: (a) amplitude decay, (b) propagation speed, (c) magnitude spectrum, and (d) temporal width of the pulse.Results show that: (a) The decay could not be considered fully consistent with the linear theory in the audible bandwidth even 1 m away from the LIB.Nevertheless, the decay is close to the linear theory with lower frequencies decaying at a slightly lower rate.(b) The pulse travels with significantly increased speed up to around 10 cm from the LIB.Otherwise, the estimated packet speed shows slightly increased speed (100%-100.25%c linear ).Moreover, the estimated instantaneous speed of the pulse peak is converging to the expected linear value around 15-20 cm away from the LIB.(c) The magnitude spectrum is mostly stable in the audible range starting from around 4 cm away from the LIB.(d) In a broadband sense, the pulse significantly elongates up to about 0.5 m away from the LIB.Based on the previous analysis, the following distances from the LIB were chosen for the validation study: • r LIB ¼ 50 mm, the pulse is nonlinear.
• r LIB ¼ 200 mm, the pulse elongation starts decelerating and the instantaneous velocity of the pulse converges to c linear .• r LIB ¼ 500 mm, the pulse elongation becomes "insignificant." • r ¼ 1000 mm, HRTF could be compared to an HRTF measured using a conventional source at 1 m from the center of the head.
Here, r LIB represents the distance from the LIB which, for the measurements/simulations of P ear in Eq. ( 1), will later correspond to the (minimal) distance from the LIB to any surface of the head.Note such distances apply only to the current lasing set-up.

B. MoRa head replica
An accurate in vivo scan of a human head was acquired using an accurate blue-light scanner (Space Spider, Artec 3D, Luxembourg).Since the pinna cavities are difficult to scan due to occlusion, 54 a cast was obtained of each pinna surface which was subsequently scanned (i.e., the surface of the inner cavities of the pinna was scanned as seen from "within" the pinna) and fused with the full head scan.The estimated geometrical accuracy of the resulting mesh is 1 mm.Such mesh was processed and 3D printed to facilitate the subsequent validation study.The resulting head replica, a Modified Ravish's head, will be referred to as the MoRa head.
The scanned mesh was first aligned such that the interaural axis passed through the center of the ear canals and formed the horizontal plane with the tip of the nose.The horizontal plane matched the xy Cartesian plane with z pointing toward the top.The ear canals were manually blocked and the concha slightly altered around the interaural axis such that the local mesh formed a vertical plane.Two aligned cylindrical holes were then cut in the mesh such that the center of the fitted microphones' diaphragms will be on the interaural axis.The design involved surrounding the microphones with vibro-isolating material.The bottom of the neck was cut with a plane parallel to the xy plane.Scanned extra hair not covered by a cap was manually removed from the mesh close to the auricles.
To reduce orientation-related errors for a validation study, well-defined reference one-dimesional (1D)/two-dimensional (2D) manifolds are needed. 55As such, 4 mm-length oval holes of 1 mm height were designed at various points on the head surface to mark the horizontal, median, and frontal reference planes.Crosshair holes were designed to mark the intersection of such planes at the top, front, and back of the head.
The MoRa head was finally split into two parts.For each part, an interior surface was specially designed to facilitate fixation and the mounting of the microphones.The top part was fixed with three screws to the lower part.The lower part was designed with a hole around the back of the neck for the microphone cables.The bottom head-plane had four holes designed to fit four M6 screws for fixation.See Fig. 2.
The designed MoRa mesh was 3D printed in a stiff nylon-based material Polyamide 12 (PA12) having a density q ¼ 1:01 gcm À3 (ASTM D792) and a tensile modulus of 1700 MPa-1800 MPa (ASTM D638), 56 using a multi-jet printer (Jet Fusion 3D 4200A, Hewlett-Packard, Palo Alto, CA).The fabricated MoRa head weighed about 1.9 kg.The reference holes were filled with white room-temperaturevulcanizing (RTV) silicone and a rubber ring of 20A hardness was fitted between the two parts for acoustic sealing.The back-neck hole was sealed with a 30A flex filament.

Printed MoRa head validation
The accuracy of the 3D-printed exterior surface was validated against the original scanned mesh.Here, the FIG. 2. (Color online) Designed MoRa head replica, vertical-polar coordinate system, and validation locations (colored balls).r LIB represents the distance to the surface of the head and not its center r.More details in Table I printed MoRa head was scanned using an identical procedure and apparatus as described in Sec.III B. The resulting mesh was cleaned of artifacts, aligned to the original mesh using a simple iterative closest point (ICP) algorithm, 57 and the two-sided Hausdorff distance 58 was calculated with an open-source tool. 59The Hausdorff error was found to be within 1 mm for the surfaces of interest with the exception of the inferior crus of the antihelix which was within 2 mm.Results are shown in Fig. 3.
The Hausdorff distance will embed different errors: 3Dprinting errors, scanning errors, mesh-alignment errors, software parameters, and expertise of the operator.Although repeated Hausdorff measures were done with similar results, error cancellation is still possible.Nevertheless, due to the complexity of the surface, it is assumed that such a scenario is quite unlikely and the accuracy of the fabricated MoRa head replica is considered satisfactory for the present study.

C. Acoustical impedance assessment of MoRa head material
Due to discrete-boundary convergence issues, 60 the present study focuses on acoustically-rigid BCs.Nevertheless, since impedance data on 3D-printed materials are missing, the surface impedance of the printed material is investigated in a measurement tube to reduce the general uncertainty.See more details about the tube measurements in the Appendix.
The 3D-printed samples were cylinders having a diameter 0.5 mm smaller than the tube diameters of f2:9; 10g cm.Measurement bias was captured by varying different sample parameters.Since the roughness (equivalently, porosity 61 ) of the surface and the mechanical strength 62,63 of the part can change based on deposition angle, each part was printed in two orientations relative to the height of the cylinder: horizontal and vertical.A total of five heights were designed for each build orientation, h 2 f5; 10; 14:5; 25; 50g mm, where 14.5 mm is the average thickness of the MoRa head mesh.Since sample mounting is known to induce measurement bias, 64,65 three different mountings were employed: unhardened polymer clay, unhardened light clay, and duct tape.The clay was applied to the edges of the sample.The top and the bottom of each printed sample was measured once, yielding a total of N m ¼ 2 Â 3 Â 5 measurements per print orientation.
No large differences in the reflection coefficient R were found for build orientation in the analyzed 0.5-6.4kHz range.Regarding the mounting type, the duct tape likely allowed for more acoustical leakage and showed slightly lower R values.The sample height did not seem to correlate with R in a consistent way.
The reflection coefficient R results were pooled across all frequencies (60 measurements/frequency), measurement conditions, and sample sizes.Results are shown in Fig. 4 (top) together with double-tailed studentized CIs for the jRj values (i.e., not standard errors) at a ¼ 0:05 significance level.The results for no sample (i.e., with the hard termination of the tube, N m ¼ 2 Â 10) are also shown.Results indicate that the employed setup contains some bias (i.e., acoustical leakage) and some variability.On average, the PA12 material shows 0.01 lower jRj: average value significantly smaller as found by a Welch's t-test (t ¼ 230.7, p ¼ 0).As such, given the accuracy of the set-up, the BC of the printed material could locally be considered as acoustically rigid for all practical purposes below 6.4 kHz.
Finally, the phase of the samples was analyzed and compared to the rigid termination of the tube.Results in Fig. 4 (middle) show that the reflection coefficient is almost linear-phase, suggesting a purely resistive BC.The phase of the boundary impedance is shown in Fig. 4 (bottom).It is reliably calculated up to about 4 kHz where it shows the typical behavior for a plane wave impinging on a rigid wall; the velocity is about 90 out of phase with the pressure, independent of frequency.Thus, provided the local reaction approximation remains valid, a large purely resistive surface impedance seems like a good model, at least below 4 kHz.

D. HRTF measurements
The measurement setup is depicted in Fig. 5(a).More hardware details can be found in the Appendix.The lasing was configured on the central axis of the end stage of the expander (see Fig. 1) such that the LIB occurred in a vertical reference plane.The LIB could move vertically (i.e., among the Cartesian z axis) in the reference plane which also passed through the horizontal translator axis, the center of the turntable, and the center of the head.The reference plane is perpendicular to the horizontal plane and includes the rotational axis of the turntable.
The measurement bias needs to be minimized since only one measurement set-up is employed; the measurement bias is only weakly captured by the stochastic processes of the LIB.The following angles were measured with a digital inclinometer Bosch GIM 60L of 0.05 claimed accuracy in the horizontal and vertical planes; the lasing optics were within 0.1 of the horizontal plane, the MoRa head replica was horizontally within inclinometer accuracy of 0.05 , while the linear translators were within 0.05 -0.1 of the horizontal and vertical planes, respectively.The exact location of the LIB was manually adjusted in the reference plane with an estimated maximal error of % 2 mm.Some potentially small (on the order of a few millimeters) heightdependent bias was observed during the acoustical assessment of the LIB. 52o explore any potential vibrational leakage from the printed MoRa head to the microphones, each microphone was fitted surrounded by a layer of viscoelastic material.The P ear measurements were repeated for three 3D-printed urethane rubber layers of different shore hardness in the following order: 60A, 20A, 50A.The microphones were positioned such that their grids were flush to the surface surrounding the ear-canal holes.Consistency in placing the microphones between measurements with different hardness was ensured with the help of a specially designed tool.For all cases, the viscoelastic material was not perfectly flush to the microphone grid; a small (i.e., <1 mm) tubular hole was created between the microphone and the walls of the ear canal [see Fig. 5(c)].Finally, some absorptive material was placed inside the MoRa head replica.
For the P ref measurements, the same microphones and signal routing apparatus used during the P ear measurements were employed.Each of the two microphones was placed at a distance from the LIB equal to the radius of each direction in the validation study (i.e., the r value in the first column of Table I).For these measurements, the microphones were placed at 90 incidence at the same height as the LIB; each microphone was mounted on an arm fixed to the horizontal translator (through the rotating table).Only the horizontal translator was used to adjust the distance from the LIB; thus, the measurement direction relative to the LIB for P ref was the same direction as for the linearity measurement points in Fig. 1.
The temperature (T) and relative humidity (RH) sensor HMP110 (Vaisala, Helsinki, Finland) and the laser front panel ICE450 were present in the anechoic chamber at about 1.5-2 m from the spark [see Fig. 5(a)].Since the laser front panel generated some clearly audible noise, the power spectral density (PSD) inside the anechoic chamber was estimated using Welch method 66   IV.VALIDATION

A. Validation measurements and directions
A total of 17 directions were initially measured.For each direction and r LIB 2 f5; 10; 20; 50g cm value, the radius r was manually evaluated such that a ball of r LIB radius does not intersect the MoRa head surface.A distance of r ¼ 1 m was also measured.To limit the analysis, only two directions are selected for the validation study: ð45 ; 0 Þ and ð270 ; 45 Þ. Due to limitations in the set-up, an elevation of / ¼ 40 was used for r ¼ 1 m and ð270 ; 45 Þ direction.Finally, only one r LIB ¼ 1 cm location was acquired at ð125; 277:2 ; 0 Þ(see Fig. 2).A measurement interval of 6 ms was used for the P ear measurement, while 3 ms was used for P ref subsequently zero-padded to 6 ms.The measurement sampling frequency was f s;m ¼ 4 MHz.
To reduce the random errors in the measurements, a total of N LIB ¼ 500 spark measurements were acquired for each term in Eq. ( 1); P ear was measured for each shore hardness totaling 500 Â 3 measurements.Due to memory limitations, the P ear measurements for each direction were split into two ensembles: the translators were homed and then 250 measurements were obtained for each direction; subsequently, the remaining N LIB ¼ 250 measurements from the same .csvfile were acquired without homing the linear translators.Thus, splitting the two measurements into two consecutive ensembles could capture both environmental drifts and drifts/ordering issues in the positioning apparatus.

Post-processing
Measurements for the validation directions were qualitatively analyzed and two conservative observation intervals t obs were chosen: 2.5 and 5 ms.This entailed two frequency resolutions Df of 400 and 200 Hz, respectively.For each P ear signal and direction, sine-squared onset/offset ramps were sampled and applied.For P ref signals, a similar window was applied with a constant onset/offset time of 0.05 ms applied around a rectangular window 52 of t LIB ¼ 91 ms.See Table I for more details.

Microphone viscoelastic treatment
To begin with, the differences between the two ensembles of N LIB ¼ 250 measurements for each direction are analyzed.An error metric similar to the validation error metric in Eq. ( 4) was used to study the differences in the averaged responses-the uncertainty intervals were calculated based on a two-tailed t-distribution based on the standard error of the mean at a ¼ 0:05 significance level.The analysis is conducted based on the post-processed HRTFs for the validation directions (see Table I).For each validation direction, the two measurement ensembles are generally within 0.5 dB of each other, independent of shore hardness of the viscoelastic material.As such, in the subsequent discussion, the two measurement ensembles will be pooled and analyzed as an ensemble of N LIB ¼ 500 measurements.
Considering the large sound pressure level (SPL) values of the LIB acoustic pulse, 52 both vibroacoustic coupling and vibrational leakage are of concern in the audible range for the present validation study.The latter can be investigated by analyzing the changes in the HRTFs with different viscoelastic materials; for instance, changes in the resonant frequencies assuming a single degree of freedom damped oscillator.Again, a metric similar to the validation error metric in Eqs. ( 4)-( 6) was used.
Figure 6 shows the differences in the HRTF magnitudes for the different viscoelastic materials used.It can be seen that the influence of the vibrational isolation material surrounding the microphone is generally small below around 17 kHz.The largest deviations are seen for ð125 mm, 277:2 ; 0 Þ: the acoustic pulse is the furthest away from the linear regime at this location.In addition, most sharp deviations below 17 kHz in Fig. 6 appear at HRTF-notch frequencies.No significant change in the HRTF features was observed below 17 kHz.As such, it is unlikely that vibrational leakage to the microphone or vibroacoustical coupling of the microphone to the field happens below 17 kHz.Consequently, to simplify the analysis and greatly reduce the measurement uncertainty, the P ear measurements for each direction will be pooled and averaged across all the N LIB ¼ 1500P ear measurements.
Although the present analysis excludes some potential measurement errors, it cannot assess the vibroacoustical coupling between the acoustic field and the MoRa head replica, which could cause violations in the local-reaction assumption of the model.

B. Validation simulations
Due to concerns of scattering from the lasing optics, a 3D mesh replica of the measurement system was created.Such a triangulated replica used the original meshes from various manufacturers and contained no cables, connecting parts (e.g., lasing cooling tubes), or the mesh floor grille of the chamber (see Fig. 5).The individual virtual objects were hierarchized in 3ds max V R (Autodesk Inc., San Rafael, CA) and forward kinematics was employed to easily replicate the positioning and orientation of the various parts of the measurement set-up.
The MoRa mesh was first processed to match the measurements; the reference holes used in orientation matching and the ear-canal holes were capped with planar surfaces.To reduce the number of triangles, the inner surfaces of the MoRa mesh were deleted.Subsequently, the MoRa mesh was initially oriented in the reference 3D Cartesian system as shown in Fig. 5.For directions with t obs ¼ 5 ms, the meshes of the set-up were mildly decimated 59 (measured mesh error 0.5 mm).
Since the relative position of the lasing optics to the MoRa head is direction dependent, a simulation domain was created for each direction in Table I that contained a point source and two blocked-meatus point-receivers.For each grid and validation location, the voxel location of each receiver was found by searching for the first air voxel from the blocked-meatus surface in the outward direction on the interaural axis.The P ref and P ear simulations were configured with the expected small-acoustics sound speed c based on the averaged T and RH measured for each validation location (ambient pressure assumed to be 1 atm).
The size of each simulation domain was set as an axisaligned rectangular box such that no domain reflection arrived from the source to the two receivers.The six distances 6d x=y=z from the source to each wall of the domain were estimated for each validation location based on the simulation interval t obs (see Table I), the direct specular reflection, and the corresponding sound speed c.For each location, after positioning the measurement set-up replica mesh, the domain bounding box was created as uniquely specified by d x=y=z from the source.Then, any objects from the set-up outside such boxes were cut using Boolean operations on meshes.An example can be seen in Fig. 5(b).
Each simulation was driven by a discrete delta sequence.All meshes of the set-up had an acoustically rigid BC.For the MoRa head mesh, two BCs were simulated: perfectly rigid (b ¼ 0) and absorbing (b abs , see Sec.IV B 1).For each direction and b value, a convergence study on six grids was conducted (see Table I).The grid refinement ratio between grids was around f s;new =f s;old ¼ 1:1 and the voxel sizes were chosen based on GPU memory limitations for the smallest grid.
For each grid i, f s;i was chosen such that for each location, the frequency resolution Df i is kept constant (see Table I) and such that the same frequency bins are sampled in the discrete Fourier transform.
The P ear and P ref simulations were post-processed with the same sampled windows as the measurements (see Sec. IV A 1).The resulting simulated HRTFs were then used to estimate, without interpolation, the magnitude of the asymptotic solution g HRTF and its uncertainty similarly to Prepelit ¸ a et al. 41 1. b abs Absorbing b abs values are used here to study the effects of the surface impedance concept for an accurate HRTF validation study.Rigorously, an uncertainty analysis is required; since the code only supports frequency-independent resistive BC, the impedance input uncertainty is treated here as epistemic uncertainty.Such uncertainty is usually treated as an interval without any associated probability density function (Ref.16, p. 544).Under the assumption that the HRTF magnitude decreases monotonically with increasing b, a tedious uncertainty analysis can be avoided; estimating the HRTFs for the maximal b abs and for b ¼ 0 should approximately yield the full uncertainty interval in the HRTFs.
The maximal value is chosen based on the impedance measurements (see Sec. III C); the U jRj;95% uncertainty from Fig. 4 (top) is first propagated (Ref.67, p. 52) to b yielding U b;95% .Then, the maximal possible value is retained, V. RESULTS

A. Qualitative comparison
Figure 7 shows the qualitative comparison of simulated and measured nf-HRTF magnitude.Except for the contralateral ear for the ð125; 277:2; 0Þ location, measurements seem to agree well with the rigid-wall simulations.Poor results for the ð125; 277:2; 0Þ location were expected due to the strong nonlinearities in the resulting pressure pulse of the LIB (see Sec. III A 1).Thus, it is surprising that results agree well for the ipsilateral ear for the source at ð125; 277:2; 0Þ.The notch mismatch for such location could be due to differences in sound speeds.
The asymptotic solutions generally show good precision below 15-17 kHz for the grid sizes DX employed.Although their center frequencies were reasonably well predicted, some notches show poor magnitude prediction: their magnitude was around 0 [e.g., the 10 kHz notch for the contralateral ear for (196, 45, 0)].Moreover, the first notch for the ipsilateral ear and ð45 ; 0 Þ validation direction is systematically mispredicted at slightly higher frequencies, suggesting some small bias in the validation process.
Considering the magnitude, the rigid BC predicts the measurements rather well up to 16-17 kHz.At higher frequencies, more absorption is present and the b abs simulations are closer to the measurements.There are three unexpected exceptions where the measured magnitude is larger than the rigid BC predictions: the peaks around 11 kHz for (196, 45, 0) and (298, 45, 0), and the notch around 12.7 kHz for (298, 45, 0).The authors do not have a good explanation for such a result.For instance, considering the 12.7 kHz notch, it does not seem to be a locally-reacting impedance mismatch since increasing the impedance pushes the magnitude even lower-vibro-acoustical coupling of nonlinearities in the resulting pressure pulse could offer a better explanation.
Figure 7 clearly shows that the signal coming from the right-ear microphone was contaminated with noise: although higher frequencies are well-predicted, both contralateral and ipsilateral right-ear measurements show poor SNR below 3 kHz.The exact source of the responsible noise is unknown.Finally, due to the large number of measurements (N LIB ¼ 1500), the uncertainty in the averaged measured nf-HRTFs is negligible.

B. Quantitative comparison
Figure 8 shows the validation results for the metric in Eqs. ( 4)-( 6); here, H 1 is the asymptotic prediction, while H 2 is the measured HRTF.The CIs for the simulated magnitudes are given by the 95% BC a pair-bootstrapped intervals, while the standard error and two-sided student t distribution are used to express the uncertainties in the measured P ear and P ref in Eq. ( 6).
Due to the poor results, the result for the contralateral ear for the ð125; 277:2; 0Þ location is excluded from the analysis.
First, the large deviations below 5 kHz seen for some locations in Fig. 8 are due to the increased measurement noise in the right microphone signal (see Fig. 7) and are disregarded in the following analysis.
The metric in Fig. 8 is strongly affected by the deeper notches of the simulation.The simulated notches are expected to have lower prediction precision and an increased SNR compared to measurements. 41Moreover, measured notches are also quite sensitive to measurement noise (see, e.g., Sec.IV A 2).As such, if such notches are not considered, the asymptotic solutions of the employed wave equation agree with the measurements within 1 dB up to around 11-12 kHz.Moreover, the agreement is within 2 dB up to 17 kHz, with small exceptions.
Above 17 kHz, despite the generally good qualitative prediction (see Fig. 7), the rigid-wall simulated HRTFs seem to be a poor predictor.Considering the results for b abs in Fig. 7, the model needs to incorporate additional losses at such frequencies.Also note the viscoelastic treatment around the microphone also affects the results at such frequencies (see Sec. IV A 2) .
Considering the two ears, the ipsilateral ear shows a slightly improved match with measurements: qualitatively, the oscillations around the ideal 0 dB error appear smaller in Fig. 8. Quantitatively, the common spectral distortion (SD) metric can be used.Excluding the results for ð125; 277:2; 0Þ at the contralateral ear, data shows an average (across frequencies and direction) absolute error on the dB scale within 1-17 kHz of 1.85 dB, an average of 1.49 dB for results at the ipsilateral ear, while 2.27 dB for the contralateral ear.The exact reason for a better match for the ipsilateral ear is unclear: both the simulation and the measurement uncertainties show similar magnitude at both contralateral and ipsilateral ears.Potential unacknowledged biases could be responsible: they could affect the contralateral ear more due to, e.g., increased sensitivity in HRTF magnitude. 68

VI. DISCUSSION
With few exceptions, the predicted asymptotic magnitudes of the employed lossless model agree within 1-2 dB  4)] and 95% CIs [Eqs.(5) and (6)] between the estimated asymptotic solution for a lossless wave model [Eq.
(2) with rigid boundaries] and the LIB HRTF measurements.Contralateral ear results for ð125; 277:2; 0Þ are not shown.up to 17 kHz with the LIB measurements-a match which can be considered acceptable for many practical purposes for HRTFs.This result holds even though the LIB-generated pressure pulses are not completely in the linear regime.From a validation perspective, this could be problematic; (linear) modeling errors could have canceled out with other unacknowledged errors in the validation process.However, since Figs.7 and 8 do not show any clear correlation between the validation error and the chosen radii r in Table I, and since the number of nonlinearities in the pressure pulses decreases with such radii (see supplemental material 52 ), the chances of cancellation of modeling-errors are considered small.Additional analyses are needed to further increase the confidence in this conclusion-e.g., employing a more conventional acoustical source at r ¼ 1 m.The most surprising result in the present study is the relatively good prediction for r LIB ¼ 5 cm for the ipsilateral ear: the acoustic pulse generated by the LIB is expected to behave non-linearly when it hits the surfaces of the head replica.One explanation could be that the pulse is weakly non-linear in the audible range even at 5 cm from the LIB.This would imply, provided the measurement errors are small, that the used non-linearity markers are poor indicators of a non-linear regime.Alternatively, the regime of the acoustic pulse might switch closer to linear after a boundary is encountered (at least for non-grazing incidence 69 ).For instance, we expect the incoming and reflected waves at a boundary not to be fully compatible with the principle of superposition.Otherwise, the authors are unaware of the frequency behavior of the pressure shock wave as it hits a boundary-most studies only address the free-field propagation, while the few studies involving boundaries involve broadband analyses. 37,69,70Such broadband results generally show increased high-frequency losses in the pressure waveform with surface roughness. 37,70Nevertheless, it is unclear then why the results are so poor for the contralateral ear at the same location; although very low frequencies are well predicted by the linear model [see Fig. 7(b)], the expected strong head-shadowing could push higher frequencies in the noise floor for both measurements and simulations [see, e.g., the large simulation CIs in Fig. 7(b)].Further investigations are required.
The present study also showed the difficulty in employing measured input-impedance values (e.g., b values).Results in Figs. 4 (top) and 7 show that: (i) the present methods to assess b values are too inaccurate for thorough validation studies, and (ii) the tube-measurements results should not be taken at their face value due to various instrumental biases.This is despite the good repeatability and reproducibility in the present tube-measurements compared to other studies. 64,71Note the presented nf-HRTFs results for b abs serve mostly as a guide for absorbing BC behavior due to various limitations in choosing b abs value in Eq. ( 7): first, the value is likely to be biased due to potential leaks in the impedance tube (see Sec. III C); second, using error propagation laws yields slightly less conservative intervals since, e.g., Dð1 À jRjÞ ( 1 À jRj is not necessarily true; finally, the asymptotic magnitude is likely slightly biased due to convergence issues in the discrete surface.
Provided the local-reactance assumption is correct, the b abs plots in Fig. 7 could also hint at the mechanism of notch generation: if the amplitude at the frequency of the notch increases while its Q-factor decreases with higher absorption at the boundary, this could suggest a destructive interference pattern with the direct sound; by contrast, if the notch deepens with increased absorption, it could indicate, e.g., a modal antiresonance or some more complicated interference.][74] This could also explain the slight frequency misprediction for some notches: modes require more time to form, allowing the speeds of the packets to converge to the expected values.
The HRTFs simulated with b abs generally show a better match in magnitude with the measurements above 17 kHz (see Fig. 7).This only indicates increased absorption at higher frequencies; the cause could be an impedance mismatch or other loss mechanisms (e.g., boundary viscous losses).Since any impedance measurement is currently highly inaccurate at such frequency range, calibration of b values might be a solution.However, other sources of losses need to be excluded first; otherwise, the calibration is done on the wrong target, and the predictive power is lost.
Based on the shape of the signed errors in Fig. 8, the validation results do not pinpoint any systematic deficiency of the employed PDE model relative to the spark nf-HRTFs up to about 17 kHz.Nevertheless, some of the notches are slightly mispredicted in both frequency and magnitude.Since most of such notches seem to be caused by destructive interference at higher frequencies, the main causes could be slight location/orientation mismatches, 52 geometrical errors, or phase delays introduced by the boundary impedance (note the impedance measurements in Fig. 4 indicate potential deviations from a rigid boundary together with illbehavior of the phase delay with increased frequency).Other boundary effects could be responsible such as nonlocal reactance.Based on the present study, it seems that the uncertainty in both simulation and measurements is larger for such notches, showing the further ill-conditioning of the notch estimation problem.
Finally, extrapolation of the current results to similar problems needs to be done with care; more validation studies are needed to confirm the present results for other pinna shapes.Extrapolating the results to in vivo HRTFs is also problematic.

VII. CONCLUSION
Employing the SRL FDTD scheme, the predicted formal nf-HRTF solutions agree in magnitude with LIB nf-HRTF measurements within about 1-2 dB up to around 17 kHz for a geometrically-validated stiff head replica named MoRa.Except for a few notches, no systematic deficiencies of the used lossless model (wave equation with rigid boundaries) could be observed relative to the measurements in such bandwidth.Above 17 kHz, increased uncertainty and real-world losses affected the results.
The present study introduced the usage of air LIB for nf-HRTF measurements.Such an acoustical source generates a highly omnidirectional (for our set-up, results showed a 4r deviation in directivity of about 1.6 dB), time-coherent, and phase-stable (our measurements show a 2r variability in group delay below about 1 ms for each direction) pressure pulse in the audible range.In the absence of unacknowledged error cancellation, the present findings showed that such pulse is weakly non-linear for the HRTF problem, such that the measured nf-HRTFs could be validated against the employed PDE model.Surprisingly, good validation results were found even for locations where the spark-generated acoustical pulse would be considered fairly non-linear.
Finally, the present study showed the difficulties in employing the surface impedance model to accurate validation studies: the errors and uncertainties in existing impedance measurements techniques are too large for reliable/ credible validation conclusions to be drawn for the sensitive HRTF problem.

ACKNOWLEDGMENTS
This research has received funding from Facebook Reality Labs.The computational resources provided by the CSC-IT Center for Science, Espoo, Finland, are also acknowledged.Pasi Karppinen from ProtoRhino Oy (Espoo, Finland) is thanked for the design and implementation of the Labview control software and help in designing the optical set-up.Ilkka Huhtakallio is thanked for help in designing and implementing the measurement set-up.Tyson Lindsay is thanked for the work and help in designing and 3D-printing the MoRa head replica and impedance tube samples.Li-Chung Chih is thanked for scanning work and help in the MoRa replica validation.Terry Cho is thanked for help and discussions.Julie Meyer is thanked for commenting on the manuscript.Jukka P€ atynen and Tapio Lokki are thanked for the tube-measurement MATLAB scripts.

FIG. 1 .
FIG. 1. (Color online) Laser setup and measurement points for the acoustical characterization.The expander separation from the laser is for visualization purposes.The scale of the measurement points is inconsistent with the other objects.Double 3D arrows represent the movement paths of the microphone relative to the LIB.

FIG. 4 .
FIG. 4. (Color online) Pooled reflection coefficient R magnitude (top), unwrapped R angle (middle), and unwrapped angle of the boundary surface impedance Z b (bottom).E is the expectation.Top plot shows 95% CI of the measured jRj values U jRj;95% (standard errors of the two means are <10 À3 dB).The two lines in the middle/bottom plots show results for each tube size.

FIG. 6 .
FIG. 6. (Color online) Differences in HRTF magnitude measured with the microphones surrounded by viscoelastic material of different hardness.Transparent fills represent 95% CIs (see text for details).Dotted lines represent the HRTFs for the right ear.

FIG. 7 .
FIG. 7. (Color online) Measured nf-HRTF magnitude and predicted asymptotic magnitude g FDTD for b ¼ 0 and b abs [see Eq. (7)].The CIs are for the mean values: based on standard error for the measurement and the BC a pair-bootstrapped 95% for the computed solution (Ref.41).The same color coding for the direction text as in Fig. 2 and Table I.CIs not propagated through the log function.(a) Validation direction 1: (45 , 0 ).(b) Validation direction 2: (270 , 45 ).

TABLE I .
Validation directions.r LIB represents the distance to the surface of the MoRa head replica.Voxel sizes DX values are given only for the P ear simulation and are rounded to 2 decimal places.The same sampling frequencies f s were used for both P ear and P ref simulations.