PROCEEDINGS of the
23rd International Congress on Acoustics
9 to 13 September 2019 in Aachen, Germany

Independent modes and dimensionality reduction of head-relatedtransfer functions based on tensor decomposition
Tong ZHAO1; Bosun XIE2
1

Acoustic lab, School of Physics and Optoelectronics, South China University of Technology, Guangzhou, China

Acoustic lab, School of Physics and Optoelectronics, South China University of Technology, Guangzhou, China

ABSTRACT
Head-related transfer functions (HRTFs) are acoustical transfer functions from a source to two ears in thefree-field and vital to spatial audio. Generally, far-field HRTFs vary with source direction, frequency andsubject. Accordingly, the dimensionality of a full set of HRTFs is very large. Linear decomposition has beenused to reduce the dimensionality of HRTFs. Due to the multiple variable-dependent characteristics of HRTFs,
in present work, tensor decomposition is applied to analyze the independent modes of HRTF variation so as
to achieve a more efficient representation of them. The measured HRTF magnitudes are first smoothed witha frequency bandwidth of 1 ERB. Then the smoothed HRTF magnitudes of 52 Chinese subject, 493 directions
for each subjects and 42 discrete frequencies (up to 14 kHz) are combined as a 3rd-order tensor and thenundergone a Tucker decomposition. The results indicate that 17 subject-related modes, 13 directional-relatedmodes and 7 frequency-related modes account for more than 95% variations of HRTFs at corresponding
dimension respectively, and thus the HRTF dataset can be reduced by truncating the modes to correspondingorders. Based on the results, it is also deduced that 17 anthropometric parameters at least are needed for theanthropometry-based HRTF customization.
Keywords: HRTF, Tensor decomposition, Dimensionality reduction

1. INTRODUCTION
Head-related transfer functions (HRTFs) are the acoustical transfer functions from a sound source
to two ears in the free field, which characterize the comprehensive filtering effect s of anatomical
structures on sound waves (1, 2). HRTFs contain sound source localization information such as
binaural time difference (ITD), binaural sound level difference (ILD), and high -frequency spectral
cue. They play a very important role on the researches of binaural hearing and applications in virtual
auditory display (VAD).
Generally, at the far-field source distance r > 1.0 m with respect to head centre, HRTFs vary ascomplicated functions of source direction and frequency. Due to the differences in the anthropometricgeometries of head, pinnae and torso etc., HRTFs also depend on individual. Due to the multiple
variable-dependent characteristics of HRTFs, the dimensionality of a full set of HRTFs is very large.
In researches and applications, it is often required to reduce the dimensionality of HRTFs. For
example, an accurate VAD needs individualized or matched HRTFs (3). Measurement and numericalcalculation are two common methods to acquire HRTFs. However, it is difficult to acquire a full setof HRTFs for each individual by measurement or numerical calculation. Alternatively, a full set ofmatched HRTFs of a new subject can be customized from an existing HRTF database andanthropometric-based measurement or from a few HRTF measurements on the ne w subject (4, 5).
Prior to customization, dimensionality reduction of the existing HRTF database is needed. In addition,a low-dimensional representation of HRTFs reduces the required memory space in practical VAD.
Actually, there exist some similarities or correlations among the HRTFs of different directions,frequencies and subjects. The dimensionality of HRTFs can be reduced by removing these correlations.
1
2

5212

371870906@qq.com
Corresponding to: phbsxie@scut.edu.cn

Linear decomposition is often used to derive the low-dimensional representation of HRTFs. Principal
component analysis (PCA), or equivalently, singular value decomposition (SVD), is a common
example of linear decomposition (5-7). By using PCA, HRTFs are decomposed into a weighted sumof spectral shape basis functions or spatial basis functions. Each ba sis function represents an
independent mode of HRTF variation in the frequency domain or spatial domain. The low -dimensional
representation of HRTFs is derived by truncating the basis functions decomposition to some precedingorders that represent the most important variance of HRTFs.
Tensor is a multi-dimensional extension of conventional vector and matrix (8). Tensor
decomposition is a multilinear modelling technique and an extension of traditional PCA or SVD
calculation. Due to the multiple variable-dependent characteristics of HRTFs, a more efficient
representation of HRTFs can be achieved by tensor decomposition. One advantage of tensor
decomposition is that the interactions between multiple variables are considered. Grindlay and
Vasilescu first presented a tensor framework for HRTF decomposition and used it to customize
individualized HRTFs from anthropometric measurement (9). Huang and Li also explored theindividualized HRTF customization using tensor decomposition (10). However, the independent
dimensionalities or modes required for HRTF reconstruction were not analyzed in detail in
aforementioned works, especially for subject-related modes.
In present work, tensor decomposition is further applied to analyse the independent modes of HRTF
variation, so as to achieve a more efficient representation of the original HRTF datasets.

