ICSV14Cairns • Australia
9-12 July, 2007

THE HEAD RELATED TRANSFER FUNCTION SIMULATION BY
FEM/IEM AND RECIPROCAL THEOREM
Hiroshi Muraoka1, Yuta Takahashi1, Tubagus Haedar1, Makoto Otani2 and Tatsuya Hirahara2
1

Research Center of Computational Mechanics, Inc.
Togoshi NI-Bldg. 1-7-1 Togoshi, Shinagawa-ku, Tokyo, 142-0041, Japanmuraoka@rccm.co.jp
2
Department of Intelligent Systems Design Engineering
Toyama Prefectural University, Imizu City, Toyama, 939-0398, Japan

Abstract
An effective method of reducing computational cost involving Head Related Transfer Function(HRTF) simulation for high frequency is to apply the reciprocal theorem. Commonly, theBoundary Element Method is used since both calculation and mesh generation effort areminimal, however it is sensitive to geometry complexity and unstable during the discretizationof the reciprocal theorem. Hence, in order to eliminate the errors and increase the efficiency ofthe outcome, this report demonstrates a new approach to Head Related Transfer Function(HRTF) simulations by coupling Finite Element Method (FEM) with Infinite Element Method(IEM). First, the validity of the reciprocal theorem was examined by comparing numericalresults along with the impulse response up to 1000Hz. Once the accuracy has been established,further studies were carried out between 10 and 20KHz to illustrate the practicability of thismethod.

1. INTRODUCTION
HRTF describes how a particular sound wave is filtered by the diffraction and reflection of
facial features before reaching the eardrum. It is quite essentially how we identify a sound and it
is affected by our individual facial characteristics (head, nose, ear, cheekbone) along with theposition and strength of the sound source. Due to computational improvement over the years,
higher frequency calculations are no longer out of reach. Although being restricted tohomogeneous problems and without uniqueness in the solution, BEM [1] has been considered
as a reference method since it offers reduction in dimensionality and the fact that it handles thesound radiation at the far field rigorously.
When retorting FEM [2], two methods are considered; one is FEM/IEM and the other is
DtN-FEM. Domain-based methods such as FEM/IEM are more suitable for solving exterior
acoustic problems than the ones based on DtN [3] due to the sparse matrices problem andefficient solution procedures such as parallel calculations (frequency and domain wise).
Coupling FEM with IEM will enable a natural extension of finite element to an unboundedICSV14 • 9-12 July 2007 • Cairns • Australia

domain. FEM was utilized, as it is not affected by geometrical complexity since it is based on
the weak variational formulation and sound hard condition is automatically accounted as anatural boundary condition. Although the FEM offers high accuracy results, it is well known
that the computational load is heavy due to the need of discretization of the whole
computational domain (large stiffness matrix). For instance, regarding high frequency
simulations, the elements around the head, sound source and within the volumetric space need
to be refined in a manner to capture the wavelength of the monopole, which culminates in a finegrid generation, increased computational time and memory usage. Integrating the reciprocal
theorem to the FEM/IEM does not only solve the dimensional issue but also facilitates theprocess of obtaining the output, such as directivity, as unlimited number of observation points
can be placed into a single calculation.

Figure 1. Two acoustic problems diagram

Consider a case in which a particular sound of frequency is generated by two unit point
sources around a solid body S. The reciprocal theorem [4] would indicate that sound pressure
obtained at B from a source at A is equal to one obtained at A from a source at B as follows
G ( x A , xB , ω ) = G ( xB , x A , ω )(1)

The functional forms of the respective velocity potentials generated by the two sources are asfollows
Sound pressure from A satisfies (∆ + k 2 )G ( x, x A , ω ) = δ ( x − xA )(2)Sound pressure from B satisfies (∆ + k 2 )G ( x, xB , ω ) = δ ( x − xB )(3)

Wave number , k =

2π fc

(4)

This theorem aids the process of computational reduction by minimizing the analysis domain,which is highlighted in this paper. This paper is organized as follows: After this introduction,
the governing equations for the weak variational formulation are presented in section 2. Section
3 outlines the problem definition and section 4 shows some numerical results. The penultimate
section discusses high frequency results and the final section summarizes the work carried outand identifies future guidance for improvements.

2ICSV14 • 9-12 July 2007 • Cairns • Australia

