% SUpDEq - Spatial Upsampling by Directional Equalization
function HRIRs_sfd = supdeq_sofa2sfd(SOFAobj, N, samplingGrid, FFToversize, transformCore, tikhEps)
This function transformes spherical HRIR measurements in SOFA format
('SimpleFreeFieldHRIR' convention) to the SH-domain by spherical Fourier
transform and writes a output struct suitable for the SUpDEq toolbox.
Output:
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 sourceDistance in m.
Input:
SOFAobj - Measurement data in form of a SOFA object
N - Transform order N
samplingGrid - Optional input [default = nan]
If a spatial sampling grid is passed, it should be a Qx3
matrix where the first column holds the azimuth, the
second the elevation, and the third the sampling weights.
In this case, the SH transform is performed with weights.
If no samplingGrid (or []) is passed, the sampling grid
according to the SOFA object is used for a spherical
Fourier transform without sampling weights. In this case,
the transformCore 'ak' is used by default.
Azimuth positions in degree. Can be scalar or vector of size(el)
(0=front, 90=left, 180=back, 270=right)
(0 points to positive x-axis, 90 to positive y-axis)
Elevations in degree. Can be scalar or vector of size(az)
(0=North Pole, 90=front, 180=South Pole)
(0 points to positive z-axis, 180 to negative z-axis)
FFToversize - FFToversize rises the FFT Blocksize [default = 1]
A FFT of the blocksize (FFToversize*NFFT) is applied
to the time domain data, where NFFT is determinded
as the next power of two of the signalSize which is
signalSize = (lastSample-firstSample).
transformCore - String to define method to be used for the spherical
Fourier transform.
'sofia - sofia_itc from SOFiA toolbox
'ak' - AKisht from AKtools
The results are exactly the same, but AKisht is faster
with big sampling grids
Default: 'sofia'
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.
Variable 'transformCore' is neglected when 'tikhEps' is
defined as the regularized least-square spherical Fourier
transform is applied directly without any third party toolbox.
Default: 0 (no Tikhonov regularization)
Dependencies: SOFiA toolbox, AKtools, SOFA API
References:
Benjamin Bernschütz: Microphone Arrays and Sound Field Decomposition
for Dynamic Binaural Recording. Ph.D. dissertation, Technical University
Berlin (2016).
Boaz Rafaely: Fundamentals of spherical array processing. In. Springer
topics in signal processing. Benesty, J.; Kellermann, W. (Eds.),
Springer, Heidelberg et al. (2015).
Franz Zotter: Analysis and synthesis of sound-radiation with spherical
arrays. Ph.D. dissertation, University of Music and Performing arts (2009).
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.
(C) 2018/2019 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 = supdeq_sofa2sfd(SOFAobj, N, samplingGrid, FFToversize, transformCore, tikhEps) 0004 % 0005 % This function transformes spherical HRIR measurements in SOFA format 0006 % ('SimpleFreeFieldHRIR' convention) to the SH-domain by spherical Fourier 0007 % transform and writes a output struct suitable for the SUpDEq toolbox. 0008 % 0009 % Output: 0010 % HRIRs_sfd - Struct with Spherical Harmonic coefficients 0011 % (SH-coefficients) for the left (Hl_nm) and right (Hr_nm) 0012 % channel/ear, absolute frequency scale f, 0013 % transform order N, FFToversize, and sourceDistance in m. 0014 % 0015 % Input: 0016 % SOFAobj - Measurement data in form of a SOFA object 0017 % N - Transform order N 0018 % samplingGrid - Optional input [default = nan] 0019 % If a spatial sampling grid is passed, it should be a Qx3 0020 % matrix where the first column holds the azimuth, the 0021 % second the elevation, and the third the sampling weights. 0022 % In this case, the SH transform is performed with weights. 0023 % If no samplingGrid (or []) is passed, the sampling grid 0024 % according to the SOFA object is used for a spherical 0025 % Fourier transform without sampling weights. In this case, 0026 % the transformCore 'ak' is used by default. 0027 % Azimuth positions in degree. Can be scalar or vector of size(el) 0028 % (0=front, 90=left, 180=back, 270=right) 0029 % (0 points to positive x-axis, 90 to positive y-axis) 0030 % Elevations in degree. Can be scalar or vector of size(az) 0031 % (0=North Pole, 90=front, 180=South Pole) 0032 % (0 points to positive z-axis, 180 to negative z-axis) 0033 % FFToversize - FFToversize rises the FFT Blocksize [default = 1] 0034 % A FFT of the blocksize (FFToversize*NFFT) is applied 0035 % to the time domain data, where NFFT is determinded 0036 % as the next power of two of the signalSize which is 0037 % signalSize = (lastSample-firstSample). 0038 % transformCore - String to define method to be used for the spherical 0039 % Fourier transform. 0040 % 'sofia - sofia_itc from SOFiA toolbox 0041 % 'ak' - AKisht from AKtools 0042 % The results are exactly the same, but AKisht is faster 0043 % with big sampling grids 0044 % Default: 'sofia' 0045 % tikhEps - Define epsilon of Tikhonov regularization if 0046 % regularization should be applied 0047 % Applying the Tikhonov regularization will always result 0048 % in a least-square fit solution for the SH transform. 0049 % Variable 'transformCore' is neglected when 'tikhEps' is 0050 % defined as the regularized least-square spherical Fourier 0051 % transform is applied directly without any third party toolbox. 0052 % Default: 0 (no Tikhonov regularization) 0053 % 0054 % Dependencies: SOFiA toolbox, AKtools, SOFA API 0055 % 0056 % References: 0057 % Benjamin Bernschütz: Microphone Arrays and Sound Field Decomposition 0058 % for Dynamic Binaural Recording. Ph.D. dissertation, Technical University 0059 % Berlin (2016). 0060 % 0061 % Boaz Rafaely: Fundamentals of spherical array processing. In. Springer 0062 % topics in signal processing. Benesty, J.; Kellermann, W. (Eds.), 0063 % Springer, Heidelberg et al. (2015). 0064 % 0065 % Franz Zotter: Analysis and synthesis of sound-radiation with spherical 0066 % arrays. Ph.D. dissertation, University of Music and Performing arts (2009). 0067 % 0068 % R. Duraiswami, D. N. Zotkin, and N. A. Gumerov, ?Interpolation and range 0069 % extrapolation of HRTFs,? in Proceedings of the IEEE International 0070 % Conference on Acoustics, Speech, and Signal Processing, 2004, pp. IV45?IV48. 0071 % 0072 % (C) 2018/2019 by JMA, Johannes M. Arend 0073 % TH Köln - University of Applied Sciences 0074 % Institute of Communications Engineering 0075 % Department of Acoustics and Audio Signal Processing 0076 0077 function HRIRs_sfd = supdeq_sofa2sfd(SOFAobj, N, samplingGrid, FFToversize, transformCore, tikhEps) 0078 0079 if nargin < 3 || isempty(samplingGrid) 0080 samplingGrid = []; 0081 end 0082 0083 if nargin < 4 || isempty(FFToversize) 0084 FFToversize = 1; 0085 end 0086 0087 if nargin < 5 || isempty(transformCore) 0088 transformCore = 'sofia'; 0089 end 0090 0091 if nargin < 6 || isempty(tikhEps) 0092 tikhEps = 0; 0093 end 0094 0095 %Check FFToversize parameter 0096 if FFToversize < 1 0097 warning('FFToversize needs to >= 1. Set to default = 1'); 0098 FFToversize = 1; 0099 end 0100 0101 %Check grid if passed 0102 weightsPassed = true; 0103 if ~isempty(samplingGrid) 0104 if size(samplingGrid,2) > size(samplingGrid,1) 0105 samplingGrid = samplingGrid'; 0106 warning('Assuming grid was passed in wrong dimensions - Grid flipped'); 0107 end 0108 0109 if size(samplingGrid,2) ~= 3 0110 weightsPassed = false; 0111 end 0112 end 0113 0114 %Set transformCore to 'ak' if no weights passed but sofia is choosen 0115 if ~weightsPassed && strcmp(transformCore,'sofia') 0116 disp('Changed transformCore to ak as no sampling weights were passed!'); 0117 transformCore = 'ak'; 0118 end 0119 0120 %Check SOFA object 0121 if ~strcmp(SOFAobj.GLOBAL_SOFAConventions,'SimpleFreeFieldHRIR') 0122 error('Function only valid for SOFA convention "SimpleFreeFieldHRIR".'); 0123 end 0124 0125 %% If sampling grid with weights passed and transformCore sofia 0126 0127 if tikhEps == 0 %Without Tikhonov Regularization 0128 0129 if ~isempty(samplingGrid) && weightsPassed && strcmp(transformCore,'sofia') 0130 0131 %Convert az and el to rad again 0132 samplingGrid(:,1:2) = samplingGrid(:,1:2) * pi / 180; 0133 0134 %Write timeData struct needed for sofia_fdt 0135 timeDataL.impulseResponses = squeeze(SOFAobj.Data.IR(:,1,:)); 0136 timeDataL.FS = SOFAobj.Data.SamplingRate; 0137 timeDataL.radius = abs(SOFAobj.ReceiverPosition(3)); 0138 timeDataL.averageAirTemp = 29.3100; %Number according to miro object HRIR_L2702 0139 0140 timeDataR.impulseResponses = squeeze(SOFAobj.Data.IR(:,2,:)); 0141 timeDataR.FS = SOFAobj.Data.SamplingRate; 0142 timeDataR.radius = abs(SOFAobj.ReceiverPosition(3)); 0143 timeDataR.averageAirTemp = 29.3100; %Number according to miro object HRIR_L2702 0144 0145 %Transform to Fourier domain with regular FFT (sofia_fdt) 0146 fftDataL = sofia_fdt(timeDataL,FFToversize); 0147 [fftDataR, ~, f] = sofia_fdt(timeDataR,FFToversize); 0148 0149 %Transform to SH-domain with spherical Fourier transform (sofia_stc) 0150 %Get SH-coefficients for left and right channel 0151 Hl_nm = sofia_stc(N,fftDataL,samplingGrid); 0152 Hr_nm = sofia_stc(N,fftDataR,samplingGrid); 0153 0154 %Write output struct 0155 HRIRs_sfd.Hl_nm = Hl_nm; 0156 HRIRs_sfd.Hr_nm = Hr_nm; 0157 HRIRs_sfd.f = f; 0158 HRIRs_sfd.N = N; 0159 HRIRs_sfd.FFToversize = FFToversize; 0160 HRIRs_sfd.sourceDistance = SOFAobj.SourcePosition(1,3); 0161 0162 disp('Transformation done with SOFiA toolbox - Sampling weights passed'); 0163 end 0164 end 0165 0166 %% If sampling grid passed with or without weights and transformCore ak 0167 0168 if tikhEps == 0 %Without Tikhonov Regularization 0169 0170 if ~isempty(samplingGrid) && strcmp(transformCore,'ak') 0171 0172 %Get IRs needed for AKsht 0173 %Arrays need to be flipped because AKsht needs vectors with [M x N] 0174 %(columns = channels) 0175 irL = squeeze(SOFAobj.Data.IR(:,1,:))'; 0176 irR = squeeze(SOFAobj.Data.IR(:,2,:))'; 0177 0178 %Zero pad to get length according to FFToversize 0179 if FFToversize > 1 0180 NFFT = 2^nextpow2(size(irL,1)); 0181 NFFT = NFFT*FFToversize; 0182 irL = [irL; zeros(NFFT-size(irL,1),size(irL,2))]; 0183 irR = [irR; zeros(NFFT-size(irR,1),size(irR,2))]; 0184 else 0185 if size(irL,1) == 2^nextpow2(size(irL,1)) 0186 NFFT = size(irL,1); 0187 else 0188 NFFT = 2^nextpow2(size(irL,1)); 0189 irL = [irL; zeros(NFFT-size(irL,1),size(irL,2))]; 0190 irR = [irR; zeros(NFFT-size(irR,1),size(irR,2))]; 0191 end 0192 end 0193 0194 %Transform to SH-domain with spherical Fourier transform (AKsht) 0195 %Get SH-coefficients for left and right channel 0196 Hl_nm = AKsht(irL,true,samplingGrid,N,'complex',SOFAobj.Data.SamplingRate); 0197 [Hr_nm, f] = AKsht(irR,true,samplingGrid,N,'complex',SOFAobj.Data.SamplingRate); 0198 0199 %Write output struct 0200 HRIRs_sfd.Hl_nm = Hl_nm; 0201 HRIRs_sfd.Hr_nm = Hr_nm; 0202 HRIRs_sfd.f = f'; 0203 HRIRs_sfd.N = N; 0204 HRIRs_sfd.FFToversize = FFToversize; 0205 HRIRs_sfd.sourceDistance = SOFAobj.SourcePosition(1,3); 0206 0207 if weightsPassed 0208 disp('Transformation done with AKTools - Sampling weights passed'); 0209 else 0210 disp('Transformation done with AKTools - No sampling weights passed'); 0211 end 0212 end 0213 end 0214 0215 %% If no sampling grid passed - Use AKsht 0216 0217 if tikhEps == 0 %Without Tikhonov Regularization 0218 0219 if isempty(samplingGrid) 0220 0221 %Get samplingGrid from SOFA object 0222 samplingGrid = SOFAobj.SourcePosition(:,1:2); 0223 0224 %Transform samplingGrid back to AKsht coordinate system 0225 samplingGrid(:,1) = mod(samplingGrid(:,1),360); 0226 samplingGrid(:,2) = 90-samplingGrid(:,2); 0227 0228 %Get IRs needed for AKsht 0229 %Arrays need to be flipped because AKsht needs vectors with [M x N] 0230 %(columns = channels) 0231 irL = squeeze(SOFAobj.Data.IR(:,1,:))'; 0232 irR = squeeze(SOFAobj.Data.IR(:,2,:))'; 0233 0234 %Zero pad to get length according to FFToversize 0235 if FFToversize > 1 0236 NFFT = 2^nextpow2(size(irL,1)); 0237 NFFT = NFFT*FFToversize; 0238 irL = [irL; zeros(NFFT-size(irL,1),size(irL,2))]; 0239 irR = [irR; zeros(NFFT-size(irR,1),size(irR,2))]; 0240 else 0241 if size(irL,1) == 2^nextpow2(size(irL,1)) 0242 NFFT = size(irL,1); 0243 else 0244 NFFT = 2^nextpow2(size(irL,1)); 0245 irL = [irL; zeros(NFFT-size(irL,1),size(irL,2))]; 0246 irR = [irR; zeros(NFFT-size(irR,1),size(irR,2))]; 0247 end 0248 end 0249 0250 %Transform to SH-domain with spherical Fourier transform (AKsht) 0251 %Get SH-coefficients for left and right channel 0252 Hl_nm = AKsht(irL,true,samplingGrid,N,'complex',SOFAobj.Data.SamplingRate); 0253 [Hr_nm, f] = AKsht(irR,true,samplingGrid,N,'complex',SOFAobj.Data.SamplingRate); 0254 0255 %Write output struct 0256 HRIRs_sfd.Hl_nm = Hl_nm; 0257 HRIRs_sfd.Hr_nm = Hr_nm; 0258 HRIRs_sfd.f = f'; 0259 HRIRs_sfd.N = N; 0260 HRIRs_sfd.FFToversize = FFToversize; 0261 HRIRs_sfd.sourceDistance = SOFAobj.SourcePosition(1,3); 0262 0263 disp('Transformation done with AKTools - Applied sampling grid from SOFA object without weights'); 0264 end 0265 end 0266 0267 %% If tikhEps is not zero - Calculate SH transform with Tikhonov regularization 0268 0269 if tikhEps ~= 0 0270 0271 sofaGrid = false; 0272 if isempty(samplingGrid) %If no sampling grid passed 0273 %Get samplingGrid from SOFA object 0274 samplingGrid = SOFAobj.SourcePosition(:,1:2); 0275 %Transform samplingGrid back to SUpDEq coordinate system 0276 samplingGrid(:,1) = mod(samplingGrid(:,1),360); 0277 samplingGrid(:,2) = 90-samplingGrid(:,2); 0278 %Set weightsPassed 0279 weightsPassed = false; 0280 sofaGrid = true; 0281 end 0282 0283 %Get IRs 0284 irL = squeeze(SOFAobj.Data.IR(:,1,:)); 0285 irR = squeeze(SOFAobj.Data.IR(:,2,:)); 0286 0287 %Zero pad to get length according to FFToversize 0288 if FFToversize > 1 0289 NFFT = 2^nextpow2(size(irL,2)); 0290 NFFT = NFFT*FFToversize; 0291 irL = [irL, zeros(size(irL,1),NFFT-size(irL,2))]; 0292 irR = [irR, zeros(size(irR,1),NFFT-size(irR,2))]; 0293 else 0294 if size(irL,2) == 2^nextpow2(size(irL,2)) 0295 NFFT = size(irL,2); 0296 else 0297 NFFT = 2^nextpow2(size(irL,2)); 0298 irL = [irL, zeros(size(irL,1),NFFT-size(irL,2))]; 0299 irR = [irR, zeros(size(irR,1),NFFT-size(irR,2))]; 0300 end 0301 end 0302 0303 %Transform HRIRs to frequency domain (HRTFs) 0304 HRTF_L = fft(irL,[],2); 0305 HRTF_R = fft(irR,[],2); 0306 %Cut at fs/2 0307 HRTF_L = HRTF_L(:,1:end/2+1); 0308 HRTF_R = HRTF_R(:,1:end/2+1); 0309 0310 %Write f 0311 f = linspace(0,SOFAobj.Data.SamplingRate/2,NFFT/2+1); 0312 0313 %Get SH functions, weights omitted 0314 [Ynm,n] = AKsh(N,[],samplingGrid(:,1),samplingGrid(:,2)); 0315 n = n'; 0316 nSH = (N+1).^2; 0317 0318 %Create diagonal matrix according to Duraiswami2004 0319 I = eye(nSH); 0320 D = diag(1 + n.*(n+1)) .* I; 0321 0322 % Inverse SH matrix for Least-Square SH transform with Tikhonov regularization 0323 YnmInvTik = (Ynm' * Ynm + tikhEps*D)^-1 * Ynm'; 0324 0325 %Get SH-coefficients for left and right channel 0326 Hl_nm = YnmInvTik*HRTF_L; 0327 Hr_nm = YnmInvTik*HRTF_R; 0328 0329 %Write output struct 0330 HRIRs_sfd.Hl_nm = Hl_nm; 0331 HRIRs_sfd.Hr_nm = Hr_nm; 0332 HRIRs_sfd.f = f; 0333 HRIRs_sfd.N = N; 0334 HRIRs_sfd.FFToversize = FFToversize; 0335 HRIRs_sfd.sourceDistance = SOFAobj.SourcePosition(1,3); 0336 HRIRs_sfd.tikhEps = tikhEps; 0337 0338 if sofaGrid 0339 disp('Transformation done with least-squares method with Tikhonov regularization - Applied sampling grid from SOFA object'); 0340 else 0341 disp('Transformation done with least-squares method with Tikhonov regularization'); 0342 end 0343 0344 end 0345 end 0346