2. PRINCIPLE of Tensor decomposition of HRTFs
In present work, spatial direction is specified by elevation -90º ≤ Ф ≤90º and azimuth 0º ≤ θ < 360º,
where Ф = -90°, 0° and 90° represents the bottom, horizontal and top direction, respectively; In thehorizontal plane, θ = 0°, 90° and 180° represents the front, right and back direction, respectively.
Suppose that there are a full set of HRTFs from S subjects, M directions for each subject, and N
discrete frequencies at each direction. The data can be HRTF magnitudes or complex -valued HRTFs
with magnitude and phase. Then the full set of HRTFs are represented by a S × M × N tensor H SuM uN
with entries H(s, m, n) =Hsmn . Where s =1,2…S, m = 1,2…M and n = 1,2…N represent the indexesof subject, direction and frequency, respectively. To decompose the HRTF efficiently, the mean HRTF
across S subjects and M directions are subtracted from the entries of H SuM uN , yielding a zero-mean
HRTF tensor H ' S uM u N with entries:

1 S M¦ ¦ H (s' , m' , n)(1)
SM s ' 1 m ' 1
s 1,2...S , m 1,2...M , n 1,2...NThe Tucker decomposition of tensor H ' S uM u N is expressed as the product of a core tensor and three
unitary matrixes (8):(2)
H'SuM uN WSuM uN us U SsuS um U Mm uM un U Nn uN
H ' sm n

H ( s, m, n) 

Where WS uM u N is a S × M × N core tensor with entries Wsmn , s = 1,2…S, m = 1,2,…M, and n = 1,2…N.
n
s
n
, m
and U Nu
are three unitary matrixes with entries U sss ' , U mmm' and U nn
respectively.
U Su
N
S U M uM'
Subscript in each matrix indicates the dimensionality of matrix. These unitary matrixes satisfy
following orthonormal properties:(3)
U sSuS U Ss uS I S uS
U mMuM U mM uM I M uM
U nNu N U nN u N I N u N

Where, the superscript “+” denotes the transpose and complex conjugation of the matrix. I is identity
matrix.
The products of core tensor and three unitary matrixes are calculated so that the entries of tensor
H ' S uM u N
are related to the entries of core tensor WS uM u N and entries of three unitary matrixes asfollows:

H 'sm nS

MN

¦ ¦ ¦Wn
U sss 'U mmm'U nn'

s 'm 'n '(4)s' 1 m' 1 m' 1s, an S × (MN) matrix
Similar to the procedures on PCA (2), to find the unitarity matrix U Su
S
1
1
is constructed. The rows of H Su
represent the mean-subtracted HRTFs for different
H Su
(MN )
(MN )

5213

1
subjects, and the columns of H Su
represent the mean-subtracted HRTFs at various directions and
(MN )
1
1
frequencies. An S×S Hermitian matrix RSu
is constructed from H Su
as,
(MN )
S

1

(5)
H S1u( MN ) H (1MN)uSNMs
The preceding S’ columns of unitarity matrix U Su
are constructed from the orthonormal eigenvectors
S
1of matrix RSuS associated with the S’ ≤ S positive eigenvalues in a descending order:
RS1uS

O1s ! Os2 ! ... ! OsS ' ! 0
Similarly, to find the unitarity matrix Um
M uM, an M × (SN) matrix H(6)2
M u(SN )is constructed. The rows

represent the mean-subtracted HRTFs at different directions, and the columns of H M2 u(SN )represent the mean-subtracted HRTFs for different subjects and at various frequencies. An M×M
Hermitian matrix RM2 uM is constructed from H M2 u(SN ) as,of H

