Get 20M+ Full-Text Papers For Less Than $1.50/day. Start a 14-Day Trial for You or Your Team.

Learn More →

Modeling of Individual HRTFs based on Spatial Principal Component Analysis

Modeling of Individual HRTFs based on Spatial Principal Component Analysis IEEE/ACM TRANSACTIONS ON AUDIO, SPEECH AND LANGUAGE PROCESSING, VOL. 28, NO. 1, 2020 1 Modeling of Individual HRTFs based on Spatial Principal Component Analysis Mengfan Zhang, Zhongshu Ge, Tiejun Liu, Xihong Wu and Tianshu Qu Abstract—Head-related transfer function (HRTF) plays an multimedia, and virtual reality. Measuring the high spatial important role in the construction of 3D auditory display. This resolution HRTFs for each potential user is difficult, so schol- paper presents an individual HRTF modeling method using deep ars basically use non-individual HRTFs. However, using non- neural networks based on spatial principal component analysis. individual HRTFs may lead to some perception errors such as The HRTFs are represented by a small set of spatial principal in-head localization, front-back confusion, and a breakdown components combined with frequency and individual-dependent weights. By estimating the spatial principal components using of elevation discrimination ability [2], [3]. Thus, attaining deep neural networks and mapping the corresponding weights to individual HRTFs is very important and urgent in virtual a quantity of anthropometric parameters, we predict individual auditory scene synthesis. HRTFs in arbitrary spatial directions. The objective and subjec- At present, experimental measuring is an accurate method tive experiments evaluate the HRTFs generated by the proposed to obtain individual HRTFs. In the past two decades, a number method, the principal component analysis (PCA) method, and the generic method. The results show that the HRTFs generated of research groups have performed HRTF measurements and by the proposed method and PCA method perform better than established HRTF databases [4], [5]. However, experimental the generic method. For most frequencies the spectral distortion measuring of individual HRTFs requires rigorous experimental of the proposed method is significantly smaller than the PCA conditions and complicated equipment and keeps the subjects method in the high frequencies but significantly larger in the low not moving during the measuring procedure, which make this frequencies. The evaluation of the localization model shows the PCA method is better than the proposed method. The subjective method hard to implement. localization experiments show that the PCA and the proposed With the development of computer technology, numerical methods have similar performances in most conditions. Both calculation can be used to obtain HRTFs. Common numerical the objective and subjective experiments show that the proposed calculation methods include the boundary element method method can predict HRTFs in arbitrary spatial directions. (BEM) [6], the finite element method (FEM) [7] and the finite Index Terms—anthropometric parameters, HRTF, individual, difference method (FDM) [8]. However, numerical calculation SPCA. methods are computationally expensive and depend on the availability of precise 3D geometry. For example, magnetic resonance imaging (MRI) is used to obtain individual mor- I. I NTRODUCTION phology. This method requires a non-trivial acquisition process HE head-related transfer function (HRTF) describes the and complicated calculation. acoustic transmission of sound waves from a sound In recent years, many researchers have concentrated on source to a listener’s binaural ears. It assumes the head modeling individual HRTFs. Brown et al. [9] separated the exists in a free field. In time domain, HRTF is called head- effects of different physiological structures on HRTFs, and related impulse response (HRIR) [1]. HRTF has been widely each effect was modeled with a low-order sub-filter. The used in virtual sound technology, room acoustics simulation, combination of all sub-filters represents an HRTF. Middle- brooks [10] used the frequency scaling method, assuming that Copyright 2020 IEEE. Personal use of this material is permitted. Permission from IEEE must be obtained for all other uses, in any current or future media, the HRTF spectral characteristics of diverse individuals are including reprinting/republishing this material for advertising or promotional similar, but the corresponding frequencies of spectral charac- purposes, creating new collective works, for resale or redistribution to servers teristics are different. Through the frequency scaling method, or lists, or reuse of any copyrighted component of this work in other works. Manuscript received June 6, 2019; revised September 17, 2019 and De- a new subject’s HRTF can be obtained. Zotkin et al. [11] cember 15, 2019; accepted January 11, 2020. Date of publication January 17, selected the HRTF data of the subject whose anthropometric 2020; date of current version February 4, 2020. This work was supported by parameters were closest to the new subject. Jin et al. [12] the National Natural Science Foundation of China under Grant 11590773, Grant 61175043, and Grant 61421062, and in part by High-performance applied the principal component analysis (PCA) to the HRTF Computing Platform of Peking University. The associate editor coordinating amplitude spectrum and anthropometric parameters separately the review of this manuscript and approving it for publication was Prof. A. and then constructed a linear mapping from the PCA weights Sarti. (Corresponding author: Tianshu Qu.) M. Zhang, Z. Ge, X. Wu, and T. Qu are with the Key Labora- of the anthropometric parameters to the PCA weights of tory on Machine Perception (Ministry of Education), Speech and Hear- HRTFs. Hu et al. [13] used back-propagation artificial neural ing Research Center, Peking University, Beijing 100871, China (e-mail: networks to map the PCA weights of HRTFs to the selected zhangmengfan@pku.edu.cn; gezhongshu@pku.edu.cn; wxh@cis.pku.edu.cn; qutianshu@pku.edu.cn). anthropometric parameters. However, this approach required T. Liu is with the State key Laboratory of Robotics, Shenyang Institute separately training neural networks for each spatial direction. of Automation, Chinese Academy of Sciences, Shenyang, Liaoning 110004, To predict high spatial resolution HRTFs, this method had to China (e-mail: tjliu@sia.cn). Digital Object Identifier 10.1109/TASLP.2020.2967539 measure a large HRTF database and train thousands of neural arXiv:1910.09484v6 [eess.AS] 6 Feb 2020 IEEE/ACM TRANSACTIONS ON AUDIO, SPEECH AND LANGUAGE PROCESSING, VOL. 28, NO. 1, 2020 2 network models. XXX As the subjective perception of spatialization is the ultimate (f) = HRTF (; '; f; s); (2) log goal, the individual HRTFs can also be obtained based on S  D s ' the listeners feedback [14]. Fink et al. [15] let subjects tune After removing the mean logarithmic HRTFs, we obtain the the PCA weights from average HRTFs to obtain individual log-magnitude function, which is called HRTF . HRTFs. This tuning procedure can reduce localization errors; log however, obtaining a customized HRTF for a subject is very time-consuming. The subjects need a lot of time to finish the HRTF (; '; f; s) = HRTF (; '; f; s) (f): (3) log log tuning part. Luo et al. [16] also used the tuning method to obtain individual HRTFs and first introduced deep learning C. SPCA autoencoders to HRTF. The autoecoder was used to perform feature reduction and obtained a better result than PCA. PCA is mathematically defined as an orthogonal linear The aim of this paper is to realize the modeling of individual transformation that transforms the data to a new coordinate HRTFs and to predict the HRTFs in arbitrary spatial directions. system, such that the greatest variance by some projection Spatial principal component analysis (SPCA) [17] or spherical of the data comes to lie on the first coordinate (called the harmonics (SHs) basis functions [18] can be used to spatially first principal component), the second greatest variance on the decompose HRTFs. In this study, we use SPCA to decom- second coordinate, and so forth [20]. pose HRTF into a weighted combination of spatial principal The traditional PCA method is generally used in the time components (SPCs). Through the deep neural network (DNN) or frequency domain of HRTFs [21], [22]. In contrast to tra- training, the SPCs in arbitrary spatial directions are estimated. ditional PCA models, SPCA is applied to the spatial domain. A small quantity of anthropometric parameters were selected The high spatial resolution HRTFs can be represented as the and mapped to the SPCA weights using neural networks. Then, weighted combination of orthonormal SPCs [17]. The SPCA we combined the predicted SPCs and the SPCA weights for applied to a HRTF is shown below. log each individual to reconstruct the HRTFs in arbitrary spatial directions. HRTF (; '; f; s) = d (f; s)W (; ') + H (; '); log q q av The rest of the paper is organized as follows. In Section II, q=1 the data preprocessing and SPCA are performed. In Section (4) III, modeling individual HRTFs based on SPCA using DNN where q is the identification of the SPC, W is the SPC is described. In Section IV, the objective experiments are that depends only on the source direction, ' is the elevation conducted and the objective error is analyzed. In Section V, the angle, and  is the azimuth angle. D is the number of subjective evaluation of the proposed approach is described. spatial directions. d (f; s) is the SPCA weight which varies In Section VI, the conclusion is presented. as function of frequency f and individual s. H is the mean av HRTF magnitude across the frequencies and subjects and log II. SPCA AND DATA P REPROCESSING can be calculated as follows: XX A. CIPIC Database H (; ') = HRTF (; '; f; s); (5) av log The HRIRs used in this paper are derived from the CIPIC N  S s f database [4], which is measured by U. C. Davis CIPIC where N and S are the total number of frequencies and Interface Laboratory. In this database, a blocked ear technique subjects respectively. is performed for 45 subjects (27 males, 16 females, 2 KE- To calculate the SPCs and SPCA weights, we combine MARs), and 1250 directions of HRTF data were measured for HRTF of all the frequencies, directions and subjects into log each subject. Sound source directions are in interaural-polar an (NS) D matrix H . Each column of H corresponds to coordinate. The database also contains up to 27 anthropometric a spatial direction, and each row of H represents the HRTF parameters for each subject. of an individual at a discrete frequency. We subtract the mean value from H to obtain H . Then we calculate the covariance B. Data Preprocessing matrix R: We first transform the raw HRIRs in CIPIC database into T R = H H ; (6) the frequency domain. Fourier transformation is applied to the where R is a D  D matrix. Its eigenvectors are extracted HRIRs to calculate the HRTFs. Then we transform the HRTFs and arranged as the eigenvalue reduced-order. Then the first into a logarithmic scale. Because a logarithmic scale is much Q eigenvectors are taken as the base vectors, i.e. the SPCs closer to our auditory perception [19]. Therefore, the base 10 which represent the values of basis functions at D discrete log-magnitude responses of HRTFs are computed. directions, and form a Q D matrix W : HRTF (; '; f; s) = 20log (jHRTF (; '; f; s)j): (1) log 10 W = [W (0); W (1); :::; W (D 1)]: (7) q q q Then the mean logarithmic HRTFs is calculated from all the HRTF . The mean function includes the direction and Each row of W corresponds to a SPC, and each column of W log subject independent spectral features shared by all HRTFs in represents the values of SPCs at a spatial direction. We name the CIPIC database. the column of W as direction vector of SPCs (DV-SPCs). IEEE/ACM TRANSACTIONS ON AUDIO, SPEECH AND LANGUAGE PROCESSING, VOL. 28, NO. 1, 2020 3 TABLE I TABLE II CUMULATIVE PERCENTAGE OF VARIANCE IN RECONSTRUCTED HRTF S ELECTED KEY ANTHROPOM ETRIC PARAMETERS MAGNITUDES IS INCREASED W ITH THE NUMBER OF SELECTED SPCS Variable Measurement Variable Measurement Cumulative percentage x head width d cavum concha height 1 1 The number of SPCs (Q) of variance (%) x * head height d cavum concha width 2 3 Left ear Right ear x head depth d fossa height 3 4 1 16.54 20.14 x shoulder width d pinna height 12 5 5 52.20 55.33 d pinna width 10 62.29 64.84 * This parameter is only used in the prediction of ITDs. 20 70.10 71.85 50 78.33 79.54 60 80.09 81.22 individuals [13], [22]. We process the HRIRs and the anthro- 80 82.93 83.98 100 85.11 86.09 pometric parameters in CIPIC database as follows. 200 91.03 91.56 First, Fourier transformation is applied to the HRIRs in the 500 97.07 97.22 CIPIC database, and the mean of the obtained HRTF data is subtracted. For each sampled direction in the CIPIC database, we apply the traditional PCA to the HRTF data of all the Since all the SPCs are orthogonal to one another, by the use subjects, of these SPCs, we can obtain the SPCA weights: d = H W ; (8) HRTF (f; s) = d (s)W (f) + H (f); (11) q q av where d is a (NS)Q matrix composed of the SPCA weights q=1 for all the individuals and frequencies, W is the transpose where W is the PC. Note that the traditional PCA needs to of W . Finally, we can predict all the HRTFs as model HRTFs in each measured spatial direction, which means we apply PCA 1250 times in total. The resultant PCs and H = dW + H ; (9) AV H depend only on the frequency, whereas the SPCs and H av av where H is an (NS)  D matrix. Each row of H AV AV obtained by SPCA depend only on spatial directions. d (s) is corresponds to the H in D directions, and all of the rows av the PCA weight which varies as function of individual s. are identical. Second, multiple linear regression analysis is used to an- If the number of the SPCs is chosen to be equal to the alyze the relationship between the PCA weights and the number of spatial directions, i.e. Q = D, then the HRTFs anthropometric parameters. can be fully represented without loss, as shown in Eq.(4). If d = s + e; (12) Q < D is selected, then only an approximate representation of q the original HRTF magnitude can be obtained. The principal where is the regression coefficient, and e is the error. component variances are the eigenvalues of the covariance ma- Then we use t statistics to identify which parameters have trix R. The cumulative percentage of variance in reconstructed a significant effect on the PCA weights. HRTF is calculated as Third, Pearson correlation coefficient analysis is introduced to measure the strength of dependence between each of the q=1 V ar =  100%; (10) D two anthropometric parameters, q=1 (s s )(s s ) i i j j where  is an eigenvalue, and V ar is the cumulative per- r = p ; (13) ij 2 2 centage of variance. V ar is increased with the selection of Q, ( (s s ) (s s ) ) i i j j as shown in Table I. We can restore more than 70% of the where s is a vector that comprises 27 anthropometric pa- total variability by selecting the first 20 SPCs. We can recover rameters, i; j = 1; 2; :::; 27 and i 6= j. r is the degree of ij more than 80% of the total variability by selecting the first 60 correlation between two anthropometric parameters s and s . i j SPCs. Prior studies reported that more than 90% of the total This procedure is used to reduce the number of parameters to variability is enough to recover the HRTF magnitudes when make the measurement more feasible. A stronger correlation applying PCA in frequency domain [13], [23]. Therefore, we indicates that one parameter could be represented by another select the first 200 SPCs to recover more than 90% of the total [23]. variability when applying PCA in spatial domain. Finally, based on the correlation between the anthropomet- ric parameters as well as the analysis of the PCA weights D. Key Anthropometric Parameters Selection and anthropometric parameters, we balance the principle of There are 27 anthropometric parameters in the CIPIC simpleness, completeness and feasibility in practice to reduce database. Measuring all the anthropometric parameters for parameters. Anthropometric parameters with large correlation each individual is a time-consuming and difficult process. are reduced while considering the theoretical and practical in- In addition, not all of these anthropometric parameters are fluence of them on HRTFs to determine which to be remained. strongly related to HRTFs. Therefore, it is necessary to select Eight anthropometric parameters that have a significant effect a small set of anthropometric parameters which are most on the PCA weights are finally selected. The quantities of strongly related to the variations in the HRTFs of different selected anthropometric parameters are the same as [13], [22]. IEEE/ACM TRANSACTIONS ON AUDIO, SPEECH AND LANGUAGE PROCESSING, VOL. 28, NO. 1, 2020 4 Fig. 2. The framework of individual HRTF modeling. separately in order to keep a confident link between the HRTF dissimilarity and the angular difference [24]. Additionally, Fig. 1. Five selected pinna parameters in our model. there exists some differences between an individual’s two ears; we therefore consider one individual as containing two observations which should be modeled separately to obtain As shown in Table II, the eight anthropometric parameters different SPCA weights. When modeling the SPCs, the H av are head width, head depth, shoulder width, cavum concha and the ITDs, which are related to spatial directions, we split height, cavum concha width, fossa height, pinna height and the data of these parameters into two parts, the front and rear pinna width, corresponding to x , x , x , d , d , d , d 1 3 12 1 3 4 5 parts. In the CIPIC database, the azimuth angles composed a and d in the CIPIC database. Note that the head height vector azi = [80;65;55;45 : 5 : 45; 55; 65; 80], and parameter x is only used in the prediction of ITDs. The first the elevation angles make up a vector elev = 45+5:625[0 : four anthropometric parameters are head and torso parameters, 49]. The front part contains elev  90 in all of the azimuth which can be measured by a caliber rule or taking pictures. angles, and the remaining angles belong to the rear part. Thus, The last five anthropometric parameters are pinna parameters, each part contains a total of 1250 sets of data, including 625 which can be obtained by photograph annotation as shown in sets of data for two ears. For the three parameters, the SPCs, Fig. 1. the H and the ITDs, we train two DNNs for each parameter av in order to obtain a desirable prediction result. III. M ODELING OF I NDIVIDUAL HRTF S To sum up, through securing a quantity of anthropomet- ric parameters, we can reconstruct the individual’s binaural A. Outline HRIRs in arbitrary spatial directions. In this paper, the magnitude spectra of HRTFs are modeled based on SPCA using neural networks. The phase of HRTFs B. Prediction of SPCA Weights for an Individual is calculated by minimum-phase reconstruction [21]. Fig. 2 depicts the framework of individual HRTF modeling. Given that the SPCA weights are a function of frequency First we model the SPCA weights, the SPCs, the H and and an individual’s morphology, we model the SPCA weights av the ITDs, respectively. The SPCA weights can be obtained by and eight key anthropometric parameters based on neural modeling the individual’s morphology, for the SPCA weights networks. After securing eight anthropometric parameters of vary as functions of anthropometric parameters and frequency. an individual, we can estimate the SPCA weights for the A quantity of key anthropometric parameters is selected in individual in all of the frequency points using neural network Section II-D to represent the human morphology, therefore, the models. SPCA weights for any individual outside the database can be A total of 37 subjects s (m = 1; 2; :::; 37) in the CIPIC estimated from a small set of anthropometric measurements. database contain all of the eight anthropometric parameters. As DNNs are used to predict the SPCs, the H and the ITDs aforementioned in Section III-A, one individual’s two ears are av in arbitrary spatial directions. Accordingly, HRTFs with high modeled separately to obtain different SPCA weights. Thus, directional resolution can be recovered by solving the Eq. we have 74 sets of anthropometric parameters; considering the (4) using the predicted SPCs, SPCA weights and H . Then total number of frequency points is N = 200, we have 14800 av the minimum-phase reconstruction method is used to generate sets of the SPCA weights. mono HRIRs [21]. Finally, binaural HRIRs are obtained using For each frequency point f (k = 1; 2; :::; N), we build estimated ITDs and the corresponding left and right mono a model for the eight anthropometric parameters and the HRIRs. corresponding SPCA weights. The specific architecture of the Due to the symmetry between front and rear HRTFs, we neural network model is shown in Fig. 3. The inputs are model the front and rear HRTFs separately. Because HRTF the eight anthropometric parameters, and the ground truth are dissimilarity increases with their angular difference, however, SPCA weights in one frequency point. So in total we build the low dissimilarity may occur for large angular difference. 101 neural network models because of the symmetry property For instance, if the two HRTFs taken from two locations are of the Fourier transformation. symmetric with respect to the interaural axis, the two HRTFs 30 subjects’ binaural data are carefully selected to guarantee will be very similar despite the large angular difference. that the data sets for each anthropometric parameter are widely Thus, it is preferred to consider the front and rear HRTFs distributed. That is to say, most of the test data are within the IEEE/ACM TRANSACTIONS ON AUDIO, SPEECH AND LANGUAGE PROCESSING, VOL. 28, NO. 1, 2020 5 subject_033 average 0 2 4 6 8 10 12 14 16 18 20 Frequency (kHz) Fig. 5. The prediction errors of SPCA weights as a function of frequency. Fig. 3. The architecture of neural network model for SPCA weights prediction. 4. As previously discussed in Section II, the first 50 SPCs can restore about 80% of the total variability, and the first 200 SPCs can reconstruct above 90% of the total variability. Therefore, when the range of the q is from 51 to 200, the total amount of variability is around 10%. Moreover, we observe that the SPCA weights are approximately equal to 0 when q is larger than 50. Thus, we only plot the curves where q is less than 50 to make comparisons between the estimated SPCA weights and the real ones. Fig. 5 shows the prediction error, which is the mean of the absolute errors across all directions, of SPCA weights for subject 033 and the average prediction error. The figure shows that the difference between the real value and the estimated value of the SPCA weights increases when frequency increases. The SPCA weights in Fig. 4. Comparison of the real value (top) and the estimated value (bottom) high frequency are more unstable than the lower ones. for the first 50 orders of SPCA weights. C. Modeling of SPCs corresponding training data ranges [13]. These data, which As discussed in Section II, W is a QD matrix composed comprise a total of 60 sets, are used in the training phase, and of the first Q SPCs. Each row of W corresponds to an SPC, 10 of them are used as the validation set to prevent over-fitting. and each column of W represents the values of SPCs at The remaining 14 sets of data are used to test the average error a spatial direction. The DV-SPCs are a function of spatial of the neural network model. The mean and variance of the directions and represented by the column of W . We model test set and validation set are normalized using the training set the first Q SPCs by predicting DV-SPCs in all sampled D statistics to have zero mean and unit variance. All the neural spatial directions, then recombining W . network models comprise a single hidden layer, a hyperbolic 2 3 tangent activation function and a hyperbolic tangent output W (0); W (1) : : : W (D 1) 1 1 1 function, and they are feedforward backpropagation neural 6 7 W (0); W (1) : : : W (D 1) 2 2 2 6 7 networks with a learning rate of 0.001. W = 6 7 : (15) . . . . . . . . 4 5 For network models of all frequencies, each network takes . . . . about 500 epochs to converge averagely. After training all the W (0); W (1) : : : W (D 1) Q Q Q neural networks, we obtain a system for predicting the SPCA The DV-SPCs are modeled with DNNs. We set the angles weights. The mean square error (MSE) is used to test the of directly ahead (' = 0 ,  = 0 ) and the directly behind prediction system: (' = 180 ,  = 0 ) as the reference direction for the front XXX ^ and rear data sets, respectively. The inputs of DNN are the e = (d (f ; s ) d (f ; s )) ; d q k m q k m Q N  S DV-SPCs in the reference direction, the target azimuth  and q m (14) the target elevation ' in degrees. We set the ground truth as where d is estimated by neural network model, and e is the the DV-SPCs in the target direction. The specific architecture q d reconstruction error of the SPCA weights. The MSE of the of the DNN for predicting the DV-SPCs is shown in Fig. 6. overall prediction, including all the SPCA weights for all the Inputing the DV-SPCs in the reference direction reduces the individuals and frequencies, is 9.25. complexity of learning process. The main task of the DNN We randomly select a subject 033 in the CIPIC database and becomes learning the difference between the DV-SPCs in the plot the comparison of the estimated value and real value of reference direction and target direction. This can effectively its left ear’s first 50 orders of SPCA weights, as shown in Fig. improve the task result. prediction error IEEE/ACM TRANSACTIONS ON AUDIO, SPEECH AND LANGUAGE PROCESSING, VOL. 28, NO. 1, 2020 6 Fig. 8. Comparison of the real value (left) and the estimated value (right) of the H . av Fig. 6. The architecture of DNN for predicting the DV-SPCs. directions can be predicted. We combine the DV-SPCs in all the D sampled directions to obtain the QD matrix W . MSE is used to calculate the reconstruction error: XX e = (W (; ) W (; )) ; (16) W q q where W is estimated by the DNN model and e is the q W reconstruction error whose value is 1:97 10 . Fig. 7 illustrates the comparison of the first, the second, the third and the fourth SPCs’ real values and estimated values. The predicted SPCs are close to the real SPCs. As discussed in Section II, the first 5 SPCs can restore more than 50% of the total variability, and our algorithm performs well. With an increment of q, the corresponding SPC gradually becomes unstable, and the difference between the predicted SPC and the real SPC widens. Fig. 7. Comparison of the real value (top row) and the estimated value (bottom row) of the first 4 SPCs (4 columns). D. Modeling of H av The modeling of the H is similar to the prediction of the av To guarantee the variability of the test data, we index the DV-SPCs. This model is also based on DNN. The input of 1250 sets of data and select the test data every four indexes. DNN is the H in the reference direction, the target azimuth av The extra directions of data are then re-indexed and we choose and the target elevation in degrees. The ground truth is the H av the validation data every five indexes. The remaining directions in the target direction. The reference directions in addition to of data are used as training data. Thus, we uniformly distribute the selection of training set, validation set and test set are the the training set, validation set and test set in the space. same as in Section III-C. The mean and variance of the test The validation set is used to control the training phase and set and validation set are normalized using the training set prevent over-fitting, and the test set is used to estimate the statistics to have zero mean and unit variance. Both of the generalization error of the modeling. The mean and variance modeled DNNs are fully connected and have three hidden of the test set and validation set were normalized using the layers. The activation function and the output function are training set statistics to have zero mean and unit variance. The set hyperbolic tangent, and the learning rate is 0.001. The DNNs were implemented by MATLAB ’s DeepLearnToolbox- architecture and the parameters are determined through the master [25]. Each DNN is fully connected and has three hidden results of many experiments. layers. Three hidden layers are chosen because it has a better It takes about 700 epochs for the DNN models to converge. performance than using only one or two hidden layers; if more After training the DNNs, we can obtain a system to predict than three hidden layers are used, the performance does not the H in arbitrary spatial directions. All the D directions in av have a greater improvement, and the neural network is also the CIPIC database are used as the target directions to predict easy to overfit. Both the activation function and the output the corresponding H . The predicted results are shown in Fig. av function are set hyperbolic tangents with the learning rate set as 0.001. These settings which were determined through The MSE is used to calculate the reconstruction error of the the results of many experiments can lead to good prediction H : av performance. XX It takes about 5000 epochs for the DNN models to converge. 2 e = (H (; ) H (; )) ; (17) H av av By training the DNN models, the DV-SPCs in arbitrary spatial IEEE/ACM TRANSACTIONS ON AUDIO, SPEECH AND LANGUAGE PROCESSING, VOL. 28, NO. 1, 2020 7 Fig. 10. Comparison of the real value (left) and the estimated value (right) of the ITD. Fig. 9. The architecture of DNN for predicting the ITDs. where T is the estimated ITD and T is the real ITD. e is the reconstruction error of the ITDs, and its value is 2:22 10 ms. where H is estimated by the DNN model and e is the av H reconstruction error whose value is 0.195. F. Recovery of Individual HRIRs To reconstruct an individual’s HRIRs, we first model the SPCA weights using the eight key anthropometric parameters. E. Modeling of ITDs Then DNN is used to model SPCs, H and ITDs respectively. av The azimuth angle and elevation angle are introduced into the As ITD not only relates to spatial directions, but also input layer of the DNN models, the DV-SPC, H and ITD in varies with different individuals, we use the anthropometric av arbitrary spatial directions can then be predicted. parameters and spatial directions to model it. The inputs of The HRTF magnitude of a new subject s in azimuth angle DNN include head width x , head height x , as well as head m 1 2 and elevation angle ' can be reconstructed by solving depth x of an individual, target azimuth and target elevation in d d the Eq. (9). The SPCA weights for the new subject can be degrees. The individual’s ITD in the target direction is taken combined into a matrix. as the ground truth, as shown in Fig. 9. x , x and x are 1 2 3 2 3 chosen by many experiments. The three parameters together d (f ; s ); d (f ; s ) : : : d (f ; s ) 1 1 m 2 1 m Q 1 m show the best results, and adding other parameters does not 6 7 d (f ; s ); d (f ; s ) : : : d (f ; s ) 1 2 m 2 2 m Q 2 m 6 7 d = ; yield an improvement. This indicates the strong correlation 6 . . . . 7 . . . . 4 5 . . . . between the ITDs and the head size [26]. d (f ; s ); d (f ; s ) : : : d (f ; s ) 1 N m 2 N m Q N m Similar to Section III-C, 625 directions of ITD data for (19) an individual are indexed and the test data are selected every the DV-SPCs in  and ' is d d four indexes. The extra directions of data are then re-indexed W = [W ( ; ' ); W ( ; ' ); :::; W ( ; ' )] ; (20) and we choose the validation data every five indexes. The 1 d d 2 d d Q d d remaining directions of data are used as training data. In the the H in  and ' is a vector of length N : av d d end, the training data of the 30 subjects, selected in Section III-B, are combined as the training set. The validation data of H = [H ( ; ' ); H ( ; ' ); :::; H ( ; ' )] ; (21) AV av d d av d d av d d the 30 subjects and the test data of the remaining 7 subjects are then the HRTF magnitude in  and ' of the new subject d d integrated as the validation set and the test set, respectively. can be restored. The mean and variance of the test set and validation set The minimum phase reconstruction method is then em- are normalized using the training set statistics to have zero ployed to the HRTF magnitudes to generate the mono HRIRs mean and unit variance. Each DNN is fully connected and has [21]. The ITDs obtained in Section III-E are used to calculate three hidden layers to yield a better result. Both the activation the binaural HRIRs. The individual’s HRIRs in arbitrary function and the output function are set hyperbolic tangents spatial directions can then be obtained. and the learning rate is 0.001. It takes about 5500 epochs for the DNN models to converge. IV. OBJECTIVE EXPERIM ENTS Fig. 10 describes the results of ITD prediction. A total A. Evaluation using spectral distortion of 1250 spatial directions for a subjects’ ITDs are plotted compared with its real value. To evaluate the effectiveness of our proposed approach, we carried out a set of objective experiments for the proposed The mean absolute error (MAE) is used to evaluate the SPCA method, the PCA method and the generic method. reconstruction error of ITDs. The generic method used the HRTFs of the CIPIC KEMAR XXX with small ears. The PCA method applied traditional PCA to e = jT (; ; s ) T (; ; s )j; (18) T m m S  D HRTFs in each sampled spatial direction, which is described IEEE/ACM TRANSACTIONS ON AUDIO, SPEECH AND LANGUAGE PROCESSING, VOL. 28, NO. 1, 2020 8 Fig. 11. The SFRSs of the real value (1st column), the generic method(2nd column), the PCA method (3rd column) and the SPCA method (4th column) 2 at frequency of 12.35 kHz of subject 163. Generic PCA SPCA 1.5 3.3 5.1 6.8 8.6 10.4 12.1 13.9 15.7 17.4 19.2 20.9 Frequency (kHz) Fig. 13. Comparison of the generic method, the PCA method and the SPCA method for the mean SD across all directions for some frequency bins and their deviations. Fig. 12. The average prediction errors of SFRSs across all frequencies and subjects in the test set of the generic method (left), the PCA method (middle) and the SPCA method (right). Fig. 14. The localization performance of the median plane using the localization model with the generic method (left), the PCA method (middle) in Section II-D. The first twelve principal components (PCs) and the SPCA method (right). for all the PCA models are selected. When we apply PCA to HRTFs, for all of the spatial directions, we can restore an average of 92.02% of the total variability for the left ear HRTF data. and 91.71% for the right ear. These percentages of the total XX variability are close to our proposed SPCA model’s selections. SD(f; s) = jH(; ; f; s) H(; ; f; s)j; (22) Then, neural network models are used to map the selected eight anthropometric parameters to the PCA weights. Thus, a where H is the magnitude response (dB) of the measured total of 1250 neural network models are trained to yield the HRTF from the CIPIC database, and H is the magnitude PCA weights in each spatial direction. This is quite arduous response (dB) of the test HRTF. in PCA method, which will bring computation issues in the Fig. 13 shows how the SD varies with frequency of the three applications of HRTF prediction. Finally, individual HRTFs method. Slight translations were performed in frequency axis can be reconstructed by combining the PCs and the PCA in order to avoid overlapping of symbols for different methods. weights. As a result, the proposed method has the lowest SD for most The Spatial Frequency Response Surfaces (SFRSs) [27], frequencies compared to the other two methods. The average a spatial-domain visualization tool for HRTFs, is used to SD of the proposed model in all the sampled directions is evaluate the three methods. Each frequency bin in the HRTF 5.54 dB, which is 1.13 dB lower than that of the generic left or right magnitude responses constructs one SFRS, where method and 0.29 dB lower than that of the PCA method. T- magnitude is plotted against azimuth and elevation. Fig. 11 test applied on the 12 frequency bins in Fig. 13 shows that shows the SFRS of the three methods compared with the real SPCA method has significant smaller SD (p < 0:001) than value of subject 163 at frequency of 12.35 kHz. Fig. 12 shows PCA method above the frequency of 6 kHz but has significant the average prediction errors of SFRS across all frequencies larger SD (p < 0:001) than PCA method blow the frequency and subjects in the test set of the three methods. The SPCA of 6 kHz, except for 3 frequency bins, 3526 Hz (p = 0:052), method has smaller and more averagely distributed prediction 10584 Hz (p = 0:157) and 21168 Hz (p = 0:073) with no errors than the other two methods. Namely, our proposed significant difference. The results indicate that our proposed SPCA method has a much more smooth spatial shape than SPCA method predicted the HRTFs in the test set well. Since PCA method, which indicates that our algorithm predict the the test sets in the modeling of DV-SPCs, H and ITDs are av HRTF in spatial domain well. different from the training spatial directions, our proposed Then a frequency dependent spectral distortion (SD) [28] model will definitely predict the HRTFs in unmeasured spatial is used as an evaluation metric between the real and the test directions. SD magnitude (dB) IEEE/ACM TRANSACTIONS ON AUDIO, SPEECH AND LANGUAGE PROCESSING, VOL. 28, NO. 1, 2020 9 TABLE III THE AVERAGE RESULTS OF THE OBJ ECTIVE EVALUATION ON MEDIAN PLANE LOCALIZATION EXPERIMENT. Confusion rate (%) Angle of PE Method Up-down Front-back error (Deg) (Deg) Generic 11.05 16.32 17.71  2.48 26.84 PCA 7.81 11.22 13.93  2.21 21.63 SPCA 12.62 14.51 15.91  2.53 24.89 (a) Interface for azimuth experiments (b) Interface for elevation experiments B. Evaluation using the localization model Fig. 15. The user interfaces for subjects to give percepted directions of For further evaluation, an auditory based localization model sounds. [29], [30] was used to evaluate the localization performance of the generic, the PCA and SPCA calculated HRIRs. The BAUMGARTNER2014 function in the AMT tool box [31] subjects’ anthropometric parameters listed in Table II. These was used in our experiment. The BAUMGARTNER2014 func- parameters were measured before the experiments. The 4 tion is a sagittal plane localization model. As a template-based variables listed in the 1st column were measured using a comparison model, it can be used to evaluate the performance ruler with an accuracy of 2 mm. The 5 variables listed in of HRTF individualization methods. The DTFs [32] were the 3rd column were measured through photographing and extracted from HRTFs [30] as the inputs of the auditory model. manually demarcating them as shown in Fig. 1. The ruler was We apply an listener-specific sensitivity threshold of 1 and set used to obtain the scale between the picture and real objects. the differential order of the spectral gradient extraction to be All parameters were measured three times then averaged to 0 to acquire a reasonable prediction error. Other settings of reduce the measurement errors. The impulse responses from the function were set as default. the headphone, used in the experiments, to the entrances We applied the localization model on the median plane to of the blocked ear canals of subjects were measured using the BK Type 4101-A binaural in-ear microphone. During the test the performance of the generic method, the PCA method experiments, headphone equalization was performed [33]. and the SPCA method. For each method, the HRIRs of all the 37 individuals used in our paper were tested. For each elevation angle, 2 runs were conducted to reduce probable A. Azimuth Localization Experiments prediction errors. Namely, for each elevation angle of the The aim of this experiment is to compare the azimuth median plane, there were 74 predicted elevation angles of the localization performance of the HRTFs generated by the target HRIRs in total. generic method, the PCA method and the SPCA method. For Fig. 14 shows the localization performance of the generic the PCA method, the phases of the HRTFs are obtained by method, the PCA method and the SPCA method of the median the minimum phase reconstruction method, and the ITDs are plane. Table III shows the statistical analysis of the experiment yielded using the modeling method of Section III-E. results. Four indices, including the up-down confusion rate, the The stimulus in this experiment was a train of eight 250-ms front-back confusion rate, the angle of errors with standard bursts of Gaussian noise (20-ms cosine-squared onset-offset deviations and the polar RMS (root-mean square) error (PE), ramps), with 300 ms of silence between the bursts. The HRIRs were used to evaluate the localization performance. For the an- of twelve azimuth angles, 0, 30, 55, 80, 125, 150, 180, 210, gle errors of the median plane, the repeated-measures ANOVA 235, 280, 305 and 330 degrees, in three elevation angles, 0, was used to verify the significance of the mean difference. 22.5 and 45 degrees, are generated by the generic method, the According to the Shapiro-Wilk test, the data of each group PCA method and the SPCA method respectively. Note that obey the normal distribution (p > 0:05). After the Mauchly’s the HRIRs of the four azimuth angles, 55, 150, 210 and 305 spherical hypothesis test, the variance covariance matrix of degrees, in the three elevation angles are not in the training set the dependent variable is not equal ( (2) = 6:279; p = and the validation set of each DNN model, i.e. Section III-C, 0:043 < 0:05). The data are corrected by Huynh-Feldt Method Section III-D and Section III-E. Then, the stimulus is filtered ( = 0:898). The mean errors of the three methods have by the HRIRs produced by the three methods to create three significant difference (F (1:80; 64:65) = 34:37; p < 0:001). kinds of sounds. Bonferroni post-hoc test shows that the SPCA and PCA A total of three azimuth localization experiments are per- methods are significantly better than the generic method (both formed, and each experiment includes three tests. The three p < 0:001). Further more, the PCA method is better than the experiments correspond to three elevation angles, 0, 22.5 and SPCA method (p < 0:001). 45 degrees, respectively. Each test is a kind of sound generated by one method, and the order of the three tests is arranged by V. SUBJ ECTIVE EXPERIM ENTS latin square design across every three subjects. Before each In this section, the azimuth localization experiments and the test, the subject is trained using the test sound of the other elevation localization experiments were conducted to evaluate eight azimuth angles, 0, 45, 90, 135, 180, 225, 270 and 315 the performance of the proposed method. In the experiments, degrees. Through listening to these sounds, the subject can three methods were used to generate the HRTFs according to build up the spatial perception for this kind of virtual sound. IEEE/ACM TRANSACTIONS ON AUDIO, SPEECH AND LANGUAGE PROCESSING, VOL. 28, NO. 1, 2020 10 TABLE IV THE AVERAGE RESULTS OF THE AZIMUTH LOCALIZATION EXPERIM ENTS. Elevation Front-back Angle of Method angle (Deg) confusion rate (%) error (Deg) Generic 28.2  9.6 18.90  4.25 0 PCA 22.4  11.1 14.84  3.63 SPCA 24.5  12.7 15.15  4.04 Generic 29.0  12.2 19.49  4.22 22.5 PCA 17.0  9.9 15.25  3.25 SPCA 25.6  13.3 14.96  3.74 180 Generic 31.8  11.9 19.79  3.30 45 PCA 25.3  10.8 15.74  3.81 SPCA 29.0  10.7 15.28  3.75 azimuth of 0 degrees is not directly ahead. The front-back confusion rate and error of angle with their standard deviations of the azimuth localization experiments are shown in Table IV. Angle of error is the angle difference between the target angle and the judgement angle. Front- 0 60 120 180 240 300 0 60 120 180 240 300 0 60 120 180 240 300 back confusion is corrected before calculating the difference. Target Angle (Deg) Target Angle (Deg) Target Angle (Deg) The repeated-measures ANOVA was used to compare the Fig. 16. Judged direction versus target direction of all subjects using the localization performance (confusion rate and error of angle) generic method (left column), the PCA method (middle column) and the SPCA of the three methods (generic, PCA and SPCA) for three method (right column) in elevations of 0 degrees (top row), 22.5 degrees elevations (0, 22.5 and 45 degrees). According to the Shapiro- (middle row) and 45 degrees (bottom row). Two oblique lines with a slope of 135 degrees correspond to the front-back confusions. Wilk test, all groups of data obey the normal distribution (P > 0:05). Through the Mauchly’s spherical hypothesis test, the variance covariance matrix of the dependent variable is After that, thirty-six binaural sounds are randomly played to dominantly equal (p > 0:05 except for the condition of error the subject by a Sennheiser HD 650 headphone through a of angle in 22.5 degrees,  (2) = 7:97; p < 0:05). The Sound Blaster sound card. The thirty-six sounds contain twelve data, which doesn’t satisfy spherical assumption, are corrected directions’ sounds, and each direction appears three times. The by Greenhouse-geisser Method ( = 0:718). For confusion subject can listen to one sound many times until he/she can rate, there is no significant difference for three methods in identify the exact perceived direction. The subject gave the elevation of 0 and 45 degrees (F (2; 34) = 2:41; p = 0:11 exact direction of each sound he/she perceived during the test and F (2; 34) = 2:10; p = 0:14), but a significant difference through an interface on a computer, which is shown in Fig. in elevation of 22.5 degrees (F (2; 34) = 8:71; p < 0:005). 15 (a). The blue dashed line in the interface is moved with Bonferroni post-hoc test for elevation of 22.5 degrees shows the mouse cursor and indicates the exact angle. After each that the PCA method has significantly smaller confusion experiment, there were five minutes for a break. rate than the SPCA and generic methods (p < 0:05 and Eighteen subjects (11 male 7 female, age from 21 to 27) p < 0:005); the confusion rates of the SPCA method and with normal hearing took part in the experiments. All of the the generic method have no significant difference (p = 0:92). experiments were performed in a sound booth (Background For error of angle, there are significant differences for all three noise level: 20.9 dBA), with no light during each test. elevations (F (2; 34) = 13:23; p < 0:001, F (1:436; 24:42) = 17:38; p < 0:001 and F (2; 34) = 16:25; p < 0:001): the Fig. 16 shows the results of the localization experiments generic method has significantly larger errors than the PCA of all eighteen subjects in three elevation angles respectively. and SPCA methods (p < 0:005 and p < 0:005 for all The judgments are plotted as a function of the coordinates of three elevations); the PCA method and SPCA method have the targets. The left column, the middle column and the right no significant difference (p = 1:0 for all three elevations). column depict the judgments using the generic method, the PCA method and the SPCA method respectively. There are 648 judgments shown in each panel, corresponding to the 54 B. Elevation Localization Experiments judgments made for each of the twelve binaural sounds. The purpose of this experiment is to compare the eleva- Note that the localization performance of the directions in tion localization performance of the HRIRs generated by the the test set of DNN modeling is as good as the directions in generic method, the PCA method and the SPCA method. the training set. The localization performances of the PCA method and the SPCA method are better than that of the In the CIPIC database, the azimuth angle and the elevation generic method. More judgments fall upon the diagonal line angle are measured in a head-centered interaural-polar coor- in the PCA and the SPCA method, which indicates a higher dinate system shown in Fig. 17 (a). However, the elevation precision of localization. The reconstruction error of ITD is angle in a spherical coordinate system is much closer to our larger than the just noticeable differences (JNDs) for some auditory perception. Therefore, the azimuth angle and the specific directions. For example, a few subjects perceived that elevation angle in the CIPIC database are represented in the Judgement Angle (Deg) Judgement Angle (Deg) Judgement Angle (Deg) IEEE/ACM TRANSACTIONS ON AUDIO, SPEECH AND LANGUAGE PROCESSING, VOL. 28, NO. 1, 2020 11 -30 -60 -90 -30 (a) (b) -60 Fig. 17. Definition of (a) interaural-polar coordinate system and (b) head- -90 -42 -19 0 19 47 70 90 -42 -19 0 19 47 70 90 -42 -19 0 19 47 70 90 related spherical coordinate system. Target Angle (Deg) Target Angle (Deg) Target Angle (Deg) Fig. 18. Judged direction versus target direction of all subjects using the head-related spherical coordinate system plotted in Fig. 17 (b). generic method (left column), the PCA method (middle column) and the SPCA The transformation formulas are as follows: method (right column) in azimuth of 30 degrees (top row) and 150 degrees (bottom row). The oblique line with a slope of 135 degrees corresponds to sin( ) = sin() cos('); the up-down confusions. (23) tan(' ) = tan(') /cos() ; TABLE V where  and ' refer to the azimuth angle and the elevation THE AVERAGE RESULTS OF THE ELEVATION LOCALIZATION angle in the head-related spherical coordinate system respec- EXPERIMENTS. 0 0 tively, and  and ' are the azimuth angle and the elevation Azimuth Up-down Angle of angle in the interaural-polar coordinate system respectively. Method angle (Deg) confusion rate (%) error (Deg) In the elevation localization experiments, we use the same Generic 16.1  12.1 18.69  6.42 stimulus in the Section V-A. The HRIRs of seven elevation 30 PCA 18.8  14.2 16.65  3.75 SPCA 20.4  10.3 15.04  4.13 angles, 90, 70, 47, 19, 0, -19 and -42 degrees, in two azimuth Generic 17.2  9.3 25.00  8.96 angles, nearly 30 and 150 degrees, are generated by the generic 150 PCA 18.0  13.0 21.86  7.86 method, the PCA method and the SPCA method respectively. SPCA 16.7  14.8 20.15  9.02 Then, the stimulus is filtered by the HRIRs produced by the three methods to create three kinds of sounds. A total of two elevation localization experiments are per- The judgments are plotted as a function of the coordinates of formed, and each experiment includes three tests. The two the targets. The left column, the middle column and the right column depict the judgments using the generic method, the experiments correspond to two azimuth angles, 30 and 150 PCA method and the SPCA method respectively. There are degrees, respectively. The order of the two experiments is 378 judgments shown in each panel, corresponding to the 54 balanced. Each test is a kind of sound generated by one judgments made to each of the seven binaural sounds. method, and the order of the three tests is arranged by latin The up-down confusion rate and error of angle with their square design across every three subjects. Before each test, the standard deviations of the elevation localization experiments subject is trained using the test sound of five elevation angles, are shown in Table V. Angle of error is the angle difference 90, 58, 30, 0 and -30 degrees. Through listening to these between the target angle and the judgement angle. Up-down sounds, the subject gradually builds up the spatial perception confusion is corrected before calculating the difference. The for this kind of virtual sound. After that, twenty-one binaural sounds are randomly played to the subject by a Sennheiser repeated-measures ANOVA was used to compare the localiza- HD 650 headphone through a Sound Blaster sound card. The tion performance (confusion rate and error of angle) of the twenty-one sounds contain seven directions’ sounds, and each three methods (generic, PCA and SPCA) for two azimuths direction appears three times. The subject can listen to one (30 and 150 degrees). According to the Shapiro-Wilk test, sound many times until he/she can identify the exact direction. all groups of data obey the normal distribution (p > 0:05). The subject gave the exact direction of each sound he/she Through the Mauchly’s spherical hypothesis test, the variance perceived during the test through an interface on a computer, covariance matrix of the dependent variable is dominantly which is shown in Fig. 15 (b). The blue dashed line in the equal (p > 0:05 except for the condition of error of angle in interface is moved with the mouse cursor and indicates the 30 degrees,  (2) = 9:11; p < 0:05). The data, which doesn’t exact angle. After each experiment, there were five minutes satisfy spherical assumption, are corrected by Greenhouse- for a break. geisser Method ( = 0:697). For confusion rates, there is Eighteen subjects (the same as in the azimuth localization no significant difference for three methods in two azimuths experiments) with normal hearing took part in the experiments. (F (2; 34) = 0:70; p = 0:51 and F (2; 34) = 0:08; p = 0:92). All experiments were performed in a sound booth (Back- For error of angle, there are significant differences for three methods in two azimuths (F (1:39; 23:7) = 3:92; p < 0:05 and ground noise level: 20.9 dBA), with no light during each test. Fig. 18 shows the results of the localization experiments F (2; 34) = 4:76; p < 0:05). With Bonferroni post-hoc test, the of all eighteen subjects in two azimuth angles respectively. generic method has significantly larger errors than the SPCA Judgement Angle (Deg) Judgement Angle (Deg) IEEE/ACM TRANSACTIONS ON AUDIO, SPEECH AND LANGUAGE PROCESSING, VOL. 28, NO. 1, 2020 12 method (p < 0:05 for two azimuths); there is no significant [7] F. Ma, J. Wu, M. Huang, W. Zhang, W. Hou, and C. Bai, “Finite element determination of the head-related transfer function,” Journal difference between the PCA method and the generic and SPCA of Mechanics in Medicine and Biology, vol. 15, no. 05, p. 1550066, methods for two azimuths (PCA VS generic: p = 0:73 and p = 0:15; PCA VS SPCA: p = 0:45 and p = 1:0). [8] T. Xiao and Q. Liu, “Finite difference computation of head-related transfer function for human hearing.” Journal of the Acoustical Society Based on both azimuth and elevation localization experi- of America, vol. 113, no. 5, pp. 2434–41, 2003. ments, we conclude that both the SPCA and PCA methods are [9] C. P. Brown and R. O. Duda, “A structural model for binaural sound superior than the generic method. For most of the localization synthesis,” IEEE transactions on speech and audio processing, vol. 6, no. 5, pp. 476–488, 1998. experiments, the SPCA method and the PCA method have [10] J. C. Middlebrooks, “Individual differences in external-ear transfer similar performances, except for the confusion rate in the functions reduced by scaling in frequency.” Journal of the Acoustical elevation of 22.5 degrees, where PCA method has significantly Society of America, vol. 106, no. 3, pp. 1480–1492, 1999. [11] D. N. Zotkin, J. Hwang, R. Duraiswaini, and L. S. Davis, “Hrtf smaller confusion rate than the SPCA method. However, personalization using anthropometric measurements,” in 2003 IEEE the proposed SPCA method predicts the HRTFs outside the Workshop on Applications of Signal Processing To Audio and Acoustics, training data well, which indicates that the SPCA method can New Paltz, NY, USA, 2003, pp. 157–160. [12] C. Jin, P. Leong, J. Leung, A. Corderoy, and S. Carlile, “Enabling in- predict the HRTFs in unmeasured spatial directions. dividualized virtual auditory space using morphological measurements,” in Proceedings of the First IEEE Pacific-Rim Conference on Multimedia, Sydney, Australia, 2000, pp. 235–238. VI. CONCLUSION [13] H. Hu, L. Zhou, H. Ma, and Z. Wu, “Hrtf personalization based on artificial neural network in individual virtual auditory space,” Applied This paper proposed the individual HRTFs modeling method Acoustics, vol. 69, no. 2, pp. 163–172, 2008. using DNN based on SPCA. By modeling the SPCs, the SPCA [14] C. Guezenoc and R. Seguier, “Hrtf individualization: A survey,” in Audio weights, the H and the ITDs respectively, we reconstruct the Engineering Society Convention 145. Audio Engineering Society, 2018. av [15] K. J. Fink and L. Ray, “Tuning principal component weights to indi- individual HRIRs in arbitrary spatial directions. The objective vidualize hrtfs,” in IEEE International Conference on Acoustics, Speech and subjective experiments are performed to evaluate the and Signal Processing, Kyoto, Japan, 2012, pp. 389–392. individual HRTFs generated by the proposed method, the [16] Y. Luo, D. N. Zotkin, and R. Duraiswami, “Virtual autoencoder based recommendation system for individualizing head-related transfer func- PCA method, and the generic method. Both the objective and tions,” in IEEE Workshop on Applications of Signal Processing to Audio subjective experiments’ results show that the proposed SPCA and Acoustics, New York, USA, 2013, pp. 1–4. method and PCA method are superior than the generic method. [17] B. Xie, “Recovery of individual head-related transfer functions from a small set of measurements.” Journal of the Acoustical Society of The spectral distortion of the SPCA method is significantly America, vol. 132, no. 1, pp. 282–294, 2012. smaller than PCA method in high frequencies but significantly [18] G. D. Romigh, D. S. Brungart, R. M. Stern, and B. D. Simpson, larger in low frequencies, except for few frequencies with no “Efficient real spherical harmonic representation of head-related transfer functions,” IEEE Journal of Selected Topics in Signal Processing, vol. 9, significant difference. The evaluation using the localization no. 5, pp. 921–930, 2015. model shows that the PCA method is better than the SPCA [19] J. O. Smith, “Techniques for digital filtering design and system identi- method. The subjective experiments show the SPCA method fication with the violin,” Ph.D. dissertation, CCRMA, Stanford, 1983. and the PCA method have similar performances, except for [20] I. Jolliffe, “Principal component analysis,” Technometrics, vol. 45, no. 3, p. 276, 2003. the confusion rate in the elevation of 22.5 degrees, where PCA [21] D. J. Kistler and F. L. Wightman, “A model of head-related transfer method has a better performance than the SPCA method. Nev- functions based on principal components analysis and minimum-phase ertheless, the results indicate that our proposed SPCA method reconstruction,” Journal of the Acoustical Society of America, vol. 91, no. 3, p. 1637, 1992. could predict the HRTFs in arbitrary spatial directions well. [22] M. Zhang, R. A. Kennedy, T. D. Abhayapala, and W. Zhang, “Statistical For the future work, we will use more HRTF data to implement method to identify key anthropometric parameters in hrtf individual- the individual HRTF models for better performance. ization,” in The Workshop on Hands-Free Speech Communication & Microphone Arrays, Edinburgh, UK, 2011, pp. 213–218. [23] X. Zeng, S. Wang, and L. Gao, “A hybrid algorithm for selecting head-related transfer function based on similarity of anthropometric REFERENCES structures,” Journal of Sound & Vibration, vol. 329, no. 19, pp. 4093– 4106, 2010. [1] J. Blauert, Spatial hearing: the psychophysics of human sound localiza- [24] R. Nicol, V. Lemaire, A. Bondu, and S. Busson, “Looking for a relevant tion. MIT press, 1997. similarity criterion for hrtf clustering: a comparative study,” in Audio [2] E. Wenzel, F. Wightman, and D. Kistler, “Localization with non- Engineering Society Convention 120. Paris, France: Audio Engineering individualized virtual acoustic display cues,” in Sigchi Conference on Society, 2006, p. 6653. Human Factors in Computing Systems, 1991, pp. 351–359. [25] R. B. Palm, “Prediction as a candidate for learning deep hierarchical [3] E. M. Wenzel, M. Arruda, D. J. Kistler, and F. L. Wightman, “Lo- models of data,” Master’s thesis, 2012. calization using nonindividualized head-related transfer functions,” The Journal of the Acoustical Society of America, vol. 94, no. 1, pp. 111– [26] V. R. Algazi, C. Avendano, and R. O. Duda, “Estimation of a spherical- 123, 1993. head model from anthropometry,” Journal of the Audio Engineering [4] V. R. Algazi, R. O. Duda, D. M. Thompson, and C. Avendano, “The Society, vol. 49, no. 6, pp. 472–479, 2001. cipic hrtf database,” in 2001 IEEE Workshop on the Applications of [27] C. I. Cheng and G. H. Wakefield, “Spatial frequency response surfaces Signal Processing to Audio and Acoustics. New Platz, NY, USA,: (sfr’s): An alternative visualization and interpolation technique for head- IEEE, 2001, pp. 99–102. related transfer functions (hrtf’s),” in Audio Engineering Society Con- [5] T. Qu, Z. Xiao, M. Gong, Y. Huang, X. Li, and X. Wu, “Distance- ference: 16th International Conference: Spatial Sound Reproduction. dependent head-related transfer functions measured with high spatial Audio Engineering Society, 1999. resolution using a spark gap,” IEEE Transactions on Audio, Speech, [28] S. Prepelia, ˘ M. Geronazzo, F. Avanzini, and L. Savioja, “Influence of and Language Processing, vol. 17, no. 6, pp. 1124–1132, 2009. voxelization on finite difference time domain simulations of head-related [6] W. Kreuzer, P. Majdak, and Z. Chen, “Fast multipole boundary element transfer functions,” The Journal of the Acoustical Society of America, method to calculate head-related transfer functions for a wide frequency vol. 139, no. 5, pp. 2489–2504, 2016. range,” Journal of the Acoustical Society of America, vol. 126, no. 3, [29] M. S. A. Zilany, I. C. Bruce, and L. H. Carney, “Updated parameters pp. 1280–90, 2009. and expanded simulation options for a model of the auditory periphery,” IEEE/ACM TRANSACTIONS ON AUDIO, SPEECH AND LANGUAGE PROCESSING, VOL. 28, NO. 1, 2020 13 The Journal of the Acoustical Society of America, vol. 135, no. 1, pp. 283–286, 2014. [30] P. Majdak, R. Baumgartner, and B. Laback, “Acoustic and non-acoustic factors in modeling listener-specific performance of sagittal-plane sound localization,” Frontiers in psychology, vol. 5, p. 319, 2014. [31] P. L. Søndergaard and P. Majdak, “The auditory modeling toolbox,” in The technology of binaural listening. Springer, 2013, pp. 33–56. [32] J. C. Middlebrooks, “Virtual localization improved by scaling nonindi- vidualized external-ear transfer functions in frequency,” The Journal of the Acoustical Society of America, vol. 106, no. 3, pp. 1493–1510, 1999. [33] B. Masiero and J. Fels, “Perceptually robust headphone equalization for binaural reproduction,” in Audio Engineering Society Convention 130. Audio Engineering Society, 2011. http://www.deepdyve.com/assets/images/DeepDyve-Logo-lg.png Electrical Engineering and Systems Science arXiv (Cornell University)

