E. Subjective experiment
VI. CONCLUSIONS
Two different numerical methods are compared in this work. One is the simple numerical method, virtual source representation. This method is motivated by the layer potential theory. The acoustic field of arbitrarily shaped radiators is described by using the principle of wave superposition in this method. The other is indirect boundary element method which is implemented by using a Vibro-Acoustic software-LMS SYSNOISE. The measurement of HRTFs is also implemented by using one simple measurement mechanism here. The goal of this investigation is to calculate the HRTFs of each azimuth and elevation efficiently. Both of these numerical methods can work and the calculated results which can catch the obvious characteristics of HRTFs. In the comparison between the numerical results and the measured HRTF databases, important characteristics of KEMAR dummy head with DB65 pinna can be captured by the two numerical methods in the frequency response
analysis. In directional response, a reasonable and clear direction was focused by the main lobe in all numerical results. Obvious ILD and ITD can also be found in the calculated HRIRs. In the process of HRTF computation, large efforts were paid in processing the ill-conditioned matrix in virtual source representation. By using the Tikhonov regularization in virtual source representation, the suitable regularization parameters can be selected. Appropriate regularization to the ill-conditioned matrix can make the numerical results close to the measured ones.
Although both the numerical results are close to the measured ones in IBEM and virtual source representation, the computation loadings are very different between these methods. The computation loadings in virtual source representation are much lighter than the IBEM in memory use and computation time. However, the errors made by using virtual source representation can still be improved by increasing the number of both virtual sources and mesh nodes and how to process the ill-conditioned matrix effectively will be still an important point in this method.
REFERENCES
1 Dolby Laboratories, Inc., “5.1-Channel production guidelines.”
(www.dolby.com)
2 J. W. Strutt (Lord Rayleigh), “On our perception of sound direction,” Philosophical Mag. 13, 214-232(1907).
3 B. Gardner and K. Martin, “HRTF measurements of KEMAR dummy-head microphone,” MIT Media Lab (1994).
( http://sound.media.mit.edu/KEMAR.html )
4 V. R. Algazi, R. O. Duda, D. M. Thompson and C. Avendano, “The CIPIC HRTF Database,” Proc. 2001 IEEE Workshop on Applications of Signal Processing to Audio and Electroacoustics, Mohonk Mountain House, New Paltz, NY, pp. 99-102, Oct. 21-24 (2001).
5 Suzuki Laboratory(Acoustic Information Systems), Tohoku University (2001).
( http://www.itakura.nuee.nagoya-u.ac.jp/ )
6 D. W. Batteau, “The role of the pinna in human localization,” Proc. Royal Society of London 168(B), 158-180(1967).
7 C. Phillip Brown and Richard O. Duda, “A Structural Model for Binaural Sound Synthesis”, IEEE TRANSACTIONS ON SPEECH AND SUDIO PROCESSING,, 6, 476-488, 1998
8 V. R. Algazi, R. O. Duda, R. Duraiswami, N. A. Gumerov, and Z. Tang,
“Approximation the head-related transfer function using simple geometric models of the head and torso”, J. Acoust .Soc. Am. 112, 2053-2064 (2002).
9 R. O. Duda, “3-D audio for HCI,” Department of Electrical Engineering, San Jose StateUniversity (2000).
( http://www.engr.sjsu.edu/~knapp/HCIROD3D/3D_home.htm )
10 Y. Kahana, “Numerical modeling of the head related transfer function”, a thesis
submitted for the degree of doctor philosophy, University of Southampton, Faculty of engineering and applied science, Institute of sound and vibration research, 2000.
11 P. A. Nelson, Y. Kahana, “Spherical harmonics, singular-value decomposition and the head-related transfer function”, J. Sound and Vib., 239, 607-637 (2001).
12 LMS SYSNOISE release notes and getting started manual.
13 LMS SYSNOISE user manual.
14 E. A. G. Shaw, “Transfrormation of sound pressure level from the free-field to the eardrum in the horizontal plane”, J. Acoust. Soc. of Am. 56, 1848-1861 (1974).
15 J. Chen, B. D. Van Veen, and K. E. Hecox, ”A Spatial feature extraction and regularization model for the head-related transfer function,” J. Acoust. Soc. Am. 97, 439-452 (1995).
16 黃嘉文, “Headphone Reproduction of 3D Immersive Spatial Sound Field”, 國立交 通大學機械工程研究所碩士論文.
17 Douglars D. Rife and John Vanderkooy, “Transfer-Function Measurement with Maximum-Length Sequences”, J. Audio Eng. Soc., 37, 419-443 (1989).
18 G. H. Koopmann, L. Song, and J. B. Fahnline, “A method for computinfg acoustic fields based on the principle of wave superposition”, J. Acoust. Soc. Am., 86, 2433-2438, (1989).
19 R. Jeans and I. C. Mathews, “The wave superposition method as a robust technique for computing acoustic fields”, J. Acoust. Soc. Am., 92, 1156-1166, (1992).
20 白明憲, 聲學理論與應用 : 主動式噪音控制, 修訂版 (全華書局, 2001).
21 L. Song, G. H. Koopmann, and J. B. Fahnline, “Numerical errors associated with the method of superposition for computing acoustic fields”, J. Acoust. Soc. Am. 89, June 1991.
22 MATLAB 6.5 Help contents, MathWorks. Inc.
23 A. Schuhmacher and J. Hald, K. B. Rasmussen, P. C. Hansen, “Sound source
reconstruction using inverse boundary element calculations”, J. Acoust. Soc. Am.
113, 114-127 (2003).
24 E. G. Williams, Fourier Acoustics, (Academic Press, 1999).
25 A. V. Oppenheim, R. W. Schafer, J. R. Buck, Discrete-time signal processing, second edition, (Prentice-Hall, 1999).
TABLES
Table I. Number of measurement and azimuth increments at each elevation
Elevation Number of
Table II. Comparison of computation loadings between IBEM and virtual source representation
Virtual source representation case 1 :
The nodes of head mesh are decimated into 2% of the original ones and the nodes of pinna mesh are also decimated into 17% of the original ones).
Virtual source representation case 2 :
The nodes of head mesh are decimated into 10% of the original ones and the nodes of pinna mesh are also decimated into 50% of the original ones).
Case
FIGURES
(a)
(b)
Fig. 1 Coordinate systems for HRTF representation
(a) Interaural-Polar coordinate system (b) Vertical-Polar coordinate system
X Y
Z
φ θ
θ φ
X Y
Z
(a)
(b)
Fig. 2 An HRTF example of KEMAR at azimuth=45∘and elevation=0∘
(a) The impulse response (HRIR). (b) The HRTF frequency response
Fig. 3 Schematic diagram of measurement mechanism
φ
θ
Loudspeaker
Turntable
Fig. 4 Photo of HRTF measurement mechanism
Fig. 5 Frequency response of measurement loudspeaker and power amplifier
Fig. 6 The program interface of HRTF measurement written by using Visual Basic
Fig. 7 The interface of Pulse multi analyzer software – Pulse Labshop
Fig. 8 Measured HRTFs post-processing flow chart Subjective measure
Free field measure Measured
Loudspeaker+Microphone Frequency Response
Measured HRTF
C om plex D ivision
Equalized HRTF HRIR Sound source
Sound source
Inverse FFT
Mic.
Fig. 9 Measured HRTFs on the horizontal plane (elevation = 0 )∘
(dash line: unequalized MIT HRTF, solid line : unequalized HRTF) (a) azimuth = 0 (b) azimuth = 30 (c) azimuth = 60∘ ∘ ∘
(d) azimuth = 90 (e) azimuth = 120 (f) azimuth = 150∘ ∘ ∘
Fig. 9 Measured HRTFs on the horizontal plane (elevation = 0 )∘
(dash line: unequalized MIT HRTF, solid line : unequalized HRTF) (g)azimuth = 180 (h) azimuth = 210 (i) azimuth = 240∘ ∘ ∘ (j)azimuth = 270 (k) azimuth = 300 (l) azimuth = 330∘ ∘ ∘
Fig. 10 Measured HRTFs on the horizontal plane (elevation = 0 )∘ (dash line: equalized MIT HRTF, solid line : equalized HRTF) (a) azimuth = 0 (b) azimuth = 30 (c) azimuth = 60∘ ∘ ∘ (d) azimuth = 90 (e) a∘ zimuth = 120 (f) azimuth = 150∘ ∘
Fig. 10 Measured HRTFs on the horizontal plane (elevation = 0 )∘ (dash line: equalized MIT HRTF, solid line : equalized HRTF) (g)azimuth = 180 (h) azimuth = 210 (i) azimuth = 240∘ ∘ ∘ (j)azimuth = 270 (k) azimuth = 30∘ 0 (l) azimuth = 330∘ ∘
Fig. 11 Indirect boundary element method definition (close body)
Fig. 12 Generic L-curve form; a plot in the log-log scale of the discrete smoothing norm versus the residual norm as a function of regularization parameter
log Lx
2less filtering
log Ax b -
2more filtering
Fig. 13 Spacing of nodes at maximum frequency. Seven points denote the six elements in one wavelength.
Field amplitude
Field variation in space 1
2 3
4
5 6
7
Fig. 14 Circular piston of radius α set in the rigid sphere. The piston vibrates with a velocityWg . The z axis passes through the center of the piston.
Baffle
Polar Cap
(θ) W w & = &
(θ) 0 w& =
α z
Fig. 15 The settings of the mesh elements in the SYSNOISE. The light color areas denote the rigid baffle and others denote the vibrating circular piston.
(a)
(b)
Fig. 16 The settings of mesh nodes in virtual source representation.
The plus sign in (a) and (b) denotes the vibrating circular piston.
Inner black points denote the virtual sources.
Fig. 17 Frequency responses of circular piston in a spherical baffle with optimum and smallτ in virtual source representation and analytical solution.
(a)
(b)
Fig. 18 L-curve of circular piston in a spherical baffle at
(a) 100 Hz. (b) Distribution of the value (γ ) of cost function.
The regularization parameter β is taken as the first singular value multiplied by the values marked on the curve.
Fig. 19 Frequency responses of circular piston in a spherical baffle calculated by using analytical solution, virtual source representation and IBEM.
(a)
(b)
Fig. 20 Directional response of circular piston in a spherical baffle at (a) 300 Hz (b) 500 Hz
(c)
(d)
Fig. 20 Directional response of circular piston in a spherical baffle at (c) 700 Hz (d) 1000 Hz
(a)
Fig. 21 The mesh settings of KEMAR dummy head using IBEM in the SYSNOISE (ITRI case).
(a) Head mesh given by ITRI
(b)
Fig. 21 The mesh setting of KEMAR dummy head using IBEM in the SYSNOISE (ITRI case).
(b) Pinna (White color areas denote the vibrating panels).
(a)
Fig. 22 The mesh settings of KEMAR dummy head using IBEM in the SYSNOISE (Kahana case).
(a) Head mesh given by Yuvi Kahana
(b)
Fig. 22 The mesh settings of KEMAR dummy head using IBEM in the SYSNOISE (Kahana case).
(b) Pinna (White areas are the vibrating panels).
(a)
(b)
Fig. 23 The distributions of KEMAR mesh nodes in virtual source representation.
(Only 2% of the original scanned head mesh nodes and 17% of the original scanned pinna mesh nodes).
(a)The right hand side of KEMAR (case 1) (b)The front of KEMAR (case 1)
(The cluster of light color points denote the virtual sources)
(c)
(d)
Fig. 23 The distributions of KEMAR mesh nodes in virtual source representation.
(Only 10% of the original scanned head mesh nodes and 50% of the original scanned pinna mesh nodes).
(c)The right hand side of KEMAR (case 2) (d)The front of KEMAR (case 2)
(The cluster of light color points denote the virtual sources)
(a)
(b)
Fig. 24 L-curve of HRTF computations using virtual source representation at (a) 100 Hz. (b) Distribution of value (γ ) of cost function.
The regularization parameter β is taken as the first singular value multiplied by the values marked on the curve. (Virtual source case 2)
(a)
(b)
Fig. 25 Calculated HRTFs on the horizontal plane by using IBEM (ITRI case).
(a) 0 azimuth degrees (b) 90 azimuth degrees
(c)
(d)
Fig. 25 Calculated HRTFs on the horizontal plane by using IBEM (ITRI case).
(c) 180 azimuth degrees (d) 270 azimuth degrees
(a)
(b)
Fig. 26 Calculated HRTFs on the horizontal plane by using IBEM (Kahana case).
(a) 0 azimuth degrees (b) 90 azimuth degrees
(c)
(d)
Fig. 26 Calculated HRTFs on the horizontal plane by using IBEM (Kahana case).
(c) 180 azimuth degrees (d) 270 azimuth degrees
(a)
(b)
Fig. 27 Calculated HRTFs on the horizontal plane by using virtual source representation (virtual source case 1).
(a) 0 azimuth degrees (b) 90 azimuth degrees
(c)
(d)
Fig. 27 Calculated HRTFs on the horizontal plane by using virtual source representation (virtual source case 1).
(c) 180 azimuth degrees (d) 270 azimuth degrees
(a)
(b)
Fig. 28 Calculated HRTFs on the horizontal plane by using virtual source representation (virtual source case 2).
(a) 0 azimuth degrees (b) 90 azimuth degrees
(c)
(d)
Fig. 28 Calculated HRTFs on the horizontal plane by using virtual source representation (virtual source case 2).
(c) 180 azimuth degrees (d) 270 azimuth degrees
(a)
(b)
Fig. 29 Directional response of measured HRTF on the horizontal plane at (a) 500 Hz (b) 2100 Hz
(c)
(d)
Fig. 29 Directional response of measured HRTF on the horizontal plane at (c) 5100 Hz (d) 9100 Hz
(a)
(b)
Fig. 30 Directional response calculated by using IBEM (Kahana case) at (a) 500 Hz (b) 2100 Hz
dash line: measured HRTF solid line: calculated HRTF (elevation = 0∘)
(c)
(d)
Fig. 30 Directional response calculated by using IBEM (Kahana case) at (c) 5100 Hz (d) 9100 Hz
dash line: measured HRTF solid line: calculated HRTF (elevation = 0∘)
(a)
(b)
Fig. 31 Directional response calculated by using IBEM (ITRI case) at (a) 500 Hz (b) 2100 Hz
dash line: measured HRTF, solid line: calculated HRTF (elevation = 0∘)
(c)
(d)
Fig. 31 Directional response calculated by using IBEM (ITRI case) at (c) 5100 Hz (d) 9100 Hz
dash line: measured HRTF, solid line: calculated HRTF (elevation = 0∘)
(a)
(b)
Fig. 32 Directional response on the horizontal plane calculated by using virtual source representation (virtual source case 1) at
(a) 500 Hz (b) 2100 Hz
(dash line: measured HRTF, solid line: calculated HRTF)
(c)
(d)
Fig. 32 Directional response on the horizontal plane calculated by using virtual source representation (virtual source case 1) at
(c) 5100 Hz (d) 9100 Hz
(dash line: measured HRTF, solid line: calculated HRTF)
(a)
(b)
Fig. 33 Directional response on the horizontal plane calculated by using virtual source representation (virtual source case 2) at
(a) 500 Hz (b) 2100 Hz
(dash line: measured HRTF, solid line: calculated HRTF)
(c)
(d)
Fig. 33 Directional response on the horizontal plane calculated by using virtual source representation (virtual source case 2) at
(c) 5100 Hz (d) 9100 Hz
(dash line: measured HRTF, solid line: calculated HRTF)
Fig. 34 Head related impulse response on the horizontal plane from the measured HRTF database.
(solid line: 90 azimuth degrees.) (dash line: 270 azimuth degrees.)
Fig. 35 Head related impulse response on the horizontal plane from the calculated HRTFs by using IBEM. (ITRI case).
Fig. 36 Head related impulse response on the horizontal plane from the calculated HRTFs by using virtual source representation (case 1).
Fig. 37 Head related impulse response on the horizontal plane from the calculated HRTFs by using virtual source representation (case 2).
(a)
(b)
Fig. 38 Sound localization results of an azimuth test using white noise (a) Measured HRTFs (b) HRTFs from MIT HRTF database
(c)
(d)
Fig. 38 Sound localization results of an azimuth test using white noise (c) Calculated HRTFs of virtual source case 1
(d) Calculated HRTFs by using IBEM (ITRI case)
Fig. 38 Sound localization results of an azimuth test using white noise (e) Calculated HRTFs of virtual source case 2
(a)
(b)
Fig. 39 The total average error and perceived angle error at each direction of sound localization
(a) Measured HRTFs
(b) HRTFs from MIT HRTF database
(c)
(d)
Fig. 39 The total average error and perceived angle error at each direction of sound localization
(c) Calculated HRTF of virtual source case 1 (d) Calculated HRTF by using IBEM (ITRI case)
(e)
Fig. 39 The total average error and perceived angle error at each direction of sound localization
(e) Calculated HRTF of virtual source case 2