2
M u(SN )1(7)
H M2 u( SN ) H (2SN )uM
SN
The preceding M’ columns of unitarity matrix U Mm uM are constructed from the orthonormal
eigenvectors of matrix RM2 uM associated with the M’ ≤ M positive eigenvalues in a descending order.
n3
Finally, to find the unitarity matrix U Nu
, an N × (SM) matrix H Nu
is constructed. The rows
(SM )N
RM2 uM3
3
of H Nu
represent the mean-subtracted HRTFs at various frequencies, and the column of H Nu
(SM )
(SM )
represents the mean-subtracted HRTFs for different subjects and at various directions. An N×N
3
3
is constructed from H Nu
as,Hermitian matrix RNu
(SM )
N

1

(8)
H N3 u( SM ) H (2SM
)u N
SM
n
The preceding N’ columns of unitarity matrix U Nu
are constructed from the orthonormalN3eigenvectors of matrix RNuN associated with the N’ ≤ N positive eigenvalues in a descending order.
Once the three unitarity matrixes are found, using the orthonormal properties of Eq.(3), the coretensor for decomposition is calculated by:(9)
WSuM uN H 'S uM uN us USs uS um UmMuM un UnNuN
or
RN3 uN

Ws ' m ' n '

S

M

N

¦¦¦ H 'sm nn*
U sss*'U mm*m'U nn'(10)

s 1 m 1 n 1

The superscript “*” denotes complex conjugation.
The nature of tensor decomposition in Eq.(2) is to decompose the mean-subtracted HRTFs into ascombination of modes of independent variables. The preceding S’ columns of matrix U Su
representSsS’ ≤ S subject-related modes, each row of matrix U Su
represents a specified subject. The precedingS
M’ columns of matrix U Mm uM represent M’ ≤ M direction-related modes, each row of matrix U Mm uMnrepresents a specified direction. The preceding N’ columns of matrix U Nu
represent N’ ≤ NNmfrequency-related modes, each row of matrix U NuN represents a specified frequency. The orthogonalcondition given by Eq.(3) indicates that each mode is independent. In addition, the core tensor
WSuM uN is only relevant to the modes of subject, direction and frequency, but independent from thevariations of direction, frequency and subject.
When the mean-subtracted HRTFs are approximately reconstructed from the modes or eigenvectors
nsassociated with the preceding eigenvaluesˈthe columns or modes of matrixes U Su
, U Mm uM and U NuNSsˆare truncated to S” < S’ ≤ S, M” < M’ ≤ M, and N” < N’ ≤ N, respectively, yielding S × S” matrix U,
S uS "

M × M” matrix Uˆ Mm uM " and N × N” matrix U Nn uN " . Then, the mean-subtracted HRTFs can beapproximated as:

Hˆ 'SuM uN | Wˆ S "u M "u N '' us Uˆ s SuS " um Uˆ m M uM " un Uˆ n N uN "

5214

(11)

Where Wˆ S "u M "u N " is a 3-order core tensor with dimensionality of S” × M” × N”, which is also
constructed from S× M × N tensor WS uM u N by truncating it to corresponding dimensionality.
Eq.(11) can also be written as:
S" M " N"nHˆ 'sm n | ¦¦¦ Ws 'm 'n 'U sss 'U mmm'U nn'(12)s' 1 m' 1 n' 1

To reconstruct the complete HRTFs, the mean in Eq.(1) should be supplemented into Eq.(11) orEq.(12).
Similar to the case of PCA (2, 5, 6), the cumulative percentage variance of the energy represented
by the preceding modes is used to evaluate the accuracy or error in HRTF reconstruction. For subjectrelated modes, it is calculated by corresponding eigenvalues:
S"

Ks

¦O

s

¦O

s

s 1
S'

(13)

u 100%

s 1

Similar equations are available to the direction-related and frequency-related modes, respectively.
The error in reconstruction can also be evaluated by relative error er(s, m, n).
| H ( s, m, n)  Hˆ ( s, m, n) |2
(14)
(dB )
er ( s, m, n) 10 log 10| H ( s, m, n) |2
Where H ( s, m, n) and Hˆ ( s, m, n) are original and reconstructed HRTFs, respectively.
It can also be evaluated by the mean relative error er av1 (m, n) across subjects:
1 S
erav1 (m, n)¦ er (s, m, n) (dB)S s1and the mean relative error er av2 across subjects, directions and frequencies:
1 S M Nerav2¦¦¦ er (s, m, n) (dB)SMN s 1 m 1 N 1

(15)

