COMSOL 6.4 - Head and Torso HRTF Computation

This tutorial model shows how to import a 3D scanned geometry of a human head and torso and compute the head related transfer function (HRTF). The scan is imported as a stl file and converted into a COMSOL geometry. The HRTF is computed using the reciprocity principle, locating the source at the ear canal entrance; this approach greatly reduces the computational cost to get a full 3D response. The acoustics are modeled using the Pressure Acoustics, Boundary Element interface of the Acoustics Module. The simulated results are compared with measured data from the actual subject and show good agreement.

The HRTF gives a complete description of the way the head and torso of an individual distorts incident sound fields. The HRTF is an important component of spatial hearing. The HRTF includes both so-called monaural and binaural cues. Binaural cues include the interaural time differences (ITD) and interaural level differences (ILD), whereas the monaural cues represent a spectral distortion of the sound that is identical for both ears, see Ref. 1. The HRTF is defined as the sound pressure level (SPL) measured at the eardrum (or the ear canal entrance as in this model) relative to the SPL when no head is present.

When virtual sound is used (or acoustic virtual reality), the HRTF is important in order to make the test subject experience a virtual sound scene. The HRTF can be measured, which can be a tedious task, or it can be simulated based on a scan of the individual. This model presents the latter approach on a scanned head geometry provided by the, Institute for Hearing Technology and Acoustics (IHTA), RWTH Aachen University, Germany, Ref. 2. The scan is of an actual individual where the facial features have been removed, while all the details of the ear geometry have been retained.

The scanned geometry (stl mesh) and measured data is with courtesy of the Teaching and Research Area of Medical Acoustics, Institute of Technical Acoustics, RWTH Aachen University, Germany. The stl mesh is licensed under the Creative Commons Attribution 4.0 International License and is provided “as is” with all warranties disclaimed, as stated in that license. See Ref. 7 and Ref. 8 for details.

Model Definition

A common approach when simulating the HRTF is to use the reciprocity principle; the source and receiver locations are reverted, Ref. 3. This means that, in the model, the source is located at the entrance of the ear canal and the evaluation is performed along a circle (or on a sphere for the full bubble) with its center in the middle of the head, between the ears. In this way, the HRTF can be deduced for all spatial directions for each frequency with just one simulation. Not using reciprocity requires solving one problem per incidence direction per frequency, which is not practical. Reciprocity is used for similar applications in Ref. 4 and Ref. 5.

The acoustic problem is modeled using the boundary element method (BEM) with the Pressure Acoustics, Boundary Element interface. This is especially efficient since the present model represents a pure radiation problem.

The imported stl mesh is depicted in Figure 1 and the COMSOL geometry generated from the stl mesh is depicted in Figure 2. Notice that the geometry has been moved and rotated (in Figure 2) to align the coordinate axis with the commonly used directions for directivity assessment. The evaluation circle for the HRTF that is used in the model, is represented in Figure 2. The evaluation is performed using a Radiation Pattern plot.

Figure 1: The imported stl mesh. The imported stl mesh is courtesy of the Teaching and Research Area of Medical Acoustics, Institute of Technical Acoustics, RWTH Aachen University, Germany. The stl mesh is licensed under the Creative Commons Attribution 4.0 International License and is provided “as is” with all warranties disclaimed, as stated in that license. See Ref. 7 and Ref. 8.

Figure 2: The generated COMSOL geometry from the stl mesh. Representation of the location for the HRTF evaluation on a circle in the horizontal plane.

Results and Discussion

The model is solved for the three frequencies where measurement data is provided (f = 1033.6 Hz, 2067.5 Hz, and 3962.1 Hz) as well as two additional higher frequencies (f = 6000 Hz and 8000 Hz).

The pressure field generated from the excitation at the ear canal entrance is depicted in Figure 3 for the frequencies f = 1033.6 Hz, 2067.5 Hz, and 3962.1 Hz. These have been selected as they coincide with the measurement data (octave band center frequencies). A unit normal velocity is assigned to the ear canal entrance.

The acoustic pressure is also depicted in a cross-section plane in Figure 4 and the corresponding sound pressure level (SPL) is depicted in Figure 5. Both are evaluated at the 4 kHz octave band center frequency. The SPL plot clearly shows the presence of notches (cancellations) for certain directions. These are more evident at the higher frequencies.

Figure 3: Acoustic pressure at the surface of the head and torso evaluated at three frequencies.

Figure 4: Acoustic pressure on the head and torso and in a cut plane.

Figure 5: Sound pressure level on the surface of the head and torso and in a cut plane.

The HRTF in the horizontal plane (the xy-plane) is depicted for the five evaluation frequencies in Figure 6. In the plot, the HRTF is normalized to 0 dB toward the front (polar angle θ = 0°). In the following three plots — Figure 7, Figure 8, and Figure 9 — the computed HRTFs are compared to measurement data; see Ref. 6 and Ref. 8. In the plots, the HRTF data has been rotated by θ0 = 4.5° to make the location of the notches match (defined by the parameter theta0).

