Presented at the 8th International Conference on Auditory Display (ICAD), Kyoto, Japan, July 2-5, 2002.In order to develop the 3D auditory display using a head phone which allow head movement, real time calculation of the head related transfer function (HRTF) is necessary. In the conventional studies, the calculations are performed approximately based on the mathematical model by regarding the head as the sphere shape. Although the boundary element method is also possible to calculate the transfer function by solving the wave equation from the accurate boundary condition of the head including the shape of the face and the ears measured by the 3D scanner, it is thought to be impossible to calculate it in real time because the amount of the calculation is too big. In this study, we discuss a new calculation method of the HRTF based on the reciprocity principle which enables significant speed-up of the calculation.The demonstrations are available at the website: http://acoust.archi.kyoto-u.ac.jp/HRTF

Ise, Shiro

Otani, M.

English

Scholarly Materials And Research @ Georgia Tech

Proceedings of the 2002 International Conference on Auditory Display, Kyoto, Japan, July 2-5, 2002REAL TIME CALCULATION OF THE HEAD RELATED TRANSFER FUNCTION BASEDON THE BOUNDARY ELEMENT METHODShiro Ise Makoto OtaniDepartment of Architecture and Architectural SystemsFaculty of Engineering, Kyoto University.Yoshida-Honmachi, Sakyo-ku, Kyoto 606-8501, Japanise@archi.kyoto-u.ac.jpABSTRACTIn order to develop the 3D auditory display using a head phonewhich allow head movement, real time calculation of the head re-lated transfer function (HRTF) is necessary. In the conventionalstudies, the calculations are performed approximately based on themathematical model by regarding the head as the sphere shape. Al-though the boundary element method is also possible to calculatethe transfer function by solving the wave equation from the ac-curate boundary condition of the head including the shape of theface and the ears measured by the 3D scanner, it is thought to beimpossible to calculate it in real time because the amount of thecalculation is too big. In this study, we discuss a new calculationmethod of the HRTF based on the reciprocity principle which en-ables significant speed-up of the calculation.The demonstrations are available at the website:http://acoust.archi.kyoto-u.ac.jp/HRTF/1. INTRODUCTIONGenerally, human obtains spatial informations through binaurallistening, which is one of the functions of the phonoreception. Par-ticularly, it is known that human perceives the direction of soundsources through the interaural sound level differences and time dif-ferences of the ear input signals [1].The head-related transfer functions (HRTFs), which are thetransfer functions between the sound source and listener’s ears,contains such informations that the interaural sound level differ-ences and time differences of the ear input signals. The spectralcues caused by these informations is said to be related to the soundlocalization. If a source signal and the HRTFs are convolved andthe resulting signal can be precisely produced at both ears, it wouldbe possible to let the listener perceive the sound source in that di-rection. In order to measure the HRTFs, we usually use a realhead or a dummy head. However, because we can not measurethe HRTFs in every directions, we have to predict the HRTFs inunmeasured direction. It is possible to interpolate HRTFs in theunmeasured direction from those in the measured directions [2].However, if the HRTFs are needed for the other distance from thesource, a number of measuring HRTFs are needed. Drastic ap-proximation, such that the head shape is approximated to sphereor spheroids [3], is also realistic if we may ignore the differencebetween the front and the back.On the other hand, it is possible to solve the wave equationstrictly by using the finite/boundary element method [4, 5]. Al-though these methods can calculate HRTFs accurately, the calcu-lation time is extremely long compared with the conventional cal-culation methods.In this study, we describes the calculation method of the HRTFsby using the conventional BEM, the Chief method [6] which min-imize the error caused at the internal eigen-frequencies and thereciprocity method which enables speeding up of the calculation.In the numerical study, these methods are compared and are dis-cussed. Furthermore, we introduce the server system which calcu-late the HRTFs based on the BEM.2. CALCULATING HRTFS BY USING THE BEMThe BEM which calculates the HRTFs is categorized as the ex-ternal problem of the boundary integral equation. As the bound-ary condition, the boundary surface is given by the shape of thedummy head. By solving the integral equation, we obtain theacoustic pressure at ears when a sound source is located at a posi-tion considered.2.1. Calculation of the acoustic pressure on the boundary sur-faceThe BEM is a numerical method which solves the algebra equationby discretizing the boundary surface of the Kirchhoff Helmholtzintegral equation expressed as follows:C(s)p(s) = jωρNk=1qkG(r′k |s)+ SG(r|s)∂p(r)∂n−p(r)∂G(r|s)∂ndSC(s) = 1/2 (s ∈ S)1 (s ∈ V )0 (s /∈ V ), (1)whereICAD02-1Proceedings of the 2002 International Conference on Auditory Display, Kyoto, Japan, July 2-5, 2002S : SurfaceV : Volumes : Position within the volume Vr′k : Point source position (k = 1 . . .N)r : Position on the surface Sp : acoustic pressureG : the Green functionqk : the volume velocity of the source kω : angular frequencyρ : the air densityn : the normal vector.In Eq.(1), the normal derivative is given as follows,∂p(r)∂n= −jωρvn(r) = −jωρyn(r),where vn and yn are the particle velocity and the acoustic admit-tance in normal, respectively.The boundary surface S is divided into M elements for dis-cretization. In this study, we suppose the constant element whichregard the physical parameters such as the acoustic pressure andthe particle velocity are constant within the discretized element.Then ,if the j-th element is set to Sj and the center of Sj is de-fined as rj , Eq.(1) can be descretized asC (s)p(s)= jωρNk=1qkG(r′k|s)−Mj=1jωρvn(rj) SjG(r|s)dS−Mj=1 Sj∂G(r|s)∂ndS+jωρyn(rj) SjG(r|s)dSp(rj),(2)wherevn(rj /∈S′)=0,yn(rj∈S′)=0.Here, if the vector r1, r2, r3, · · · , rM which are located onthe surface S is substituted for s in Eq.(2), we can solve them assimultaneous equations of many unknown pr i(i = 1 · · ·M). Asri is on the surface S, that is C(s) = 1/2, Eq.(2) can be writtenas12p(ri) = jωρNk=1qkG(r′k|ri)−Mj=1jωρvn(rj)Gij−Mj=1Gnij +jωρyn(rj)Gij	p(rj),whereGij = SjG(r|ri)dS,Gnij = Sj∂G(r|ri)∂ndS.By transposing p(ri) in the third term in the right side of Eq.(3),we can express simply as the matrix equation.12IM + Gn + jωρGYP = gA − jωρGV (3)whereGn=Gnij ∈ CM×MG = [Gij] ∈ CM×MY = diag(yn(r1), yn(r2), yn(r3), · · ·, yn(rM )) ∈ CM×MP = [p(r1), p(r2), p(r3), · · · , p(rM )]T ∈ CMg =G(r′k |ri) ∈ CM×NA = [A1, A2, A3, · · · , AN ]T ∈ CNV = [v(r1), v(r2), v(r3), · · · , v(rM )]T ∈ CMTherefore, if the boundary conditions such as the source am-plitude, the vibration velocity and the acoustic admittance on theboundary surface are given and if (IM/2 + Gn + jωρGY) is theregular matrix, the acoustic pressure P̂ on the surface is given byP̂ =12IM + Gn + jωρGY−1(gA − jωρGV) . (4)2.2. Calculation of the acoustic pressure within the volumeIn the same manner of deriving Eq.(eq:discrete1), the acousticpressure within the volume V (p(s ∈ V )) can be written asp(s) = gsA − jωρGsV − (Gns + jωρGsY) P̂ (5)wheregs =G(r′k|s) ∈ CN ,Gs = SjG(rj |s)ds∈ CM ,Gns = Sj∂G(rj |s)∂nds∈ CM .Therefore, after obtaining the acoustic pressure on the surface S,that is P̂, by solving Eq.(4), we can obtain the acoustic pressure atany position within the volume V .2.3. Reciprocity principleSuppose a sound fields driven by a sound source A with the volumevelocity v located at ra within the closed surface S, the sound fieldis expressed as the non-homogeneous wave equation as follows.(∇2 + k2)pa(r) = −jωρ0v · δ(r − ra) (6)In the same manner, if we suppose another sound field driven bya sound source B with the volume velocity v located at r b withinthe closed surface S, the sound field is also expressed as the non-homogeneous wave equation as follows.(∇2 + k2)pb(r) = −jωρ0v · δ(r − rb) (7)By substituting rb for r in Eq.(6) and ra for r in Eq.(7) and com-paring both, we can easily derive that pb(ra) = pa(rb). The reci-procity principle in the wave equation can be transformed into theintegral equation aspa(s) = jωρvG(ra|s)+ G′(r|s), (8)pb(s) = jωρvG(rb|s)+ G′(r|s), (9)ICAD02-2Proceedings of the 2002 International Conference on Auditory Display, Kyoto, Japan, July 2-5, 2002whereG′(r|s) = SG(r|s)∂p(r)∂n−p(r)∂G(r|s)∂ndS.If ra is the sound source position and r b is the ear position, theHRTF can be obtained by calculating both p a(rb) and pb(ra).Because pb(ra) is the acoustic pressure at the source position r awhen the sound source is located at the ear position r b, HRTF canbe solved by exchanging the sound source position and the ear po-sition. This fact makes great merit in the HRTFs calculation. Ifwe calculate the pressure signal P̂ with the sound source locate atthe ear position rb based on Eq. 4 once, the HRTF at any sourceposition can be obtained by calculating Eq. 5.In the calculation of the BEM, the most part of the calculationis occupied with the first calculation which calculates P̂ by solv-ing Eq. 4. The second calculation finishes with a instant. However,because the source point is close to the boundary surface, the ac-curacy of solving P̂ can reduce. If the problem of the internaleigen-frequency mentioned below coincide it, further degradationof the accuracy is expected.2.4. Cheif methodIn solving Eq. (4), if (1/2IM + Gn + jωρGY) is the regularmatrix, the pressure on the boundary surface, that is P̂, have anunique solution. It can, however, not be regular at all the frequen-cies. In particular, at the high frequencies, Eq.(4) tends to havemany eigen-frequencies and not to have any unique solution. Inthe internal problem of the integral equation, it is physically nat-ural because of the resonance inside of the closed volume. It is,however, unnatural to cause the same phenomena outside of theclosed volume. We do not usually observe the resonance outsideof the closed volume. This problem can be solved by changing thematrix. The Chief method is one of the solution by adding obser-vation points where the acoustic pressure becomes zero, i.e. insidethe volume in the external problem. In this method, first s /∈ Vare chosen and are substituted for s in Eq.(1). Next, the new equa-tions are attached to Eq. (4). Because the equations are more thanthe number of the variable which is the acoustic pressure on thesurface Sj(j = 1 · · ·M), we must solve by using the least meansquare method. If we choose tildeM points of s /∈ V and defineas r1,r2, , · · · , rM , the equation is given as follows:0 = jωρNk=1qkG(r′k|ri)− Mj=1jωρvn(rj) Gij−Mj=1 Gnij +jωρyn(rj) Gij p(rj),whereGij =  SjG(r|ri)dS,Gnij =  Sj∂G(r|ri)∂ndS.By expressing the matrix equation, it becomesGn + jωρ GYP = gA − jωρ GV, (10)whereGn=  Gnij ∈ CM×M ,G =  Gij ∈ CM×M ,g = G(r′k|ri) ∈ CM×N .(11)By attaching Eq.(10) to Eq.(3), it is expressed asHP = L, (12)whereH = 12IM + Gn + jωρGYGn + jωρ GY∈ C(M+M)×M ,L =ggA − jωρGGV ∈ C(M+M). (13)In this case, as Eq.(12) is the over-determined system, P does nothave any solutions. However, we can obtain P which has the leasterror of the equation by solvingP̂ = (HT H)−1HT L. (14)3. COMPUTER SIMULATION3.1. Boundary condition of head shapeWe prepared triangle elements group which shapes the surface ofthe dummy head. First, the 3D coordinate data was measuredby using 3D scanner (Cyberware 3030RGB/I). Next, the trianglemesh data was calculated by using the triangle mesh generator [7].Although the measured coordinate data of the dummy head sur-face include the body of the head torso simulator, we use only thehead part of the data in this study.It is said that the element size should be less than 1/4 to 1/6of the wave length to obtain enough accuracy in the calculationfor the BEM. Therefore, we made upper limit of the element sizewhen the surface is divided into the triangle element. Besides,because the ear pinna could affect the accuracy of the calculation,we controlled the upper limitation size near the ear pinna to besmall enough. After determining these condition, the boundaryelements are made as shown in Fig.1.The number of the triangle element in this model is 5014. Theaverage size of the element is about 6.8mm. Although the inlet ofthe external ear canal is closed in this model, it can be regarded thatit is the same condition as the HRTF measured under the conditionof the external ear canal closed which is said to include the spatialinformation [8].The sound source is locate at an angle of front, right 30 degreeand a distance of 1.4m from the center of the head torso. We cal-culate the acoustic pressure at the position of left and right ear. Inorder to confirm the accuracy of the computation, we calculate theacoustic pressure inside of the head, too.3.2. HRTF using the conventional BEMAs the conventional method, we solve Eq.(4) first and Eq.(5) nextto obtain the HRTFs. The frequency range of the calculation is43-12015Hz at intervals of 43Hz. Fig.2 shows the HRTFs in theICAD02-3Proceedings of the 2002 International Conference on Auditory Display, Kyoto, Japan, July 2-5, 2002Figure 1: Computer modelfrequency domain and the solid line and the dashed line indicatethe S.P.L. at the right ear and left ear, respectively. The dotted lineis the error level.0 2 4 6 8 10 12−120−100−80−60−40−200frequency [kHz]relative SPL [dB]right earleft earinternal pointFigure 2: HRTF in the frequency domain by the conventional BEMAs shown in Fig.3, it can be seen that the amplitude of theright ear is higher than that of the left ear and the wave arrived atthe right ear is earlier than that at the left ear. time is 45 hoursusing Xeon 2.5GHz.3.3. HRTF calculation using the reciprocity principleFig.4 shows the HRTFs calculated by the method based on thereciprocity principle, that is s in Eq.(5) is set as the source positionand r′k in Eq.(4) is set as the ear position. Fig.5 is the HRTFs inthe time domain, respectively.The computation time for the first calculation is 45 hours andthat for the second calculation is 2 seconds. The difference be-0 2 4 6 8 10−0.10.1time [ms] amplituderight earleft earFigure 3: HRTF in the time domain by the conventional BEM0 2 4 6 8 10 12−120−100−80−60−40−200frequency [kHz]relative SPL [dB]right earleft earFigure 4: HRTF in the frequency domain by the reciprocitymethodtween the peak and the dip in Fig.4 are less than the conventionalmethod. The noise reduction can be seen before the peak of theHRTF waveform in Fig.5.3.4. HRTF calculation using Chief methodFig.6 shows the HRTFs calculated by the Chief method, that isp(s) in Eq.(5) is calculated by substituting the acoustic pressure P̂on the surface S given by Eq.(14).The error level in the frequency domain in Fig.6 does notincreases compared with the conventional method illustrated inFig.2. The computation time is 116 hours using Xeon 2.5GHz.ICAD02-4Proceedings of the 2002 International Conference on Auditory Display, Kyoto, Japan, July 2-5, 20020 2 4 6 8 10−0.10.1time [ms] amplituderight earleft earFigure 5: HRTF in the time domain by the reciprocity method0 2 4 6 8 10 12−120−100−80−60−40−200frequency [kHz]relative SPL [dB]right earleft earinternal pointFigure 6: HRTF in the frequency domain by the Chief method4. DEVELOPMENT OF REAL TIME HRTF SERVERWe developed the computation server as shown in Fig.8 which cal-culate the HRTFs quickly. The URL address of it is as follows:http://acoust.archi.kyoto-u.ac.jp/HRTF/The computation time is 2 seconds for one position of thesound source.At the moment, it can not be said the perfect real time server,but comparing the computation time in the conventional BEM, thatis 45 hours, it can be said that the revolutionary speed up of thecalculation achieved. We are going to open this server for anyapplication which will uses HRTF after confirming the copyrightconcerning the triangle mesh generator [7].0 2 4 6 8 10−0.10.1time [ms] amplituderight earleft earFigure 7: HRTF in the time domain by the Chief method5. REFERENCES[1] Jens Blauert, Spatial Hearing (The MIT Press, Cambridge,Massachusetts, London, England 1997)[2] “Interpolation of the head related transfer function on thehorizontal plane,” Takanori Nishino, Shoji Kajita, KazuyaTakeda ,and Fumitada Itakura, J. Acoust. Soc. Jpn.., 55, 91-99 (1999) (in Japanese).[3] “Calculation of acoustic coefficients between two ears usingspheroids,” Kiyoshi Sugiyama, Tatsuya Sakaguchi, ShigeakiAoki ,and Ikuichiro Kinoshita, J. Acoust. Soc. Jpn.., 51, 117-122 (1995) (in Japanese).[4] “Boundary Element Simulation of HRTFs and Sound FieldsProduced by Virtual Acoustic Imaging Systems,” Yuvi Ka-hana, Philip A. Nelson, M. Petyt, and Sunghoon Choi, Au-dio Engineering Society 105th Convention, pre-print 4817(1998)[5] M. Otani and S. Ise, “Numerical calculation of the headrelated transfer function by using the boundary elementmethod”. In Proc. of WESTPRAC VII, pages 305–308, 2000.[6] H.A.Schenck , “Improved integral formulation for acousticradiation problems”. Journal of the Acoustical Society ofAmerica, 44,41-58 ,(1968)[7] J.Schoeberl, “NETGEN - An advancing front 2D/3D-meshgenerator based on abstract rules”, Computing And Visuali-sation in Science, 1997[8] H. Mφller, “Fundamentals of Binaural Technology”,AppliedAcoustics, 36, 171-218, (1992)ICAD02-5Proceedings of the 2002 International Conference on Auditory Display, Kyoto, Japan, July 2-5, 2002Figure 8: HRTF serverICAD02-6

Real time calculation of the head related transfer function based on the boundary element method

https://smartech.gatech.edu/bitstream/1853/51372/1/IseOtani2002.pdf

Real time calculation of the head related transfer function based on the boundary element method

Abstract

Similar works

Full text

Available Versions

Scholarly Materials And Research @ Georgia Tech