Spatial Principal Component Analysis of Head-Related Transfer Functions
using their complex logarithm with unwrapping of phase
Shouichi Takane (1)(1) Department of Information and Computer Science, Faculty of Systems Science and Technology, Akita Prefectural University,

84-4 Ebinokuchi, Tsuchiya, Yurihonjo, Akita, 015-0055 Japan, E-mail: takane@akita-pu.ac.jp

Abstract
Head-Related Transfer Function (HRTF) changes depending on sound source position and has strong individuality. In order to compactly represent the spatial change of the HRTF in acceptable accuracy, the Spatial Principal
Component Analysis (SPCA) has been adopted in many researches. Accuracy of the reconstructed HRTFs with
the SPCA is different depending on which domain is used for the SPCA. In this article, a new procedure for the
SPCA of the HRTFs using complex logarithm of the HRTFs with unwrapping of phase was proposed. Moreover,the performance of the procedure was evaluated via the SPCA of a single HATS. The results showed that the
proposed SPCA procedure has the potential to reconstruct the HRTFs with the smaller distortion in frequency
domain using the smaller number of principal component than those in the conventional ones.
Keywords: Spatial Principal Component Analysis, Head-Related Transfer Function, Complex logarithm,
Unwrapping of phase

1 INTRODUCTION
Head-Related Transfer Function (HRTF) is deﬁned as an acoustic transfer function from sound acquired at a
center point when a listener is absent to that acquired at the listener’s ear[1]. The HRTF varies due to thesound source position and has strong individuality in both objective and subjective senses. Some systems havebeen proposed in order to synthesize the sound at the listener’s ears based on the synthesis of binaural soundsignals involving his/her HRTFs [2, 3, 4, 5, 6]. In order to achieve the accurate synthesis in these systems,
the HRTFs are ideally acquired in all directions around the listener and optimized for him/her. While a study
considering the effecient sampling scheme of the HRTF measurement exists[7], data size of such set of HRTFs
may become numerous. Therefore a requirement arises that data size of his/her HRTFs acquired in all directions
must be as compact as possible with their synthesis accuracy accomplished to some extent.
One of effective methods for the compact representation of the HRTFs is the Principal Component Analysis(PCA)[8, 9], or the spatial feature extraction method[10, 11, 12]. All of them have their theoretical basis on the
PCA or Singular Value Decomposition (SVD) of matrices assembled by using the HRTFs or their corresponding
impulse responses called the Head-Related Impulse Responses (HRIRs). In such kind of methods, the spatial
variation of HRTFs is modeled by using small number of principal components or eigenvectors. In this paper,
this method is called the Spatial PCA (SPCA) of the HRTFs/HRIRs, as Xie did[11].
Although the researches on the compact representation of HRTFs based on the SPCA exist [8, 9, 10, 11, 12],
there are some differences in which domain was used for the SPCA among these studies. Kistler et al. used thelogarithm of the HRTF amplitude[9], and Xie used the linear amplitude of the HRTF[11]. In these two studies,
the HRTFs were assumed to have only minimum-phase component and the phase component were reconstructedby using the Hilbert transform[13]. Chen et al. used the complex-valued frequency spectrum of the HRTFs[10],and Wu et al. used the HRIRs[12]. No minimum-phase approximation was required in these studies. The
SPCA can be successively and commonly adopted with the use of each parameter, and the modelling accuracyof the HRTFs was regarded sufﬁcient in both objective and subjective viewpoints in all of these studies, withlinear combination of relatively small number of principal component. However, it is known that difference inthe modelling accuracy by choosing each of the four domains. Liang et al. compared the modelling accuracy3048