Modeling of Individual HRTFs based on Spatial Principal Component Analysis

Loading next page...
 
/lp/arxiv-cornell-university/modeling-of-individual-hrtfs-based-on-spatial-principal-component-d9K2snYDYV

References (34)

ISSN
2329-9290
eISSN
ARCH-3348
DOI
10.1109/TASLP.2020.2967539
Publisher site
See Article on Publisher Site

Abstract

IEEE/ACM TRANSACTIONS ON AUDIO, SPEECH AND LANGUAGE PROCESSING, VOL. 28, NO. 1, 2020 1 Modeling of Individual HRTFs based on Spatial Principal Component Analysis Mengfan Zhang, Zhongshu Ge, Tiejun Liu, Xihong Wu and Tianshu Qu Abstract—Head-related transfer function (HRTF) plays an multimedia, and virtual reality. Measuring the high spatial important role in the construction of 3D auditory display. This resolution HRTFs for each potential user is difficult, so schol- paper presents an individual HRTF modeling method using deep ars basically use non-individual HRTFs. However, using non- neural networks based on spatial principal component analysis. individual HRTFs may lead to some perception errors such as The HRTFs are represented by a small set of spatial principal in-head localization, front-back confusion, and a breakdown components combined with frequency and individual-dependent weights. By estimating the spatial principal components using of elevation discrimination ability [2], [3]. Thus, attaining deep neural networks and mapping the corresponding weights to individual HRTFs is very important and urgent in virtual a quantity of anthropometric parameters, we predict individual auditory scene synthesis. HRTFs in arbitrary spatial directions. The objective and subjec- At present, experimental measuring is an accurate method tive experiments evaluate the HRTFs generated by the proposed to obtain individual HRTFs. In the past two decades, a number method, the principal component analysis (PCA) method, and the generic method. The results show that the HRTFs generated of research groups have performed HRTF measurements and by the proposed method and PCA method perform better than established HRTF databases [4], [5]. However, experimental the generic method. For most frequencies the spectral distortion measuring of individual HRTFs requires rigorous experimental of the proposed method is significantly smaller than the PCA conditions and complicated equipment and keeps the subjects method in the high frequencies but significantly larger in the low not moving during the measuring procedure, which make this frequencies. The evaluation of the localization model shows the PCA method is better than the proposed method. The subjective method hard to implement. localization experiments show that the PCA and the proposed With the development of computer technology, numerical methods have similar performances in most conditions. Both calculation can be used to obtain HRTFs. Common numerical the objective and subjective experiments show that the proposed calculation methods include the boundary element method method can predict HRTFs in arbitrary spatial directions. (BEM) [6], the finite element method (FEM) [7] and the finite Index Terms—anthropometric parameters, HRTF, individual, difference method (FDM) [8]. However, numerical calculation SPCA. methods are computationally expensive and depend on the availability of precise 3D geometry. For example, magnetic resonance imaging (MRI) is used to obtain individual mor- I. I NTRODUCTION phology. This method requires a non-trivial acquisition process HE head-related transfer function (HRTF) describes the and complicated calculation. acoustic transmission of sound waves from a sound In recent years, many researchers have concentrated on source to a listener’s binaural ears. It assumes the head modeling individual HRTFs. Brown et al. [9] separated the exists in a free field. In time domain, HRTF is called head- effects of different physiological structures on HRTFs, and related impulse response (HRIR) [1]. HRTF has been widely each effect was modeled with a low-order sub-filter. The used in virtual sound technology, room acoustics simulation, combination of all sub-filters represents an HRTF. Middle- brooks [10] used the frequency scaling method, assuming that Copyright 2020 IEEE. Personal use of this material is permitted. Permission from IEEE must be obtained for all other uses, in any current or future media, the HRTF spectral characteristics of diverse individuals are including reprinting/republishing this material for advertising or promotional similar, but the corresponding frequencies of spectral charac- purposes, creating new collective works, for resale or redistribution to servers teristics are different. Through the frequency scaling method, or lists, or reuse of any copyrighted component of this work in other works. Manuscript received June 6, 2019; revised September 17, 2019 and De- a new subject’s HRTF can be obtained. Zotkin et al. [11] cember 15, 2019; accepted January 11, 2020. Date of publication January 17, selected the HRTF data of the subject whose anthropometric 2020; date of current version February 4, 2020. This work was supported by parameters were closest to the new subject. Jin et al. [12] the National Natural Science Foundation of China under Grant 11590773, Grant 61175043, and Grant 61421062, and in part by High-performance applied the principal component analysis (PCA) to the HRTF Computing Platform of Peking University. The associate editor coordinating amplitude spectrum and anthropometric parameters separately the review of this manuscript and approving it for publication was Prof. A. and then constructed a linear mapping from the PCA weights Sarti. (Corresponding author: Tianshu Qu.) M. Zhang, Z. Ge, X. Wu, and T. Qu are with the Key Labora- of the anthropometric parameters to the PCA weights of tory on Machine Perception (Ministry of Education), Speech and Hear- HRTFs. Hu et al. [13] used back-propagation artificial neural ing Research Center, Peking University, Beijing 100871, China (e-mail: networks to map the PCA weights of HRTFs to the selected zhangmengfan@pku.edu.cn; gezhongshu@pku.edu.cn; wxh@cis.pku.edu.cn; qutianshu@pku.edu.cn). anthropometric parameters. However, this approach required T. Liu is with the State key Laboratory of Robotics, Shenyang Institute separately training neural networks for each spatial direction. of Automation, Chinese Academy of Sciences, Shenyang, Liaoning 110004, To predict high spatial resolution HRTFs, this method had to China (e-mail: tjliu@sia.cn). Digital Object Identifier 10.1109/TASLP.2020.2967539 measure a large HRTF database and train thousands of neural arXiv:1910.09484v6 [eess.AS] 6 Feb 2020 IEEE/ACM TRANSACTIONS ON AUDIO, SPEECH AND LANGUAGE PROCESSING, VOL. 28, NO. 1, 2020 2 network models. XXX As the subjective perception of spatialization is the ultimate (f) = HRTF (; '; f; s); (2) log goal, the individual HRTFs can also be obtained based on S  D s ' the listeners feedback [14]. Fink et al. [15] let subjects tune After removing the mean logarithmic HRTFs, we obtain the the PCA weights from average HRTFs to obtain individual log-magnitude function, which is called HRTF . HRTFs. This tuning procedure can reduce localization errors; log however, obtaining a customized HRTF for a subject is very time-consuming. The subjects need a lot of time to finish the HRTF (; '; f; s) = HRTF (; '; f; s) (f): (3) log log tuning part. Luo et al. [16] also used the tuning method to obtain individual HRTFs and first introduced deep learning C. SPCA autoencoders to HRTF. The autoecoder was used to perform feature reduction and obtained a better result than PCA. PCA is mathematically defined as an orthogonal linear The aim of this paper is to realize the modeling of individual transformation that transforms the data to a new coordinate HRTFs and to predict the HRTFs in arbitrary spatial directions. system, such that the greatest variance by some projection Spatial principal component analysis (SPCA) [17] or spherical of the data comes to lie on the first coordinate (called the harmonics (SHs) basis functions [18] can be used to spatially first principal component), the second greatest variance on the decompose HRTFs. In this study, we use SPCA to decom- second coordinate, and so forth [20]. pose HRTF into a weighted combination of spatial principal The traditional PCA method is generally used in the time components (SPCs). Through the deep neural network (DNN) or frequency domain of HRTFs [21], [22]. In contrast to tra- training, the SPCs in arbitrary spatial directions are estimated. ditional PCA models, SPCA is applied to the spatial domain. A small quantity of anthropometric parameters were selected The high spatial resolution HRTFs can be represented as the and mapped to the SPCA weights using neural networks. Then, weighted combination of orthonormal SPCs [17]. The SPCA we combined the predicted SPCs and the SPCA weights for applied to a HRTF is shown below. log each individual to reconstruct the HRTFs in arbitrary spatial directions. HRTF (; '; f; s) = d (f; s)W (; ') + H (; '); log q q av The rest of the paper is organized as follows. In Section II, q=1 the data preprocessing and SPCA are performed. In Section (4) III, modeling individual HRTFs based on SPCA using DNN where q is the identification of the SPC, W is the SPC is described. In Section IV, the objective experiments are that depends only on the source direction, ' is the elevation conducted and the objective error is analyzed. In Section V, the angle, and  is the azimuth angle. D is the number of subjective evaluation of the proposed approach is described. spatial directions. d (f; s) is the SPCA weight which varies In Section VI, the conclusion is presented. as function of frequency f and individual s. H is the mean av HRTF magnitude across the frequencies and subjects and log II. SPCA AND DATA P REPROCESSING can be calculated as follows: XX A. CIPIC Database H (; ') = HRTF (; '; f; s); (5) av log The HRIRs used in this paper are derived from the CIPIC N  S s f database [4], which is measured by U. C. Davis CIPIC where N and S are the total number of frequencies and Interface Laboratory. In this database, a blocked ear technique subjects respectively. is performed for 45 subjects (27 males, 16 females, 2 KE- To calculate the SPCs and SPCA weights, we combine MARs), and 1250 directions of HRTF data were measured for HRTF of all the frequencies, directions and subjects into log each subject. Sound source directions are in interaural-polar an (NS) D matrix H . Each column of H corresponds to coordinate. The database also contains up to 27 anthropometric a spatial direction, and each row of H represents the HRTF parameters for each subject. of an individual at a discrete frequency. We subtract the mean value from H to obtain H . Then we calculate the covariance B. Data Preprocessing matrix R: We first transform the raw HRIRs in CIPIC database into T R = H H ; (6) the frequency domain. Fourier transformation is applied to the where R is a D  D matrix. Its eigenvectors are extracted HRIRs to calculate the HRTFs. Then we transform the HRTFs and arranged as the eigenvalue reduced-order. Then the first into a logarithmic scale. Because a logarithmic scale is much Q eigenvectors are taken as the base vectors, i.e. the SPCs closer to our auditory perception [19]. Therefore, the base 10 which represent the values of basis functions at D discrete log-magnitude responses of HRTFs are computed. directions, and form a Q D matrix W : HRTF (; '; f; s) = 20log (jHRTF (; '; f; s)j): (1) log 10 W = [W (0); W (1); :::; W (D 1)]: (7) q q q Then the mean logarithmic HRTFs is calculated from all the HRTF . The mean function includes the direction and Each row of W corresponds to a SPC, and each column of W log subject independent spectral features shared by all HRTFs in represents the values of SPCs at a spatial direction. We name the CIPIC database. the column of W as direction vector of SPCs (DV-SPCs). IEEE/ACM TRANSACTIONS ON AUDIO, SPEECH AND LANGUAGE PROCESSING, VOL. 28, NO. 1, 2020 3 TABLE I TABLE II CUMULATIVE PERCENTAGE OF VARIANCE IN RECONSTRUCTED HRTF S ELECTED KEY ANTHROPOM ETRIC PARAMETERS MAGNITUDES IS INCREASED W ITH THE NUMBER OF SELECTED SPCS Variable Measurement Variable Measurement Cumulative percentage x head width d cavum concha height 1 1 The number of SPCs (Q) of variance (%) x * head height d cavum concha width 2 3 Left ear Right ear x head depth d fossa height 3 4 1 16.54 20.14 x shoulder width d pinna height 12 5 5 52.20 55.33 d pinna width 10 62.29 64.84 * This parameter is only used in the prediction of ITDs. 20 70.10 71.85 50 78.33 79.54 60 80.09 81.22 individuals [13], [22]. We process the HRIRs and the anthro- 80 82.93 83.98 100 85.11 86.09 pometric parameters in CIPIC database as follows. 200 91.03 91.56 First, Fourier transformation is applied to the HRIRs in the 500 97.07 97.22 CIPIC database, and the mean of the obtained HRTF data is subtracted. For each sampled direction in the CIPIC database, we apply the traditional PCA to the HRTF data of all the Since all the SPCs are orthogonal to one another, by the use subjects, of these SPCs, we can obtain the SPCA weights: d = H W ; (8) HRTF (f; s) = d (s)W (f) + H (f); (11) q q av where d is a (NS)Q matrix composed of the SPCA weights q=1 for all the individuals and frequencies, W is the transpose where W is the PC. Note that the traditional PCA needs to of W . Finally, we can predict all the HRTFs as model HRTFs in each measured spatial direction, which means we apply PCA 1250 times in total. The resultant PCs and H = dW + H ; (9) AV H depend only on the frequency, whereas the SPCs and H av av where H is an (NS)  D matrix. Each row of H AV AV obtained by SPCA depend only on spatial directions. d (s) is corresponds to the H in D directions, and all of the rows av the PCA weight which varies as function of individual s. are identical. Second, multiple linear regression analysis is used to an- If the number of the SPCs is chosen to be equal to the alyze the relationship between the PCA weights and the number of spatial directions, i.e. Q = D, then the HRTFs anthropometric parameters. can be fully represented without loss, as shown in Eq.(4). If d = s + e; (12) Q < D is selected, then only an approximate representation of q the original HRTF magnitude can be obtained. The principal where is the regression coefficient, and e is the error. component variances are the eigenvalues of the covariance ma- Then we use t statistics to identify which parameters have trix R. The cumulative percentage of variance in reconstructed a significant effect on the PCA weights. HRTF is calculated as Third, Pearson correlation coefficient analysis is introduced to measure the strength of dependence between each of the q=1 V ar =  100%; (10) D two anthropometric parameters, q=1 (s s )(s s ) i i j j where  is an eigenvalue, and V ar is the cumulative per- r = p ; (13) ij 2 2 centage of variance. V ar is increased with the selection of Q, ( (s s ) (s s ) ) i i j j as shown in Table I. We can restore more than 70% of the where s is a vector that comprises 27 anthropometric pa- total variability by selecting the first 20 SPCs. We can recover rameters, i; j = 1; 2; :::; 27 and i 6= j. r is the degree of ij more than 80% of the total variability by selecting the first 60 correlation between two anthropometric parameters s and s . i j SPCs. Prior studies reported that more than 90% of the total This procedure is used to reduce the number of parameters to variability is enough to recover the HRTF magnitudes when make the measurement more feasible. A stronger correlation applying PCA in frequency domain [13], [23]. Therefore, we indicates that one parameter could be represented by another select the first 200 SPCs to recover more than 90% of the total [23]. variability when applying PCA in spatial domain. Finally, based on the correlation between the anthropomet- ric parameters as well as the analysis of the PCA weights D. Key Anthropometric Parameters Selection and anthropometric parameters, we balance the principle of There are 27 anthropometric parameters in the CIPIC simpleness, completeness and feasibility in practice to reduce database. Measuring all the anthropometric parameters for parameters. Anthropometric parameters with large correlation each individual is a time-consuming and difficult process. are reduced while considering the theoretical and practical in- In addition, not all of these anthropometric parameters are fluence of them on HRTFs to determine which to be remained. strongly related to HRTFs. Therefore, it is necessary to select Eight anthropometric parameters that have a significant effect a small set of anthropometric parameters which are most on the PCA weights are finally selected. The quantities of strongly related to the variations in the HRTFs of different selected anthropometric parameters are the same as [13], [22]. IEEE/ACM TRANSACTIONS ON AUDIO, SPEECH AND LANGUAGE PROCESSING, VOL. 28, NO. 1, 2020 4 Fig. 2. The framework of individual HRTF modeling. separately in order to keep a confident link between the HRTF dissimilarity and the angular difference [24]. Additionally, Fig. 1. Five selected pinna parameters in our model. there exists some differences between an individual’s two ears; we therefore consider one individual as containing two observations which should be modeled separately to obtain As shown in Table II, the eight anthropometric parameters different SPCA weights. When modeling the SPCs, the H av are head width, head depth, shoulder width, cavum concha and the ITDs, which are related to spatial directions, we split height, cavum concha width, fossa height, pinna height and the data of these parameters into two parts, the front and rear pinna width, corresponding to x , x , x , d , d , d , d 1 3 12 1 3 4 5 parts. In the CIPIC database, the azimuth angles composed a and d in the CIPIC database. Note that the head height vector azi = [80;65;55;45 : 5 : 45; 55; 65; 80], and parameter x is only used in the prediction of ITDs. The first the elevation angles make up a vector elev = 45+5:625[0 : four anthropometric parameters are head and torso parameters, 49]. The front part contains elev  90 in all of the azimuth which can be measured by a caliber rule or taking pictures. angles, and the remaining angles belong to the rear part. Thus, The last five anthropometric parameters are pinna parameters, each part contains a total of 1250 sets of data, including 625 which can be obtained by photograph annotation as shown in sets of data for two ears. For the three parameters, the SPCs, Fig. 1. the H and the ITDs, we train two DNNs for each parameter av in order to obtain a desirable prediction result. III. M ODELING OF I NDIVIDUAL HRTF S To sum up, through securing a quantity of anthropomet- ric parameters, we can reconstruct the individual’s binaural A. Outline HRIRs in arbitrary spatial directions. In this paper, the magnitude spectra of HRTFs are modeled based on SPCA using neural networks. The phase of HRTFs B. Prediction of SPCA Weights for an Individual is calculated by minimum-phase reconstruction [21]. Fig. 2 depicts the framework of individual HRTF modeling. Given that the SPCA weights are a function of frequency First we model the SPCA weights, the SPCs, the H and and an individual’s morphology, we model the SPCA weights av the ITDs, respectively. The SPCA weights can be obtained by and eight key anthropometric parameters based on neural modeling the individual’s morphology, for the SPCA weights networks. After securing eight anthropometric parameters of vary as functions of anthropometric parameters and frequency. an individual, we can estimate the SPCA weights for the A quantity of key anthropometric parameters is selected in individual in all of the frequency points using neural network Section II-D to represent the human morphology, therefore, the models. SPCA weights for any individual outside the database can be A total of 37 subjects s (m = 1; 2; :::; 37) in the CIPIC estimated from a small set of anthropometric measurements. database contain all of the eight anthropometric parameters. As DNNs are used to predict the SPCs, the H and the ITDs aforementioned in Section III-A, one individual’s two ears are av in arbitrary spatial directions. Accordingly, HRTFs with high modeled separately to obtain different SPCA weights. Thus, directional resolution can be recovered by solving the Eq. we have 74 sets of anthropometric parameters; considering the (4) using the predicted SPCs, SPCA weights and H . Then total number of frequency points is N = 200, we have 14800 av the minimum-phase reconstruction method is used to generate sets of the SPCA weights. mono HRIRs [21]. Finally, binaural HRIRs are obtained using For each frequency point f (k = 1; 2; :::; N), we build estimated ITDs and the corresponding left and right mono a model for the eight anthropometric parameters and the HRIRs. corresponding SPCA weights. The specific architecture of the Due to the symmetry between front and rear HRTFs, we neural network model is shown in Fig. 3. The inputs are model the front and rear HRTFs separately. Because HRTF the eight anthropometric parameters, and the ground truth are dissimilarity increases with their angular difference, however, SPCA weights in one frequency point. So in total we build the low dissimilarity may occur for large angular difference. 101 neural network models because of the symmetry property For instance, if the two HRTFs taken from two locations are of the Fourier transformation. symmetric with respect to the interaural axis, the two HRTFs 30 subjects’ binaural data are carefully selected to guarantee will be very similar despite the large angular difference. that the data sets for each anthropometric parameter are widely Thus, it is preferred to consider the front and rear HRTFs distributed. That is to say, most of the test data are within the IEEE/ACM TRANSACTIONS ON AUDIO, SPEECH AND LANGUAGE PROCESSING, VOL. 28, NO. 1, 2020 5 subject_033 average 0 2 4 6 8 10 12 14 16 18 20 Frequency (kHz) Fig. 5. The prediction errors of SPCA weights as a function of frequency. Fig. 3. The architecture of neural network model for SPCA weights prediction. 4. As previously discussed in Section II, the first 50 SPCs can restore about 80% of the total variability, and the first 200 SPCs can reconstruct above 90% of the total variability. Therefore, when the range of the q is from 51 to 200, the total amount of variability is around 10%. Moreover, we observe that the SPCA weights are approximately equal to 0 when q is larger than 50. Thus, we only plot the curves where q is less than 50 to make comparisons between the estimated SPCA weights and the real ones. Fig. 5 shows the prediction error, which is the mean of the absolute errors across all directions, of SPCA weights for subject 033 and the average prediction error. The figure shows that the difference between the real value and the estimated value of the SPCA weights increases when frequency increases. The SPCA weights in Fig. 4. Comparison of the real value (top) and the estimated value (bottom) high frequency are more unstable than the lower ones. for the first 50 orders of SPCA weights. C. Modeling of SPCs corresponding training data ranges [13]. These data, which As discussed in Section II, W is a QD matrix composed comprise a total of 60 sets, are used in the training phase, and of the first Q SPCs. Each row of W corresponds to an SPC, 10 of them are used as the validation set to prevent over-fitting. and each column of W represents the values of SPCs at The remaining 14 sets of data are used to test the average error a spatial direction. The DV-SPCs are a function of spatial of the neural network model. The mean and variance of the directions and represented by the column of W . We model test set and validation set are normalized using the training set the first Q SPCs by predicting DV-SPCs in all sampled D statistics to have zero mean and unit variance. All the neural spatial directions, then recombining W . network models comprise a single hidden layer, a hyperbolic 2 3 tangent activation function and a hyperbolic tangent output W (0); W (1) : : : W (D 1) 1 1 1 function, and they are feedforward backpropagation neural 6 7 W (0); W (1) : : : W (D 1) 2 2 2 6 7 networks with a learning rate of 0.001. W = 6 7 : (15) . . . . . . . . 4 5 For network models of all frequencies, each network takes . . . . about 500 epochs to converge averagely. After training all the W (0); W (1) : : : W (D 1) Q Q Q neural networks, we obtain a system for predicting the SPCA The DV-SPCs are modeled with DNNs. We set the angles weights. The mean square error (MSE) is used to test the of directly ahead (' = 0 ,  = 0 ) and the directly behind prediction system: (' = 180 ,  = 0 ) as the reference direction for the front XXX ^ and rear data sets, respectively. The inputs of DNN are the e = (d (f ; s ) d (f ; s )) ; d q k m q k m Q N  S DV-SPCs in the reference direction, the target azimuth  and q m (14) the target elevation ' in degrees. We set the ground truth as where d is estimated by neural network model, and e is the the DV-SPCs in the target direction. The specific architecture q d reconstruction error of the SPCA weights. The MSE of the of the DNN for predicting the DV-SPCs is shown in Fig. 6. overall prediction, including all the SPCA weights for all the Inputing the DV-SPCs in the reference direction reduces the individuals and frequencies, is 9.25. complexity of learning process. The main task of the DNN We randomly select a subject 033 in the CIPIC database and becomes learning the difference between the DV-SPCs in the plot the comparison of the estimated value and real value of reference direction and target direction. This can effectively its left ear’s first 50 orders of SPCA weights, as shown in Fig. improve the task result. prediction error IEEE/ACM TRANSACTIONS ON AUDIO, SPEECH AND LANGUAGE PROCESSING, VOL. 28, NO. 1, 2020 6 Fig. 8. Comparison of the real value (left) and the estimated value (right) of the H . av Fig. 6. The architecture of DNN for predicting the DV-SPCs. directions can be predicted. We combine the DV-SPCs in all the D sampled directions to obtain the QD matrix W . MSE is used to calculate the reconstruction error: XX e = (W (; ) W (; )) ; (16) W q q where W is estimated by the DNN model and e is the q W reconstruction error whose value is 1:97 10 . Fig. 7 illustrates the comparison of the first, the second, the third and the fourth SPCs’ real values and estimated values. The predicted SPCs are close to the real SPCs. As discussed in Section II, the first 5 SPCs can restore more than 50% of the total variability, and our algorithm performs well. With an increment of q, the corresponding SPC gradually becomes unstable, and the difference between the predicted SPC and the real SPC widens. Fig. 7. Comparison of the real value (top row) and the estimated value (bottom row) of the first 4 SPCs (4 columns). D. Modeling of H av The modeling of the H is similar to the prediction of the av To guarantee the variability of the test data, we index the DV-SPCs. This model is also based on DNN. The input of 1250 sets of data and select the test data every four indexes. DNN is the H in the reference direction, the target azimuth av The extra directions of data are then re-indexed and we choose and the target elevation in degrees. The ground truth is the H av the validation data every five indexes. The remaining directions in the target direction. The reference directions in addition to of data are used as training data. Thus, we uniformly distribute the selection of training set, validation set and test set are the the training set, validation set and test set in the space. same as in Section III-C. The mean and variance of the test The validation set is used to control the training phase and set and validation set are normalized using the training set prevent over-fitting, and the test set is used to estimate the statistics to have zero mean and unit variance. Both of the generalization error of the modeling. The mean and variance modeled DNNs are fully connected and have three hidden of the test set and validation set were normalized using the layers. The activation function and the output function are training set statistics to have zero mean and unit variance. The set hyperbolic tangent, and the learning rate is 0.001. The DNNs were implemented by MATLAB ’s DeepLearnToolbox- architecture and the parameters are determined through the master [25]. Each DNN is fully connected and has three hidden results of many experiments. layers. Three hidden layers are chosen because it has a better It takes about 700 epochs for the DNN models to converge. performance than using only one or two hidden layers; if more After training the DNNs, we can obtain a system to predict than three hidden layers are used, the performance does not the H in arbitrary spatial directions. All the D directions in av have a greater improvement, and the neural network is also the CIPIC database are used as the target directions to predict easy to overfit. Both the activation function and the output the corresponding H . The predicted results are shown in Fig. av function are set hyperbolic tangents with the learning rate set as 0.001. These settings which were determined through The MSE is used to calculate the reconstruction error of the the results of many experiments can lead to good prediction H : av performance. XX It takes about 5000 epochs for the DNN models to converge. 2 e = (H (; ) H (; )) ; (17) H av av By training the DNN models, the DV-SPCs in arbitrary spatial IEEE/ACM TRANSACTIONS ON AUDIO, SPEECH AND LANGUAGE PROCESSING, VOL. 28, NO. 1, 2020 7 Fig. 10. Comparison of the real value (left) and the estimated value (right) of the ITD. Fig. 9. The architecture of DNN for predicting the ITDs. where T is the estimated ITD and T is the real ITD. e is the reconstruction error of the ITDs, and its value is 2:22 10 ms. where H is estimated by the DNN model and e is the av H reconstruction error whose value is 0.195. F. Recovery of Individual HRIRs To reconstruct an individual’s HRIRs, we first model the SPCA weights using the eight key anthropometric parameters. E. Modeling of ITDs Then DNN is used to model SPCs, H and ITDs respectively. av The azimuth angle and elevation angle are introduced into the As ITD not only relates to spatial directions, but also input layer of the DNN models, the DV-SPC, H and ITD in varies with different individuals, we use the anthropometric av arbitrary spatial directions can then be predicted. parameters and spatial directions to model it. The inputs of The HRTF magnitude of a new subject s in azimuth angle DNN include head width x , head height x , as well as head m 1 2 and elevation angle ' can be reconstructed by solving depth x of an individual, target azimuth and target elevation in d d the Eq. (9). The SPCA weights for the new subject can be degrees. The individual’s ITD in the target direction is taken combined into a matrix. as the ground truth, as shown in Fig. 9. x , x and x are 1 2 3 2 3 chosen by many experiments. The three parameters together d (f ; s ); d (f ; s ) : : : d (f ; s ) 1 1 m 2 1 m Q 1 m show the best results, and adding other parameters does not 6 7 d (f ; s ); d (f ; s ) : : : d (f ; s ) 1 2 m 2 2 m Q 2 m 6 7 d = ; yield an improvement. This indicates the strong correlation 6 . . . . 7 . . . . 4 5 . . . . between the ITDs and the head size [26]. d (f ; s ); d (f ; s ) : : : d (f ; s ) 1 N m 2 N m Q N m Similar to Section III-C, 625 directions of ITD data for (19) an individual are indexed and the test data are selected every the DV-SPCs in  and ' is d d four indexes. The extra directions of data are then re-indexed W = [W ( ; ' ); W ( ; ' ); :::; W ( ; ' )] ; (20) and we choose the validation data every five indexes. The 1 d d 2 d d Q d d remaining directions of data are used as training data. In the the H in  and ' is a vector of length N : av d d end, the training data of the 30 subjects, selected in Section III-B, are combined as the training set. The validation data of H = [H ( ; ' ); H ( ; ' ); :::; H ( ; ' )] ; (21) AV av d d av d d av d d the 30 subjects and the test data of the remaining 7 subjects are then the HRTF magnitude in  and ' of the new subject d d integrated as the validation set and the test set, respectively. can be restored. The mean and variance of the test set and validation set The minimum phase reconstruction method is then em- are normalized using the training set statistics to have zero ployed to the HRTF magnitudes to generate the mono HRIRs mean and unit variance. Each DNN is fully connected and has [21]. The ITDs obtained in Section III-E are used to calculate three hidden layers to yield a better result. Both the activation the binaural HRIRs. The individual’s HRIRs in arbitrary function and the output function are set hyperbolic tangents spatial directions can then be obtained. and the learning rate is 0.001. It takes about 5500 epochs for the DNN models to converge. IV. OBJECTIVE EXPERIM ENTS Fig. 10 describes the results of ITD prediction. A total A. Evaluation using spectral distortion of 1250 spatial directions for a subjects’ ITDs are plotted compared with its real value. To evaluate the effectiveness of our proposed approach, we carried out a set of objective experiments for the proposed The mean absolute error (MAE) is used to evaluate the SPCA method, the PCA method and the generic method. reconstruction error of ITDs. The generic method used the HRTFs of the CIPIC KEMAR XXX with small ears. The PCA method applied traditional PCA to e = jT (; ; s ) T (; ; s )j; (18) T m m S  D HRTFs in each sampled spatial direction, which is described IEEE/ACM TRANSACTIONS ON AUDIO, SPEECH AND LANGUAGE PROCESSING, VOL. 28, NO. 1, 2020 8 Fig. 11. The SFRSs of the real value (1st column), the generic method(2nd column), the PCA method (3rd column) and the SPCA method (4th column) 2 at frequency of 12.35 kHz of subject 163. Generic PCA SPCA 1.5 3.3 5.1 6.8 8.6 10.4 12.1 13.9 15.7 17.4 19.2 20.9 Frequency (kHz) Fig. 13. Comparison of the generic method, the PCA method and the SPCA method for the mean SD across all directions for some frequency bins and their deviations. Fig. 12. The average prediction errors of SFRSs across all frequencies and subjects in the test set of the generic method (left), the PCA method (middle) and the SPCA method (right). Fig. 14. The localization performance of the median plane using the localization model with the generic method (left), the PCA method (middle) in Section II-D. The first twelve principal components (PCs) and the SPCA method (right). for all the PCA models are selected. When we apply PCA to HRTFs, for all of the spatial directions, we can restore an average of 92.02% of the total variability for the left ear HRTF data. and 91.71% for the right ear. These percentages of the total XX variability are close to our proposed SPCA model’s selections. SD(f; s) = jH(; ; f; s) H(; ; f; s)j; (22) Then, neural network models are used to map the selected eight anthropometric parameters to the PCA weights. Thus, a where H is the magnitude response (dB) of the measured total of 1250 neural network models are trained to yield the HRTF from the CIPIC database, and H is the magnitude PCA weights in each spatial direction. This is quite arduous response (dB) of the test HRTF. in PCA method, which will bring computation issues in the Fig. 13 shows how the SD varies with frequency of the three applications of HRTF prediction. Finally, individual HRTFs method. Slight translations were performed in frequency axis can be reconstructed by combining the PCs and the PCA in order to avoid overlapping of symbols for different methods. weights. As a result, the proposed method has the lowest SD for most The Spatial Frequency Response Surfaces (SFRSs) [27], frequencies compared to the other two methods. The average a spatial-domain visualization tool for HRTFs, is used to SD of the proposed model in all the sampled directions is evaluate the three methods. Each frequency bin in the HRTF 5.54 dB, which is 1.13 dB lower than that of the generic left or right magnitude responses constructs one SFRS, where method and 0.29 dB lower than that of the PCA method. T- magnitude is plotted against azimuth and elevation. Fig. 11 test applied on the 12 frequency bins in Fig. 13 shows that shows the SFRS of the three methods compared with the real SPCA method has significant smaller SD (p < 0:001) than value of subject 163 at frequency of 12.35 kHz. Fig. 12 shows PCA method above the frequency of 6 kHz but has significant the average prediction errors of SFRS across all frequencies larger SD (p < 0:001) than PCA method blow the frequency and subjects in the test set of the three methods. The SPCA of 6 kHz, except for 3 frequency bins, 3526 Hz (p = 0:052), method has smaller and more averagely distributed prediction 10584 Hz (p = 0:157) and 21168 Hz (p = 0:073) with no errors than the other two methods. Namely, our proposed significant difference. The results indicate that our proposed SPCA method has a much more smooth spatial shape than SPCA method predicted the HRTFs in the test set well. Since PCA method, which indicates that our algorithm predict the the test sets in the modeling of DV-SPCs, H and ITDs are av HRTF in spatial domain well. different from the training spatial directions, our proposed Then a frequency dependent spectral distortion (SD) [28] model will definitely predict the HRTFs in unmeasured spatial is used as an evaluation metric between the real and the test directions. SD magnitude (dB) IEEE/ACM TRANSACTIONS ON AUDIO, SPEECH AND LANGUAGE PROCESSING, VOL. 28, NO. 1, 2020 9 TABLE III THE AVERAGE RESULTS OF THE OBJ ECTIVE EVALUATION ON MEDIAN PLANE LOCALIZATION EXPERIMENT. Confusion rate (%) Angle of PE Method Up-down Front-back error (Deg) (Deg) Generic 11.05 16.32 17.71  2.48 26.84 PCA 7.81 11.22 13.93  2.21 21.63 SPCA 12.62 14.51 15.91  2.53 24.89 (a) Interface for azimuth experiments (b) Interface for elevation experiments B. Evaluation using the localization model Fig. 15. The user interfaces for subjects to give percepted directions of For further evaluation, an auditory based localization model sounds. [29], [30] was used to evaluate the localization performance of the generic, the PCA and SPCA calculated HRIRs. The BAUMGARTNER2014 function in the AMT tool box [31] subjects’ anthropometric parameters listed in Table II. These was used in our experiment. The BAUMGARTNER2014 func- parameters were measured before the experiments. The 4 tion is a sagittal plane localization model. As a template-based variables listed in the 1st column were measured using a comparison model, it can be used to evaluate the performance ruler with an accuracy of 2 mm. The 5 variables listed in of HRTF individualization methods. The DTFs [32] were the 3rd column were measured through photographing and extracted from HRTFs [30] as the inputs of the auditory model. manually demarcating them as shown in Fig. 1. The ruler was We apply an listener-specific sensitivity threshold of 1 and set used to obtain the scale between the picture and real objects. the differential order of the spectral gradient extraction to be All parameters were measured three times then averaged to 0 to acquire a reasonable prediction error. Other settings of reduce the measurement errors. The impulse responses from the function were set as default. the headphone, used in the experiments, to the entrances We applied the localization model on the median plane to of the blocked ear canals of subjects were measured using the BK Type 4101-A binaural in-ear microphone. During the test the performance of the generic method, the PCA method experiments, headphone equalization was performed [33]. and the SPCA method. For each method, the HRIRs of all the 37 individuals used in our paper were tested. For each elevation angle, 2 runs were conducted to reduce probable A. Azimuth Localization Experiments prediction errors. Namely, for each elevation angle of the The aim of this experiment is to compare the azimuth median plane, there were 74 predicted elevation angles of the localization performance of the HRTFs generated by the target HRIRs in total. generic method, the PCA method and the SPCA method. For Fig. 14 shows the localization performance of the generic the PCA method, the phases of the HRTFs are obtained by method, the PCA method and the SPCA method of the median the minimum phase reconstruction method, and the ITDs are plane. Table III shows the statistical analysis of the experiment yielded using the modeling method of Section III-E. results. Four indices, including the up-down confusion rate, the The stimulus in this experiment was a train of eight 250-ms front-back confusion rate, the angle of errors with standard bursts of Gaussian noise (20-ms cosine-squared onset-offset deviations and the polar RMS (root-mean square) error (PE), ramps), with 300 ms of silence between the bursts. The HRIRs were used to evaluate the localization performance. For the an- of twelve azimuth angles, 0, 30, 55, 80, 125, 150, 180, 210, gle errors of the median plane, the repeated-measures ANOVA 235, 280, 305 and 330 degrees, in three elevation angles, 0, was used to verify the significance of the mean difference. 22.5 and 45 degrees, are generated by the generic method, the According to the Shapiro-Wilk test, the data of each group PCA method and the SPCA method respectively. Note that obey the normal distribution (p > 0:05). After the Mauchly’s the HRIRs of the four azimuth angles, 55, 150, 210 and 305 spherical hypothesis test, the variance covariance matrix of degrees, in the three elevation angles are not in the training set the dependent variable is not equal ( (2) = 6:279; p = and the validation set of each DNN model, i.e. Section III-C, 0:043 < 0:05). The data are corrected by Huynh-Feldt Method Section III-D and Section III-E. Then, the stimulus is filtered ( = 0:898). The mean errors of the three methods have by the HRIRs produced by the three methods to create three significant difference (F (1:80; 64:65) = 34:37; p < 0:001). kinds of sounds. Bonferroni post-hoc test shows that the SPCA and PCA A total of three azimuth localization experiments are per- methods are significantly better than the generic method (both formed, and each experiment includes three tests. The three p < 0:001). Further more, the PCA method is better than the experiments correspond to three elevation angles, 0, 22.5 and SPCA method (p < 0:001). 45 degrees, respectively. Each test is a kind of sound generated by one method, and the order of the three tests is arranged by V. SUBJ ECTIVE EXPERIM ENTS latin square design across every three subjects. Before each In this section, the azimuth localization experiments and the test, the subject is trained using the test sound of the other elevation localization experiments were conducted to evaluate eight azimuth angles, 0, 45, 90, 135, 180, 225, 270 and 315 the performance of the proposed method. In the experiments, degrees. Through listening to these sounds, the subject can three methods were used to generate the HRTFs according to build up the spatial perception for this kind of virtual sound. IEEE/ACM TRANSACTIONS ON AUDIO, SPEECH AND LANGUAGE PROCESSING, VOL. 28, NO. 1, 2020 10 TABLE IV THE AVERAGE RESULTS OF THE AZIMUTH LOCALIZATION EXPERIM ENTS. Elevation Front-back Angle of Method angle (Deg) confusion rate (%) error (Deg) Generic 28.2  9.6 18.90  4.25 0 PCA 22.4  11.1 14.84  3.63 SPCA 24.5  12.7 15.15  4.04 Generic 29.0  12.2 19.49  4.22 22.5 PCA 17.0  9.9 15.25  3.25 SPCA 25.6  13.3 14.96  3.74 180 Generic 31.8  11.9 19.79  3.30 45 PCA 25.3  10.8 15.74  3.81 SPCA 29.0  10.7 15.28  3.75 azimuth of 0 degrees is not directly ahead. The front-back confusion rate and error of angle with their standard deviations of the azimuth localization experiments are shown in Table IV. Angle of error is the angle difference between the target angle and the judgement angle. Front- 0 60 120 180 240 300 0 60 120 180 240 300 0 60 120 180 240 300 back confusion is corrected before calculating the difference. Target Angle (Deg) Target Angle (Deg) Target Angle (Deg) The repeated-measures ANOVA was used to compare the Fig. 16. Judged direction versus target direction of all subjects using the localization performance (confusion rate and error of angle) generic method (left column), the PCA method (middle column) and the SPCA of the three methods (generic, PCA and SPCA) for three method (right column) in elevations of 0 degrees (top row), 22.5 degrees elevations (0, 22.5 and 45 degrees). According to the Shapiro- (middle row) and 45 degrees (bottom row). Two oblique lines with a slope of 135 degrees correspond to the front-back confusions. Wilk test, all groups of data obey the normal distribution (P > 0:05). Through the Mauchly’s spherical hypothesis test, the variance covariance matrix of the dependent variable is After that, thirty-six binaural sounds are randomly played to dominantly equal (p > 0:05 except for the condition of error the subject by a Sennheiser HD 650 headphone through a of angle in 22.5 degrees,  (2) = 7:97; p < 0:05). The Sound Blaster sound card. The thirty-six sounds contain twelve data, which doesn’t satisfy spherical assumption, are corrected directions’ sounds, and each direction appears three times. The by Greenhouse-geisser Method ( = 0:718). For confusion subject can listen to one sound many times until he/she can rate, there is no significant difference for three methods in identify the exact perceived direction. The subject gave the elevation of 0 and 45 degrees (F (2; 34) = 2:41; p = 0:11 exact direction of each sound he/she perceived during the test and F (2; 34) = 2:10; p = 0:14), but a significant difference through an interface on a computer, which is shown in Fig. in elevation of 22.5 degrees (F (2; 34) = 8:71; p < 0:005). 15 (a). The blue dashed line in the interface is moved with Bonferroni post-hoc test for elevation of 22.5 degrees shows the mouse cursor and indicates the exact angle. After each that the PCA method has significantly smaller confusion experiment, there were five minutes for a break. rate than the SPCA and generic methods (p < 0:05 and Eighteen subjects (11 male 7 female, age from 21 to 27) p < 0:005); the confusion rates of the SPCA method and with normal hearing took part in the experiments. All of the the generic method have no significant difference (p = 0:92). experiments were performed in a sound booth (Background For error of angle, there are significant differences for all three noise level: 20.9 dBA), with no light during each test. elevations (F (2; 34) = 13:23; p < 0:001, F (1:436; 24:42) = 17:38; p < 0:001 and F (2; 34) = 16:25; p < 0:001): the Fig. 16 shows the results of the localization experiments generic method has significantly larger errors than the PCA of all eighteen subjects in three elevation angles respectively. and SPCA methods (p < 0:005 and p < 0:005 for all The judgments are plotted as a function of the coordinates of three elevations); the PCA method and SPCA method have the targets. The left column, the middle column and the right no significant difference (p = 1:0 for all three elevations). column depict the judgments using the generic method, the PCA method and the SPCA method respectively. There are 648 judgments shown in each panel, corresponding to the 54 B. Elevation Localization Experiments judgments made for each of the twelve binaural sounds. The purpose of this experiment is to compare the eleva- Note that the localization performance of the directions in tion localization performance of the HRIRs generated by the the test set of DNN modeling is as good as the directions in generic method, the PCA method and the SPCA method. the training set. The localization performances of the PCA method and the SPCA method are better than that of the In the CIPIC database, the azimuth angle and the elevation generic method. More judgments fall upon the diagonal line angle are measured in a head-centered interaural-polar coor- in the PCA and the SPCA method, which indicates a higher dinate system shown in Fig. 17 (a). However, the elevation precision of localization. The reconstruction error of ITD is angle in a spherical coordinate system is much closer to our larger than the just noticeable differences (JNDs) for some auditory perception. Therefore, the azimuth angle and the specific directions. For example, a few subjects perceived that elevation angle in the CIPIC database are represented in the Judgement Angle (Deg) Judgement Angle (Deg) Judgement Angle (Deg) IEEE/ACM TRANSACTIONS ON AUDIO, SPEECH AND LANGUAGE PROCESSING, VOL. 28, NO. 1, 2020 11 -30 -60 -90 -30 (a) (b) -60 Fig. 17. Definition of (a) interaural-polar coordinate system and (b) head- -90 -42 -19 0 19 47 70 90 -42 -19 0 19 47 70 90 -42 -19 0 19 47 70 90 related spherical coordinate system. Target Angle (Deg) Target Angle (Deg) Target Angle (Deg) Fig. 18. Judged direction versus target direction of all subjects using the head-related spherical coordinate system plotted in Fig. 17 (b). generic method (left column), the PCA method (middle column) and the SPCA The transformation formulas are as follows: method (right column) in azimuth of 30 degrees (top row) and 150 degrees (bottom row). The oblique line with a slope of 135 degrees corresponds to sin( ) = sin() cos('); the up-down confusions. (23) tan(' ) = tan(') /cos() ; TABLE V where  and ' refer to the azimuth angle and the elevation THE AVERAGE RESULTS OF THE ELEVATION LOCALIZATION angle in the head-related spherical coordinate system respec- EXPERIMENTS. 0 0 tively, and  and ' are the azimuth angle and the elevation Azimuth Up-down Angle of angle in the interaural-polar coordinate system respectively. Method angle (Deg) confusion rate (%) error (Deg) In the elevation localization experiments, we use the same Generic 16.1  12.1 18.69  6.42 stimulus in the Section V-A. The HRIRs of seven elevation 30 PCA 18.8  14.2 16.65  3.75 SPCA 20.4  10.3 15.04  4.13 angles, 90, 70, 47, 19, 0, -19 and -42 degrees, in two azimuth Generic 17.2  9.3 25.00  8.96 angles, nearly 30 and 150 degrees, are generated by the generic 150 PCA 18.0  13.0 21.86  7.86 method, the PCA method and the SPCA method respectively. SPCA 16.7  14.8 20.15  9.02 Then, the stimulus is filtered by the HRIRs produced by the three methods to create three kinds of sounds. A total of two elevation localization experiments are per- The judgments are plotted as a function of the coordinates of formed, and each experiment includes three tests. The two the targets. The left column, the middle column and the right column depict the judgments using the generic method, the experiments correspond to two azimuth angles, 30 and 150 PCA method and the SPCA method respectively. There are degrees, respectively. The order of the two experiments is 378 judgments shown in each panel, corresponding to the 54 balanced. Each test is a kind of sound generated by one judgments made to each of the seven binaural sounds. method, and the order of the three tests is arranged by latin The up-down confusion rate and error of angle with their square design across every three subjects. Before each test, the standard deviations of the elevation localization experiments subject is trained using the test sound of five elevation angles, are shown in Table V. Angle of error is the angle difference 90, 58, 30, 0 and -30 degrees. Through listening to these between the target angle and the judgement angle. Up-down sounds, the subject gradually builds up the spatial perception confusion is corrected before calculating the difference. The for this kind of virtual sound. After that, twenty-one binaural sounds are randomly played to the subject by a Sennheiser repeated-measures ANOVA was used to compare the localiza- HD 650 headphone through a Sound Blaster sound card. The tion performance (confusion rate and error of angle) of the twenty-one sounds contain seven directions’ sounds, and each three methods (generic, PCA and SPCA) for two azimuths direction appears three times. The subject can listen to one (30 and 150 degrees). According to the Shapiro-Wilk test, sound many times until he/she can identify the exact direction. all groups of data obey the normal distribution (p > 0:05). The subject gave the exact direction of each sound he/she Through the Mauchly’s spherical hypothesis test, the variance perceived during the test through an interface on a computer, covariance matrix of the dependent variable is dominantly which is shown in Fig. 15 (b). The blue dashed line in the equal (p > 0:05 except for the condition of error of angle in interface is moved with the mouse cursor and indicates the 30 degrees,  (2) = 9:11; p < 0:05). The data, which doesn’t exact angle. After each experiment, there were five minutes satisfy spherical assumption, are corrected by Greenhouse- for a break. geisser Method ( = 0:697). For confusion rates, there is Eighteen subjects (the same as in the azimuth localization no significant difference for three methods in two azimuths experiments) with normal hearing took part in the experiments. (F (2; 34) = 0:70; p = 0:51 and F (2; 34) = 0:08; p = 0:92). All experiments were performed in a sound booth (Back- For error of angle, there are significant differences for three methods in two azimuths (F (1:39; 23:7) = 3:92; p < 0:05 and ground noise level: 20.9 dBA), with no light during each test. Fig. 18 shows the results of the localization experiments F (2; 34) = 4:76; p < 0:05). With Bonferroni post-hoc test, the of all eighteen subjects in two azimuth angles respectively. generic method has significantly larger errors than the SPCA Judgement Angle (Deg) Judgement Angle (Deg) IEEE/ACM TRANSACTIONS ON AUDIO, SPEECH AND LANGUAGE PROCESSING, VOL. 28, NO. 1, 2020 12 method (p < 0:05 for two azimuths); there is no significant [7] F. Ma, J. Wu, M. Huang, W. Zhang, W. Hou, and C. Bai, “Finite element determination of the head-related transfer function,” Journal difference between the PCA method and the generic and SPCA of Mechanics in Medicine and Biology, vol. 15, no. 05, p. 1550066, methods for two azimuths (PCA VS generic: p = 0:73 and p = 0:15; PCA VS SPCA: p = 0:45 and p = 1:0). [8] T. Xiao and Q. Liu, “Finite difference computation of head-related transfer function for human hearing.” Journal of the Acoustical Society Based on both azimuth and elevation localization experi- of America, vol. 113, no. 5, pp. 2434–41, 2003. ments, we conclude that both the SPCA and PCA methods are [9] C. P. Brown and R. O. Duda, “A structural model for binaural sound superior than the generic method. For most of the localization synthesis,” IEEE transactions on speech and audio processing, vol. 6, no. 5, pp. 476–488, 1998. experiments, the SPCA method and the PCA method have [10] J. C. Middlebrooks, “Individual differences in external-ear transfer similar performances, except for the confusion rate in the functions reduced by scaling in frequency.” Journal of the Acoustical elevation of 22.5 degrees, where PCA method has significantly Society of America, vol. 106, no. 3, pp. 1480–1492, 1999. [11] D. N. Zotkin, J. Hwang, R. Duraiswaini, and L. S. Davis, “Hrtf smaller confusion rate than the SPCA method. However, personalization using anthropometric measurements,” in 2003 IEEE the proposed SPCA method predicts the HRTFs outside the Workshop on Applications of Signal Processing To Audio and Acoustics, training data well, which indicates that the SPCA method can New Paltz, NY, USA, 2003, pp. 157–160. [12] C. Jin, P. Leong, J. Leung, A. Corderoy, and S. Carlile, “Enabling in- predict the HRTFs in unmeasured spatial directions. dividualized virtual auditory space using morphological measurements,” in Proceedings of the First IEEE Pacific-Rim Conference on Multimedia, Sydney, Australia, 2000, pp. 235–238. VI. CONCLUSION [13] H. Hu, L. Zhou, H. Ma, and Z. Wu, “Hrtf personalization based on artificial neural network in individual virtual auditory space,” Applied This paper proposed the individual HRTFs modeling method Acoustics, vol. 69, no. 2, pp. 163–172, 2008. using DNN based on SPCA. By modeling the SPCs, the SPCA [14] C. Guezenoc and R. Seguier, “Hrtf individualization: A survey,” in Audio weights, the H and the ITDs respectively, we reconstruct the Engineering Society Convention 145. Audio Engineering Society, 2018. av [15] K. J. Fink and L. Ray, “Tuning principal component weights to indi- individual HRIRs in arbitrary spatial directions. The objective vidualize hrtfs,” in IEEE International Conference on Acoustics, Speech and subjective experiments are performed to evaluate the and Signal Processing, Kyoto, Japan, 2012, pp. 389–392. individual HRTFs generated by the proposed method, the [16] Y. Luo, D. N. Zotkin, and R. Duraiswami, “Virtual autoencoder based recommendation system for individualizing head-related transfer func- PCA method, and the generic method. Both the objective and tions,” in IEEE Workshop on Applications of Signal Processing to Audio subjective experiments’ results show that the proposed SPCA and Acoustics, New York, USA, 2013, pp. 1–4. method and PCA method are superior than the generic method. [17] B. Xie, “Recovery of individual head-related transfer functions from a small set of measurements.” Journal of the Acoustical Society of The spectral distortion of the SPCA method is significantly America, vol. 132, no. 1, pp. 282–294, 2012. smaller than PCA method in high frequencies but significantly [18] G. D. Romigh, D. S. Brungart, R. M. Stern, and B. D. Simpson, larger in low frequencies, except for few frequencies with no “Efficient real spherical harmonic representation of head-related transfer functions,” IEEE Journal of Selected Topics in Signal Processing, vol. 9, significant difference. The evaluation using the localization no. 5, pp. 921–930, 2015. model shows that the PCA method is better than the SPCA [19] J. O. Smith, “Techniques for digital filtering design and system identi- method. The subjective experiments show the SPCA method fication with the violin,” Ph.D. dissertation, CCRMA, Stanford, 1983. and the PCA method have similar performances, except for [20] I. Jolliffe, “Principal component analysis,” Technometrics, vol. 45, no. 3, p. 276, 2003. the confusion rate in the elevation of 22.5 degrees, where PCA [21] D. J. Kistler and F. L. Wightman, “A model of head-related transfer method has a better performance than the SPCA method. Nev- functions based on principal components analysis and minimum-phase ertheless, the results indicate that our proposed SPCA method reconstruction,” Journal of the Acoustical Society of America, vol. 91, no. 3, p. 1637, 1992. could predict the HRTFs in arbitrary spatial directions well. [22] M. Zhang, R. A. Kennedy, T. D. Abhayapala, and W. Zhang, “Statistical For the future work, we will use more HRTF data to implement method to identify key anthropometric parameters in hrtf individual- the individual HRTF models for better performance. ization,” in The Workshop on Hands-Free Speech Communication & Microphone Arrays, Edinburgh, UK, 2011, pp. 213–218. [23] X. Zeng, S. Wang, and L. Gao, “A hybrid algorithm for selecting head-related transfer function based on similarity of anthropometric REFERENCES structures,” Journal of Sound & Vibration, vol. 329, no. 19, pp. 4093– 4106, 2010. [1] J. Blauert, Spatial hearing: the psychophysics of human sound localiza- [24] R. Nicol, V. Lemaire, A. Bondu, and S. Busson, “Looking for a relevant tion. MIT press, 1997. similarity criterion for hrtf clustering: a comparative study,” in Audio [2] E. Wenzel, F. Wightman, and D. Kistler, “Localization with non- Engineering Society Convention 120. Paris, France: Audio Engineering individualized virtual acoustic display cues,” in Sigchi Conference on Society, 2006, p. 6653. Human Factors in Computing Systems, 1991, pp. 351–359. [25] R. B. Palm, “Prediction as a candidate for learning deep hierarchical [3] E. M. Wenzel, M. Arruda, D. J. Kistler, and F. L. Wightman, “Lo- models of data,” Master’s thesis, 2012. calization using nonindividualized head-related transfer functions,” The Journal of the Acoustical Society of America, vol. 94, no. 1, pp. 111– [26] V. R. Algazi, C. Avendano, and R. O. Duda, “Estimation of a spherical- 123, 1993. head model from anthropometry,” Journal of the Audio Engineering [4] V. R. Algazi, R. O. Duda, D. M. Thompson, and C. Avendano, “The Society, vol. 49, no. 6, pp. 472–479, 2001. cipic hrtf database,” in 2001 IEEE Workshop on the Applications of [27] C. I. Cheng and G. H. Wakefield, “Spatial frequency response surfaces Signal Processing to Audio and Acoustics. New Platz, NY, USA,: (sfr’s): An alternative visualization and interpolation technique for head- IEEE, 2001, pp. 99–102. related transfer functions (hrtf’s),” in Audio Engineering Society Con- [5] T. Qu, Z. Xiao, M. Gong, Y. Huang, X. Li, and X. Wu, “Distance- ference: 16th International Conference: Spatial Sound Reproduction. dependent head-related transfer functions measured with high spatial Audio Engineering Society, 1999. resolution using a spark gap,” IEEE Transactions on Audio, Speech, [28] S. Prepelia, ˘ M. Geronazzo, F. Avanzini, and L. Savioja, “Influence of and Language Processing, vol. 17, no. 6, pp. 1124–1132, 2009. voxelization on finite difference time domain simulations of head-related [6] W. Kreuzer, P. Majdak, and Z. Chen, “Fast multipole boundary element transfer functions,” The Journal of the Acoustical Society of America, method to calculate head-related transfer functions for a wide frequency vol. 139, no. 5, pp. 2489–2504, 2016. range,” Journal of the Acoustical Society of America, vol. 126, no. 3, [29] M. S. A. Zilany, I. C. Bruce, and L. H. Carney, “Updated parameters pp. 1280–90, 2009. and expanded simulation options for a model of the auditory periphery,” IEEE/ACM TRANSACTIONS ON AUDIO, SPEECH AND LANGUAGE PROCESSING, VOL. 28, NO. 1, 2020 13 The Journal of the Acoustical Society of America, vol. 135, no. 1, pp. 283–286, 2014. [30] P. Majdak, R. Baumgartner, and B. Laback, “Acoustic and non-acoustic factors in modeling listener-specific performance of sagittal-plane sound localization,” Frontiers in psychology, vol. 5, p. 319, 2014. [31] P. L. Søndergaard and P. Majdak, “The auditory modeling toolbox,” in The technology of binaural listening. Springer, 2013, pp. 33–56. [32] J. C. Middlebrooks, “Virtual localization improved by scaling nonindi- vidualized external-ear transfer functions in frequency,” The Journal of the Acoustical Society of America, vol. 106, no. 3, pp. 1493–1510, 1999. [33] B. Masiero and J. Fels, “Perceptually robust headphone equalization for binaural reproduction,” in Audio Engineering Society Convention 130. Audio Engineering Society, 2011.

Journal

Electrical Engineering and Systems SciencearXiv (Cornell University)

Published: Oct 21, 2019

There are no references for this article.