Notice that the model assumes that the head and torso is located in free space. This is consistent with the measured data where the floor reflections have been removed. A time-windowed approach is used, so that the reflections will not affect the data. In the measurements, there is a maximum length of the impulse response of about 330 samples at a 44,100 Hz sampling rate.

The COMSOL models results agree well with the measurement data. There are some general small discrepancies; these can be due to head movement during the measurements. In the current measurement setup at RWTH, head movement is tracked and compensated for. A larger discrepancy can be seen in the 1 kHz plot toward 30° (in Figure 7). Shoulder reflections are typically seen at around 1.5 kHz so a slight under- or overestimate of the shoulder size, when generating the head scan, could introduce this error.

Figure 6: Comparison of the normalized HRTF evaluated at the three frequencies.

Figure 7: Comparison of the simulated HRTF with measured data at 1033 Hz.

Figure 8: Comparison of the simulated HRTF with measured data at 2066 Hz.

Figure 9: Comparison of the simulated HRTF with measured data at 3962 Hz.

References

1. M. Vorländer, Auralization, Springer, 2008.

2. Web link: www.akustik.rwth-aachen.de/

3. A.D. Pierce, Acoustics: An Introduction to Its Physical Principles and Applications, Acoustical Society of America, 1991.

4. Z. Conrad, “Hats Off to the Boundary Element Method,” IEEE Spectrum Multiphysics Simulation, October 2018, p. 30, web link: www.comsol.com/offers/multiphysics-simulation-2018.

5. M.H. Jensen, “Improving the Performance of Hearing Aids Using Acoustic Simulations,” COMSOL Conference 2009, web link: www.comsol.com/paper/improving-the-performance-of-hearing-aids-using-acoustic-simulations-7227.

6. Web link: www.akustik.rwth-aachen.de/cms/Technische-Akustik/Forschung/~lxfd/Downloads/lidx/1/.

7. H.S. Braren and J. Fels, “A High-Resolution Individual 3D Adult Head and Torso Model for HRTF Simulation and Validation: 3D Data,” web link (DOI): doi.org/10.18154/RWTH-2020-06760

8. H.S. Braren and J. Fels, “A High-Resolution Individual 3D Adult Head and Torso Model for HRTF Simulation and Validation: HRTF Measurement,” web link (DOI): doi.org/10.18154/RWTH-2020-06761

Acoustics_Module/Tutorials,_Pressure_Acoustics/head_torso_hrtf

Modeling Instructions

From the menu, choose .

New

In the window, click

Model Wizard

In the window, click

In the tree, select > > .

Click .

Click

In the tree, select > .

Click

Global Definitions

Parameters 1

In the window, under click .

In the window for , locate the section.

Click

Browse to the model’s Application Libraries folder and double-click the file head_torso_hrtf_parameters.txt.

Interpolation 1 (int1)

In the toolbar, click

and choose > .

In the window for , locate the section.

From the list, choose .

Click

Browse to the model’s Application Libraries folder and double-click the file head_torso_hrtf_measured.txt.

Locate the section. In the table, enter the following settings:


Columns	Type	Settings
Column 2	Function values	Function name=int1a
Column 3	Function values	Function name=int1b
Column 4	Function values	Function name=int1
Column 5	Function values	Function name=int1c
Column 6	Function values	Function name=int1d
Column 7	Function values	Function name=int1e

In the table, click to select the cell at row number 1 and column number 1.

In the text field, type rad.

In the table, click to select the cell at row number 2 and column number 1.

In the text field, type HRTF_1033_real.

In the text field, type Pa.

In the table, click to select the cell at row number 3 and column number 1.

In the text field, type HRTF_1033_imag.

In the text field, type Pa.

In the table, click to select the cell at row number 4 and column number 1.

In the text field, type HRTF_2067_real.

In the text field, type Pa.

In the table, click to select the cell at row number 5 and column number 1.

In the text field, type HRTF_2067_imag.

In the text field, type Pa.

In the table, click to select the cell at row number 6 and column number 1.

In the text field, type HRTF_3962_real.

In the text field, type Pa.

In the table, click to select the cell at row number 7 and column number 1.

In the text field, type HRTF_3962_imag.

In the text field, type Pa.

Locate the section. From the list, choose .

From the list, choose .

Locate the section. Click

Analytic 1 (an1)

In the toolbar, click

and choose > .

In the window for , type p1033 in the text field.

Locate the section. In the text field, type HRTF_1033_real(theta+theta0)+i*HRTF_1033_imag(theta+theta0).

In the text field, type theta.

Click to expand the section. Select the checkbox.

In the text field, type 2*pi.

Locate the section. In the table, enter the following settings:


Argument	Unit
theta	rad

In the text field, type Pa.

Click to expand the section. Select the checkbox.

Locate the section. In the table, enter the following settings:


Plot	Argument	Lower limit	Upper limit	Fixed value	Unit
√	theta	-pi	3*pi	0	rad

