# A Review on Head-Related Transfer Function Generation for Spatial Audio
## 6 Interpolation
### Nearest Neighbour Interpolation  
The nearest neighbour interpolation techniques are based on calculating the interpolated HRIR \(\hat{h}(n)\) as the combination of the nearest measured impulse responses as follows:  
\[\hat{h}(n)=\sum_{i=1}^{N}w_{i}h_{i}(n), \tag{19}\]  
where \(N\) is the number of the measured impulse responses \(h_{i}\) considered in the interpolation, and \(w_{i}\) are the weights that can be estimated in several ways. The most simple interpolation technique is linear interpolation [171], in which two impulse responses are involved (\(N=2\)), and the weights are computed as  
\[w_{1}=\frac{x_{2}-\hat{x}}{x_{2}-x_{1}},\quad\text{and}\quad w_{2}=\frac{\hat{ x}-x_{1}}{x_{2}-x_{1}}, \tag{20}\]  
where \(\hat{x}\) is the position of the interpolated point, and \(x_{1}\) and \(x_{2}\) are the positions of the two measured impulse responses. In the HRTF dataset, where the measurements depend on the azimuth and the elevation, \(x\) resembles the azimuth or the elevation angles, and HRIRs should have the same elevation or the same azimuth, respectively. In 1993, [PERSON] and [PERSON] employed linear interpolation for a perceptual study conducted on individual HRTFs [171]. In 2018, [PERSON] used linear interpolation as part of the entire algorithm applied to binaural room responses, considering a horizontal measurement distribution [33]. This approach was also used by [PERSON] et al. in 2020, who considered a spherical distribution [172]. Finally, in 2023, [PERSON] et al. applied linear interpolation in combination with a frequency-warping approach to interpolate HRTF for non-linear 3D audio systems [27]. Frequency-warping is applied to reduce the resolution at high frequencies, improving binaural reproduction performance. One widespread technique is bilinear interpolation, as presented by [PERSON] and [PERSON] in 1999 [183]. In this case, \(N=4\) impulse responses are involved, and the four weights are calculated as follows:  
\[w_{1}=(1-c_{\theta})(1-c_{\phi}),\ \ \ w_{2}=c_{\theta}(1-c_{\phi}),\ \ \ w_{3}=c_{\phi}c_{\phi},\ \ \text{and}\ \ \ w_{4}=(1-c_{\theta})c_{\phi}, \tag{21}\]  
where the coefficients \(c_{\theta}\) and \(c_{\phi}\) are computed as  
\[c_{\theta}=\frac{C_{\theta}}{\theta_{\mathrm{s}}}=\frac{\theta\mathrm{mod} \theta_{\mathrm{s}}}{\theta_{\mathrm{s}}},\ \ \ \text{and}\ \ \ c_{\phi}=\frac{C_{\phi}}{\phi_{\mathrm{s}}}=\frac{\phi\mathrm{mod} \phi_{\mathrm{s}}}{\phi_{\mathrm{s}}}, \tag{22}\]  
where \(\theta_{\mathrm{s}}\) and \(\phi_{\mathrm{s}}\) are the angular steps, and \(C_{\phi}\) and \(C_{\phi}\) are the relative positions for azimuth and elevation, respectively, as shown in Figure 6a. The angular steps \(\theta_{\mathrm{s}}\) and \(\phi_{\mathrm{s}}\) must be fixed, so the positions of the dataset should be homogeneously distributed on a grid on the sphere surface. In 2002, [PERSON] et al. improved the bilinear technique by adding a new auxiliary inter-positional transfer function (IPTF) and reducing the number of points needed to three, resulting in a triangular interpolation in which the weights are calculated as follows [25; 184]  
\[w_{1}=1-w_{2}-w_{3},\ \ \ \ w_{2}=\frac{\Delta\theta_{\mathrm{A}}-w_{3}\Delta \theta_{\mathrm{AC}}}{\Delta\theta_{\mathrm{s}}},\ \ \ \ \ \text{and}\ \ \ \ \ w_{3}=\frac{\Delta\phi}{\Delta\phi_{\mathrm{s}}}, \tag{23}\]  
where the angular distances are determined as \(\Delta\phi=\phi-\phi_{\mathrm{A}}\), \(\Delta\theta_{\mathrm{s}}=\phi_{\mathrm{C}}-\phi_{\mathrm{A}}\), \(\Delta\theta_{\mathrm{s}}=\theta_{\mathrm{B}}-\theta_{\mathrm{A}}\), \(\Delta\theta_{\mathrm{A}}=\theta-\thetaAnother technique that can be used for calculating weights is vector base amplitude panning (VBAP), which was introduced by [PERSON] in 1997 [PERSON] (1997) for loudspeaker reproduction and used by [PERSON] et al. in 2011 [PERSON] et al. (2011) for HRTF interpolation. Three-dimensional VBAP involves \(N=3\) measured impulse responses, and the weights are obtained as follows:  
\[ \begin{bmatrix}w_{1}\\ w_{2}\\ w_{3}\end{bmatrix} =\mathbf{L}^{-1}\mathbf{p}=