Comparative evaluation of interpolation methods for the directivity of musical instruments

Measurements of the directivity of acoustic sound sources must be interpolated in almost all cases, either for spatial upsampling to higher resolution representations of the data, for spatial resampling to another sampling grid, or for use in simulations of sound propagation. The performance of different interpolation techniques applied to sparsely sampled directivity measurements depends on the sampling grid used but also on the radiation pattern of the sources themselves. Therefore, we evaluated three established approaches for interpolation from a low-resolution sampling grid using high-resolution measurements of a representative sample of musical instruments as a reference. The smallest global error on average occurs for thin plate pseudo-spline interpolation. For interpolation based on spherical harmonics (SH) decomposition, the SH order and the spatial sampling scheme applied have a strong and difficult to predict influence on the quality of the interpolation. The piece-wise linear, spherical triangular interpolation provides almost as good results as the first-order spline approach, albeit with on average 20 times higher computational effort. Therefore, for spatial interpolation of sparsely sampled directivity measurements of musical instruments, the thin plate pseudo-spline method applied to absolute-valued data is recommended and, if necessary, a subsequent modeling of the phase.


Introduction
The first studies on the specific sound radiation characteristics of the human voice were conducted as early as the late 1930s [1], while systematic investigations of the directivity of musical instruments began 30 years later [2]. The radiation patterns of acoustic sound sources such as speakers, singers or musical instruments are commonly measured in anechoic environments with the source centered in an enclosing spherical microphone array.
For a comprehensive analysis of the directivity of 40 human speakers a nearly full spherical array was used, measured sequentially at 253 positions [3]. With respect to the singing voice, the radiation characteristics of 8 *Correspondence: david.ackermann@tu-berlin.de 1 Audio Communication Group, Technische Universität Berlin, Einsteinufer 17c, 10587 Berlin, Germany Full list of author information is available at the end of the article opera singers [4] and 15 trained singers [5] were determined in the horizontal plane, measured at 9 and 13 positions, respectively. A higher spatial resolution was used for measurements of a professional male singer using an adjustable semi-circular microphone array with 24 receivers [6]. For a recent review of research on the sound radiation of singing voices, see [7].
The directivities of eight musical instruments were measured using 64 microphones [8], 22 microphones were used for a measurement of 22 instruments [9]. A recently generated database for 14 instruments and a speaker contains radiation patterns measured at 2522 positions on a sphere [10]. However, these data contain only (third-)octave band directivities, which limits their use for research purposes. The most comprehensive public database was collected for 41 modern and historic instruments measured with 32 microphones and contains single tones within the playable range of each instrument and directivities computed from the stationary parts of these tones [11,12].
The spatial resolutions of the available directivity measurements of acoustic sound sources thus depend on the technology used and differs greatly from each other. At the same time, many source directivity based applications such as room acoustical simulations require either continuous or higher resolution data. And even if the application uses a discrete spatial representation, the sampling grid required is usually different from that used in the measurement. The measured data must therefore be interpolated or resampled.
In the common polar representation, the measured values are usually linearly interpolated (cf. [9]) and occasionally also smoothed in addition (cf. [5]). For 3D balloon plots showing the spherical radiation pattern for single frequencies or frequency bands, a common method is to decompose the sound pressure measurements into spherical harmonic (SH) basis functions followed by spatial oversampling on the surface of a sphere. The resampled grid is sometimes linearly interpolated at the end for visual display (cf. [8]).
The accuracy of room acoustical simulations was shown to strongly depend on the directivity of the sound source and thus also on the quality of the chosen interpolation method. The angular resolution of the directivity affects the simulated room impulse response (RIR) and several other room acoustic parameters up to a spherical harmonics (SH) order of N = 10 [13] even if the information incorporated in higher order components may no longer be perceptually relevant, at least at larger distances from the source [14].
Not only the required resolution and the required sampling scheme, but also the required physical information inherent in the radiation pattern depends on the subsequent application. In wave-based simulations, such as the Boundary Element Method (BEM) and the Finite Element Method (FEM), it appears beneficial to have both, a continuous magnitude and phase response of the source [15]. In simulations based on geometrical acoustics [16] that combine image sources and stochastic ray tracing to compute early reflections and the late reverberation [17,18], a complex-valued description of the source directivity might be beneficial for the image source part, while the phase response is spurious in an energy-histogram-based ray tracing approach.
Moreover, if directivities are calculated from the steady part of played tones, the phase spectrum may be subject to fluctuations, especially if the source in the center of the measurement system is not completely spatially fixed, causing a fluctuating excess phase that renders phase information practically useless (Fig. 1). To account for this, Zagala & Zotter [19] suggested to iteratively optimize the sign of the absolute magnitude response prior to SH interpolation to minimize the mean squared error (MSE) between the input and interpolated data. Ahrens & Bilbao [20] chose to make the magnitude response minimum phase to avoid excess phase and to get directivity more easily decomposed into SH impulse responses applicable to time domain room acoustical simulations. However, both studies did not investigate the general suitability of SH for interpolating the magnitude response of the directivity.
The question about whether and how to interpolate directivities with phase has been successfully addressed for head-related transfer functions. A variety of techniques either pre-align the entire or a high-frequency portion of the impulse responses, or manipulate the corresponding phase to improve magnitude interpolation (cf. [21] and ( [22], Chapter 4.11) for an overview). These methods either reconstruct the phase after interpolation or are justified by the irrelevance of inter-aural phase at high frequencies [23] and rely on the relation of the frequency-domain directivity to a short, impulse-shaped time-domain representation. Hereby, they do not apply to the directional spectra of musical instruments, as exemplified in Fig. 1.
The aim of this study is thus to evaluate the suitability of established methods for interpolating the magnitude response of sparsely sampled directivities of musical instruments. For this purpose, high-resolution measurements of four different musical instruments, whose technical construction and radiation characteristics cover a wide range of natural sound sources, were selected. The data were sub-sampled at 32 sparse grid points, interpolated to a high-resolution grid, and evaluated against the measured reference. Note that in this paper we use the term "interpolation" for any kind of continuous approximation of discrete spatial radiation patterns, no matter whether the grid points are precisely reproduced by this approximation or not.

