- Email: [email protected]
A boundary element phantom head for binaural hearing simulation
Engineering Analysis with Boundary Elements 30 (2006) 309–314 www.elsevier.com/locate/enganabound
Technical note
A boundary element phantom head for binaural hearing simulation Y. Kagawa *, S. Ohnogi, L. Chai Department of Electronics and Information Systems, Akita Prefectural University, 84-4 Tsuchiya-Ebinokuchi, Honjo, Akita 015-0055, Japan Received 17 May 2005; received in revised form 8 September 2005; accepted 26 October 2005 Available online 26 January 2006
Abstract A numerical phantom head with boundary element modeling is proposed and demonstrated for binaural hearing investigation. The solution is involved with the scattering and diffraction. The boundary integral expression corresponding to the Helmholtz equation is numerically solved by the help of the boundary element discretization. Two possible formulations are considered; one is the formulation for the total wave field, and another is the formulation in which the wave field consists of the incident wave and the scattered wave. The solutions of the two methods are found to be completely coinciding. The authors believe that as it provides great flexibility the modeling paves the way to the application not only to the field distribution evaluation but also to the listening experiments with the digital filter whose transfer characteristics is based on the present numeral modeling. q 2005 Elsevier Ltd. All rights reserved. Keywords: Boundary element modeling; Directional sensation; Physical phantom
1. Introduction Important factors for the directional sensation in hearing are due to the differences in sound pressure and phase (or time delay) between a pair of ears, which is known as the binaural effect. The difference in pressure and phase depends on the difference between the distances from a source position to the each position of ears. As there is a head between the two ears, diffraction and scattering due to its presence cannot be ignored. This effect, which is sometimes defined as the head transfer function, is usually investigated by using a physical phantom head with two microphones provided at the ears positions [1]. The present paper proposes and demonstrates a numerical phantom using a boundary element modeling. The boundary element method is a finite element version of the boundary integral equation method, and becomes increasingly popular also for acoustic field analysis [2], since its first introduction under the name of boundary element method by Birebia [3]. Kagawa et al. demonstrated the use of the boundary element method for the unbounded field problems [4], and for the vocal tract modeling with radiation field [5] and the vocal tract shape identification problems [6].
* Corresponding author. Tel./fax: C184 27 2092. E-mail address: [email protected] (Y. Kagawa).
0955-7997/$ - see front matter q 2005 Elsevier Ltd. All rights reserved. doi:10.1016/j.enganabound.2005.10.005
The advantages of the boundary element modeling include the dimensional reduction from three-dimensional space into the two dimensional surface in which the object surface is only taken into account with the inclusion of the infinite radiation space. On the other hand, there are some disadvantages. The discretized algebraic equation becomes full. Though its dimension of the coefficient matrix is much smaller than that of the finite element counterpart, which is banded, the computation is not always cheaper. The solution is sometimes involved with fictitious resonance, as the interior and exterior of the surface in consideration is only distinguished by the sign of the directional normal taken to the boundary surface in the formulation. The fictitious resonance should carefully be removed. Another disadvantage is that since wave number is scatteredly present in the coefficient matrix, the standard solver for the eigenvalue problems, which is not involved in our present problem, cannot be used. 2. Boundary elements for acoustics 2.1. Total wave formulation We consider a field exterior to the head–body model shown in Fig. 1. The Helmholtz equation is the governing equation for the steady-state wave field, which is V2 p C k2 p Z dðrKrs Þ ðin domain UÞ and the boundary condition is
(1)
310
Y. Kagawa et al. / Engineering Analysis with Boundary Elements 30 (2006) 309–314
Fig. 2. A half sphere scatterer. Fig. 1. Sound scattering.
vp Z q^ vn
ðin the surface GÞ
The incident wave is defined as (2)
where p is sound pressure, k(Zu/c,u,: angler frequency and c: sound speed) is wave number. r is the position vector and rs is the source position vector. (^) indicates the value prescribed ^ 0 for the rigid surface. and qZ The boundary integral expression corresponding is ð ð vj vp Ci pi C p dG Z j p dGKjðrKrs Þ (3) vn vn G
G
where piZp(ri) and Ci (Z1/2 on a smooth surface) is a constant that depends on the solid angle at point i on the surface. j is the fundamental solution of the Helmholtz equation, which is given by jkR
jZ
e 4pR
(4)
(7)
The Helmholtz equation for the scattered wave is 2 ðsÞ
V p C k2 pðsÞ Z 0
ðin domain UÞ
(8)
and the boundary condition is vpðsÞ vpðiÞ ^ h qðsÞ Z qK vn vn
ðon the surface GÞ
The boundary integral expression is ð ð ð vj ðsÞ vpðsÞ vpðiÞ p dG Z j p dG ZK j p dG C vn vn vn
ci pðsÞ i
G
G
(9)
(10)
G
and the boundary element expression to match is ½HfpðsÞ g Z ½GfqðsÞ g (s)
where R is the distance between the source point and the point of consideration (RZjrKri)). For the boundary elements, the surface is here divided into triangular elements. Proper interpolation functions are assumed for the sound pressure and its derivatives or flux in each element. With the compatibility conditions are imposed at the connecting nodes, one has the dicretized expression of the form ½Hfpg Z ½Gfqg C fbg
pðiÞ Z p0 ejkR
(11)
(s)
where {p } and {q } are the pressure and flux vectors for the scattered wave. 3. Numerical verification for a sphere 3.1. Sound pressure on a rigid sphere for plane wave incidence The program has already been developed for linear triangular elements, found in our previous work [5,6], which
(5)
[H],[G] are the matrices to match the corresponding pressure and flux vectors {p} and {q}, and is the source vector {b}. The equation is solved for the nodes under the boundary condition prescribed, from whose values on the boundary the pressure at an arbitrary point in space is calculated. 2.2. Scattered wave formulation Total pressure p consists of incident wave p(i) and scattered wave p(s) p Z pðiÞ C pðsÞ
(6)
Fig. 3. The sound pressure on the surface in different angle for the plane wave incidence. (The lines indicate the analytical solutions).
Y. Kagawa et al. / Engineering Analysis with Boundary Elements 30 (2006) 309–314
where po is the sound pressure at the center of the sphere and pa is the pressure at an arbitrary point on the surface, jn and nn are, respectively, spherical Bessal and spherical Neumann function, and Pn Legendre function. By virtue of the symmetry, a half sphere is only considered for the calculation for which the element division is illustrated in Fig. 2. The numbers of the total elements and the associated nodes are, respectively, 1104 and 579. The boundary element numerical solutions are compared with the analytical solutions for the sound pressure on the ^ 0Þ, when a sphere is illuminated by the rigid surface ðqZ plane wave. Fig. 3 shows the results for the case kaZ4(aZ0.1 m), in which the average element length to the wavelength is about one-twelfth. The two numerical solutions well agree with the analytical except for large ka.
Fig. 4. Sound pressure on the sphere surface for point source.
is here used with a little modification. The fictitious resonance associated with the interior space are suppressed by providing the constraint of several null points as in our previous cases [6,7]. The fictitious resonances are moved out of the frequency range of interest. Fig. 2 shows a rigid sphere used for the examinations, for which the analytical solution is available [8]. That is
j
311
N pðsÞ C pðiÞ p 1 X jnC1 ð2n C 1Þ Pn ðcos fÞ j Zj a j Z 2 2 0 po po k a nZ0 jn ðkaÞKjnn0 ðkaÞ
(12)
3.2. Sound pressure distribution over a rigid sphere for a point source Fig. 4 shows the frequency characteristics of the pressure at the various points on the space surface. The pressure distribution, which shows the diffraction effect due to the presence of a sphere, is illustrated for various ka in Fig. 5, when it is illuminated by a point source placed 1 m apart. The sphere is a simplest head model for binaural effect consideration. The pressure is referred to the value at the center of the sphere. For the evaluation, both algorithms, Eqs. (5) and (11) are used, which require comparable computation time and give the same
Fig. 5. Pressure distribution on the sphere surface.
312
Y. Kagawa et al. / Engineering Analysis with Boundary Elements 30 (2006) 309–314
results. The results given in Fig. 4 are somewhat different from those in Fig. 3 as the distance between the sphere and the source point is not very long (R/aZ10). It is found that as R/a becomes larger, the results approach to those of Fig. 3 as expected. 4. Head–body model The binaural effect is then examined for a more elaborate head–body model as shown in Fig. 1. Fig. 6 shows the boundary elements. This is the similar model to that used in our previous work [5], in which the vocal tract is removed. The numbers of the elements and the corresponding nodes are 1744 and 881. Fig. 7 shows a plane that includes a pair of ears and the source point. The frequency characteristics are evaluated at the positions corresponding to the ears positions referring to the pressure at the center of the head. The results are shown in Fig. 8, which shows that as the source position rotates the pressure difference between the pressure at the right ear, pR and that of the left ear pL becomes significant. Fluctuation in the frequency characteristic as approaching to 2 kHz could not be trusted, since the average element length to the wavelength approaches one- quarter there. Fig. 9 illustrates the field distribution of equi-pressure and equi-phase contours for different frequencies and different source positions. The distributions also demonstrate that the sound pressure and phase differences depend on the source position and they become significant as the frequency increases.
Fig. 7. A plane that includes a pair of ears and source.
The approach provides a numerical phantom head alternative to the physical phantom head. The transfer functions thus obtained can be realized in an on-line digital filter, which can be used not only for binaural hearing but also for auralization sensation experiments.
5. Final remarks With the boundary element modeling, the problem is solved for a head–body model as an exterior acoustic field, and the binaural effect is well demonstrated.
Fig. 6. A head–body model and surface elements.
Fig. 8. Pressure at each ear.
Y. Kagawa et al. / Engineering Analysis with Boundary Elements 30 (2006) 309–314
Fig. 9. Pressure and phase in exterior field.
313
314
Y. Kagawa et al. / Engineering Analysis with Boundary Elements 30 (2006) 309–314
The paper was partly presented at JSST Computational Electromagnetic and Electronics Symposium [9]. References [1] For examples; Brual & Karjaerr Catalog. http://www.bksv.com and CAS communication analysis system, HEAD Acoustics GmbHhttp://www. head-acoustics.de. [2] See; Internet websitehttp://www.boundary-element-method.com. [3] Brebbia CA. The boundary element mehtod for engineers. London: Pentech Press; 1978. [4] Kagawa Y. Finite/boundary element methods for unbounded field problems. Tokyo: Science-sha; 1983 [in Japanese].
[5] Kagawa Y, Shimoyama R, Yamabuchi T, Murai T, Takarada T. Boundary element models of the vocal tract and radiation field and their response characteristics. J Sound Vib 1992;157(3):385–403. [6] Kagawa Y, Ohtani Y, Shinoyama R. Vocal tract shape identification from formant frequency spectra—a simulation using three-dimensional boundary element models. J Sound Vib 1997;203(4):581–96. [7] Schenck HA. Improved integral formulation for acoustic radiation problems. J Acoust Soc Am 1967;44:41–58. [8] Miida Y. Acoustical engineering. Tokyo: Shokodo; 1987 [in Japanese]. [9] Ohnogi S, Kagawa Y, Murai T, Yamabuchi T. Boundary element simulation of sound scattering/diffraction by human head. The 11th computation electromagnetic and electronics, JSST; 1990. p. 293–8.
Technical note
A boundary element phantom head for binaural hearing simulation Y. Kagawa *, S. Ohnogi, L. Chai Department of Electronics and Information Systems, Akita Prefectural University, 84-4 Tsuchiya-Ebinokuchi, Honjo, Akita 015-0055, Japan Received 17 May 2005; received in revised form 8 September 2005; accepted 26 October 2005 Available online 26 January 2006
Abstract A numerical phantom head with boundary element modeling is proposed and demonstrated for binaural hearing investigation. The solution is involved with the scattering and diffraction. The boundary integral expression corresponding to the Helmholtz equation is numerically solved by the help of the boundary element discretization. Two possible formulations are considered; one is the formulation for the total wave field, and another is the formulation in which the wave field consists of the incident wave and the scattered wave. The solutions of the two methods are found to be completely coinciding. The authors believe that as it provides great flexibility the modeling paves the way to the application not only to the field distribution evaluation but also to the listening experiments with the digital filter whose transfer characteristics is based on the present numeral modeling. q 2005 Elsevier Ltd. All rights reserved. Keywords: Boundary element modeling; Directional sensation; Physical phantom
1. Introduction Important factors for the directional sensation in hearing are due to the differences in sound pressure and phase (or time delay) between a pair of ears, which is known as the binaural effect. The difference in pressure and phase depends on the difference between the distances from a source position to the each position of ears. As there is a head between the two ears, diffraction and scattering due to its presence cannot be ignored. This effect, which is sometimes defined as the head transfer function, is usually investigated by using a physical phantom head with two microphones provided at the ears positions [1]. The present paper proposes and demonstrates a numerical phantom using a boundary element modeling. The boundary element method is a finite element version of the boundary integral equation method, and becomes increasingly popular also for acoustic field analysis [2], since its first introduction under the name of boundary element method by Birebia [3]. Kagawa et al. demonstrated the use of the boundary element method for the unbounded field problems [4], and for the vocal tract modeling with radiation field [5] and the vocal tract shape identification problems [6].
* Corresponding author. Tel./fax: C184 27 2092. E-mail address: [email protected] (Y. Kagawa).
0955-7997/$ - see front matter q 2005 Elsevier Ltd. All rights reserved. doi:10.1016/j.enganabound.2005.10.005
The advantages of the boundary element modeling include the dimensional reduction from three-dimensional space into the two dimensional surface in which the object surface is only taken into account with the inclusion of the infinite radiation space. On the other hand, there are some disadvantages. The discretized algebraic equation becomes full. Though its dimension of the coefficient matrix is much smaller than that of the finite element counterpart, which is banded, the computation is not always cheaper. The solution is sometimes involved with fictitious resonance, as the interior and exterior of the surface in consideration is only distinguished by the sign of the directional normal taken to the boundary surface in the formulation. The fictitious resonance should carefully be removed. Another disadvantage is that since wave number is scatteredly present in the coefficient matrix, the standard solver for the eigenvalue problems, which is not involved in our present problem, cannot be used. 2. Boundary elements for acoustics 2.1. Total wave formulation We consider a field exterior to the head–body model shown in Fig. 1. The Helmholtz equation is the governing equation for the steady-state wave field, which is V2 p C k2 p Z dðrKrs Þ ðin domain UÞ and the boundary condition is
(1)
310
Y. Kagawa et al. / Engineering Analysis with Boundary Elements 30 (2006) 309–314
Fig. 2. A half sphere scatterer. Fig. 1. Sound scattering.
vp Z q^ vn
ðin the surface GÞ
The incident wave is defined as (2)
where p is sound pressure, k(Zu/c,u,: angler frequency and c: sound speed) is wave number. r is the position vector and rs is the source position vector. (^) indicates the value prescribed ^ 0 for the rigid surface. and qZ The boundary integral expression corresponding is ð ð vj vp Ci pi C p dG Z j p dGKjðrKrs Þ (3) vn vn G
G
where piZp(ri) and Ci (Z1/2 on a smooth surface) is a constant that depends on the solid angle at point i on the surface. j is the fundamental solution of the Helmholtz equation, which is given by jkR
jZ
e 4pR
(4)
(7)
The Helmholtz equation for the scattered wave is 2 ðsÞ
V p C k2 pðsÞ Z 0
ðin domain UÞ
(8)
and the boundary condition is vpðsÞ vpðiÞ ^ h qðsÞ Z qK vn vn
ðon the surface GÞ
The boundary integral expression is ð ð ð vj ðsÞ vpðsÞ vpðiÞ p dG Z j p dG ZK j p dG C vn vn vn
ci pðsÞ i
G
G
(9)
(10)
G
and the boundary element expression to match is ½HfpðsÞ g Z ½GfqðsÞ g (s)
where R is the distance between the source point and the point of consideration (RZjrKri)). For the boundary elements, the surface is here divided into triangular elements. Proper interpolation functions are assumed for the sound pressure and its derivatives or flux in each element. With the compatibility conditions are imposed at the connecting nodes, one has the dicretized expression of the form ½Hfpg Z ½Gfqg C fbg
pðiÞ Z p0 ejkR
(11)
(s)
where {p } and {q } are the pressure and flux vectors for the scattered wave. 3. Numerical verification for a sphere 3.1. Sound pressure on a rigid sphere for plane wave incidence The program has already been developed for linear triangular elements, found in our previous work [5,6], which
(5)
[H],[G] are the matrices to match the corresponding pressure and flux vectors {p} and {q}, and is the source vector {b}. The equation is solved for the nodes under the boundary condition prescribed, from whose values on the boundary the pressure at an arbitrary point in space is calculated. 2.2. Scattered wave formulation Total pressure p consists of incident wave p(i) and scattered wave p(s) p Z pðiÞ C pðsÞ
(6)
Fig. 3. The sound pressure on the surface in different angle for the plane wave incidence. (The lines indicate the analytical solutions).
Y. Kagawa et al. / Engineering Analysis with Boundary Elements 30 (2006) 309–314
where po is the sound pressure at the center of the sphere and pa is the pressure at an arbitrary point on the surface, jn and nn are, respectively, spherical Bessal and spherical Neumann function, and Pn Legendre function. By virtue of the symmetry, a half sphere is only considered for the calculation for which the element division is illustrated in Fig. 2. The numbers of the total elements and the associated nodes are, respectively, 1104 and 579. The boundary element numerical solutions are compared with the analytical solutions for the sound pressure on the ^ 0Þ, when a sphere is illuminated by the rigid surface ðqZ plane wave. Fig. 3 shows the results for the case kaZ4(aZ0.1 m), in which the average element length to the wavelength is about one-twelfth. The two numerical solutions well agree with the analytical except for large ka.
Fig. 4. Sound pressure on the sphere surface for point source.
is here used with a little modification. The fictitious resonance associated with the interior space are suppressed by providing the constraint of several null points as in our previous cases [6,7]. The fictitious resonances are moved out of the frequency range of interest. Fig. 2 shows a rigid sphere used for the examinations, for which the analytical solution is available [8]. That is
j
311
N pðsÞ C pðiÞ p 1 X jnC1 ð2n C 1Þ Pn ðcos fÞ j Zj a j Z 2 2 0 po po k a nZ0 jn ðkaÞKjnn0 ðkaÞ
(12)
3.2. Sound pressure distribution over a rigid sphere for a point source Fig. 4 shows the frequency characteristics of the pressure at the various points on the space surface. The pressure distribution, which shows the diffraction effect due to the presence of a sphere, is illustrated for various ka in Fig. 5, when it is illuminated by a point source placed 1 m apart. The sphere is a simplest head model for binaural effect consideration. The pressure is referred to the value at the center of the sphere. For the evaluation, both algorithms, Eqs. (5) and (11) are used, which require comparable computation time and give the same
Fig. 5. Pressure distribution on the sphere surface.
312
Y. Kagawa et al. / Engineering Analysis with Boundary Elements 30 (2006) 309–314
results. The results given in Fig. 4 are somewhat different from those in Fig. 3 as the distance between the sphere and the source point is not very long (R/aZ10). It is found that as R/a becomes larger, the results approach to those of Fig. 3 as expected. 4. Head–body model The binaural effect is then examined for a more elaborate head–body model as shown in Fig. 1. Fig. 6 shows the boundary elements. This is the similar model to that used in our previous work [5], in which the vocal tract is removed. The numbers of the elements and the corresponding nodes are 1744 and 881. Fig. 7 shows a plane that includes a pair of ears and the source point. The frequency characteristics are evaluated at the positions corresponding to the ears positions referring to the pressure at the center of the head. The results are shown in Fig. 8, which shows that as the source position rotates the pressure difference between the pressure at the right ear, pR and that of the left ear pL becomes significant. Fluctuation in the frequency characteristic as approaching to 2 kHz could not be trusted, since the average element length to the wavelength approaches one- quarter there. Fig. 9 illustrates the field distribution of equi-pressure and equi-phase contours for different frequencies and different source positions. The distributions also demonstrate that the sound pressure and phase differences depend on the source position and they become significant as the frequency increases.
Fig. 7. A plane that includes a pair of ears and source.
The approach provides a numerical phantom head alternative to the physical phantom head. The transfer functions thus obtained can be realized in an on-line digital filter, which can be used not only for binaural hearing but also for auralization sensation experiments.
5. Final remarks With the boundary element modeling, the problem is solved for a head–body model as an exterior acoustic field, and the binaural effect is well demonstrated.
Fig. 6. A head–body model and surface elements.
Fig. 8. Pressure at each ear.
Y. Kagawa et al. / Engineering Analysis with Boundary Elements 30 (2006) 309–314
Fig. 9. Pressure and phase in exterior field.
313
314
Y. Kagawa et al. / Engineering Analysis with Boundary Elements 30 (2006) 309–314
The paper was partly presented at JSST Computational Electromagnetic and Electronics Symposium [9]. References [1] For examples; Brual & Karjaerr Catalog. http://www.bksv.com and CAS communication analysis system, HEAD Acoustics GmbHhttp://www. head-acoustics.de. [2] See; Internet websitehttp://www.boundary-element-method.com. [3] Brebbia CA. The boundary element mehtod for engineers. London: Pentech Press; 1978. [4] Kagawa Y. Finite/boundary element methods for unbounded field problems. Tokyo: Science-sha; 1983 [in Japanese].
[5] Kagawa Y, Shimoyama R, Yamabuchi T, Murai T, Takarada T. Boundary element models of the vocal tract and radiation field and their response characteristics. J Sound Vib 1992;157(3):385–403. [6] Kagawa Y, Ohtani Y, Shinoyama R. Vocal tract shape identification from formant frequency spectra—a simulation using three-dimensional boundary element models. J Sound Vib 1997;203(4):581–96. [7] Schenck HA. Improved integral formulation for acoustic radiation problems. J Acoust Soc Am 1967;44:41–58. [8] Miida Y. Acoustical engineering. Tokyo: Shokodo; 1987 [in Japanese]. [9] Ohnogi S, Kagawa Y, Murai T, Yamabuchi T. Boundary element simulation of sound scattering/diffraction by human head. The 11th computation electromagnetic and electronics, JSST; 1990. p. 293–8.