(16)

3. HRTF data and pre-processing
The aforementioned tensor decomposition method is applicable to complex-valued HRTFs, linear
or logarithmic HRTF magnitudes or even head-related impulse responses (HRIRs) in the time domain.
Here, the tensor decomposition method is applied to linear HRTF magnitudes.
The measured HRTF database used for analysis contains 52 Chinese subjects (26 male and 26
female) and 493 source directions at source distance of 1.5 m for each subject (11). The measured
directions are uniformly distributed in each elevation. Table 1 lists the number of measured directions
and azimuthal intervals in each elevation. The data are binaural HRIRs measured using the blocked ear-canal technique (12). The length of each HRIR is 512 points at the sampling frequency of 44.1
kHz.
Table 1 – Numbers and intervals of azimuthal HRIR measurement in each elevation
Elevation Ф

-30°

-15°

0°

15°

30°

45°

60°

75°

90°

72

72

72

72

72

72

36

24

1

5°

5°

5°

5°

5°

5°

10°

15°

360°

Number of
azimuthalmeasurement
Interval of azimuthal
measurement Δθ

The measured HRIRs are first padded zeros, and undergone discrete time Fourier transform (DFT),
yielding HRTF magnitudes of 8192 points. To adapt the auditory resolution of frequency, the 8192 point HRTF magnitudes are smoothed above 500 Hz by using a moving frequency window with abandwidth of 1 ERB (below 500 Hz, smoothing is not needed), yielding the 128 -point HRTF

5215

magnitude in the linear frequency scale. Due to the conjugate symmetry of HRTFs along discrete
frequency axes, HRTF magnitudes at the preceding 42 discrete frequencies in the linear scale are used
for the tensor decomposition. Preceding 42 discrete frequencies correspond to a frequency range from
0 Hz to 14 kHz at a sampling frequency of 44.1 kHz and thus include most important localization cuesand other perceivable attributes. Moreover, because the measured HRTFs exhibit larger measured
error above 14 kHz, the data at this frequency range should not be included in the analysis.
Therefore, aforementioned pre-processing yields a tensor of HRTF magnitudes for S = 52 subjects,
M = 493 directions and N = 42 discrete frequencies, which is used as the original HRTF magnitudes
for subsequent tensor decomposition.

4. RESULTS
Because the results of left-ear HRTFs and right-ear HRTFs are similar, the results of left-ear HRTFs
are mainly discussed here. Table 2 lists the relationship between the cumulative percentage varianceof energy (Var-s) and the number S” of subject-related modes (basis vectors). As seen, the Var-s
increases with the number S” of modes. The first mode accounts for nearly 80% of the subject-related
variance. Therefore, there exists basic similarity among the HRTF magnitudes of different subjects.
The contributions of other modes descend with the order. The preceding 7 and 17 modes a ccount for
approximately 90 % and 95 % of subject-related variation, respectively.
Table 2 – Cumulative percentage variance of energy for various numbers of subject-related modes (left-ear
HRTFs)
S”

1

2

4

7

17

31

39

Var-s (%)

79.7

82.4

86.9

90.2

95.2

98.1

99.1

Similarly, Table 3 lists the relationship between the cumulative percentage variance of energy (Var m) and the number M” of subject-related modes. Var-m also increases with the number M” of modes.
The first mode accounts for nearly 70% of the direction-related variance. Therefore, there exists basic
similarity among the HRTF magnitudes of different directions. The contribu tions of other modes
descend with the order. The preceding 6 and 13 modes account for approximate 90 % and 95 % of
directions-related variation, respectively.
Table 3 – Cumulative percentage variance of energy for various numbers of direction-related modes (leftear HRTFs)
M”

1

3

4

6

13

28

49

Var-m(%)

69.5

81.3

85.9

90.1

95.3

98.1

99.0

Table 4 lists the relationship between the cumulative percentage variance of energy (Var-n) and thenumber N” of frequency-related modes. Var-n also increases with the number N” of modes. The firstmode accounts for approximate 73% of the frequency-related variance. Therefore, there exists
correlation among the HRTF magnitudes of different frequencies. The contributions of other modes
descend with the order. The preceding 4 and 7 modes account for approximate 90 % and 95 %directions-related variation, respectively.
Table 4 – Cumulative percentage variance of energy for various numbers of frequency-related modes (leftear HRTFs)