Background
A plethora of interpolation techniques for real-valued scattered data exist that make different assumptions about the distribution of the discrete set of known data points [24]. Because the quality of the interpolation depends on how well these assumptions are fulfilled, the performance of the interpolation methods considerably depends on the specific application. Simple techniques include discontinuous nearest-neighbor interpolation, as well as continuous linear and natural neighbor interpolation. More commonly used are advanced concepts such as deterministic inverse distance weighted or spline interpolation [25], as well as kriging [26]-a stochastic technique from the field of geostatistics that minimizes the spatial variance between the value to be estimated and the ambient measurements. An essential tool for data fitting and interpolation in the field of computer aided geometric design (CAGD) are barycentric coordinates defined on spherical triangles, which can be used to define the associated spherical Bernstein-Bézier polynomials for constructing piece-wise functional and parametric surfaces [27]. For acoustical sound sources, a decomposition into SH basis functions has become particularly popular [28][29][30], since it not only allows for a synthesis of the radiation pattern in virtual acoustic reality [31], but also for a decomposition of the room impulse response into SH-based spatial components [32]. In case of an order-limited directivity, SH interpolation is physically correct.
Based on the above review, we selected three interpolation approaches for the detailed evaluation. SH interpolation was included because of its widespread use in musical acoustics. Spline interpolation was chosen because it is superior to inverse distance weighting and kriging if only a small number of sample points are available [33,34]. The spherical triangular interpolation technique corresponds to a piece-wise degree-1 barycentric spherical Bernstein-Bézier polynomial interpolation; in audio technology it is commonly employed in three-dimensional vector based amplitude panning (VBAP) as introduced by Pulkki [35] for robust virtual sound source positioning [22].

