Finite-difference time-domain simulations: Verification on head-related transfer functions of a rigid sphere model

: This paper presents a veriﬁcation procedure for ﬁnite-difference time-domain-simulated head-related transfer functions (HRTFs) from a simpliﬁed model of a human head, a sphere. The analytic solution required by the code veriﬁcation is computed with the multipole reexpansion technique and used to estimate convergence rates. A solution veriﬁcation process follows in which asymptotic predictions are computed. For the HRTFs considered and employed grids, results show that the convergence rates attain the expected ﬁrst-order accuracy at lower frequencies, after which scattered estimates are observed. Results also reveal that the asymptotic predictions are accurate up to 10kHz, after which bias is observed. V C 2022 Author(s). All article content, except where otherwise noted, is licensed under a Creative Commons Attribution (CC BY) license (http://creativecommons.org/licenses/by/4.0/) .


Introduction
The finite-difference time-domain (FDTD) method has been and continues to be widely used to solve the acoustic wave equation in various application areas (e.g., underwater acoustics, 1 room acoustics 2 ).To ensure that the simulation results can be trusted, the reliability and the accuracy of the implemented method and computed results should be assessed.Such an assessment can be done via a verification procedure, which usually comprises a code verification followed by a solution verification. 3The purpose of the code verification is to assess the correctness of the discrete mathematical model implementation, while the solution verification aims at quantifying the numerical errors in the computed solution for a given problem.In practice, the code verification procedure typically consists in estimating convergence rates and requires a known analytic solution.The solution verification, which is done in the absence of a known analytic solution, usually consists in either estimating the discretization error or predicting asymptotic solutions (i.e., computing asymptotic predictions) from a series of simulations.
One application for which the use of FDTD simulations is particularly attractive, as an alternative to measurements, is the prediction of head-related transfer functions (HRTFs).Verification in that context is therefore relevant, especially considering that the "operational validity" 4 of a simulation method is context-dependent.Moreover, studies focused on verification tasks in the context of HRTF predictions are scarce. 5In Ref. 6, FDTD-simulated HRTFs of a single sphere were only visually compared to analytic solutions, and no convergence rate was estimated nor asymptotic solution predicted.Reference 7 included a convergence rate assessment for another sphere model but did not perform a solution verification.Reference 5 presented finitedifference predicted asymptotic solutions for pinna-related transfer functions.However, the code verification in Ref. 5 was done in the absence of the voxelization error.As such, the case of the single sphere model is relevant since code verification is possible as there exists an analytic solution for the problem of scattering from such a model.It also provides a case for verification where the voxelization error is present.Last, such a model, as it represents the simplest approximation of a human head, gives a foundation on what to expect for more complex geometries.

Rigid sphere model
This paper considers a sphere of radius 8.25 cm with rigid (specific admittance ¼ 0) boundaries.The spherical coordinates (r, h, /) using the convention from Ref. 8 were adopted, where r is the radial distance (or radius), h is the colatitude, and / is the azimuthal angle defined counterclockwise (as seen from the top of the sphere) and whose range is 0 / < 360 .Seven positions on the sphere surface were selected for investigation.Relative to the center of the sphere, taken as the center of the coordinate system, these positions were at a colatitude h ¼ 90 and azimuthal angles / ranging from 0 to 180 with 30 increments.The frequency of analysis varied from 250 Hz to 20 kHz with 125 Hz increments.A monopole source was positioned at two radial distances of 82.5 cm and 1.65 m from the sphere center at (h, /) ¼ (90 , 0 ).The rationale behind the choice of the source distances was that near-(within 1 m distance from the sphere, as in Ref. 9) and far-field differences could be investigated.The positions of investigation on the sphere and source positions are shown in Fig. 1(a).

Analytic solution
The problem of acoustic wave scattering from a sphere can be solved using numerical computational methods such as FDTD.Another set of methods uses analytic expressions whose accuracy is typically controlled by the truncation order of the (infinite) summations involved.The computational technique identified as the multipole reexpansion technique, 10 which is based on the T-matrix method, 11 belongs to the second set of methods and is chosen here to provide the analytic solution of the problem of scattering from the rigid sphere model described in Sec.2.1.The interested reader is referred to Refs. 10 and 12 for a thorough description of the method and definitions of the underlying mathematical functions used in the computations.The multipole reexpansion method was preferred over the method from Ref. 13 as rigorous error bounds were found in Ref. 12, whereas Ref. 13 did not provide such error bounds.
In the present context of HRTF predictions, the multipole reexpansion technique is used to compute the function H given in Eq. ( 1), which defines the HRTFs of the single sphere, where f denotes the frequency, and P S is the Fourier-transformed pressure on the sphere at the position indicated by (r, h, /).P inc denotes the Fourier-transformed pressure of the incident wave measured in the free field at the sphere center.
To ensure that the results from the implementation of the multipole reexpansion method by the present paper's authors were in line with previously published results from the authors who developed the method, comparisons with Fig. 3 from Ref. 10 were carried out.These comparisons demonstrated good agreement between the two implementation results (results not shown).Furthermore, the case of the scattering from the rigid sphere described in Sec.2.1 was compared with the results from Ref. 13.For this latter comparison, a truncation order of 64 was employed (the rationale behind this choice is explained in the following paragraph) for the multipole reexpansion technique, while the series from Ref. 13 was truncated when the ratio of the subsequent terms in the series was smaller than MATLAB's machine epsilon in double precision e % 2.2 Â 10 À16 .The maximum absolute difference between the two implementation results was around 10 À12 .
As previously mentioned, the computation of the analytic solution requires the truncation of the summations involved to a finite order.For the problem of scattering from a sphere, the truncation number required to provide computations with a prescribed accuracy can be evaluated rigorously. 12Using the error bound Eq. (9.1.42)from Ref. 12, the truncation number was found to be 64 to attain a prescribed error of (k max =4pÞ Â e in the incident and the scattered fields, separately, across all frequencies and the two source positions, where k max represents the maximum wavenumber considered.The final error bound on 20 log 10 ðH analytic Þ was estimated using interval arithmetic, assuming 2(k max =4pÞ Â e absolute error in jP S j and ðk max =4pÞ Â e in jP inc j: the error was evaluated to be below 1.3 Â 10 À11 dB across all conditions.The computed analytic solution with such prescribed accuracy is shown in Fig. 1(b) for the seven HRTF directions and for the near-field source.

Convergence rate estimation
As mentioned in Sec. 1, code verification is concerned with the assessment of convergence rates.Essentially, this implies analyzing the discretization error behavior by comparing the observed (from FDTD simulations) with the expected (i.e., the formal order of the employed scheme) rate at which the discretization error decreases as a function of spatial grid spacing. 5Note that a first-order convergence rate is expected due to the first-order voxelization error 14 at the boundaries as well as the nearest neighbor error (i.e., the computations are done at voxel centers, and no spatial interpolation is done; see, e.g., Table 1 from Ref. 5).
The basic equation to evaluate the discretization error 15,16 is the asymptotic expansion given in Eq. ( 2), where cð f Þ represents the coefficient of the principal error term, which is independent of the spatial grid spacing X, HOT denotes the higher-order terms, and q is the convergence rate.H analytic and H i are the functions H computed using the multipole reexpansion technique and simulated on the ith FDTD spatial grid spacing X i , respectively.Considering that it is possible to estimate the discretization error in a least squares sense if the number of FDTD grids is greater than or equal to four, 15 a regression model based on Eq. ( 2) (with the HOT neglected) linearized using a logarithm function was used for estimating q.More specifically, the slope resulting from the least squares fit of the regression model gave the q estimates.Additionally, the weights introduced in Appendix B.1 from Ref. 15 were also applied to the linear regression model.
Alternatively, by neglecting the HOT of Eq. ( 2), the convergence rate q can be estimated as follows: where s ¼ X iþ1 =X i is the FDTD spatial grid refinement ratio relating the spacings from two consecutive FDTD spatial grids.H i and H iþ1 are the functions H simulated on the FDTD spatial grid spacings X i and X iþ1 , respectively, with The results from these two convergence rate estimation formulas are presented in Sec.2.5.

Simulations
An FDTD solver using graphics processing units (GPUs) and a message-passing interface (i.e., using parallel processing) was utilized to run the simulations.The FDTD solver was previously verified on a free-field computational domain (i.e., no scattering object, thus without the voxelization error), using the methods described in Ref. 5, as well as for free-field propagating spherical wavefronts emanating from a point source (results not shown).Here, the code verification analysis from Ref. 5 is extended to when the voxelization error is present.The simulations were run on 10 or 12 nodes containing 8 or 16 GPUs (Tesla V100, Nvidia, Santa Clara, CA) each.The standard rectilinear scheme of the FDTD method was used at its stability limit (i.e., with a Courant number set to k ¼ 1= ffiffi ffi 3 p ).The simulations were run in single and double precision (see the results in Secs.2.5.1 and 2.5.2, respectively).The FDTD code computes the update equations at each air voxel within a finite-volume interpretation 17 on a uniform and structured grid (see, e.g., Ref. 5 for an explicit update).
The voxelization was made using a surface-conservative voxelization algorithm. 18To ensure that no receiver (i.e., a grid node at which the pressure is captured) ended up on a solid voxel of the voxelized sphere mesh, a radius larger than that of the sphere was considered to position the receivers at the center of air voxels of the FDTD spatial grid.The radius r i for the receiver positions was determined using Eq. ( 4), where r sphere ¼ 8.25 cm is the sphere radius, X i is the FDTD spatial grid spacing, and e is machine epsilon in double precision.Note that r i in Eq. ( 4) converges to r sphere (within e) as the spatial grid spacing X i gets smaller.
A soft source with a Gaussian pulse [GðtÞ ¼ e ðÀðtÀlÞ 2 =2r 2 Þ , where r 2 ¼ 0.034 Â 10 À8 and l ¼ 1.8 Â 10 À4 ] as the driving function was employed to excite the FDTD grids.A 4 m Â 4 m Â 4 m rigid box was designed around the sphere model.The simulation time was set to 8 ms such that the scattered sound field from the sphere was fully captured by the receivers without including any reflection from the box.The FDTD spatial grid refinement ratio s ¼ X iþ1 =X i , where X i < X iþ1 , was chosen to be 1.1 (minimum recommended in practice) since including a large number of grids is beneficial in the context of code verification tasks 15 (see also Sec. 2.3).
The speed of sound c for the air voxels was chosen such that the temperature was 20 C and the atmospheric pressure was 101.325 kPa in dry air (i.e., with density 1.2041 kg/m) for all simulations and thus approximately equaled 343.4 m/s.This allowed the study of 19 spatial grid spacings X ranging from approximately 0.75 to 4.20 mm (corresponding to temporal sampling frequencies f s ranging from 793 125 to 142 750 Hz) while keeping the simulation time constant throughout the series of simulations.Also, by keeping the simulation time constant, the frequency resolution defined as the ratio of the temporal sampling frequency of the simulations f s over the number of time steps N was kept constant between grids and equaled 125 Hz.Because the present FDTD solver requires an integer sampling frequency f s , keeping a constant frequency resolution while using a refinement ratio of s ¼ 1.1 was not always possible.Consequently, f s was readjusted to keep the frequency resolution constant, which led to X values rounded to two decimals.For the same reason, the refinement ratio between two consecutive grids was not exactly 1.1 but varied between 1.1 6 4.8 Â 10 À4 after re-adjusting f s .
The discrete Fourier transform of each of the 19 FDTD-simulated Gaussian pulse responses was computed using a fast Fourier transform (FFT) algorithm with an FFT size equal to N. The FFT was computed separately for the responses generated in the presence of the sphere model to obtain P S and for those generated in the free field to obtain P inc .Finally, P S was divided by P inc as in Eq. ( 1) to obtain H.

Single precision
The two convergence rate estimation formulas presented in Sec.2.3 were utilized with the 19 FDTD simulations run in single precision described in Sec.2.4.The 95% confidence intervals (CIs) on the q estimates from the linear regression model based on Eq. ( 2) with the HOT neglected were also computed using the bias-corrected and accelerated (BC a ) bootstrap method (with 5000 replicates) introduced in Ref. 19.The results for two positions on the sphere and the near-field source position are shown in Figs.2(a)-2(d).The results for the second source position were similar (not shown).
As can be seen in Figs.2(a) and 2(b), the q estimates from Eq. (3) are highly scattered, which suggests either that the asymptotic range was not attained or that Eq. ( 3) is not a reliable indicator of asymptotic range for the used models (i.e., the used FDTD scheme including the voxelization error for solving the wave equation).The extent of the scatter seen in Figs.2(a) and 2(b) was similar for the rest of the HRTF directions.The scatter in the q estimates was, however, considerably reduced at lower frequencies when the linear regression model was applied.As can also be seen in Figs.2(c) and 2(d), applying the weights from Ref. 15 to the linear regression model did not further improve the q estimates.When also comparing Fig. 2(c) with Fig. 2(d), the extent of the scatter seems to vary depending on the HRTF direction considered when the linear regression model is applied.Another observation concerns the 95% CIs on the q estimates, which increase as a function of scatter in the data.Despite the scattered estimates seen at higher frequencies, the expected firstorder accuracy is nonetheless attained at least up to 1125 Hz for most of the HRTF directions (slightly worse q estimates were found in that lower frequency range for the far-field source for / ¼ 60 ), which implies that the used FDTD solver behaves as expected in that lower frequency range.Also, note that the extent of this lower frequency range slightly varied across HRTF directions (e.g., it was up to 2 kHz for the near-field source for / ¼ 30 ).

Double precision
To verify if using double rather than single precision improved the q estimates, the 19 FDTD simulations run in double precision described in Sec.2.4 were first compared with the 19 FDTD simulations run in single precision.The maximum of the absolute error between the 19 H functions (converted to the dB scale) run in single and double precision across all conditions (the two source positions and the seven positions of investigation) was less than 6.9 Â 10 À5 dB up to 10 kHz and less than 1.9 Â 10 À3 dB for higher frequencies.
Similarly to Sec. 2.5.1, the q estimates from the linear regression model based on Eq. ( 2) (with the HOT neglected) were recomputed using the simulations run using double precision.Visual inspection of the results (not shown) revealed the q estimated from the simulations run in single and double precision were almost identical across the entire 250 Hz-20 kHz frequency range for the HRTF directions and two source positions considered.

Spatial interpolation
To verify if avoiding location-related inconsistencies improved the convergence rate estimates, second-order accurate trilinear spatial interpolation was applied to both the source and receiver positions for each of the 19 FDTD grids run in single precision.The interpolation was done for the near-field source position only.Unlike the previous simulation series where the receivers were positioned using Eq. ( 4) (see Sec. 2.4), here the interpolated receiver positions followed Eq. ( 4) solely calculated with the coarsest spatial grid spacing of the simulation series, that is, X i % 4.20 mm (corresponding to a sampling frequency of f s ¼ 142 750 Hz).Similarly, the analytic solution was recomputed for the same receiver positions located outside of the sphere surface.
The q estimates from Eq. ( 3) and from the linear regression model based on Eq. ( 2) (with the HOT neglected) were recomputed.The q estimated using Eq. ( 3) for the non-interpolated and the interpolated data were similar (results not shown).However, the q estimated using the linear regression model based on Eq. ( 2) (with the HOT neglected) were generally improved for the interpolated data, in comparison to the non-interpolated data.This can be seen for two HRTF directions by comparing Figs.2(c) and 2(d) with Figs.2(e) and 2(f).From Figs. 2(e) and 2(f), it can be seen that not only are the q estimates brought closer to the first-order accuracy, but the 95% CIs are also smaller across a wider frequency range in comparison to Figs. 2(c) and 2(d).The frequency range in which the q estimates attained the first-order accuracy was widened due to interpolation for all the HRTF directions except for / ¼ 150 and 180 , by 6375 Hz at most (for / ¼ 60 ).For all the HRTF directions except for / ¼ 120 , the 95% CIs were smaller due to interpolation for at most and at least 98% and 58% of the entire 250 Hz-20 kHz frequency range, respectively.It is recalled that the convergence of the voxelized geometry is of first order. 5As such, the first-order accuracy observed for the interpolated data implies that the voxelization error is dominant.

Solution verification
As stated in Sec. 1, solution verification, which is preceded by code verification, is concerned with the quantification of the numerical errors in the computed solution.Common techniques focus on the quantification of the discretization error, which is usually assumed dominant (Ref.3, p. 286), by employing a series of simulations run with different spatial grid spacings related to each other by a grid refinement ratio.Alternatively and as in the present work, asymptotic predictions can be computed using the same kind of simulation series.Here, although an accurate analytic solution is available, the solution verification procedure is only applied to the FDTD-simulated results, which implies assuming that the analytic solution is unknown.The outputs of the solution verification procedure can then be compared with the analytic solution, effectively assessing the accuracy of the employed solution verification procedure.Fig. 2. Convergence rates q estimated using Eq. for positions (a) / ¼ 0 and (b) / ¼ 180 for the non-interpolated data.Convergence rates q were estimated using the non-weighted (N-W) and weighted (W) linear regression model based on Eq. ( 2) with the HOT neglected and the 95% CIs from the BC a bootstrap method from Ref. 19 for the non-interpolated data and for positions (c) / ¼ 0 and (d) / ¼ 180 .(e) and (f) exhibit the same q estimates as (c) and (d), respectively, but for the interpolated data (see Sec. 2.5.3).In (a)-(f), the expected first-order rate together with the typical tolerance of 610% for convergence rates (Ref.3, p. 326) are also indicated.The source was located at a radial distance of 82.5 cm from the sphere center.
The solution verification was performed using the 19 FDTD simulations run in single precision described in Sec.2.4.The asymptotic prediction was computed from the intercept resulting from the least squares fit of the first-order asymptotic model expressed in Eq. ( 5), 5

H i ðr
where H asymptotic ðr; h; /; f Þ denotes the sought asymptotic prediction and bð f Þ represents the coefficient of the principal error term.H i denotes the function H simulated on the FDTD spatial grid spacing X i .
It is worth mentioning that the first-order asymptotic model was assumed to fail when the magnitude of the asymptotic prediction was below machine epsilon in single precision (%1.2 Â 10 À7 ).The 95% CIs on the asymptotic predictions were estimated using the BC a bootstrap method 19 with 5000 replicates.Both non-weighted and weighted asymptotic predictions from the FDTD-simulated grids as well as the 95% CIs are shown (in dB scale) in Figs.3(a)-3(c) for the near-field source position and three HRTF directions.Results (not shown) were similar for the far-field source position and the other HRTF directions.Figure 3 shows that the 95% CIs increase as a function of frequency independently of the HRTF direction, which is in line with the results reported in Ref. 5. For all the conditions considered (source distance, HRTF direction), the weighted asymptotic prediction had slightly smaller 95% CIs (averaged across the entire 250 Hz-20 kHz frequency range) than the non-weighted asymptotic prediction.From Fig. 3, it can also be seen that the asymptotic predictions successfully recover the analytic solution for the employed grids up to about 10 kHz, for which the CIs are within 3.0 dB, after which prediction bias is observed together with increased CIs.As such, larger CIs suggest a mediocre least squares fit of the first-order asymptotic model and can be useful to indicate low accuracy of the asymptotic prediction and large deviations from the asymptotic range.Although it is difficult to extrapolate to more complex acoustical scenarios since the results depend on the formal solution and the set of employed grids fX i g, a >3.0 dB threshold for the CIs could provide a conservative baseline indicator for asymptotic prediction bias.
The solution verification procedure was re-applied to the 19 FDTD simulations run in double precision.Both non-weighted and weighted asymptotic predictions as well as their 95% CIs (not shown) were almost identical to these computed using the FDTD simulations run in single precision.This implies that in the present modeling and simulation context, the solution verification procedure is largely unaffected by round-off errors.
The asymptotic predictions were also computed for the interpolated data.While failure of the first-order asymptotic model was observed for the positions / ¼ 0 around the notch at 13 kHz [shown in Fig. 3(d)] and / ¼ 30 , the Fig. 3. Individual FDTD responses with respective spatial grid spacing X, non-weighted (N-W) and weighted (W) asymptotic predictions computed from the first-order asymptotic model from Eq. ( 5), and the 95% CIs for the non-interpolated data and for positions (a) / ¼ 0 , (b) / ¼ 120 , (c) / ¼ 180 .Figure 3(d) exhibits the same asymptotic predictions as in Fig. 3(a) but for the interpolated data (notice the failure of the first-order asymptotic model between 11 250 and 17 250 Hz).In (a)-(d), the source was located at a radial distance of 82.5 cm from the sphere center.The analytic solution is solely used here for visual reference, and the uncertainty through the logarithm function was not propagated.
results for the other HRTF directions (not shown) were similar to those of the non-interpolated data.Figure 3(d) suggests that such a failure is due to the inclusion of computations far from the asymptotic range in the asymptotic model: notice the computations in light gray that cause prediction bias.This indicates that the solution verification procedure is mostly insensitive to interpolation errors in the present context.
Based on the results from the solution verification, it is recommended that using a series of grids, providing that the grids X i are sufficiently small relative to the curvature and size of the input geometry (i.e., grids are in the asymptotic range), to compute asymptotic predictions be preferred over the simulation of individual FDTD responses.

Conclusion
This paper presented a verification procedure for FDTD-simulated HRTFs from a rigid sphere model.The results from the code verification showed scattered convergence rate estimates that were improved when a linear regression model was used as the estimation procedure.Spatial interpolation further improved the convergence rate estimates by bringing them closer to the expected first-order accuracy.Using double instead of single precision to run the FDTD simulations did not seem to improve the convergence rate estimates.The extent of the scatter in the convergence rates estimated using the linear regression model was found to vary with the HRTF direction considered.The results from the solution verification showed that the 95% CIs on the asymptotic predictions increased as a function of frequency regardless of the condition considered (source distance, HRTF direction).This also indicates that the bias in the asymptotic predictions increased as the confidence in the predictions decreased with frequency.For all the conditions, the weighted asymptotic prediction had slightly smaller 95% CIs (averaged across the entire frequency range considered) than the non-weighted asymptotic prediction.For the numerical modeling method employed and the present problem at hand, interpolation and double precision did not seem to improve the accuracy of the asymptotic predictions and, as such, did not seem to be required for the assessed solution verification procedure.It was also shown that the solution verification procedure is able to recover the analytic solution when the CIs are reasonably small and the computed solutions are not largely far from the asymptotic range.In the present work, this was observed for frequencies up to 10 kHz for all the HRTF directions considered.The results of the verification procedure were similar for the near-and far-field source distances.

Fig. 1 .
Fig. 1.(a) Top view of the rigid sphere model with the positions of investigation (orange dots) located at a colatitude h ¼ 90 and azimuthal angles / ranging from 0 to 180 with 30 increments and the source positions (black dots).The indicated dimensions are not to scale.(b)The log-magnitude of H analytic as a function of frequency for the seven positions of investigation on the sphere surface and for the near-field source position.H analytic denotes the function H from Eq. (1) computed using the multipole reexpansion technique (as implemented by the present paper's authors).The computations were done in double precision using a truncation order of 64.