Pinna-related transfer functions and lossless wave equation using finite-difference methods

A common approach when employing discrete mathematical models is to assess the reliability and credibility of the computation of interest through a process known as solution veriﬁcation. Present-day computed head-related transfer functions (HRTFs) seem to lack robust and reliable assessments of the numerical errors embedded in the results which makes validation of wave-based models difﬁcult. This process requires a good understanding of the involved sources of error which are systematically reviewed here. The current work aims to quantify the pinna-related high-frequency computational errors in the context of HRTFs and wave-based simulations with ﬁnite-difference models. As a prerequisite for solution veriﬁcation, code veriﬁcation assesses the reliability of the proposed implementation. In this paper, known and manufactured formal solutions are used and tailored for the wave equation and frequency-independent boundary conditions inside a rectangular room of uniform acoustic wall-impedance. Asymptotic estimates for pinna acoustics are predicted in the frequency domain based on regression models and a convergence study on sub-millimeter grids. Results show an increasing uncertainty with frequency and a signiﬁcant


I. INTRODUCTION
The head-related transfer functions (HRTFs) represent the acoustical paths from a stationary sound source to the ears of the listener under free field conditions, encoding the anatomical filtering effects which embed the auditory cues used for localizing sound sources.
Unless cue remapping occurs, 1 HRTFs are usually not perceptually transferable. 2This is mainly due to the high degree of personalization of the human pinna. 3Currently, the main two methods for obtaining an estimate of individualized HRTFs are measurements and simulations.5][6][7] Although HRTF measurements can nowadays be done within a relatively short time, [8][9][10] they are still, in many ways, impractical: high accuracy could be challenging, specialized equipment prevents scalability, while subjects undergo a generally tedious measurement session.
Wave-based numerical simulations are an alternative which could eventually offer greater flexibility when compared to measurements.Presently, the boundary element method (BEM) [11][12][13][14][15] and the finite difference time domain (FDTD) [16][17][18] methods are the most common HRTF simulation methods.Despite the many attractive properties of the HRTF numerical simulations, their validity and reliability in the full audible frequency bandwidth has not yet been established using strong validation studies as those defined in the verification and validation (V&V) literature. 19,20 Electronic mail: sebastian.prepelita@aalto.fi b Current address: Hefio Ltd., Otakaari 5 A, 02150 Espoo, Finland.
The most challenging HRTF phenomena to be simulated are the pinna effects which occur above 3 kHz. 22The pinnarelated transfer function (PRTF) provides the isolated information of the pinna acoustics and can be used to focus the analysis to such a challenging problem.
Any computation will contain approximations and specific errors which could render the result unreliable.Present HRTF computations seem to lack a strong formal assessment of reliability and credibility of the computed solutions.Consequently, it not only becomes difficult to assess the predictive quality of the wave-based model employed, but also wrong conclusions could be reached.The present study aims to quantify the quality of computed wave-based PRTF/ HRTF solutions with focus on FDTD methods.Formal solutions and their precision are estimated by asymptotically extrapolating the results from computations on multiple submillimeter grids.In turn, the uncertainty between HRTF measurements and simulations can be quantified to the complementary benefit of both approaches. 23his work is the first part of a comprehensive V&V study.Validation with measurements of the present PRTF results will be part of a separate validation study.

II. BACKGROUND
Linear acoustics is described by a linear and wellposed continuous problem (in the Hadamard sense).Thus, any consistent and stable approximation will theoretically converge to the same solution (invoking the Lax-Richtmyer equivalence theorem), independent of the mathematical formulation of the continuous model and of the discrete model employed (e.g., BEM, FDTD).Depending on the formulation and properties of the continuous solution, any computational implementation and execution will involve model-specific deviations from the formal behavior.In practice, such computational errors could interact in a complex manner thus requiring empirical reliability assessment of the final computations.
Here, a simple explicit and stable 24 finite-difference method is chosen which is amenable to parallelization.Its simplicity and extensive study facilitates the understanding of various error sources present in a final computation.Moreover, the finite-difference literature is extensive and offers formal methods/models to study the usually-dominant (Ref.20, p. 286) discretization error.Such methods can be transferred from other fields and advantageously employed to empirically estimate the quality of the computations, given the complex PRTF problem at hand.

