% SUpDEq - Spatial Upsampling by Directional Equalization
script supdeq_demo
This script presents the SUpDEq processing and a basic technical
evaluation of the final (de-equalized) dense HRTF dataset
Reference:
C. Pörschmann*, J.M. Arend*, and F. Brinkmann, "Directional Equalization
of Sparse Head-Related Transfer Function Sets for Spatial Upsampling,"
IEEE/ACM Trans. Audio, Speech, Lang. Process., vol. 27, no. 6,
pp. 1060?1071, 2019, *These authors contributed equally to this work.
(C) 2018 by JMA, Johannes M. Arend
TH Köln - University of Applied Sciences
Institute of Communications Engineering
Department of Acoustics and Audio Signal Processing0001 %% SUpDEq - Spatial Upsampling by Directional Equalization 0002 % 0003 % script supdeq_demo 0004 % 0005 % This script presents the SUpDEq processing and a basic technical 0006 % evaluation of the final (de-equalized) dense HRTF dataset 0007 % 0008 % Reference: 0009 % C. Pörschmann*, J.M. Arend*, and F. Brinkmann, "Directional Equalization 0010 % of Sparse Head-Related Transfer Function Sets for Spatial Upsampling," 0011 % IEEE/ACM Trans. Audio, Speech, Lang. Process., vol. 27, no. 6, 0012 % pp. 1060?1071, 2019, *These authors contributed equally to this work. 0013 % 0014 % (C) 2018 by JMA, Johannes M. Arend 0015 % TH Köln - University of Applied Sciences 0016 % Institute of Communications Engineering 0017 % Department of Acoustics and Audio Signal Processing 0018 0019 %% (1) - Load reference HRTF dataset (stored as SH-coefficients in spherical harmonics domain) 0020 %See small script "get_HRIRs_sfd_N35.m" in materials folder to see how the 0021 %dataset was obtained based on the "HRIR_L2702.sofa" HRIR dataset of a 0022 %Neumann KU100 dummy head. The reference dataset is only needed for 0023 %comparison with the final de-equalized dataset... 0024 referenceHRTFdataset = importdata('HRIRs_sfd_N35.mat'); 0025 0026 %% (2) - Load sparse HRTF dataset (stored as FFT-bins in Fourier domain) 0027 %This is just an example. You could use any sparse HRTF dataset here! Also 0028 %have a look at the script "get_sparse_HRTF_set.m" in materials folder to 0029 %see how the sparse datasets are obtained. 0030 0031 %Here, we use a sparse HRTF daset with only 38 sample points (Lebedev grid) 0032 %The sampling grid as well as the highste possible transform oder N are 0033 %required for the equalization! It is saved in the struct here, but it 0034 %could also be passed separately to the function "subdeq_eq" 0035 load('sparseHRTFdataset_L38.mat'); 0036 0037 %Optionally, a sparse HRIR dataset in SOFA format can be loaded. With the 0038 %function supdeq_sofa2hrtf, the SOFA file will be transformed to 0039 %frequency domain and saved in a struct as used by the SUpDEq toolbox 0040 %(compare to the structs 'sparseHRTFdataset'). FFToversize and Nmax have to 0041 %be defined manually now, because SOFA does not support such meta data. 0042 %Furthermore, the sampling grid in SOFA does not support sampling weights, 0043 %which can be useful for spherical Fourier / SH transform. Thus, optionally 0044 %a sampling grid with weigths can be passed, which will than be stored in 0045 %the struct. Otherwise, the sampling grid stored in the SOFA file will be 0046 %written to the HRTF struct. 0047 SOFAinput = false; 0048 if SOFAinput 0049 %Load sparse HRIR dataset in SOFA format 0050 sparseHRIRdataset_SOFA = SOFAload('sparseHRIRdataset_L38.sofa'); 0051 %Transform to sparseHRTFdataset struct with pre-defined samplingGrid 0052 %(Lebedev grid with 38 nodes here), Nmax = 4, and FFToversize = 4. 0053 sparseHRTFdataset = supdeq_sofa2hrtf(sparseHRIRdataset_SOFA,4,supdeq_lebedev(38),4); 0054 end 0055 0056 %% (3) - Get equalization dataset (SH-coefficients) 0057 %The eqDataset describes the sound pressure distribution on a sphere 0058 %Use defaults: N = 35, earDistance = 0.165m, NFFT = 512, fs = 48000; 0059 eqDataset = supdeq_getEqDataset; 0060 0061 %% (4) - Perform equalization 0062 %Here, the sparse HRTF dataset is equalized with the eqDataset. The 0063 %equalized HRTF are transformed to the SH-domain again with the maximal 0064 %order N which is possible with the sparse sampling grid. 0065 %N and the sparse sampling grid are part of the sparseHRTFdataset struct 0066 sparseSamplingGrid = sparseHRTFdataset.samplingGrid; 0067 Nsparse = sparseHRTFdataset.Nmax; 0068 0069 eqHRTFdataset = supdeq_eq(sparseHRTFdataset,eqDataset,Nsparse,sparseSamplingGrid); 0070 0071 %% (5) - Perform de-equalization 0072 %Here, the sparse equalized HRTF dataset is de-equalized with the 0073 %deqDataset. This is done on a dense spatial sampling grid. The results is a 0074 %dense HRTF/HRIR dataset. In this case, deqDataset is the same as the 0075 %eqDataset... 0076 0077 %First, define dense spatial sampling grid. Here, we use the lebedev grid 0078 %with 2702 points again (same as the reference HRIR dataset). 0079 %The highest stable grid order here is N = 44. However, we use N = 35 for the 0080 %spherical Fourier transform of the de-equalized HRTFs. 0081 denseSamplingGrid = supdeq_lebedev(2702); 0082 Ndense = 35; 0083 0084 %Perform de-equalization. Apply head and tail window (8 and 32 samples 0085 %respectively) to de-equalized HRIRs/HRTFs. 0086 [denseHRTFdataset, denseHRIRdataset, denseHRTFdataset_sh] = supdeq_deq(eqHRTFdataset, eqDataset, Ndense, denseSamplingGrid,[8,32]); 0087 0088 %% (6) - Optional: Save as SOFA object 0089 %Use defaults: fs = 48000, earDistance = 0.165m, sourceDistance = 3.0m 0090 denseHRIRdataset_SOFA = supdeq_writeSOFAobj(denseHRIRdataset.HRIR_L,denseHRIRdataset.HRIR_R,denseSamplingGrid); 0091 0092 %% (7) - Optional: Plot HRIRs 0093 %Get HRIRs from reference dataset and de-equalized dense dataset. In this 0094 %example, we chose a lateral source, because most differences between the 0095 %reference and the de-equalized HRIRs occure at the contralateral ear. 0096 azPlot = 90; 0097 elPlot = 90; 0098 [hrir_ref(:,1),hrir_ref(:,2)] = supdeq_getArbHRIR(referenceHRTFdataset,[azPlot,elPlot]); 0099 [hrir_deq(:,1),hrir_deq(:,2)] = supdeq_getArbHRIR(denseHRTFdataset_sh,[azPlot,elPlot]); 0100 0101 %Plot left and right channel of reference and deq set in two plots with 0102 %32 x FFT oversampling 0103 supdeq_plotIR(hrir_ref(:,1),hrir_deq(:,1),[],[],32); 0104 supdeq_plotIR(hrir_ref(:,2),hrir_deq(:,2),[],[],32); 0105 0106 %% (8) - Optional: Listen to the results 0107 %Here, you can listen to a short drums sequence, spatialized with a 0108 %reference HRIR and a de-equalized HRIR 0109 0110 %Load drums test signal 0111 [testSignal,fsTestSignal] = audioread('drums.wav'); 0112 0113 %Define az and el for playback 0114 azPlayback = 30; 0115 elPlayback = 90; 0116 0117 %Convolve with reference HRIR and play back 0118 supdeq_listen(referenceHRTFdataset,testSignal,[azPlayback,elPlayback]); 0119 %Wait for length of test-signal 0120 pause(length(testSignal)/fsTestSignal+0.05); 0121 %Convolve with de-equalized HRIR and play back 0122 supdeq_listen(denseHRTFdataset_sh,testSignal,[azPlayback,elPlayback]); 0123 0124 %% (9) - Optional: Technical Evaluation 0125 0126 %Perform test of sagittal plane localization with default spatial sampling 0127 %grid 0128 spLocalization = supdeq_testSPlocalization(denseHRTFdataset_sh, referenceHRTFdataset); 0129 0130 %Perform test of horizontal plane localizatoin with default spatial 0131 %sampling grid 0132 hpLocalization = supdeq_testHPlocalization(denseHRTFdataset_sh, referenceHRTFdataset); 0133 0134 %Calculate spectral differences with default spatial sampling grid and 0135 %make simple plot 0136 specDifferences = supdeq_calcSpectralDifference(denseHRTFdataset_sh, referenceHRTFdataset,[],[],false,true); 0137 0138 %Calculate spectral differences in auditory bands with default spatial 0139 %sampling grid and make simple plot 0140 specDifferences_auditoryBands = supdeq_calcSpectralDifference(denseHRTFdataset_sh, referenceHRTFdataset,[],[],true,true); 0141 0142 %Calculate energy per order N and make simple plot 0143 energyPerOrder = supdeq_calcEnergyPerOrder(denseHRTFdataset_sh,true);