# A Review on Head-Related Transfer Function Generation for Spatial Audio
## 6 Interpolation
### Functional Models for HRTF Interpolation  
The first attempt at HRTF interpolation was functional representation, i.e., the calculation of a mathematical model that depends on frequency and direction, derived from a set of HRTF measurements. The first functional approach was proposed in 1967 by [PERSON], who modelled the external ear as a three-channel, two-delay system [174]. The same model was then used in [175; 176] to investigate localisation cues. In 1986, [PERSON] proposed a model based on a 16-channel filterbank, in which the filter parameters are related to the shape of the external ear [177]. Based on [PERSON]'s idea [177], [PERSON] and [PERSON] developed a model that applies proper filters and delays to reproduce the effect of the head, shoulders, and pinna [179; 180]. The model proposed in [179] consists of different modules obtained by combining properly parameterised delays and filters; each component of the model affects at least one of the three spatial dimensions, i.e., azimuth \(\theta\), elevation \(\phi\), and distance \(r\). Based on this idea, in 2007, [PERSON] et al. [211] proposed a structure dependent on \(\theta\) and \(\phi\) by generalising the continuous variable digital delay (CVDD) method introduced in [212]. The system performs polynomial interpolation, optimising the coefficients inside a given region of the sphere [211]. In 1992, [PERSON] et al. modelled the external ear as a multisensor broadband beamformer [178]. This approach is limited by the high computational complexity and numerical instability. Aiming to reduce the computational complexity of the model, [PERSON] introduced principal components analysis (PCA) in 1987 [186]. [PERSON] and [PERSON] extended this method in 1992 demonstrating the efficiency of representing a low-order HRTF [185]. PCA was also applied by [PERSON], who considered different HRTF measurements captured by in-ear microphones [187]. HRTF functions can also be simplified also by modelling low-order IIR filters [189; 190]. In this case, starting from the measured HRTF, the optimal coefficients of the IIR filter of any desired order that represents the HRTF are calculated [189]. Moreover, interpolation techniques, such as bilinear or triangular interpolation, can be applied to the IIR coefficients [190]. The methods discussed above can reduce the order of the HRTF, decreasing the computational cost. However, the weights are only optimised for measured directions, and the final model is not described by a continuous function. In 1995, [PERSON] et al. proposed a model based on Karhunen-Loeve expansion that defines HRTFs in the complex frequency domain as a series of orthogonal eigentransfer functions (EFs) [192]. This solution offers a continuous function that can predict the phase and magnitude information of unmeasured directions. Another more widespread and effective method that guarantees a continuous and orthogonal representation is the spherical harmonics decomposition (SHD). Following the Helmoltz reciprocity principle, SHD models the HRTF spectrum as the weighted sum of spherical harmonics [198]. In this way, the spectrum of a generic HRTF function in the spatial domain can be represented as follows:  
\[H(\theta,\phi,f)=\sum_{n=0}^{N-1}\sum_{m=-n}^{n}c_{n}^{m}(f)Y_{n}^{m}(\theta, \phi), \tag{26}\]  
where \(N\) is the truncation order of the SHD, \(c_{n}^{m}(f)\) are the SH coefficients, and \(Y_{n}^{m}(\theta,\phi)\) are the spherical harmonic basis function, defined as follows [213]:  
\[Y_{n}^{m}(\theta,\phi)=(-1)^{m}\sqrt{\frac{(2n+1)(n-|m|)!}{4\pi(n+|m|)!}}P_{n} ^{|m|}(\cos\theta)e^{im\phi}, \tag{27}\]  
where \(n=0,1,..,N\), \(m=-n,..,n\), and \(P_{n}^{m}(\cdot)\) is the associated Legendre function. Starting with the measured HRTFs, the SH coefficients can be analytically calculated by inverting Equation (26), or by applying an optimisation technique, as the least square (LS) method [200]. Once \(c_{n}^{m}(f)\) are computed, the HRTF of any position can be calculated using Equation (26).