# Binaural Auralization of Room Acoustics with a Highly Scalable Wave-Based Acoustics Simulation
## 2 Theory
### Expansion Coefficient Calculation via Spatial Differentiation of the FE Shape Function  
A novel method for computing expansion coefficients \(A_{l,m}\) of spherical harmonics via differentiation of the shape function is presented. The method uses the relation between \(a_{l,m}\) and the spatial differentiation of sound pressure noted in the literature [31]. According to the literature, we define a differentiation operator \(D_{l,m}\) to replace \(\gamma_{x}\), \(\gamma_{y}\) and \(\gamma_{z}\) in \(Y_{l,m}\) respectively with partial differential operators \(\partial/\partial x\), \(\partial/\partial y\) and \(\partial/\partial z\). Table 1 presents examples of correspondence between \(Y_{l,m}\) and \(D_{l,m}\) up to the first order. By applying the differential operator \(D_{l,m}\) to \(p(\mathbf{r},\omega)\) in Equation (23) yields the following relation:  
\[D_{l,m}p(\mathbf{r},\omega)=(\mathrm{jk})^{l}\sum_{l^{\prime}=0}^{\infty}\sum_{m^ {\prime}=-l^{\prime}}^{l^{\prime}}a_{l^{\prime},m^{\prime}}(\omega)\int_{0}^ {2\pi}\int_{0}^{\pi}Y_{l,m}Y_{l^{\prime},m^{\prime}}\mathrm{e}^{\mathrm{j}(kx \gamma_{x}+ky\gamma_{y}+kz\gamma_{z})}\sin\theta\mathrm{d}\theta\mathrm{d}\phi. \tag{29}\]  
From the equation presented above, the following relation between \(a_{l,m}\) and spatial differentiation of sound pressure can be derived as  
\[a_{l,m}(\omega)=(\frac{1}{\mathrm{jk}})^{l}D_{l,m}p(\omega), \tag{30}\]  
with a simple assumption of \(\mathbf{r}=0\), and with orthonormality of the spherical harmonics function, which is given as  
\[\int_{0}^{2\pi}\int_{0}^{\pi}Y_{l,m}Y_{l^{\prime},m^{\prime}}\mathrm{d}\theta \mathrm{d}\phi= \begin{cases}1&(l=l^{\prime}\text{ and }m=m^{\prime}),\\ 0&(l\
eq l^{\prime}\text{ or }m\
eq m^{\prime}).\end{cases}  \tag{31}\]  
In the literature [31], Equation (30) is transformed into a time-domain expression via inverse Fourier transform. Then, the time-domain expression is solved with multiplication of a low-order finite difference. However, finite difference approximation needs a regular computational grid. Consequently, it is unsuitable for FEM, which uses nonuniform FEs. In contrast to the approach presented earlier, this study directly evaluates \(a_{l,m}\) in the frequency domain with spatial differentiation of finite element shape functions. The present approach is available for arbitrarily-shaped FEs. With the present approach, the computed nodal RIRs on an eight-node hexahedral FE are first Fourier transformed to obtain a frequency-domain sound pressure vector \(\mathbf{p}(\omega)\) \begin{table}
\begin{tabular}{c c c} \hline \hline \((l,m)\) & \(\mathbf{Y_{l,m}}\) & \(\mathbf{D_{l,m}}\) \\ \hline \((0,0)\) & \(\sqrt{\frac{1}{4\pi}}\) & \(\sqrt{\frac{1}{4\pi}}\) \\ \((1,-1)\) & \(\sqrt{\frac{3}{4\pi}}\gamma_{y}\) & \(\sqrt{\frac{3}{4\pi}}\frac{\partial}{\partial y}\) \\ \((1,0)\) & \(\sqrt{\frac{3}{4\pi}}\gamma_{z}\) & \(\sqrt{\frac{3}{4\pi}}\frac{\partial}{\partial z}\) \\ \((1,1)\) & \(\sqrt{\frac{3}{4\pi}}\gamma_{x}\) & \(\sqrt{\frac{3}{4\pi}}\frac{\partial}{\partial x}\) \\ \hline \hline \end{tabular}
\end{table} 
Table 1: Correspondence between \(Y_{l,m}\) and \(D_{l,m}\) up to the first order. With \(\mathbf{p}(\omega)\), the spatial derivative of sound pressure at an arbitrary position \(\mathbf{r}_{a}\) in the FE defined by \(p(\mathbf{r}_{a},\,\omega)\) is calculated as  
\[ \begin{pmatrix}\frac{\partial p(\mathbf{r}_{a},\omega)}{\partial x}\\ \frac{\partial p(\mathbf{r}_{a},\omega)}{\partial y}\\ \frac{\partial p(\mathbf{r}_{a},\omega)}{\partial z}\end{pmatrix} = \begin{pmatrix} \frac{\partial N}{\partial x}\\ \frac{\partial N}{\partial y}\\ \frac{\partial N}{\partial z}\end{pmatrix} \mathbf{p}(\omega), \tag{32}\]  
where \(\mathbf{N}\) is the shape function vector. Because this study uses isoparametric elements represented using the local coordinate system (\(\xi\), \(\eta\), \(\zeta\)), the following chain rule is applied to Equation (32):  
\[