between the linear and the logarithmic magnitude of the HRTFs in the horizontal plane, assuming that theminimum phase approximation was satisfactory[14]. Takane compared the reconstruction accuracy of the HRTFs
among the results of the SPCA with each of the four domains[15]. When the accuracy was compared in thesame number of principal components, both researches reached the same conclusion that the SPCA using thelinear amplitudes of the HRTFs has the relatively better accuracy than that using the logarithmic ampitude under
the assumption that the minimum phase approximation is accepted. Takane also concluded, in the case that theminimum phase approximation is not accepted, that the SPCA using the complex frequency spectrum of the
HRTFs brings about better accuracy. In other words, the more compact representation of the spatial variation of
the HRTFs is possible with the SPCA using the complex frequency spectrum of the HRTFs.
On the other side, Takane also mentioned that the accuracy in frequency domain, observed from the Spectral
Distortion (SD) between the reconstructed and the original HRTFs, is exceptionally small when the logarithmic
amplitude was chosen as the domain of the SPCA[15]. This property has a potential that the HRTFs are reconstructed by using the relatively small number of principal components when only the accuracy in frequency
domain is focused. Since cues in the frequency domain are well known to be important[17, 18], the SPCA
using the logarithmic domain may have the property desirable for the sound localization. Based on such consideration, yet another SPCA using the complex logarithm of the HRTFs is proposed in this article. Complex
logarithm of the HRTFs is used for assembling the covariance matrix, and obtained principal component arealso complex. It should be noted that the imaginary component of the complex logarithm of the HRTFs must
be unwrapped. Effectiveness of the proposed method is demonstrated via the modeling of a set of HRTFs
covering all directions. The proposed procedure does not require the minimum-phase approximation, and themost attractive feature of the proposed procedure is its lower (and better) distortion characteristics in frequency
domain. This may come from the usage of the logarithm of the HRTFs. Kistler et al. also used the logarithmof the HRTFs under the minimum-phase assumption[9], but they did not report such feature since they did notcompare the results of the SPCA using various domains.

2 MODELING OF THE HRTFs/HRIRs BASED ON SPCA
2.1 Outline of the SPCA
The following is general procedures for the SPCA of the HRTFs/HRIRs:
1. Average of a certain set of M vectors gm (m = 1,· · · , M) is calculated as follows:
1 g .
M m=1 m
M

gav =

(1)

2. Covariance matrix R is obtained by calculating the following equation:
 

1 
g − gav · gm − gav H .
M m=1 mM

R=

(2)

Note that H indicates the Hermitian transpose. Size of the matrix R is N × N, where N is denoted as thelength of the vector gm .
3. The matrix R is decomposed into N pairs of principal component (eigenvectors) and eigenvalues by solving the following eigenvalue problem:
R · qk = λ k · qk .
(3)
As a result, a set of the eigenvalues and principal components, λ k and qk (k = 1,· · · , N ), is obtained. Notethat λ k (k = 1,· · · , N ) are sorted from the largest to the smallest, i. e., λ 1 ≥ λ 2 ≥ λ 3 ≥ · · · ≥ λ N , and theprincipal components, qk , are also arranged to the corresponding eigenvalues.

3049

4. By using the matrix, Q, with qk in its column vector, the weighting vector, wm , corresponding to them-th vector gm is calculated as follows:
wm = QH gm − gav .

(4)

As a result of the PCA, the weight wm is approximated by using q1 ∼ q K (1 ≤ K ≤ N ) as follows:(wm ) K = QHK gm − gav ,(5)where Q K is a matrix with its column vectors q1 ∼ q K . The length of the vector (wm ) K becomes K.
The vectors and matrices are assumed to have complex values in the procedures stated above. If their values
are real, the Hermitian transpose is changed to the simple tranpose T . The value m reﬂects the sound source
position, and also the subject if the HRTFs/HRIRs of multiple subjects are used for assembling the covariancematrix.
The parameters are reconstructed by using the principal components as follows:
gm K = Q K · (wm ) K + gav .

(6)

The computed parameter gm K becomes approximation when K < N, but this may have acceptable accuracy
in principle when the Cummulative Proportion of Variance (CPV) R2 (K ) is close to 1.0. The CPV is deﬁned
by using the eigenvalues of the covariance matrix, λ k (k = 1,· · · , N ), as follows:K
k=1
R (K ) =
N


λk
,

2

(7)

λk

k=1

