% SUpDEq - Spatial Upsampling by Directional Equalization
function [shiftedHRTF_L, shiftedHRTF_R, ffShiftedHRTFdataset] = supdeq_dvf(HRTFdataset, newDistance, samplingGrid ,shiftToFarField, adjustAmplitude, radius, earPosition, N_HRTF)
This function applies a distance shift to HRTFs using distance variation
functions (DVFs). As it works with HRTFs in spherical harmonics domain,
acoustic parallax effects (high frequency parallax effects induced by
the pinna which are not covered by the rigid sphere transfer functions)
can be embedded into the distance shifted HRTFs.
Output:
shiftedHRTF_L/R - Shifted HRTF_L/R in frequency domain (single sided complex spectrum)
ffShiftedHRTFdataset - Struct with far-field shifted HRTF dataset in SH
domain according to input HRTFdataset. Only if
'shiftToFarField' was applied
Input:
HRTFdataset - Struct with the HRTF dataset as
SH-coefficients for the left (Hl_nm) and
right (Hr_nm) channel/ear, absolute frequency scale f,
transform order N, FFToversize and sourceDistance
If field 'sourceDistance' does not exist, the
function assumes far-field measurements and option
'shiftToFarField' will always be set to 'false'
newDistance - New distance of HRTF dataset in m
samplingGrid - Spatial sampling grid (Q x 2 matrix) for the distance shifted HRTFs,
where the first column holds the azimuth and the second
the elevation (both in degree).
Azimuth in degree (0=front, 90=left, 180=back, 270=right)
(0 points to positive x-axis, 90 to positive y-axis)
Elevations in degree (0=North Pole, 90=front, 180=South Pole)
(0 points to positive z-axis, 180 to negative z-axis)
shiftToFarField - Boolean with true / false. Option to shift
measurements, which strictly speaking were not measured in
the far-field (e.g., at d = 3.25 m as our Neumann KU100 set)
and thus still have acoustic parallax effects, to far field before
shifting to new distance.
Default: false
Note: If set to false, only the inverse SH
transform will be applied to get HRTFs for the
new distance according to the sampling grid.
Thus, processing will be much faster. If set to
true, various forward and inverse SH transform
are gonne be performed on a high order grid, so
processing takes longer.
adjustAmplitude - Boolen with true / false. Option to adjust the HRTF level
according to the 1/r law.
Default: false --> HRTF have roughly the same
amplitude after distance shift
radius - Define radius of rigid sphere / HRTFs in m
Default: 0.0875
earPosition - 4 x 1 row vector describing the position of the ears in
spherical coordinates in degree [azL, elL, azR, elR]
Default: [90, 90, 270, 90] (left-right symmetrical)
N_HRTF - Spherical harmonics order for HRTF processing. If
not specified, Nmax from the HRTFdataset struct will be applied (mostly
the best choice)
Default: Nmax from HRTFdataset struct
Dependencies: -
References:
[1] A. Kan, C. Jin, and A. van Schaik, “A psychophysical evaluation of
near-field head-related transfer functions synthesized using a distance
variation function,” J. Acoust. Soc. Am., vol. 125, no. 4, pp. 2233–2242, 2009.
[2] S. Spagnol, E. Tavazzi, and F. Avanzini, “Distance rendering and
perception of nearby virtual sound sources with a near-field filter model,”
Appl. Acoust., vol. 115, pp. 61–73, 2017.
[3] D. Romblom and B. Cook, “Near-Field Compensation for HRTF Processing,”
in Proceedings of the 125th AES Convention, San Francisco, USA, 2008, pp. 1–6.
(C) 2020 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 [shiftedHRTF_L, shiftedHRTF_R, ffShiftedHRTFdataset] = supdeq_dvf(HRTFdataset, newDistance, samplingGrid ,shiftToFarField, adjustAmplitude, radius, earPosition, N_HRTF) 0004 % 0005 % This function applies a distance shift to HRTFs using distance variation 0006 % functions (DVFs). As it works with HRTFs in spherical harmonics domain, 0007 % acoustic parallax effects (high frequency parallax effects induced by 0008 % the pinna which are not covered by the rigid sphere transfer functions) 0009 % can be embedded into the distance shifted HRTFs. 0010 % 0011 % Output: 0012 % shiftedHRTF_L/R - Shifted HRTF_L/R in frequency domain (single sided complex spectrum) 0013 % ffShiftedHRTFdataset - Struct with far-field shifted HRTF dataset in SH 0014 % domain according to input HRTFdataset. Only if 0015 % 'shiftToFarField' was applied 0016 % 0017 % Input: 0018 % HRTFdataset - Struct with the HRTF dataset as 0019 % SH-coefficients for the left (Hl_nm) and 0020 % right (Hr_nm) channel/ear, absolute frequency scale f, 0021 % transform order N, FFToversize and sourceDistance 0022 % If field 'sourceDistance' does not exist, the 0023 % function assumes far-field measurements and option 0024 % 'shiftToFarField' will always be set to 'false' 0025 % newDistance - New distance of HRTF dataset in m 0026 % samplingGrid - Spatial sampling grid (Q x 2 matrix) for the distance shifted HRTFs, 0027 % where the first column holds the azimuth and the second 0028 % the elevation (both in degree). 0029 % Azimuth in degree (0=front, 90=left, 180=back, 270=right) 0030 % (0 points to positive x-axis, 90 to positive y-axis) 0031 % Elevations in degree (0=North Pole, 90=front, 180=South Pole) 0032 % (0 points to positive z-axis, 180 to negative z-axis) 0033 % shiftToFarField - Boolean with true / false. Option to shift 0034 % measurements, which strictly speaking were not measured in 0035 % the far-field (e.g., at d = 3.25 m as our Neumann KU100 set) 0036 % and thus still have acoustic parallax effects, to far field before 0037 % shifting to new distance. 0038 % Default: false 0039 % Note: If set to false, only the inverse SH 0040 % transform will be applied to get HRTFs for the 0041 % new distance according to the sampling grid. 0042 % Thus, processing will be much faster. If set to 0043 % true, various forward and inverse SH transform 0044 % are gonne be performed on a high order grid, so 0045 % processing takes longer. 0046 % adjustAmplitude - Boolen with true / false. Option to adjust the HRTF level 0047 % according to the 1/r law. 0048 % Default: false --> HRTF have roughly the same 0049 % amplitude after distance shift 0050 % radius - Define radius of rigid sphere / HRTFs in m 0051 % Default: 0.0875 0052 % earPosition - 4 x 1 row vector describing the position of the ears in 0053 % spherical coordinates in degree [azL, elL, azR, elR] 0054 % Default: [90, 90, 270, 90] (left-right symmetrical) 0055 % N_HRTF - Spherical harmonics order for HRTF processing. If 0056 % not specified, Nmax from the HRTFdataset struct will be applied (mostly 0057 % the best choice) 0058 % Default: Nmax from HRTFdataset struct 0059 % 0060 % Dependencies: - 0061 % 0062 % References: 0063 % [1] A. Kan, C. Jin, and A. van Schaik, “A psychophysical evaluation of 0064 % near-field head-related transfer functions synthesized using a distance 0065 % variation function,” J. Acoust. Soc. Am., vol. 125, no. 4, pp. 2233–2242, 2009. 0066 % 0067 % [2] S. Spagnol, E. Tavazzi, and F. Avanzini, “Distance rendering and 0068 % perception of nearby virtual sound sources with a near-field filter model,” 0069 % Appl. Acoust., vol. 115, pp. 61–73, 2017. 0070 % 0071 % [3] D. Romblom and B. Cook, “Near-Field Compensation for HRTF Processing,” 0072 % in Proceedings of the 125th AES Convention, San Francisco, USA, 2008, pp. 1–6. 0073 % 0074 % (C) 2020 by JMA, Johannes M. Arend 0075 % TH Köln - University of Applied Sciences 0076 % Institute of Communications Engineering 0077 % Department of Acoustics and Audio Signal Processing 0078 0079 function [shiftedHRTF_L, shiftedHRTF_R, ffShiftedHRTFdataset] = supdeq_dvf(HRTFdataset, newDistance, samplingGrid, shiftToFarField, adjustAmplitude, radius, earPosition, N_HRTF) 0080 0081 if nargin < 2 || isempty(newDistance) 0082 error('Please specify new distance in m'); 0083 end 0084 0085 if nargin < 3 || isempty(samplingGrid) 0086 error('Please define sampling grid'); 0087 end 0088 0089 if nargin < 4 || isempty(shiftToFarField) 0090 shiftToFarField = false; 0091 end 0092 0093 if nargin < 5 || isempty(adjustAmplitude) 0094 adjustAmplitude = false; 0095 end 0096 0097 if nargin < 6 || isempty(radius) 0098 radius = 0.0875; 0099 end 0100 0101 if nargin < 7 || isempty(earPosition) 0102 earPosition = [90, 90, 270, 90]; 0103 end 0104 0105 if nargin < 8 || isempty(N_HRTF) 0106 if isfield(HRTFdataset,'N') 0107 N_HRTF = HRTFdataset.N; 0108 end 0109 if isfield(HRTFdataset,'Nmax') 0110 N_HRTF = HRTFdataset.Nmax; 0111 end 0112 end 0113 0114 if isfield(HRTFdataset,'sourceDistance') 0115 md = HRTFdataset.sourceDistance; 0116 else 0117 md = 100; %Assume far field 0118 shiftToFarField = false; 0119 end 0120 0121 NFFT = length(HRTFdataset.f)*2-2; 0122 fs = HRTFdataset.f(end)*2; 0123 ffShiftedHRTFdataset = nan; 0124 0125 %% Get various grids 0126 0127 %Omit weights of samplingGrid if defined 0128 samplingGrid = samplingGrid(:,1:2); 0129 0130 %Get sg for STF 0131 sgSTF = samplingGrid; sgSTF(:,2) = 90-sgSTF(:,2); 0132 sgSTF_nd = sgSTF; sgSTF_nd(:,3) = newDistance; 0133 sgSTF_md = sgSTF; sgSTF_md(:,3) = md; 0134 ffDistance = 100; 0135 sgSTF_ff = sgSTF; sgSTF_ff(:,3) = ffDistance; 0136 0137 earPositionSTF = earPosition; earPositionSTF(2) = 90-earPositionSTF(2); earPositionSTF(4) = 90-earPositionSTF(4); 0138 0139 %Get parallax grids of new distance 0140 [sg_nd_para_l,sg_nd_para_r] = supdeq_parallax(samplingGrid,newDistance,radius,earPosition); 0141 0142 %% Get rigid sphere transfer functions for new distance 0143 0144 stf_nd = AKsphericalHead(sgSTF_nd,earPositionSTF,false,radius,newDistance,100,NFFT,fs); 0145 stf_nd_L = AKboth2singleSidedSpectrum(fft(squeeze(stf_nd(:,:,1)))); 0146 stf_nd_R = AKboth2singleSidedSpectrum(fft(squeeze(stf_nd(:,:,2)))); 0147 0148 0149 %% If dataset should be shifted to far field first before shifting to new distance 0150 %This compensates parallax effects which are already in the measurement 0151 0152 if shiftToFarField 0153 0154 disp('Shifting measured HRTFs first to far field'); 0155 0156 %Get Lebedev sampling grid according to N 0157 lebN = [1,2,3,4,5,6,7,8,9,10,11,12,13,14,15,17,20,23,26,29,32,35,38,41,44,47,50,53,56,59,62,65]; 0158 0159 %Get closest N 0160 [~,cId] = min(abs(N_HRTF-lebN)); 0161 N_HRTF_leb = lebN(cId); 0162 0163 if N_HRTF_leb < N_HRTF 0164 cId = cId + 1; 0165 N_HRTF = lebN(cId); 0166 else 0167 N_HRTF = N_HRTF_leb; 0168 end 0169 0170 %Get sampling grid according to N_HRTF for internal sh processing 0171 sg_int = supdeq_lebedev([],N_HRTF); 0172 %Omit weights 0173 sg_int = sg_int(:,1:2); 0174 0175 %Get parallax grid for measurement distance 0176 [sg_md_para_l,sg_md_para_r] = supdeq_parallax(sg_int,md,radius,earPosition); 0177 0178 %Get STFs for far field position on dense grid 0179 ffDistance = 100; 0180 sgSTF_int_ff = sg_int; sgSTF_int_ff(:,2) = 90-sgSTF_int_ff(:,2); 0181 sgSTF_int_ff(:,3) = ffDistance; 0182 stf_int_ff = AKsphericalHead(sgSTF_int_ff,earPositionSTF,false,radius,ffDistance,100,NFFT,fs); 0183 stf_int_ff_L= AKboth2singleSidedSpectrum(fft(squeeze(stf_int_ff(:,:,1)))); 0184 stf_int_ff_R= AKboth2singleSidedSpectrum(fft(squeeze(stf_int_ff(:,:,2)))); 0185 0186 %Get STFs for measurement position on dense grid 0187 sgSTF_int_md = sg_int; sgSTF_int_md(:,2) = 90-sgSTF_int_md(:,2); 0188 sgSTF_int_md(:,3) = md; 0189 stf_int_md = AKsphericalHead(sgSTF_int_md,earPositionSTF,false,radius,md,100,NFFT,fs); 0190 stf_int_md_L= AKboth2singleSidedSpectrum(fft(squeeze(stf_int_md(:,:,1)))); 0191 stf_int_md_R= AKboth2singleSidedSpectrum(fft(squeeze(stf_int_md(:,:,2)))); 0192 0193 %Calculate DVFs for left and right ear, Target = far 0194 DVF_int_L = stf_int_ff_L./stf_int_md_L; 0195 DVF_int_R = stf_int_ff_R./stf_int_md_R; 0196 0197 %Get parallaxe grid HRTFs of measured dataset for measurement position 0198 [hrtf_L] = supdeq_getArbHRTF(HRTFdataset,sg_md_para_l,'DEG',0,'ak'); 0199 [~,hrtf_R] = supdeq_getArbHRTF(HRTFdataset,sg_md_para_r,'DEG',1,'ak'); 0200 0201 %Apply DVF to shift to far-field 0202 ff_hrtf_L = (DVF_int_L.').*hrtf_L; 0203 ff_hrtf_R = (DVF_int_R.').*hrtf_R; 0204 0205 %Transform far-field shifted HRTF set to SH domain again for further 0206 %processing 0207 ffShiftedHRTFdataset = supdeq_hrtf2sfd(ff_hrtf_L,ff_hrtf_R,N_HRTF,sg_int,fs,'ak'); 0208 0209 %Get some parameters 0210 ffShiftedHRTFdataset.FFToversize = HRTFdataset.FFToversize; 0211 ffShiftedHRTFdataset.sourceDistance = 100; 0212 ffShiftedHRTFdataset.ffShifted = true; 0213 end 0214 0215 %% Apply DVFs 0216 0217 disp('Shifting HRTFs to new position'); 0218 0219 if shiftToFarField 0220 %Calculate DVFs for left and right ear, Target = close 0221 %Work with far-field distance STFs 0222 stf_ff = AKsphericalHead(sgSTF_ff,earPositionSTF,false,radius,ffDistance,100,NFFT,fs); 0223 stf_ff_L = AKboth2singleSidedSpectrum(fft(squeeze(stf_ff(:,:,1)))); 0224 stf_ff_R = AKboth2singleSidedSpectrum(fft(squeeze(stf_ff(:,:,2)))); 0225 0226 DVF_L = stf_nd_L./stf_ff_L; 0227 DVF_R = stf_nd_R./stf_ff_R; 0228 0229 %Get parallaxe grid HRTFs of far-field dataset for new position 0230 [hrtf_L] = supdeq_getArbHRTF(ffShiftedHRTFdataset,sg_nd_para_l,'DEG',0,'ak'); 0231 [~,hrtf_R] = supdeq_getArbHRTF(ffShiftedHRTFdataset,sg_nd_para_r,'DEG',1,'ak'); 0232 0233 else 0234 %Calculate DVFs for left and right ear, Target = close 0235 %Work with measurement distance STFs 0236 stf_md = AKsphericalHead(sgSTF_md,earPositionSTF,false,radius,md,100,NFFT,fs); 0237 stf_md_L = AKboth2singleSidedSpectrum(fft(squeeze(stf_md(:,:,1)))); 0238 stf_md_R = AKboth2singleSidedSpectrum(fft(squeeze(stf_md(:,:,2)))); 0239 0240 DVF_L = stf_nd_L./stf_md_L; 0241 DVF_R = stf_nd_R./stf_md_R; 0242 0243 %Get parallaxe grid HRTFs of far-field dataset for new position 0244 [hrtf_L] = supdeq_getArbHRTF(HRTFdataset,sg_nd_para_l,'DEG',0,'ak'); 0245 [~,hrtf_R] = supdeq_getArbHRTF(HRTFdataset,sg_nd_para_r,'DEG',1,'ak'); 0246 end 0247 0248 %Apply DVFs 0249 shiftedHRTF_L = (DVF_L.').*hrtf_L; 0250 shiftedHRTF_R = (DVF_R.').*hrtf_R; 0251 0252 if adjustAmplitude 0253 0254 %Set referenceDistance to 1 if md was not defined 0255 if md == 100 0256 refDistance = 1; 0257 else 0258 refDistance = md; 0259 end 0260 0261 shiftedHRTF_L = shiftedHRTF_L*(refDistance/newDistance); 0262 shiftedHRTF_R = shiftedHRTF_R*(refDistance/newDistance); 0263 0264 end 0265 0266 disp('Done...'); 0267 0268 end 0269