A. PRTFs
The PRTFs are formally defined in the frequency domain as the following transfer function: PRTFðr 1 ; XÞ ¼ P ear ðr 1 ; XÞ P ref ðr 1 ; XÞ ; where r 1 2 R 3 denotes the continuous location of the source, P ear ðXÞ represents the Fourier transform of the captured pressure at a fixed location of interest inside a pinna cavity, and P ref ðXÞ is the Fourier transform of a recorded pressure with the pinna absent.
At the boundary surface(s) @D, the following frequencyindependent resistive boundary condition (BC) of local reaction is employed for the pressure field: where n represents the outward (pointing, locally, from D into any solid in contact with D) normal to the surface boundary at point r b 2 @D, and b is the specific acoustic admittance at the boundary (Ref.

C. Discrete models
In this study, an FDTD mathematical model is used to discretize Eq. (2).Rigorously, the update discretizes the integral form of Eq. ( 2), however for the uniform Cartesian grid used, the discretization of the strong form in Eq. ( 2) is found 24 for the interior update.
The second-order accurate standard rectilinear (SRL) explicit update is obtained on the interior and employed on a cubic Cartesian grid with voxel size of DX.A frequency independent impedance can be set for the locally reacting staircased boundaries using an "upwind" approximation for the pressure field at the boundary. 28The code supports three types of sources: soft, hard, and transparent. 29The explicit FDTD update can be derived using finite volume (FV) methods 24 and can compactly be written for pressure and a soft source as where P i represents the discretized pressure field at the 3D Cartesian index i, j indexes the neighboring air voxels for P i ; k C ¼ cDT=DX is the Courant number (c is the sound speed), N i 2 f0; 1; 2; …; 6g represents the total number of solid (i.e., boundary) voxels adjacent to P i [for N i ¼ 6 the sum in Eq. ( 4) is assumed to be 0], and F i represents the local forcing applied at the center of voxel i.In Eq. ( 4), for computer memory considerations, a unique specific acoustic admittance b i is assumed for all boundary surfaces of each voxel i.
For simplicity and unambiguous definition of the grid spacing DX, the boundaries are discretized on the same uniform 3D Cartesian grid resulting in stair-stepped boundaries which are obtained from the mesh by using a voxelization algorithm.The adopted voxelization algorithm is the conservative 30 one since, although not ideal, it is well-defined and convenient to use.

Distributed computing model
The simulation domain was split against the Cartesian z axis: D is partitioned in domains distributed across computing nodes, while each domain is partitioned in subdomains across Graphics Processing Unit (GPU) devices available at each compute node.See more details about the implementation in Appendix B.
To avoid extra difficulties for a potential validation study such as convergence rates and correspondence with reality, the PRTF domain will be chosen without a commonly-employed perfectly matched layer absorbing bounding boundary.Consequently, the PRTF domain size is scaled such that no reflection from the room walls arrives at the receiver locations.

D. Simulation errors
Depending on modeling and computational choices, there are multiple sources of error, 31 which could coexist in a simulated result and affect the computations.Such errors make it highly unlikely that the formal (i.e., asymptotic, as DX !0) solution is computed in practice: given a discrete model and its computational implementation, the discrete solution will at best be within a ball Bðpðr; tÞ; eÞ of radius e proportional to the machine epsilon away from the unique formal solution p(r,t).As such, it is useful to identify and understand the behavior of the main sources of error for the model and problem at hand, also considering that error cancellation could occur (Ref.20, p. 382) in a computation when multiple sources of error are present, together with other unacknowledged 23 errors.
Table I summarizes the following sources of error: Discretization error is caused by the approximation of the continuous operators in the employed PDE in Eqs. ( 2) and ( 3) by discrete operators.
A first-order truncation error is found for the BC which propagates over interior solutions (Ref.33, p. 159).For example, the local truncation error for a velocity v b used in the flux at a boundary orthogonal to the x axis is where p b is the pressure at the boundary surface b, P j is the pressure at the center of a neighboring voxel j, b b is the specific acoustic admittance at b, and HOT stands for 19 "higher order terms" in DX.
Analysis of the global discretization error is generally making use of a series expansion of the truncation error.For the employed scheme and temporarily ignoring the initial conditions and the boundary modeling, as DX,DT !0, an FDTD-simulated pressure e p r 1 ðtÞ at r ¼ r 1 2 R 3 can be expressed as (Ref.33  33 ) while the temporal dependence of x;r 1 ðtÞ is given by p r 1 ðtÞ. 34Thus, even a discrete delta sequence can be used as the excitation signal without aliasing of the discretetime Fourier transform (DTFT) of the error terms E t=x ðxÞ ¼ DTFTf t=x ðtÞg and HOT(x) in Eq. ( 6).Asymptotic range is attained when the deviation of the discrete solution from the formal solution behaves according to the expected leading 19,35 space and/or time truncation error terms.In the following, it is assumed the discrete solutions are stepped at stable values of k C with negligible global interactions (e.g., cancellation) between the spatial and temporal truncation errors.
Round-off errors emerge due to the finite precision representation of numbers in the FDTD solver.A priori analysis shows divergence of results: assuming an additive round-off error, a theoretical O(1/DX) behavior is found for an ordinary differential equation (see, e.g., May and Noye 36 ).In practice, the behavior of the round-off error is less understood. 37nitial conditions could affect both the convergence rates and the magnitude of the error terms (see, e.g., Chap. 10 in the work by Strikwerda 38 ).Here, two main aspects are of interest: (i) spectrum inconsistencies [e.g., a discrete delta has a continuous temporal gain of O(DT) and spatial gain 39 of O(DX 3 )]; and (ii) the pollution error 40 [e.g., the local principal spatial error grows as (Ref.84, p. 5) Oðx 4 DX 2 Þ].The latter is avoided by using "smooth enough" (i.e., low-passed) sampled initial conditions. 41oundary modeling errors are related to the nontrivial grid-generation problem and emerge here due to the voxelized boundary surface f @DðDXÞ.The voxelization process will induce: (1) a volume convergence of OðDX 3 Þ; (2) a first-order 42,43 distance convergence non-parallel to the boundary @D behaving as OðxDXÞ; 44 (3) spurious local "excess scattering" 45 of high-frequency waves that basically pollute the solution; and (4) unlikely convergence of the voxelized surface area 46 f @DðDXÞ for a complex geometry such as a human pinna.
Such errors are scheme-dependent, 47 generally converge non-smoothly in DX, and could negatively impact the observed convergence rates even below first-order. 48,49ocation-related inconsistencies will cause a firstorder error due to source/receiver locations at voxel centers.For a point source/receiver (Ref.19, p. 173) in the FV interpretation of the scheme, 50 an unpractical odd grid refinement ratio r ¼ DX old =DX new of r ¼ 2k þ 1 with k 2 N Ã would be required.Alternatively, distributed sources/receivers could be employed, or interpolation used as post-processing.

III. METHODS
Verification procedures are broadly concerned with the precision and/or accuracy of actual computations. 23In the V&V literature, they involve two main procedures: • code verification: assesses whether the computational model is a faithful implementation of the discrete mathematical model; 23 • solution verification: assesses the accuracy, relative to the formal solution, of computed solutions by the/a previously verified computational model.
Both procedures imply empirically studying how refining the grid (i.e., DX; DT) affects the computed solutions.
The formal (time, space) order of accuracy of the scheme will initially be denoted by (g, q).More details on mathematical notation can be found in Appendix A.

A. Code verification
The order of accuracy test is usually considered the most difficult to satisfy code verification procedure for PDE simulations. 35It involves comparing the observed order of accuracy (i.e., the rate by which the discretization error empirically decreases with DX) with the formal order of the employed scheme (i.e., the expected rate of decrease in error with DX emerging from a truncation error analysis).This is usually 19,35,51 achieved by comparing computations against band-and energy-limited continuous solutions: (i) simple exact solutions (ESs) or (ii) manufactured solutions (MSs) (i.e., solutions which do not satisfy the given PDE).
To avoid voxelization errors (see Sec. II D), the utilized continuous solutions are chosen inside an axis-aligned rectangular room of dimensions ½L x ; L y ; L z T with uniform admittance on all the walls.
For the comparison with the analytic solutions, the discretized l 2 DX -norm is used to compute the global error where N ¼ N x N y N z ; DV ¼ ðDXÞ 3 , m indexes all the voxels residing in the interior air volume DX Á ½0; x=y=z DX; M ¼ DX½iþ0:5;jþ0:5;k þ0:5 T ; P t n m is the computed pressure field at time t n ¼ nDT ¼ t sim , and p is the manufactured or ES.
The initial conditions are used to drive the simulations and are applied at all points M of the interior air volume and are set based on the continuous solutions: pðM; 0Þ and pðM; DTÞ.For the MS, the forcing terms are implemented as soft sources and are evaluated in the center of the voxels from 2DT onwards.See more details in Appendix B.

ES
To avoid possible difficulties (Ref.19, p. 82) arising from using existing power series solutions in rectangular rooms, 52,53 a first simple approach to code verification is to assume the walls are perfectly reflecting for which an eigenvalue solution exists (Ref.26, p. 571).The ES was chosen as the following closed-form eigenfunction: where No analytic forcing term exists for the ES [i.e., f ES ðr; tÞ ¼ 0].

MS
Although rigid-wall PRTF computations are assessed in the present work, results of code verification for various boundary impedance values are presented to (i) have full code verification coverage for the solver, and (ii) examine a toy problem in which multiple sources of error having different asymptotic behaviors are present.
A solution is manufactured inside a rectangular room for b 6 ¼ 0. The method of MSs solves the PDE problem backwards; 51 the forcing term in Eq. ( 2) is determined by applying the differential operators of the PDE from Eq. ( 2) to the a priori chosen solution.Moreover, to test the boundary update, the MS should satisfy the BCs.The MS is chosen in the form where P 0 2 R represents the amplitude, k d 2 R can be viewed as the "decay" 53 or "attenuation" 52 constant, k x=y=z 2 R is the Cartesian components of the wavenumber, and / x=y=z 2 R are the phase components of the spatial modes.If the BCs are satisfied, then where d x=y=z 2 Z and a uniform b is assumed for all walls.
For the chosen MS to satisfy both Eqs.(10a) and (10b), only some combinations of L x=y=z ; / x=y=z and d x=y=z are possible.By replacement in the wave Eq. ( 2), the MS from Eq. ( 9) yields an analytic forcing term Boundaries and interior test case.Here, the initial conditions are set for all grid points based on Eq. ( 9), while the discrete forcing terms are calculated based on Eqs. ( 4) and (11).
Interior test case.The code verification program is also error prone especially for the MS, for example, if any of the continuous functions are evaluated at an incorrect time sample, the error becomes proportional to DT, yielding a first-order observed accuracy.Consequently, the code verification is generally divided into more manageable tests: a second-order accuracy is expected for the MS if the code update is not exercised near the boundary voxels.To achieve this, a Dirichlet BC is usually used 51 at the edges of the domain-here, a hard source was employed instead of the soft source at each boundary voxel m b surrounding the interior domain.Such source imposes a pressure field which is evaluated from the MS p MS ðM b ; t n Þ.For this scenario, the l 2 DX -norm in Eq. ( 7) is restricted to the DX Á ½1; N x À 1Þ Â ½1; N y À 1Þ Â ½1; N z À 1Þ volume.

Simulations
To further simplify the problem, a cubic room is considered (i.e., L x=y=z ¼ L) together with d x=y=z ¼ d for both the manufactured and ESs.In addition, / x=y=z ¼ p=4 is chosen for the MS.
To keep the observed error within an asymptotic range and within reasonable values relative to single precision, d is fixed at 1, k C ¼ 0:5; In addition, the grid is halved (at constant k C and c) starting from DX max ¼ 0:16 [m] yielding the coarsest interior grid having N x=y=z ¼ 8 voxels per Cartesian dimension.The grid spacing used in the refinement study is: DX 2 f16; 8; 4; 2; 1g [cm].The cubic room mesh was specifically designed for the used grids with double walls such that the voxelized domain D remains consistent between grids.
The grid halving is also convenient since the error can be evaluated at exactly the same time t sim : the number of time samples N t doubles by each halving of the grid for fixed k c and c.Thus, no temporal interpolation is needed.Initial conditions excite the code-verification simulations which are initialized at t ¼ 0 and t ¼ DT from the formal solution and ran until a fixed time t sim % 6:59 [ms].t sim is chosen such that it is divided by the coarsest sampling time t sim ¼ N t DT max with N t 2 N Ã and such that it is of similar order with PRTF simulation times (see Sec. III B 5).The final number of time steps including the two initial conditions is N t 2 f28; 56; 112; 224; 448g [steps], since the initial conditions were given as input, N t À 1 steps are actually computed for each grid.
Large-scale and distributed simulations could present additional challenges.In the present framework, a reasonable Message Passing Interface (MPI) code verification should satisfy the following aspects: (i) each type of partition (e.g., bottom, middle, top) in the MPI-partitioning configuration will be part of the simulation and (ii) the FDTD update CUDA kernel is exercised in each domain and subdomain.Here, the former condition is satisfied when at least three MPI nodes are entailed in the code verification procedure; for the latter, each subdomain must contain at least two z-slices (due to the implementation).
For large values of DX, not all simulations can be computed utilizing all the available devices per MPI node.As such, for each node the number of utilized devices is lowered prior to the domain partitioning such that a subdomain of two z-slices is assigned to each of the utilized devices.For example, each utilized node contained 16 available GPU devices which meant the following number of devices per node was used: f2; 3; 7; 13; 16g [nodes] corresponding to DX 2 f16; 8; 4; 2; 1g [cm].
The code verification program has four main inputs: (1) a list of grid spacings; (2) the number of MPI nodes; (3) the number of utilized GPU devices per node; (4) a triangulated mesh geometry that represents the described cubic room which is further voxelized for each grid spacing DX i .
Consequently, the MPI domain and subdomain partitioning together with the voxelization process are also tested to some small extent during code verification.

Code verification results and discussion
The observed errors are shown in Figs.1(a) and 1(b).Table II displays q obs from an ordinary least squares fit on the error curves in Fig. 1 based on the kk 2 ¼ p DX q obs regression model-linearized using the logarithm function.
As expected, the results for the ES are very close to second-order.A similar, yet not so clear, trend is seen in Fig. 1(a) for the MS and hard sources on the edges of the domain (second test case in Sec.III A 2) which provides evidence for the correctness of the code verification procedure.When the boundary update is also tested for the MS scenario (first test case in Sec.III A 2), the observed order of accuracy generally drops to first-order.
Figure 1(b) also shows that the discretization error introduced at the domain boundary can be weak enough at t sim : for small values of b and certain grid spacings, the secondorder principal error term from the interior update dominates the global error and the observed order of accuracy is 2. As such (in this case, smooth) transitions between different sources of error which have different smooth asymptotic behavior could be observed at some t sim and for a range of DX, as one source of error starts dominating over the others.8)] while the other lines are for the MS in Eq. ( 9).The slopes give the convergence rates: the solid black and gray lines are shown as visual anchors for a first-and secondorder error convergence, respectively.(a) Interior and (b) boundary.TABLE II.Ordinary least squares fit on the error curves in Fig. 1. b represents the specific acoustic admittance, while R 2 is the coefficient of determination rounded to three decimal places.Finally, it can also be seen in Fig. 1(a) that the round-off error starts predominating for DX 2:0 [cm] and certain b values: the convergence rate appears to change to zero-order when the (global) error reaches a value around 0:5 Â 10 À7 , close to the machine epsilon of single precision for P 0 ¼ 1.
The proposed code verification has some limitations.Ideally, the implementation of the voxelization process should be assessed independently of the discretization error through relevant @D-related convergence studies.The authors are unaware of rigorous procedures to do so.Moreover, the present code verification exercise did not include the interpolation algorithms but their results were cross-checked against default software packages.Finally, the calculation of the l 2 DX norms is based on the full interiorpressure domain which is retrieved from the GPU memory at the end of each simulation.Consequently, the code for individual receivers was not considered.Nevertheless, such code involves simple memory transfers at each time step.
In summary, our results show that the utilized code behaves as expected and could be used for a posteriori error quantification.Although code verification cannot be complete, 23 the current results are considered satisfactory because the implementation of the FDTD update code, the initial conditions, two types of sources, the domain partitioning, and the MPI communication were assessed.