Click

Analytic 2 (p2)

Right-click and choose .

In the window for , type p2067 in the text field.

Locate the section. In the text field, type HRTF_2067_real(theta+theta0)+i*HRTF_2067_imag(theta+theta0).

Click

Analytic 3 (p3)

Right-click and choose .

In the window for , type p3962 in the text field.

Locate the section. In the text field, type HRTF_3962_real(theta+theta0)+i*HRTF_3962_imag(theta+theta0).

Click

In the window, right-click and choose > .

Mesh Part 1

In the window for , locate the section.

Select the checkbox.

From the list, choose .

Import 1

In the window, under > > click .

In the window for , locate the section.

Click

Browse to the model’s Application Libraries folder and double-click the file head_torso_hrtf_scan.stl.

From the list, choose .

Locate the section. In the text field, type 180.

Add a attribute to mirror the STL mesh in the xy-plane.

Transform 1

In the toolbar, click

and choose .

In the window for , locate the section.

From the list, choose .

In the text field, type -1.

Click

Now, cut the surface mesh (using a cylinder) in order to create selections for the entrance of the ear canal.

Partition with Cylinder 1

In the toolbar, click

and choose .

In the window for , locate the section.

In the text field, type 2.8.

Locate the section. In the text field, type -1.3.

In the text field, type 0.6.

Locate the section. From the list, choose .

Click

Intersect with Plane 1

In the toolbar, click

and choose .

In the window for , locate the section.

From the list, choose .

In row , set to -95.

In row , set to 10.

In row , set to -100.

In row , set to -95.

Leave all other components at their default zero value.

Click

Finalize

In the window, right-click and choose .

The final mesh part created from the imported .stl file, of the scanned head and torso, should look like the image in Figure 1. Use the mouse to rotate, zoom, and move the geometry in the graphics window.

Geometry 1

In the window, expand the > node, then click .

Import 1 (imp1)

In the toolbar, click

In the window for , locate the section.

From the list, choose .

Click

Form Union (fin)

In the window, click .

In the window for , click

Remove Details 1 (rmd1)

In the toolbar, click

and choose .

In the window for , click

In the section the number of details removed can be seen.

Use the mouse to rotate, zoom, and move the geometry to see it from the front. The geometry should look like the image in Figure 2.

Add Material

In the toolbar, click

to open the window.

Go to the window.

In the tree, select > .

Click the button in the window toolbar.

In the toolbar, click

to close the window.

Materials

Air (mat1)

In the window for , locate the section.

From the list, choose .

Pressure Acoustics, Boundary Elements (pabe)

In the window, under click .

In the window for , locate the section.

From the list, choose .

Normal Velocity 1

In the toolbar, click

and choose .

Select Boundary 6 only.

In the window for , locate the section.

In the vn text field, type 1.

Mesh 1

Free Triangular 1

In the toolbar, click

and choose .

Size

In the window, click .

In the window for , locate the section.

Click the button.

Locate the section. In the text field, type min(20[mm],lam0/4).

In the text field, type 3[mm].

In the text field, type 2.

Free Triangular 1

In the window, click .

In the window for , locate the section.

From the list, choose .

Size 1

Right-click and choose .

In the window for , locate the section.

Click the button.

Locate the section.

Select the checkbox. In the associated text field, type lam0/4.

Locate the section. In the list box, select .

Click

Select Boundaries 1 and 2 only.

In the window, right-click and choose .

The mesh should look like the image below, here meshed to resolve a frequency of 4000 Hz. You can change the parameter f0 and build the mesh again, to see how it looks at different frequencies.

Definitions (comp1)

Before solving the model, add a variable theta that defines the horizontal polar angle. The variable is used when postprocessing the measured HRTF data.

Variables 1

In the window, expand the > node.

Right-click > and choose .

In the window for , locate the section.

In the table, enter the following settings:


Name	Expression	Unit	Description
theta	atan2(y,x)	rad	Horizontal polar angle

Study 1

Step 1: Frequency Domain

In the window, under click .

In the window for , locate the section.

In the text field, type f0.

Parametric Sweep

In the toolbar, click

In the window for , locate the section.

Click

Solve the model for the three frequencies used for comparison to the measurements and two additional higher frequencies.

In the table, enter the following settings:


Parameter name	Parameter value list	Parameter unit
f0 (Frequency)	{1033.6, 2067.2, 3962.1, 6000, 8000}	Hz

In the toolbar, click

Results

Grid 3D 1

In the window, expand the > node, then click .

In the window for , locate the section.

Find the subsection. In the text field, type -0.3.

In the text field, type 0.3.

Find the subsection. In the text field, type -0.5.

In the text field, type 0.5.

Find the subsection. In the text field, type -0.5.

In the text field, type 0.5.

Click to expand the section. Now, set the resolution of the evaluation grid. This resolution should be adequate to resolve the wave pattern being visualized (the wavelength). If the resolution is too poor the results will appear wrong.

In the text field, type 40.

In the text field, type 120.