2. GOVERNING EQUATIONS
Finite Element Method provides calculation stability compared with Boundary ElementMethod, as it is not affected by the complexity of the model structure. The FEM is used todiscretize the bounded domain between the head model and the infinite boundary [5]. HRTF isa problem of solving scattering sound waves by a head model (scattering body) from a unit
source at a certain point in space. In order to visualize this, the Helmholtz equation including
non-homogeneous term creating the incident wave must be solved as an outgoing wave atinfinity.

Figure 2. Acoustic boundary domain∇2 p + k 2 p = gin Ω(5)∇p ⋅ n = 0on Γ(6)1
∂p+ikp = O( 2 )
r
∂r

(7)

Ω represents an unbounded domain and Γ its boundary. Forcing function is represented by g ,n is the outward unit normal vector on the boundary and r is the reference distance.
ACTRAN’s conjugated infinite element provides both high accuracy and flexibility inmodeling radiation problems. Based on multipole expansions in spherical or spheroidal or
ellipsoidal coordinates system, the infinite element can accommodate any geometry andarbitrary radial interpolation orders while remaining stable [6].

3. PROBLEM DESCRIPTION
Evaluation step for the reciprocal theorem were divided in two folds; case A had the sound source
inclined 45 degrees to the right at 0.5m with the recipient point being next to the right ear, whereas case
B had the exact opposite settings as illustrated below

Figure 3. Mesh structure and sound source/recipient point location

3

ICSV14 • 9-12 July 2007 • Cairns • Australia

Mesh comprised of 70000 nodes with 330000 tetrahedral elements, suitable for analysis up to 1000Hz
and the head model was obtained from MRI scan (shown below), which had a hight ofapproximately 0.25m with a 0.1m radius.

Figure 4. Head model and boundary condition of the room

To demonstrate the fact that the efficiency will not be lowered by the diffraction effect within
the computational domain, both cases were tested again with a reflective wall set on the left
hand side of the room. Simulations were carried out for frequencies between 50Hz to 1000Hz in
the goal of obtaining the impulse response and outline any numerical variation between thecases.

4. NUMERICAL COMPARISON
Efficiency of the reciprocal theory was first tested utilizing two separate cases, one by placing amonopole away from the head model having the receiving point next to the right ear and thesecond case exactly the opposite to the first. Comparison was carried out across a frequency
band of 50 to 1000Hz at 50Hz interval

Figure 5. Sound propagation patterns for case A and B at 1000Hz
It can clearly be observed that the pressure waves are generated from the monopole soundsource and the head diffracts its patterns. The overall sound strength within the room differs dueto the position of the sound source; one being close to the FE/IE boundary and the other close to
the head model. Although the outcome from a contour point of view differed, analysing real andimaginary figures at each recipient point respectively across the frequency range at 50Hz
interval were very closely matched as shown in the following page (figure 6). The following
step was to determine the level of deviation present in the numerical solution.

4ICSV14 • 9-12 July 2007 • Cairns • Australia

Figure 6. Real and Imaginary comparison between case A and B (Pa)
Quantitative error estimation of the result reiterated that error margin rises after reaching a certain
frequency but up to 600Hz (approximately 10 element per wavelength), the solution is fairly stable. This
relation is due to the fact that the grid has been created to contain 6 elements per wavelength andreaching the optimum frequency tends to affect the accuracy of the numerical solution as seen for both
amplitude and phase error margin in figure 7.

Amplitude error (eA ) =AmplitudeNo reciprocal − Amplitudewith reciprocal
AmplitudeNo reciprocal

Phase error (eP ) = ⎡⎣φNo reciprocal − φwith reciprocal ⎤⎦ *180 / π(8)(9)

Figure 7. Amplitude and Phase error estimate for non-reflective caseFor further verification of the stability of the solution, the left hand side of the room was set
as a reflective wall. Adding the extra boundary condition affected both real and imaginaryvalues as shown below (figure 8). A small irregularity is seen around the 800Hz region but as in
the previous non-reflective case, most of the points are identical.

Figure 8. Real and Imaginary comparison between case A and B with reflective wall (Pa)5ICSV14 • 9-12 July 2007 • Cairns • Australia

Error estimate study displayed a significant instability of the numerical solution at 800Hz for both
amplitude and phase values. Again, similar reasoning regarding element size was deducted for thisproblem, hence we believe refining the grid further or increasing the radial interpolation order could
diminish the error magnitude.

Figure 9. Amplitude and Phase error estimate for reflective wall case
The corresponding impulse response between the default and reflective wall case with thesound source inclined at 45 degrees toward the right are shown below for practical purpose. It isclear that the sound wave reaches the right ear first in the default case, as the source is closer toit compared to the left ear.