B. Solution verification of PRTF simulations
Solution verification is concerned with the assessment of numerical reliability of a computed solution: since asymptotic solutions are neither known for general problems nor can easily be computed in practice (see Sec. II D), the accuracy and/or precision of the result relative to the formal solution needs to be quantified.
Solution verification usually entails an a posteriori quantification of the dominant error for a particular computation.Most common in PDE problems, a technique based on Richardson extrapolation is used (Ref.20, p. 284).This technique is reliably applied when the computed solution is in the asymptotic range.Thus, the asymptotic range of the computed PRTFs needs to be assessed such that the proper error estimate is used.

Source/receiver interpolation
To avoid location-induced errors (see Sec. II D), for each grid DX i , second-order accurate trilinear interpolation is used at each simulated temporal sample and calculated using a tensor product of three linear Lagrange univariate interpolants. 54ubsequent to the receiver interpolation, the principle of reciprocity is re-invoked: each simulation can now be viewed as having the source at an interpolated location and the receiver at the center of an air voxel r s .Running multiple simulations with different source locations displaced relative to r s by m DX with m 2 Z 3 , the source location can thus be interpolated to a desired continuous location.
To facilitate the described interpolation and to minimize the number of simulations, the continuous interpolated source location r ear was chosen as the center of the source voxel on the grid on which the source was furthest from the input surface @D. Farther was quantified in the outward direction on an imaginary "inter-aural" axis.
Due to the frequency-domain behavior of the used interpolant, 55 the monotonic increase with frequency in discretization error is expected to be preserved.