Spherical harmonics interpolation
If the sound pressure on the surface of a sphere is sampled with a finite number of microphones, spherical Fourier coefficients can be calculated from the measured values, which can then be used to estimate the sound pressure function on the entire measuring surface [36]. The limited number of sample points results in an order-limited sound pressure function on the measurement surface.
Thus, the spherical function f (θ, φ) (θ = azimuth, φ = colatitude) is represented by a weighted sum of a finite set of orthogonal base functions: where N ∈ N indicates the spherical harmonics order and f nm are the considered weights of the corresponding spherical harmonics where P m n (·) are the associated Legendre functions, (·)! represents the factorial function, m ∈ Z specifies the function degree, and n ∈ N the order of the function. Consequently, the Fourier coefficients f nm completely describe the order-constrained function f (θ, φ) on the entire sphere and their determination is yet sufficient for a correct SH interpolation.
By sampling the sound pressure function f (θ, φ) with a Q channel spherical microphone array, the samples p q = f (θ q , φ q ) are given at the positions (θ q , φ q ) of the respective microphones for q ∈ {1, 2, ..., Q} = N Q . In matrix form Eq. 1 can be written as where the matrix Y of dimensions Q × (N + 1) 2 is given by and the vector f =[ p 1 , . . . , p Q ] T contains the Q sound pressure measurements at position (θ q , φ q ) for q ∈ N Q . For the rare scenario, when the number of microphones Q matches the spherical harmonics order N, i.e. Q = (N+1) 2 , under consideration of perfectly distributed measuring points [37] and thus a well-conditioned full-rank matrix Y, Eq. 3 can be solved with the inverse of matrix Y: For Q > (N + 1) 2 an over-determined system of linear equations results which can be solved through best fit, in the least-squares sense, by taking the Moore-Penrose inverse of Y and thus seeking a solution f nm that minimizes the energy of the error: with Y † = (Y H Y) −1 Y H and · denoting the Euclidean norm. For functions that are not order-limited, errors occur due to spatial aliasing and f = Yf nm and consequently f (θ q , φ q ) = p q [38]. For Q < (N + 1) 2 , the system of equations is underdetermined and Eq. 3 provides infinitely many solutions. In this case the Moore-Penrose inverse of the matrix Y seeks a solution f nm with minimum Euclidean norm, i.e. with minimal wave-spectral power f nm 2 ( [29], p. 79): To interpolate samples of the sound pressure measurements on a sphere, the calculated weights of the spherical harmonics can be used in the inverse spherical Fourier transform from Eq. 1 and arbitrary points between the samples can be estimated. The values at the sampling positions (θ q , φ q ) for q ∈ N Q can be reproduced exactly if the order N is sufficiently high. In the case of underdetermined systems, however, notches occur between the sample points due to the chosen constraint of minimum wave-spectral power and therefore even order-limited functions can no longer be represented accurately.
An indication for the numerical accuracy of SH interpolation based on matrix inversion (Eq. 5) is the condition number κ of Y N . A large condition number indicates that small changes in the measured sound pressures f could lead to large changes in the Fourier coefficient matrix f nm . The solution of the linear system of equations is thus highly sensitive to errors and noise in the input data. While κ = 1 is ideal, a system with κ > 3.5 is considered as ill-conditioned [39]. The condition number depends on the chosen spatial sampling scheme and the SH order N.

Thin plate pseudo-spline interpolation
The thin plate pseudo-spline solution [40,41] allows the regularized interpolation of sparsely distributed measurements on the sphere with closed-form expressions that make this approach well suited for numerical computation. The aim is to find a smooth function f (θ, φ), where the values for f (θ q , φ q ) should be as close as possible to the measured values p q while containing minimum bending energy on the surface of the sphere S. An interpolating (A) or smoothing (B) thin plate pseudo-spline can therefore be obtained by seeking the solution to one of the following problems: for (A) or with the option of regularization for (B), where λ ≥ 0 denotes the tuning parameter and J k (f ) is defined by withf (11) and A solution of the two problems given by Eqs. 8 and 9 is obtained with where P n are the associated Legendre polynomials and z denotes the cosine of the spherical angle γ between the two arguments of the kernel function with The spline order M ∈ N determines the derivability of the solution from Eq. 13. We define the spline order as M = 2k − 2, and corresponding splines are continuous up to the (M − 1)th derivative, so they are called C M−1 smooth.
A closed-form expression for the reproducing kernel R(θ, φ; θ q , φ q ), suitable for numerical computation, is given by with (17) and 2k − 2 = M. A recursive evaluation of q 2k−2 (z) for k = 3 2 , 2, 5 2 , ..., 6 can be found in ( [40,41], Tab. 1), as well as the determination of the coefficients c and d from Eq. 13 in matrix form 1 : where If the measured values are noisy, it can be advantageous to regularize the interpolation in order to suppress outliers; the tuning parameter λ > 0 will smooth the estimated function f on the surface of the sphere. Due to the low noise measurement data used for this study, smoothing of the estimation function did not improve the quality of the interpolation (c.f. Section 5), therefore the thin plate pseudo-splines were performed without regularization (λ = 0).

Piece-wise linear, spherical triangular interpolation
The entire set of Q microphone positions (θ q , φ q ) can be equivalently expressed as a 3 × Q matrix containing its three-dimensional unit direction vectors Using the Quickhull algorithm [42] vertex index triplets v l =[ v 1l , v 2l , v 3l ] are obtained to describe a set of triangular facets that span the convex hull of the vertices stored in U. Any arbitrary unit direction vector u can be represented by the non-negative spherical barycentric/area coordinates g =[ g 1 , g 2 , g 3 ] T of the vertices U l of the lth triangle, where g i ≥ 0 and i g i ≥ 1. Note that the required allpositive spherical barycentric coordinates are only found if a suitable spherical triangle l is selected from the convex hull, which will then contain u. While the spherical barycentric coordinates g reproduce the direction u, spherical triangular interpolation uses the corresponding planar barycentric coordinatesg i = g i j g j [27] to linearly interpolate the values measured at the microphones of the triangle l by their weighted average, At the boundaries, this interpolation exactly reproduces the values at the triangle vertices and linearly interpolates the value pairs along any edge of the lth triangle. Because neighboring triangles share edges and vertices, interpolation across triangles is continuous. There is no condition for the first-order derivatives, therefore this interpolation is C 0 smooth.

