output converted filters, improve comments

%    the United Nations Convention on Contracts on the International Sales of Goods.

% 

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

function SHD_2_ROM( sofa_file, rom_c_file )

% SHD_2_ROM( sofa_file, rom_c_file )

function IR_cldfb = SHD_2_ROM( rom_c_file, sofa_file, py_path )

% SHD_2_ROM( rom_c_file, sofa_file, py_path )

%

% converts Sphere-sampled Head Related Impulse Responses (HRIRs) to the Spherical Harmonics domain (SHD)  

% to the Complex Low Delay Filter Bank (CLDFB) domain to c-code ROM tables.

%   

%   HRIRs --> SHD --> CLDFB --> ROM

%   HRIRs --> SHD: Using ./generate_HOA_HRIRs_MOD_lens.m

%   SHD --> CLDFB: Least quares optimization

%   SHD --> CLDFB: Least squares optimization

%   CLDFB --> ROM: Using ../param_bin/writeData3L.m

[thispath,~,~] = fileparts(mfilename('fullpath'));

thispath = [thispath,filesep];

%% Load SHD HRIRs

if exist('sofa_file','var') && ~isempty(sofa_file)

if 0 % exist('sofa_file','var') && exist('py_path','var')

   % Note: this path not yet functional!

   % requires

   % python -m pip install sofar

   % python -m pip install numpy

   %py_path = 'C:\Users\xxxx\AppData\Local\Programs\Python\Python39\python.exe';

   %sofa_file = fullfile(thispath,'..','HRIRs_sofa',filesep,'HRIR_128_48000_dolby_SBA3.sofa'); % 'HRIR_128_Meth5_IRC_53_Q10_symL_Itrp1_48000.sofa'

   [sofa_path,sofa_name] = fileparts(sofa_file);

   py_path = 'C:\Users\hmund\AppData\Local\Programs\Python\Python39\python.exe';

   %sofa_path = fullfile(thispath,'..','HRIRs_sofa',filesep);

   %sofa_name = 'HRIR_128_Meth5_IRC_53_Q10_symL_Itrp1_48000.sofa';%'HRIR_128_48000_dolby_SBA3.sofa';

   % convert HRIRs to SHD

   IR = generate_HOA_HRIRs_MOD_lens(1, py_path, sofa_path, sofa_name, 128);

   % convert sphere-sampled HRIRs to SHD HRIRs

   hrir_len = 128;

   IR = generate_HOA_HRIRs_MOD_lens(1, py_path, sofa_path, sofa_name, hrir_len);

else

   % load pre-converted SHD 

   % load pre-converted SHD HRIRs

   in_file = [thispath,'HRIR_128_48000_dolby_SBA3.mat'];

   disp(['Loading Spherical Harmonics Domain HRIRs from ',in_file])

   load(in_file,'IR');

end

[~,num_ears,num_ch] = size(IR);

%% SHD -> CLDFB

%% SHD -> CLDFB via least squares error optimization

num_cldfb_taps = 3;

IR_cldfb = zeros(60,num_cldfb_taps,num_ears,num_ch);

eval_flag = 0; % optional, requires signal processing toolbox

IR_cldfb = zeros(60,num_cldfb_taps,num_ears,num_ch); % 60 frequency bands

eval_flag = 0; % optional, = 1 requires signal processing toolbox (fftfilt)

for pos = 1:num_ch

   disp(['Channel ',num2str(pos),'/',num2str(num_ch)])

   for ear = 1:num_ears

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 = permute(IR_cldfb, [3 1 4 2]);

order_int = (round(sqrt(size(IR_cldfb,3))-1));

IR_cldfb_rom = permute(IR_cldfb, [3 1 4 2]); % after permute: [ears(2), bands(60), chans(16), taps(3)]

order_int = (round(sqrt(size(IR_cldfb_rom,3))-1));

if order_int == 1

    order = 'FOA';

else

fprintf(fid, ['const float FASTCONV_' order '_latency_s = %10.9ff;\n'], latency_s);

addpath(fullfile(thispath,'..','param_bin')); % for writeData3L()

writeData3L(fid, ['const float leftHRIRReal_' order '[BINAURAL_CONVBANDS][' order '_CHANNELS][BINAURAL_NTAPS_SBA]'], real( squeeze(IR_cldfb(1,1:max_band,:,:))));

writeData3L(fid, ['const float leftHRIRImag_' order '[BINAURAL_CONVBANDS][' order '_CHANNELS][BINAURAL_NTAPS_SBA]'], imag( squeeze(IR_cldfb(1,1:max_band,:,:))));

writeData3L(fid, ['const float rightHRIRReal_' order '[BINAURAL_CONVBANDS][' order '_CHANNELS][BINAURAL_NTAPS_SBA]'], real( squeeze(IR_cldfb(2,1:max_band,:,:))));

writeData3L(fid, ['const float rightHRIRImag_' order '[BINAURAL_CONVBANDS][' order '_CHANNELS][BINAURAL_NTAPS_SBA]'], imag( squeeze(IR_cldfb(2,1:max_band,:,:))));

writeData3L(fid, ['const float leftHRIRReal_' order '[BINAURAL_CONVBANDS][' order '_CHANNELS][BINAURAL_NTAPS_SBA]'], real( squeeze(IR_cldfb_rom(1,1:max_band,:,:))));

writeData3L(fid, ['const float leftHRIRImag_' order '[BINAURAL_CONVBANDS][' order '_CHANNELS][BINAURAL_NTAPS_SBA]'], imag( squeeze(IR_cldfb_rom(1,1:max_band,:,:))));

writeData3L(fid, ['const float rightHRIRReal_' order '[BINAURAL_CONVBANDS][' order '_CHANNELS][BINAURAL_NTAPS_SBA]'], real( squeeze(IR_cldfb_rom(2,1:max_band,:,:))));

writeData3L(fid, ['const float rightHRIRImag_' order '[BINAURAL_CONVBANDS][' order '_CHANNELS][BINAURAL_NTAPS_SBA]'], imag( squeeze(IR_cldfb_rom(2,1:max_band,:,:))));

fclose(fid);

rmpath(fullfile(thispath,'..','param_bin'));

   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; % some good margin for shifts while also integer number of strides

   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;

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')

legend('time domain','CLDFB domain','difference')

title('Filtered noise')

function [h, D, S, L] = get_cldfb_filter()