where N is the total number of component, equals to the length of the vector gm (m = 1,· · · , M).
2.2 Domains used for the SPCA
As mentioned in the Introduction, some kinds of parameters are subject to the assembly of the covariance matrix. Kistler et al. used the logarithm of the HRTF amplitude[9], and Xie used the amplitude of the HRTF[11].
In these studies, the minimum-phase properties of HRTFs are assumed. Chen et al. used the complex HRTF
spectrum[10]. While these studies extracted the parameters from the HRTFs, Wu et al. used the HRIRs[12]. Yetanother parameter, complex logarithm of the HRTFs, is newly listed in this paper. According to these studies
and the new proposal, such differences are treated as the modelling domains in this paper. It is called domain“I” when the HRIR is used for the PCA, domain “C” when the complex HRTF spectrum is used, domain“F” when the amplitude spectrum of the HRTF is used, domain “L” when the real logarithm of the amplitude
spectrum of the HRTF is used, domain “CL” when the complex logarithm of the HRTF is used together withunwrapping, and domain “CLN” when the complex logarithm of the HRTF is used without unwrapping. Thesedeﬁnitions are summarized in Table 1.
Considering the SPCA in each domain, it is obvious in the domains C, F, L, CL and CLN that the frequency
spectrum has the following symmetric relations for k = 1,· · · , N, where N indicates the length of a HRTF:∗
Hm (k) = Hm
(N − k),(8)∗log Hm (k) = log Hm(N − k),(9)

Am (k) = Am (N − k),3050

L m (k) = L m (N − k).

(10)Table 1. Correspondence between modelling conditions and symbols.
DomainHRIR
HRTF
Amplitude of HRTF
Real log-amplitude of HRTF
Complex log-amplitude of HRTF (w/ unwrapping)Complex log-amplitude of HRTF (w/o unwrapping)Symbol
I
C
F
L
CL
CLN

In Eqs. (8) to (10), Hm (k) is the k-th component of Hm , Am (k) = |Hm (k)|, L m (k) = log Am (k), and ∗ indicatesthe conjugation. This means that the vector length is almost halved in these domains.
When the covariance matrices assembled in the domains I, C, F, L, CL and CLN are respectively denoted as
R(I) , R(C) , R(F) , R(C) , R(CL) , and R(CLN) , size of R(I) is N × N, while that of R (C) , R (F) , R (CL) and R (CLN)
is (N/2 + 1) × (N/2 + 1). The domains C, F, L, CL and CLN have the advantage in the smaller size of thecovariance matrix, but the data in the frequency domain become complex in the domains C, CL and CLN.
Therefore data size per one component of the covariance matrix becomes twice in these domains. Using theamount of data of R(I) as reference, that of R(C) , R(CL) and R(CLN) is about 1/2, and that of R(F) and R(L) isaround 1/4.
2.3 The SPCA using complex logarithm of HRTFs
The m-th HRIR and HRTF are respectively expressed as hm and Hm , and Hm is further decomposed into its
amplitude and phase components as follows:
Hm = Am exp { jΘm } ,(11)where Am and Θm indicates vectors stored the amplitude and phase component of the m-th HRTF, respectively.
The domains I, C, F, L, CL and CLN mean that hm , Hm , Am , Lm ≡ log Am (m = 1,· · · , M), log Hm and Log Hm
are respectively used for the SPCA. Obviously the following equation stands in the domain CLN:=

Log Hmlog Am + jΘm ,(12)where LogH denotes the complex logarithm and its imaginary part of is in the range [−π,π). When the above
equation is calculated in computer, the imaginary part of this equation is automatically truncated in the range[−π,π), leading to its discontinuities along the frequency. This must be affected to the reconstruction of HRTF
by using the SPCA. Therefore, this imaginary part is unwrapped in the domain CL as follows:log Hm=log Am + jW [Θm ] ,(13)where W [·] is unwrapping function.
2.4 Calculation of HRTFs from the SPCA in various domains
When the m-th vector is reconstructed by using
domains  the ﬁrst
K principal (F)  components
 (L)   in theI, C, F, L, CL(C)(CL)
and CLN, they are respectively denoted as h(I)and LogH(CLN).
m K , H m K , A m K , L m K , logH mm
KKIn order to obtain the corresponding reconstructed HRIR and HRTF from the SPCA results in each domain, thelength of the vectors are almost doubled with the relations in Eqs. (8)∼(10) except for the domain I. Then the
FFT or IFFT is used with the additional
complexexponentials
 computation of the real and for the domains
 the

(CL)
(CL)
L, CL and CLN. For example, h(CL)isobtainedby:calculating
H
from
logH
, then the length
m
m
m
K
K
K
 (CL) 
