# A Review on Head-Related Transfer Function Generation for Spatial Audio
## 4 Physical Modelling
### Boundary Element Method
BEM is a numerical method for solving the integral Equation (14) by means of a discretisation of the boundary \(\Sigma\) with a mesh (see Section 3 on the acquisition of meshes). In order to have high-quality simulations, the mesh should fulfil some requirements (see [99]):  
* As a rough rule, for resolving a wavelength \(\lambda\), the diameter \(d\) of elements of the mesh should be contained 6 to 10 times in \(\lambda\). For instance, at a frequency of 6000 Hz and with a speed of sound of 343 m/s, the wavelength \(\lambda\) is approximately 5.7 cm, and the diameter \(d\) should not exceed a few millimetres. For a discussion of this requirement, see [100; 101]. In order to resolve the geometry of the pinna, the elements should not exceed the length of 1 to 2 mm, depending on the local curvature of the surface. * All the triangles are close to being equilateral. For planar quadrilateral meshes, the requirement is that the lengths of the edges of each element are as similar as possible. * At each vertex, the number of concurring triangles is as close as possible to six. * The rate of change of the area of the triangles across the mesh is gradual. Once a high-quality mesh is available, the functional operators of Equation (15) in the Kirchhoff-Helmholtz integral Equation (14) can be discretised by associating discrete values of the fields to discrete elements of the mesh. A mesh provides a partition of the boundary \(\Sigma\) into \(N\) non-overlapping elements \(\Sigma_{j}\), where \(j\) is a discrete index. The discretisation of the field \(u\) and its normal derivative \(q\) in Equation (14) are based on the choice of shape functions \(S_{j}\) associated with the elements \(\Sigma_{j}\), which are typically piece-wise polynomials. A field \(\phi\) on \(\Sigma\), either \(u\) or its normal derivative \(q\), is discretised as follows:  
\[\phi(\mathbf{x})=\sum_{j=1}^{N}S_{j}(\mathbf{x})\phi_{j}. \tag{16}\]  
When a set of basis functions \(S_{j}\) has been chosen, there are several ways to obtain a discretised system of linear equations from the Kirchhoff-Helmholtz Equation (14). Here, we mention the following two methods:  
* In the collocation schema, a set of collocation points \(x_{j}\) is chosen and Equation (14), with fields \(u\) and \(q\) expanded as in Equation (16), is required to hold at the collocation points. * In the Galerkin schema, the problem is written in a weak form using a set of test functions, usually the same as the shape functions. For the collocation schema, implementation is simpler and, for this reason, is frequently considered the first choice. When the observation point in Equation (14) is on the surface \(\Sigma\), the coefficient \(c(\mathbf{r})\) is equal to \(\frac{1}{2}\), and the result of the discretisation process is one of the following linear problems, where the first equation refers to the collocation schema and the second one to the Galerkin schema:  
\[\frac{1}{2}\mathbf{u}_{\text{sc}}=G_{\text{c}}q-H_{\text{c}}\mathbf{u}_{\text {sc}},\qquad\quad\frac{1}{2}B_{\text{g}}\mathbf{u}_{\text{sc}}=G_{\text{g}} \mathbf{q}-H_{\text{g}}\mathbf{u}_{\text{sc}}, \tag{17}\]  
where \(\mathbf{u}_{\text{sc}}\) is the \(N\)-dimensional vector of the pressure values associated with the \(N\) elements of the mesh, \(\mathbf{q}\) is the corresponding vector associated with the normal derivatives of the pressure field, \(G_{\text{c}}\) and \(H_{\text{c}}\) are square matrices of dimension \(N\) corresponding to the operators \(L\) and \(M\) in Equation (15) in the collocation schema, \(G_{\text{g}}\) and \(H_{\text{g}}\) are the corresponding matrices for the Galerkin schema, and \(B_{\text{g}}\) is a real-valued, sparse mass matrix. The computation of the components of these matrices requires numerically integrating shape functions on the surface \(\Sigma\). In both schemata, the integral operators in Equation (15) must be evaluated, by linearity, on the shape functions, which can be achieved via numerical methods.