3D pinna and PRTF simulation domain
The employed pinna mesh is a triangulated scan of a cast done on an otologically normal living human.The focus of this study is the blocked-meatus PRTFs.Consequently, an ear-blocking was designed to rigidly fit inside the ear canal and encapsulate a miniature microphone.In the present simulations, the microphone is replaced by a cylinder.
The origin O of the PRTF coordinate system was chosen along the axis of an M8 nut, in the back plane of the nut which is considered a sagittal plane (see Fig. 2 for more details).The central axis of the bolt was defined to be parallel to the interaural y axis.The up axis was chosen as the Cartesian z-axis yielding the x axis as the front direction.
For the ear simulation P ear in Eq. ( 1), the pinna, the earblocking, and the nut meshes were placed inside a domain bounding-box which was chosen such that no reflections would reach the receiver locations during a considered time window of 6.25 [ms].The corresponding direction-of-arrivals are given as (h,/) pairs with h 2 ½À90 ; 90 and / 2 ½À90 ; 90 representing azimuth and elevation angles, respectively (see Fig. 2).The reference simulation to estimate P ref in Eq. ( 1) only contained the bounding-box and the source was placed and interpolated at a location close to the origin of the pinna coordinate system-see Sec.III B 1.
In Secs.III B 3 and III B 4, the PRTFs are computed using the principle of reciprocity (Ref.25, p. 197) with the source at the blocked meatus location r ear and at r ref O with the ear absent.As such, the receivers will be placed at a point in the far-field r 1 relative to r ref.

Discretization error analysis in the frequency domain
Although not a new domain for error analysis, 56 a relevant and advantageous approach to the PRTF problem is to evaluate the convergence and the discretization error in the frequency domain: as formally defined, the PRTFs effectively become steady-state in the frequency domain after a minimal simulated time.
There are several advantages to studying the PRTF discretization error in the frequency domain.First of all, the PRTF problem becomes insensitive to the bandwidth of the initial conditions (see Sec. II D): the norm of the initial data will cancel out due to the pressure division in Eq. ( 1)-assuming a sampled excitation signal and a linear grid transfer function.Consequently, the used source would not require volumetric scaling of the amplitude for a monopole 29,57 or spatial distribution for a volumetric source (Ref.19, p. 173).Moreover, this also eliminates a type of input error commonly referred to as system excitation (Ref.20, p. 558) for a potential validation study.Similarly, the source/receiver interpolationinduced errors could cancel out during the frequency division.In addition, using only a subset of the Fourier coefficients effectively smooths the solution and initial conditions.
Employing the discrete Fourier transform (DFT) sampling of the DTFT spectrum, Eq. ( 6) can be written as e P r 1 x; DT; DX The simulated PRTF can be expressed in the frequency domain using Eqs.( 1), (12), and a geometric series expansion 58 as where HOT ¼ OðDX qþ2 Þ while C T ½x and C X ½x(þC 0 X ½x) are the new temporal and spatial principal error terms, respectively.
Here, "h-extrapolation" 35 is used which involves the quantification of the discretization error employing successive grid refinements using a grid refinement ratio r defined (Ref.19, p. 109) for structured grids as where f si is the sampling frequency for grid i.A minimum of r !1:1 is recommended in practice (Ref.19, p. 124) for a simple grid topology (Ref.20, p. 311).
In the absence of round-off divergence, the error of lowest order will dominate as DX !0, yielding an observed order of accuracy q obs in Eq. ( 13) with a corresponding principal error term C PRTF , where the Courant grid ratio k C is used to exchange DT to DX.
Evaluating the observed convergence rate q obs can be a good indicator to whether asymptotic range is attained: asymptotic range can be considered if q obs ¼ minðg; qÞ, with a typical tolerance 59 of, e.g., 610% (Ref.20, p. 326).
Based on the point Fourier coefficients, the observed order of accuracy can be determined as where g PRTF i ½x represents the simulated PRTF on grid DX i with DX 1 < DX 2 < DX 3 and consistent grid refinements r ¼ r 12 ¼ r 23 .
To study the convergence based on the PRTF magnitude, a power series expansion of f ðxÞ ¼ ffiffiffiffiffiffiffiffiffiffi ffi 1 þ x p can be used together with Eq. ( 15), (17) Ignoring the HOT in Eq. ( 17) yields Equation ( 17) is also valid on the decibel scale-since results were similar, they are omitted in the present work.Equation ( 18) has the advantage over Eq. (16) in that it is insensible to grid-dependent phase errors for fixed x.
In practice, the derivation in Eqs. ( 12)-( 16) tacitly assumes that the calculated discrete Fourier coefficients converge to the continuous ones at least as fast as q obs.This is true, within machine epsilon, 60 for the usual midpoint/trapezoidal calculation of the DFT coefficients on equispaced temporal grids and band limited continuous functions sampled below Nyquist (see, e.g., Corollaries 2.3 and 5.3 in the work by Trefethen and Weideman 61 ).
To estimate the discretization error, a minimal of two separate grids i and j are usually used to compute an asymptotic error estimate E ij of order O(DX qþ1 ) (see, e.g., Chap. 5 in the work by Roache 19 ).Since E ij is approximate, a safety 19,62 factor acting as a "coverage factor" C f !1.0 is usually introduced to conservatively scale E ij , e.g., the grid convergence index (GCI) (Ref.19, pp.112-114).The C f notation is chosen to avoid confusion with the sampling frequency f s .
For example, an uncertainty band for the error of interest E dB ½x ¼ 20 log 10 j g PRTF i ½xj À 20 log 10 jPRTF½xj can be obtained using the power rule (Ref.63, p. 18) in relative error estimation together with the monotonicity of the logarithm function to include the commonly-used GCI, where q is usually based on the observed convergence of the magnitude in Eq. ( 18) if some acceptable 59 asymptotic range is observed.

