% SUpDEq - Spatial Upsampling by Directional Equalization
function [interpHRTFset, HRTF_interp_L, HRTF_interp_R] = supdeq_interpHRTF(HRTFset, ipSamplingGrid, ppMethod, ipMethod, mc, headRadius, tikhEps, limitMC, mcKnee, mcMinPhase, limFade)
This function performs HRTF interpolation of the input HRTFset according
to the passed interpolation sampling grid. Various interpolation and
time-alignment (pre-processing) methods can be applied. Furthermore,
post-interpolation magnitude correction for time-aligned interpolation
(MCA) can be applied to further improve the interpolation results.
Output:
interpHRTFset - Struct with the interpolated HRTF for the
left (HRTF_L) and right (HRTF_R) ear, absolute
frequency scale f, various interpolation parameters
according to the applied method, and FFToversize
HRTF_L/R_ip - Interpolated complex HRTFs (single-sided spectrum)
Input:
HRTFset - Struct with (sparse) HRTF dataset in frequency
domain (single-sided complex spectra). Struct needs to
provide HRTF_L/HRTF_R, the samplingGrid, and Nmax
ipSamplingGrid - Spatial sampling grid (Q x 2 or Q x 3 matrix),
defining the interpolation points, where the first
column holds the azimuth, the second the elevation
(both in degree), and optionally the third the
sampling weights.
Azimuth in degree
(0=front, 90=left, 180=back, 270=right)
(0 points to positive x-axis, 90 to positive y-axis)
Elevation in degree
(0=North Pole, 90=front, 180=South Pole)
(0 points to positive z-axis, 180 to negative z-axis)
ppMethod - String defining the pre/post-processing applied
to the HRTFs before/after interpolation. Besides
'MagPhase', all of the pre-processing methods are
time-alignment approaches.
'SUpDEq' - SUpDEq method [3][5]
'SUpDEq_Lim' - Modified SUpDeq method where the eqDataset is limited to frequencies above the spatial aliasing frequency fA [5]
'SUpDEq_AP' - Modified SUpDEq method applying rigid sphere allpass filter (only the phase components of the rigid sphere transfer functions)
'SUpDEq_Lim_AP' - Modified SUpDeq method where the eqDataset is limited as above (SUpDEq_Lim) and only the phase components are used (SUpDEq_AP)
'PC' - Phase-Correction (open sphere transfer functions) [4][5]
'OBTA' - Onset-Based Time-Alignment [1][2][5]
'MagPhase' - Split HRTF to magnitude and unwrapped phase - No time-alignment
'None' - Perform interpolation without any pre/post-processing
Default: 'SUpDEq'
ipMethod - String defining the ipMethod
'SH' - Spherical harmonics interpolation [1-5]
'NN' - Natural neighbor interpolation with Voronoi weights [6]
'Bary' - Barycentric interpolation [7]
Default: 'SH'
mc - Define maximum boost in dB if magnitude correction according to [8] should be applied to interpolated HRTFs in a further postprocessing step
Set to nan if magnitude correction should not be applied
Set to inf if no limiting should be applied, i.e., unlimited boost
Default: inf
headRadius - Head radius in m. Required for time-alignment methods SUpDEq [3][5] and PC [4][5]
Default: 0.0875
tikhEps - Define epsilon of Tikhonov regularization if regularization should be applied
Applying the Tikhonov regularization will always result in a least-square fit
solution for the SH transform. Only relevant if ipMethod is 'SH'.
Default: 0 (no Tikhonov regularization)
limitMC - Boolean to set magnitude correction filters to 0 dB below spatial aliasing frequency fA,
assuming that (SH) interpolation below fA is physically correct.
Default: true (correction filter operate only above fA)
Only applies if 'mc' is not nan
mcKnee - Knee in dB of the soft-limiting applied to the correction filters. Soft-limiting is applied according to maximum boost 'mc'
Only applies if 'mc' is not nan or inf
Default: 0 dB (no knee)
mcMinPhase - Boolean to design correction filters as minimum-phase filters. If set to false, correction filters are zero-phase filters
Default: true (minimum-phase filters) as this provides a better performance in terms of ITD errors in combination with classical SUpDEq processing
When applying other preprocessing methods, using zero-phase filters can be more appropriate.
limFade - String to define whether to apply a linear fade upwards from aliasing frequency fA to fAt = fA + 1/3 Oct, or downwards, i.e., from fA to fAt = fA - 1/3 Oct
'fadeUp' - upwards
'fadeDown' - downwards
Default: 'fadeDown'
Dependencies: SUpDEq toolbox, AKtools, SOFiA toolbox, SFS toolbox, TriangleRayIntersection
References: (Pre-Processing methods)
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.
Spatial Upsampling by Directional Equalization (SUpDEq)
[3] 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)
[4] 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.
Comparison of time-alignment methods & perceptual evaluation
[5] J. M. Arend, F. Brinkmann, and C. Pörschmann, “Assessing Spherical
Harmonics Interpolation of Time-Aligned Head-Related Transfer Functions,”
J. Audio Eng. Soc., vol. 69, no. 1/2, pp. 104–117, 2021.
Spherical harmonics interpolation
[1-5]
Natural-Neighbor Interpolation
[6] R. Sibson, “A brief description of natural neighbor interpolation,”
in Interpreting Multivariate Data, V. Barnett, Chichester, England: John Wiley & Sons, 1981, pp. 21–36.
Good read about NN interpolation: https://github.com/gwlucastrig/Tinfour/wiki/Introduction-to-Natural-Neighbor-Interpolation
Barycentric Coordinates / Interpolation
[7] P. Shirley and S. Marschner, Fundamentals of Computer Graphics, 3rd ed. Boca Raton, FL: Taylor & Francis, 2009.
Good read about Barycentric interpolation: https://www.scratchapixel.com/lessons/3d-basic-rendering/ray-tracing-rendering-a-triangle/barycentric-coordinates
Magnitude-Corrected and Time-Aligned Interpolation (MCA)
[8] J. M. Arend, C. Pörschmann, S. Weinzierl, F. Brinkmann,
"Magnitude-Corrected and Time-Aligned Interpolation of Head-Related Transfer Functions,"
(Manuscript submitted for publication).
(C) 2020-2022 by JMA, Johannes M. Arend
TU Berlin, Audio Communication Group
TH Köln, Institute of Communications Engineering0001 %% SUpDEq - Spatial Upsampling by Directional Equalization 0002 % 0003 % function [interpHRTFset, HRTF_interp_L, HRTF_interp_R] = supdeq_interpHRTF(HRTFset, ipSamplingGrid, ppMethod, ipMethod, mc, headRadius, tikhEps, limitMC, mcKnee, mcMinPhase, limFade) 0004 % 0005 % This function performs HRTF interpolation of the input HRTFset according 0006 % to the passed interpolation sampling grid. Various interpolation and 0007 % time-alignment (pre-processing) methods can be applied. Furthermore, 0008 % post-interpolation magnitude correction for time-aligned interpolation 0009 % (MCA) can be applied to further improve the interpolation results. 0010 % 0011 % Output: 0012 % interpHRTFset - Struct with the interpolated HRTF for the 0013 % left (HRTF_L) and right (HRTF_R) ear, absolute 0014 % frequency scale f, various interpolation parameters 0015 % according to the applied method, and FFToversize 0016 % HRTF_L/R_ip - Interpolated complex HRTFs (single-sided spectrum) 0017 % 0018 % Input: 0019 % HRTFset - Struct with (sparse) HRTF dataset in frequency 0020 % domain (single-sided complex spectra). Struct needs to 0021 % provide HRTF_L/HRTF_R, the samplingGrid, and Nmax 0022 % ipSamplingGrid - Spatial sampling grid (Q x 2 or Q x 3 matrix), 0023 % defining the interpolation points, where the first 0024 % column holds the azimuth, the second the elevation 0025 % (both in degree), and optionally the third the 0026 % sampling weights. 0027 % Azimuth in degree 0028 % (0=front, 90=left, 180=back, 270=right) 0029 % (0 points to positive x-axis, 90 to positive y-axis) 0030 % Elevation in degree 0031 % (0=North Pole, 90=front, 180=South Pole) 0032 % (0 points to positive z-axis, 180 to negative z-axis) 0033 % ppMethod - String defining the pre/post-processing applied 0034 % to the HRTFs before/after interpolation. Besides 0035 % 'MagPhase', all of the pre-processing methods are 0036 % time-alignment approaches. 0037 % 'SUpDEq' - SUpDEq method [3][5] 0038 % 'SUpDEq_Lim' - Modified SUpDeq method where the eqDataset is limited to frequencies above the spatial aliasing frequency fA [5] 0039 % 'SUpDEq_AP' - Modified SUpDEq method applying rigid sphere allpass filter (only the phase components of the rigid sphere transfer functions) 0040 % 'SUpDEq_Lim_AP' - Modified SUpDeq method where the eqDataset is limited as above (SUpDEq_Lim) and only the phase components are used (SUpDEq_AP) 0041 % 'PC' - Phase-Correction (open sphere transfer functions) [4][5] 0042 % 'OBTA' - Onset-Based Time-Alignment [1][2][5] 0043 % 'MagPhase' - Split HRTF to magnitude and unwrapped phase - No time-alignment 0044 % 'None' - Perform interpolation without any pre/post-processing 0045 % Default: 'SUpDEq' 0046 % ipMethod - String defining the ipMethod 0047 % 'SH' - Spherical harmonics interpolation [1-5] 0048 % 'NN' - Natural neighbor interpolation with Voronoi weights [6] 0049 % 'Bary' - Barycentric interpolation [7] 0050 % Default: 'SH' 0051 % mc - Define maximum boost in dB if magnitude correction according to [8] should be applied to interpolated HRTFs in a further postprocessing step 0052 % Set to nan if magnitude correction should not be applied 0053 % Set to inf if no limiting should be applied, i.e., unlimited boost 0054 % Default: inf 0055 % headRadius - Head radius in m. Required for time-alignment methods SUpDEq [3][5] and PC [4][5] 0056 % Default: 0.0875 0057 % tikhEps - Define epsilon of Tikhonov regularization if regularization should be applied 0058 % Applying the Tikhonov regularization will always result in a least-square fit 0059 % solution for the SH transform. Only relevant if ipMethod is 'SH'. 0060 % Default: 0 (no Tikhonov regularization) 0061 % limitMC - Boolean to set magnitude correction filters to 0 dB below spatial aliasing frequency fA, 0062 % assuming that (SH) interpolation below fA is physically correct. 0063 % Default: true (correction filter operate only above fA) 0064 % Only applies if 'mc' is not nan 0065 % mcKnee - Knee in dB of the soft-limiting applied to the correction filters. Soft-limiting is applied according to maximum boost 'mc' 0066 % Only applies if 'mc' is not nan or inf 0067 % Default: 0 dB (no knee) 0068 % mcMinPhase - Boolean to design correction filters as minimum-phase filters. If set to false, correction filters are zero-phase filters 0069 % Default: true (minimum-phase filters) as this provides a better performance in terms of ITD errors in combination with classical SUpDEq processing 0070 % When applying other preprocessing methods, using zero-phase filters can be more appropriate. 0071 % limFade - String to define whether to apply a linear fade upwards from aliasing frequency fA to fAt = fA + 1/3 Oct, or downwards, i.e., from fA to fAt = fA - 1/3 Oct 0072 % 'fadeUp' - upwards 0073 % 'fadeDown' - downwards 0074 % Default: 'fadeDown' 0075 % 0076 % Dependencies: SUpDEq toolbox, AKtools, SOFiA toolbox, SFS toolbox, TriangleRayIntersection 0077 % 0078 % References: (Pre-Processing methods) 0079 % Onset-Based Time-Alignment (OBTA) 0080 % [1] M. J. Evans, J. A. S. Angus, and A. I. Tew, "Analyzing head-related transfer function measurements 0081 % using surface spherical harmonics," 0082 % J. Acoust. Soc. Am., vol. 104, no. 4, pp. 2400-2411, 1998. 0083 % [2] F. Brinkmann and S. Weinzierl, "Comparison of head-related transfer functions pre-processing 0084 % techniques for spherical harmonics decomposition," 0085 % in Proceedings of the AES Conference on Audio for Virtual and Augmented Reality, 2018, pp. 1-10. 0086 % 0087 % Spatial Upsampling by Directional Equalization (SUpDEq) 0088 % [3] C. Pörschmann, J. M. Arend, and F. Brinkmann, 0089 % "Directional Equalization of Sparse Head-Related Transfer Function Sets for Spatial Upsampling," 0090 % IEEE/ACM Trans. Audio, Speech, Lang. Process., vol. 27, no. 6, pp. 1060-1071, 2019. 0091 % 0092 % Phase-Correction (PC) 0093 % [4] Z. Ben-Hur, D. Lou Alon, R. Mehra, and B. Rafaely, 0094 % "Efficient Representation and Sparse Sampling of Head-Related Transfer Functions 0095 % Using Phase-Correction Based on Ear Alignment," 0096 % IEEE/ACM Trans. Audio, Speech, Lang. Process., vol. 27, no. 12, pp. 2249-2262, 2019. 0097 % 0098 % Comparison of time-alignment methods & perceptual evaluation 0099 % [5] J. M. Arend, F. Brinkmann, and C. Pörschmann, “Assessing Spherical 0100 % Harmonics Interpolation of Time-Aligned Head-Related Transfer Functions,” 0101 % J. Audio Eng. Soc., vol. 69, no. 1/2, pp. 104–117, 2021. 0102 % 0103 % Spherical harmonics interpolation 0104 % [1-5] 0105 % 0106 % Natural-Neighbor Interpolation 0107 % [6] R. Sibson, “A brief description of natural neighbor interpolation,” 0108 % in Interpreting Multivariate Data, V. Barnett, Chichester, England: John Wiley & Sons, 1981, pp. 21–36. 0109 % Good read about NN interpolation: https://github.com/gwlucastrig/Tinfour/wiki/Introduction-to-Natural-Neighbor-Interpolation 0110 % 0111 % Barycentric Coordinates / Interpolation 0112 % [7] P. Shirley and S. Marschner, Fundamentals of Computer Graphics, 3rd ed. Boca Raton, FL: Taylor & Francis, 2009. 0113 % Good read about Barycentric interpolation: https://www.scratchapixel.com/lessons/3d-basic-rendering/ray-tracing-rendering-a-triangle/barycentric-coordinates 0114 % 0115 % Magnitude-Corrected and Time-Aligned Interpolation (MCA) 0116 % [8] J. M. Arend, C. Pörschmann, S. Weinzierl, F. Brinkmann, 0117 % "Magnitude-Corrected and Time-Aligned Interpolation of Head-Related Transfer Functions," 0118 % (Manuscript submitted for publication). 0119 % 0120 % (C) 2020-2022 by JMA, Johannes M. Arend 0121 % TU Berlin, Audio Communication Group 0122 % TH Köln, Institute of Communications Engineering 0123 0124 function [interpHRTFset, HRTF_L_ip, HRTF_R_ip] = supdeq_interpHRTF(HRTFset, ipSamplingGrid, ppMethod, ipMethod, mc, headRadius, tikhEps, limitMC, mcKnee, mcMinPhase, limFade) 0125 0126 0127 if nargin < 3 || isempty(ppMethod) 0128 ppMethod = 'SUpDEq'; 0129 end 0130 0131 if nargin < 4 || isempty(ipMethod) 0132 ipMethod = 'SH'; 0133 end 0134 0135 if nargin < 5 || isempty(mc) 0136 mc = inf; 0137 end 0138 0139 if nargin < 6 || isempty(headRadius) 0140 headRadius = 0.0875; 0141 end 0142 0143 if nargin < 7 || isempty(tikhEps) 0144 tikhEps = 0; 0145 end 0146 0147 if nargin < 8 || isempty(limitMC) 0148 limitMC = true; 0149 end 0150 0151 if nargin < 9 || isempty(mcKnee) 0152 mcKnee = 0; 0153 end 0154 0155 if nargin < 10 || isempty(mcMinPhase) 0156 mcMinPhase = true; 0157 end 0158 0159 if nargin < 11 || isempty(limFade) 0160 limFade = 'fadeDown'; 0161 end 0162 0163 %% Get some required variables 0164 0165 fs = HRTFset.f(end)*2; 0166 f = HRTFset.f; 0167 NFFT = length(HRTFset.f)*2-2; 0168 samplingGrid = HRTFset.samplingGrid; 0169 HRTF_L = HRTFset.HRTF_L; 0170 HRTF_R = HRTFset.HRTF_R; 0171 c = 343; 0172 0173 %Get HRIRs - Required for threshold detection and mc 0174 %Get mirror spectrum 0175 HRIR_L = AKsingle2bothSidedSpectrum(HRTF_L.'); 0176 HRIR_R = AKsingle2bothSidedSpectrum(HRTF_R.'); 0177 %Transform back to time domain 0178 HRIR_L = real(ifft(HRIR_L)); 0179 HRIR_R = real(ifft(HRIR_R)); 0180 hrirLength = size(HRIR_L,1)/HRTFset.FFToversize; 0181 0182 %% Perform preprocessing 0183 0184 switch ppMethod 0185 case 'SUpDEq' 0186 disp('Preprocessing: SUpDEq') 0187 %Apply SUpDEq - Equalize HRTFs 0188 0189 %Get eqDataset at N = 44 0190 eqDataset = supdeq_getEqDataset(44,2*headRadius,NFFT,fs); 0191 0192 %Equalization 0193 [eqTF_L,eqTF_R] = supdeq_getEqTF(eqDataset,samplingGrid,'DEG',2,'ak',false); 0194 pHRTF_L = HRTF_L./eqTF_L; 0195 pHRTF_R = HRTF_R./eqTF_R; 0196 0197 case 'SUpDEq_Lim' 0198 disp('Preprocessing: SUpDEq_Lim') 0199 %Apply SUpDEq - Equalize HRTFs 0200 0201 %Get eqDataset at N = 44 0202 eqDataset = supdeq_getEqDataset(44,2*headRadius,NFFT,fs); 0203 0204 %Limit eqDataset (Small improvement explained in [5] leading to better 0205 %results below the spatial aliasing frequency fA) 0206 eqDataset = supdeq_limitEqDataset(eqDataset,HRTFset.Nmax,headRadius); 0207 0208 %Equalization 0209 [eqTF_L,eqTF_R] = supdeq_getEqTF(eqDataset,samplingGrid,'DEG',2,'ak',false); 0210 pHRTF_L = HRTF_L./eqTF_L; 0211 pHRTF_R = HRTF_R./eqTF_R; 0212 0213 case 'SUpDEq_AP' 0214 disp('Preprocessing: SUpDEq_AP') 0215 %Apply SUpDEq_AP - Equalize HRTFs with rigid sphere allpass functions 0216 0217 %Get eqDataset at N = 44 0218 eqDataset = supdeq_getEqDataset(44,2*headRadius,NFFT,fs); 0219 0220 %Equalization 0221 [eqTF_L,eqTF_R] = supdeq_getEqTF(eqDataset,samplingGrid,'DEG',2,'ak',true); 0222 pHRTF_L = HRTF_L./eqTF_L; 0223 pHRTF_R = HRTF_R./eqTF_R; 0224 0225 case 'SUpDEq_Lim_AP' 0226 disp('Preprocessing: SUpDEq_Lim_AP') 0227 %Apply SUpDEq_AP - Equalize HRTFs with limited rigid sphere allpass functions 0228 0229 %Get eqDataset at N = 44 0230 eqDataset = supdeq_getEqDataset(44,2*headRadius,NFFT,fs); 0231 0232 %Limit eqDataset (Small improvement explained in [5] leading to better 0233 %results below the spatial aliasing frequency fA) 0234 eqDataset = supdeq_limitEqDataset(eqDataset,HRTFset.Nmax,headRadius); 0235 0236 %Equalization 0237 [eqTF_L,eqTF_R] = supdeq_getEqTF(eqDataset,samplingGrid,'DEG',2,'ak',true); 0238 pHRTF_L = HRTF_L./eqTF_L; 0239 pHRTF_R = HRTF_R./eqTF_R; 0240 0241 case 'PC' 0242 disp('Preprocessing: PC') 0243 %Apply Phase-Correction - Equalize HRTFs with open sphere allpass functions 0244 0245 %Get wavenumber k 0246 k = 2*pi*f/c; k = k.'; 0247 0248 %Transform sampling grid to radiant 0249 sg = samplingGrid * pi / 180; 0250 0251 %Get phase correction term for left/right ear according to Eq. 13 in [4] 0252 cosThetaL = cos(sg(:,2)')*cos(pi/2) + sin(sg(:,2)')*sin(pi/2) .* cos(sg(:,1)'-pi/2); %Left ear with -pi/2 0253 phaseCorL = exp(-1j*headRadius * k .* cosThetaL); 0254 cosThetaR = cos(sg(:,2)')*cos(pi/2) + sin(sg(:,2)')*sin(pi/2) .* cos(sg(:,1)'+pi/2); %Right ear with +pi/2 0255 phaseCorR = exp(-1j*headRadius * k .* cosThetaR); 0256 0257 %Apply to HRTFs 0258 pHRTF_L = HRTF_L .* phaseCorL.'; %Just transpose, no conjugate complex 0259 pHRTF_R = HRTF_R .* phaseCorR.'; 0260 0261 case 'OBTA' 0262 disp('Preprocessing: OBTA') 0263 %Apply Onset-Based Time-Alignment with fractional delay filters 0264 0265 %Estimate TOA on low-passed and 10 times upsampled HRIRs according to [2] 0266 toa_L = AKonsetDetect(HRIR_L, 10, -20, 'rel', [3e3 fs]); 0267 toa_R = AKonsetDetect(HRIR_R, 10, -20, 'rel', [3e3 fs]); 0268 0269 %Time-Align HRIRs (remove TOAs) with fractional delay 0270 pHRIR_L = AKfractionalDelayCyclic(HRIR_L, -toa_L); 0271 pHRIR_R = AKfractionalDelayCyclic(HRIR_R, -toa_R); 0272 0273 %Transform to Fourier domain with same FFToversize as defined in struct 0274 pHRTF_L = AKboth2singleSidedSpectrum (fft(pHRIR_L)).'; 0275 pHRTF_R = AKboth2singleSidedSpectrum (fft(pHRIR_R)).'; 0276 0277 case 'MagPhase' 0278 disp('Preprocessing: MagPhase') 0279 0280 pHRTF_L = abs(HRTF_L); 0281 pHRTF_R = abs(HRTF_R); 0282 pHRTF_L_ph = unwrap(angle(HRTF_L).').'; 0283 pHRTF_R_ph = unwrap(angle(HRTF_R).').'; 0284 0285 case 'None' 0286 disp('Preprocessing: None') 0287 0288 %Simply copy HRTF to pHRTF_L - Easier for further processing 0289 pHRTF_L = HRTF_L; 0290 pHRTF_R = HRTF_R; 0291 end 0292 0293 %Save intermediate result (pre-processed HRTFs) in output struct 0294 interpHRTFset.p.pHRTF_L = pHRTF_L; 0295 interpHRTFset.p.pHRTF_R = pHRTF_R; 0296 if strcmp(ppMethod,'OBTA') 0297 interpHRTFset.p.toa_L = toa_L; 0298 interpHRTFset.p.toa_R = toa_R; 0299 end 0300 if strcmp(ppMethod,'MagPhase') 0301 interpHRTFset.p.pHRTF_L_ph = pHRTF_L_ph; 0302 interpHRTFset.p.pHRTF_R_ph = pHRTF_R_ph; 0303 end 0304 0305 %% Perform interpolation 0306 0307 switch ipMethod 0308 case 'SH' 0309 disp('HRTF Interpolation: SH') 0310 0311 %Transform preprocessed HRTFs to SH domain 0312 %Tikhonov regularization can be applied if weights in sampling grid are erased and tikhEps ~= 0 0313 pHRTF_L_nm = AKsht(pHRTF_L.',false,samplingGrid,HRTFset.Nmax,'complex',fs,false,'complex',tikhEps); 0314 pHRTF_R_nm = AKsht(pHRTF_R.',false,samplingGrid,HRTFset.Nmax,'complex',fs,false,'complex',tikhEps); 0315 0316 %Apply SH interpolation to ipSamplingGrid 0317 pHRTF_L_ip = AKisht(pHRTF_L_nm,false,ipSamplingGrid(:,1:2),'complex').'; 0318 pHRTF_R_ip = AKisht(pHRTF_R_nm,false,ipSamplingGrid(:,1:2),'complex').'; 0319 0320 %Also interpolate TOAs if ppMethod = 'OBTA' 0321 if strcmp(ppMethod,'OBTA') 0322 0323 %Transform TOAs to SH domain 0324 %Tikhonov regularization can be applied if weights in sampling grid are erased and tikhEps ~= 0 0325 toa_L_nm = AKsht(toa_L,false,samplingGrid,HRTFset.Nmax,'complex',fs,false,'real',tikhEps); 0326 toa_R_nm = AKsht(toa_R,false,samplingGrid,HRTFset.Nmax,'complex',fs,false,'real',tikhEps); 0327 0328 %Apply SH interpolation to ipSamplingGrid 0329 toa_L_ip = AKisht(toa_L_nm, false, ipSamplingGrid(:,1:2), 'complex', false, false, 'real'); 0330 toa_R_ip = AKisht(toa_R_nm, false, ipSamplingGrid(:,1:2), 'complex', false, false, 'real'); 0331 0332 end 0333 0334 %Also interpolate Phase if ppMethod = 'MagPhase' 0335 if strcmp(ppMethod,'MagPhase') 0336 0337 %Transform Phase to SH domain 0338 %Tikhonov regularization can be applied if weights in sampling grid are erased and tikhEps ~= 0 0339 pHRTF_L_ph_nm = AKsht(pHRTF_L_ph.',false,samplingGrid,HRTFset.Nmax,'complex',fs,false,'complex',tikhEps); 0340 pHRTF_R_ph_nm = AKsht(pHRTF_R_ph.',false,samplingGrid,HRTFset.Nmax,'complex',fs,false,'complex',tikhEps); 0341 0342 %Apply SH interpolation to ipSamplingGrid 0343 pHRTF_L_ph_ip = AKisht(pHRTF_L_ph_nm,false,ipSamplingGrid(:,1:2),'complex').'; 0344 pHRTF_R_ph_ip = AKisht(pHRTF_R_ph_nm,false,ipSamplingGrid(:,1:2),'complex').'; 0345 0346 end 0347 0348 case 'NN' 0349 disp('HRTF Interpolation: NN') 0350 0351 %Transform grids to cartesian coordinates 0352 [samplingGrid_cart(:,1), samplingGrid_cart(:,2), samplingGrid_cart(:,3)] = sph2cart(samplingGrid(:,1)/360*2*pi, pi/2-samplingGrid(:,2)/360*2*pi,1); 0353 [ipSamplingGrid_cart(:,1), ipSamplingGrid_cart(:,2), ipSamplingGrid_cart(:,3)] = sph2cart(ipSamplingGrid(:,1)/360*2*pi, pi/2-ipSamplingGrid(:,2)/360*2*pi,1); 0354 0355 %NN interpolate preprocessed HRTFs to ipSamplingGrid 0356 pHRTF_L_ip = zeros(length(ipSamplingGrid),length(f)); 0357 pHRTF_R_ip = zeros(length(ipSamplingGrid),length(f)); 0358 0359 for nPoint = 1:size(ipSamplingGrid,1) 0360 0361 [idx_voronoi,w_voronoi] = findvoronoi(samplingGrid_cart,ipSamplingGrid_cart(nPoint,:)); 0362 0363 for n = 1:length(idx_voronoi) 0364 pHRTF_L_ip(nPoint,:) = pHRTF_L_ip(nPoint,:) + w_voronoi(n)*pHRTF_L(idx_voronoi(n),:); 0365 pHRTF_R_ip(nPoint,:) = pHRTF_R_ip(nPoint,:) + w_voronoi(n)*pHRTF_R(idx_voronoi(n),:); 0366 end 0367 end 0368 0369 %Also interpolate TOAs if ppMethod = 'OBTA' 0370 if strcmp(ppMethod,'OBTA') 0371 0372 toa_L_ip = zeros(1,length(ipSamplingGrid)); 0373 toa_R_ip = zeros(1,length(ipSamplingGrid)); 0374 0375 for nPoint = 1:size(ipSamplingGrid,1) 0376 0377 [idx_voronoi,w_voronoi] = findvoronoi(samplingGrid_cart,ipSamplingGrid_cart(nPoint,:)); 0378 0379 for n = 1:length(idx_voronoi) 0380 toa_L_ip(1,nPoint) = toa_L_ip(1,nPoint) + w_voronoi(n)*toa_L(1,idx_voronoi(n)); 0381 toa_R_ip(1,nPoint) = toa_R_ip(1,nPoint) + w_voronoi(n)*toa_R(1,idx_voronoi(n)); 0382 end 0383 end 0384 end 0385 0386 %Also interpolate Phase if ppMethod = 'MagPhase' 0387 if strcmp(ppMethod,'MagPhase') 0388 0389 pHRTF_L_ph_ip = zeros(length(ipSamplingGrid),length(f)); 0390 pHRTF_R_ph_ip = zeros(length(ipSamplingGrid),length(f)); 0391 0392 for nPoint = 1:size(ipSamplingGrid,1) 0393 0394 [idx_voronoi,w_voronoi] = findvoronoi(samplingGrid_cart,ipSamplingGrid_cart(nPoint,:)); 0395 0396 for n = 1:length(idx_voronoi) 0397 pHRTF_L_ph_ip(nPoint,:) = pHRTF_L_ph_ip(nPoint,:) + w_voronoi(n)*pHRTF_L_ph(idx_voronoi(n),:); 0398 pHRTF_R_ph_ip(nPoint,:) = pHRTF_R_ph_ip(nPoint,:) + w_voronoi(n)*pHRTF_R_ph(idx_voronoi(n),:); 0399 end 0400 end 0401 end 0402 0403 0404 case 'Bary' 0405 disp('HRTF Interpolation: Bary') 0406 0407 %Transform grids to cartesian coordinates 0408 [samplingGrid_cart(:,1), samplingGrid_cart(:,2), samplingGrid_cart(:,3)] = sph2cart(samplingGrid(:,1)/360*2*pi, pi/2-samplingGrid(:,2)/360*2*pi,1); 0409 [ipSamplingGrid_cart(:,1), ipSamplingGrid_cart(:,2), ipSamplingGrid_cart(:,3)] = sph2cart(ipSamplingGrid(:,1)/360*2*pi, pi/2-ipSamplingGrid(:,2)/360*2*pi,1); 0410 0411 %Get simplified convex hull (could also be done after delaunay-triangulation) 0412 convHullIdx = convhull(samplingGrid_cart(:,1), samplingGrid_cart(:,2), samplingGrid_cart(:,3), 'simplify', true); 0413 0414 %Get HRTFs / directions assigned to convex hull triangles 0415 convHullHRTF_A = samplingGrid_cart(convHullIdx(:,1),:); 0416 convHullHRTF_B = samplingGrid_cart(convHullIdx(:,2),:); 0417 convHullHRTF_C = samplingGrid_cart(convHullIdx(:,3),:); 0418 0419 %Interpolate with barycentric weights 0420 %Right triangle of convex hull picked by ray-triangle intersection test (intersectLinePolygon3d from geom3D toolbox should work too) 0421 %Function returns barycentric weights / coordinates u and v of the intersection point -> w = 1-u-v 0422 %P = w*A + u*B + v*C 0423 orig = [0 0 0]; 0424 pHRTF_L_ip = zeros(length(ipSamplingGrid),length(f)); 0425 pHRTF_R_ip = zeros(length(ipSamplingGrid),length(f)); 0426 0427 for nPoint = 1:size(ipSamplingGrid,1) 0428 0429 [intersect, ~, u, v, ~] = TriangleRayIntersection(orig, ipSamplingGrid_cart(nPoint,:), convHullHRTF_A, convHullHRTF_B, convHullHRTF_C, 'border', 'inclusive'); 0430 intersectIdx = find(intersect, 1, 'first'); %Evaluate first intersection point 0431 u = u(intersectIdx); v = v(intersectIdx); w = 1-u-v; 0432 bw = [w u v]; %Barycentric weights 0433 idx_bw = convHullIdx(intersectIdx,:); %Indizes of HRTFs of convex hull / triangles 0434 0435 for n = 1:3 %Always 3 because of triangulization 0436 pHRTF_L_ip(nPoint,:) = pHRTF_L_ip(nPoint,:) + bw(n)*pHRTF_L(idx_bw(n),:); 0437 pHRTF_R_ip(nPoint,:) = pHRTF_R_ip(nPoint,:) + bw(n)*pHRTF_R(idx_bw(n),:); 0438 end 0439 end 0440 0441 %Also interpolate TOAs if ppMethod = 'OBTA' 0442 if strcmp(ppMethod,'OBTA') 0443 0444 toa_L_ip = zeros(1,length(ipSamplingGrid)); 0445 toa_R_ip = zeros(1,length(ipSamplingGrid)); 0446 0447 for nPoint = 1:size(ipSamplingGrid,1) 0448 0449 [intersect, ~, u, v, ~] = TriangleRayIntersection(orig, ipSamplingGrid_cart(nPoint,:), convHullHRTF_A, convHullHRTF_B, convHullHRTF_C, 'border', 'inclusive'); 0450 intersectIdx = find(intersect, 1, 'first'); %Evaluate first intersection point 0451 u = u(intersectIdx); v = v(intersectIdx); w = 1-u-v; 0452 bw = [w u v]; %Barycentric weights 0453 idx_bw = convHullIdx(intersectIdx,:); %Indizes of HRTFs of convex hull / triangles 0454 0455 for n = 1:3 %Always 3 because of triangulization 0456 toa_L_ip(1,nPoint) = toa_L_ip(1,nPoint) + bw(n)*toa_L(1,idx_bw(n)); 0457 toa_R_ip(1,nPoint) = toa_R_ip(1,nPoint) + bw(n)*toa_R(1,idx_bw(n)); 0458 end 0459 end 0460 end 0461 0462 %Also interpolate Phase if ppMethod = 'MagPhase' 0463 if strcmp(ppMethod,'MagPhase') 0464 0465 pHRTF_L_ph_ip = zeros(length(ipSamplingGrid),length(f)); 0466 pHRTF_R_ph_ip = zeros(length(ipSamplingGrid),length(f)); 0467 0468 for nPoint = 1:size(ipSamplingGrid,1) 0469 0470 [intersect, ~, u, v, ~] = TriangleRayIntersection(orig, ipSamplingGrid_cart(nPoint,:), convHullHRTF_A, convHullHRTF_B, convHullHRTF_C, 'border', 'inclusive'); 0471 intersectIdx = find(intersect, 1, 'first'); %Evaluate first intersection point 0472 u = u(intersectIdx); v = v(intersectIdx); w = 1-u-v; 0473 bw = [w u v]; %Barycentric weights 0474 idx_bw = convHullIdx(intersectIdx,:); %Indizes of HRTFs of convex hull / triangles 0475 0476 for n = 1:3 %Always 3 because of triangulization 0477 pHRTF_L_ph_ip(nPoint,:) = pHRTF_L_ph_ip(nPoint,:) + bw(n)*pHRTF_L_ph(idx_bw(n),:); 0478 pHRTF_R_ph_ip(nPoint,:) = pHRTF_R_ph_ip(nPoint,:) + bw(n)*pHRTF_R_ph(idx_bw(n),:); 0479 end 0480 end 0481 end 0482 0483 end 0484 0485 %Save intermediate result (interpolated pre-processed HRTFs) in output struct 0486 interpHRTFset.p.pHRTF_L_ip = pHRTF_L_ip; 0487 interpHRTFset.p.pHRTF_R_ip = pHRTF_R_ip; 0488 if strcmp(ppMethod,'OBTA') 0489 interpHRTFset.p.toa_L_ip = toa_L_ip; 0490 interpHRTFset.p.toa_R_ip = toa_R_ip; 0491 end 0492 if strcmp(ppMethod,'MagPhase') 0493 interpHRTFset.p.pHRTF_L_ph_ip = pHRTF_L_ph_ip; 0494 interpHRTFset.p.pHRTF_R_ph_ip = pHRTF_R_ph_ip; 0495 end 0496 0497 %% Perform postprocessing 0498 0499 switch ppMethod 0500 case 'SUpDEq' 0501 disp('Postprocessing: SUpDEq') 0502 %Apply SUpDEq - De-Equalize interpolated HRTFs 0503 0504 %Get De-Equalization functions 0505 [deqTF_L,deqTF_R] = supdeq_getEqTF(eqDataset,ipSamplingGrid,'DEG',2,'ak',false); 0506 0507 %De-Equalize 0508 HRTF_L_ip = pHRTF_L_ip.*deqTF_L; 0509 HRTF_R_ip = pHRTF_R_ip.*deqTF_R; 0510 0511 case 'SUpDEq_Lim' 0512 disp('Postprocessing: SUpDEq_Lim') 0513 %Apply SUpDEq - De-Equalize interpolated HRTFs 0514 0515 %Get De-Equalization functions of limited eqDataset 0516 [deqTF_L,deqTF_R] = supdeq_getEqTF(eqDataset,ipSamplingGrid,'DEG',2,'ak',false); 0517 0518 %De-Equalize 0519 HRTF_L_ip = pHRTF_L_ip.*deqTF_L; 0520 HRTF_R_ip = pHRTF_R_ip.*deqTF_R; 0521 0522 case 'SUpDEq_AP' 0523 disp('Postprocessing: SUpDEq_AP') 0524 %Apply SUpDEq_AP - De-Equalize interpolated HRTFs with rigid sphere allpass functions 0525 0526 %Get De-Equalization functions 0527 [deqTF_L,deqTF_R] = supdeq_getEqTF(eqDataset,ipSamplingGrid,'DEG',2,'ak',true); 0528 0529 %De-Equalize 0530 HRTF_L_ip = pHRTF_L_ip.*deqTF_L; 0531 HRTF_R_ip = pHRTF_R_ip.*deqTF_R; 0532 0533 case 'SUpDEq_Lim_AP' 0534 disp('Postprocessing: SUpDEq_Lim_AP') 0535 %Apply SUpDEq_Lim_AP - De-Equalize interpolated HRTFs with limited rigid sphere allpass functions 0536 0537 %Get De-Equalization functions of limited eqDataset 0538 [deqTF_L,deqTF_R] = supdeq_getEqTF(eqDataset,ipSamplingGrid,'DEG',2,'ak',true); 0539 0540 %De-Equalize 0541 HRTF_L_ip = pHRTF_L_ip.*deqTF_L; 0542 HRTF_R_ip = pHRTF_R_ip.*deqTF_R; 0543 0544 case 'PC' 0545 disp('Postprocessing: PC') 0546 %Apply Inverse Phase-Correction 0547 0548 %Transform interpolation sampling grid to radiant 0549 sg = ipSamplingGrid * pi / 180; 0550 0551 %Get inverse phase correction term for left/right ear according to Eq. 13 in [4] 0552 cosThetaL = cos(sg(:,2)')*cos(pi/2) + sin(sg(:,2)')*sin(pi/2) .* cos(sg(:,1)'-pi/2); 0553 %No minus before 1j to get inverse phase term. Could also be achieved by division instead of multiplication 0554 phaseCorL = exp(1j*headRadius * k .* cosThetaL); 0555 cosThetaR = cos(sg(:,2)')*cos(pi/2) + sin(sg(:,2)')*sin(pi/2) .* cos(sg(:,1)'+pi/2); 0556 %No minus before 1j to get inverse phase term. Could also be achieved by division instead of multiplication 0557 phaseCorR = exp(1j*headRadius * k .* cosThetaR); 0558 0559 %Apply to interpolated HRTFs 0560 HRTF_L_ip = pHRTF_L_ip .* phaseCorL.'; %Just transpose, no conjugate complex 0561 HRTF_R_ip = pHRTF_R_ip .* phaseCorR.'; 0562 0563 case 'OBTA' 0564 disp('Postprocessing: OBTA') 0565 %Re-Insert interpolated TOAs with fractional delays 0566 0567 %Get interpolated HRIRs 0568 pHRIR_L_ip = AKsingle2bothSidedSpectrum(pHRTF_L_ip.'); 0569 pHRIR_R_ip = AKsingle2bothSidedSpectrum(pHRTF_R_ip.'); 0570 pHRIR_L_ip = real(ifft(pHRIR_L_ip)); 0571 pHRIR_R_ip = real(ifft(pHRIR_R_ip)); 0572 0573 %Re-Insert TOAs with fractional delays 0574 HRIR_L_ip = AKfractionalDelayCyclic(pHRIR_L_ip, toa_L_ip); 0575 HRIR_R_ip = AKfractionalDelayCyclic(pHRIR_R_ip, toa_R_ip); 0576 0577 %Transform back to frequency domain 0578 HRTF_L_ip = AKboth2singleSidedSpectrum(fft(HRIR_L_ip)).'; 0579 HRTF_R_ip = AKboth2singleSidedSpectrum(fft(HRIR_R_ip)).'; 0580 0581 case 'MagPhase' 0582 disp('Postprocessing: MagPhase') 0583 0584 %Combine magnitude and phase of interpolated HRTFs again 0585 HRTF_L_ip = pHRTF_L_ip.*exp(1i*pHRTF_L_ph_ip); 0586 HRTF_R_ip = pHRTF_R_ip.*exp(1i*pHRTF_R_ph_ip); 0587 0588 case 'None' 0589 disp('Postprocessing: None') 0590 0591 %Simply copy pHRTF_L_ip to HRTF_L_ip - Easier for further processing 0592 HRTF_L_ip = pHRTF_L_ip; 0593 HRTF_R_ip = pHRTF_R_ip; 0594 end 0595 0596 %% Correct magnitude if mc ~= nan 0597 0598 if ~isnan(mc) 0599 0600 disp('MC postprocessing'); 0601 0602 %Save intermediate result (interpolated HRTFs without mc) in output struct 0603 interpHRTFset.p.HRTF_L_ip_noMC = HRTF_L_ip; 0604 interpHRTFset.p.HRTF_R_ip_noMC = HRTF_R_ip; 0605 0606 %Get input HRTFs in auditory bands (ERB filter) 0607 fLowErb = 50; 0608 [cl,ferb] = AKerbErrorPersistent(HRIR_L(1:hrirLength,:),AKdirac(hrirLength),[fLowErb fs/2],fs); 0609 [cr] = AKerbErrorPersistent(HRIR_R(1:hrirLength,:),AKdirac(hrirLength),[fLowErb fs/2],fs); 0610 0611 %Interpolate ERB filters to ipSamplingGrid with respective ipMethod 0612 switch ipMethod 0613 case 'SH' 0614 0615 %Transform ERB filters to SH domain 0616 %Tikhonov regularization can be applied if weights in sampling grid are erased and tikhEps ~= 0 0617 cl_nm = AKsht(cl, false, samplingGrid, HRTFset.Nmax, 'complex', fs, false, 'real',tikhEps); 0618 cr_nm = AKsht(cr, false, samplingGrid, HRTFset.Nmax, 'complex', fs, false, 'real',tikhEps); 0619 0620 %Apply SH interpolation to ipSamplingGrid 0621 cl_ip = AKisht(cl_nm, false, ipSamplingGrid(:,1:2), 'complex', false, false, 'real'); 0622 cr_ip = AKisht(cr_nm, false, ipSamplingGrid(:,1:2), 'complex', false, false, 'real'); 0623 0624 case 'NN' 0625 0626 %Transform grids to cartesian coordinates 0627 [samplingGrid_cart(:,1), samplingGrid_cart(:,2), samplingGrid_cart(:,3)] = sph2cart(samplingGrid(:,1)/360*2*pi, pi/2-samplingGrid(:,2)/360*2*pi,1); 0628 [ipSamplingGrid_cart(:,1), ipSamplingGrid_cart(:,2), ipSamplingGrid_cart(:,3)] = sph2cart(ipSamplingGrid(:,1)/360*2*pi, pi/2-ipSamplingGrid(:,2)/360*2*pi,1); 0629 0630 %NN interpolate ERB filters to ipSamplingGrid 0631 cl_ip = zeros(length(ferb),length(ipSamplingGrid)); 0632 cr_ip = zeros(length(ferb),length(ipSamplingGrid)); 0633 0634 for nPoint = 1:size(ipSamplingGrid,1) 0635 0636 [idx_voronoi,w_voronoi] = findvoronoi(samplingGrid_cart,ipSamplingGrid_cart(nPoint,:)); 0637 0638 for n = 1:length(idx_voronoi) 0639 cl_ip(:,nPoint) = cl_ip(:,nPoint) + w_voronoi(n)*cl(:,idx_voronoi(n)); 0640 cr_ip(:,nPoint) = cr_ip(:,nPoint) + w_voronoi(n)*cr(:,idx_voronoi(n)); 0641 end 0642 end 0643 0644 case 'Bary' 0645 0646 %Transform grids to cartesian coordinates 0647 [samplingGrid_cart(:,1), samplingGrid_cart(:,2), samplingGrid_cart(:,3)] = sph2cart(samplingGrid(:,1)/360*2*pi, pi/2-samplingGrid(:,2)/360*2*pi,1); 0648 [ipSamplingGrid_cart(:,1), ipSamplingGrid_cart(:,2), ipSamplingGrid_cart(:,3)] = sph2cart(ipSamplingGrid(:,1)/360*2*pi, pi/2-ipSamplingGrid(:,2)/360*2*pi,1); 0649 0650 %Get simplified convex hull (could also be done after delaunay-triangulation) 0651 convHullIdx = convhull(samplingGrid_cart(:,1), samplingGrid_cart(:,2), samplingGrid_cart(:,3), 'simplify', true); 0652 0653 %Get HRTFs / directions assigned to convex hull triangles 0654 convHullHRTF_A = samplingGrid_cart(convHullIdx(:,1),:); 0655 convHullHRTF_B = samplingGrid_cart(convHullIdx(:,2),:); 0656 convHullHRTF_C = samplingGrid_cart(convHullIdx(:,3),:); 0657 0658 %Interpolate with barycentric weights 0659 %Right triangle of convex hull picked by ray-triangle intersection test (intersectLinePolygon3d from geom3D toolbox should work too) 0660 %Function returns barycentric weights / coordinates u and v of the intersection point -> w = 1-u-v 0661 %P = w*A + u*B + v*C 0662 orig = [0 0 0]; 0663 cl_ip = zeros(length(ferb),length(ipSamplingGrid)); 0664 cr_ip = zeros(length(ferb),length(ipSamplingGrid)); 0665 0666 for nPoint = 1:size(ipSamplingGrid,1) 0667 0668 [intersect, ~, u, v, ~] = TriangleRayIntersection(orig, ipSamplingGrid_cart(nPoint,:), convHullHRTF_A, convHullHRTF_B, convHullHRTF_C, 'border', 'inclusive'); 0669 intersectIdx = find(intersect, 1, 'first'); %Evaluate first intersection point 0670 u = u(intersectIdx); v = v(intersectIdx); w = 1-u-v; 0671 bw = [w u v]; %Barycentric weights 0672 idx_bw = convHullIdx(intersectIdx,:); %Indizes of HRTFs of convex hull / triangles 0673 0674 for n = 1:3 %Always 3 because of triangulization 0675 cl_ip(:,nPoint) = cl_ip(:,nPoint) + bw(n)*cl(:,idx_bw(n)); 0676 cr_ip(:,nPoint) = cr_ip(:,nPoint) + bw(n)*cr(:,idx_bw(n)); 0677 end 0678 end 0679 end 0680 0681 %Save interpolated ERBs in 3D-array for further processing 0682 c_ip(:,:,1) = cl_ip; 0683 c_ip(:,:,2) = cr_ip; 0684 0685 %Get interpolated HRTFs in auditory bands (ERB filter) 0686 HRIR_L_ip = AKsingle2bothSidedSpectrum(HRTF_L_ip.'); 0687 HRIR_R_ip = AKsingle2bothSidedSpectrum(HRTF_R_ip.'); 0688 HRIR_ip(:,:,1) = real(ifft(HRIR_L_ip)); 0689 HRIR_ip(:,:,2) = real(ifft(HRIR_R_ip)); 0690 HRIR_ip = HRIR_ip(1:hrirLength,:,:); 0691 %Save in 3D-array for further processing 0692 c_hrir_ip = AKerbErrorPersistent(HRIR_ip,AKdirac(size(HRIR_ip,1)),[fLowErb fs/2],fs); 0693 0694 %Get correction filters for all HRTFs 0695 corrFilt = c_ip-c_hrir_ip; 0696 %Spline interpolate correction filters to 0-fs/2 0697 corrFilt_l = AKinterpolateSpectrum( squeeze(corrFilt(:,:,1)), ferb, NFFT, {'nearest' 'spline' 'nearest'}, fs); 0698 corrFilt_r = AKinterpolateSpectrum( squeeze(corrFilt(:,:,2)), ferb, NFFT, {'nearest' 'spline' 'nearest'}, fs); 0699 0700 %Limit to f > fA (set to 0 below fA) if desired 0701 if limitMC 0702 0703 %Get spatial alasing frequency 0704 fA = HRTFset.Nmax*c/2/pi/headRadius; 0705 0706 disp(['MC postprocessing limited to f > fA, with fA = ',num2str(round(fA,2)),'Hz']); 0707 0708 if strcmp(limFade,'fadeDown') 0709 0710 disp('MC fade downwards'); 0711 0712 %Get third-octave below fA 0713 fAt = fA/2^(1/3); 0714 0715 %Set to 0 (logarithmic) below fAt and apply linear fade from fAt to fA 0716 [~,fA_bin] = min(abs(HRTFset.f-fA)); 0717 [~,fAt_bin] = min(abs(HRTFset.f-fAt)); 0718 ramp = linspace(0,1,abs(fA_bin-fAt_bin)+1); 0719 corrFilt_l(1:fAt_bin-1,:) = 0; 0720 corrFilt_l(fAt_bin:fA_bin,:) = corrFilt_l(fAt_bin:fA_bin,:).*ramp.'; 0721 corrFilt_r(1:fAt_bin-1,:) = 0; 0722 corrFilt_r(fAt_bin:fA_bin,:) = corrFilt_r(fAt_bin:fA_bin,:).*ramp.'; 0723 end 0724 0725 if strcmp(limFade,'fadeUp') 0726 0727 disp('MC fade upwards'); 0728 0729 %Get third-octave above fA 0730 fAt = fA*2^(1/3); 0731 0732 %Set to 0 (logarithmic) below fA and apply linear fade from fA to fAt 0733 [~,fA_bin] = min(abs(HRTFset.f-fA)); 0734 [~,fAt_bin] = min(abs(HRTFset.f-fAt)); 0735 ramp = linspace(0,1,abs(fAt_bin-fA_bin)+1); 0736 corrFilt_l(1:fA_bin-1,:) = 0; 0737 corrFilt_l(fA_bin:fAt_bin,:) = corrFilt_l(fA_bin:fAt_bin,:).*ramp.'; 0738 corrFilt_r(1:fA_bin-1,:) = 0; 0739 corrFilt_r(fA_bin:fAt_bin,:) = corrFilt_r(fA_bin:fAt_bin,:).*ramp.'; 0740 end 0741 0742 end 0743 0744 %Transform to linear values 0745 corrFilt_l = 10.^(corrFilt_l/20); 0746 corrFilt_r = 10.^(corrFilt_r/20); 0747 0748 if ~isinf(mc) 0749 %Apply soft-limiting to correction filters according to mc if not inf 0750 disp(['MC maximum boost: ',num2str(mc),'dB']) 0751 disp(['MC soft-limiting knee: ',num2str(mcKnee),'dB']) 0752 0753 corrFilt_l_lim = AKsoftLimit(corrFilt_l, mc, mcKnee,[0 fs/2], fs, true); 0754 corrFilt_r_lim = AKsoftLimit(corrFilt_r, mc, mcKnee,[0 fs/2], fs, true); 0755 0756 else %No soft-limiting / Unlimited gain 0757 0758 disp('MC maximum boost: inf') 0759 0760 corrFilt_l_lim = corrFilt_l; 0761 corrFilt_r_lim = corrFilt_r; 0762 0763 end 0764 0765 %Design minimum phase filters and use for correction instead of zero 0766 %phase filters (if mcMinPhase = false) 0767 if mcMinPhase 0768 0769 disp('MC phase: minimum'); 0770 0771 corrFilt_l_lim = AKsingle2bothSidedSpectrum(corrFilt_l_lim); 0772 corrFilt_r_lim = AKsingle2bothSidedSpectrum(corrFilt_r_lim); 0773 0774 corrFilt_l_lim = real(ifft(corrFilt_l_lim)); 0775 corrFilt_r_lim = real(ifft(corrFilt_r_lim)); 0776 0777 corrFilt_l_lim = AKphaseManipulation(corrFilt_l_lim,fs,'min',1,0); 0778 corrFilt_r_lim = AKphaseManipulation(corrFilt_r_lim,fs,'min',1,0); 0779 0780 %Go back to frequency domain 0781 corrFilt_l_lim = AKboth2singleSidedSpectrum(fft(corrFilt_l_lim)); 0782 corrFilt_r_lim = AKboth2singleSidedSpectrum(fft(corrFilt_r_lim)); 0783 0784 else 0785 0786 disp('MC phase: zero'); 0787 0788 end 0789 0790 %Write in 3D array 0791 corrFilt_lim(:,:,1) = corrFilt_l_lim; 0792 corrFilt_lim(:,:,2) = corrFilt_r_lim; 0793 0794 %Save intermediate results (correction filter) in output struct 0795 interpHRTFset.p.corrFilt_lim = corrFilt_lim; 0796 0797 %Apply magnitude correction filters to HRTFs 0798 HRTF_L_ip = HRTF_L_ip.*corrFilt_lim(:,:,1).'; 0799 HRTF_R_ip = HRTF_R_ip.*corrFilt_lim(:,:,2).'; 0800 0801 end 0802 0803 %% Write result in struct 0804 0805 %Get final HRIRs 0806 HRIR_L_ip = AKsingle2bothSidedSpectrum(HRTF_L_ip.'); 0807 HRIR_R_ip = AKsingle2bothSidedSpectrum(HRTF_R_ip.'); 0808 HRIR_L_ip = real(ifft(HRIR_L_ip)); 0809 HRIR_R_ip = real(ifft(HRIR_R_ip)); 0810 %Cut zeros depending on FFToversize 0811 HRIR_L_ip = HRIR_L_ip(1:hrirLength,:); 0812 HRIR_R_ip = HRIR_R_ip(1:hrirLength,:); 0813 0814 interpHRTFset.HRTF_L = HRTF_L_ip; 0815 interpHRTFset.HRTF_R = HRTF_R_ip; 0816 interpHRTFset.HRIR_L = HRIR_L_ip; 0817 interpHRTFset.HRIR_R = HRIR_R_ip; 0818 interpHRTFset.f = HRTFset.f; 0819 interpHRTFset.fs = fs; 0820 interpHRTFset.FFToversize = HRTFset.FFToversize; 0821 interpHRTFset.samplingGrid = ipSamplingGrid; 0822 interpHRTFset.headRadius = headRadius; 0823 interpHRTFset.ppMethod = ppMethod; 0824 interpHRTFset.ipMethod = ipMethod; 0825 interpHRTFset.headRadius = headRadius; 0826 if tikhEps ~= 0 0827 interpHRTFset.tikhEps = tikhEps; 0828 end 0829 if ~isnan(mc) 0830 interpHRTFset.mc = mc; 0831 interpHRTFset.mcKnee = mcKnee; 0832 interpHRTFset.mcMinPhase = mcMinPhase; 0833 if limitMC 0834 interpHRTFset.limitMC = 1; 0835 interpHRTFset.limitMC_fA = fA; 0836 interpHRTFset.limitMC_fAt = fAt; 0837 interpHRTFset.limitMC_fade = limFade; 0838 end 0839 end 0840 disp('Done...'); 0841 0842 end 0843