% SUpDEq - Spatial Upsampling by Directional Equalization
function HRIRs_sfd_ext = supdeq_rangeExt(HRIRs_sfd, r1, regMethod, ampLimit, compTimeShift, compEnergy, rmin, r0)
This function performs range extrapolation of HRTFs in SH-domain. This
way, near-field HRTFs can be obtained from far-field HRTFs. In brief,
radial propagation functions are applied to the SH-coefficients to
calculate the range-extrapolated SH-coefficients for the new distance r1.?
As the radial functions (in case of an exterior problem the hankel
quotient) increase dramatically towards lower kr, a regularization method
like for example soft-limiting or the Tikhonov regularization needs to be
applied. See references for more information about range extrapolation
and regularization approaches.
Output:
HRIRs_sfd_ext - Struct with Spherical Harmonic coefficients
(SH-Coefficients) for the left (Hl_nm) and right (Hr_nm)
channel/ear of the range-extrapolated HRTFs, absolute
frequency scale f, transform order N, FFToversize,
(new range-extrapolated) source distance r1,
non-regularized radial functions (hankel quotient),
regularized (and applied) hankel quotient, and info
about applied regularization and range extrapolation
Input:
HRIRs_sfd - Struct with Spherical Harmonic coefficients
(SH-Coefficients) for the left (Hl_nm) and right (Hr_nm)
channel/ear, absolute frequency scale f,
transform order N, FFToversize, and source distance r0
r1 - New range-extrapolated source distance r1 in m
regMethod - Regularization method (string parameter) for the
radial function (the hankel quotient). The following
methods are available:
'none' - no regularization
'disc' - simply discard contributions where the
radial filter magnitude exceeds ampLimit
'disc_kr' - discard contributions at higher orders
based on kr_min. Results are pretty
similar to normal 'disc'
'hardLimit'- set a hard limit for all radial filter
components at ampLimit instead of
dropping them
'softLimit'- soft-limit the components
'Tikh' - apply Tikhonov regularization
See Bernschütz2011 and Rettberg2014 for an overview
of most regularization methods. Implementation of most
of the regularization methods according to Rettberg2014
Regularization 'disc_kr' implemented according to
Duraiswami2004 and Pollow2012
Default: 'softLimit'
ampLimit - Amplitude limit of radial functions in dB. Not needed
for all regularization methods, but since the default
regularization is 'softLimiting', the amplitude limit
can be freely chosen here.
Default: 20 dB
compTimeShift - true/false - compensate time shift induced by
range-extrapolation with an inverse phase term
Default: true
compEnergy - true/false - compensate the in-/decrease in energy
averaged over all orders and modes
(single value compensation)
Default: true
rmin - Smallest radius in m, meaning the radius of the
smallest sphere enclosing the listener?s head
with all its significant scattering sources. This
defines the truncation order N = [krmin] when
regMethod 'disc_kr' according to Duraiswami2004 and
Pollow2012 is used. Otherwise, rmin has no influence.
Default: 0.1
r0 - Actual source distance in m input HRIRs_sfd.
The HRIRs_sfd structs have a field with the actual
source distance and the function will use the value
in this field as r0. However, if this field is
missing, r0 can be given as an argument.
Default: [] (use value in HRIRs_sfd.sourceDistance)
Dependencies: AKtools
References:
R. Duraiswami, D. N. Zotkin, and N. A. Gumerov, ?Interpolation and range
extrapolation of HRTFs,? in Proceedings of the IEEE International
Conference on Acoustics, Speech, and Signal Processing, 2004, pp. IV45-IV48.
M. Pollow, K.-V. Nguyen, O. Warusfel, T. Carpentier, M. Müller-Trapet,
M. Vorländer, and M. Noisternig, ?Calculation of Head-Related Transfer
Functions for Arbitrary Field Points Using Spherical Harmonics
Decomposition,? Acta Acust. united Ac., vol. 98, no. 1, pp. 72?82, 2012.
T. Rettberg and S. Spors, ?Time-Domain Behaviour of Spherical Microphone
Arrays at High Orders,? in Proceedings of the 40th DAGA, 2014, pp. 598?599.
B. Bernschutz, C. Pörschmann, S. Spors, and S. Weinzierl, ?Soft-Limiting
der modalen Amplitudenverstärkung bei sphärischen Mikrofonarrays im
Plane Wave Decomposition Verfahren,?
in Proceedings of the 37th DAGA, 2011, pp. 661?662.
Boaz Rafaely: Fundamentals of spherical array processing. In. Springer
topics in signal processing. Benesty, J.; Kellermann, W. (Eds.),
Springer, Heidelberg et al. (2015).
Benjamin Bernschütz: Microphone Arrays and Sound Field Decomposition
for Dynamic Binaural Recording. Ph.D. dissertation, Technical University
Berlin (2016).
(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 % function HRIRs_sfd_ext = supdeq_rangeExt(HRIRs_sfd, r1, regMethod, ampLimit, compTimeShift, compEnergy, rmin, r0) 0004 % 0005 % This function performs range extrapolation of HRTFs in SH-domain. This 0006 % way, near-field HRTFs can be obtained from far-field HRTFs. In brief, 0007 % radial propagation functions are applied to the SH-coefficients to 0008 % calculate the range-extrapolated SH-coefficients for the new distance r1.? 0009 % As the radial functions (in case of an exterior problem the hankel 0010 % quotient) increase dramatically towards lower kr, a regularization method 0011 % like for example soft-limiting or the Tikhonov regularization needs to be 0012 % applied. See references for more information about range extrapolation 0013 % and regularization approaches. 0014 % 0015 % Output: 0016 % HRIRs_sfd_ext - Struct with Spherical Harmonic coefficients 0017 % (SH-Coefficients) for the left (Hl_nm) and right (Hr_nm) 0018 % channel/ear of the range-extrapolated HRTFs, absolute 0019 % frequency scale f, transform order N, FFToversize, 0020 % (new range-extrapolated) source distance r1, 0021 % non-regularized radial functions (hankel quotient), 0022 % regularized (and applied) hankel quotient, and info 0023 % about applied regularization and range extrapolation 0024 % 0025 % Input: 0026 % HRIRs_sfd - Struct with Spherical Harmonic coefficients 0027 % (SH-Coefficients) for the left (Hl_nm) and right (Hr_nm) 0028 % channel/ear, absolute frequency scale f, 0029 % transform order N, FFToversize, and source distance r0 0030 % r1 - New range-extrapolated source distance r1 in m 0031 % regMethod - Regularization method (string parameter) for the 0032 % radial function (the hankel quotient). The following 0033 % methods are available: 0034 % 'none' - no regularization 0035 % 'disc' - simply discard contributions where the 0036 % radial filter magnitude exceeds ampLimit 0037 % 'disc_kr' - discard contributions at higher orders 0038 % based on kr_min. Results are pretty 0039 % similar to normal 'disc' 0040 % 'hardLimit'- set a hard limit for all radial filter 0041 % components at ampLimit instead of 0042 % dropping them 0043 % 'softLimit'- soft-limit the components 0044 % 'Tikh' - apply Tikhonov regularization 0045 % See Bernschütz2011 and Rettberg2014 for an overview 0046 % of most regularization methods. Implementation of most 0047 % of the regularization methods according to Rettberg2014 0048 % Regularization 'disc_kr' implemented according to 0049 % Duraiswami2004 and Pollow2012 0050 % Default: 'softLimit' 0051 % ampLimit - Amplitude limit of radial functions in dB. Not needed 0052 % for all regularization methods, but since the default 0053 % regularization is 'softLimiting', the amplitude limit 0054 % can be freely chosen here. 0055 % Default: 20 dB 0056 % compTimeShift - true/false - compensate time shift induced by 0057 % range-extrapolation with an inverse phase term 0058 % Default: true 0059 % compEnergy - true/false - compensate the in-/decrease in energy 0060 % averaged over all orders and modes 0061 % (single value compensation) 0062 % Default: true 0063 % rmin - Smallest radius in m, meaning the radius of the 0064 % smallest sphere enclosing the listener?s head 0065 % with all its significant scattering sources. This 0066 % defines the truncation order N = [krmin] when 0067 % regMethod 'disc_kr' according to Duraiswami2004 and 0068 % Pollow2012 is used. Otherwise, rmin has no influence. 0069 % Default: 0.1 0070 % r0 - Actual source distance in m input HRIRs_sfd. 0071 % The HRIRs_sfd structs have a field with the actual 0072 % source distance and the function will use the value 0073 % in this field as r0. However, if this field is 0074 % missing, r0 can be given as an argument. 0075 % Default: [] (use value in HRIRs_sfd.sourceDistance) 0076 % 0077 % Dependencies: AKtools 0078 % 0079 % References: 0080 % R. Duraiswami, D. N. Zotkin, and N. A. Gumerov, ?Interpolation and range 0081 % extrapolation of HRTFs,? in Proceedings of the IEEE International 0082 % Conference on Acoustics, Speech, and Signal Processing, 2004, pp. IV45-IV48. 0083 % 0084 % M. Pollow, K.-V. Nguyen, O. Warusfel, T. Carpentier, M. Müller-Trapet, 0085 % M. Vorländer, and M. Noisternig, ?Calculation of Head-Related Transfer 0086 % Functions for Arbitrary Field Points Using Spherical Harmonics 0087 % Decomposition,? Acta Acust. united Ac., vol. 98, no. 1, pp. 72?82, 2012. 0088 % 0089 % T. Rettberg and S. Spors, ?Time-Domain Behaviour of Spherical Microphone 0090 % Arrays at High Orders,? in Proceedings of the 40th DAGA, 2014, pp. 598?599. 0091 % 0092 % B. Bernschutz, C. Pörschmann, S. Spors, and S. Weinzierl, ?Soft-Limiting 0093 % der modalen Amplitudenverstärkung bei sphärischen Mikrofonarrays im 0094 % Plane Wave Decomposition Verfahren,? 0095 % in Proceedings of the 37th DAGA, 2011, pp. 661?662. 0096 % 0097 % Boaz Rafaely: Fundamentals of spherical array processing. In. Springer 0098 % topics in signal processing. Benesty, J.; Kellermann, W. (Eds.), 0099 % Springer, Heidelberg et al. (2015). 0100 % 0101 % Benjamin Bernschütz: Microphone Arrays and Sound Field Decomposition 0102 % for Dynamic Binaural Recording. Ph.D. dissertation, Technical University 0103 % Berlin (2016). 0104 % 0105 % (C) 2018 by JMA, Johannes M. Arend 0106 % TH Köln - University of Applied Sciences 0107 % Institute of Communications Engineering 0108 % Department of Acoustics and Audio Signal Processing 0109 0110 function HRIRs_sfd_ext = supdeq_rangeExt(HRIRs_sfd, r1, regMethod, ampLimit, compTimeShift, compEnergy, rmin, r0) 0111 0112 if nargin < 1 || isempty(HRIRs_sfd) 0113 error('Please specifiy HRIRs_sfd!'); 0114 end 0115 0116 if nargin < 2 || isempty(r1) 0117 error('Please specifiy r1!'); 0118 end 0119 0120 if nargin < 3 || isempty(regMethod) 0121 regMethod = 'softLimit'; 0122 end 0123 0124 if nargin < 4 || isempty(ampLimit) 0125 ampLimit = 20; 0126 end 0127 0128 if nargin < 5 || isempty(compTimeShift) 0129 compTimeShift = 1; 0130 end 0131 0132 if nargin < 6 || isempty(compEnergy) 0133 compEnergy = 1; 0134 end 0135 0136 if nargin < 7 || isempty(rmin) 0137 rmin = 0.1; 0138 end 0139 0140 if nargin < 7 || isempty(rmin) 0141 rmin = 0.1; 0142 end 0143 0144 if nargin < 8 || isempty(r0) 0145 if isfield(HRIRs_sfd,'sourceDistance') 0146 r0 = HRIRs_sfd.sourceDistance; 0147 else 0148 error('Please specifiy r0!'); 0149 end 0150 end 0151 0152 %% Get kr and N 0153 0154 f = HRIRs_sfd.f; 0155 c = 343; %Set speed of sound as default with 343 m/s 0156 k = 2*pi*f/c; 0157 kr0 = k*r0; 0158 kr1 = k*r1; 0159 N = HRIRs_sfd.N; 0160 0161 %% Get hankel quotient for range extrapolation 0162 0163 %Get hankel functions 0164 hn_r0 = zeros(N+1,length(f)); 0165 hn_r1 = zeros(N+1,length(f)); 0166 for kk = 0:N 0167 hn_r0(kk+1,:) = AKshRadial(kr0,'hankel',2,kk); 0168 hn_r1(kk+1,:) = AKshRadial(kr1,'hankel',2,kk); 0169 end 0170 0171 %Set 0 Hz bin to 1 0172 hn_r0(:,1) = 1; 0173 hn_r1(:,1) = 1; 0174 0175 %Calculate hankel quotient 0176 hn_ext = hn_r1./hn_r0; 0177 0178 %Compensate time shift if intended 0179 if compTimeShift 0180 t = ((r0-r1)/c); 0181 timeShift = exp((0-1i) * complex((2*pi*f*t)+0i)); 0182 for kk = 0:N 0183 hn_ext(kk+1,:) = hn_ext(kk+1,:).*timeShift; 0184 end 0185 end 0186 hn_ext_noReg = hn_ext; 0187 0188 %% Apply regularization to hankel quotient 0189 0190 %Implementation according to Rettberg2014 0191 %https://github.com/spatialaudio/sfa-numpy/blob/master/micarray/modal/radial.py 0192 %Function "regularize" 0193 0194 if ~strcmp(regMethod,'disc_kr') %Regularization disc_kr implemented in a different way... 0195 a0_lin = 10^(ampLimit/20); 0196 %Get values where abs(hn_ext) > maxGain 0197 idx = abs(hn_ext) > a0_lin; 0198 0199 %Get regularisation function hn_reg depending on method 0200 if strcmp(regMethod,'none') 0201 hn_reg = ones(size(hn_ext,1),size(hn_ext,2)); 0202 0203 elseif strcmp(regMethod,'disc') 0204 hn_reg = ones(size(hn_ext,1),size(hn_ext,2)); 0205 hn_reg(idx) = 0; 0206 0207 elseif strcmp(regMethod,'hardLimit') 0208 hn_reg = ones(size(hn_ext,1),size(hn_ext,2)); 0209 hn_reg(idx) = a0_lin ./ abs(hn_ext(idx)); 0210 0211 elseif strcmp(regMethod,'softLimit') 0212 scaling = pi/2; 0213 hn_reg = a0_lin ./ abs(hn_ext); 0214 hn_reg = 2/pi * atan(scaling*hn_reg); 0215 0216 elseif strcmp(regMethod,'Tikh') 0217 a0_lin = sqrt(a0_lin/2); 0218 alpha = (1 - sqrt(1 - 1/(a0_lin^2))) / (1 + sqrt(1 - 1/(a0_lin^2))); 0219 hn_reg = 1 ./ (1 + alpha^2 * abs(hn_ext).^2); 0220 end 0221 0222 %Apply hn_reg to hn_ext 0223 hn_ext = hn_ext .* hn_reg; 0224 end 0225 0226 0227 %% Get hankel quotient matrix with size of SH-coefficients matrix 0228 0229 %Set nan values at kr = 0 to 0 0230 hn_ext(isnan(hn_ext)) = 0; 0231 0232 hn_ext_matrix = zeros(size(HRIRs_sfd.Hl_nm,1),size(HRIRs_sfd.Hl_nm,2)); 0233 modeCounter = 1; 0234 for n = 1:N+1 0235 %Get "corner"-modes 0236 modeLeft = n^2-modeCounter+1; 0237 modeRight = n^2; 0238 0239 %Define modes vector 0240 modes = modeLeft:modeRight; 0241 0242 for mcr = 1:length(modes) 0243 hn_ext_matrix(modes(mcr),:) = hn_ext(n,:); 0244 end 0245 modeCounter = modeCounter + 2; 0246 end 0247 0248 %% Apply hankel quotient matrix to SH-coefficients 0249 0250 Hl_nm_r0 = HRIRs_sfd.Hl_nm; 0251 Hr_nm_r0 = HRIRs_sfd.Hr_nm; 0252 Hl_nm_r1 = zeros(size(Hl_nm_r0,1),size(Hl_nm_r0,2)); 0253 Hr_nm_r1 = zeros(size(Hr_nm_r0,1),size(Hr_nm_r0,2)); 0254 0255 if strcmp(regMethod,'disc_kr') %With disc_kr regularization according to Duraiswami2004 0256 krmin = k*rmin; 0257 Nh = round(krmin)+1; %Duraiswami2004, Eq.8 0258 Nh(Nh>N) = N; 0259 0260 for bin = 1:size(Hl_nm_r0,2) 0261 Nbin = Nh(bin); 0262 range = (Nbin+1)^2; 0263 Hl_nm_r1(1:range,bin) = Hl_nm_r0(1:range,bin) .* hn_ext_matrix(1:range,bin); 0264 Hr_nm_r1(1:range,bin) = Hr_nm_r0(1:range,bin) .* hn_ext_matrix(1:range,bin); 0265 end 0266 0267 else %With any other regularization already applied above to hn_ext 0268 0269 for kk = 1:size(hn_ext_matrix,1) 0270 Hl_nm_r1(kk,:) = Hl_nm_r0(kk,:) .* hn_ext_matrix(kk,:); 0271 Hr_nm_r1(kk,:) = Hr_nm_r0(kk,:) .* hn_ext_matrix(kk,:); 0272 end 0273 end 0274 0275 %% Apply energy compensation if intended 0276 0277 if compEnergy 0278 %Get average energy over all orders and modes of SH-coefficients 0279 [~,eL_r0] = AKshEnergy(Hl_nm_r0); 0280 [~,eR_r0] = AKshEnergy(Hr_nm_r0); 0281 [~,eL_r1] = AKshEnergy(Hl_nm_r1); 0282 [~,eR_r1] = AKshEnergy(Hr_nm_r1); 0283 0284 %Apply energy-compensation factor to SH-coefficients 0285 Hl_nm_r1 = Hl_nm_r1 * sqrt(eL_r0/eL_r1); 0286 Hr_nm_r1 = Hr_nm_r1 * sqrt(eR_r0/eR_r1); 0287 end 0288 0289 %% Write new output struct with additional info on range extrapolation 0290 0291 HRIRs_sfd_ext = HRIRs_sfd; 0292 HRIRs_sfd_ext.Hl_nm = Hl_nm_r1; 0293 HRIRs_sfd_ext.Hr_nm = Hr_nm_r1; 0294 HRIRs_sfd_ext.sourceDistance = r1; 0295 HRIRs_sfd_ext.rangeExt.r0 = r0; 0296 HRIRs_sfd_ext.rangeExt.r1 = r1; 0297 HRIRs_sfd_ext.rangeExt.regMethod = regMethod; 0298 HRIRs_sfd_ext.rangeExt.ampLimit = ampLimit; 0299 HRIRs_sfd_ext.rangeExt.compTimeShift = compTimeShift; 0300 HRIRs_sfd_ext.rangeExt.compEnergy = compEnergy; 0301 HRIRs_sfd_ext.rangeExt.rmin = rmin; 0302 HRIRs_sfd_ext.rangeExt.hn_ext = hn_ext; 0303 HRIRs_sfd_ext.rangeExt.hn_ext_noReg = hn_ext_noReg; 0304 0305 end 0306