Regression models in convergence
In complex practical applications, the estimated q obs values could lose precision and show the so-called scatter in the data. 64Potential causes of scatter in q obs include: error cancellation, 20 the round-off error on finer grids, 64 the grid 20,64 and boundary 18 modelling errors, and the actual function of the continuous solutions. 65,66To account for such behavior, Ec ¸a and Hoekstra 64 introduced a least squares analysis of convergence.If the regression model is used for prediction, it is important to pick the correct model.From a regression perspective, both a priori (e.g., knowledge of error behavior) and a posteriori (e.g., hypothesis testing on model parameters) information can be used to correctly choose the regression model.
Based on the error analysis in Sec.II D, the lowest convergence rate is first-order due to the voxelization process.However, if any PRTF feature is linked to volumes of the pinna cavities, the third-order volumetric convergence of the voxelization process will be seen as second-order due to the update scheme.Finally, first-order convergence can also be observed if any non-rigid BC in the computation domain is close enough to the receivers (such that, e.g., the firstorder error propagates to the receiver without strong 1/jrj decay from the propagation on the continuous solutions).Similar to Ec ¸a and Hoekstra, 64 the following linear regression models based on magnitude are considered: • first-order asymptotic model, • second-order convergence, • combined first-and second-order type, where jPRTF r 1 ½xj represents the sought asymptotic magnitude at location r 1 .In addition, Ec ¸a and Hoekstra 64 also proposed the following non-linear regression model: which can also be used to quantify q obs from the leastsquares fit of the parameter q.There are different ways to estimate q in Eq. ( 21), e.g., using non-linear regression, solving non-linear equations, 64 or through the Box-Tidwell 67 procedure.
The above models assume that DX is not too small such that round-off divergence is observed in which case error/asymptotic estimation becomes more difficult.
By utilizing the regression models proposed by Ec ¸a and Hoekstra, 64 the asymptotic solution PRTF r 1 ½x becomes a random variable as a parameter of a probabilistic model.As such, one can quantify the uncertainty of the asymptotic solution through the derivation of confidence intervals (CIs) as long as the probabilistic models are used appropriately.The asymptotic solution can be predicted from the regression models based on the mean output of each model as where E is the expected value operator.In practice, jPRTF r 1 ½xj in Eq. ( 22) is considered the intercept yielded by the least-squares fit of the regression models.).Although (i) is satisfied due to the uniform and structured grid used, it is difficult to give strong arguments for (ii)-(iv) given the problem at hand.The usual visual techniques to assess normality are difficult to apply to the PRTF problem due to the large amount of data.One alternative is to use statistical tests of normality, which might not be very reliable for a small number of grids N g .
For each frequency, homoscedasticity (iv) is likely violated. 18As such, the weights introduced by Ec ¸a and Hoekstra 64 for weighted least squares models could be helpful.
Non-parametric methods-The alternative to the Gaussian-based CIs is to use non-parametric estimation methods.One such technique is bootstrapping 70 which, in its basic form, only relies on the residual independence (ii) assumption.Such assumption may still be difficult to motivate for minute changes in DX.In the context of regression, two main bootstrapping methods exist (Ref.68, p. 494): bootstrapping pairs and bootstrapping residuals.Bootstrapping pairs is generally used whenever there is increased doubt about the assumed regression model.There are also different methods to estimate the CI: the bias-corrected and accelerated (BC a ) method is chosen which corrects for bias and skewness in the bootstrap samples and is "second-order" accurate [i.e., estimate error is proportional with O(1/N g )]-see Chap.14 from Efron and Tibshirani. 70uch uncertainty estimates are different than the estimates proposed by Ec ¸a and Hoekstra 64 based on the safety factor C f .Nevertheless, if required, C f can be included in the analysis as an appropriate coverage factor.

PRTF simulations
A total of 10 directions are simulated (see Fig. 2) on N g ¼ 6 grids and a constant grid refinement ratio of with i 2 f1; 2; …; N g À 1g.To directly apply Eqs. ( 16)-( 19), the same discrete frequencies f k ¼ x k ð2pÞ À1 must be sampled with the DFT on each grid.This implies a consistent frequency resolution Df ¼ 1=t obs , where t obs is the observation interval where N i 2 N Ã .To obtain consistency of t obs , the sampling frequency f si on grid i was adjusted incrementally from an initial value that satisfied the grid refinement ratio r i;iþ1 ¼ 1:1 in Eq. ( 14) until Eq. ( 25) was satisfied for t obs ¼ 6:25 [ms].Since t obs and k C were fixed while DX changed given an integer sampling frequency f si 2 N Ã enforced by the solver, it is not always possible to perfectly preserve both r ij and Df.
The chosen frequency resolution is Df ¼ 160 [Hz].A sub-millimeter voxel is desirable for increased precision of the calculated FDTD PRTFs at higher frequencies. 18The smallest voxel is limited by the available resources.The voxel sizes used vary from 0.95 to 0.59 [mm]-Table III shows the details of the convergence study.
The surface boundary impedance for the pinna, earcanal blocking, and surface of the microphone was set as fully rigid b ¼ 0. For all simulations, a temperature of 294.78 [K] was used which entailed a dry-air sound speed (Ref.71, p. 121) of about 344.34 [m/s].

