Bilbao, Stefan

Webb, Craig

English

Three dimensional ﬁnite difference time domain schemes can be used as an approach to spatial audio simulation. By embedding a model of the human head in a 3D computational space, such simulations can emulate binaural sound localisation. This approach normally relies on using high sample rates to give ﬁnely detailed models, and is computationally intensive.This paper examines the use of head models within audio rate FDTD schemes, ranging from 176.4 down to 44.1 kHz. Using GPU computing with Nvidia’s CUDA architecture, simulations can be accelerated many times over a serial computation in C. This allows efﬁcient, dynamic simulations to be produced where sounds can be moved around during the runtime. Sound examples have been generated by placing a personalised head model inside an anechoic cube. At the lowest sample rate, 44.1 kHz, localisation is clear in the horizontal plane but much less so in the other dimensions. At 176.4, there is far greater three dimensional depth, with perceptible front to back, and some vertical movement

Edinburgh Research Explorer

     Edinburgh Research Explorer                                      Binaural Simulations Using Audio Rate FDTD Schemes andCUDACitation for published version:Webb, C & Bilbao, S 2012, 'Binaural Simulations Using Audio Rate FDTD Schemes and CUDA'. inProceedings of the 15th International Conference on Digital Audio Effects. York, UK.Link:Link to publication record in Edinburgh Research ExplorerDocument Version:Publisher final version (usually the publisher pdf)Published In:Proceedings of the 15th International Conference on Digital Audio EffectsPublisher Rights Statement:©Webb, C., & Bilbao, S. (2012). BINAURAL SIMULATIONS USING AUDIO RATE FDTD SCHEMES ANDCUDA. In Proceedings of the 15th International Conference on Digital Audio Effects. York, UK.General rightsCopyright for the publications made accessible via the Edinburgh Research Explorer is retained by the author(s)and / or other copyright owners and it is a condition of accessing these publications that users recognise andabide by the legal requirements associated with these rights.Take down policyThe University of Edinburgh has made every reasonable effort to ensure that Edinburgh Research Explorercontent complies with UK legislation. If you believe that the public display of this file breaches copyright pleasecontact openaccess@ed.ac.uk providing details, and we will remove access to the work immediately andinvestigate your claim.Download date: 20. Feb. 2015Proc. of the 15th Int. Conference on Digital Audio Effects (DAFx-12), York, UK , September 17-21, 2012BINAURAL SIMULATIONS USING AUDIO RATE FDTD SCHEMES AND CUDACraig J. Webb∗Acoustics Group/Edinburgh Parallel Computing CentreUniversity of Edinburgh, UKC.J.Webb-2@sms.ed.ac.ukStefan BilbaoAcoustics GroupUniversity of Edinburgh, UKsbilbao@staffmail.ed.ac.ukABSTRACTThree dimensional finite difference time domain schemes canbe used as an approach to spatial audio simulation. By embeddinga model of the human head in a 3D computational space, such sim-ulations can emulate binaural sound localisation. This approachnormally relies on using high sample rates to give finely detailedmodels, and is computationally intensive.This paper examines the use of head models within audio rateFDTD schemes, ranging from 176.4 down to 44.1 kHz. UsingGPU computing with Nvidia’s CUDA architecture, simulationscan be accelerated many times over a serial computation in C.This allows efficient, dynamic simulations to be produced wheresounds can be moved around during the runtime. Sound exampleshave been generated by placing a personalised head model insidean anechoic cube. At the lowest sample rate, 44.1 kHz, localisa-tion is clear in the horizontal plane but much less so in the otherdimensions. At 176.4, there is far greater three dimensional depth,with perceptible front to back, and some vertical movement.1. INTRODUCTIONIn virtual acoustics applications, direct numerical methods are in-creasingly becoming a viable approach. Full wave-based tech-niques can give highly detailed and accurate results, and the finitedifference time domain method (FDTD) is particularly suitable foracceleration using GPU computing [1]. Such an approach allowsthe possibility of 3D localization, giving binaural results when us-ing stereo headphones.Research into the simulation of head-related transfer functions(HRTFs) has been conducted using techniques such as boundaryelement methods [2], and indeed finite difference methods [3].This research has focused on calculating accurate HRTFs, basedon generic or individual scans of the human head. In the case ofFDTD, very high sample rates are used (up to 500 kHz) in order todiscretize the models at a fine level of detail [4].This paper investigates 3D localization within audio rate FDTDsimulations, from 44.1 to 176.4 kHz. Rather than computing staticHRTFs based on impulse responses, sound can be localized withina fully dynamic virtual acoustic environment. These simulationsallow the movement of sound sources during runtime, which isdifficult to achieve using impulse response rendering, requiring in-terpolation.Section 2 outlines the discretization of head geometry over auniform grid. Section 3 describes the FDTD scheme employed,and Section 4 its implementation using CUDA. Results are pre-sented in Section 5.∗ This work was supported by the European Research Council, underGrant StG-2011-279068-NESS2. 3D HEAD MODELThe first stage in creating a binaural FDTD simulation is to obtainan accurate model of the human head. In order to minimise theinter-personal effect of HRTFs in testing, a model was constructedfrom scans of the primary author’s own head. This was producedusing a stereo-photographic process, giving 3D scans from variousangles which are then combined to give a complete model1. Theresult is a dense point cloud of coordinates indicating the surfaceof the head.With a view to parallelization in CUDA, the point cloud wasvoxelized into a 3D data grid. The grid density is set according tothe stability condition for the FDTD scheme (see Section 3) whichdictates a minimum spacing X for a given audio sample rate Fs.(In this article, a simple scheme is employed, so the relationshipto be satisfied is X ≥ √3c/Fs, where c is the wave speed in air.)The data was voxelized using the following algorithm:1. Iterate through the point cloud list, scaling the coordinatesto the required grid spacing and rounding them to the near-est point in the 3D data grid.2. Identify all points on or outside the surface boundary, bytracing rays from each face of the minimum parallelepipedsurrounding the head. What is left is considered to be theinterior of the head.3. Calculate the correct FDTD scheme coefficients for all pointsnot inside the surface.The voxelization for a grid spacing of 3.4mm (176.4 kHz)shows detailed resolution (Figure 1). Facial features and pinnadefinition are readily apparent.Figure 1: Voxelization at 3.4mm for 176.4 kHz, with a magnifiedview of the left pinna.1Stereo-photographic scanning provided by Prof. R. Fisher and S. Mc-Donagh, School of Informatics, University of Edinburgh.DAFX-1Proc. of the 15th Int. Conference on Digital Audio Effects (DAFx-12), York, UK , September 17-21, 2012A sample rate of 44.1 kHz requires a grid spacing of at least13.5mm, Figure 2 shows the result. Whilst the basic shape of thehead is still clear, the detail of the pinna is almost completelylost. At 88.2 kHz, the grid spacing of 6.8mm gives some defi-Figure 2: Voxelization at 13.5mm for 44.1kHz.nition of the shape of the pinna, certainly better than at 13.5mm.The effect of the different voxelizations was tested by embeddingthe head models into a 3D space, using a single variable FDTDscheme. This incorporates basic reflective boundaries on the headsurface with loss. The rectangular cuboid containing the head datais placed inside a larger space with numerically absorbing bound-ary conditions, as shown in Figure 3.Figure 3: Simulation setup in anechoic FDTD cube.Sound sources were input into the scheme during runtime, andthen dynamically moved around inside the simulation space.3. FINITE DIFFERENCE SCHEMEAcoustic propagation is modelled using the second order 3D waveequation:∂2u∂t2= c2∇2u (1)Here, u(x, y, z, t) is the target acoustical field quantity, c is thewave speed in air, and ∇2 is the 3D Laplacian. The simplest pos-sible FDTD scheme is of the “seven point" variety []:un+1l,m,p = (2− 6λ2)unl,m,p + λ2Snl,m,p − un−1l,m,p (2)where λ = cTX(T = Fs is the time step and X the grid spacing),and Snl,m,p is the sum of the six direct neighbouring points in un.The stability condition is derived using von Neumann analysis [5],leading to the constraint mentioned above:X ≥ cT√3 (3)Two types of boundary condition were used. First, the surface ofthe head is assumed to be rigid and a lossy condition is used:∂u∂t= cβn · ∇u (4)where n is the normal to the surface at the head boundary, and βis a simple loss coefficient. The integer 6 from the interior updateis reduced for each of the neighbouring points of un which fallinside the boundary, and Snl,m,p adjusted accordingly. It is thisinteger value which is stored in the head definition data.To create an anechoic boundary around the cube, an absorb-ing condition is required. Whilst a perfectly matched layer is thestandard approach from electro-magnetics [6], a computationallycheaper method can be employed based on the Enquist Majda con-dition [7]. If n is a normal to a planar boundary, then the conditioncan be written as: (∂∂t+ cn · ∇)qu = 0 (5)where q determines the order of the approximation. Although thisgives some minor reflections for high frequency waves which aretangential to the boundary, informal listening tests suggest that theconditions are reasonably transparent at audio sample rates, evenwith an order as low as q = 2.For a dynamic source location, appropriately sampled audiosignals were summed into the grid at each time step using tri-linearinterpolation. The input position was calculated along a 3D pathdefined by numerous points inside the space. Tri-linear interpola-tion provides a simple, but effective, method for moving the input,with minimal perceptible artefacts.4. IMPLEMENTATION USING CUDAInitial prototyping was performed in MATLAB. However, even atthe lowest sample rate of 44.1kHz and for a cube of less than 1m3,a serial C code takes around a minute for one second of simulation.At 176.4 kHz, the computation is 256 times larger and takes hours.Using GPU computing and Nvidia’s CUDA architecture, one cangreatly reduce the runtime [8]. Porting the code is simplified usingthe latest CUDA version 4.1 and the FERMI GPU cards.CUDA version 4 enables a three dimensional grid for the threadblocks [9], allowing a direct mapping of the 3D data to block in-dices. Using a FERMI architecture card, complex shared memorycoding is unnecessary as the caching system provides the samefunction [10]. Each thread can read data directly from global mem-ory, and setting the device to prefer the L1 cache gives a furtherimprovement in performance.Achieving optimal performance still relies on global memorycoalescing. This is a particular issue for the exterior boundaryupdates of the cube, where the faces and edges necessarily havedifferent orientations and their data is not always contiguous inmemory. Considerable speedup was obtained by rearranging theleft and right faces into temporary arrays to facilitate coalescing(Figure 4). Each Z-layer of these non-contiguous faces was loadedDAFX-2Proc. of the 15th Int. Conference on Digital Audio Effects (DAFx-12), York, UK , September 17-21, 2012Figure 4: Non-contiguous memory on left and right faces.into linear memory, which can then be accessed in a contiguousmanner when applying the update to the face boundaries. Thisrearranging of data can also be masked by performing it duringthe data loading for the interior points.For the interior and head data updates, a two stage approachwas used. An initial kernel performs an update over the entireinterior space, followed by a second kernel which corrects the gridfor the head definition. A further optimisation was possible for the44.1 and 88.2 kHz simulations, as the head definition data couldbe stored as single byte chars in constant memory.5. RESULTSTesting was performed on a single Nvidia Tesla C2050 card at dou-ble precision, for an outer cube of width 86 cm. Table 1 shows theresulting computation times for one second of simulation. MovingSample Rate (kHz) Total Grid points Time (sec)44.1 262,144 6.088.2 2,097,152 70.3176.4 16,777,216 1,360.6Table 1: Computation times for 1 second of simulation.from 88.2 to 176.4 kHz takes nearly twenty times longer, as thehead definition data is read from global memory.Initial audio tests used dry recordings of a matchbox and alsoa short vocal sample. During a six second simulation, this audiosignal was moved around the head on a flat circular path and alsoan angled path as shown in Figure 5. The left and right readoutpositions were taken from individual grid points on the boundaryof the head, approximately at the entrance to the ear canal. Notethat the ear canal itself was not modelled due to the coarse natureFigure 5: Source movement used for initial testing.of the grid used. Audio files are available at http://www2.ph.ed.ac.uk/~s0956654/Site/VirtualRoomAcoustics.htmlStarting at 44.1 kHz, movement around the azimuth angle ofthe horizontal plane is immediately apparent and the position canbe judged. What is less clear is the location either in front or be-hind the head, and variation in elevation angle is only slightly dis-cernible at the front. At 88.2 kHz, the result is similar to the previ-ous ones, but azimuth variation is more clearly defined. At 176.4kHz there is a distinct change and spatialisation effects are moreeasily discernible. Azimuth variation is very clear, and there issense of the movement from front to back, along with some effectsof changes in elevation. At all sample rates, the head definitiongives spatial effects which are clearly absent from the simulationsusing just an empty space.Figure 6: Semi-circular path at 40cm from centre of head.Analysing the system using impulse response simulations canbe performed to derive the head-related impulse response (HRIR)and HRTF for various source positions. Initial tests were per-formed using impulses at five degree increments around a semi-circlular path (Figure 6). Results are shown in Figures 7 and 8.6. CONCLUSIONS3D audio localization is of course based on various features: objec-tive factors such as the type of sound source used, and also neces-sarily subjective factors such as the individual nature of the HRTF.The purpose of this initial study was to examine if any form oflocalization was perceivable using audio rate FDTD schemes andbasic boundary conditions.At the lower sample rates of 44.1 and 88.2 kHz, only the vari-ation in the horizontal plane is clearly defined. Movement in el-evation and judging back to front is difficult. The audio resultsat 176.4 kHz are much improved, giving clear localisation in theDAFX-3Proc. of the 15th Int. Conference on Digital Audio Effects (DAFx-12), York, UK , September 17-21, 2012Figure 7: HRIR for left ear at 176 kHz.Figure 8: HRTF for left ear at 176 kHz.horizontal plane as well as perceptible back to front movement.Changes in elevation are apparent but harder to accurately place,and would benefit from the inclusion of the upper torso in themodel.Running simulations at sample rates greater than 44.1 kHz isfar less efficient, increasing the workload by x16 for each dou-bling of the rate. For large scale virtual acoustics modelling, this isclearly not advantageous, both in terms of computation times andindeed memory usage. Further development using more refinedboundary conditions for defining and voxelizing the head modelshould allow the effects of the 176.4 kHz simulation to be heardat the lower sample rates. An open question is the level of detailnecessary in the representation of the pinna. Future work will beconcerned with perceptual testing, with regard to the issues men-tioned above.7. REFERENCES[1] L. Savioja, D. Manocha, and M. Lin, “Use of GPUs in roomacoustic modeling and auralization,” in Proc. Int. Symposiumon Room Acoustics, Melbourne, Australia, Aug 2010.[2] Y. Kahana and P. Nelson, “Boundary element simulationsof the transfer functions of human heads and baffled pinnaeusing accurate geometric hearing,” J. Sound and Vibration,vol. 300, pp. 552–579, 2007.[3] P. Mokhtari, H. Takemoto, R. Nishimura, and H. Kato,“Computer simulation of HRTFs for personalization of 3Daudio,” in Proc. International Symposium on Universal Com-munications, Osaka, Japan, Dec 2008.[4] M. Nakazawa and A. Nishikata, “Develompent of a soundlocalization system with tube earphone using human headmodel with ear canal,” IEICE Trans. Fundamentals, vol.E88-A(12), pp. 3584–3592, Dec 2005.[5] J. Strikwerda, Finite Difference Schemes and Partial Dif-ferential Equations, Wadsworth and Brooks/Cole AdvancedBooks and Software, Pacific Grove, California, 1989.[6] J.P. Berenger, “Three-dimensional perfectly matched layerfor the absorption of electromagnetic waves,” J. Comp. Phys,vol. 127(2), pp. 363–379, Sept 1996.[7] B. Enquist and A. Majda, “Absorbing boundary conditionsfor the numerical evaluation of waves,” Mathematics of Com-putation, vol. 31(139), pp. 629–651, 1997.[8] C.J. Webb and S. Bilbao, “Computing room acoustics withCUDA - 3D FDTD schemes with boundary losses and vis-cosity,” in Proc. of the IEEE Int. Conf. on Acoustics, Speechand Signal Processing, Prague, Czech Republic, May 2011.[9] Nvidia Corp, “CUDA C Programming Guide v4,” Available: http://developer.nvidia.com/cuda-toolkit-41, May, 2011.[10] Nvidia Corp, “Fermi compute arch. whitepaper,” Available: http://www.nvidia.com/object/fermiarchitecture.html, Oct,2010.DAFX-4

Absorbing boundary conditions forthenumericalevaluationofwaves,”

Boundary element simulations of the transfer functions of human heads and bafﬂed pinnae using accurate geometric hearing,”

Computer simulation of HRTFs for personalization of 3D audio,” inProc.InternationalSymposiumonUniversalCommunications,

Computing room acoustics with CUDA - 3D FDTD schemes with boundary losses and viscosity,”

Develompent of a sound localization system with tube earphone using human head model with ear canal,”

Fermi compute arch. whitepaper,” Available : http://www.nvidia.com/object/fermiarchitecture.html,

Finite Difference Schemes and Partial Differential Equations,

Three-dimensional perfectly matched layer for the absorption of electromagnetic waves,”

Use of GPUs in room acoustic modeling and auralization,” in

Binaural Simulations Using Audio Rate FDTD Schemes and CUDA

Binaural Simulations Using Audio Rate FDTD Schemes and CUDA

Abstract

Similar works

Full text

Available Versions

Edinburgh Research Explorer