Robustness and bias
Robustness is often measured by observing the range of amplifications that stochastic perturbations linearly superimposed with the input data can undergo. Due to linearity, it is insightful and common practice to observe changes that uncorrelated Gaussian noise as the only input f = N undergoes, which we adopt to analyze the robustness of the three above-mentioned interpolation methods. We consider the 32 nodes of a pentakis dodecahedron as directional sampling for the input data f, which is interpolated using the 2520 nodes of a Chebyshev-type quadrature [43], yielding the 2520 output valuesf. Figure 2 shows a statistical analysis of the two ratios RMS{f} RMS{f} and max{f} max{f} , analyzing these ratios for 1000 independent instances of a random input vector f.
Regarding changes between RMS values from input to output, we observe that SH methods for N = {7, 8}, Spline for all M = {1, 2, 3}, and TrI produce a bias towards smaller output RMS values of 2 dB or more for stochastic input. For TrI, it is understandable that within any triangle, three uncorrelated inputs get averaged linearly, therefore the output RMS gets reduced by stochastic instead of additive interference. For SH interpolation with N = 4, this reduction only happens sparsely. The implicit minimization of the Euclidian norm for N ≥ 5 minimizes the output RMS value, and therefore causes the observable bias towards lower RMS. This minimization might be optimistically regarded as an increase in robustness between the sampling nodes for N ≤ 6, but it also implies a general decrease in magnitude there, a bias causing dips between the observed samples when interpolating omnidirectional directivities. Spline methods process constant inputs separately, therefore, it is reasonable to assume that the observed reduction in output RMS rather displays increased robustness to stochastic perturbation. For the chosen spatial sampling scheme, all methods appear robust enough to avoid enlarged output RMS values.
As a more critical test, SH-based interpolation exhibits the largest differences between maxima in the interpolated output compared to those in the input, with around ±3 dB for N = {4, 5, 6, 7}, while the settings SH 8 and Spline 3 behave reasonably. Rigorously, TrI as a linear interpolation is capable of precisely avoiding enlarged output maxima, and the same benefit is observed for Spline 1.

High-resolution reference directivities
For an objective evaluation of the estimation accuracy of a spatial interpolation method based on finite samples on a measurement surface, a high-resolution reference is required. This could theoretically be an analytical function sampled at the evaluation points. However, since the quality of the interpolation also depends significantly on the properties of the pattern to be interpolated, a diverse sample of musical instruments was chosen for which high-resolution measurements were made. The sample included a trombone, a violin, a flue pipe and a bassoon, and thus different types of sound production, different physical principles of sound radiation, and different sizes and geometries of the radiator. To achieve high reliability, the excitation of the instruments was automated, the instruments were rotated by a computer-controlled 3D loudspeaker measurement system ELF, and the sound pressure measurements were obtained on a dense spatial sampling grid. The measurements were conducted in the anechoic chamber of the OWL University of Applied Sciences and Arts in Lemgo. A 1/2" free-field equalized BK4190 cartridge was used as measurement microphone, placed 2 m from the sound source. The trombone, a member of the brass family, is a relatively small and straightforward sound source, with the bell being the only port from which sound energy is emitted. The directional dependence of sound radiation of brass instruments is rotationally largely symmetric with respect to the center axis of the bell, which is also the main radiation direction. With increasing frequency the main lobe of the directivity constricts, resulting in a more directional sound radiation. To determine the directivity of the trombone, the shortened instrument (without slide) was artificially excited with a sine sweep signal of the order 16 (2 16 samples ≈ 1.4 s @48 kHz sampling rate), emitted by a horn driver directly attached to the small end of the bell [44]. An equal angle sampling grid with an angular resolu-tion of 5 • in azimuth and colatitude was chosen and thus 2522 unique positions were measured.
The sound radiation for a violin, a string instrument, is partially determined by the parallel plates of the instruments' body which vibrate locally in different amplitudes and phases. Particularly at low frequencies the sound is additionally radiated through the characteristic open fholes that build a Helmholtz resonator in connection with the air cavity of the body. With a source extension of about 40 cm it is a medium-sized instrument. However, the two vibrating plates with a different local phasing cause interferences and therefore a distinct directional characteristic in the far-field. The directivity of a violin was measured exemplarily for the open A string (f 0 = 440 Hz), applying a repeatable bowing machine for excitation [45] and utilizing an equal angle grid with an angular resolution of 6 • for azimuth and colatitude, yielding the sound pressure at 1742 positions on a sphere.
The sound of a flue pipe of an organ is radiated through the mouth as well as through the open end of the resonator. The two spatially separated partial sound sources thus produce frequency-dependent directional characteristics, which get more complex with increasing frequency following the characteristics of a dipole and corresponding to the sound radiation behavior of a flute [2]. To measure the directivity of the flue pipe, a horn driver has been attached directly to its toe hole which was artificially excited, again with a sine sweep of order 16. An equal angle sampling grid within an angular resolution of 5 • was chosen, again yielding 2522 positions. The pipe used has a length of 51.3 cm with a diameter of 4.8 cm and a fundamental frequency of f 0 = 280 Hz [46].
The bassoon, a woodwind instrument of the double reed family, has a bell and numerous extended tone holes distributed irregularly across a long, bent corpus. The openings act as secondary sound sources depending on the fingering. The superposition of their radiated sound fields can cause a relatively complex directivity in the far- field. The sound radiation of a bassoon, fingered for the note Eb3 (f 0 = 156 Hz), was measured with an equal angle sampling distribution within a horizontal and vertical angular resolution of 5 • , applying a repeatable artificial excitation [47]. Accordingly, the data was acquired at 2522 positions on a sphere.

