Initial commit for TD renderer modeling

function Gen_Hrf_IVAS_Binary(dataSpec, modSpec)

%% Config

filePath = modSpec.folderMod; % Path for outputs

modelVersion = 3;

useModel = 0;           % Indicates model format, '0' for Bspline model

useITD = 1;             % use ITD model

always_48k_ITD = 1;     % Use 48 kHz sampling for ITD model also for lower sampling rates, resampling done in real-time in C

resAzim = 1;            % Resolution in degrees for sampled basis functions

resElev = 4;            % Resolution in degreed for sampled basis functions

%% Load model

modelName = [dataSpec.dataBase, '_', ...

        dataSpec.subjId, '_', modSpec.method, '_D_Model.mat'];

fn = fullfile(modSpec.folderMod, modelName);

load(fn,'mod');

mod_hrf_org = mod.mod.interp.hrf;

mod_itd_org = mod.mod.interp.itd;

fs_orig = 48000;

%% Generate paramters for model

for fs = [48000 32000 16000]    

    fs_khz = fs/1000;

    outputBinaryFileName = sprintf('hrfilter_model_v%03d_%dkHz.bin',modelVersion,fs_khz);

    % HR filters, sample basis functions

    mod_hrf = mod_hrf_org;

    mod_itd = mod_itd_org;

    % Resample if needed

    if fs ~= fs_orig

        k_orig = size(mod_hrf.WL{1},2);

        mod_hrf.WL{1} = (fs_orig/fs)*resample(mod_hrf.WL{1}',fs,fs_orig)';

        mod_hrf.WR{1} = (fs_orig/fs)*resample(mod_hrf.WR{1}',fs,fs_orig)';

        k = size(mod_hrf.WL{1},2);

        resample_factor = fs/fs_orig;

        if useITD && ~always_48k_ITD

            % Rescale such that the number of samples corresponds to the correct delay.

            mod_itd.W = mod_itd.W * resample_factor;

        end

    end

    % Azimuth

    M = length(mod_hrf.azimKSeq);

    N = sum(unique([mod_hrf.azimBfNum{:}])~=1);

    unique_azimKSeq = cell(N,1);

    num_unique_splines = 0;

    azimShapeIdx = -1*ones(M,1);

    azimShapeSampFactor = -1*ones(M,1);

    a_num_points = -1*ones(M,1);

    a_range = cell(N,1);

    azimBfVec = cell(N,1);

    % Find basis function shapes

    [azimBfNum_sorted, azimBfNum_sorted_idx] = sort([mod_hrf.azimBfNum{:}]);

    last_L = -1; % previous number of basis functions

    for m = 1:M

        L = azimBfNum_sorted(m);

        sm = azimBfNum_sorted_idx(m);

        if(L == 1)

            azimShapeIdx(sm) = -1;

            azimShapeSampFactor(sm) = -1;

            last_L = L;

        elseif( L == last_L )

            azimShapeIdx(sm) = num_unique_splines-1; % -1 for C code indexing

            azimShapeSampFactor(sm) = 1;

        elseif (last_L > 1 && rem(L,last_L)==0)

            azimShapeIdx(sm) = num_unique_splines-1; % -1 for C code indexing

            azimShapeSampFactor(sm) = L/last_L;

        else

            num_unique_splines = num_unique_splines+1;

            a_knot_interval = (mod_hrf.azimKSeq{sm}(end)-mod_hrf.azimKSeq{sm}(1))/mod_hrf.azimBfNum{sm}; % deg between knot points

            a_num_points(num_unique_splines) = ceil(a_knot_interval / resAzim);

            a_step = a_knot_interval/a_num_points(num_unique_splines); % sampling in degrees

            a_range{num_unique_splines} = mod_hrf.azimKSeq{sm}(1):a_step:mod_hrf.azimKSeq{sm}(end);

            i = 1;

            for a = a_range{num_unique_splines}

                azimBfVec{num_unique_splines}(i,:) = BSplineSampVec( mod_hrf.azimBf{sm}, mod_hrf.azimKmSeq{sm}, a );

                i = i+1;

            end

            unique_azimKSeq{num_unique_splines} = mod_hrf.azimKSeq{sm};

            azimShapeIdx(sm) = num_unique_splines-1; % -1 for C code indexing

            azimShapeSampFactor(sm) = 1;

            last_L = L;

        end

    end

    % Elevation

    e_knot_interval = (mod_hrf.elevKSeq{1}(end)-mod_hrf.elevKSeq{1}(1))/(mod_hrf.elevBfNum{1}-3); % deg between knot points

    e_num_points = ceil(e_knot_interval / resElev);

    e_step = e_knot_interval/e_num_points; % sampling in degrees

    e_range = mod_hrf.elevKSeq{1}(1):e_step:mod_hrf.elevKSeq{1}(end);

    j=1;

    elevBfVec = -1*ones(length(e_range),size(mod_hrf.elevBf{1},3));

    for e = e_range

        elevBfVec(j,:) = BSplineSampVec( mod_hrf.elevBf{1}, mod_hrf.elevKmSeq{1}, e );

        j = j+1;

    end

    % HR azimuth shapes

    azimSplineShape = cell(num_unique_splines,1);

    len_a_shapes = cell(num_unique_splines,1);

    for n = 1:num_unique_splines

        b = azimBfVec{n}(:,2); 

        azimSplineShape{n} = b(1:2*a_num_points(n)+1);

        len_a_shapes{n} = 2*a_num_points(n)+1;

    end

    % HR elevation shapes

    elevSplineShape = [];

    b1 = elevBfVec(:,1); 

    elevSplineShape{1} = b1(1:(e_num_points+1)); 

    b2 = elevBfVec(:,2);

    elevSplineShape{2} = b2(1:(e_num_points*2+1)); 

    b3 = elevBfVec(:,3);

    elevSplineShape{3} = b3(1:(e_num_points*3+1)); 

    b4 = elevBfVec(:,4); 

    elevSplineShape{4} = b4((1:e_num_points*4+1)); 

    elevSplineShape{4} = elevSplineShape{4}(1:(end+1)/2); % symmetric

    len_e = [length(elevSplineShape{1}) length(elevSplineShape{2}) length(elevSplineShape{3}) length(elevSplineShape{4})];

    start_e = [0, len_e(1), sum(len_e(1:2)), sum(len_e(1:3))];

    elevSplineShape_all = [elevSplineShape{1};elevSplineShape{2};elevSplineShape{3};elevSplineShape{4}];

    % ITD, generate basis functions

    % Azimuth

    resAzim = 1; % resolution, degrees

    a_knot_interval_ITD = (mod_itd.azimKSeq{2}(end)-mod_itd.azimKSeq{2}(1))/((mod_itd.azimBfNum{2}+1)/2-3); % deg between knot points

    a_num_points_ITD = ceil(a_knot_interval_ITD / resAzim);

    a_step_ITD = a_knot_interval_ITD/a_num_points_ITD; % sampling in degrees

    a_range_ITD = mod_itd.azimKSeq{2}(1):a_step_ITD:mod_itd.azimKSeq{2}(end);

    azimBfVecITD = -1*ones(length(a_range_ITD),size(mod_itd.azimBf{2},3));

    i = 1;

    for a = a_range_ITD

        azimBfVecITD(i,:) = BSplineSampVecITD( mod_itd.azimBf{2}, mod_itd.azimKmSeq{2}, a );

        i = i+1;

    end

    % Elevation

    e_knot_interval_ITD = (mod_itd.elevKSeq(end)-mod_itd.elevKSeq(1))/(mod_itd.elevBfNum-3); % deg between knot points

    e_num_points_ITD = ceil(e_knot_interval_ITD / resElev);

    e_step_ITD = e_knot_interval_ITD/e_num_points_ITD; % sampling in degrees

    e_range_ITD = mod_itd.elevKSeq(1):e_step_ITD:mod_itd.elevKSeq(end);

    elevBfVecITD = -1*ones(length(e_range_ITD),size(mod_itd.elevBf,3));

    j=1;

    for e = e_range_ITD

        elevBfVecITD(j,:) = BSplineSampVecITD( mod_itd.elevBf, mod_itd.elevKmSeq, e );

        j = j+1;

    end

    % ITD azimuth shapes

    azimSplineShapeITD = [];

    b1 = azimBfVecITD(:,1); 

    azimSplineShapeITD{1} = b1(1:(a_num_points_ITD+1)); 

    b2 = azimBfVecITD(:,2); 

    azimSplineShapeITD{2} = b2(1:(a_num_points_ITD*2+1));

    b3 = azimBfVecITD(:,3); 

    azimSplineShapeITD{3} = b3(1:(a_num_points_ITD*3+1)); 

    b4 = azimBfVecITD(:,4); 

    azimSplineShapeITD{4} = b4((1:a_num_points_ITD*4+1)); 

    azimSplineShapeITD{4} = azimSplineShapeITD{4}(1:(end+1)/2); % symmetric

    len_a_ITD = [length(azimSplineShapeITD{1}) length(azimSplineShapeITD{2}) length(azimSplineShapeITD{3}) length(azimSplineShapeITD{4})];

    start_a_ITD = [0, len_a_ITD(1), sum(len_a_ITD(1:2)), sum(len_a_ITD(1:3))];

    azimSplineShapeITD_all = [azimSplineShapeITD{1};azimSplineShapeITD{2};azimSplineShapeITD{3};azimSplineShapeITD{4}];

    % ITD elevation shapes

    elevSplineShapeITD = [];

    b1 = elevBfVecITD(:,1); 

    elevSplineShapeITD{1} = b1(1:(e_num_points_ITD+1)); 

    b2 = elevBfVecITD(:,2);

    elevSplineShapeITD{2} = b2(1:(e_num_points_ITD*2+1)); 

    b3 = elevBfVecITD(:,3); 

    elevSplineShapeITD{3} = b3(1:(e_num_points_ITD*3+1));

    b4 = elevBfVecITD(:,4); 

    elevSplineShapeITD{4} = b4((1:e_num_points_ITD*4+1));

    elevSplineShapeITD{4} = elevSplineShapeITD{4}(1:(end+1)/2); % symmetric

    len_e_ITD = [length(elevSplineShapeITD{1}) length(elevSplineShapeITD{2}) length(elevSplineShapeITD{3}) length(elevSplineShapeITD{4})];

    start_e_ITD = [0, len_e_ITD(1), sum(len_e_ITD(1:2)), sum(len_e_ITD(1:3))];

    elevSplineShapeITD_all = [elevSplineShapeITD{1};elevSplineShapeITD{2};elevSplineShapeITD{3};elevSplineShapeITD{4}];

    %% Write to binary file

    % Open file

    fclose('all');

    if(exist(filePath, 'dir')==0)

      mkdir(filePath);

    end

    fileID = fopen(fullfile(filePath,outputBinaryFileName),'w');

    % Header

    % Format for file:

    % 0 = BSpline model

    fwrite(fileID, useModel, 'short');

    % ITD model active/inactive:

    % 1 = ITD model is used

    % 0 = ITD model is not used

    fwrite(fileID, useITD, 'short');

    % The sampling frequency in kHz of the HR filter set:

    fwrite(fileID, fs_khz, 'short');

    % General - model-specific parts

    fwrite(fileID, size(mod_hrf.elevBf{1}, 1), 'short');            % N = 4 i.e. coefficients for cubic including zeroth order

    fwrite(fileID, size(mod_hrf.WR{1}, 2), 'short');                % K, filter length

    % Elevation model structure

    elevDim2 = size(mod_hrf.elevBf{1}, 2);

    elevDim3 = size(mod_hrf.elevBf{1}, 3);    

    fwrite(fileID, elevDim2, 'short');                              % elevDim2

    fwrite(fileID, elevDim3, 'short');                              % elevDim3 = P

    fwrite(fileID, mod_hrf.elevKSeq{1}, 'float');                   % length = elevDim3-2

    % Azimuth model structure

    azim_start_idx = 0;

    for i = 1:elevDim3

        azimDim2 = size(mod_hrf.azimBf{i}, 2);

        azimDim3 = size(mod_hrf.azimBf{i}, 3);

        fwrite(fileID, azimDim2, 'short');                          % azimDim2

        fwrite(fileID, azimDim3, 'short');                          % azimDim3 = Q

        fwrite(fileID, azim_start_idx, 'short');                    % start azim index per elevation

        azim_start_idx = azim_start_idx + azimDim3;

        fwrite(fileID, mod_hrf.azimKSeq{i},       'float');         % length = azimDim3+1

    end

    % Weights

    fwrite(fileID, size(mod_hrf.WL{1},1), 'short');                 % (P*Q)

    fwrite(fileID, mod_hrf.WL{1}, 'float');                         % (P*Q) by K

    fwrite(fileID, mod_hrf.WR{1}, 'float');                         % (P*Q) by K

    % Azimuth basis functions

    fwrite(fileID, num_unique_splines, 'short');                    % number of unique spline functions

    for i = 1:num_unique_splines 

        fwrite(fileID, len_a_shapes{i}, 'short');                   % length of azimuth shape

        fwrite(fileID, azimSplineShape{i}, 'float');                % azimuth shape

        fwrite(fileID, a_num_points(i), 'short');                   % samples between knot points

    end

    fwrite(fileID, azimShapeIdx, 'short');                          % indices for spline functions to use

    fwrite(fileID, azimShapeSampFactor, 'short');                   % decimation factor for spline functions

    % Elevation basis functions

    fwrite(fileID, len_e, 'short');                                 % length of elevation shapes

    fwrite(fileID, start_e, 'short');                               % start idx (C indexing) of elevation shapes

    fwrite(fileID, length(elevSplineShape_all), 'short');           % total length elevation shapes

    fwrite(fileID, elevSplineShape_all, 'float');                   % elevation shapes

    fwrite(fileID, e_num_points, 'short');                          % samples between knot points

    % If ITD model is used, parameters are stored in 2nd part of the same file

    if useITD

        % General

        fwrite(fileID, size(mod_itd.elevBf, 1), 'short');           % N = 4 i.e. coefficients for cubic including zeroth order

        %fwrite(fileID, size(mod_itd.W, 2), 'short');               % K = 1 always for ITD, so do not need to write.

        % Elevation model structure

        elevDim2 = size(mod_itd.elevBf, 2);

        elevDim3 = size(mod_itd.elevBf, 3);        

        fwrite(fileID, elevDim2, 'short');                          % elevDim2

        fwrite(fileID, elevDim3, 'short');                          % elevDim3 = P

        fwrite(fileID, mod_itd.elevKSeq,       'float');            % length = elevDim3-2

        % Azimuth model structure

        azimDim2 = size(mod_itd.azimBf{2}, 2);

        azimDim3 = size(mod_itd.azimBf{2}, 3);        

        fwrite(fileID, azimDim2, 'short');                          % azimDim2

        fwrite(fileID, azimDim3, 'short');                          % azimDim3 = Q

        fwrite(fileID, mod_itd.azimKSeq{2}, 'float');               % length = azimDim3+1

        % Weights

        fwrite(fileID, size(mod_itd.W,1), 'short');                 % (P*Q)

        fwrite(fileID, mod_itd.W, 'float');                         % (P*Q) by K

        % Azimuth basis functions

        fwrite(fileID, len_a_ITD, 'short');                         % length of azimuth shapes

        fwrite(fileID, start_a_ITD, 'short');                       % start idx (C indexing) of azimuth shapes

        fwrite(fileID, length(azimSplineShapeITD_all), 'short');    % total length azimuth shapes

        fwrite(fileID, azimSplineShapeITD_all, 'float');            % azimuth shape

        fwrite(fileID, a_num_points_ITD, 'short');                  % samples between knot points

        % Elevation basis functions

        fwrite(fileID, len_e_ITD, 'short');                         % length of elevation shapes

        fwrite(fileID, start_e_ITD, 'short');                       % start idx (C indexing) of elevation shapes

        fwrite(fileID, length(elevSplineShapeITD_all), 'short');    % total length elevation shapes

        fwrite(fileID, elevSplineShapeITD_all, 'float');            % elevation shapes

        fwrite(fileID, e_num_points_ITD, 'short');                  % samples between knot points

    end

    % Close file

    fclose(fileID);

    fprintf("Wrote model parameters to: %s\n",fullfile(filePath,outputBinaryFileName));

end % fs loop

end % function

 No newline at end of file

function [dataSpec, modSpec] = HrfModBsp_Config(dataSpec)

    %% HRF directory

    myPath = dataSpec.hrfDir;

    %% settings for modeling individual dataset

    % modeling method 

    modSpec.method='BspTdFir'; 

    % number of DFT points

    modSpec.nfft = 512;

    % enable/disable regularization

    modSpec.flag.reg = 1; 

    % B-spline polynomial order

    modSpec.hrf.pAng = 3;

    modSpec.itd.PAng = 3;

    % for delay correction

    modSpec.upSampFac = 1;  % Upsampling factor for delay correction

    modSpec.maxDel = 5;     % Number of samples shifted (upsampled by modSpec.upSampFac)

    modSpec.numIter = 20;

    modSpec.trainSetSpec = 'TRAINSET_SPAT-AREA';

    switch dataSpec.dataBase

        case 'Orange'

            % separate ITD model

            dataSpec.truncType = 'D';

            dataSpec.onsetDelayInit = 1;

            % Knot sequence specification

            modSpec.hrf.elevKSeq{1, 1} = -90:15:90;

            azimKSeq = 0:10:360;

            N = length(modSpec.hrf.elevKSeq{1, 1});

            modSpec.hrf.azimKSeq{1,1} = [0, 360];

            for n = 2:N+1

                modSpec.hrf.azimKSeq{1,n} = azimKSeq;

            end

            modSpec.hrf.azimKSeq{1,N+2} = [0, 360];

            modSpec.itd.elevKSeq = -90:12:90;

            modSpec.itd.azimKSeq = 0:10:180;            

            % desired sample rate

            modSpec.fs = 48000;   

            % regularization factor

            modSpec.hrf.lambda = 0.5;

            modSpec.itd.lambda = 0.2; 

            modSpec.enThr  = 0.3;

            % folder to save the model

            myPath = fullfile(myPath, 'Orange_53');

            dataSpec.folderNameOnset = ...

                fullfile(myPath, 'OnsetDelay');

            modSpec.folderMod = myPath;            

    end

%% check if the required folder exsits, if not, create one

% the folder for model

if ~isfolder(modSpec.folderMod)

    mkdir(modSpec.folderMod)

end

 No newline at end of file