almost doubled, and calculating the IFFT of the resultant vector.
of Hm
K

3051

3 RELATION BETWEEN ACCURACY AND NUMBER OF PRINCIPAL COMPONENTS
A database of HRIRs of KEMAR HATS (Head And Torso Simulator) presented by Media lab. of MIT[19] wasused. This database involves 710 pairs of HRIRs (total: 1420) with sampling frequency of 44.1 kHz. This value
corresponds to M in Eqs. (1) and (2).
3.1 Conditions of analysis
The initial delay in each response was extracted, then 256 sample points were taken as the data for the analysis, windowing with latter half of 512-points Blackman-Harris window function adjusting its peak at that of
the HRIR. The SPCA was executed by constructing the covariance matrices from HRIRs (called as domainI), HRTFs (domain C), amplitude of HRTFs (domain F), log-amplitude of HRTFs (domain L), complex logamplitude of HRTFs with unwrapping (domain CL) and that without unwrapping (domain CLN). The HRTFs/HRIRs in all directions (710 directions×2 ears) were used, and the average vector (Eq. (1)) and the covariancematrix (Eq. (2)) were calculated in each domain.

Cumulative proportion of variance

3.2 Cumulative Proportion of Variance (CPV)
When the PCA is utilized for some data, the cumulative proportion of variance R2 (K ) deﬁned as Eq. (7)1
is used for the reference indicating how much variancein data is covered by using the ﬁrst K principal com0.9
ponents. Change of the cumulative proportion of vari0.8
ance in each domain is plotted in Fig. 1. It is foundout from this ﬁgure that the CPV is monotonically in0.7
creased and converges to 1.0 as the number of component is increased in all domain. Among these six
0.6
domains, the cases CLN and CL respectively have theslowest and the fastest increase as the number of prin0.5
cipal components increases. This indicates that the the0.4
smallest number of principal components can cover acertain proportion of variance in the domain CL and
0.3
the largest number is required for the domain CLN.
0
5
10
15
20
25
30
35
40
Number of component(s)Since the difference between the domain CL and CLN
is whether the unwrapping of phase was executed or
Fig. 1. Change in cumulative proportion of variancenot, the difference in Fig. 1 certainly comes from the(CPV) with number of component in PCA. Lines andunwrapping processing.
marks for the cases C, L, CL and CLN are relatively
The previous studies had various reference values in thicker and larger than those for the cases I and F.
the proportion of variance, more than which the corresponding principal components are omitted. Kistler et
al. set this value to 0.90[9], Chen et al. and Wu et al. set this value to 0.999[10, 12], and Xie set this to around
0.98[11]. Since the amount of data and analyzing purposes are different among these studies, direct comparison
among these values is impossible. However, it can be said that the value more than 0.9 is set in these studies.
After these studies, the least numbers of component to cover four values of the CPV, 0.90, 0.95, 0.99 and 0.999
are indicated in Table 2. Seeing Table 2, the domain CLN has the numbers of component largest in all of theCPV values. Since the dimension of the covariance matrix for this domain is 129 × 129, around 81% (105/129)
of principal component are required to cover 0.999 of the variance in the domain CLN. On the contrary, thedomain CL has the smallest values among six domains when the CPV value is 0.90 or 0.95, and the domain C
has the smallest value when the CPV value is 0.99 or 0.999. From these comparison, it can be stated that thedomains C and CL have almost the same characteristics. However, the domain C requires number of component
less than that for the domain CL when the CPV value is raised to 0.999.

3052

I
C
F
L
CL
CLN

Table 2. The least number of component to cover the cumulative proportion of variance in each case.
Case
I
C
F
L
CL
CLN

0.90
8
4
5
6
2
29

Variance
0.95 0.99
10
20
6
11
7
14
11
32
4
12
47
84

0.999
39
20
31
78
40
105

3.3 Comparison of reconstruction accuracy
The CPV of the covariance matrix is an effective criterion for the coverage of variance with a certain numberof component derived from the PCA. However, Takane indicated that the direct comparison of the CPV against
a certain number of principal component among each domain is not possible[15]. Therefore, the reconstruction
accuracy obtained for each domain is compared in this section. For this purpose, Spectral Distortion (SD) and
Signal-to-Distortion Ratio (SDR) is used.
The SD is deﬁned as standard deviation in log-amplitude of two frequency spectra as follows:Nf A(k) 1  [dB],SD [A, A]
20 log10 (14)^ =N f k=1
 A(k) 
