# A Real-Time Sound Field Rendering Processor
## 2 HO-FDTD Algorithm
### Boundary Condition  
A reflective boundary can be modeled as a locally reacting surface by assuming that a wave does not propagate along with the boundary surface, and the acoustical behavior is affected by the sound pressure and particle velocity perpendicular to the boundary surface. If a sound wave travels in a positive axis \((x,y,z)\) direction, the boundary impedance \(Z\) is represented by the sound pressure and the particle vibration through Equation (7) [23; 24]. \[Z=\frac{P}{U} \tag{7}\]  
Here, \(U\) is the particle velocity component perpendicular to the boundary. Differentiating both sides of Equation (7) and substituting \(U\) by the momentum conservation equation of wave propagation, the boundary conditions are obtained in terms of sound pressure [25]. \[\frac{\partial P}{\partial t}=-c\xi\
abla P \tag{8}\]  
where \(\xi=Z/\rho c\) is the normalized boundary impedance. For a rectangular sound space, boundary grids are classified into interior grids of a boundary, edges, and corners according to their position. Different formulas are applied to update sound pressures of different types of boundary grids because their conditions are different. For example, for the interior grids of right boundary, Equation (9) is derived by applying the centered finite difference method on Equation (8) and assuming the normalized boundary impedances of all boundaries are \(\xi\). \[\frac{P^{n+1}(i,j,k)-P^{n-1}(i,j,k)}{2\Delta t}=-c\xi\frac{P^{n}(i+1,j,k)-P^{n }(i-1,j,k)}{2\Delta x} \tag{9}\]  
By rearranging the terms in Equation (9) and introducing the parameter \(\chi\), Equation (10) is derived to express a virtual point \(P^{n}(i+1,j,k)\), which lies outside of the sound space [25]. \[P^{n}(i+1,j,k)=P^{n}(i-1,j,k)+\frac{1}{\chi^{\xi}}(P^{n-1}(i,j,k)-P^{n+1}(i,j,k) \tag{10}\]Substituting the related items in Equation (5) by Equation (10), then  
\[ \begin{array}{l}P^{n+1}(i,j,k)=\{\chi^{2}(2P_{i-1,j,k}^{n}(i-1,j,k)+P^{n}(i,j+1, k)+P^{n}(i,j-1,k)+P^{n}(i,j,k+1)\\ +P^{n}(i,j,k-1))+2(1-3\chi^{2})P^{n}(i,j,k)+(\frac{\chi}{\xi}-1)P_{i,j,k}^{n-1} (i,j,k)]/(\frac{\chi}{\xi}+1)\end{array}  \tag{11}\]  
If the reflection factor \(R\) is defined as \((\xi-1)/(\xi+1)\) and \(\chi\) is assumed to be \(1/2\), Equation (11) is changed to  
\[ \begin{array}{l}P^{n+1}(i,j,k)=\frac{1+R}{2(3+R)}[2P^{n}(i-1,j,k)+P^{n}(i,j+1,k)+P^{n}(i,j-1,k)\\ +P^{n}(i,j,k+1)+P^{n}(i,j,k-1)+2P^{n}(i,j,k)]-\frac{3R+1}{3+R}P^{n-1}(i,j,k) \end{array}  \tag{12}\]  
Equation (12) consists of two parts, one is the sum associated with the sound pressures of a grid and its neighbor grids at the time step \(n\), and another corresponds to the sound pressure of a grid at the time step \(n-1\). Compared with Equation (6), except the multiplicands, Equation (12) only replaces the sound pressure of the virtual grid \(P^{n}(i+1,j,k)\) by the sound pressure of the neighbor grid \(P^{n}(i-1,j,k)\) in the sum part. Moreover, for the interior grids on other boundaries, the multiplicands of two parts are the same while just substituting the sound pressure of the virtual grid by the sound pressure of the related neighbor grid in the summation. For example, when a grid is on the interior of the left boundary, the updated equation is achieved through substituting the sound pressure of the virtual grid \(P^{n}(i-1,j,k)\) with the sound pressure of the neighbor grid \(P^{n}(i+1,j,k)\) in Equation (12). The sum part is, therefore, changed from \((2P^{n}(i-1,j,k)+P^{n}(i,j+1,k)+P^{n}(i,j-1,k)+P^{n}(i,j,k+1)+P^{n}(i,j,k-1)+2 P^{n}(i,j,k))\) to \((2P^{n}(i+1,j,k)+P^{n}(i,j+1,k)+P^{n}(i,j-1,k)+P^{n}(i,j,k+1)+P^{n}(i,j,k-1)+2 P^{n}(i,j,k))\). The similar derivation procedure can be applied to edges and corners by using different boundary conditions. For example, when grids are on edges, which are intersections of two boundary planes, two boundary conditions are satisfied simultaneously. Expressions for two virtual points are consequently required to be derived.