PRTF results and discussion
To avoid ripple distortions in the PRTF spectra due to "edge effects" between DFT periods, a continuous Gaussian window was sampled at t n ¼ nDT i (n 2 N) and applied to both the interpolated simulated temporal signals in Eq. ( 1), where t 0 ¼ t obs =2 and the parameter a was chosen such that the amplitude of the window had negligible values at the beginning and end of the observation interval: wðt obs Þ % 1 Â 10 À7 .Due to the frequency division in Eq. ( 1), the DT scaling in the DFT cancels out and the PRTFs were computed with the usual fast Fourier transform (FFT) routines which are working on normalized temporal grids (Ref.72, p. 3).A discrete delta sequence was used as the input to the driving soft source for both simulations in Eq. ( 1).No circular shifts were applied to the computed g PRTF.In the following, unless otherwise mentioned, the results are presented for the interpolated PRTFs (both source and receiver-see Sec.III B 1): interpolation was done on the center of the source voxel from grid 5.The analysis is limited to 24 [kHz] yielding a total of N freqs ¼ 151 positive frequency bins.
Assessing asymptotic range.Since an a posteriori error estimation in Eq. ( 19) may be problematic if the asymptotic range is not attained, 35 a first analysis of the computation of the observed convergence rate q obs is performed.Q obs can be computed either based on three grids [Eqs.( 16) and ( 18)] or estimated based on the non-linear regression model in Eq. ( 21).Results for one typical direction are presented in Figs. 3 and 4. Results show similar behavior for the other directions.Undefined results in R such as the logarithm of a negative number are excluded from the plots.
The three-grid estimation of q obs is shown in Fig. 3.In this case, since 6 grids are available for each frequency bin, q obs can be estimated 4 times based on three consecutive grids (i.e., fj; j þ 1; j þ 2g; j ¼ 1…4 grid triplets) using r ¼ e r and twice on every other grid (i.e., f1; 3; 5g; f2; 4; 6g grid triplets) using r ¼ e r 2 in Eqs. ( 16) and ( 18).In the calculation of q obs in Eqs. ( 16) and ( 18), the ideal value of e r ¼ 1:1 was used and not the empirical values in Table III-results TABLE III.Grids used in the PRTF convergence study.Interpolation done on source/receiver location from grid 5.The voxel sizes DX and N voxels are rounded to two decimal places for their units.N voxels represents the total number of voxels used in simulations excluding the halos on each MPI node.N Nodes represents the number of MPI nodes-see Appendix B. r i;iþ1 also contains the deviation from the chosen e r in Eq. (24).N DFT represents the number of temporal samples used to compute the DFT (i.e., the simulation interval in samples) such that a Df frequency resolution is attained in the DFT.N freqs represents the number of positive frequencies used in the analysis, while N dirs represents the total number of PRTF directions.Grid  3. (Color online) Estimated observed order of accuracy q obs for direction ðÀ90 ; 0 Þ based on magnitude [Eq.( 18)] and on complex-values [Eq.( 16)].q obs ½e r uses three consecutive grids while q obs ½e r 2 uses three alternate grids.Horizontal dashed lines show expected first-and second-order q obs embedded in 610[%] intervals.
were almost identical when q obs was numerically determined 35 for r 12 6 ¼ r 23 .Results show a high amount of scatter with unreasonable values of q obs , for instance, the divergence is worse than OðDX À10 Þ for some frequency bins.Data for PRTFs directions close to the dispersion-free grid diagonal [i.e., ðÀ45 ; 40:6 Þ and ð45 ; 40:6 Þ] did not exhibit improved results.Moreover, the discrimination ratio 64 does not show any consistent bandwidth of monotonic convergence.
It is unclear why such unreliable results are obtained for the three-grid estimates of q obs .Lower frequencies might not show convergence due to the round-off noise floor: for instance, the 160 Hz frequency bin is sampled more than 2000 times in a wavelength.The grid size DX might still be too large for higher frequencies where solution could be affected by pollution error even at the current voxel sizes employed, poor precision in the PRTF data due to the voxelization process, and small refinement ratios r.
Too small r are expected to cause too little changes in the computed solution (Ref.20, p. 331) and would render the obtained data unreliable for a convergence study.Figure 3 shows that computing q obs with the larger grid refinement ratio of e r 2 ¼ 1:21 only provides marginal improvements: the range of values for q obs diminishes but not to a useful extent.Thus, a much larger refinement ratio r might be needed which conflicts with the limitations in the methods: memory limitations for lower DX and increased solution scattering 18 for larger DX.Finally, a clear convergence trend is seen in the PRTF magnitudes as DX is decreased: see the convergence of the "concha depth" 73 mode around 5 [kHz] in Fig. 5.
Thus, estimating q obs does not seem to be wellconditioned for the current problem and models: the random error in the voxelization process coupled with numerical sensitivity might render the three-grid method unreliable.Since solutions depend on the boundary, such voxelizationinduced "boundary noise" could translate into slight changes in the asymptotic solutions.The expected strong sensitivity of the PRTFs to the pinna shape further aggravates such changes.Thus, it is possible that computed solutions from different grids belong to neighboring analytic solutions 47 and render the estimation of q obs unusable at a single point in space.Spatially averaged solutions might be less sensitive to such small boundary changes.
Consider now the estimation of q obs based on the nonlinear regression model.Figure 4 shows the estimated q obs based on Eq. ( 21).Results are shown for both weighted [see Eq. ( 23)] and non-weighted regression models.Compared to Fig. 3, results show less erratic results but are still unreliable: many frequencies far away from the direct current (DC) one (i.e., x ¼ 0) show divergence with other frequencies showing unreasonably high convergence rates.For many frequency bins, results are similar for different methods of calculation for the non-weighted models.Similar to the last example in the work by Ec ¸a and Hoekstra, one reason for such results could be the small number of grids. 64Non-linear regression was also applied to Eq. ( 21) but results were too sensitive to the starting values which were poorly estimated based on the three-grid calculations-see Fig. 3.
Based on the results in Figs. 3 and 4, it would be difficult to estimate a bandwidth where asymptotic range is attained.Thus, the commonly used two-grid error estimates such as Eq. ( 19) will yield unreliable results and cannot be used on the current simulated data with a high degree of confidence.
Estimating the asymptotic solution.Results of the predicted asymptotic solution by the linear regression models are shown in Figs. 5 and 6 for the ðÀ90 ; 0 Þ direction.Results are similar for the other directions.
CI estimation-Notice first that although the plots in Figs. 5 and 6 are on the decibel scale, the uncertainty in the CI is not propagated through the logarithm function [see, e.g., Baird (Ref.63, p. 19)].Figure 5 shows the behavior of the three types of estimates for the CIs for two regression models: first-order and weighted second-order.A number of  21) and different calculation methods for direction ðÀ90 ; 0 Þ. y-axis is truncated to ½À10; 10. Green dashed lines show expected first-and second-order q obs embedded in 610[%] intervals.
B ¼ 5000 bootstrap replicates are used for the non-parametric estimates. 74Results show that the two-tailed parametric CI estimates are usually within the nonparametric ones.Moreover, the tightest estimates are found by employing BC a on bootstrapping residuals while the largest estimates are expressed through the BC a estimates on bootstrapping pairs.Such results on the CIs generally extrapolate to all the considered linear regression models: no strong qualitative deviations were observed.Since no normality of residuals is required for the non-parametric CIs, they are to be preferred to the parametric CIs.If more conservative estimates of region of the asymptotic solution are preferred, CIs are better estimated with bootstrapping pairs.
Asymptotic predictions-The asymptotic predictions of the linear regression models in Eqs.(20a)-(20c) are shown in Fig. 6(a) together with the most conservative CI estimates.A first result is that the first-and second-order model in Eq. (20c) shows poor results: predictions contain PRTF features which are not present in any computed solution-see Fig. 5.Moreover, the model predicts unphysical discontinuities or negative magnitudes.Unreasonably large CIs also suggest a poor fit.Finally, a posteriori significance tests for the two regressors DX and DX 2 from model (20c) showed significance in less than 2% of cases for both weighted and non-weighted models.The significance was such that, e.g., H 0 : C 00 12 ¼ 0 was rejected (Ref.68, p. 84) for a p-value 0:005 with the data pooled over all frequency bins and directions.Such results are considered poor enough to show the ill-fit of such regression model.Consequently, the first-and second-order model in Eq. (20c) is excluded from further However, it is not so to exclude one of the remaining two linear regression models: both models predict similar asymptotic solutions which are consistent with the computed solutions on the six utilized grids.One way to assess the fit of the regression models is to explore the coefficient of determination R 2 .For all the N freqs Â N dirs ¼ 1510 data points, R 2 for the first-order model in Eq. ( 20a) is larger than for the second-order model in Eq. (20b) in 49.21% cases.For the weighted linear models, 50.73% of the data points indicate a better fit for the firstorder model based on R 2 .Thus, neither of the two models can be excluded based on aggregated results.A frequencydependent analysis of R 2 shows that the first-order model generally fits the data better up to about 8.5 [kHz] and other sporadic frequency bands such as ½11:5 À 12:5 [kHz] and ½19 À 20 [kHz].Such results are seen for both weighted and non-weighted linear models-results not shown.
In regression analysis, model selection strategies based on "domain knowledge" (Ref.68, pp.285-286) are favored to descriptive data-driven selection.The non-interpolated PRTFs can be used as a model selection criterion: the PRTFs will converge in a first-order fashion due to source/receiver inconsistencies (see Sec. II D) assuming other errors such as the pollution error are not dominating.The source/receiver location for the non-interpolated PRTFs are calculated in a nearest-neighbor fashion, at a voxel close to the continuous location given by grid 5.
Figures 6(a) and 6(b) both contain the first-order prediction of the asymptotic solution for the non-interpolated PRTFs together with the pair-bootstrapped CIs.The firstorder regression model in Eq. (20a) was used for the noninterpolated PRTFs.It can be seen that for the concha depth mode, the second-order model predicts an inconsistent solution: the CIs are disjoint.This result is consistent with the frequency-dependent analysis of the coefficient of determination R 2 .Results suggest that, at low frequencies, the simulations are dominated by first-order errors originating in the voxelization process.
For the other PRTF directions in the study, in addition to the general misprediction of the second-order model of the asymptotic amplitude of the concha depth mode, its CIs do not overlap with the non-interpolated ones at various 1-2 kHz wide frequency bands (results not shown).Moreover, although not statistically significant, the predicted asymptotic solution from Eq. ( 22) is more consistent with the first-order model than the second-order one: one could see small deviations in the peaks and notches of the PRTF spectra relative to the non-interpolated PRTF.
It should be mentioned that: (i) for direction ð45 ; 0 Þ the CIs overlap for the three regression models: first-order, second-order, and non-interpolated first-order; (ii) compared to the second-order regression model, the first-order model fails to predict one PRTF notch for direction ð90 ; 40:6 Þ: the asymptotic value of the magnitude is negative.
Figure 6(b) shows the weighted results.Results are highly similar with the non-weighted ones in Fig. 6(a).Both the asymptotic solutions and the CIs in the weighted case undergo small corrections: the asymptotic solutions change mostly less than 0.5 [dB] for the first-and second-order models while the CIs are generally smaller for the weighted models.The 0.5 [dB] value is estimated based on a histogram of the difference for the data pooled across all frequencies and directions-histogram not shown.All the previous results from non-interpolated regression models are found adequate also for the weighted regression models.Since the precision of computed PRTFs increases with smaller voxel, 18 the weighted models are preferred over non-weighted models and the small corrections are considered improvements in the accuracy of predictions.
Interpolation impact-The effect of interpolation is also seen: using interpolated PRTFs improves the precision of the estimated asymptotic solutions.Although not shown, the largest effect on the precision is for interpolation within the concha volume (here, source interpolation) as compared to interpolation in the far-field.Moreover, using the weights in Eq. ( 23) appears to partly correct the prediction bias in the asymptotic solution in Eq. ( 22) of the non-interpolated PRTFs relative to the interpolated PRTFs.Consequently, since the boundary @D gives the lowest order of convergence for the employed discrete models, a weighted least-squares could be used without the need to interpolate the solution if lower precision is acceptable.
Peak/Notch estimation-Of interest for the PRTF problem are the frequency locations of the peaks and notches and their convergence.Such PRTF features were coarsely extracted in a nearest-neighbor fashion and shown in Fig. 7 for the same ðÀ90 ; 0 Þ direction.Here, only the first-order model is considered.Despite the OðDf Þ error in frequency estimation, it can be seen that such features undergo corrections of various degrees.The amount of correction was also found to be direction-dependent.Moreover, the individual computations could display spurious features compared to the asymptotic solution, mostly below 20 [kHz].Similarly, the asymptotic solution could show extra features above 20 [kHz]-nevertheless, such conclusion is weak since the uncertainty is high at such frequencies.It is difficult to say based on the present results whether any PRTF feature converges non-linearly in DX which might indicate a volumetric dependence.