Interpolation of microphone array measurements
In the first step, the high-resolution reference data was sub-sampled at 32 microphone positions used in the Berlin-Aachen database of musical instruments [12]. The 32 sampling points are located at the vertices of a pentakis dodecahedron (Fig. 3) and were chosen as one possible sampling scheme to evaluate the interpolation techniques under realistic conditions. Except at the two poles, the 32 positions of the sparse grid are not contained in the reference grids. To account for this, the final sparse grid was generated from the closest positions of the highresolution reference grids. The resulting grid is called sample grid in the following and diverges from ideal sparse grid by 1.7 • /1.0 • on average, with a maximum deviation of 2.4 • /2.6 • . For the SH orders examined in this paper, the condition numbers of the ideal sparse grid are small due to the equal distribution of the sampling points (κ = 1.0 for N ≤ 4; κ < 2.0 for N ≥ 6). The large condition number of κ = 1.24 × 10 16 for N = 5, however, indicates that this SH order should not be used for the selected grid. The condition numbers increased only slightly due to using the nearest neighbors for the sample grid, i.e., κ < 2.0 for N ≤ 4 and N ≥ 6 still holds.
The interpolation functions were sampled at the corresponding high-resolution reference grid points, allowing a direct comparison between the interpolation result and the reference, i.e., the measured directivity.
The interpolation was done on the magnitude responses, i.e., the phase information was neglected for two reasons. First, interpolation of the phase spectrum is very susceptible to noisy data and errors can occur, particularly at high frequencies [49]. Second, natural sound sources, in contrast to artificial sound sources,  do not have a stationary phase response neither for a certain frequency nor for a certain radiation direction [20], which could be used for room acoustics simulation without further ado. Figure 1 shows a spectrogram for the note A4 (440 Hz) played by a trumpet (recording taken from [12]). While the amplitude of the fundamental and the overtones are almost constant during the observed time window, the phase of the trumpet signal fluctuates strongly and cannot be determined unambiguously. A proposal on how absolute-valued interpolated directivity patterns can be used for wave-based simulation methods is presented in Section 5.

Global error measure
For the evaluation of the interpolation algorithms, a global single-number error measure is proposed. To describe the mathematical accuracy of the interpolation, the difference of the sound pressure levels of interpolated directivity and reference directivity averaged over all directions could be used, in the way Arend et al. have done to describe the accuracy in interpolation of head-related transfer functions (HRTFs) [21]. However, to describe the acoustic effect of erroneous excitation of the sound field caused by an incorrect directivity, it seems important to consider the sound power radiated (incorrectly or correctly) and not level differences, which correspond to larger differences in power at high levels than at low levels. As a physically meaningful measure we therefore propose to calculate the sound power erroneously radiated with respect to the direction due to the interpolation error and to relate this to the total radiated sound power. To obtain such as measure, the relative error in radiated sound power can be calculated as the summed area weighted relative differences of the squared sound pressures of interpolation p 2 (θ r , φ r ) and reference p 2 (θ r , φ r ) over the R directions (θ r , φ r ) of the reference grid for r ∈ {1, 2, ..., R}, related to the summed area weighted squared sound pressure of the reference: where w a are the normalized area weights of the reference grid with R r=1 w a (θ r , φ r ) = 1. (24) Note that an error of = 0.5 states that the interpolated directivity emits 50% of the sound power in incorrect directions, whereas a value of = 0 shows that the radiation pattern of the interpolated fully matches the reference.

