diff --git a/scripts/binauralRenderer_interface/matlab_hrir_generation_scripts/README.md b/scripts/binauralRenderer_interface/matlab_hrir_generation_scripts/README.md new file mode 100644 index 0000000000000000000000000000000000000000..f6eb2cfad6ea65309e27c337e672b3fe2945b238 --- /dev/null +++ b/scripts/binauralRenderer_interface/matlab_hrir_generation_scripts/README.md @@ -0,0 +1,7 @@ + +### 2023-07-17 +Entry point to convert SD HRIRs to SHD HRIRs is the convert_SD2SHD_HRIRs.m script. +Entry point to convert SHD HRIRs to CLDFB domain HRIRs is SHD_2_ROM.m. + +Python 3.9.x needs to be installed with the sofar python module. The convert_SD2SHD_HRIRs.m script needs the path to this python. +Matlab 2020 has seen some issues with using python, so use a newer version if possible. If you are unable to use a newer matlab version, try loading in the OS default python, as this works sometimes. diff --git a/scripts/binauralRenderer_interface/matlab_hrir_generation_scripts/SHD_2_ROM.m b/scripts/binauralRenderer_interface/matlab_hrir_generation_scripts/SHD_2_ROM.m new file mode 100644 index 0000000000000000000000000000000000000000..3fb9a9e8ea4e9932cd9ca5307787c1d131af5322 --- /dev/null +++ b/scripts/binauralRenderer_interface/matlab_hrir_generation_scripts/SHD_2_ROM.m @@ -0,0 +1,132 @@ +%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% +% +% (C) 2022-2023 IVAS codec Public Collaboration with portions copyright Dolby International AB, Ericsson AB, +% Fraunhofer-Gesellschaft zur Foerderung der angewandten Forschung e.V., Huawei Technologies Co. LTD., +% Koninklijke Philips N.V., Nippon Telegraph and Telephone Corporation, Nokia Technologies Oy, Orange, +% Panasonic Holdings Corporation, Qualcomm Technologies, Inc., VoiceAge Corporation, and other +% contributors to this repository. All Rights Reserved. +% +% This software is protected by copyright law and by international treaties. +% The IVAS codec Public Collaboration consisting of Dolby International AB, Ericsson AB, +% Fraunhofer-Gesellschaft zur Foerderung der angewandten Forschung e.V., Huawei Technologies Co. LTD., +% Koninklijke Philips N.V., Nippon Telegraph and Telephone Corporation, Nokia Technologies Oy, Orange, +% Panasonic Holdings Corporation, Qualcomm Technologies, Inc., VoiceAge Corporation, and other +% contributors to this repository retain full ownership rights in their respective contributions in +% the software. This notice grants no license of any kind, including but not limited to patent +% license, nor is any license granted by implication, estoppel or otherwise. +% +% Contributors are required to enter into the IVAS codec Public Collaboration agreement before making +% contributions. +% +% This software is provided "AS IS", without any express or implied warranties. The software is in the +% development stage. It is intended exclusively for experts who have experience with such software and +% solely for the purpose of inspection. All implied warranties of non-infringement, merchantability +% and fitness for a particular purpose are hereby disclaimed and excluded. +% +% Any dispute, controversy or claim arising under or in relation to providing this software shall be +% submitted to and settled by the final, binding jurisdiction of the courts of Munich, Germany in +% accordance with the laws of the Federal Republic of Germany excluding its conflict of law rules and +% the United Nations Convention on Contracts on the International Sales of Goods. +% +%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% +function IR_cldfb = SHD_2_ROM( py_path, rom_c_file, sofa_file, ambi_order, hrir_len ) +% SHD_2_ROM( python_path, rom_c_file, sofa_file, ambisonics_order, hrir_length ) +% +% - converts sphere-sampled Head Related Impulse Responses (HRIRs) given in sofa_file +% to the Spherical Harmonics domain (SHD) using generate_HOA_HRIRs_MOD_lens.m +% - converts SHD HRIRs to Complex Low Delay Filter Bank (CLDFB) domain using fir_to_cldfb_fir.m +% - writes CLDFB HRIRs to c-code ROM tables. +[thispath,~,~] = fileparts(mfilename('fullpath')); +thispath = [thispath,filesep]; + +%py_path = 'C:\Users\xxxx\AppData\Local\Programs\Python\Python39\python.exe'; % may look like this +if ~exist('sofa_file','var') || isempty(sofa_file) + sofa_file = fullfile(thispath,'..','HRIRs_sofa','HRIR_128_Meth5_IRC_53_Q10_symL_Itrp1_48000.sofa'); +end +if ~exist('ambi_order','var') + ambi_order = 3; +end +if ~exist('hrir_len','var') + hrir_len = 128; +end +%% convert sphere-sampled HRIRs to SHD HRIRs +% requires: +% python -m pip install sofar +% python -m pip install numpy + +% convert sphere-sampled HRIRs to SHD HRIRs +write_out_sofa = 0; +[sofa_path,sofa_name, sofa_ext] = fileparts(sofa_file); +IR = generate_HOA_HRIRs_MOD_lens(ambi_order, py_path, sofa_path, [sofa_name,sofa_ext], hrir_len, write_out_sofa); + +%% SHD -> CLDFB via least squares error optimization +[~,num_ears,num_ch] = size(IR); +num_cldfb_taps = 3; +IR_cldfb = zeros(60,num_cldfb_taps,num_ears,num_ch); % 60 frequency bands +eval_flag = 0; % optional, = 1 requires signal processing toolbox (fftfilt) +legacy_flag = 1; % = 1 used to indicate slightly too short buffers as used to generate tested coefficients +for pos = 1:num_ch + disp(['Channel ',num2str(pos),'/',num2str(num_ch)]) + for ear = 1:num_ears + IR_cldfb(:,:,ear,pos) = fir_to_cldfb_fir( IR(:,ear,pos), num_cldfb_taps, eval_flag, legacy_flag ); + end +end + +%% CLDFB -> ROM +latency_s = 2.083333333333333e-05; % No added latency from conversion method +max_band = 50; % Compute 60 bands, but only use 50 in ROM table + +IR_cldfb_rom = permute(IR_cldfb, [3 1 4 2]); % after permute: [ears(2), bands(60), chans(16), taps(3)] +IR_cldfb_rom = IR_cldfb_rom(:,1:max_band,:,:); +if ambi_order == 1 + order = 'FOA'; +else + order = ['HOA' int2str(ambi_order)]; +end +if ~exist('rom_c_file','var') || isempty(rom_c_file) + rom_c_file = [thispath,'ivas_rom_binauralRenderer_',order,'.c']; % fullfile(thispath,'..','..','..','lib_rend',['ivas_rom_binauralRenderer_',order,'.c']); +end + +fid = fopen(rom_c_file,'wt'); + +fprintf(fid, ['const float FASTCONV_' order '_latency_s = %10.9ff;\n'], latency_s); +write_one_ear( fid, ['const float leftHRIRReal_' order '[BINAURAL_CONVBANDS][' order '_CHANNELS][BINAURAL_NTAPS_SBA]='], real(IR_cldfb_rom(1,:,:,:))); +write_one_ear( fid, ['const float leftHRIRImag_' order '[BINAURAL_CONVBANDS][' order '_CHANNELS][BINAURAL_NTAPS_SBA]='], imag(IR_cldfb_rom(1,:,:,:))); +write_one_ear( fid, ['const float rightHRIRReal_' order '[BINAURAL_CONVBANDS][' order '_CHANNELS][BINAURAL_NTAPS_SBA]='], real(IR_cldfb_rom(2,:,:,:))); +write_one_ear( fid, ['const float rightHRIRImag_' order '[BINAURAL_CONVBANDS][' order '_CHANNELS][BINAURAL_NTAPS_SBA]='], imag(IR_cldfb_rom(2,:,:,:))); + +fclose(fid); + +function write_one_ear( fid, first_line, IR_cldfb_rom) +IR_cldfb_rom = squeeze(IR_cldfb_rom); +[num_bands,num_chans, num_taps] = size(IR_cldfb_rom); +num_spaces = 4; +num_spaces_cur = 0; +fprintf(fid,[first_line,'\n{\n']); +num_spaces_cur = num_spaces_cur + num_spaces; +for band = 1:num_bands + fprintf(fid,[blanks(num_spaces_cur),'{\n']); + num_spaces_cur = num_spaces_cur + num_spaces; + for chan = 1:num_chans + fprintf(fid,[blanks(num_spaces_cur),'{']); + for tap = 1:num_taps + if tap == num_taps + fprintf(fid,'%+ff',IR_cldfb_rom(band,chan,tap)); + else + fprintf(fid,'%+ff, ',IR_cldfb_rom(band,chan,tap)); + end + end + if chan == num_chans + fprintf(fid,'}\n'); + else + fprintf(fid,'},\n'); + end + end + num_spaces_cur = num_spaces_cur - num_spaces; + if band == num_bands + fprintf(fid,[blanks(num_spaces_cur),'}\n']); + else + fprintf(fid,[blanks(num_spaces_cur),'},\n']); + end +end +fprintf(fid,'};\n\n'); diff --git a/scripts/binauralRenderer_interface/matlab_hrir_generation_scripts/convert_SD2SHD_HRIRs.m b/scripts/binauralRenderer_interface/matlab_hrir_generation_scripts/convert_SD2SHD_HRIRs.m new file mode 100644 index 0000000000000000000000000000000000000000..dceee36e4a98e4894b4e319f4a264b7f475b00f3 --- /dev/null +++ b/scripts/binauralRenderer_interface/matlab_hrir_generation_scripts/convert_SD2SHD_HRIRs.m @@ -0,0 +1,64 @@ +%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% +% +% (C) 2022-2023 IVAS codec Public Collaboration with portions copyright Dolby International AB, Ericsson AB, +% Fraunhofer-Gesellschaft zur Foerderung der angewandten Forschung e.V., Huawei Technologies Co. LTD., +% Koninklijke Philips N.V., Nippon Telegraph and Telephone Corporation, Nokia Technologies Oy, Orange, +% Panasonic Holdings Corporation, Qualcomm Technologies, Inc., VoiceAge Corporation, and other +% contributors to this repository. All Rights Reserved. +% +% This software is protected by copyright law and by international treaties. +% The IVAS codec Public Collaboration consisting of Dolby International AB, Ericsson AB, +% Fraunhofer-Gesellschaft zur Foerderung der angewandten Forschung e.V., Huawei Technologies Co. LTD., +% Koninklijke Philips N.V., Nippon Telegraph and Telephone Corporation, Nokia Technologies Oy, Orange, +% Panasonic Holdings Corporation, Qualcomm Technologies, Inc., VoiceAge Corporation, and other +% contributors to this repository retain full ownership rights in their respective contributions in +% the software. This notice grants no license of any kind, including but not limited to patent +% license, nor is any license granted by implication, estoppel or otherwise. +% +% Contributors are required to enter into the IVAS codec Public Collaboration agreement before making +% contributions. +% +% This software is provided "AS IS", without any express or implied warranties. The software is in the +% development stage. It is intended exclusively for experts who have experience with such software and +% solely for the purpose of inspection. All implied warranties of non-infringement, merchantability +% and fitness for a particular purpose are hereby disclaimed and excluded. +% +% Any dispute, controversy or claim arising under or in relation to providing this software shall be +% submitted to and settled by the final, binding jurisdiction of the courts of Munich, Germany in +% accordance with the laws of the Federal Republic of Germany excluding its conflict of law rules and +% the United Nations Convention on Contracts on the International Sales of Goods. +% +%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% +function convert_SD2SHD_HRIRs(python_path, sofa_path, sofa_file, IR_size) + +data_struct = struct.empty(3,0); + +sr = [48000, 32000, 16000]; +sr_short = [48, 32, 16]; +sr_dft_size = [240, 160, 80]; + +% FOA +data_struct(1).IR_data = generate_HOA_HRIRs_MOD_lens(1, python_path, sofa_path, sofa_file, IR_size, 0); +data_struct(1).HOA_name = 'FOA'; +data_struct(1).n_HOA_ch = 4; +data_struct(1).sr = sr; +data_struct(1).sr_short = sr_short; +data_struct(1).sr_dft_size = sr_dft_size; + +% HOA2 +data_struct(2).IR_data = generate_HOA_HRIRs_MOD_lens(2, python_path, sofa_path, sofa_file, IR_size, 0); +data_struct(2).HOA_name = 'HOA2'; +data_struct(2).n_HOA_ch = 9; +data_struct(2).sr = sr; +data_struct(2).sr_short = sr_short; +data_struct(2).sr_dft_size = sr_dft_size; + +% HOA3 +data_struct(3).IR_data = generate_HOA_HRIRs_MOD_lens(3, python_path, sofa_path, sofa_file, IR_size, 0); +data_struct(3).HOA_name = 'HOA3'; +data_struct(3).n_HOA_ch = 16; +data_struct(3).sr = sr; +data_struct(3).sr_short = sr_short; +data_struct(3).sr_dft_size = sr_dft_size; + +generate_rom_tables(data_struct) diff --git a/scripts/binauralRenderer_interface/matlab_hrir_generation_scripts/fir_to_cldfb_fir.m b/scripts/binauralRenderer_interface/matlab_hrir_generation_scripts/fir_to_cldfb_fir.m new file mode 100644 index 0000000000000000000000000000000000000000..5043f4f2c7577e11b070a33c18ae8bcffb6bef5f --- /dev/null +++ b/scripts/binauralRenderer_interface/matlab_hrir_generation_scripts/fir_to_cldfb_fir.m @@ -0,0 +1,304 @@ +%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% +% +% (C) 2022-2023 IVAS codec Public Collaboration with portions copyright Dolby International AB, Ericsson AB, +% Fraunhofer-Gesellschaft zur Foerderung der angewandten Forschung e.V., Huawei Technologies Co. LTD., +% Koninklijke Philips N.V., Nippon Telegraph and Telephone Corporation, Nokia Technologies Oy, Orange, +% Panasonic Holdings Corporation, Qualcomm Technologies, Inc., VoiceAge Corporation, and other +% contributors to this repository. All Rights Reserved. +% +% This software is protected by copyright law and by international treaties. +% The IVAS codec Public Collaboration consisting of Dolby International AB, Ericsson AB, +% Fraunhofer-Gesellschaft zur Foerderung der angewandten Forschung e.V., Huawei Technologies Co. LTD., +% Koninklijke Philips N.V., Nippon Telegraph and Telephone Corporation, Nokia Technologies Oy, Orange, +% Panasonic Holdings Corporation, Qualcomm Technologies, Inc., VoiceAge Corporation, and other +% contributors to this repository retain full ownership rights in their respective contributions in +% the software. This notice grants no license of any kind, including but not limited to patent +% license, nor is any license granted by implication, estoppel or otherwise. +% +% Contributors are required to enter into the IVAS codec Public Collaboration agreement before making +% contributions. +% +% This software is provided "AS IS", without any express or implied warranties. The software is in the +% development stage. It is intended exclusively for experts who have experience with such software and +% solely for the purpose of inspection. All implied warranties of non-infringement, merchantability +% and fitness for a particular purpose are hereby disclaimed and excluded. +% +% Any dispute, controversy or claim arising under or in relation to providing this software shall be +% submitted to and settled by the final, binding jurisdiction of the courts of Munich, Germany in +% accordance with the laws of the Federal Republic of Germany excluding its conflict of law rules and +% the United Nations Convention on Contracts on the International Sales of Goods. +% +%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% +function F = fir_to_cldfb_fir( target_fir, num_cldfb_taps, eval_flag, legacy_flag ) +% F = fir_to_cldfb_fir( target_fir, num_cldfb_taps, eval_flag, legacy_flag ) +% +% computes complex-valued FIR coefficients to be applied in the CLDFB analysis domain +% which approximate time domain filtering with given target FIR coefficients +% in a least squares sense. + +[pFilter, D, S, L] = get_cldfb_filter; +cldfb_delay = D - S + 1; % processing delay given stride +pFilt = sqrt(S)/L * pFilter(:); +N = length(pFilt); +%legacy_flag = 1; % = 1 used to indicate slightly too short buffers as used to generate tested coefficients + +% Filter Bank analysis / synthesis filters +cldfb_mod_mat = exp(1i*((0:N-1)'-D/2) * ((0:L-1)+0.5) * pi/L ); % N x L modulation matrix with alternating sign periodicity of 2*L in time (row) direction +cldfb_ansyn_filters = pFilt .* cldfb_mod_mat; % all (L) analysis/synthesis filters + +target_fir = target_fir(:); +fir_length = length(target_fir); + +%% brute force least squares +if legacy_flag + num_slots = ceil(fir_length/S) + 10; % used to generate tested IVAS coefficients in CLDFB HRIR ROM with slightly too short buffers + inp_len = S * num_slots; + len = inp_len - cldfb_delay; + idx_opt = cldfb_delay + (1:len); + pos_offset = 200; +else + num_slots = ceil((fir_length + 2 * N - 1) / S); % convolution with target fir, analysis + synthesis filters (of length N each) + len = S * num_slots; + inp_len = len; + idx_opt = (1:len); % use all output samples in the optimization, including the filter bank delay + pos_offset = 0; +end + +% Big trial and error +Nbins = num_cldfb_taps * L; +Vreal = zeros(S*len, Nbins); % all real for all shifts concatenated for each mask bin, +Vimag = zeros(S*len, Nbins); % then all imag +r = zeros(S*len,1); % right hand side all shifts concatenated + +x = zeros(inp_len,1); +for pos = 1:S + % input dirac pulse + x(:) = 0; + + % target impulse response y0 for input pulse at pos + out = x; + out(pos+pos_offset+(0:fir_length-1)) = target_fir; + y0 = [zeros(cldfb_delay,1);out(1:end-cldfb_delay)]; + r(len*(pos-1)+(1:len)) = y0(idx_opt); + fprintf('.'); + + % CLDFB analysis for impulse input + x(pos+pos_offset) = 1; + X = zeros(L,num_slots); + frm_idx = (1:N)'; + x_tmp = [zeros(N-S,1);x]; + for slt = 1:num_slots + X(:,slt) = sum(cldfb_ansyn_filters.*flipud(x_tmp(frm_idx))).'; + frm_idx = frm_idx + S; + end + for cldfb_bin = 1:Nbins + cldfb_band = rem(cldfb_bin-1,L)+1; + cldfb_lag = fix((cldfb_bin-1)/L); + + % real-valued filter contribution of a single tap + frm_idx = cldfb_lag * S + (1:N)'; + y_tmp = zeros( num_slots * S + N - S,1); + for slt = cldfb_lag+1:num_slots + y_tmp(frm_idx) = y_tmp(frm_idx) + real(X(cldfb_band,slt-cldfb_lag) * cldfb_ansyn_filters(:,cldfb_band)); + frm_idx = frm_idx + S; + end + Vreal(len*(pos-1)+1:len*pos,cldfb_bin) = y_tmp(idx_opt); + + % imaginary-valued filter contribution of a single tap + frm_idx = cldfb_lag * S + (1:N)'; + y_tmp(:) = 0; + for slt = cldfb_lag+1:num_slots + y_tmp(frm_idx) = y_tmp(frm_idx) + real(1i*X(cldfb_band,slt-cldfb_lag) * cldfb_ansyn_filters(:,cldfb_band)); + frm_idx = frm_idx + S; + end + Vimag(len*(pos-1)+1:len*pos,cldfb_bin) = y_tmp(idx_opt); + end % bin +end % pos +fprintf('\n'); + +% solve lsq +V = [Vreal,Vimag]; +M = V'*V; +Mrel = M + 1e-8*norm(M)*eye(size(M)); +c = Mrel\(V'*r); + +% map back to mask locations +F = zeros(L,num_cldfb_taps); +for cldfb_bin = 1:Nbins + F(cldfb_bin) = c(cldfb_bin)+1i*c(Nbins+cldfb_bin); +end + +if eval_flag + % evaluation (not needed for table generation, needs signal processing toolbox) + get_snr( target_fir, F, cldfb_ansyn_filters, cldfb_delay, S ); +end + +function get_snr( target_fir, cldfb_firs, cldfb_ansyn_filters, cldfb_delay, S ) + +rng(0); +num_samples = 5*48000; +x = randn(num_samples,1); + +% Filter noise in time domain +y1 = fftfilt(target_fir,x); + +% Filter noise in CLDFB domain +[N,L] = size(cldfb_ansyn_filters); +num_slots = ceil((num_samples + cldfb_delay) / S); +% analysis +X = zeros(L,num_slots); +frm_idx = (1:N)'; +x_tmp = [zeros(N-S,1);x;zeros(cldfb_delay,1)]; +for slt = 1:num_slots + X(:,slt) = sum(cldfb_ansyn_filters.*flipud(x_tmp(frm_idx))).'; + frm_idx = frm_idx + S; +end +% filter +X = fftfilt(cldfb_firs.',X.').'; +% synthesis +frm_idx = (1:N)'; +y_tmp = zeros( num_slots * S + N - S,1); +for slt = 1:num_slots + y_tmp(frm_idx) = y_tmp(frm_idx) + real(sum(X(:,slt).' .* cldfb_ansyn_filters,2)); + frm_idx = frm_idx + S; +end +y2 = y_tmp(cldfb_delay+(1:num_samples)); + +plot(y1), hold on, plot(y2,'r--'), plot(y1-y2,'k'), hold off +set(gca,'xlim',[1,num_samples]) +fprintf('SNR: %.1f dB \n', 10*log10(sum(y1.^2)/sum((y1-y2).^2))); +legend('time domain','CLDFB domain','difference') +title('Filtered noise') + +function [h, D, S, L] = get_cldfb_filter() +% const float LDQMF_60[] in \lib_com\rom_com.c, line 5219 +S = 60; % stride +L = 60; % frequency bands +D = 240 + S - 1; % system delay +h = [ + 0.0000953076, 0.0001230230, 0.0001279964, 0.0001260533, 0.0001211219 + 0.0001122123, 0.0001010860, 0.0000876540, 0.0000719883, 0.0000545472 + 0.0000352143, 0.0000145686, -0.0000074264, -0.0000303788, -0.0000539205 + -0.0000782743, -0.0001028838, -0.0001275374, -0.0001520015, -0.0001760167 + -0.0001997108, -0.0002226708, -0.0002446725, -0.0002655797, -0.0002852145 + -0.0003034996, -0.0003203036, -0.0003356283, -0.0003493345, -0.0003614030 + -0.0003719004, -0.0003807641, -0.0003881051, -0.0003939842, -0.0003985357 + -0.0004019095, -0.0004041938, -0.0004056677, -0.0004065430, -0.0004069925 + -0.0004072535, -0.0004075877, -0.0004083676, -0.0004098394, -0.0004122990 + -0.0004160839, -0.0004214063, -0.0004285777, -0.0004378651, -0.0004495422 + -0.0004637682, -0.0004806494, -0.0005003878, -0.0005231378, -0.0005489803 + -0.0005777747, -0.0006095612, -0.0006443121, -0.0006813223, -0.0007226231 + -0.0007722576, -0.0008268412, -0.0008839625, -0.0009417336, -0.0010004630 + -0.0010601623, -0.0011206097, -0.0011817788, -0.0012432419, -0.0013045983 + -0.0013656860, -0.0014260965, -0.0014855355, -0.0015435946, -0.0015999591 + -0.0016543545, -0.0017062968, -0.0017554691, -0.0018015467, -0.0018441341 + -0.0018829798, -0.0019177221, -0.0019480695, -0.0019736972, -0.0019943134 + -0.0020097434, -0.0020197174, -0.0020240925, -0.0020226294, -0.0020152442 + -0.0020017736, -0.0019820682, -0.0019561697, -0.0019240153, -0.0018855907 + -0.0018409232, -0.0017900462, -0.0017330211, -0.0016699535, -0.0016009507 + -0.0015261442, -0.0014456788, -0.0013597424, -0.0012685407, -0.0011722331 + -0.0010710671, -0.0009652392, -0.0008549765, -0.0007405236, -0.0006221444 + -0.0005001140, -0.0003745670, -0.0002458634, -0.0001142541, 0.0000199491 + 0.0001564174, 0.0002949402, 0.0004350246, 0.0005769439, 0.0007203126 + -0.0008803223, -0.0010328424, -0.0011841310, -0.0013346316, -0.0014848098 + -0.0016343417, -0.0017832819, -0.0019316213, -0.0020790498, -0.0022252349 + -0.0023701149, -0.0025136294, -0.0026556554, -0.0027960713, -0.0029348312 + -0.0030717771, -0.0032068293, -0.0033399195, -0.0034709862, -0.0035999804 + -0.0037267797, -0.0038513308, -0.0039736414, -0.0040935921, -0.0042111278 + -0.0043262239, -0.0044388464, -0.0045489701, -0.0046565188, -0.0047614835 + -0.0048637423, -0.0049632201, -0.0050599808, -0.0051539382, -0.0052450863 + -0.0053333500, -0.0054187514, -0.0055012843, -0.0055808770, -0.0056575472 + -0.0057313135, -0.0058021732, -0.0058701355, -0.0059352517, -0.0059975707 + -0.0060571772, -0.0061141332, -0.0061685541, -0.0062205540, -0.0062703062 + -0.0063179093, -0.0063635921, -0.0064075105, -0.0064498796, -0.0064908965 + -0.0065308069, -0.0065698619, -0.0066083665, -0.0066466411, -0.0066849431 + -0.0067233290, -0.0067621553, -0.0068021296, -0.0068436749, -0.0068870094 + -0.0069324085, -0.0069801519, -0.0070305937, -0.0070840055, -0.0071406048 + -0.0072006541, -0.0072644479, -0.0073321410, -0.0074039386, -0.0074799177 + -0.0075602704, -0.0076450342, -0.0077342330, -0.0078278277, -0.0079257628 + -0.0080279401, -0.0081341872, -0.0082442267, -0.0083577875, -0.0084744738 + -0.0085938899, -0.0087156557, -0.0088391500, -0.0089637861, -0.0090888245 + -0.0092134504, -0.0093367994, -0.0094579896, -0.0095760096, -0.0096898535 + -0.0097982995, -0.0099003557, -0.0099947909, -0.0100801717, -0.0101551116 + -0.0102182031, -0.0102678994, -0.0103026126, -0.0103207529, -0.0103206923 + -0.0103006857, -0.0102590285, -0.0101939747, -0.0101036867, -0.0099863587 + -0.0098401112, -0.0096632261, -0.0094537362, -0.0092098210, -0.0089295702 + -0.0086111929, -0.0082527259, -0.0078523541, -0.0074084769, -0.0069190590 + 0.0063841688, 0.0057985946, 0.0051621343, 0.0044734711, 0.0037309236 + 0.0029329660, 0.0020781513, 0.0011651339, 0.0001925042, -0.0008409545 + -0.0019364181, -0.0030950012, -0.0043176264, -0.0056051607, -0.0069584334 + -0.0083780792, -0.0098646237, -0.0114185056, -0.0130400723, -0.0147295250 + -0.0164868534, -0.0183120724, -0.0202049762, -0.0221651513, -0.0241921283 + -0.0262852497, -0.0284437388, -0.0306666382, -0.0329528190, -0.0353010744 + -0.0377098918, -0.0401776619, -0.0427025780, -0.0452826768, -0.0479161367 + -0.0506004691, -0.0533332452, -0.0561118126, -0.0589331910, -0.0617944039 + -0.0646922663, -0.0676232576, -0.0705836788, -0.0735698044, -0.0765774846 + -0.0796026587, -0.0826408863, -0.0856874809, -0.0887378305, -0.0917868018 + -0.0948293805, -0.0978601947, -0.1008738130, -0.1038645208, -0.1068264544 + -0.1097536832, -0.1126400903, -0.1154794544, -0.1182654947, -0.1209914312 + -0.1236500666, -0.1262341589, -0.1287376434, -0.1311538219, -0.1334753036 + -0.1356947273, -0.1378047168, -0.1397978216, -0.1416664869, -0.1434033662 + -0.1450008005, -0.1464512348, -0.1477471888, -0.1488809884, -0.1498452872 + -0.1506324410, -0.1512351334, -0.1516460329, -0.1518578976, -0.1518635303 + -0.1516559124, -0.1512281001, -0.1505732536, -0.1496847868, -0.1485562176 + -0.1471813470, -0.1455538720, -0.1436681300, -0.1415183097, -0.1390990764 + -0.1364052594, -0.1334318966, -0.1301742792, -0.1266280264, -0.1227891371 + -0.1186537445, -0.1142183766, -0.1094799563, -0.1044355705, -0.0990828425 + -0.0934195668, -0.0874440819, -0.0811550021, -0.0745511875, -0.0676321834 + -0.0603975877, -0.0528475679, -0.0449828543, -0.0368040986, -0.0283128861 + -0.0195108838, -0.0104003223, -0.0009837818, 0.0087356847, 0.0187546927 + 0.0290693250, 0.0396753438, 0.0505682528, 0.0617432520, 0.0731955394 + -0.0849232078, -0.0969054326, -0.1091460735, -0.1216373071, -0.1343720406 + -0.1473424733, -0.1605402082, -0.1739567965, -0.1875831038, -0.2014097124 + -0.2154271752, -0.2296251506, -0.2439934313, -0.2585212290, -0.2731975317 + -0.2880111337, -0.3029502928, -0.3180032372, -0.3331578076, -0.3484017253 + -0.3637222052, -0.3791064322, -0.3945416212, -0.4100143015, -0.4255111217 + -0.4410185516, -0.4565227628, -0.4720100164, -0.4874662757, -0.5028775334 + -0.5182296634, -0.5335084200, -0.5486994982, -0.5637886524, -0.5787616372 + -0.5936041474, -0.6083019376, -0.6228409410, -0.6372069120, -0.6513859630 + -0.6653640866, -0.6791275144, -0.6926627755, -0.7059561610, -0.7189947963 + -0.7317654490, -0.7442554235, -0.7564523220, -0.7683438063, -0.7799182534 + -0.7911639810, -0.8020697832, -0.8126249313, -0.8228194118, -0.8326428533 + -0.8420860767, -0.8511404991, -0.8597975969, -0.8680517077, -0.8758881092 + -0.8832823634, -0.8902196884, -0.8967157602, -0.9027729034, -0.9083824754 + -0.9135394692, -0.9182395935, -0.9224776030, -0.9262499809, -0.9295535684 + -0.9323854446, -0.9347436428, -0.9366261959, -0.9380323887, -0.9389615655 + -0.9394137263, -0.9393896461, -0.9388904572, -0.9379178882, -0.9364743829 + -0.9345626831, -0.9321863055, -0.9293491840, -0.9260557890, -0.9223110080 + -0.9181203246, -0.9134896994, -0.9084255695, -0.9029349089, -0.8970250487 + -0.8907034993, -0.8839784265, -0.8768582940, -0.8693521619, -0.8614694476 + -0.8532197475, -0.8446131349, -0.8356599212, -0.8263708353, -0.8167568445 + -0.8068289757, -0.7965991497, -0.7860788107, -0.7752800584, -0.7642148733 + -0.7528960109, -0.7413358092, -0.7295469642, -0.7175422311, -0.7053351402 + -0.6929380894, -0.6803644896, -0.6676273942, -0.6547405124, -0.6417166591 + -0.6285686493, -0.6153115034, -0.6019562483, -0.5885198116, -0.5750215650 + 0.5615197420, 0.5478612781, 0.5341838002, 0.5204906464, 0.5067980289 + 0.4931168854, 0.4794588387, 0.4658361673, 0.4522601366, 0.4387422502 + 0.4252935350, 0.4119254053, 0.3986486793, 0.3854739666, 0.3724119067 + 0.3594728410, 0.3466667533, 0.3340034485, 0.3214924335, 0.3091430366 + 0.2969639599, 0.2849639654, 0.2731511295, 0.2615332901, 0.2501178682 + 0.2389119864, 0.2279221565, 0.2171545923, 0.2066148520, 0.1963084787 + 0.1862401515, 0.1764142811, 0.1668347418, 0.1575049609, 0.1484276950 + 0.1396053135, 0.1310400218, 0.1227331311, 0.1146853194, 0.1068974212 + 0.0993694067, 0.0921007246, 0.0850901082, 0.0783365741, 0.0718384907 + 0.0655927584, 0.0595967993, 0.0538481586, 0.0483424664, 0.0430756323 + 0.0380428955, 0.0332404599, 0.0286619961, 0.0242999699, 0.0201510899 + 0.0162059069, 0.0124559226, 0.0088928537, 0.0054926532, 0.0023052765 + -0.0005515143, -0.0030201224, -0.0052712574, -0.0073737046, -0.0093160523 + -0.0111072771, -0.0127562135, -0.0142635731, -0.0156361461, -0.0168790054 + -0.0179969221, -0.0189934950, -0.0198726747, -0.0206398536, -0.0212980714 + -0.0218509119, -0.0223025978, -0.0226570386, -0.0229178313, -0.0230882075 + -0.0231725387, -0.0231746566, -0.0230979007, -0.0229462404, -0.0227237809 + -0.0224345829, -0.0220820960, -0.0216706358, -0.0212045144, -0.0206875466 + -0.0201238506, -0.0195175279, -0.0188730918, -0.0181944817, -0.0174855441 + -0.0167510118, -0.0159947462, -0.0152208358, -0.0144332750, -0.0136361914 + -0.0128338682, -0.0120294262, -0.0112272501, -0.0104311826, -0.0096443929 + -0.0088709844, -0.0081134979, -0.0073764324, -0.0066623385, -0.0059733889 + -0.0053142183, -0.0046856776, -0.0040914025, -0.0035321070, -0.0030089030 + -0.0025271603, -0.0020749648, -0.0016621647, -0.0012705614, -0.0008115423 +]; +h = h.'; +h = h(:); +h(1*120+(1:120)) = -h(1*120+(1:120)); +h(3*120+(1:120)) = -h(3*120+(1:120)); +h(:) = h(end:-1:1); \ No newline at end of file diff --git a/scripts/binauralRenderer_interface/matlab_hrir_generation_scripts/generate_HOA_HRIRs_MOD_lens.m b/scripts/binauralRenderer_interface/matlab_hrir_generation_scripts/generate_HOA_HRIRs_MOD_lens.m new file mode 100644 index 0000000000000000000000000000000000000000..ff42cf79c45cb1e02022e153d34463ae72e1555e --- /dev/null +++ b/scripts/binauralRenderer_interface/matlab_hrir_generation_scripts/generate_HOA_HRIRs_MOD_lens.m @@ -0,0 +1,141 @@ +%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% +% +% (C) 2022-2023 IVAS codec Public Collaboration with portions copyright Dolby International AB, Ericsson AB, +% Fraunhofer-Gesellschaft zur Foerderung der angewandten Forschung e.V., Huawei Technologies Co. LTD., +% Koninklijke Philips N.V., Nippon Telegraph and Telephone Corporation, Nokia Technologies Oy, Orange, +% Panasonic Holdings Corporation, Qualcomm Technologies, Inc., VoiceAge Corporation, and other +% contributors to this repository. All Rights Reserved. +% +% This software is protected by copyright law and by international treaties. +% The IVAS codec Public Collaboration consisting of Dolby International AB, Ericsson AB, +% Fraunhofer-Gesellschaft zur Foerderung der angewandten Forschung e.V., Huawei Technologies Co. LTD., +% Koninklijke Philips N.V., Nippon Telegraph and Telephone Corporation, Nokia Technologies Oy, Orange, +% Panasonic Holdings Corporation, Qualcomm Technologies, Inc., VoiceAge Corporation, and other +% contributors to this repository retain full ownership rights in their respective contributions in +% the software. This notice grants no license of any kind, including but not limited to patent +% license, nor is any license granted by implication, estoppel or otherwise. +% +% Contributors are required to enter into the IVAS codec Public Collaboration agreement before making +% contributions. +% +% This software is provided "AS IS", without any express or implied warranties. The software is in the +% development stage. It is intended exclusively for experts who have experience with such software and +% solely for the purpose of inspection. All implied warranties of non-infringement, merchantability +% and fitness for a particular purpose are hereby disclaimed and excluded. +% +% Any dispute, controversy or claim arising under or in relation to providing this software shall be +% submitted to and settled by the final, binding jurisdiction of the courts of Munich, Germany in +% accordance with the laws of the Federal Republic of Germany excluding its conflict of law rules and +% the United Nations Convention on Contracts on the International Sales of Goods. +% +%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% +function IR_data = generate_HOA_HRIRs_MOD_lens(order, python_path, sofa_path, sofa_file_name, ir_len, write_out_sofa) + % HRIR convertor - Takes sphere sampled HRIRs and converts them to + % HOA HRIRs. + % + % order - HOA order to be converted to. + % python_path - path that points to the python binary that has the + % sofar python module installed. sofar is needed to read the .sofa + % files that contain the HRIRs. + % sofa_path - path to the directory that contains the sofa files to be + % converted. + % sofa_file - file name of the HRTFs to be converted + % ir_len - length of the IRs to be used. + % write_out_sofa - boolean to enable/disable writing out of sofa + % + % Typical usage: + % generate_HOA_HRIRs_MOD_lens(1, "~/.pyenv/versions/3.8.16/bin/python", '~/git/ivas-pc-testfiles/sofa-files/', 'HRIR_128_48000.sofa', 128) + % + +% Pointing to the correct python binary +py_env = pyenv; +if ~strcmp(py_env.Executable, python_path) + pyenv('Version', python_path); +end + +% Load in the support coefs +load('hrtf_support_coefs.mat', 'hrtf_support_coefs'); +rmsSphere = hrtf_support_coefs(order).rmsSphere; +LR_Odd = hrtf_support_coefs(order).LR_Odd; +XYZ_to_Pan = hrtf_support_coefs(order).XYZ_to_Pan; +AllPass = hrtf_support_coefs(order).AP; + +% Choose a hi-res set of positions to sample the input HRTFs +Vs_hi_res = load("sphere_packing_2562.mat"); +Vs_hi_res = Vs_hi_res.Vs_hi_res; +N = 512; + +% Fetch the HRTFs, and figure out the ITD for every direction +% name = ['HRIR_', num2str(ir_len), '_Meth5_IRC_53_Q10_symL_Itrp1_48000']; +H = hrtf_library_loader(); +H.readSOFA(char(fullfile(sofa_path, sofa_file_name))); + +IRs_hi_res = H.XYZ_to_IR(Vs_hi_res); +FRs_hi_res = m_dft(IRs_hi_res, N); % freq x ears x Vs + +FRs_hi_res_minP = mag2min_phase(FRs_hi_res); + +Excess_Phase = squeeze(unwrap(diff(angle(FRs_hi_res), 1, 2) - ... + diff(angle(FRs_hi_res_minP), 1, 2))); + +bin1200 = ceil( 1200/24000*N ); +ITD_hi_res = Excess_Phase(bin1200,:)' / ((bin1200-0.5)/N*24000*2*pi); + +MaxDel = max(ITD_hi_res, [], 'all'); + +% Create 2 ears +Ear_dels_hi_res = (repmat(ITD_hi_res, 1, 2) .* [0.5 -0.5]) + 0.5*MaxDel .* ... + (1 - 2/pi*cos(ITD_hi_res * pi/2 / MaxDel)); + +MRs_hi_res = abs(FRs_hi_res_minP); +% Generate permutation +[~, perm] = ismembertol(... + Vs_hi_res', Vs_hi_res'.*[1,-1,1], ... + 1e-4, "ByRows",true); +MRs_hi_res(:,2,:) = MRs_hi_res(:,2,perm); + +New_FreqResp_L = mag2min_phase(squeeze(mean(MRs_hi_res, 2))) .* ... + m_dft(get_allpass_IRs(AllPass, Ear_dels_hi_res(:, 1) * 48000), N, 1); + +% Create solving weights +weights = abs(New_FreqResp_L); +weights(weights < 0.1) = 0.1; +weights = 1 ./ (sqrt(sqrt(weights))); + +% Solve to compute the HOA frequency responses. +[m, ~] = size(XYZ_to_Pan); +FreqResp_HOA = zeros(m, N); +for k=1:N + Aw = New_FreqResp_L(k,:) .* weights(k,:); + Bw = XYZ_to_Pan .* weights(k,:); + + FreqResp_HOA(:,k) = Aw * pinv(Bw, 0); +end + +FreqResp_HOA_abs2 = real(FreqResp_HOA.*conj(FreqResp_HOA)); +FreqResp_HOA = FreqResp_HOA .* mag2min_phase(((FreqResp_HOA_abs2' * rmsSphere.^2) .^ (-0.5)), 1).'; +% Convert back to IRs +IR_HOA = m_idft(FreqResp_HOA.', [], 1); + +IR_HOA = cat(3, IR_HOA, (IR_HOA .* (1-2*LR_Odd)')); +% Put matrix dimensions in the right order +IR_HOA = permute(IR_HOA, [3, 1, 2]); + +% Get the IRs to the right length +IR_HOA = IR_HOA(:,1:ir_len,:) .* sin(interp1([0,150/192*ir_len,ir_len+1],[1,1,0]*pi/2, 1:ir_len)); + +IR = permute(IR_HOA, [2, 1, 3]); +HOAformat_str = ['HOA',num2str(order),'S']; +save(fullfile(erase(sofa_file_name, ".sofa") + "_converted_" + HOAformat_str + ".mat"), "IR") + +if (write_out_sofa == 1) + H.updateSOFA(IR); + H.writeSOFA(fullfile(erase(sofa_file_name, ".sofa") + "_converted_" + HOAformat_str + ".sofa")); +end + +IR_data = IR; + + + + + diff --git a/scripts/binauralRenderer_interface/matlab_hrir_generation_scripts/generate_rom_tables.m b/scripts/binauralRenderer_interface/matlab_hrir_generation_scripts/generate_rom_tables.m new file mode 100644 index 0000000000000000000000000000000000000000..4f0e959c6d783ae3c0b997909064d00de209573f --- /dev/null +++ b/scripts/binauralRenderer_interface/matlab_hrir_generation_scripts/generate_rom_tables.m @@ -0,0 +1,253 @@ +%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% +% +% (C) 2022-2023 IVAS codec Public Collaboration with portions copyright Dolby International AB, Ericsson AB, +% Fraunhofer-Gesellschaft zur Foerderung der angewandten Forschung e.V., Huawei Technologies Co. LTD., +% Koninklijke Philips N.V., Nippon Telegraph and Telephone Corporation, Nokia Technologies Oy, Orange, +% Panasonic Holdings Corporation, Qualcomm Technologies, Inc., VoiceAge Corporation, and other +% contributors to this repository. All Rights Reserved. +% +% This software is protected by copyright law and by international treaties. +% The IVAS codec Public Collaboration consisting of Dolby International AB, Ericsson AB, +% Fraunhofer-Gesellschaft zur Foerderung der angewandten Forschung e.V., Huawei Technologies Co. LTD., +% Koninklijke Philips N.V., Nippon Telegraph and Telephone Corporation, Nokia Technologies Oy, Orange, +% Panasonic Holdings Corporation, Qualcomm Technologies, Inc., VoiceAge Corporation, and other +% contributors to this repository retain full ownership rights in their respective contributions in +% the software. This notice grants no license of any kind, including but not limited to patent +% license, nor is any license granted by implication, estoppel or otherwise. +% +% Contributors are required to enter into the IVAS codec Public Collaboration agreement before making +% contributions. +% +% This software is provided "AS IS", without any express or implied warranties. The software is in the +% development stage. It is intended exclusively for experts who have experience with such software and +% solely for the purpose of inspection. All implied warranties of non-infringement, merchantability +% and fitness for a particular purpose are hereby disclaimed and excluded. +% +% Any dispute, controversy or claim arising under or in relation to providing this software shall be +% submitted to and settled by the final, binding jurisdiction of the courts of Munich, Germany in +% accordance with the laws of the Federal Republic of Germany excluding its conflict of law rules and +% the United Nations Convention on Contracts on the International Sales of Goods. +% +%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% +function generate_rom_tables(data_struct) + + h_file_name = 'ivas_rom_binaural_crend_head.h'; + c_file_name = 'ivas_rom_binaural_crend_head.c'; + tab = ' '; + + copyright_str = string(join({ + '/******************************************************************************************************' + '' + ' (C) 2022-2023 IVAS codec Public Collaboration with portions copyright Dolby International AB, Ericsson AB,' + ' Fraunhofer-Gesellschaft zur Foerderung der angewandten Forschung e.V., Huawei Technologies Co. LTD.,' + ' Koninklijke Philips N.V., Nippon Telegraph and Telephone Corporation, Nokia Technologies Oy, Orange,' + ' Panasonic Holdings Corporation, Qualcomm Technologies, Inc., VoiceAge Corporation, and other' + ' contributors to this repository. All Rights Reserved.' + '' + ' This software is protected by copyright law and by international treaties.' + ' The IVAS codec Public Collaboration consisting of Dolby International AB, Ericsson AB,' + ' Fraunhofer-Gesellschaft zur Foerderung der angewandten Forschung e.V., Huawei Technologies Co. LTD.,' + ' Koninklijke Philips N.V., Nippon Telegraph and Telephone Corporation, Nokia Technologies Oy, Orange,' + ' Panasonic Holdings Corporation, Qualcomm Technologies, Inc., VoiceAge Corporation, and other' + ' contributors to this repository retain full ownership rights in their respective contributions in' + ' the software. This notice grants no license of any kind, including but not limited to patent' + ' license, nor is any license granted by implication, estoppel or otherwise.' + '' + ' Contributors are required to enter into the IVAS codec Public Collaboration agreement before making' + ' contributions.' + '' + ' This software is provided "AS IS", without any express or implied warranties. The software is in the' + ' development stage. It is intended exclusively for experts who have experience with such software and' + ' solely for the purpose of inspection. All implied warranties of non-infringement, merchantability' + ' and fitness for a particular purpose are hereby disclaimed and excluded.' + '' + ' Any dispute, controversy or claim arising under or in relation to providing this software shall be' + ' submitted to and settled by the final, binding jurisdiction of the courts of Munich, Germany in' + ' accordance with the laws of the Federal Republic of Germany excluding its conflict of law rules and' + ' the United Nations Convention on Contracts on the International Sales of Goods.' + '' + '*******************************************************************************************************/' + '' + '/* clang-format off */' + '' + '/*-------------------------------------------------------------------------' + ' * Binaural rendering related ROM tables' + ' *------------------------------------------------------------------------*/' + '' + '/* Binaural rendering data set based on HRIRs */' + '/* Tables generated by scripts/binauralRenderer_interface/generate_cren_ivas_tables.c, see mixer_conv_sofa_to_rom_table_converter_readme.txt */' + '/* Can be replaced by your own generated HRIR or BRIR tables */' + '' + '' + }, newline)); + + %% Write out the header + fileID = fopen(h_file_name,'w'); + fprintf(fileID,'%s', copyright_str); + + h_file_content = string(join({ + '' + '' + '#ifndef _IVAS_ROM_BINAURAL_CREND_HEAD_' + '#define _IVAS_ROM_BINAURAL_CREND_HEAD_' + '' + '#include "ivas_cnst.h"' + '' + '' + }, newline)); + fprintf(fileID,'%s', h_file_content); + + for hoa_idx=1:length(data_struct) + HOA_name = data_struct(hoa_idx).HOA_name; + n_HOA_ch = data_struct(hoa_idx).n_HOA_ch; + + h_file_content = string(join({ + sprintf('\n/********************** CRendBin_%s_HRIR **********************/', HOA_name) + '' + sprintf('extern float CRendBin_%s_HRIR_latency_s;', HOA_name) + '' + }, newline)); + fprintf(fileID,'%s', h_file_content); + for sr_idx=1:length(data_struct(hoa_idx).sr) + sr = data_struct(hoa_idx).sr(sr_idx); + sr_short = data_struct(hoa_idx).sr_short(sr_idx); + sr_dft_size = data_struct(hoa_idx).sr_dft_size(sr_idx); + + fprintf(fileID, '\n/* Sample Rate = %d */\n\n', sr); + fprintf(fileID, 'extern int16_t CRendBin_%s_HRIR_max_num_iterations_%dkHz;\n', HOA_name, sr_short); + fprintf(fileID, 'extern uint16_t CRendBin_%s_HRIR_num_iterations_%dkHz[%d][BINAURAL_CHANNELS];\n', HOA_name, sr_short, n_HOA_ch); + fprintf(fileID, 'extern uint16_t CRendBin_%s_HRIR_num_iterations_diffuse_%dkHz[BINAURAL_CHANNELS];\n', HOA_name, sr_short); + fprintf(fileID, 'extern uint16_t CRendBin_%s_HRIR_pIndex_frequency_max_%dkHz[%d][BINAURAL_CHANNELS][1];\n', HOA_name, sr_short, n_HOA_ch); + fprintf(fileID, 'extern uint16_t CRendBin_%s_HRIR_index_frequency_max_diffuse_%dkHz;\n', HOA_name, sr_short); + fprintf(fileID, 'extern float CRendBin_%s_HRIR_inv_diffuse_weight_%dkHz[%d];\n', HOA_name, sr_short, n_HOA_ch); + fprintf(fileID, 'extern uint16_t *CRendBin_%s_HRIR_pIndex_frequency_max_diffuse_%dkHz[BINAURAL_CHANNELS];\n', HOA_name, sr_short); + fprintf(fileID, 'extern float CRendBin_%s_HRIR_coeff_re_%dkHz[%d][BINAURAL_CHANNELS][%d];\n', HOA_name, sr_short, n_HOA_ch, sr_dft_size); + fprintf(fileID, 'extern float CRendBin_%s_HRIR_coeff_im_%dkHz[%d][BINAURAL_CHANNELS][%d];\n', HOA_name, sr_short, n_HOA_ch, sr_dft_size); + fprintf(fileID, 'extern float *CRendBin_%s_HRIR_coeff_diffuse_re_%dkHz[BINAURAL_CHANNELS];\n', HOA_name, sr_short); + fprintf(fileID, 'extern float *CRendBin_%s_HRIR_coeff_diffuse_im_%dkHz[BINAURAL_CHANNELS];\n', HOA_name, sr_short); + end + end + fprintf(fileID,'%s\n', '#endif /* _IVAS_ROM_BINAURAL_CREND_HEAD_ */'); + fclose(fileID); + + %% Write out the c file + fileID = fopen(c_file_name,'w'); + fprintf(fileID,'%s', copyright_str); + + c_file_content = string(join({ + '' + '' + '#include ' + '#include "ivas_cnst.h"' + '' + '/* clang-format off */' + '' + '#define WMC_TOOL_SKIP' + '' + '' + '' + }, newline)); + fprintf(fileID,'%s', c_file_content); + + for hoa_idx=1:length(data_struct) + HOA_name = data_struct(hoa_idx).HOA_name; + n_HOA_ch = data_struct(hoa_idx).n_HOA_ch; + + % First Sample rate must be 48k + assert(data_struct(hoa_idx).sr(1) == 48000) + + % Force Latency of HRIRs to 1 sample to retain bit exactness + latency = 1 / data_struct(hoa_idx).sr(1); + latency = latency + 0.000000001; + + fprintf(fileID,'\n/********************** CRendBin_%s_HRIR **********************/\n\n', HOA_name); + fprintf(fileID,'const float CRendBin_%s_HRIR_latency_s = %10.9ff;\n', HOA_name, latency); + + % 48k data + resampled_data{1} = data_struct(hoa_idx).IR_data; + % Resample to 32k + resampled_data{2} = resample(data_struct(hoa_idx).IR_data, 2, 3, 128, 'Dimension', 1); + % Resample to 16k + resampled_data{3} = resample(data_struct(hoa_idx).IR_data, 1, 3, 128, 'Dimension', 1); + for sr_idx=1:length(data_struct(hoa_idx).sr) + sr = data_struct(hoa_idx).sr(sr_idx); + sr_short = data_struct(hoa_idx).sr_short(sr_idx); + sr_dft_size = data_struct(hoa_idx).sr_dft_size(sr_idx); + + FR = m_dft(resampled_data{1}, data_struct(hoa_idx).sr_dft_size(sr_idx), 1); + + fprintf(fileID, '\n'); + fprintf(fileID, '/* Sample Rate = %d */\n', sr); + fprintf(fileID, '\n'); + fprintf(fileID, 'const int16_t CRendBin_%s_HRIR_max_num_iterations_%dkHz = 1;\n', HOA_name, sr_short); + fprintf(fileID, 'const uint16_t CRendBin_%s_HRIR_num_iterations_%dkHz[%d][BINAURAL_CHANNELS]={', HOA_name, sr_short, n_HOA_ch); + c_file_content = ''; + for ch_idx=1:n_HOA_ch + c_file_content = [c_file_content, '{1, 1}, ']; + end + c_file_content(end-1:end) = []; + fprintf(fileID, '%s', c_file_content); + fprintf(fileID, ' };\n'); + fprintf(fileID, 'const uint16_t CRendBin_%s_HRIR_num_iterations_diffuse_%dkHz[BINAURAL_CHANNELS] = {0, 0};\n', HOA_name, sr_short); + fprintf(fileID, 'const uint16_t CRendBin_%s_HRIR_pIndex_frequency_max_%dkHz[%d][BINAURAL_CHANNELS][1]={', HOA_name, sr_short, n_HOA_ch); + c_file_content = ''; + for ch_idx=1:n_HOA_ch + c_file_content = [c_file_content '{']; + for ear_idx=1:2 + c_file_content = [c_file_content, sprintf('{%d},', sr_dft_size)]; + end + c_file_content(end) = []; + c_file_content = [c_file_content '},']; + end + c_file_content(end) = []; + fprintf(fileID, '%s', c_file_content); + fprintf(fileID, '};\n'); + fprintf(fileID, 'const uint16_t CRendBin_%s_HRIR_index_frequency_max_diffuse_%dkHz = 0;\n', HOA_name, sr_short); + fprintf(fileID, 'const float CRendBin_%s_HRIR_inv_diffuse_weight_%dkHz[%d]={', HOA_name, sr_short, n_HOA_ch); + c_file_content = ''; + for ch_idx=1:data_struct(hoa_idx).n_HOA_ch + c_file_content = [c_file_content, sprintf('%ff, ', 0.0)]; + end + c_file_content(end-1:end) = []; + fprintf(fileID, '%s', c_file_content); + fprintf(fileID, '};\n'); + fprintf(fileID, 'const uint16_t *CRendBin_%s_HRIR_pIndex_frequency_max_diffuse_%dkHz[BINAURAL_CHANNELS]={NULL,NULL};\n', HOA_name, sr_short); + + % Real Part + write_out_HRIR(fileID, tab, n_HOA_ch, 2, sr_dft_size, FR, @real, 're', HOA_name, sr_short); + % Imag Part + write_out_HRIR(fileID, tab, n_HOA_ch, 2, sr_dft_size, FR, @imag, 'im', HOA_name, sr_short); + + fprintf(fileID, 'const float *CRendBin_%s_HRIR_coeff_diffuse_re_%dkHz[BINAURAL_CHANNELS]={NULL,NULL};\n', HOA_name, sr_short); + fprintf(fileID, 'const float *CRendBin_%s_HRIR_coeff_diffuse_im_%dkHz[BINAURAL_CHANNELS]={NULL,NULL};\n', HOA_name, sr_short); + end + end + fprintf(fileID, '\n#undef WMC_TOOL_SKIP\n\n'); + + fclose(fileID); + +end + +function write_out_HRIR(fileID, tab, n_hoa_ch, n_ears, dft_size, freq_resp, func_handle, func_name, hoa_name, sr_short) + fprintf(fileID,'%s\n', sprintf('const float CRendBin_%s_HRIR_coeff_%s_%dkHz[%d][BINAURAL_CHANNELS][%d]={', hoa_name, func_name, sr_short, n_hoa_ch, dft_size)); + for hoa_idx=1:n_hoa_ch + fprintf(fileID,'%s%s\n', tab, '{'); + for ear_idx=1:n_ears + fprintf(fileID,'%s%s%s', tab, tab, '{'); + for bin_idx=1:dft_size + fprintf(fileID,'%ff, ', func_handle(freq_resp(bin_idx, ear_idx, hoa_idx))); + if(mod(bin_idx, 96) == 0) + fprintf(fileID,'\n%s%s', tab, tab); + end + end + fseek(fileID, -2, "cof"); + fprintf(fileID,'%s,\n', '}'); + end + fseek(fileID, -2, "cof"); + fprintf(fileID,'\n%s%s,\n', tab, '}'); + end + fseek(fileID, -2, "cof"); + fprintf(fileID,'\n%s;\n', '}'); +end + + diff --git a/scripts/binauralRenderer_interface/matlab_hrir_generation_scripts/get_allpass_IRs.m b/scripts/binauralRenderer_interface/matlab_hrir_generation_scripts/get_allpass_IRs.m new file mode 100644 index 0000000000000000000000000000000000000000..ed4ef76faf21d52a47e055bc9169f758ea812679 --- /dev/null +++ b/scripts/binauralRenderer_interface/matlab_hrir_generation_scripts/get_allpass_IRs.m @@ -0,0 +1,74 @@ +%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% +% +% (C) 2022-2023 IVAS codec Public Collaboration with portions copyright Dolby International AB, Ericsson AB, +% Fraunhofer-Gesellschaft zur Foerderung der angewandten Forschung e.V., Huawei Technologies Co. LTD., +% Koninklijke Philips N.V., Nippon Telegraph and Telephone Corporation, Nokia Technologies Oy, Orange, +% Panasonic Holdings Corporation, Qualcomm Technologies, Inc., VoiceAge Corporation, and other +% contributors to this repository. All Rights Reserved. +% +% This software is protected by copyright law and by international treaties. +% The IVAS codec Public Collaboration consisting of Dolby International AB, Ericsson AB, +% Fraunhofer-Gesellschaft zur Foerderung der angewandten Forschung e.V., Huawei Technologies Co. LTD., +% Koninklijke Philips N.V., Nippon Telegraph and Telephone Corporation, Nokia Technologies Oy, Orange, +% Panasonic Holdings Corporation, Qualcomm Technologies, Inc., VoiceAge Corporation, and other +% contributors to this repository retain full ownership rights in their respective contributions in +% the software. This notice grants no license of any kind, including but not limited to patent +% license, nor is any license granted by implication, estoppel or otherwise. +% +% Contributors are required to enter into the IVAS codec Public Collaboration agreement before making +% contributions. +% +% This software is provided "AS IS", without any express or implied warranties. The software is in the +% development stage. It is intended exclusively for experts who have experience with such software and +% solely for the purpose of inspection. All implied warranties of non-infringement, merchantability +% and fitness for a particular purpose are hereby disclaimed and excluded. +% +% Any dispute, controversy or claim arising under or in relation to providing this software shall be +% submitted to and settled by the final, binding jurisdiction of the courts of Munich, Germany in +% accordance with the laws of the Federal Republic of Germany excluding its conflict of law rules and +% the United Nations Convention on Contracts on the International Sales of Goods. +% +%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% +function IR = get_allpass_IRs(AllPass, dels) + XSet = poly2XSet(AllPass.proto_poly, (dels(:)'-16)*AllPass.upper_freq/AllPass.proto_bw); + + IR = makeIRs(XSet, AllPass.upper_freq); +end + +function XSet = poly2XSet(p, dels_samples) + XSet = zeros(size(p,1), numel(dels_samples)); + for k = 1:size(p,1) + XSet(k,:) = polyval(p(k,:),dels_samples(:)'/32); + end +end + +function irs = makeIRs(X,bw) + irs = map_poles2irs( map2poles(X,bw) ); +end + +function p = map2poles(X,bw) + p = map_2_s_poles(X) * bw; + for k = 1:size(p(:,:),2) + p(:,k) = bilinear(p(:,k),p(:,k),1,48000,bw); + end +end + +function irs = map_poles2irs(p) + irs = zeros(512,size(p(:,:),2)); + for k = 1:size(irs,2) + [~,A] = zp2tf(p(:,k),p(:,k),1); + irs(:,k) = filter(fliplr(A),A, [1;zeros(511,1)]); + end +end + +function p = map_2_s_poles(X) + Order = size(X,1); + if Order==0 + p=[]; + elseif Order==1 + p=-2*pi*X; + else + ang = atan(X(2,:))/2+pi/4; + p = [-2*pi*X(1,:).*exp(1i*[1;-1]*ang) ; map_2_s_poles(X(3:end,:))]; + end +end \ No newline at end of file diff --git a/scripts/binauralRenderer_interface/matlab_hrir_generation_scripts/hrtf_library_loader.m b/scripts/binauralRenderer_interface/matlab_hrir_generation_scripts/hrtf_library_loader.m new file mode 100644 index 0000000000000000000000000000000000000000..f08016036250fdeff472706adb532cbbb6c771f0 --- /dev/null +++ b/scripts/binauralRenderer_interface/matlab_hrir_generation_scripts/hrtf_library_loader.m @@ -0,0 +1,478 @@ +%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% +% +% (C) 2022-2023 IVAS codec Public Collaboration with portions copyright Dolby International AB, Ericsson AB, +% Fraunhofer-Gesellschaft zur Foerderung der angewandten Forschung e.V., Huawei Technologies Co. LTD., +% Koninklijke Philips N.V., Nippon Telegraph and Telephone Corporation, Nokia Technologies Oy, Orange, +% Panasonic Holdings Corporation, Qualcomm Technologies, Inc., VoiceAge Corporation, and other +% contributors to this repository. All Rights Reserved. +% +% This software is protected by copyright law and by international treaties. +% The IVAS codec Public Collaboration consisting of Dolby International AB, Ericsson AB, +% Fraunhofer-Gesellschaft zur Foerderung der angewandten Forschung e.V., Huawei Technologies Co. LTD., +% Koninklijke Philips N.V., Nippon Telegraph and Telephone Corporation, Nokia Technologies Oy, Orange, +% Panasonic Holdings Corporation, Qualcomm Technologies, Inc., VoiceAge Corporation, and other +% contributors to this repository retain full ownership rights in their respective contributions in +% the software. This notice grants no license of any kind, including but not limited to patent +% license, nor is any license granted by implication, estoppel or otherwise. +% +% Contributors are required to enter into the IVAS codec Public Collaboration agreement before making +% contributions. +% +% This software is provided "AS IS", without any express or implied warranties. The software is in the +% development stage. It is intended exclusively for experts who have experience with such software and +% solely for the purpose of inspection. All implied warranties of non-infringement, merchantability +% and fitness for a particular purpose are hereby disclaimed and excluded. +% +% Any dispute, controversy or claim arising under or in relation to providing this software shall be +% submitted to and settled by the final, binding jurisdiction of the courts of Munich, Germany in +% accordance with the laws of the Federal Republic of Germany excluding its conflict of law rules and +% the United Nations Convention on Contracts on the International Sales of Goods. +% +%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% +classdef hrtf_library_loader < handle + % HRTF Library Loader and Interpolator + % + % Typical usage: + % H1 = hrtf_library_loader('CIPIC_137_cleaned', 48000); + % H1.setLen(64); + % IR_L = H1.getHRTF_L([1;1;0]); + % IR_R = H1.getHRTF_R([1;1;0]); + + properties + Len=0; % The currently selected HRTF length + Window=[]; % The window used to set the HRTF length + Lib_Name = []; % The name of the library (as supplied to the class constructor + Lib_SampleRate = []; % The sample rate of the library + SourceRadius = []; % The source distance for HRTFs in this library + PlaybackRadius = []; % The radius for HRTFs during playback + Info = {''}; % A cell array of strings with info about this library + end + properties (Access=protected) + % Define the set of uni-vectors for the HRTF directions in a discrete Library + % If our library is Discrete, we will always specifiy a set of unit-vectors: + Discrete_UnitVectors = []; + % A Discrete Library is generally defined by a set of HRTF Impulse Responses + Discrete_HRTFs = []; + % Define the HRTFs in terms of GroupDelay and Power spectrum. + % Sometimes we will convert the impulse responses into + % GroupDelay/Magnitude format, to allow a better mode of interpolation. + % Sometimes, our library already contains these GD_Mag responses, so we + % don't need to generate them in this class + Discrete_HRTFs_GD_Mag = []; + + % Sometimes we aim to get a smoother group-delay (as a function of + % angle) by using a simplified model. If the following member variable + % is not empty, then we will use this model to compute the group-delay + GD_Model=[]; + + % For the GD_Mag interpolation method, we need to define the FFT length + GD_Mag_FFTLen=[]; + % For the GD_Mag interpolation method, we need to define the freq bands + GD_Mag_Band2FR=[]; + % For the GD_Mag interpolation method, low-pass filter applied to imp-resps + GD_Mag_LoPass=[]; + + % Keep the most recntly computed HRTFs: + % To save us the pain of recalcuating a whole lot of stuff when we + % make consecutive calls the get the Left and Right HRTFs, we store the + % most recent HRTFs, and their associated unit-vectors. + Last_HRTFs=[]; + % Remember the unit vectors for the most recent HRTFs computed + Last_UV=[]; + + Sofa_File = []; + end + + methods + + function obj = hrtf_library_loader() + end + + function readSOFA(obj, Lib_Name) + + % Read in the sofa file using 'sofar' python module + sofa_file = py.sofar.read_sofa(Lib_Name); + + % This looks like a library built from a bunch of discrete + % locations in free field + % Let's figure out the source locations used in the HRTF library: + Pos = hrtf_library_loader.convert_numpy_2darray(sofa_file.SourcePosition)'; + Units = strtrim(strsplit(string(sofa_file.SourcePosition_Units), ',')); + + assert( any(strcmpi(Units{1}, {'degree','radian'})), 'Unknown units'); + if strcmpi(Units{1},'degree'), Pos(1,:)=Pos(1,:)*pi/180; end + assert( any(strcmpi(Units{2}, {'degree','radian'})), 'Unknown units'); + if strcmpi(Units{2},'degree'), Pos(2,:)=Pos(2,:)*pi/180; end + assert( any(strcmpi(Units{3}, {'metre','meter','inch'})), 'Unknown units'); + if strcmpi(Units{3},'inch' ), Pos(3,:)=Pos(3,:)*0.0254; end + Pos = Pos(3,:) .* ... + [cos(Pos(2,:)).*[cos(Pos(1,:));sin(Pos(1,:))];sin(Pos(2,:))]; + + Radii = sqrt(sum(Pos.*Pos,1)); + XYZ = make_unit_vectors(Pos); + NumLoc = size(XYZ,2); + + % Now, get the impulse repsonses: + assert( isprop(sofa_file, 'GLOBAL_DataType'), 'Expected field: GLOBAL_DataType'); + Data.IR = hrtf_library_loader.convert_numpy_3darray(sofa_file.Data_IR); + Data.SamplingRate = double(sofa_file.Data_SamplingRate); + Data.SamplingRate_Units = string(sofa_file.Data_SamplingRate_Units); + Data.Delay = hrtf_library_loader.convert_numpy_2darray(sofa_file.Data_Delay); + switch lower(string(sofa_file.GLOBAL_DataType)) + case 'fir' + assert( size(Data.IR,2)>=2, 'Expecting 2 receivers (ears)'); + if size(Data.IR,2)>2 + Data.IR = Data.IR(:,1:2,:); + Data.Delay = Data.Delay(:,1:2); + end + assert( size(Data.IR,1)==NumLoc, 'IR is incorrect size'); + IRLen = size(Data.IR,3); + assert( all(diff(Data.Delay,1,2)==0), ... + 'Non-zero inter-aural delay offset not (yet) implemented'); + HRTF_L = reshape( Data.IR(:,1,:), NumLoc, IRLen)'; + HRTF_R = reshape( Data.IR(:,2,:), NumLoc, IRLen)'; + Fs = Data.SamplingRate; + assert(strcmp(Data.SamplingRate_Units,'hertz'), 'unknown samplerate units'); + otherwise + error(['unknown GLOBAL_DataType = ', hrtf.GLOBAL_DataType]); + end + + % Normalise for unity gain around 1kHz + FFTLen=1024; + Mid_Bin_1K=round(FFTLen*1000/Fs); + FR_L=fft(HRTF_L,FFTLen); + FR_R=fft(HRTF_R,FFTLen); + MeanG=sqrt(mean(mean(abs([FR_L(Mid_Bin_1K+(-1:1),:),FR_R(Mid_Bin_1K+(-1:1),:)]).^2,1),2)); + HRTF_L = HRTF_L/MeanG; + HRTF_R = HRTF_R/MeanG; + + H = struct('Info',struct()); + H.SourceRadius = mean(Radii); + H.UnitVectors = XYZ; + H.HRTF_L = HRTF_L; + H.HRTF_R = HRTF_R; + + [~,NameRoot,~] = fileparts(Lib_Name); + H.Info.Name = NameRoot; + H.Info.Source = 'SOFA(Description Goes Here)'; + H.Info.SpatialSize = sprintf('%d discrete directions',NumLoc); + H.Info.FilterSize = sprintf('%d-tap FIR',IRLen); + H.Info.Scaling = 'Normalised to unity gain at 1kHz'; + + obj.process_lib(H, Data.SamplingRate); + obj.Sofa_File = sofa_file; + end + + function updateSOFA(obj, data) + right_shape = permute(data, [3, 2, 1]); + obj.Sofa_File.Data_IR = py.numpy.array(right_shape); + d = size(right_shape); + obj.Sofa_File.SourcePosition = py.numpy.array(zeros(d(1), 3)); + end + + function writeSOFA(obj, file_name) + py.sofar.write_sofa(file_name, obj.Sofa_File); + end + + function IR = XYZ_to_IR( this, XYZ ) + IR = permute(cat(3,this.getHRTF_L(XYZ),this.getHRTF_R(XYZ)),[1,3,2]); + end + + end % methods + + methods (Access=private) + + function process_lib(obj, hrtf_lib, SampleRate) + + MatData = hrtf_lib; + + % Hopefully, this file contains an 'Info' field + if isfield(MatData, 'Info') + obj.Info=MatData.Info; + obj.Lib_Name=obj.Info.Name; + else + obj.Lib_Name=hrtf_lib; + end + + if ~isempty(SampleRate) && isfield(MatData, 'FSample') && SampleRate~=MatData.FSample + error('Sample rate set in Library (%d) does not match requested rate (%d)', ... + MatData.FSample, SampleRate); + end + if isfield(MatData, 'SourceRadius') + obj.SourceRadius = MatData.SourceRadius; + else + fprintf('Warning - no source radius specififed in library\n'); + obj.SourceRadius = 10.0; + end + obj.PlaybackRadius = obj.SourceRadius; + if isfield(MatData, 'FSample') + obj.Lib_SampleRate=MatData.FSample; + elseif ~isempty(SampleRate) + obj.Lib_SampleRate=SampleRate; + else + obj.Lib_SampleRate=48000; % Default, if sample rate not specified anywhere else + end + + obj.Discrete_UnitVectors = MatData.UnitVectors; + if isfield(MatData, 'HRTF_L') + obj.Discrete_HRTFs = cat(3, MatData.HRTF_L, MatData.HRTF_R); + elseif isfield(MatData, 'Band_FR_L') && isfield(MatData, 'GD_Model') + obj.GD_Model = MatData.GD_Model; + obj.buildBand2FR(); + obj.Discrete_HRTFs_GD_Mag = cat(3, ... + [zeros(1,size(MatData.Band_FR_L,2));MatData.Band_FR_L], ... + [zeros(1,size(MatData.Band_FR_L,2));MatData.Band_FR_R]); + else + error('Discrete HRTF library has some data missing'); + end + + obj.buildGDMag(); + obj.Last_UV=[]; + obj.setLen(512); + end + + function setLen(obj,Len) + % SETLEN Set the length of the HRTFs generated by this object + % Example: + % H.setLen(48); + assert(isscalar(Len) && (Len>1) && (Len<4096)); + obj.Len=Len; + obj.Window = sin(min(1,(Len:-1:1)'*4/Len)*pi/2).^2; + obj.Last_UV=[]; + end + + function HRTFs = getHRTF_L(obj,varargin) + % GETHRTF_L Get one or more Left-ear HRTFs from the library + % Given [N] direction-of-arrival vectors, returns an array + % of size [Len]x[N], containing the N left-ear HRTFs. + % Examples: + % HRTFs = H.getHRTF_L(Az,El) % Az and El are both Nx1 or 1xN vectors (in radians) + % HRTFs = H.getHRTF_L( Angles ) % Angles is a 2xN array of Az,El column pairs (in radians) + % HRTFs = H.getHRTF_L(X, Y, Z) % X,Y,Z are all 1xN or Nx1 vectors + % HRTFs = H.getHRTF_L(Vects) % Vects is a 3xN array of X,Y,Z column vectors + UnitVecs = make_unit_vectors(varargin{:}); + if any(size(UnitVecs)~=size(obj.Last_UV)) || ... + any(UnitVecs(:)~=obj.Last_UV(:)) + obj.Interpolate_HRTFs(UnitVecs); + end + HRTFs = obj.Last_HRTFs(:,:,1); + end + + function HRTFs = getHRTF_R(obj,varargin) + % GETHRTF_R Get one or more Right-ear HRTFs from the library + % Given [N] direction-of-arrival vectors, returns an array + % of size [Len]x[N], containing the N right-ear HRTFs. + % Examples: + % HRTFs = H.getHRTF_R(Az,El) % Az and El are both Nx1 or 1xN vectors (in radians) + % HRTFs = H.getHRTF_R( Angles ) % Angles is a 2xN array of Az,El column pairs (in radians) + % HRTFs = H.getHRTF_R(X, Y, Z) % X,Y,Z are all 1xN or Nx1 vectors + % HRTFs = H.getHRTF_R(Vects) % Vects is a 3xN array of X,Y,Z column vectors + UnitVecs = make_unit_vectors(varargin{:}); + if any(size(UnitVecs)~=size(obj.Last_UV)) || ... + any(UnitVecs(:)~=obj.Last_UV(:)) + obj.Interpolate_HRTFs(UnitVecs); + end + HRTFs = obj.Last_HRTFs(:,:,2); + end + + function buildBand2FR(obj) + obj.GD_Mag_FFTLen=2^ceil(log2(obj.Lib_SampleRate)-6.8); + while ~isempty(obj.Discrete_HRTFs) && obj.GD_Mag_FFTLen<2*size(obj.Discrete_HRTFs,1) + obj.GD_Mag_FFTLen = obj.GD_Mag_FFTLen*2; + end + obj.GD_Mag_Band2FR=hrtf_library_loader.MakeBand2FR(obj.GD_Mag_FFTLen,obj.Lib_SampleRate); + + [B,A]=butter(2,0.9); + obj.GD_Mag_LoPass = freqz(B,A,(0:obj.GD_Mag_FFTLen/2)'/obj.GD_Mag_FFTLen*2*pi); + end + + function buildGDMag(obj) + obj.buildBand2FR(); + FR2Band = pinv(obj.GD_Mag_Band2FR); + obj.Discrete_HRTFs_GD_Mag=zeros( ... + size(obj.GD_Mag_Band2FR,2)+1,size(obj.Discrete_HRTFs,2),2); + for k=1:size(obj.Discrete_HRTFs,2) + HRTF_L = obj.Discrete_HRTFs(:,k,1); + HRTF_R = obj.Discrete_HRTFs(:,k,2); + for p=1:2 + % Get mag resp in each band + FR = fft(HRTF_L,obj.GD_Mag_FFTLen); + Mag = FR2Band * (abs(FR(1:end/2+1)).^2); + Mag(1)=abs(FR(1)).^2; % special case for DC + % Now, get the group-delay (excess phase at 1200Hz) + % Let's use mdft for the GD calculation + N = obj.GD_Mag_FFTLen/2; + GD_Bin = 1+floor(1200*2*N/obj.Lib_SampleRate); + FR = m_dft(HRTF_L, N); + Ph = unwrap(angle(FR)); + PhM = -imag(m_hilbert(log(max(1e-4,abs(FR))))); + GD = (PhM(GD_Bin)-Ph(GD_Bin))/pi*N/(GD_Bin-0.5)/obj.Lib_SampleRate; + obj.Discrete_HRTFs_GD_Mag(:,k,p)=[GD;Mag]; + HRTF_L=HRTF_R; + end + obj.Discrete_HRTFs_GD_Mag(1,k,:)=obj.Discrete_HRTFs_GD_Mag(1,k,:)- ... + min(obj.Discrete_HRTFs_GD_Mag(1,k,:)); + end + obj.Last_UV=[]; + end + + function IR=gdMag_to_IR(obj,GDMag) + GD=GDMag(1,:)*obj.Lib_SampleRate + 10; + Mag=GDMag(2:end,:); + if (size(Mag,1)1e-20 + tmpV=tmpV/sqrt(sum(tmpV.^2)); + tmpV=tmpV-v*(tmpV'*v); + IndSubset = IndSubset(tmpV'*obj.Discrete_UnitVectors(:,IndSubset)>0); + end + end + W2=zeros(size(obj.Discrete_UnitVectors,2),1); + W2(GotSet)=lsqnonneg([obj.Discrete_UnitVectors(:,GotSet);0*obj.Discrete_UnitVectors(1,GotSet)+1],[UnitVecs(:,k);1]); + + TempL(:,k)=obj.Discrete_HRTFs_GD_Mag(:,GotSet,1)*W2(GotSet); + TempR(:,k)=obj.Discrete_HRTFs_GD_Mag(:,GotSet,2)*W2(GotSet); + end + end + end % private methods + + methods(Static=true) + + function mat = convert_numpy_2darray(np_array) + py_list = np_array.tolist(); + % List to cell + matCell = cell(py_list)'; + shape = cell(np_array.shape); + shape_mat = zeros(1, length(shape)); + for i=1:length(shape) + shape_mat(i) = double(shape{i}); + end + mat = zeros(shape_mat); + for i = 1:length(matCell) + if isa(matCell{i}, 'py.list') + mat(i,:) = cell2mat(cell(matCell{i})); + end + end + end + + function mat = convert_numpy_3darray(np_array) + py_list = np_array.tolist(); + % List to cell + matCell = cell(py_list)'; + shape = cell(np_array.shape); + shape_mat = zeros(1, length(shape)); + for i=1:length(shape) + shape_mat(i) = double(shape{i}); + end + mat = zeros(shape_mat); + for i = 1:length(matCell) + if isa(matCell{i}, 'py.list') + for j = 1:length(matCell{i}) + if isa(matCell{i}, 'py.list') + mat(i,j,:) = cell2mat(cell(matCell{i}{j})); + end + end + end + end + end + + function [Band2FR,FCs]=MakeBand2FR(FFTLen,FSample) + % Compute the standard band filters for freq domain representations + % + % [Band2FR,FCs]=MakeBand2FR(FFTLen,FSample) + % Builds a column vector, FCs, of size NBx1 (where NB equals the + % number of bands) indicating the center-frequency of the + % bands. These center frequencies are spaced at 200Hz intervals at + % low frequencies and 1/12 octave intervals at higher frequencies. + % + % Band2FR is a real matrix of size (FFTLen/2+1)*NB, where each column + % contains the band-pass filter coefficients, and where the following + % identity is true: + % sum(Band2FR,2)==ones(FFTLen/2+1,1) + % + Band2FR=[]; + FCs=[]; + + TransLen=512; + Tmp=cumsum(sinc((-TransLen/2:TransLen/2)*4/TransLen))'; + Trans=1-Tmp(1:end-1)/Tmp(end-1); + + LastLP=0; + LastBin=1; + while 1 + NewBin = max( LastBin+200/FSample*FFTLen, ... + 1+(LastBin-1)*(2^(1/12))); % 200Hz or 1/12 oct steps + if NewBin>=FFTLen/2, break; end + NewLP=interp1( [-1e10;(0.5:TransLen)'/TransLen*4-1.5;1e10] , ... + [1;Trans;0] , ((1:FFTLen/2+1)'-LastBin)/(NewBin-LastBin) ); + FCs = [FCs ; (LastBin-1)*FSample/FFTLen]; + Band2FR = [Band2FR NewLP(:)-LastLP(:)]; + LastLP=NewLP; + LastBin=NewBin; + end + FCs = [FCs ; (LastBin-1)*FSample/FFTLen]; + Band2FR = [Band2FR 1-LastLP(:)]; + end + end % static methods +end diff --git a/scripts/binauralRenderer_interface/matlab_hrir_generation_scripts/hrtf_support_coefs.mat b/scripts/binauralRenderer_interface/matlab_hrir_generation_scripts/hrtf_support_coefs.mat new file mode 100644 index 0000000000000000000000000000000000000000..a5d6a0ec5fef997a21e039483a557fda2644ce92 --- /dev/null +++ b/scripts/binauralRenderer_interface/matlab_hrir_generation_scripts/hrtf_support_coefs.mat @@ -0,0 +1,3 @@ +version https://git-lfs.github.com/spec/v1 +oid sha256:c8fdb6bc599e4b8b122c4def021d506ddc9c57b82000608cce742a24b3e368a3 +size 395637 diff --git a/scripts/binauralRenderer_interface/matlab_hrir_generation_scripts/m_dft.m b/scripts/binauralRenderer_interface/matlab_hrir_generation_scripts/m_dft.m new file mode 100644 index 0000000000000000000000000000000000000000..b4c04c07edff268e185f67ee6d39736a4e66a4b2 --- /dev/null +++ b/scripts/binauralRenderer_interface/matlab_hrir_generation_scripts/m_dft.m @@ -0,0 +1,76 @@ +%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% +% +% (C) 2022-2023 IVAS codec Public Collaboration with portions copyright Dolby International AB, Ericsson AB, +% Fraunhofer-Gesellschaft zur Foerderung der angewandten Forschung e.V., Huawei Technologies Co. LTD., +% Koninklijke Philips N.V., Nippon Telegraph and Telephone Corporation, Nokia Technologies Oy, Orange, +% Panasonic Holdings Corporation, Qualcomm Technologies, Inc., VoiceAge Corporation, and other +% contributors to this repository. All Rights Reserved. +% +% This software is protected by copyright law and by international treaties. +% The IVAS codec Public Collaboration consisting of Dolby International AB, Ericsson AB, +% Fraunhofer-Gesellschaft zur Foerderung der angewandten Forschung e.V., Huawei Technologies Co. LTD., +% Koninklijke Philips N.V., Nippon Telegraph and Telephone Corporation, Nokia Technologies Oy, Orange, +% Panasonic Holdings Corporation, Qualcomm Technologies, Inc., VoiceAge Corporation, and other +% contributors to this repository retain full ownership rights in their respective contributions in +% the software. This notice grants no license of any kind, including but not limited to patent +% license, nor is any license granted by implication, estoppel or otherwise. +% +% Contributors are required to enter into the IVAS codec Public Collaboration agreement before making +% contributions. +% +% This software is provided "AS IS", without any express or implied warranties. The software is in the +% development stage. It is intended exclusively for experts who have experience with such software and +% solely for the purpose of inspection. All implied warranties of non-infringement, merchantability +% and fitness for a particular purpose are hereby disclaimed and excluded. +% +% Any dispute, controversy or claim arising under or in relation to providing this software shall be +% submitted to and settled by the final, binding jurisdiction of the courts of Munich, Germany in +% accordance with the laws of the Federal Republic of Germany excluding its conflict of law rules and +% the United Nations Convention on Contracts on the International Sales of Goods. +% +%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% +function Y=m_dft(X,Len,Dim) +% MDFT - Modified transform, for use in convolution +% +% Y=M_DFT(X,Len) +% +% Input vector should be Real, and length should be '2*Len'. Output vector +% will be Complex, and length with be 'Len' +% +% Note that M_DFT, and M_IDFT, use a frequency bin re-ording (not normally +% used in DSP implementations) to keep the frequency bins in monotonic +% frequency order. +% +% See also M_IDFT + +% Sort out which dimension we operate along: +if nargin<3 + Dim=[]; +end +[X,shift_perm,nshifts] = shiftdata(X,Dim); +dimensions = size(X); +X=X(:,:); + +% Figure out what length we use using: +if nargin<2 || isempty(Len) + Len=size(X,1)/2; +end +Len=double(ceil(Len)); + +% If the user-selected length is longer than data, zero-pad it +if size(X,1)<(2*Len) + X(2*Len,1)=0; +end + +% Do the MDFT: +Y= fft((X(1:Len,:) - X(Len+1:2*Len,:)*1i) .* ... + repmat(exp(-(0:Len-1)'*1i*pi/Len/2),1,size(X,2)) ,Len,1); +Y(round(size(Y,1)/2+1):end,:) = conj (Y(round(size(Y,1)/2+1):end,:) ); + +Tmp=[1:(Len+1)/2 ; Len:-1:(Len+1)/2]; +Y=Y(Tmp(1:Len),:); + +% restore the original dimensions +dimensions(1) = size(Y,1); +Y = reshape(Y, [size(Y,1), dimensions(2:end)]); +Y = unshiftdata(Y,shift_perm,nshifts); diff --git a/scripts/binauralRenderer_interface/matlab_hrir_generation_scripts/m_hilbert.m b/scripts/binauralRenderer_interface/matlab_hrir_generation_scripts/m_hilbert.m new file mode 100644 index 0000000000000000000000000000000000000000..ba4179dd0036eb9bb9d4949292cf68b1b0977813 --- /dev/null +++ b/scripts/binauralRenderer_interface/matlab_hrir_generation_scripts/m_hilbert.m @@ -0,0 +1,44 @@ +%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% +% +% (C) 2022-2023 IVAS codec Public Collaboration with portions copyright Dolby International AB, Ericsson AB, +% Fraunhofer-Gesellschaft zur Foerderung der angewandten Forschung e.V., Huawei Technologies Co. LTD., +% Koninklijke Philips N.V., Nippon Telegraph and Telephone Corporation, Nokia Technologies Oy, Orange, +% Panasonic Holdings Corporation, Qualcomm Technologies, Inc., VoiceAge Corporation, and other +% contributors to this repository. All Rights Reserved. +% +% This software is protected by copyright law and by international treaties. +% The IVAS codec Public Collaboration consisting of Dolby International AB, Ericsson AB, +% Fraunhofer-Gesellschaft zur Foerderung der angewandten Forschung e.V., Huawei Technologies Co. LTD., +% Koninklijke Philips N.V., Nippon Telegraph and Telephone Corporation, Nokia Technologies Oy, Orange, +% Panasonic Holdings Corporation, Qualcomm Technologies, Inc., VoiceAge Corporation, and other +% contributors to this repository retain full ownership rights in their respective contributions in +% the software. This notice grants no license of any kind, including but not limited to patent +% license, nor is any license granted by implication, estoppel or otherwise. +% +% Contributors are required to enter into the IVAS codec Public Collaboration agreement before making +% contributions. +% +% This software is provided "AS IS", without any express or implied warranties. The software is in the +% development stage. It is intended exclusively for experts who have experience with such software and +% solely for the purpose of inspection. All implied warranties of non-infringement, merchantability +% and fitness for a particular purpose are hereby disclaimed and excluded. +% +% Any dispute, controversy or claim arising under or in relation to providing this software shall be +% submitted to and settled by the final, binding jurisdiction of the courts of Munich, Germany in +% accordance with the laws of the Federal Republic of Germany excluding its conflict of law rules and +% the United Nations Convention on Contracts on the International Sales of Goods. +% +%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% +function x = m_hilbert(xr) +% M_HILBERT - compute the hilbert transform assuming the signal is even and +% periodic + +[xr,shift_dim_n] = shiftdim(xr); +dimensions = size(xr); +xr=xr(:,:); + +tmp = hilbert([xr;conj(flipud(xr))]); +x=tmp(1:end/2,:); + +x = reshape(x, [size(x,1), dimensions(2:end)]); +x = shiftdim(x, -shift_dim_n); diff --git a/scripts/binauralRenderer_interface/matlab_hrir_generation_scripts/m_idft.m b/scripts/binauralRenderer_interface/matlab_hrir_generation_scripts/m_idft.m new file mode 100644 index 0000000000000000000000000000000000000000..85141c44f297a2fad8fd97a9f1f85b763e53dc3b --- /dev/null +++ b/scripts/binauralRenderer_interface/matlab_hrir_generation_scripts/m_idft.m @@ -0,0 +1,72 @@ +%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% +% +% (C) 2022-2023 IVAS codec Public Collaboration with portions copyright Dolby International AB, Ericsson AB, +% Fraunhofer-Gesellschaft zur Foerderung der angewandten Forschung e.V., Huawei Technologies Co. LTD., +% Koninklijke Philips N.V., Nippon Telegraph and Telephone Corporation, Nokia Technologies Oy, Orange, +% Panasonic Holdings Corporation, Qualcomm Technologies, Inc., VoiceAge Corporation, and other +% contributors to this repository. All Rights Reserved. +% +% This software is protected by copyright law and by international treaties. +% The IVAS codec Public Collaboration consisting of Dolby International AB, Ericsson AB, +% Fraunhofer-Gesellschaft zur Foerderung der angewandten Forschung e.V., Huawei Technologies Co. LTD., +% Koninklijke Philips N.V., Nippon Telegraph and Telephone Corporation, Nokia Technologies Oy, Orange, +% Panasonic Holdings Corporation, Qualcomm Technologies, Inc., VoiceAge Corporation, and other +% contributors to this repository retain full ownership rights in their respective contributions in +% the software. This notice grants no license of any kind, including but not limited to patent +% license, nor is any license granted by implication, estoppel or otherwise. +% +% Contributors are required to enter into the IVAS codec Public Collaboration agreement before making +% contributions. +% +% This software is provided "AS IS", without any express or implied warranties. The software is in the +% development stage. It is intended exclusively for experts who have experience with such software and +% solely for the purpose of inspection. All implied warranties of non-infringement, merchantability +% and fitness for a particular purpose are hereby disclaimed and excluded. +% +% Any dispute, controversy or claim arising under or in relation to providing this software shall be +% submitted to and settled by the final, binding jurisdiction of the courts of Munich, Germany in +% accordance with the laws of the Federal Republic of Germany excluding its conflict of law rules and +% the United Nations Convention on Contracts on the International Sales of Goods. +% +%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% +function Y=m_idft(X,Len,Dim) +% M_IDFT - Modified transform, for use in convolution +% +% Y=M_IDFT(X,Len) +% +% Input vector should be Complex, and length should be 'Len'. Output vector +% will be Real, and length with be '2*Len' +% +% Note that M_DFT, and M_IDFT, use a frequency bin re-ording (not normally +% used in DSP implementations) to keep the frequency bins in monotonic +% frequency order. +% +% See also M_DFT + +% Sort out which dimension we operate along: +if nargin<3, + Dim=[]; +end; +[X,shift_perm,nshifts] = shiftdata(X,Dim); +dimensions = size(X); +X=X(:,:); + +% Figure out what length we use using: +if nargin<2 || isempty(Len), + Len=size(X,1); +end; +if size(X,1) 10 ) + fprintf(fid_source,'\n'); + fprintf(fid_source,repmat(' ',1,indent*3)); + end + counter=1; + for C = 1:indices(3) + fprintf(fid_source,'%+ff',real(data(A,B,C))); + if C < indices(3) + if mod(counter,10) == 0 + fprintf(fid_source,','); + else + fprintf(fid_source,', '); + end + end + if mod(counter,10) == 0 && counter ~= indices(3) + fprintf(fid_source,'\n'); + fprintf(fid_source,repmat(' ',1,indent*3)); + end + counter = counter+1; + end + if( indices(3) > 10 ) + fprintf(fid_source,'\n'); + fprintf(fid_source,repmat(' ',1,indent*2)); + end + fprintf(fid_source,'}'); + if B < indices(2) + fprintf(fid_source,','); + end + end + fprintf(fid_source,'\n'); + fprintf(fid_source,repmat(' ',1,indent)); + fprintf(fid_source,'}'); + if A < indices(1) + fprintf(fid_source,','); + end +end +fprintf(fid_source,'\n};\n\n'); + +end % function + diff --git a/scripts/thirdPartyLegalNotices/sofar.txt b/scripts/thirdPartyLegalNotices/sofar.txt new file mode 100644 index 0000000000000000000000000000000000000000..132366686a2aa57ac714f28155504b0cc73830b3 --- /dev/null +++ b/scripts/thirdPartyLegalNotices/sofar.txt @@ -0,0 +1,21 @@ +Applies to sofar + +Copyright (c) [2021] [The pyfar developers] + +Permission is hereby granted, free of charge, to any person obtaining a copy +of this software and associated documentation files (the "Software"), to deal +in the Software without restriction, including without limitation the rights +to use, copy, modify, merge, publish, distribute, sublicense, and/or sell +copies of the Software, and to permit persons to whom the Software is +furnished to do so, subject to the following conditions: + +The above copyright notice and this permission notice shall be included in all +copies or substantial portions of the Software. + +THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR +IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY, +FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE +AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER +LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM, +OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE +SOFTWARE. \ No newline at end of file