IV. DISCUSSION
Results in Figs. 3 and 4 show the difficulty in estimating the computational errors in practical problems with complex geometries: asymptotic range does not seem to be attained even on the submillimeter grids used.Depending on frequency, various sources of error and their interaction could be responsible for such behavior-see Secs.III B 6 and II D. Given the present models, one solution would be to increase the computational accuracy: using double precision for lower frequencies and smaller DX such that asymptotic range is reached within some bandwidth.However, such solution is presently unreachable due to restrictions in computing memory: the simulation on grid 1 already uses more than 100 Â 10 9 voxels.
Alternatively, regression models will not only predict the asymptotic solution, but will also reflect the quality of such predictions: larger CIs of the asymptotic solution will signal stronger departures from the expected asymptotic regime-equivalently, a reliability decrease of the computed solution.However, to avoid both prediction bias and incorrect CIs, the asymptotic behavior of the errors present in the computations needs to be understood such that the correct regression model is applied.
Since the mathematical theory predicts that the discretization error is a bias in the solution, statistical estimates for computed solutions are sometimes discredited. 23,75evertheless, the discontinuity in the convergence of the boundary @D (which may cause, e.g., jumps in the error asymptotes 47 as DX !0) motivates the random treatment of the convergence.Moreover, the commonly used techniques in error estimation were shown to be difficult to apply for the present problem and associated parameters.The estimated asymptotic solutions also show that a two-grid error estimate such as the GCI would have been unreasonably large for the PRTF/HRTF problem.Additionally, the asymptotic solution was estimated based on the asymptotic behavior of the discretization error of lowest order which should persist, as DX !0, until round-off error will cause divergence from the true solution.This also motivates predictions using least squares outside the observed grid interval: roundoff divergence would dominate the results as DX approaches 0. Finally, no noticeable divergence is qualitatively observed at low frequencies.
Nevertheless, such regression models are not perfect.Potential caveats in faulty model assumptions were mentioned in Sec.III B 4: unless the independence assumption in the residuals is strongly violated (see, e.g., Sec.II D), the given estimates and CIs are expected to have a reasonable degree of reliability.III and the asymptotic solution for the first-order model in Eq. (20a).Marker size is proportional to the log of the prominence in dB.
Although the fit at higher frequencies is not much better than a second-order model, the first-order regression model is still considered representative for the used discrete models since: (i) it better matches the non-interpolated predictions, (ii) there is no previous data supporting volumetric ear modes (most previous studies only consider distances 17,76 ), and (iii) it accounts for the smallest sources of error (see Sec. II D).Besides, larger CIs could signal a poor fit.
An alternative would be to use different regression models depending on frequency but this could cause discontinuous asymptotic solutions.Since pollution error could be of concern at higher frequencies and in the spirit of more conservative estimates, the CIs calculated with the nonparametric bootstrapping BC a -pairs are advised.
The present paper only deals with the quality of finitedifference PRTF calculations and not with their validity in predicting real-world PRTFs.Such model legitimacy needs to be assessed in a separate validation study where FDTD predictions along with their estimated computational errors are compared with real-world measurements.

