% SUpDEq - Spatial Upsampling by Directional Equalization
function eqDataset = supdeq_getEqDataset(N, earDistance, NFFT, fs, waveType, sourceDistance, earPosition, sphereOffset, transformCore, swSettings)
This function returns the sound pressure distribution of an incident
plane wave for different sound incidence directions on the left ear
position and on the right ear position of a sphere which models the
human head.
Output:
eqDataset - Struct with SH-coefficients describing the sound
incidence at the letf/right ear position on the sphere
(Hl_nm/Hr_nm). Length of the SH-coefficients is NFFT/2+1
(single sided spectrum). Important input variables are
stored in the struct.
Input:
N - Transform order N
Default: 35
earDistance - Distance between both ear positions in m (radius*2)
Default: 0.165
NFFT - FFT size of the SH-coefficients to be returned
Default: 512
fs - Sampling rate
Default: 48000
waveType - Boolean for plane wave (0) or spherical wave (1)
Default: 0
sourceDistance- Distance (in meter) between source and sphere.
Only considered if waveType = 1 (spherical wave)
Default: 1
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)
sphereOffset - 3 x 1 row vector describing the offset of the sphere from
the origin in m [xd, yd, zd], with x shift on "depth-axis"
of the sphere, y shift on "width-axis" of the sphere,
and z shift on "height-axis" of the sphere. See matlab
function "cart2sph" for describtion of the coordinate.
Positive values
system.
Default: [] - No offset
Combinations of earPosition and sphereOffset are stored
in the field "earPosition" in the output struct
transformCore - String to define method to be used for the rigid spher
transfer function synthesis
'sofia - sofia_wgc from SOFiA toolbox
'ak' - AKsphericalHead from AKtools
The results are almost exactly the same
Default: 'sofia'
swSettings - Settings only for sourceType = 1 (spherical wave),
indicated by integer 0,1,2,or 3.
Reference distance of spherical wave is set to 1.00m for
both transform cores!
0 - Compensate time shift according to distance to
reference distance (1.00m), but maintain level changes.
1 - Compensate time shift and compensate level changes
according to 1/r law (to level at r0 = 1m). The 1/r
compensation is not necessarily accurate enough
in the near-field!
2 - Maintain time shift but compensate level changes
3 - Maintain time shift and maintain level changes (no
compensation)
Default: 1
IMPORTANT NOTE: ALWAYS USE THE SAME TRANSFORM-CORE FOR EQUALIZATION AND
DE-EQUALIZATION. AS THE LEVEL AND PHASE PROPERTIES OF THE MODELS ARE
DIFFERENT, MIXING THEM UP RESULTS IN INCORRECT RESULTS!
Dependencies: SOFiA toolbox
References:
Benjamin Bernschütz: Microphone Arrays and Sound Field Decomposition
for Dynamic Binaural Recording. Ph.D. dissertation, Technical University
Berlin (2016).
R. O. Duda and W. L. Martens, ?Range dependence of the response of a
spherical head model,? J. Acoust. Soc. Am., vol. 104, no. 5,
pp. 3048?3058, 1998.
(C) 2018 by JMA, Johannes M. Arend
CP, Christoph Pörschmann
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 eqDataset = supdeq_getEqDataset(N, earDistance, NFFT, fs, waveType, sourceDistance, earPosition, sphereOffset, transformCore, swSettings) 0004 % 0005 % This function returns the sound pressure distribution of an incident 0006 % plane wave for different sound incidence directions on the left ear 0007 % position and on the right ear position of a sphere which models the 0008 % human head. 0009 % 0010 % Output: 0011 % eqDataset - Struct with SH-coefficients describing the sound 0012 % incidence at the letf/right ear position on the sphere 0013 % (Hl_nm/Hr_nm). Length of the SH-coefficients is NFFT/2+1 0014 % (single sided spectrum). Important input variables are 0015 % stored in the struct. 0016 % 0017 % Input: 0018 % N - Transform order N 0019 % Default: 35 0020 % earDistance - Distance between both ear positions in m (radius*2) 0021 % Default: 0.165 0022 % NFFT - FFT size of the SH-coefficients to be returned 0023 % Default: 512 0024 % fs - Sampling rate 0025 % Default: 48000 0026 % waveType - Boolean for plane wave (0) or spherical wave (1) 0027 % Default: 0 0028 % sourceDistance- Distance (in meter) between source and sphere. 0029 % Only considered if waveType = 1 (spherical wave) 0030 % Default: 1 0031 % earPosition - 4 x 1 row vector describing the position of the ears in 0032 % spherical coordinates in degree [azL, elL, azR, elR] 0033 % Default: [90, 90, 270, 90] (left-right symmetrical) 0034 % sphereOffset - 3 x 1 row vector describing the offset of the sphere from 0035 % the origin in m [xd, yd, zd], with x shift on "depth-axis" 0036 % of the sphere, y shift on "width-axis" of the sphere, 0037 % and z shift on "height-axis" of the sphere. See matlab 0038 % function "cart2sph" for describtion of the coordinate. 0039 % Positive values 0040 % system. 0041 % Default: [] - No offset 0042 % Combinations of earPosition and sphereOffset are stored 0043 % in the field "earPosition" in the output struct 0044 % transformCore - String to define method to be used for the rigid spher 0045 % transfer function synthesis 0046 % 'sofia - sofia_wgc from SOFiA toolbox 0047 % 'ak' - AKsphericalHead from AKtools 0048 % The results are almost exactly the same 0049 % Default: 'sofia' 0050 % swSettings - Settings only for sourceType = 1 (spherical wave), 0051 % indicated by integer 0,1,2,or 3. 0052 % Reference distance of spherical wave is set to 1.00m for 0053 % both transform cores! 0054 % 0 - Compensate time shift according to distance to 0055 % reference distance (1.00m), but maintain level changes. 0056 % 1 - Compensate time shift and compensate level changes 0057 % according to 1/r law (to level at r0 = 1m). The 1/r 0058 % compensation is not necessarily accurate enough 0059 % in the near-field! 0060 % 2 - Maintain time shift but compensate level changes 0061 % 3 - Maintain time shift and maintain level changes (no 0062 % compensation) 0063 % Default: 1 0064 % 0065 % IMPORTANT NOTE: ALWAYS USE THE SAME TRANSFORM-CORE FOR EQUALIZATION AND 0066 % DE-EQUALIZATION. AS THE LEVEL AND PHASE PROPERTIES OF THE MODELS ARE 0067 % DIFFERENT, MIXING THEM UP RESULTS IN INCORRECT RESULTS! 0068 % 0069 % Dependencies: SOFiA toolbox 0070 % 0071 % References: 0072 % Benjamin Bernschütz: Microphone Arrays and Sound Field Decomposition 0073 % for Dynamic Binaural Recording. Ph.D. dissertation, Technical University 0074 % Berlin (2016). 0075 % 0076 % R. O. Duda and W. L. Martens, ?Range dependence of the response of a 0077 % spherical head model,? J. Acoust. Soc. Am., vol. 104, no. 5, 0078 % pp. 3048?3058, 1998. 0079 % 0080 % (C) 2018 by JMA, Johannes M. Arend 0081 % CP, Christoph Pörschmann 0082 % TH Köln - University of Applied Sciences 0083 % Institute of Communications Engineering 0084 % Department of Acoustics and Audio Signal Processing 0085 0086 function eqDataset = supdeq_getEqDataset(N, earDistance, NFFT, fs, waveType, sourceDistance, earPosition, sphereOffset, transformCore, swSettings) 0087 0088 if nargin < 1 || isempty(N) 0089 N = 35; 0090 end 0091 0092 if nargin < 2 || isempty(earDistance) 0093 earDistance = 0.165; 0094 end 0095 0096 if nargin < 3 || isempty(NFFT) 0097 NFFT = 512; 0098 end 0099 0100 if nargin < 4 || isempty(fs) 0101 fs = 48000; 0102 end 0103 0104 if nargin < 5 || isempty(waveType) 0105 waveType = 0; 0106 sourceDistance = 1; 0107 end 0108 0109 if nargin < 6 || isempty(sourceDistance) 0110 sourceDistance = 1; 0111 end 0112 0113 if nargin < 7 || isempty(earPosition) 0114 earPosition = [90, 90, 270, 90]; 0115 end 0116 0117 if nargin < 8 || isempty(sphereOffset) 0118 sphereOffset = []; 0119 end 0120 0121 if nargin < 9 || isempty(transformCore) 0122 transformCore = 'sofia'; 0123 end 0124 0125 if nargin < 10 || isempty(swSettings) 0126 swSettings = 1; 0127 end 0128 0129 %% Get SH-Coefficients with sofia_wgc 0130 0131 %Define required parameters 0132 radius = earDistance/2; 0133 if sourceDistance < radius 0134 error('Invalid source distance. Source must be outside the radius (earDistance/2)'); 0135 end 0136 ac = 2; %Array configuration, 2 - Rigid sphere with pressure transducers 0137 c = 343; %Speed of sound in m/s - Set to 343 here 0138 %Time delay in s - Set to 0 for plane waves or, to compensate for the time shift induced by spherical wave generation, 0139 %to negative time delay according to swSettings 0140 if waveType == 0 0141 delay = 0; 0142 else %waveType == 1 0143 0144 %This way, sofia and ak IRs are at the same position when reference in ak is set to r0 = 1m. Applied to all sofia data as a global delay 0145 %0.5141 was determined empirically with several tests... 0146 globalDelay = -1*((1-0.5141)/c); 0147 globalDelay2 = -1*((0.5141)/c); 0148 0149 if swSettings == 0 %Compensate time shift but maintain level changes 0150 delay = (-1*(sourceDistance/c)) - globalDelay2; 0151 elseif swSettings == 1 %Compensate time shift and compensate level changes (Default) 0152 delay = (-1*(sourceDistance/c)) - globalDelay2; 0153 levelComp = sourceDistance/1; %Reference at r0 = 1m 0154 elseif swSettings == 2 %Maintain time shift but compensate level changes 0155 delay = globalDelay; %Just the global delay to align sofia and ak 0156 levelComp = sourceDistance/1; %Reference at r0 = 1m 0157 elseif swSettings == 3 %Maintain time shift and maintain level changes (no compensation) 0158 delay = globalDelay; %Just the global delay to align sofia and ak 0159 end 0160 0161 end 0162 0163 %Convert earPosition to radiant 0164 appliedEarPosition = earPosition*pi/180; 0165 0166 %Shift earPosition if sphereOffset applied 0167 if ~isempty(sphereOffset) 0168 azL = appliedEarPosition(1); 0169 elL = appliedEarPosition(2); 0170 azR = appliedEarPosition(3); 0171 elR = appliedEarPosition(4); 0172 0173 [xL,yL,zL] = sph2cart(azL,pi/2-elL,radius); 0174 [xR,yR,zR] = sph2cart(azR,pi/2-elR,radius); 0175 0176 xL = xL - sphereOffset(1); %Use minus because sphere offset is inverse to ear offset 0177 xR = xR - sphereOffset(1); 0178 yL = yL - sphereOffset(2); 0179 yR = yR - sphereOffset(2); 0180 zL = zL - sphereOffset(3); 0181 zR = zR - sphereOffset(3); 0182 0183 [azL, elL, ~] = cart2sph(xL,yL,zL); 0184 azL = mod(azL,2*pi); 0185 elL = pi/2 - elL; 0186 [azR, elR, ~] = cart2sph(xR,yR,zR); 0187 azR = mod(azR,2*pi); 0188 elR = pi/2 - elR; 0189 0190 appliedEarPosition = [azL,elL,azR,elR]; 0191 end 0192 0193 %Calculate SH-Coefficients 0194 if strcmp(transformCore,'sofia') 0195 eqDataset.Hl_nm = sofia_wgc(N, radius, ac, fs, NFFT, appliedEarPosition(1), appliedEarPosition(2), delay,c,waveType,sourceDistance); 0196 eqDataset.Hr_nm = sofia_wgc(N, radius, ac, fs, NFFT, appliedEarPosition(3), appliedEarPosition(4), delay,c,waveType,sourceDistance); 0197 0198 %Apply level compensation if chosen with respective swSettings 0199 if waveType == 1 && (swSettings == 1 || swSettings == 2) 0200 eqDataset.Hl_nm = eqDataset.Hl_nm*levelComp; 0201 eqDataset.Hr_nm = eqDataset.Hr_nm*levelComp; 0202 end 0203 0204 0205 elseif strcmp(transformCore,'ak') 0206 0207 %Transform coordinates first 0208 appliedEarPositionAK = appliedEarPosition*180/pi; 0209 appliedEarPositionAK(2) = 90 - appliedEarPositionAK(2); 0210 appliedEarPositionAK(4) = 90 - appliedEarPositionAK(4); 0211 0212 %Get Lebedev sampling grid according to N 0213 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]; 0214 lebNumPoints = [6, 14, 26, 38, 50, 74, 86, 110, 146, 170, 194, 230, 266, 302, 350, 434, 590, 770, 974, 1202, 1454, 1730, 2030, 2354, 2702, 3074, 3470, 3890, 4334, 4802, 5294, 5810]; 0215 0216 if N > lebN(end) 0217 error('N = 65 is the highest spatial order possible with transformcore ak. Use sofia instead'); 0218 end 0219 0220 %Get closest N and add + 1 to be safe 0221 [~,cId] = min(abs(N-lebN)); 0222 if lebN(cId) < N 0223 cId = cId + 1; 0224 end 0225 0226 %Get sampling grid 0227 samplingGridAK = supdeq_lebedev(lebNumPoints(cId)); 0228 samplingGridAKSH = samplingGridAK; 0229 0230 %Transform coordinates again 0231 samplingGridAK(:,2) = 90-samplingGridAK(:,2); 0232 samplingGridAK = samplingGridAK(:,1:2); 0233 0234 %Define SH order for AKsphericalHead 0235 shOrderAK = 100; %Maybe set to lower order? 0236 0237 %Define distance - if wavetype = 0 (plane wave), set distance to 100 m 0238 if waveType == 0 0239 distanceAK = ones(size(samplingGridAK,1),1) * 100; 0240 samplingGridAK(:,3) = distanceAK; 0241 %%Define referenceDistance 0242 refDistance = distanceAK(1); 0243 end 0244 0245 if waveType == 1 % wavetype = 1 (spherical wave) 0246 distanceAK = ones(size(samplingGridAK,1),1) * sourceDistance; 0247 samplingGridAK(:,3) = distanceAK; 0248 0249 %Define parameters according to swSettings 0250 if swSettings == 0 %Compensate time shift but maintain level changes 0251 refDistance = distanceAK(1); %Setting reference distance and source distance to the same value leads to time- and level-compensation 0252 levelComp = 1/sourceDistance; %Reference at r0 = 1m 0253 elseif swSettings == 1 %Compensate time shift and compensate level changes (Default) 0254 refDistance = distanceAK(1); %Setting reference distance and source distance to the same value leads to time- and level-compensation 0255 elseif swSettings == 2 %Maintain time shift but compensate level changes 0256 refDistance = 1; 0257 levelComp = sourceDistance/1; %Reference at r0 = 1m 0258 elseif swSettings == 3 %Maintain time shift and maintain level changes (no compensation) 0259 refDistance = 1; %Set reference to 1 m and apply to compensation/shift 0260 end 0261 0262 end 0263 0264 %Get IRs 0265 irsAK = AKsphericalHead(samplingGridAK,appliedEarPositionAK,false,radius,refDistance,shOrderAK,NFFT,fs,c); 0266 irsAK_L = squeeze(irsAK(:,:,1)); 0267 irsAK_R = squeeze(irsAK(:,:,2)); 0268 0269 %Apply level compensation if chosen with respective swSettings 0270 if waveType == 1 && (swSettings == 0 || swSettings == 2) 0271 irsAK_L = irsAK_L*levelComp; 0272 irsAK_R = irsAK_R*levelComp; 0273 end 0274 0275 %Transform IRs to SH domain with given order N 0276 eqDataset.Hl_nm = AKsht(irsAK_L,true,samplingGridAKSH,N,'complex',fs); 0277 eqDataset.Hr_nm = AKsht(irsAK_R,true,samplingGridAKSH,N,'complex',fs); 0278 0279 end 0280 eqDataset.f = linspace(0,fs/2,NFFT/2+1); 0281 eqDataset.N = N; 0282 eqDataset.earDistance = earDistance; 0283 eqDataset.radius = radius; 0284 eqDataset.waveType = waveType; 0285 eqDataset.sourceDistance = sourceDistance; 0286 eqDataset.c = c; 0287 if isempty(sphereOffset) 0288 eqDataset.earPosition = earPosition; 0289 end 0290 if ~isempty(sphereOffset) 0291 eqDataset.inputEarPosition = earPosition; 0292 eqDataset.sphereOffset = sphereOffset; 0293 eqDataset.appliedEarPosition = appliedEarPosition*180/pi; 0294 end 0295 %if waveType == 1 0296 % eqDataset.tsComp = tsComp; 0297 %end 0298 eqDataset.transformCore = transformCore; 0299 0300 end 0301