Results
Since we investigated tonal instruments, we restricted the analysis to fundamentals and overtones. In addition, we discarded tones whose spatially averaged energy was more than − 40 dB below the tone with the maximum average energy separately for each instrument and played note. Thus, only the radiation pattern of the fundamental at f 0 = 280 Hz and the first five overtones f 1 , ..., f 5 of the flue pipe were examined. For the violin the fundamental f 0 = 440 Hz and the nine first overtones f 1 , ..., f 9 were evaluated, and for the bassoon the fundamental at f 0 = 156 Hz and the first five overtones f 1 , ..., f 5 . The trombone was a special case exhibiting an unnatural overtone series due to the shortened instrument and the artificial excitation with a horn driver. Therefore only the first five resonance frequencies at 966 Hz, 2280 Hz, 4882 Hz, 7212 Hz, and 9179 Hz are considered in the following.
The results of the spatial interpolations of the bassoon's directivity for the second overtone at f 2 = 468 Hz are shown in Fig. 3, along with the related global error measures . Distinct differences emerge between the interpolation methods. The spherical triangular interpolation approach (TrI) shows the lowest reproduction error for this directivity with = 0.36, with the two distinctive indentations in the θ = 90 • ; φ = 90 • and 270 • ; 90 • radiation directions estimated quite accurately despite the small number of sample points.
For the spline interpolation, the global error ≈ 0.4 is almost constant across the order M. Differences between the orders can be seen mainly in the 90 • ; 90 • radiation direction, where the notch is well reproduced by the spline interpolation of order M = 1, however at the expense of larger errors in the transition areas between high and low radiation, which are better interpolated with order M = 2.
For SH interpolation of orders N = {4, 5}, the largest errors occur in the notch regions at 90 • ; 90 • and 270 • ; 90 • . As the SH order increases, these errors disappear, but indentations between the sample points become visible, most pronounced at N = 8. The lowest global error is obtained with order N = 7 and = 0.37, whereas the maximum error of = 0.49 occurs when interpolating with the theoretical optimal SH order of N = 4, considering the 32 points of the sample grid.
For other radiation patterns such as the bassoon's third overtone at f 3 = 624 Hz (Fig. 4), the interpolation methods used behave somewhat differently. At this slightly higher frequency, the spline interpolation of order M = 1 is superior with = 0.49. SH interpolations of orders N = {6, 7} again perform better than the reproduction with N = {4, 5} as well as N = 8, where indentations between the sample points become again visible. While for the second overtone triangular interpolation was superior to all other methods investigated, it now performs slightly worse than spline interpolation. An overview of global errors for all examined orders and musical instruments is shown in Fig. 5, where the individual error distributions contain between 5 and 10 partials, depending on the instrument. To assess the benefit of using musical instrument directivities from small sample grids, was additionally calculated for the trivial assumption of an omnidirectional directivity using the mean radiated energy over all directions.
Taking the median of the distribution over all 27 partial tones of the examined instruments analyzed as a measure of quality of the interpolation methods, the spline approach with order M = 1 shows the best result, closely followed by the triangular interpolation and the spline interpolations with orders M = {2, 3}.
The largest reproduction errors occur at the lowest and the highest examined SH order, i.e., at N = {1, 8}, whereas  the best SH interpolation is achieved with order N = 7. A detailed list of results is shown in Table 1.

Discussion
The directivities of four different tonal musical instruments were sub-sampled at 32 almost equally distributed points and interpolated using spherical triangulation, spherical splines and spherical harmonics. A global error measure was proposed to assess the quality of the various interpolation techniques.

