% SUpDEq - Spatial Upsampling by Directional Equalization
function [HRIR_L, HRIR_R] = supdeq_getArbHRIR(HRIRs_sfd,samplingGrid,mode,channel,transformCore)
This function returns a HRIR / HRTF for an arbitrary direction, based on
interpolation in the SH-domain and inverse spherical Fourier
transform.
Output:
HRIR_L/R - Result of ISFT and IFFT in time domain. HRIR for L/R
channel at each spatial sampling point
Input:
HRIRs_sfd - Struct with HRIRs in spatial Fourier domain (sfd) / SH-domain
samplingGrid - Spatial sampling grid (Q x 2 matrix), 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)
mode - 'DEG' or 'RAD' (radiant or degree).
Default: 'DEG'
channel - Integer with 0, 1, or 2
0 - Only left channel
1 - Only right channel
Default: 2 - Both channels (stereo)
transformCore - String to define method to be used for the inverse
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'
Dependencies: SOFiA toolbox, AKtools
Reference:
Bernschütz, B. (2016). Microphone Arrays and Sound Field Decomposition
for Dynamic Binaural Recording. TU Berlin.
(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 [HRIR_L, HRIR_R] = supdeq_getArbHRIR(HRIRs_sfd,samplingGrid,mode,channel,transformCore) 0004 % 0005 % This function returns a HRIR / HRTF for an arbitrary direction, based on 0006 % interpolation in the SH-domain and inverse spherical Fourier 0007 % transform. 0008 % 0009 % Output: 0010 % HRIR_L/R - Result of ISFT and IFFT in time domain. HRIR for L/R 0011 % channel at each spatial sampling point 0012 % 0013 % Input: 0014 % HRIRs_sfd - Struct with HRIRs in spatial Fourier domain (sfd) / SH-domain 0015 % samplingGrid - Spatial sampling grid (Q x 2 matrix), where the first 0016 % column holds the azimuth and the second the 0017 % elevation (both in degree). 0018 % Azimuth in degree (0=front, 90=left, 180=back, 270=right) 0019 % (0 points to positive x-axis, 90 to positive y-axis) 0020 % Elevations in degree (0=North Pole, 90=front, 180=South Pole) 0021 % (0 points to positive z-axis, 180 to negative z-axis) 0022 % mode - 'DEG' or 'RAD' (radiant or degree). 0023 % Default: 'DEG' 0024 % channel - Integer with 0, 1, or 2 0025 % 0 - Only left channel 0026 % 1 - Only right channel 0027 % Default: 2 - Both channels (stereo) 0028 % transformCore - String to define method to be used for the inverse 0029 % spherical Fourier transform. 0030 % 'sofia - sofia_itc from SOFiA toolbox 0031 % 'ak' - AKisht from AKtools 0032 % The results are exactly the same, but AKisht is faster 0033 % with big sampling grids 0034 % Default: 'sofia' 0035 % 0036 % Dependencies: SOFiA toolbox, AKtools 0037 % 0038 % Reference: 0039 % Bernschütz, B. (2016). Microphone Arrays and Sound Field Decomposition 0040 % for Dynamic Binaural Recording. TU Berlin. 0041 % 0042 % (C) 2018 by JMA, Johannes M. Arend 0043 % TH Köln - University of Applied Sciences 0044 % Institute of Communications Engineering 0045 % Department of Acoustics and Audio Signal Processing 0046 0047 function [HRIR_L, HRIR_R] = supdeq_getArbHRIR(HRIRs_sfd,samplingGrid,mode,channel,transformCore) 0048 0049 if isfield(HRIRs_sfd,'FFToversize') 0050 FFToversize = HRIRs_sfd.FFToversize; 0051 else 0052 FFToversize = 1; 0053 end 0054 0055 if nargin < 3 || isempty(mode) 0056 mode = 'DEG'; 0057 end 0058 0059 if nargin < 4 || isempty(channel) 0060 channel = 2; 0061 end 0062 0063 if nargin < 5 || isempty(transformCore) 0064 transformCore = 'sofia'; 0065 end 0066 0067 %Get angles 0068 az = samplingGrid(:,1); 0069 el = samplingGrid(:,2); 0070 0071 %SOFiA needs RAD 0072 if strcmp(transformCore,'sofia') 0073 if strcmp(mode,'DEG') 0074 az = az*pi/180; 0075 el = el*pi/180; 0076 end 0077 end 0078 0079 %AK needs DEG 0080 if strcmp(transformCore,'ak') 0081 if strcmp(mode,'RAD') 0082 az = az*180/pi; 0083 el = el*180/pi; 0084 end 0085 end 0086 0087 %Check range 0088 if sum(az < 0) > 0 || sum(el < 0) > 0 0089 error('Only positive azimuth/elevation values allowed!'); 0090 end 0091 0092 %% Perform transform with sofia_itc 0093 0094 if strcmp(transformCore,'sofia') 0095 0096 if channel == 0 %Only left channel 0097 0098 %Perform inverse spherical Fourier transform 0099 Hl = sofia_itc(HRIRs_sfd.Hl_nm, [az el]); 0100 0101 %Get mirror spectrum 0102 HlTwoSided = conj([Hl(:,:), conj(fliplr(Hl(:,2:end-1)))])'; 0103 0104 %Transform back to time domain 0105 HRIR_L = real(ifft(HlTwoSided)); 0106 0107 %Cut zeros depending on FFToversize 0108 finalLength = size(HRIR_L,1)/FFToversize; 0109 HRIR_L = HRIR_L(1:finalLength,:); 0110 HRIR_R = nan; 0111 0112 elseif channel == 1 %Only right channel 0113 0114 %Perform inverse spherical Fourier transform 0115 Hr = sofia_itc(HRIRs_sfd.Hr_nm, [az el]); 0116 0117 %Get mirror spectrum 0118 HrTwoSided = conj([Hr(:,:), conj(fliplr(Hr(:,2:end-1)))])'; 0119 0120 %Transform back to time domain 0121 HRIR_R = real(ifft(HrTwoSided)); 0122 0123 %Cut zeros depending on FFToversize 0124 finalLength = size(HRIR_R,1)/FFToversize; 0125 HRIR_L = nan; 0126 HRIR_R = HRIR_R(1:finalLength,:); 0127 0128 elseif channel == 2 %Default - Stereo 0129 0130 %Perform inverse spherical Fourier transform 0131 Hl = sofia_itc(HRIRs_sfd.Hl_nm, [az el]); 0132 Hr = sofia_itc(HRIRs_sfd.Hr_nm, [az el]); 0133 0134 %Get mirror spectrum 0135 HlTwoSided = conj([Hl(:,:), conj(fliplr(Hl(:,2:end-1)))])'; 0136 HrTwoSided = conj([Hr(:,:), conj(fliplr(Hr(:,2:end-1)))])'; 0137 0138 %Transform back to time domain 0139 HRIR_L = real(ifft(HlTwoSided)); 0140 HRIR_R = real(ifft(HrTwoSided)); 0141 0142 %Cut zeros depending on FFToversize 0143 finalLength = size(HRIR_L,1)/FFToversize; 0144 HRIR_L = HRIR_L(1:finalLength,:); 0145 HRIR_R = HRIR_R(1:finalLength,:); 0146 0147 end 0148 end 0149 0150 %% Perform transform with AKisht 0151 0152 if strcmp(transformCore,'ak') 0153 0154 if channel == 0 %Only left channel 0155 0156 %Perform inverse spherical Fourier transform 0157 HRIR_L = AKisht(HRIRs_sfd.Hl_nm,true,[az el],'complex'); 0158 0159 %Cut zeros depending on FFToversize 0160 finalLength = size(HRIR_L,1)/FFToversize; 0161 HRIR_L = HRIR_L(1:finalLength,:); 0162 HRIR_R = nan; 0163 0164 elseif channel == 1 %Only right channel 0165 0166 %Perform inverse spherical Fourier transform 0167 HRIR_R = AKisht(HRIRs_sfd.Hr_nm,true,[az el],'complex'); 0168 0169 %Cut zeros depending on FFToversize 0170 finalLength = size(HRIR_R,1)/FFToversize; 0171 HRIR_L = nan; 0172 HRIR_R = HRIR_R(1:finalLength,:); 0173 0174 elseif channel == 2 %Default - Stereo 0175 0176 %Perform inverse spherical Fourier transform 0177 HRIR_L = AKisht(HRIRs_sfd.Hl_nm,true,[az el],'complex'); 0178 HRIR_R = AKisht(HRIRs_sfd.Hr_nm,true,[az el],'complex'); 0179 0180 %Cut zeros depending on FFToversize 0181 finalLength = size(HRIR_L,1)/FFToversize; 0182 HRIR_L = HRIR_L(1:finalLength,:); 0183 HRIR_R = HRIR_R(1:finalLength,:); 0184 0185 end 0186 end 0187 0188 end 0189 0190 0191 0192