where A and A are the frequency amplitude spectrum of the original and the estimated responses, respectively,
and A(k) and A(k) are the k-th component of the vectors A and A, respectively. Here the value of N f is set
so that it gets closest to the frequency of 20 kHz. The closer A is to A, the smaller the SD is.
Signal-to-Distortion Ratio (SDR) is deﬁned as the level difference between the energy of the original impulse
response and that of the error as follows:
 
||h||
SDR h, h = 10 log10 h − h

[dB],

(15)

where ||·|| indicates the Euclid norm of the vector. The closer h is to h, the larger the SDR is.
As for the SDRs in the domains F and L, the original impulse responses are ones constructed with its minimum phase component, which are different from the ones in the other domains. The SDRs in each case werecomputed as how much the reconstructed impulse response differs from the desired one. It is noted that thispoint is related to the calculation of SDRs, and not to that of SDs, since the SD is computed by using only
magnitude of the reconstructed and the original HRTFs.
Accuracy in frequency domain with number of component set to K is computed in a certain domain X(X=I,C,F,L,CL,CLN), the overall average of SD and SDR weres calculated by using the following equations:M  2
1  AvSD(X,K ) =SD Am , A(X)
,
(16)
m K
M m=1
AvSDR(X,K )

=

M
 
⎧
⎫
⎪ 1  SDR hm , h(X)⎪m K /10 ⎬10
10 log10 ⎨
⎪M
⎪.
⎩ m=1
⎭

(17)

Changes of the average SD and SDR with number of component(s) K in each domain are plotted in Fig. 2 andFig. 3, respectively.

3053

From Fig. 2, it is found out that the least average SD is obtained for the domain L when the number ofcomponent is less than 12. When the number of component is more than 12, the domains L and CL havealmost the same average SD, and they acquire the lowest values of average SD among all domains. They
correspond to the domains that the logarithm of the HRTFs is used for assembling the covariance matrix. The
domain C has the SD smaller and almost equal to the domains C and CL when the number of component isincreased. The average SD remains relatively large for the domains I and F.
From Fig. 3, the highest SDR is obtained for the domain C over all numbers of principal component. The
characteristics of that for the domain CL is similar to that of the domain L, and the domain CLN has the lowest
average SDR over all numbers of principal components. As a result, the proposed SPCA method does not havethe best characteristics, compared with the other domains. One possible reason for this incorrespondence in thereconstruction accuracies between frequency and time domains is that the local difference in frequency domain
may bring about a relatively large difference in time domain. However, it is found that the acceptable accuracy
is achieved for the larger number of component in all domains except the domain CLN. At least the proposed
method is appropriate in the sense that the reconstruction accuracy for the domain CL is in the similar level tothat for all domains.
40

10

I
C
F
L
CL
CLN

30

Average SDR [dB]

Average SD [dB]

8

35

6

4

25
20
15
10

2

0

I
C
F
L
CL
CLN

5
0
0

5

10

15

20

25

30

35

40

Number of component(s)Fig. 2. Change of average SD with number of component(s) in six domains, respectively expressed as
AvSD(I,K ), AvSD(C,K ), AvSD(F,K ), AvSD(L,K ),AvSD(CL,K ) and AvSD(CLN,K ) (K = 1,· · · ,40).

0

5

10

15

20

25

30

35

40

Number of component(s)Fig. 3. Change of average SDR with number ofcomponent(s) in six cases, respectively expressedas AvSDR(I,K ), AvSDR(C,K ), AvSDR(F,K ),AvSDR(L,K ), AvSDR(CL,K ) and AvSDR(CLN,K )(K = 1,· · · ,40).

