Merge branch '450-Deliver-Matlab-scripts-for-generation-of-SHD-HRIRs-from-SD-HRIRs' into 'main'

### 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.

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

%

%   (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');

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

%

%   (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)

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

%

%   (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;