Comparison of interpolation algorithms
It is obvious that the absolute quality of all interpolation methods strongly depends on the acoustic size of the sound source and the resulting complexity of the radiation pattern. Acoustically small sound sources like the trombone bell can be relatively precisely interpolated already with 32 measuring points. With a median value for the difference between interpolation and reference of < 0.3 for the spline and the triangulation approach and for the SH interpolation with N = {5, 6, 7}, more than 70% of the sound power is radiated in exactly the right direction on average using the interpolated directivities. For extended sources with more complex radiation patterns, however, the far-field directivities can be increasingly poorly estimated with a sparse sampling grid.
Since the spatial frequencies of acoustic sound sources can generally be assumed not to be limited, SH decomposition based on a finite number of microphones on the surface of a sphere allows a correct reproduction of this function only up to a cutoff frequency of kr < N, where k denotes the wave number and r the radius of the microphone array. Above this frequency errors occur due to spatial aliasing [50] (depending on the SH order of the function and the sampling grid) and series truncation [39]. The almost equally distributed Q = 32 points of the applied sampling grid support a maximum SH order of N = 4 to solve the linear system of equation through best fit in the least-squares sense (Eq. 6). Taking into account the measuring distance of 2 m between the microphone and the sound source, the reference directivity can be correctly reconstructed up to a cutoff frequency of only f c ≈ 108 Hz without aliasing and truncation errors. All partial tones investigated in this study, and all musical frequency components in general, are above this frequency. Closest to this frequency is the fundamental of the bassoon at f 0 = 156 Hz, which can be reconstructed most accurately with SH order N = 4, and a global error of = 0.03. In this case, the interpolation error increases with SH order because the system of equations becomes increasingly under-determined. As a consequence, the minimization of the wave-spectral power associated with the Moore-Penrose inversion of the SH transformation (Eq. 7) entails an increasingly poor interpolation between the sampling points.
Radiation patterns for frequencies well above the cutoff frequency of the microphone array, such as for the second harmonic of the bassoon at f 2 = 468 Hz, not only show larger errors in general, but can provide smaller errors for orders N ≥ 5. In these cases, the minimization of the wave-spectral power between the 32 sample points (c.f. Fig. 3), can lead to smaller errors due to a better tradeoff between truncation and aliasing errors. Figure 6 shows the global error for all six investigated directivities of the bassoon across frequency.
At this point it seems worthwhile to take a closer look at the chosen sampling grid. An analysis showed that an exact reproduction of the 32 magnitude values at the sampling points can only be achieved with SH orders of N ≥ 6 for all radiation patterns investigated in this study and contained in [12]. Finding the best SH order for interpolating sparsely sampled directivities can thus also be interpreted as optimizing the trade-off between the desired exact reproduction of the magnitude values at the sampling points and the undesired minimization of the wave-spectral energy which causes notches between the sampling points and increases with increasing SH order. For the selected sampling grid, the optimum appears on average at an SH order of N = 7. The triangular and spline interpolation, however, both show not only smaller median errors but also smaller 25th and 75th percentiles than SH interpolation. Since the spline technique was applied without regularization (λ = 0), the 32 sample points are always correctly reconstructed regardless of the order M. The same applies to triangular interpolation, where all sample points are reproduced correctly as well. Even at low frequencies, where SH interpolation with N = 4 estimates the reference almost physically correct, the triangular and spline interpolation perform comparably well (Fig. 6). Interestingly, the errors for the spline interpolation increase with increasing order. This may be explained by the fact that splines are piece-wise defined functions f that are continuous at the sampling points up to the (M −1)th derivative, i.e. f is C (M−1) smooth. The smoothness of the splines thus increases with increasing order M, whereas the smoothness of SH increases with decreasing order. In both cases some degree of non-smoothness increases the accuracy of the interpolation for the selected sample grid.

Generalization to different sample grids
In the first instance, the results of this study only hold for the selected sample grid with almost equally distributed 32 points. Even though this is likely to be a typical design and close to designs used for other measurements (Section 1), we provide a way to check the error for other designs as well. For finding the best interpolation algorithm for a specific sample grid, we provide the Matlab tool SourceInterp.m [51]. It calculates the global error for all interpolation methods evaluated in this paper, based on the publicly available high-resolution Bassoon radiation patterns [52] for the note F3 (f 0 = 175 Hz) and will be extended for other instrument directivities once they are made publicly available.

