supdeq_compareTimeAlignmentMethods

% SUpDEq - Spatial Upsampling by Directional Equalization

 script supdeq_compareTimeAlignmentMethods
 
 In this script, various preprocessing / time-alignment methods [1,2,3,5] are implemented
 and compared to the SUpDEq method [4]. It is the basis for the journal publication [6].
 In brief, the sparse HRTF set gets time-aligned with the respective method, interpolated
 to a dense HRTF set in SH domain, and then HRTFs get reconstructed with the respective
 inverse alignment. The methods are described and compared in more detail in [6].

 Dependencies: SUpDEq toolbox with SOFiA and AKtools
               https://github.com/AudioGroupCologne/SUpDEq

 References:
 Onset-Based Time-Alignment (OBTA)
 [1] M. J. Evans, J. A. S. Angus, and A. I. Tew, "Analyzing head-related transfer function measurements 
 using surface spherical harmonics," 
 J. Acoust. Soc. Am., vol. 104, no. 4, pp. 2400-2411, 1998.
 [2] F. Brinkmann and S. Weinzierl, "Comparison of head-related transfer functions pre-processing 
 techniques for spherical harmonics decomposition," 
 in Proceedings of the AES Conference on Audio for Virtual and Augmented Reality, 2018, pp. 1-10.

 Frequency-Dependenten Time-Alignment (FDTA)
 [3] M. Zaunschirm, C. Schoerkhuber, and R. Hoeldrich, 
 "Binaural rendering of Ambisonic signals by HRIR time alignment and a diffuseness constraint," 
 J. Acoust. Soc. Am., vol. 143, no. 6, pp. 3616-3627, 2018.

 Spatial Upsampling by Directional Equalization (SUpDEq)
 [4] 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.

 Phase-Correction (PC)
 [5] Z. Ben-Hur, D. Lou Alon, R. Mehra, and B. Rafaely, 
 "Efficient Representation and Sparse Sampling of Head-Related Transfer Functions 
 Using Phase-Correction Based on Ear Alignment," 
 IEEE/ACM Trans. Audio, Speech, Lang. Process., vol. 27, no. 12, pp. 2249-2262, 2019.

 [6] J.M. Arend, F. Brinkmann, and C. Pörschmann, 
 "Assessing Spherical Harmonics Interpolation of Time-Aligned Head-Related
 Transfer Functions", Journal of the Audio Engineering Society (JAES)

 (C) 2020 by JMA, Johannes M. Arend
             Technische Hochschule Köln
             University of Applied Sciences
             Institute of Communications Engineering
             Department of Acoustics and Audio Signal Processing

%