V. CONCLUSION
Individual computations are not necessarily a faithful representation of wave-based formal solutions when complex geometries are involved: they could embed errors which should not be assumed negligible.The present work showed that the current state in HRTF/PRTF simulation literature which involve a single computation could yield erroneous results: for each input parameter such as b, single-grid computed solutions could be biased relative to the formal solution.Reliability assessment of the computational result should always be performed.
The employed FDTD solver was verified using relevant methods from the V&V literature-a mandatory step prior to the investigation of the error in the simulated HRTFs/PRTFs.
The error in computed PRTFs for the ipsilateral ear was assessed for a blocked-meatus location through a convergence study on sub-millimeter grids.Since no clear asymptotic range was observed, the commonly employed techniques based on Richardson extrapolation were found inapplicable for the used grids and models.Estimates of the asymptotic solution using weighted 64 regression models fitted on multiple grids were found more adequate.The method also provides CIs for the asymptotic solutions which could provide reliability information for any validation study.The non-parametric bias-corrected and accelerated CI estimation method based on bootstrapping pairs provided more conservative estimates with better justified assumptions.
The estimated asymptotic solution was found less reliable at high frequencies and at notches in the PRTF features even for the unprecedented grid spacing used.Moreover, the proposed method requires a thorough understanding of the sources of error in order to employ the correct regression model: prediction bias is otherwise expected.In the present case, a first-order regression model was found adequate.It was also found that, due to the stair-stepped boundaries, source and receiver interpolation does not appear to be mandatory.However, interpolation improves the precision of the estimated asymptotic solution.
Evaluation of the magnitude of the computational errors or asymptotic estimation of the solution such as given in the present work are a prerequisite for a correct assessment of modeling errors and assumptions in wave-based HRTF simulations.Due to complex interactions among numerical errors even for a simple wave model such as the SRL scheme, using HRTF/PRTF simulations as a reliable blackbox prediction tool seems presently unfounded: no singlegrid computation is necessarily correct in the entire audible range even when employing submillimeter voxels; error assessment is recommended.The proposed procedure is one such evaluation method.
Parentheses are used for arguments in the case of continuous signals, while square brackets are used for discrete signals.X is used for the angular frequency [rad/s] for continuous temporal signals, while x is for the angular frequency of discrete-time signals.The set of natural numbers is denoted by N (with N Ã :¼ N À f0g), integer numbers by Z and R denotes the set of real numbers.
Contrary to the notation in the PDE literature, independent variables as subscripts to dependent variables are used as naming conventions: p t 6 ¼ @p=@t while @ t p ¼ @p=@t.

v b ðtÞ À b b P j ðtÞ ¼ b b DX 2 @
p b ðtÞ @x þ b b Á HOTðtÞ;

FIG. 1 .
FIG. 1. (Color online) Observed global spatial l 2 DX error norm jjjj 2 at t ¼ 6.59 [ms].L ¼ 1.28 [m], P 0 ¼ 1; k C ¼ 0:5, [1,1,1] mode, 3 MPI nodes.The line for a zero specific acoustic boundary admittance b correspond to the ES case [i.e., rigid box in Eq. (8)] while the other lines are for the MS in Eq. (9).The slopes give the convergence rates: the solid black and gray lines are shown as visual anchors for a first-and secondorder error convergence, respectively.(a) Interior and (b) boundary.

FIG. 2 .
FIG. 2. (Color online) Used Cartesian and spherical coordinate systems (left) with the computed PRTF directions (right).The directions in the first column are for the left ear of a hypothetical head.

FIG. 7 .
FIG. 7. (Color online) Peaks (top) and notches (bottom) on each grid in TableIIIand the asymptotic solution for the first-order model in Eq. (20a).Marker size is proportional to the log of the prominence in dB.

TABLE I .
Relevant sources of error for the employed discrete models and problem of interest, relative to the formal solution of the PDE model in Eqs.(2) and (3) assuming the triangulated mesh as the true geometry.Results obtained by applying a priori bounds (Ref.36) within a multiplicative separation of variables frame (Ref.34) coupled with the ordering DX Àj ( DX Àðjþ1Þ as DX !0 with j 2 N Ã . a b Due to the frequency division in Eq. (1) and the frequency-analysis for computed PRTFs, temporal interpolation is not considered (N=A) in the present work.See more details in the text, Sec.II D.
OðDT 4 ; DX 2q Þ.For fixed r 1 ; t;r 1 ðtÞ satisfies a differential equation that depends on the update and p r 1 ðtÞ (see, e.g., Secs.III 9 and III 10 in Hairer et al.
Parametric methods-The CIs can be estimated based on the usual covariance of the inputs and variance of the residuals [see, e.g., Sec.3.4 from the work by Montgomery et al. (Ref.68, p. 95)].Such CIs are reliable when the regression model assumptions are satisfied (Ref.68, p. 122 and Ref. 69, p. 28): (i) DX is measured without error, (ii) the residuals of the regression model are independent of each other, and (iii) the residuals are normally distributed with (iv) constant variance.For example, gross violations of normality could have a serious effect on the estimated CIs (Ref.68, p. 139