Comparison to SH interpolation with iterative sign retrieval
As detailed earlier, only the absolute magnitude response was interpolated due to the stochastic nature of the phase information of measured natural sound sources. In case of SH interpolation, however, this approach increases the required SH order. To counteract this, iterative semidefinite relaxation methods to find a suitable real-valued sign to an absolute-valued radiation pattern prior to SH interpolation were proposed by Zagala and Zotter [19]. To assess the quality of the triangular and spline interpolation with respect to the proposed sign-retrieval algorithm, we replicated the benchmark from [19]. Therefore, 50 radiation patterns were created with randomly generated standard normally distributed real-valued SH coefficients f nm ∈ R up to SH order N = 3. The absolute unsigned values of these radiation patterns were evaluated at Q = 64 extremal sampling points (cf. [37]) according to Eq. 3 and used as input for triangular interpolation, the firstorder spline interpolation and third order SH interpolation. Finally, the area weighted Mean Square Error (MSE) between the analytical reference f and the interpolation resultf was calculated for 2522 sampling points of the reference grid used above. The distribution of the MSE across the 50 radiation patterns is shown in Fig. 7. On average, the spline interpolation and the triangular method reconstruct the random directivity patterns only slightly better than the SDR-based sign retrieval with common-sign regions algorithm (median of 0.05 vs. 0.07). However, the dispersion of the error is considerably lower for the spline approach, with a 75th percentile of 0.06 compared to 0.12 for SDR rgn+dbl. It should also be noted that the spline method reconstructs one radiation pattern on average 30 faster than the SDR rgn+dbl algorithm and 26 times faster than triangular interpolation (using Matlab on a PC with Intel Core i5-6400 CPU @ 2.70 GHz and 16 GB RAM, Fig. 7). Note that the triangular interpolation AKsphTriInterp() is already optimized by searching only sub-lists of triangles extending into the octant of any regarded interpolated coordinate. The potential speed-up factor has an upper limit of 8 and practically reached 6 with the required control flow and overlap of the sub-lists.

Combination with models for the phase response
Time domain acoustical simulations and low-order SH decompositions may benefit from complex-valued directivities with both magnitude and phase information. In this regard, the interpolation methods discussed up to this point might not be optimal, yet. The unsigned (zero phase) SH interpolation causes acausal impulse responses that are symmetrical around t = 0 and the sign-retrieval algorithm generates an non-continuous phase response that might also deteriorate the corresponding time signals. In addition, due to treating each frequency separately, it remains to be clarified based on which criteria to smooth the phase responses across frequencies. A direct integration of a numerically derived absolute-value directivity into time domain simulations such as the finite difference time domain (FDTD) method was presented by Ahrens & Bilbao using an analytical point source phase response and SH interpolation [49]. However, our results suggest that spline interpolation is more suitable for interpolating sparsely sampled magnitude responses regardless of the phase. Therefore, for interpolating sparsely resolved directivity measurements of musical instruments for the use in time domain simulations, we recommend to first interpolate the magnitude response to a high-resolution grid using first-order splines, followed by a subsequent SH expansion according to [49]. Note, however, that a physically correct extrapolation of the sound field within the near field of the sound source is not possible even with this approach, despite subsequently added phase information.

Generalization to different musical instruments
The results of the current study first only apply to the four instruments studied, and cannot simply be generalized to all musical instruments or acoustical sound sources in general. Considering the instruments of the classical orchestra, however, all fundamental radiation mechanisms such as the vibrating soundboard (violin), the piston-like air vibration in the bell of brass instruments (trombone), the complex multi-source radiation of woodwind instruments (bassoon), and the air jet excitation of flute instruments (flue pipe) are represented, at least those with stationary excitation, i.e., not including percussion instruments.
It is noticeable that piece-wise interpolation methods such as the spherical spline and spherical triangular interpolation seem to provide a better reconstruction of complex radiation patterns in cases where the sampling theorem (kr < N) is not met. This behavior is largely consistent for the 27 patterns analyzed, all of them at frequencies violating the spatial sampling theorem and representing different source types and different degrees of complexity.
We thus expect that a similar result would also be obtained also for directional characteristics measured in the presence of a player. Players as reflecting or diffracting object next to the instrument have an influence on the directivity of a musical instrument [11]. For largely omnidirectional sources they make the radiation pattern more complex, in specific cases, one could imagine it to be smoothened, for highly directional sources such as the trombone at high frequencies, the influence will be negligible. We thus see no reason why interpolation methods for radiation patterns with players should behave fundamentally differently from those without players observed in the current study.

Conclusion
The performance of different interpolation techniques applied to sparsely sampled directivity measurements of acoustic sound sources depends on the sampling grid used but also on the radiation pattern of the sources themselves. Therefore, we evaluated three established approaches for interpolation from a low-resolution sampling grid using high-resolution measurements of a representative sample of musical instruments as a reference. The smallest global error on average occurs for thin plate pseudo-spline interpolation, with order 1 performing slightly better than orders 2 and 3. For interpolation based on spherical harmonics (SH) decomposition, the SH order and the spatial sampling scheme applied have a strong influence on the quality of the interpolation that is difficult to predict in individual cases. The piece-wise linear, spherical triangular interpolation finally provides almost as good results as the first-order spline approach, albeit with on average 20 times higher computational effort. Therefore, for spatial interpolation of sparsely sampled directivity measurements of musical instruments, the thin plate pseudo-spline method applied to absolutevalued data with order M = 1 is recommended and, if necessary, a subsequent modeling of the phase.