Figure 10. Impulse response for cases with sound source away from the head (Pa)The effect of the reflective wall is seen on the right graph as the second sound wave reflectedby the wall reaches the left ear first. Culmination of the results so far proves the validity of thereciprocal theory hence, the analysis moved on to higher frequency, which were deemed
impractical without heavy machine power.

5. HIGH FREQUENCY STUDY
The head model was placed inside a sphere of 0.17m radius with a monopole set next to the left
ear at a frequency range of 10 to 20 KHz. The element size had to be restrained to one sixth of
the wavelength (2.8mm) hence a dense grid of approximately 300000 nodes and 1.6million
elements had to be computed. As a reference, this case was operated on a machine with
Dual-Core Opteron 2218 (2.6GHz) processor with 8GB memory and took approximately 50
minutes for one frequency calculation. Directivity and the dip characteristics of sound pressurelevel at 1.5 meter radius around the head model (horizontal plane to the ears) were determined
across the relevant frequency. Figure 11 demonstrates how the FE/IE boundaries were divided
for the sound source and the recipient point, where the computational domain (FE) is kept to aminimum. IE governs the outer boundary and observation points were placed at a 5 degrees

6

ICSV14 • 9-12 July 2007 • Cairns • Australiainterval on the circle with 1.5 meter radius from the center coordinate (total of 72 observation
points).

Figure 11. Domain boundary and observation perimeterAs the monopole is placed near the left ear, only a weak portion of sound wave reaches theright hand side, which results in low sound pressure level compared to the left region being
illustrated in figure 12. Also, decreasing the wavelength increased the diffraction effect.

Figure 12. Sound pressure Level and Real Part contour at 20KHz
The directivity reiterated that the strongest sound region is constantly around the front left region
regardless of the frequency, thus the HRTF variation will always be higher between 180-360 degrees
compared to the 0-90 degrees plane.

Figure 13 Directivity and HRTF variation characteristics

7

ICSV14 • 9-12 July 2007 • Cairns • Australia

6. CONCLUSIONS AND FURTHER WORKThe reciprocal theory application was successful in terms of numerical accuracy and it was notaffected by the addition of extra reflective property. Also, the sound propagation pattern and theimpulse response suggested that the coupling with FEM produced logical results by eliminating
the potential frequency dependent calculation instability. Complexity of the analysis structure
was negligible (facial characteristics) and as unlimited monitor points (recipient points) could
be placed within a single calculation, determining the directivity is simple. Overall, thecoupling analysis managed to reduce computational time considerably by minimizing theanalysis domain without affecting the efficiency of the result. The next step would be tocompare simulation efficiency to that of experimental to further validate the feasibility ofFEM/IEM coupling with reciprocal theory. Calculation at higher order of frequency andaddition of extra boundary condition such as impedance could be of interest.

REFERENCES[1]

M. Otani and S. Ise, “A fast calculation method of the head-related transfer functions for multiple
source points based on the boundary element method”, Acoust. Sci. & Tech 24, 5 (2003).

[2]

Free-Field-Technologies-S.A., Actran 2006 Aeroacoustic Solutions: Actran/TM and Actran/LA User’s Manual , Belgium, 2006.

[3]

Y. Naka, A.A. Oberai and B.G Shinn-Cunningham “ The finite element method with theDirichlet-to Neumann map for sound-hard rectangular rooms” Proc. 18th Int. Congr. Acoustics,
Vol. IV, pp. 2477-2480 (2004).

[4]
[5]
[6]

M.S. Howe, Theory of Vortex Sound, CAMBRIDGE UNIVERSITY PRESS, 1992.

[7]

M. Otani and S. Ise, “Head Related Transfer Function Analysis and Application”, Journal ofAcoustical Society of Japan 61, 404-409 (2005)[8]B.F.G. Katz, “Boundary element method calculation of individual head-related transfer function. I.
Rigid model calculation”, Journal of Acoustical Society of America 110, 2440-2448 (2001)[9]

B.F.G. Katz, “Boundary element method calculation of individual head-related transfer function.
II. Impedance effects and comparison to real measurements”, Journal of Acoustical Society ofAmerica 110, 2449-2455 (2001)

F. Ihlenburg, Finite Element Analysis of Acoustic Scattering, Springer, 1998
R.J. Astley, G.J. Macaulay and J.P. Coyette, “Mapped wave envelope elements for acoustic
radiation and scattering”, JSV, 170 (1), 97-118 (1994)

8