4 CONCLUDING REMARKSIn this paper, the SPCA of the HRTFs using complex logarithm of the HRTFs with unwrapping of phase wasproposed. Although there exist the articles in which the logarithm of the HRTFs is applied to the SPCA, theirphase component are assumed to have minimum-phase. On the contrary, the proposed procedure requires nosuch assumption. Moreover, the performance of this method was examined via the SPCA of a single HATS.
Summarizing the above results, the proposed SPCA procedure has the potential to reconstruct the HRTFs withthe smaller SD using the smaller number of principal components than those in the other domain.
The reconstruction accuracy was evaluated by mainly using the SD in this paper, but it is also necessary to
conﬁrm the above mentioned feature via subjective experiment. This is regarded as one of the important future
work concerning the proposed procedure in this article.

3054

ACKNOWLEDGMENTS
A part of this study was supported by JSPS KAKENHI (JP17K00242).

REFERENCES[1] J. Blauert, Spatial Hearing (MIT Press, 1983).
[2] M. Otani and T. Hirahara, “A dynamic auditory display: its design, performance, and problems in HRTF
switching,” Proc. Japan-China Joint Conf. Acoust., SS-1-3, 2007.
[3] S. Yairi, Y. Iwaya and Y. Suzuki, “Estimation of detection threshold of system latency of virtual auditory
display,” Applied Acoustics, 68(8), 851-863, 2007.
[4] J. D. Miller, “Slab: a software-based real-time virtual acoustic environment rendering system,” Proc. 2001
ICAD, 279-280, 2001.
[5] S. Takane, S. Takahashi, Y. Suzuki, T. Miyajima, “Elementary real-time implementation of a virtual acoustic display based on ADVISE,” Acoust. Sci. & Tech., 24(5), 304-310, 2003.
[6] D. R. Begault, 3-D sound for virtual reality and multimedia (AP Professional, 1994).
[7] W. Zhang, M. Zhang, R. A. Kennedy and T. D. Abhayapala, “On high-resolution Head-Related Transfer
Function Measurements: An efﬁcient sampling scheme,” IEEE Trans. Audio, Speech & Lang. Process.,20(2), 575-584, 2012.
[8] W. L. Martens, “Principal component analysis and resynthesis of spectral cues to perceived direction,”
Proc. ICMC, 274-281, 1987.
[9] D. J. Kistler and F. L. Wightmann, “A model of head-related transfer functions based on principal component analysis and minimum-phase reconstruction,” J. Acoust. Soc. Am., 91(3), 1637-1647, 1992.
[10] J. Chen, B. D. Van Veen and K. E. Hecox, “A spatial feature extraction and regularization model for thehead-related transfer functions,” J. Acoust. Soc. Am., 97(1), 439-452, 1995.
[11] Bo-Sun Xie, “Recovery of individual head-related transfer functions from a small set of measurements,”
J. Acoust. Soc. Am., 132(1), 282-294, 2012.
[12] Z. Wu, F. H. Y. Chan, F. K. Lam, J. C. K. Chan, “A time domain binaural model based on spatial feature
extraction for the head-related transfer functions,” J. Acoust. Soc. Am., 102(4), 2211-2218, 1995.
[13] A. V. Oppenheim and R. W. Schafer, Discrete-time signal processing, third edition (Prentice Hall, 2010).
[14] Z. Liang, B. Xie and X. Zhong, “Comparison of principal component analysis of linear and logarithmic
magnitude of Head-Related Transfer Functions,” Proc. 2nd IEEE International Congress on Image andSignal Processing, 2009.
[15] S. Takane, “Effect of domain selection for compact representation of spatial variation of Head-Related
Transfer Function in all directions based on Spatial Principal Components Analysis,” Applied Acoustics,
101(1), 64-77, 2016.
[16] T. Nishino, S. Kajita, K. Takeda and F. Itakura, “Interpolating head related transfer functions in medianplane,” Proc. 1999 IEEE Workshop on Applications of Signal Processing to Audio and Acoustics, New
Paltz, NY, 17-20, 1999.
[17] F. Asano, Y. Suzuki and T. Sone, “Role of spectral cues in median plane localization,” J. Acoust. Soc. Am.,
88(1), 159-168, 1990.
[18] K. Iida, M. Itoh, A. Itagaki and M. Morimoto, “Median plane localization using a parametric model of thehead-related transfer function based on spectral cues,” Applied Acoustics, 68, 835-850, 2007.
[19] W. G. Gardner and K. D. Martin, “HRTF measurements of a KEMAR,” J. Acoust. Soc. Am., 97, 39073908, 1995.

3055