0001 %% SUpDEq - Spatial Upsampling by Directional Equalization
0002 %
0003 % script supdeq_compareTimeAlignmentMethods
0004 %
0005 % In this script, various preprocessing / time-alignment methods [1,2,3,5] are implemented
0006 % and compared to the SUpDEq method [4]. It is the basis for the journal publication [6].
0007 % In brief, the sparse HRTF set gets time-aligned with the respective method, interpolated
0008 % to a dense HRTF set in SH domain, and then HRTFs get reconstructed with the respective
0009 % inverse alignment. The methods are described and compared in more detail in [6].
0010 %
0011 % Dependencies: SUpDEq toolbox with SOFiA and AKtools
0012 %               https://github.com/AudioGroupCologne/SUpDEq
0013 %
0014 % References:
0015 % Onset-Based Time-Alignment (OBTA)
0016 % [1] M. J. Evans, J. A. S. Angus, and A. I. Tew, "Analyzing head-related transfer function measurements
0017 % using surface spherical harmonics,"
0018 % J. Acoust. Soc. Am., vol. 104, no. 4, pp. 2400-2411, 1998.
0019 % [2] F. Brinkmann and S. Weinzierl, "Comparison of head-related transfer functions pre-processing
0020 % techniques for spherical harmonics decomposition,"
0021 % in Proceedings of the AES Conference on Audio for Virtual and Augmented Reality, 2018, pp. 1-10.
0022 %
0023 % Frequency-Dependenten Time-Alignment (FDTA)
0024 % [3] M. Zaunschirm, C. Schoerkhuber, and R. Hoeldrich,
0025 % "Binaural rendering of Ambisonic signals by HRIR time alignment and a diffuseness constraint,"
0026 % J. Acoust. Soc. Am., vol. 143, no. 6, pp. 3616-3627, 2018.
0027 %
0028 % Spatial Upsampling by Directional Equalization (SUpDEq)
0029 % [4] C. Pörschmann, J. M. Arend, and F. Brinkmann,
0030 % "Directional Equalization of Sparse Head-Related Transfer Function Sets for Spatial Upsampling,"
0031 % IEEE/ACM Trans. Audio, Speech, Lang. Process., vol. 27, no. 6, pp. 1060-1071, 2019.
0032 %
0033 % Phase-Correction (PC)
0034 % [5] Z. Ben-Hur, D. Lou Alon, R. Mehra, and B. Rafaely,
0035 % "Efficient Representation and Sparse Sampling of Head-Related Transfer Functions
0036 % Using Phase-Correction Based on Ear Alignment,"
0037 % IEEE/ACM Trans. Audio, Speech, Lang. Process., vol. 27, no. 12, pp. 2249-2262, 2019.
0038 %
0039 % [6] J.M. Arend, F. Brinkmann, and C. Pörschmann,
0040 % "Assessing Spherical Harmonics Interpolation of Time-Aligned Head-Related
0041 % Transfer Functions", Journal of the Audio Engineering Society (JAES)
0042 %
0043 % (C) 2020 by JMA, Johannes M. Arend
0044 %             Technische Hochschule Köln
0045 %             University of Applied Sciences
0046 %             Institute of Communications Engineering
0047 %             Department of Acoustics and Audio Signal Processing
0048 %
0049 %%
0050 clear all;
0051 close all;
0052 
0053 %Chose SH order of sparse HRTF dataset (Lebedev sampling scheme)
0054 sparseOrder = 2;
0055 %Get sparse HRTF set
0056 [sparseHRTFdataset,sparseHRIRdataset] = supdeq_getSparseDataset(supdeq_lebedev([],sparseOrder),sparseOrder,44);
0057 %Get fs
0058 fs = sparseHRTFdataset.f(end)*2;
0059 
0060 %Define dense sampling grid for upsampling
0061 denseSamplingGrid = supdeq_lebedev(2702);
0062 %Load reference dataset from HD
0063 referenceHRTFdataset = importdata('HRIRs_sfd_N44.mat');
0064 %Get Ndense, which in this case is the N = 44 (order of
0065 %supdeq_lebedev(2702)
0066 Ndense = referenceHRTFdataset.N;
0067 
0068 %Calculate optimal radius for KU100 according to Algazi et al. (2001)
0069 rKU100 = supdeq_optRadius(0.155,0.25,0.20,'Algazi'); 
0070 %Define c
0071 c = 343;
0072 
0073 %% Perform Onset-Based Time-Alignment (OBTA)
0074 %Implementation according to [2]
0075 
0076 %Estimate TOA on low-passed and 10 times upsampled HRIRs according to [2]
0077 toa_L = AKonsetDetect(sparseHRIRdataset.HRIR_L, 10, -20, 'rel', [3e3 fs]);
0078 toa_R = AKonsetDetect(sparseHRIRdataset.HRIR_R, 10, -20, 'rel', [3e3 fs]);
0079 
0080 %Transform TOAs to SH domain
0081 toa_L_nm = AKsht(toa_L, false, sparseHRIRdataset.samplingGrid, sparseHRIRdataset.Nmax, 'complex', fs, false, 'real');
0082 toa_R_nm = AKsht(toa_R, false, sparseHRIRdataset.samplingGrid, sparseHRIRdataset.Nmax, 'complex', fs, false, 'real');
0083 
0084 %Time-Align HRIRs (remove TOAs) with fractional delay
0085 HRIR_L_OBTA = AKfractionalDelayCyclic(sparseHRIRdataset.HRIR_L, -toa_L);
0086 HRIR_R_OBTA = AKfractionalDelayCyclic(sparseHRIRdataset.HRIR_R, -toa_R);
0087 
0088 %Transform to SH domain with Nmax according to sparse dataset
0089 HRTF_L_OBTA = AKboth2singleSidedSpectrum (fft(HRIR_L_OBTA));
0090 HRTF_R_OBTA = AKboth2singleSidedSpectrum (fft(HRIR_R_OBTA));
0091 obtaHRTFdataset = supdeq_hrtf2sfd(HRTF_L_OBTA.',HRTF_R_OBTA.',sparseHRTFdataset.Nmax,sparseHRTFdataset.samplingGrid,fs,'ak');
0092 
0093 %Upsample to dense grid by inverse SH transform
0094 [HRIR_L_Dense_OBTA, HRIR_R_Dense_OBTA] = supdeq_getArbHRIR(obtaHRTFdataset,denseSamplingGrid,'DEG',2,'ak');
0095 
0096 %Get interpolated TOAs by inverse SH transform
0097 toa_L_Dense = AKisht(toa_L_nm, false, denseSamplingGrid(:,1:2), 'complex', false, false, 'real');
0098 toa_R_Dense = AKisht(toa_R_nm, false, denseSamplingGrid(:,1:2), 'complex', false, false, 'real');
0099 
0100 %Re-Insert TOAs with fractional delays
0101 HRIR_L_Dense_OBTA = AKfractionalDelayCyclic(HRIR_L_Dense_OBTA, toa_L_Dense);
0102 HRIR_R_Dense_OBTA = AKfractionalDelayCyclic(HRIR_R_Dense_OBTA, toa_R_Dense);
0103 
0104 %Transform dense HRIR dataset to SH domain (apply zero padding according to
0105 % reference) at high order Ndense
0106 HRTF_L_Dense_OBTA = AKboth2singleSidedSpectrum (fft(HRIR_L_Dense_OBTA,size(HRIR_L_Dense_OBTA,1)*referenceHRTFdataset.FFToversize));
0107 HRTF_R_Dense_OBTA = AKboth2singleSidedSpectrum (fft(HRIR_R_Dense_OBTA,size(HRIR_R_Dense_OBTA,1)*referenceHRTFdataset.FFToversize));
0108 denseHRTFdataset_obta = supdeq_hrtf2sfd(HRTF_L_Dense_OBTA.',HRTF_R_Dense_OBTA.',Ndense,denseSamplingGrid,fs,'ak');
0109 denseHRTFdataset_obta.FFToversize = referenceHRTFdataset.FFToversize;
0110 
0111 %Clear everything besides required variables and final dense HRTF set in SH domain
0112 clearvars -EXCEPT denseHRTFdataset_obta c denseSamplingGrid fs Ndense referenceHRTFdataset rKU100 sparseHRIRdataset sparseHRTFdataset sparseOrder
0113 
0114 %% Perform Frequency-Dependente Time-Alignment (FDTA)
0115 
0116 %Transform sampling grid as required with azimuth between -pi to pi and
0117 %elevation -pi/2 to pi/2
0118 nPoints = size(sparseHRTFdataset.samplingGrid,1);
0119 az = sparseHRTFdataset.samplingGrid(:,1);
0120 az(az>180) = (360-az(az>180))*-1;
0121 el = 90-sparseHRTFdataset.samplingGrid(:,2);
0122 azRad = az * pi / 180;
0123 elRad = el * pi / 180;
0124 
0125 % Calculate the relative TOAs for each grid point according to Eq. 16 in [3]
0126 tl = cos(elRad).*sin(azRad)*rKU100*(c^-1); 
0127 tr = -tl;
0128 
0129 %Get normalized cut on frequency at f = 1500Hz
0130 fVec =  sparseHRTFdataset.f;
0131 fCut = 1500;
0132 [~,fCutBin] = min(abs(fVec-fCut));
0133 
0134 %Design allpass filter for each grid point according to Eq.26 in [3]
0135 wc = 2*pi*fVec(fCutBin);
0136 apl = zeros(nPoints,length(fVec));
0137 apl(:,1:fCutBin-1) = 1;
0138 apr = zeros(nPoints,length(fVec));
0139 apr(:,1:fCutBin-1) = 1;
0140 for id = 1:nPoints
0141     for kk = fCutBin:size(apl,2)
0142         apl(id,kk) = exp(-1j*(2*pi*fVec(kk)-wc)*tl(id));
0143         apr(id,kk) = exp(-1j*(2*pi*fVec(kk)-wc)*tr(id));
0144     end
0145 end
0146 
0147 %Apply allpass to HRTFs
0148 HRTF_L_FDTA = sparseHRTFdataset.HRTF_L .* apl;
0149 HRTF_R_TAZ = sparseHRTFdataset.HRTF_R .* apr;
0150 
0151 %Transform to SH domain with Nmax according to sparse dataset
0152 fdtaHRTFdataset = supdeq_hrtf2sfd(HRTF_L_FDTA,HRTF_R_TAZ,sparseHRTFdataset.Nmax,sparseHRTFdataset.samplingGrid,fs,'ak');
0153 
0154 %Upsample to dense grid by inverse SH transform
0155 [HRTF_L_Dense_FDTA, HRTF_R_Dense_FDTA] = supdeq_getArbHRTF(fdtaHRTFdataset,denseSamplingGrid,'DEG',2,'ak');
0156 
0157 %Add relative TOAs for dense sampling grid to interpolated HRTFs
0158 nPointsDense = size(denseSamplingGrid,1);
0159 azDense = denseSamplingGrid(:,1);
0160 azDense(azDense>180) = (360-azDense(azDense>180))*-1;
0161 elDense = 90-denseSamplingGrid(:,2);
0162 azDenseRad = azDense * pi / 180;
0163 elDenseRad = elDense * pi / 180;
0164 
0165 %Design allpass filter for each grid point of dense grid according to Eq.26 in [3]
0166 tlDense = cos(elDenseRad).*sin(azDenseRad)*rKU100*(c^-1);
0167 trDense = -tlDense;
0168 
0169 aplDense = zeros(nPointsDense,length(fVec));
0170 aplDense(:,1:fCutBin-1) = 1;
0171 aprDense = zeros(nPointsDense,length(fVec));
0172 aprDense(:,1:fCutBin-1) = 1;
0173 for id = 1:nPointsDense
0174     for kk = fCutBin:size(apl,2)
0175         %No minus before 1j to get inverse phase term. Could also be achieved by division instead of multiplication
0176         aplDense(id,kk) = exp(1j*(2*pi*fVec(kk)-wc)*tlDense(id)); 
0177         aprDense(id,kk) = exp(1j*(2*pi*fVec(kk)-wc)*trDense(id));
0178     end
0179 end
0180 
0181 %Add relative TOAs of dense sampling grid to interpolated HRTFs
0182 HRTF_L_Dense_FDTA = HRTF_L_Dense_FDTA .* aplDense;
0183 HRTF_R_Dense_FDTA = HRTF_R_Dense_FDTA .* aprDense;
0184 
0185 %Transform dense HRTF dataset to SH domain
0186 denseHRTFdataset_fdta = supdeq_hrtf2sfd(HRTF_L_Dense_FDTA,HRTF_R_Dense_FDTA,Ndense,denseSamplingGrid,fs,'ak');
0187 denseHRTFdataset_fdta.FFToversize = referenceHRTFdataset.FFToversize;
0188 
0189 %Clear everything besides required variables and final dense HRTF set in SH domain
0190 clearvars -EXCEPT denseHRTFdataset_fdta denseHRTFdataset_obta c denseSamplingGrid fs Ndense referenceHRTFdataset rKU100 sparseHRIRdataset sparseHRTFdataset sparseOrder
0191 
0192 %% Perform Spatial Upsampling by Directional Equalization (SUpDEq)
0193 
0194 %Get eqDataset
0195 eqDataset = supdeq_getEqDataset(Ndense,2*rKU100,length(referenceHRTFdataset.f)*2-2,fs);
0196 
0197 %Limit eqDataset (Small improvement explained in [6] leading to better
0198 %results below the spatial aliasing frequency fA)
0199 eqDataset = supdeq_limitEqDataset(eqDataset,sparseHRTFdataset.Nmax,eqDataset.radius);
0200 
0201 %Equalization
0202 [eqHRTFdataset, HRTF_equalized_L, HRTF_equalized_R] = supdeq_eq(sparseHRTFdataset,eqDataset,sparseHRTFdataset.Nmax,sparseHRTFdataset.samplingGrid);
0203 
0204 %De-Equalization
0205 [~,~,denseHRTFdataset_supdeq] = supdeq_deq(eqHRTFdataset, eqDataset, Ndense, denseSamplingGrid);
0206 
0207 %Clear everything besides required variables and final dense HRTF set in SH domain
0208 clearvars -EXCEPT denseHRTFdataset_supdeq denseHRTFdataset_fdta denseHRTFdataset_obta c denseSamplingGrid fs Ndense referenceHRTFdataset rKU100 sparseHRIRdataset sparseHRTFdataset sparseOrder
0209 
0210 %% Perform Phase-Correction (PC)
0211 
0212 %Get wavenumber k
0213 k  = 2*pi*sparseHRTFdataset.f/c; k = k.';
0214 
0215 %Transform sampling grid to radiant
0216 sg = sparseHRTFdataset.samplingGrid * pi / 180;
0217 
0218 %Get phase correction term for left/right ear according to Eq. 13 in [5]
0219 cosThetaL = cos(sg(:,2)')*cos(pi/2) + sin(sg(:,2)')*sin(pi/2) .* cos(sg(:,1)'-pi/2); %Left ear with -pi/2
0220 phaseCorL = exp(-1j*rKU100 * k .* cosThetaL);
0221 cosThetaR = cos(sg(:,2)')*cos(pi/2) + sin(sg(:,2)')*sin(pi/2) .* cos(sg(:,1)'+pi/2); %Right ear with +pi/2
0222 phaseCorR = exp(-1j*rKU100 * k .* cosThetaR);
0223 
0224 %Apply to HRTFs
0225 HRTF_L_PC = sparseHRTFdataset.HRTF_L .* phaseCorL.'; %Just flip, no conjugate complex
0226 HRTF_R_PC = sparseHRTFdataset.HRTF_R .* phaseCorR.';
0227 
0228 %Transform to SH domain with Nmax according to sparse dataset
0229 pcHRTFdataset = supdeq_hrtf2sfd(HRTF_L_PC,HRTF_R_PC,sparseHRTFdataset.Nmax,sparseHRTFdataset.samplingGrid,fs,'ak');
0230 
0231 %Upsample to dense grid by inverse SH transform
0232 [HRTF_L_Dense_PC, HRTF_R_Dense_PC] = supdeq_getArbHRTF(pcHRTFdataset,denseSamplingGrid,'DEG',2,'ak');
0233 
0234 %Transform dense sampling grid to radiant
0235 denseSamplingGridRad = denseSamplingGrid * pi / 180;
0236 
0237 %Add relative TOAs of dense sampling grid to HRTFs (inverse phase
0238 %correction)
0239 cosThetaDenseL = cos(denseSamplingGridRad(:,2)')*cos(pi/2) + sin(denseSamplingGridRad(:,2)')*sin(pi/2) .* cos(denseSamplingGridRad(:,1)'-pi/2);
0240 %No minus before 1j to get inverse phase term. Could also be achieved by division instead of multiplication
0241 phaseCorDenseL = exp(1j*rKU100 * k .* cosThetaDenseL);
0242 cosThetaDenseR = cos(denseSamplingGridRad(:,2)')*cos(pi/2) + sin(denseSamplingGridRad(:,2)')*sin(pi/2) .* cos(denseSamplingGridRad(:,1)'+pi/2);
0243 %No minus before 1j to get inverse phase term. Could also be achieved by division instead of multiplication
0244 phaseCorDenseR = exp(1j*rKU100 * k .* cosThetaDenseR);
0245 HRTF_L_Dense_PC = HRTF_L_Dense_PC .* phaseCorDenseL.'; %Just flip, no conjugate complex
0246 HRTF_R_Dense_PC = HRTF_R_Dense_PC .* phaseCorDenseR.';
0247 
0248 %Transform dense HRTF dataset to SH domain
0249 denseHRTFdataset_pc = supdeq_hrtf2sfd(HRTF_L_Dense_PC,HRTF_R_Dense_PC,Ndense,denseSamplingGrid,fs,'ak');
0250 denseHRTFdataset_pc.FFToversize = referenceHRTFdataset.FFToversize;
0251 
0252 %Clear everything besides required variables and final dense HRTF set in SH domain
0253 clearvars -EXCEPT denseHRTFdataset_pc denseHRTFdataset_supdeq denseHRTFdataset_fdta denseHRTFdataset_obta c denseSamplingGrid fs Ndense referenceHRTFdataset rKU100 sparseHRIRdataset sparseHRTFdataset sparseOrder
0254 
0255 %%
0256 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
0257 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
0258 %
0259 %Some plots for physical evaluation. See more detailed results in [6] and
0260 %in supplementary material of [6].
0261 %
0262 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
0263 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
0264 
0265 %% Spectral difference over the entire sphere (2702 Lebedev grid)
0266 
0267 supdeq_calcSpectralDifference(denseHRTFdataset_obta,referenceHRTFdataset,[],[],[],true);
0268 titleString = ['OBTA - N=',num2str(sparseOrder)];
0269 title(titleString);
0270 ylim([0 6]);
0271 %saveas(gcf,['SpecDiff_OBTA_N_',num2str(sparseOrder),'.jpg']);
0272 
0273 supdeq_calcSpectralDifference(denseHRTFdataset_fdta,referenceHRTFdataset,[],[],[],true);
0274 titleString = ['FDTA - N=',num2str(sparseOrder)];
0275 title(titleString);
0276 ylim([0 6]);
0277 %saveas(gcf,['SpecDiff_FDTA_N_',num2str(sparseOrder),'.jpg']);
0278 
0279 supdeq_calcSpectralDifference(denseHRTFdataset_supdeq,referenceHRTFdataset,[],[],[],true);
0280 titleString = ['SUpDEq - N=',num2str(sparseOrder)];
0281 title(titleString);
0282 ylim([0 6]);
0283 %saveas(gcf,['SpecDiff_SUpDEq_N_',num2str(sparseOrder),'.jpg']);
0284 
0285 supdeq_calcSpectralDifference(denseHRTFdataset_pc,referenceHRTFdataset,[],[],[],true);
0286 titleString = ['PC - N=',num2str(sparseOrder)];
0287 title(titleString);
0288 ylim([0 6]);
0289 %saveas(gcf,['SpecDiff_PC_N_',num2str(sparseOrder),'.jpg']);
0290 
0291 %Clear everything besides required variables and final dense HRTF set in SH domain
0292 clearvars -EXCEPT denseHRTFdataset_pc denseHRTFdataset_supdeq denseHRTFdataset_fdta denseHRTFdataset_obta c denseSamplingGrid fs Ndense referenceHRTFdataset rKU100 sparseHRIRdataset sparseHRTFdataset sparseOrder
0293 
0294 %% Spectral difference in region around the contralateral (only horizontal plane)
0295 
0296 testSamplingGrid = [(250:290).',ones(41,1)*90];
0297 
0298 supdeq_calcSpectralDifference(denseHRTFdataset_obta,referenceHRTFdataset,testSamplingGrid,[],[],true);
0299 titleString = ['OBTA - Lateral - N=',num2str(sparseOrder)];
0300 title(titleString);
0301 ylim([0 15]);
0302 %saveas(gcf,['SpecDiff_Lateral_OBTA_N_',num2str(sparseOrder),'.jpg']);
0303 
0304 supdeq_calcSpectralDifference(denseHRTFdataset_fdta,referenceHRTFdataset,testSamplingGrid,[],[],true);
0305 titleString = ['FDTA - Lateral - N=',num2str(sparseOrder)];
0306 title(titleString);
0307 ylim([0 15]);
0308 %saveas(gcf,['SpecDiff_Lateral_FDTA_N_',num2str(sparseOrder),'.jpg']);
0309 
0310 supdeq_calcSpectralDifference(denseHRTFdataset_supdeq,referenceHRTFdataset,testSamplingGrid,[],[],true);
0311 titleString = ['SUpDEq - Lateral - N=',num2str(sparseOrder)];
0312 title(titleString);
0313 ylim([0 15]);
0314 %saveas(gcf,['SpecDiff_Lateral_SUpDEq_N_',num2str(sparseOrder),'.jpg']);
0315 
0316 supdeq_calcSpectralDifference(denseHRTFdataset_pc,referenceHRTFdataset,testSamplingGrid,[],[],true);
0317 titleString = ['PC - Lateral - N=',num2str(sparseOrder)];
0318 title(titleString);
0319 ylim([0 15]);
0320 %saveas(gcf,['SpecDiff_Lateral_PC_N_',num2str(sparseOrder),'.jpg']);
0321 
0322 %Clear everything besides required variables and final dense HRTF set in SH domain
0323 clearvars -EXCEPT denseHRTFdataset_pc denseHRTFdataset_supdeq denseHRTFdataset_fdta denseHRTFdataset_obta c denseSamplingGrid fs Ndense referenceHRTFdataset rKU100 sparseHRIRdataset sparseHRTFdataset sparseOrder
0324 
0325 %% Get some HRIRs (lateral)
0326 
0327 %Define direction
0328 az = 90;
0329 el = 90;
0330 
0331 %Get HRIRs
0332 hrir_ref = supdeq_getArbHRIR(referenceHRTFdataset,[az,el]);
0333 hrir_obta = supdeq_getArbHRIR(denseHRTFdataset_obta,[az,el]);
0334 hrir_fdta = supdeq_getArbHRIR(denseHRTFdataset_fdta,[az,el]);
0335 hrir_supdeq = supdeq_getArbHRIR(denseHRTFdataset_supdeq,[az,el]);
0336 hrir_pc = supdeq_getArbHRIR(denseHRTFdataset_pc,[az,el]);
0337 
0338 %Plot HRIRs in comparison to reference
0339 supdeq_plotIR(hrir_ref,hrir_obta,[],[],32)
0340 supdeq_plotIR(hrir_ref,hrir_fdta,[],[],32)
0341 supdeq_plotIR(hrir_ref,hrir_supdeq,[],[],32)
0342 supdeq_plotIR(hrir_ref,hrir_pc,[],[],32)
0343 
0344 %Clear everything besides required variables and final dense HRTF set in SH domain
0345 clearvars -EXCEPT denseHRTFdataset_pc denseHRTFdataset_supdeq denseHRTFdataset_fdta denseHRTFdataset_obta c denseSamplingGrid fs Ndense referenceHRTFdataset rKU100 sparseHRIRdataset sparseHRTFdataset sparseOrder
0346 
0347 %% Get some HRIRs (contralateral)
0348 
0349 %Define direction
0350 az = 270;
0351 el = 90;
0352 
0353 %Get HRIRs
0354 hrir_ref = supdeq_getArbHRIR(referenceHRTFdataset,[az,el]);
0355 hrir_obta = supdeq_getArbHRIR(denseHRTFdataset_obta,[az,el]);
0356 hrir_fdta = supdeq_getArbHRIR(denseHRTFdataset_fdta,[az,el]);
0357 hrir_supdeq = supdeq_getArbHRIR(denseHRTFdataset_supdeq,[az,el]);
0358 hrir_pc = supdeq_getArbHRIR(denseHRTFdataset_pc,[az,el]);
0359 
0360 %Plot HRIRs in comparison to reference
0361 supdeq_plotIR(hrir_ref,hrir_obta,[],[],32)
0362 supdeq_plotIR(hrir_ref,hrir_fdta,[],[],32)
0363 supdeq_plotIR(hrir_ref,hrir_supdeq,[],[],32)
0364 supdeq_plotIR(hrir_ref,hrir_pc,[],[],32)
0365 
0366 %Clear everything besides required variables and final dense HRTF set in SH domain
0367 clearvars -EXCEPT denseHRTFdataset_pc denseHRTFdataset_supdeq denseHRTFdataset_fdta denseHRTFdataset_obta c denseSamplingGrid fs Ndense referenceHRTFdataset rKU100 sparseHRIRdataset sparseHRTFdataset sparseOrder
0368 
0369 %% Get circular grids for ITDs/ILDs
0370 
0371 %Get circ dataset with 1° in azimuth
0372 circGrid(1:360,1) = 0:359;
0373 circGrid(1:360,2) = 90;
0374 
0375 %Get circ grid HRTFs
0376 [circGrid_ref_L,circGrid_ref_R] = supdeq_getArbHRIR(referenceHRTFdataset,circGrid,'DEG',2,'ak'); 
0377 [circGrid_obta_L,circGrid_obta_R] = supdeq_getArbHRIR(denseHRTFdataset_obta,circGrid,'DEG',2,'ak'); 
0378 [circGrid_fdta_L,circGrid_fdta_R] = supdeq_getArbHRIR(denseHRTFdataset_fdta,circGrid,'DEG',2,'ak'); 
0379 [circGrid_supdeq_L,circGrid_supdeq_R] = supdeq_getArbHRIR(denseHRTFdataset_supdeq,circGrid,'DEG',2,'ak'); 
0380 [circGrid_pc_L,circGrid_pc_R] = supdeq_getArbHRIR(denseHRTFdataset_pc,circGrid,'DEG',2,'ak'); 
0381 
0382 %% Get ITDs (takes some time)
0383 
0384 th = -10; 
0385 for id = 1:size(circGrid,1)
0386     
0387     ITD_ref(id) = supdeq_calcITD([circGrid_ref_L(:,id),circGrid_ref_R(:,id)],fs,th);
0388     ITD_obta(id) = supdeq_calcITD([circGrid_obta_L(:,id),circGrid_obta_R(:,id)],fs,th);
0389     ITD_fdta(id) = supdeq_calcITD([circGrid_fdta_L(:,id),circGrid_fdta_R(:,id)],fs,th);
0390     ITD_supdeq(id) = supdeq_calcITD([circGrid_supdeq_L(:,id),circGrid_supdeq_R(:,id)],fs,th);
0391     ITD_pc(id) = supdeq_calcITD([circGrid_pc_L(:,id),circGrid_pc_R(:,id)],fs,th);
0392     
0393 end
0394 
0395 figure;
0396 lineWidth = 1.2;
0397 plot(ITD_ref,'LineWidth',lineWidth)
0398 hold on;
0399 plot(ITD_obta,'LineWidth',lineWidth); plot(ITD_fdta,'LineWidth',lineWidth); plot(ITD_supdeq,'LineWidth',lineWidth); plot(ITD_pc,'LineWidth',lineWidth);
0400 legend('Ref','OBTA','FDTA','SUpDEq','PA','Location','SouthEast');
0401 xlabel('Azimuth in degrees'); xlim([0 360]);
0402 ylabel('ITD in ms'); 
0403 
0404 %% Get broadband ILDs
0405 
0406 for id = 1:size(circGrid,1)
0407     
0408     ILD_ref(id) = supdeq_calcILD([circGrid_ref_L(:,id),circGrid_ref_R(:,id)]);
0409     ILD_obta(id) = supdeq_calcILD([circGrid_obta_L(:,id),circGrid_obta_R(:,id)]);
0410     ILD_fdta(id) = supdeq_calcILD([circGrid_fdta_L(:,id),circGrid_fdta_R(:,id)]);
0411     ILD_supdeq(id) = supdeq_calcILD([circGrid_supdeq_L(:,id),circGrid_supdeq_R(:,id)]);
0412     ILD_pc(id) = supdeq_calcILD([circGrid_pc_L(:,id),circGrid_pc_R(:,id)]);
0413     
0414 end
0415 
0416 figure;
0417 lineWidth = 1.2;
0418 plot(ILD_ref,'LineWidth',lineWidth)
0419 hold on;
0420 plot(ILD_obta,'LineWidth',lineWidth); plot(ILD_fdta,'LineWidth',lineWidth); plot(ILD_supdeq,'LineWidth',lineWidth); plot(ILD_pc,'LineWidth',lineWidth);
0421 legend('Ref','OBTA','FDTA','SUpDEq','PC','Location','NorthEast');
0422 xlabel('Azimuth in degrees'); xlim([0 360]);
0423 ylabel('ILD in dB');