% SUpDEq - Spatial Upsampling by Directional Equalization
function HRIRs_sfd = supdeq_hrtf2sfd(HRTF_L, HRTF_R, N, samplingGrid, fs, transformCore, tikhEps)
This function transformes HRTFs to the SH-domain by spherical Fourier transform.
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, and FFToversize
Input:
HRTF_L/R - HRTF for left and right ear as [n x m] samples, with n =
number of channels / sampling points and m = FFT bins
N - Transform order N
samplingGrid - Spatial sampling grid. Can 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 SOFiA
functions "sofia_fdt" and "sofia_stc" as well as "AKsht"
can be used
If no weights are passed (Qx2 matrix), the function "AKsht"
is used providing a least-squares solution of the spherical
Fourier transform.
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)
fs - Sampling rate. Only needed to write frequency vector or
if AKsht is used.
Default: 48000
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. Depending on the sampling grids, weights are
applied or not.
Default: 0 (no Tikhonov regularization)
Dependencies: SOFiA toolbox, AKtools
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_hrtf2sfd(HRTF_L, HRTF_R, N, samplingGrid, fs, transformCore, tikhEps) 0004 % 0005 % This function transformes HRTFs to the SH-domain by spherical Fourier transform. 0006 % 0007 % Output: 0008 % HRIRs_sfd - Struct with Spherical Harmonic coefficients 0009 % (SH-Coefficients) for the left (Hl_nm) and right (Hr_nm) 0010 % channel/ear, absolute frequency scale f, 0011 % transform order N, and FFToversize 0012 % 0013 % Input: 0014 % HRTF_L/R - HRTF for left and right ear as [n x m] samples, with n = 0015 % number of channels / sampling points and m = FFT bins 0016 % N - Transform order N 0017 % samplingGrid - Spatial sampling grid. Can be a Qx3 matrix where the first 0018 % column holds the azimuth, the second the elevation, and 0019 % the third the sampling weights. In this case, the SOFiA 0020 % functions "sofia_fdt" and "sofia_stc" as well as "AKsht" 0021 % can be used 0022 % If no weights are passed (Qx2 matrix), the function "AKsht" 0023 % is used providing a least-squares solution of the spherical 0024 % Fourier transform. 0025 % Azimuth positions in degree. Can be scalar or vector of size(el) 0026 % (0=front, 90=left, 180=back, 270=right) 0027 % (0 points to positive x-axis, 90 to positive y-axis) 0028 % Elevations in degree. Can be scalar or vector of size(az) 0029 % (0=North Pole, 90=front, 180=South Pole) 0030 % (0 points to positive z-axis, 180 to negative z-axis) 0031 % fs - Sampling rate. Only needed to write frequency vector or 0032 % if AKsht is used. 0033 % Default: 48000 0034 % transformCore - String to define method to be used for the spherical 0035 % Fourier transform. 0036 % 'sofia - sofia_itc from SOFiA toolbox 0037 % 'ak' - AKisht from AKtools 0038 % The results are exactly the same, but AKisht is faster 0039 % with big sampling grids 0040 % Default: 'sofia' 0041 % tikhEps - Define epsilon of Tikhonov regularization if 0042 % regularization should be applied 0043 % Applying the Tikhonov regularization will always result 0044 % in a least-square fit solution for the SH transform. 0045 % Variable 'transformCore' is neglected when 'tikhEps' is 0046 % defined as the regularized least-square spherical Fourier 0047 % transform is applied directly without any third party 0048 % toolbox. Depending on the sampling grids, weights are 0049 % applied or not. 0050 % Default: 0 (no Tikhonov regularization) 0051 % 0052 % Dependencies: SOFiA toolbox, AKtools 0053 % 0054 % References: 0055 % Benjamin Bernschütz: Microphone Arrays and Sound Field Decomposition 0056 % for Dynamic Binaural Recording. Ph.D. dissertation, Technical University 0057 % Berlin (2016). 0058 % 0059 % Boaz Rafaely: Fundamentals of spherical array processing. In. Springer 0060 % topics in signal processing. Benesty, J.; Kellermann, W. (Eds.), 0061 % Springer, Heidelberg et al. (2015). 0062 % 0063 % Franz Zotter: Analysis and synthesis of sound-radiation with spherical 0064 % arrays. Ph.D. dissertation, University of Music and Performing arts (2009). 0065 % 0066 % R. Duraiswami, D. N. Zotkin, and N. A. Gumerov, ?Interpolation and range 0067 % extrapolation of HRTFs,? in Proceedings of the IEEE International 0068 % Conference on Acoustics, Speech, and Signal Processing, 2004, pp. IV45?IV48. 0069 % 0070 % (C) 2018/2019 by JMA, Johannes M. Arend 0071 % TH Köln - University of Applied Sciences 0072 % Institute of Communications Engineering 0073 % Department of Acoustics and Audio Signal Processing 0074 0075 function HRIRs_sfd = supdeq_hrtf2sfd(HRTF_L, HRTF_R, N, samplingGrid, fs, transformCore, tikhEps) 0076 0077 if nargin < 5 || isempty(fs) 0078 fs = 48000; 0079 end 0080 0081 if nargin < 6 || isempty(transformCore) 0082 transformCore = 'sofia'; 0083 end 0084 0085 if nargin < 7 || isempty(tikhEps) 0086 tikhEps = 0; 0087 end 0088 0089 %Check passed grid for weights 0090 weightsPassed = true; 0091 if size(samplingGrid,2) ~= 3 0092 %warning('No sampling weights passed. Using AKsht to perform SHT with pseudo-inverse SH-matrix') 0093 weightsPassed = false; 0094 end 0095 0096 %% If sampling grid with weights passed - use sofia_stc or AKsht 0097 0098 if tikhEps == 0 %Without Tikhonov Regularization 0099 0100 if weightsPassed 0101 0102 if strcmp(transformCore,'sofia') 0103 0104 %Convert az and el to rad again 0105 samplingGrid(:,1:2) = samplingGrid(:,1:2) * pi / 180; 0106 0107 %Transform to SH-domain with spherical Fourier transform (sofia_stc) 0108 %Get SH-coefficients for left and right channel 0109 Hl_nm = sofia_stc(N,HRTF_L,samplingGrid); 0110 Hr_nm = sofia_stc(N,HRTF_R,samplingGrid); 0111 0112 %Write output struct 0113 HRIRs_sfd.Hl_nm = Hl_nm; 0114 HRIRs_sfd.Hr_nm = Hr_nm; 0115 HRIRs_sfd.f = linspace(0,fs/2,size(HRTF_L,2)); 0116 HRIRs_sfd.N = N; 0117 HRIRs_sfd.FFToversize = 1; 0118 0119 disp('Transformation done with SOFiA toolbox - Sampling weights passed'); 0120 0121 elseif strcmp(transformCore,'ak') 0122 0123 %Transform to SH-domain with spherical Fourier transform (AKsht) 0124 %Get SH-coefficients for left and right channel 0125 Hl_nm = AKsht(HRTF_L.',false,samplingGrid,N,'complex',fs); 0126 Hr_nm = AKsht(HRTF_R.',false,samplingGrid,N,'complex',fs); 0127 0128 %Write output struct 0129 HRIRs_sfd.Hl_nm = Hl_nm; 0130 HRIRs_sfd.Hr_nm = Hr_nm; 0131 HRIRs_sfd.f = linspace(0,fs/2,size(HRTF_L,2)); 0132 HRIRs_sfd.N = N; 0133 HRIRs_sfd.FFToversize = 1; 0134 0135 disp('Transformation done with AKTools - Sampling weights passed'); 0136 0137 end 0138 end 0139 end 0140 0141 %% If no weights passed - Use always AKsht 0142 if tikhEps == 0 %Without Tikhonov Regularization 0143 if ~weightsPassed 0144 0145 %Transform to SH-domain with spherical Fourier transform (AKsht) 0146 %Get SH-coefficients for left and right channel 0147 Hl_nm = AKsht(HRTF_L.',false,samplingGrid,N,'complex',fs); 0148 Hr_nm = AKsht(HRTF_R.',false,samplingGrid,N,'complex',fs); 0149 0150 %Write output struct 0151 HRIRs_sfd.Hl_nm = Hl_nm; 0152 HRIRs_sfd.Hr_nm = Hr_nm; 0153 HRIRs_sfd.f = linspace(0,fs/2,size(HRTF_L,2)); 0154 HRIRs_sfd.N = N; 0155 HRIRs_sfd.FFToversize = 1; 0156 0157 disp('Transformation done with AKTools - No sampling weights passed'); 0158 end 0159 end 0160 0161 %% If tikhEps is not zero - Calculate SH transform with Tikhonov regularization 0162 0163 if tikhEps ~= 0 0164 0165 %Get SH functions, weights omitted 0166 [Ynm,n] = AKsh(N,[],samplingGrid(:,1),samplingGrid(:,2)); 0167 n = n'; 0168 nSH = (N+1).^2; 0169 0170 %Create diagonal matrix according to Duraiswami2004 0171 I = eye(nSH); 0172 D = diag(1 + n.*(n+1)) .* I; 0173 0174 % Inverse SH matrix for Least-Square SH transform with Tikhonov regularization 0175 YnmInvTik = (Ynm' * Ynm + tikhEps*D)^-1 * Ynm'; 0176 0177 %Get SH-coefficients for left and right channel 0178 Hl_nm = YnmInvTik*HRTF_L; 0179 Hr_nm = YnmInvTik*HRTF_R; 0180 0181 %Write output struct 0182 HRIRs_sfd.Hl_nm = Hl_nm; 0183 HRIRs_sfd.Hr_nm = Hr_nm; 0184 HRIRs_sfd.f = linspace(0,fs/2,size(HRTF_L,2)); 0185 HRIRs_sfd.N = N; 0186 HRIRs_sfd.FFToversize = 1; 0187 HRIRs_sfd.tikhEps = tikhEps; 0188 0189 disp('Transformation done with least-squares method with Tikhonov regularization'); 0190 0191 end 0192 0193 end 0194