5216

N”

1

2

3

4

7

12

16

Var-n(%)

73.3

80.9

87.1

90.2

95.4

98.1

99.1

Figure1 – Mean erav1(m, n) across subjects for the left-ear HRTF magnitudes in the horizontal plane
Calculations indicate that, for construction with S” = 17 subjected-related modes, the mean relative
error er av2 across subjects, directions and frequencies in Eq. (16) are -14.5 dB for the left-ear HRTF
magnitudes and -13.6 dB for the right-ear HRTF magnitudes. In most cases, the mean relative error
er av1(m, n) across subjects is less than -10 dB, which is basically sufficient for the HRTF reconstruction.
Errors mainly occur at high frequencies and source directions contralateral to the ear. Actually, theattenuation and interference caused by head shadow reduce the contralateral HRTF magnitudes(therefore reduce the signal-to-noise ratio) and make them variate complicatedly at high frequency,
which makes the reconstruction of contralateral HRTF magnitudes difficult. However, thecontributions of contralateral HRTFs to localization are less than those of ipsilateral HRTFs (13), theerrors in contralateral HRTF reconstruction are not so significant as those of ipsilateral HRTFs.
As an example, Figure 1 illustrate the mean relative error er av1(m, n) across subjects for the leftear HRTF magnitudes in 72 horizontal direction (Ф = 0°) and up to frequency of 14 kHz. Calculationfrom the results in Figure 1 indicate that the mean erav1(m, n) across frequencies are -19.1 dB at theipsilateral direction (θ = 265°, Ф = 0°) and -5.1 dB at the contralateral direction (θ = 85º, Ф = 0°).

5. Discussion
The results in Section 4 indicate that 4 and 7 frequency-related modes are needed to represent 90 %and 95 % frequency-related variation of HRTF magnitudes, respectively; 6 and 13 direction-relatedmodes are needed to represent 90 % and 95 % directions-related variation of HRTF magnitudes.
Kistler and Wightman applied PCA to the logarithmic HRTF magnitudes of two ears for S = 10 humansubjects at M = 265 directions and N = 150 discrete frequencies (ranging from 0.2 to 15 kHz). The
results indicated that 5 spectral shape basis functions or frequency-related modes were needed torepresent 90% frequency-related variation in logarithmic HRTF magnitudes (6). Xie applied spatial
principal component analysis (SPCA) to the HRTF magnitudes of two ears for S = 52 human subjects
at M = 493 source directions and N = 59 discrete frequencies (from DC to 20 kHz) (5). The results
indicated that 16 spatial basis functions (direction-related modes) were needed to represent more than
95 % direction-related variation in HRTF magnitudes. Therefore, the results in present work arecomparable to previous works.
Xie et al. analyzed the correlations among HRTF magnitudes of different subjects and applied acluster analysis to them (14). The HRTF magnitudes used in that analysis were identical to those usedin present work. The results indicated that HRTF magnitudes for most of 52 subjects could beclassified into seven clusters. Some singleton clusters were also observed on a few subjects, whichreflected the diversity in individualized HRTF magnitudes. The results in Section 4 also indicate that
7 subject-related modes are needed to represent 90 % subject-related variations in HRTF magnitudes.
To represent the diversity in HRTF magnitudes more accurately, more subject-related modes areneeded. For example, 17 subject-related modes are needed to represent 95 % subject-related variations
in HRTF magnitudes. Therefore, the results of subject-related modes analysis in present work are also
comparable to previous works.
It should be noted that the subject-related variations of HRTF magnitudes from a Chinese-based
database are analyzed in present work. Due to the statistical difference among the HRTF magnitud esof different nationalities (such as Chinese and western subjects), it is expected that more subject -5217

related modes are needed to represent the diversity of HRTF database including subjects of multiple
nationalities. From this point of view, it is beneficial to analyze the subject-related modes for the
HRTFs of each special nationality separately.
In anthropometry-based HRTF customization, individualized HRTFs are approximately estimatedfrom anthropometric measurements by linear or nonlinear mapping technique (2, 4). The mapping
methods are derived from the statistical analysis on a baseline HRTF and anthropometric measurement
database with sufficient number of samples. To make the mapping efficiently, decomposition anddimensionality reduction of HRTFs are needed. Because 17 subject-related modes are needed torepresent 95 % subject-related variations in HRTF magnitudes, it is deduced that 17 anthropometricparameters at least are needed for the mapping in anthropometry-based HRTF customization. Grindlay
and Vasilescu suggested a tensor decomposition of HRTFs and anthropometry-based HRTF
customization method (9). They chose 5 subject-related modes to represent the subject-related
variations in HRTF, which seems not enough.

6. CONCLUSIONS
Far-field HRTFs vary with subject, direction and frequency. Tensor decomposition is an effect tool
for analyzing the independent modes of HRTF variation and reducing the dimensionality of HRTFs.
Analysis on the HRTF magnitudes of 52 Chinese subjects, 493 directions and up to a frequency of 14
kHz indicates that 17 subject-related modes, 13 directional-related modes and 7 frequency-relatedmodes account for more than 95% variations of HRTFs at corresponding dimension respectively. The
reconstructed individualized HRTF magnitudes from 17 subject-related modes exhibit a reasonableaccuracy but some errors occur at high-frequency above 10 kHz and source direction contralateral toear. Based on the results, it is also deduced that 17 anthropometric parameters at least are needed forthe anthropometry-based HRTF customization, and the number of anthropometric parameters used for
HRTF customization in some previous was insufficient. The future works include HRTF customization,
simplifying HRTF data and signal processing in virtual auditory display using tensor decomposition.

ACKNOWLEDGEMENTS
This work was supported by the National Natural Science Foundation of China (11674105) andState Key Lab of Subtropical Building Science, South China University of Technology.

REFERENCES
1. Blauert J. Spatial hearing: the psychophysics of human sound localization. MIT press; 1997.
2. Xie B. Head-related transfer function and virtual auditory display. J. Ross Publishing; 2013 Jul 23.
3. Wenzel EM, Arruda M, Kistler DJ, Wightman FL. Localization using nonindividualized headȋrelatedtransfer functions. The Journal of the Acoustical Society of America. 1993 Jul;94(1):111-23.
4. Jin C, Leong P, Leung J, Corderoy A, Carlile S. Enabling individualized virtual auditory space usingmorphological measurements. In Proceedings of the First IEEE Pacific-Rim Conference on Multimedia(2000 International Symposium on Multimedia Information Processing) 2000 Dec (pp. 235-238).
5. Xie BS. Recovery of individual head-related transfer functions from a small set of measurements. The
Journal of the Acoustical Society of America. 2012 Jul;132(1):282-94.
6. Kistler DJ, Wightman FL. A model of head-related transfer functions based on principal components
analysis and minimum-phase reconstruction. The Journal of the Acoustical Society of America. 1992
Mar;91(3):1637-47.
7. Larcher V, Warusfel O, Jot JM, Guyard J. Study and comparison of efficient methods for 3-d audio
spatialization based on linear decomposition of HRTF data. In Audio Engineering Society Convention
108 2000 Feb 1. Audio Engineering Society.
8. Cichocki A, Mandic D, De Lathauwer L, Zhou G, Zhao Q, Caiafa C, Phan HA. Tensor decompositions
for signal processing applications: From two-way to multiway component analysis. IEEE Signal
Processing Magazine. 2015 Mar;32(2):145-63.
9. Grindlay G, Vasilescu MA. A multilinear (tensor) framework for HRTF analysis and synthesis. In2007
IEEE International Conference on Acoustics, Speech and Signal Processing-ICASSP'07 2007 Apr 15(Vol. 1, pp. I-161). IEEE.
10. Huang Q, Li L. Modeling individual HRTF tensor using high-order partial least squares. EURASIP
Journal on Advances in Signal Processing. 2014 Dec 1;2014(1):58.

5218

11. Xie B, Zhong X, Rao D, Liang Z. Head-related transfer function database and its analyses. Science in
China Series G: Physics, Mechanics and Astronomy. 2007 Jun 1;50(3):267-80.
12. Møller H. Fundamentals of binaural technology. Applied acoustics. 1992 Jan 1;36(3-4):171-218.
13. Hofman P, Van Opstal A. Binaural weighting of pinna cues in human sound localization. Experimental
brain research. 2003 Feb 1;148(4):458-70.
14. Xie B, Zhong X, He N. Typical data and cluster analysis on head-related transfer functions from Chinese
subjects. Applied Acoustics. 2015 Jul 1;94:1-3.

5219