Commit 6a49f7f3 authored by Archit Tamarapu's avatar Archit Tamarapu
Browse files

- consolidate FastConv scripts into one entry point to generate ROM tables or binary data 

- modify SD/SHD_2_ROM.m to return a struct with required data
- modify hrtf_library_loader to return normalised IRs for SD
- add script for generating FastConv BRIR tables
parent e2f53203

scripts/binauralRenderer_interface/fastconv/README.txt

0 → 100644
+3 −0
Original line number Diff line number Diff line 
WIP:

The script to generate all FastConv Renderer tables is generate_tables_for_fastconv.m

Options can be set in the script to generate either or both ROM tables and a Binary file.
 
 No newline at end of file

scripts/binauralRenderer_interface/fastconv/cldfb_prototype.mat

0 → 100644LFS
+3 −0
Original line number Diff line number Diff line 
version https://git-lfs.github.com/spec/v1

oid sha256:386bbd3a02c3b95280fff7bd8f5240ee60f8858889cad0f6147e930f5f2aa7e2

size 1246

scripts/binauralRenderer_interface/fastconv/generate_BRIR_CLDFB_FASTCONV.m

0 → 100644
+369 −0
Original line number Diff line number Diff line 
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

%

%   (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 FastConv_SD_BRIR = generate_BRIR_CLDFB_FASTCONV(sofa_file, max_band)

%% generate_BRIR_CLDFB_FASTCONV(rom_c_file, sofa_file): 

% script for getting the binaural room impulse response coefficients in CLDFB domain

% - loads sphere-samples BRIRs given in sofa_file (must match requested speaker positions!)

% - converts SD BRIRs to Complex Low Delay Filter Bank (CLDFB) domain using td2cldfb.m

% - truncates CLDFB BRIRs in frequency and time

% - writes CLDFB BRIRs to c-code ROM tables

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

   sofa_file = fullfile(thispath,'..','BRIRs_sofa','IIS_BRIR_officialMPEG_Combined.sofa');

end

if ~exist('max_band', 'var') || isempty(max_band)

    max_band = 45;

end



%% Load CLDFB protopyte

load('cldfb_prototype.mat');



%% Load HRTFs

addpath("../matlab_hrir_generation_scripts");

sofaData = hrtf_library_loader();

sofaData.readSOFA(char(sofa_file));

ls_struct = get_ls_layout_config('Combined');



IR = permute(sofaData.Data.IR(:,:,:), [3 2 1]);



% match required loudspeaker positions

brirAziRad = sofaData.PosSpherical(1, :).';

brirEleRad = sofaData.PosSpherical(2, :).';

refAziRad = deg2rad(ls_struct.azi)';

refEleRad = deg2rad(ls_struct.ele)';

refVectors = [cos(refAziRad).*cos(refEleRad) sin(refAziRad).*cos(refEleRad) sin(refEleRad)];

sofaVectors = [cos(brirAziRad).*cos(brirEleRad) sin(brirAziRad).*cos(brirEleRad) sin(brirEleRad)];

for ch = 1:length(refAziRad)

    [~, maxIndex] = max(sofaVectors*refVectors(ch,:)');

    chSelect(ch) = maxIndex;

end

IR = IR(:, :, chSelect);



% Resample IRs if needed

fs = sofaData.Lib_SampleRate;

if fs~=48000

    hrir_in = permute(IR,[2 3 1]);

    IR = permute(resamp_hrir(hrir_in,fs),[3 1 2]);

end



%% TD -> CLDFB

clear IR_cldfb;



[BRIR_cldfb, rev_param] = compute_brir_parameters(IR,prototype,max_band);

BRIR_cldfb = permute(BRIR_cldfb, [3 1 4 2]);

BRIR_cldfb = BRIR_cldfb(:,1:max_band,:,:);



FastConv_SD_BRIR.IR = BRIR_cldfb;

FastConv_SD_BRIR.rev_param = rev_param;

end





function [hrir,Nhrir] = resamp_hrir(hrir_in,fs)



Nhrir_in = size(hrir_in,3);

nEar = size(hrir_in,1);

nDir = size(hrir_in,2);



if fs ~= 48000

    

    r  = fs/48000;

    fs = 48000;

    

    Nhrir = Nhrir_in/r;

    

    hrir = zeros([nEar,nDir,Nhrir]);

    data(1:Nhrir_in) = 0.0;

    

    for iEar = 1:nEar

       for iDir = 1:nDir

           data(1:Nhrir_in) = hrir_in(iEar,iDir,:);

           data(1:Nhrir) = decimate(data,r,'fir');

           hrir(iEar,iDir,1:Nhrir) = data(1:Nhrir);

       end

    end

    

else

    

    hrir = hrir_in;

    Nhrir = Nhrir_in;



end



end





function [IR_cldfb,reverb_parameters] = compute_brir_parameters(IR,prototype,max_band)

% Step 1: Compute the propogation time and remove the propogation time

in_struct.propogationTime = 0;

in_struct.fs = 48000;

[IR,in_struct] = compensate_prop_time(IR,in_struct);

in_struct.THR_DE = -20;

in_struct.max_band = max_band;

in_struct.max_index = 273; % Needs to be changed according to IRs



% Step 2: TD -> CLDFB

IR_cldfb = [];

for chIdx = 1:size(IR,3)

    for outChIdx = 1:size(IR,2)

        IR_cldfb(:,:,outChIdx,chIdx) = td2cldfb( IR(:,outChIdx,chIdx), prototype, 60 );

    end

end



% Step 3: Compute EDC (Energy Decay Curve)

in_struct.num_bands = size(IR_cldfb,1);

in_struct.timeSlots = size(IR_cldfb,2);

assert(in_struct.timeSlots > in_struct.max_index)

in_struct.num_channels = size(IR_cldfb,4);

EDC = get_EDC(IR_cldfb,in_struct);



% Step 4: Determine Mixing Time

mixingTime = getMixingTime(EDC, in_struct);



% Compute the time-slot where the IR is truncated

in_struct.NFft = floor((mixingTime));

NFilter = max(in_struct.NFft(1:in_struct.max_band));



% Variable order Filtering

for idx = 1:5

    in_struct.NFft(idx) = NFilter;

end

for idx = 6:10

    in_struct.NFft(idx) = ceil(0.6*NFilter);

end

for idx = 11:20

    in_struct.NFft(idx) = ceil(0.5*NFilter);

end

for idx = 21:30

    in_struct.NFft(idx) = ceil(0.4*NFilter);

end

for idx = 31:in_struct.max_band

    in_struct.NFft(idx) = ceil(0.3*NFilter);

end



% Step 5: Compute reverb parameters

% compute rt60 and energy

reverb_parameters = get_Reverb_parameters(IR_cldfb,in_struct);

reverb_parameters.NFilter = NFilter;



end



function [IR_out,in_struct] = compensate_prop_time(IR_in,in_struct)



% Get the min index

min_idx = zeros(1,size(IR_in,3));

for chIdx = 1:size(IR_in,3)

    [~,idx] = max(abs(IR_in(:,:,chIdx)));

    min_idx(chIdx) = min(idx);

end



prop_time = min(min_idx) - 15;

in_struct.propogationTime = prop_time;

in_struct.latency_s = (min(min_idx) - 1 - prop_time)/in_struct.fs + 31/in_struct.fs;



IR_out = [];

for chIdx = 1:size(IR_in,3)

    IR_out(:,:,chIdx) = IR_in(prop_time:end,:,chIdx);

end



end



function EDC = get_EDC(in,in_struct)

    bands = in_struct.num_bands;

    timeSlots = in_struct.timeSlots;

    numChannels = in_struct.num_channels;

    

    % EDC: Energy decay curve

    EDC = zeros(bands,timeSlots,2,numChannels);

    for chIdx = 1:numChannels

        tmp = in(:,:,:,chIdx);

        for idx = 1:2

            for bandIdx = 1:bands

                sum_1 = sum(abs(tmp(bandIdx,1:end,idx)).^2);

                for timeIdx = 1:timeSlots

                    sum_2 = sum(abs(tmp(bandIdx,timeIdx:end,idx)).^2);

                    EDC(bandIdx,timeIdx,idx,chIdx) = 10*log10(sum_2/sum_1);

                end

            end

        end

    end

end



function mixingTime = getMixingTime(EDC, in_struct)

THR_DE = in_struct.THR_DE;

bands = in_struct.num_bands;

numChannels = in_struct.num_channels;

N = size(EDC,2);

mixingTime = zeros(1,bands);



for idx = 1:bands

    mixingTime(idx) = 0;

    for chIdx = 1:numChannels

        THR = EDC(idx,1,1,chIdx) + THR_DE;

        n = 2;

        while( EDC(idx,n,1,chIdx) > THR )

            n = n + 1;

            if(n > N)

                break;

            end

        end

        mixingTime(idx) = mixingTime(idx) + (n+1);

    end

    for chIdx = 1:numChannels

        THR = EDC(idx,1,2,chIdx) + THR_DE;

        n = 2;

        while( EDC(idx,n,2,chIdx) > THR )

            n = n + 1;

            if(n > N)

                break;

            end

        end

        mixingTime(idx) = mixingTime(idx) + (n+1);

    end

    mixingTime(idx) = mixingTime(idx)/(2*numChannels);

end

end



function RT60 = computeRT60(in,s_idx,e_idx,in_struct)

% define t based on the length of reverb tail

end_idx = length(in);

s_i = (s_idx*in_struct.num_bands)/in_struct.fs;

e_i = (e_idx*in_struct.num_bands)/in_struct.fs;

t = linspace(s_i,e_i,end_idx);



% compute EDC (Energy Decay curve) of the reverb tail

EDC = zeros(1,end_idx);

sum2 = sum(abs(in(1:end_idx)).^2);

for idx = 1:end_idx

    sum1 = sum(abs(in(idx:end_idx)).^2);

    EDC(idx) = 10*log10(sum1/sum2);

end



% compute index, i5 and i35

% i5 is the time required for EDC to reduce by 5dB

% i35 is the time required for EDC to reduce by 35dB

% Based on i5 and i35, we compute the T30 and estimate T60

index = find(EDC==0,1,'last');



i0 = index;

i5 = i0+1;

while EDC(i5) > (EDC(i0)-5)

    i5 = i5+1;

end



i35 = i5+1;

while ((EDC(i35) > (EDC(i0) -35)) && (i35 < end_idx))

    i35 = i35 + 1;

end



% Compute T(30) and T(60)

if (i35 > floor(length(EDC) * 0.95))

    i35 = round(length(EDC) * 0.95);

    amp = EDC(i35);

    fact = -65 / amp;

else

    fact = 2;

end

if (i5) > length(EDC) * 0.5

    i5 = 1;

end



RT30 = t(i35) - t(i5);

RT60 = fact * RT30;

end



function reverbParameters = get_Reverb_parameters(IR_CLDFB,input_struct)



% start the analysis for computing reverb parameters

ENRG = zeros(input_struct.num_bands,1);

RT60 = zeros(input_struct.num_bands,1);

ct_num = 0;

max_index = input_struct.max_index;





disp('Late Reverberation Analysis Started')

for chIdx = 1:input_struct.num_channels

    disp(['Analyzing channel ' num2str(chIdx) '/' num2str(input_struct.num_channels)])

    for outChIdx = 1:2

        for bandIdx = 1:input_struct.num_bands

            % the computation of reverb parameters starts with some overlap

            s_idx = input_struct.NFft(bandIdx)-5;

            

            % Determine the truncation point based on noise floor estimates

            end_idx = get_truncation_point(IR_CLDFB(bandIdx,s_idx:max_index,outChIdx,chIdx));

            e_idx = s_idx + end_idx;

            if(e_idx > max_index)

                e_idx = max_index;

            end



            % Extract the reverb tail that is used for energy/RT60

            % computation

            input_signal = IR_CLDFB(bandIdx,s_idx:e_idx,outChIdx,chIdx);

            

            % Estimate RT60 and Energy

            ENRG_tmp = sum(abs(input_signal.^2));

            RT60_tmp = computeRT60(input_signal,s_idx,e_idx,input_struct);

            

            ENRG(bandIdx) = ENRG(bandIdx) + ENRG_tmp;

            RT60(bandIdx) = RT60(bandIdx) + RT60_tmp;

        end

    end

    ct_num = ct_num + 2;

end



reverbParameters.latency_s = input_struct.latency_s;

reverbParameters.nrgLr = ENRG/ct_num;

reverbParameters.rt60 = RT60/ct_num;

reverbParameters.kAna = input_struct.num_bands;

end



function end_index = get_truncation_point(in)

end_idx = length(in);



% Compute the energy in log domain

log_energy = 10*log10(abs(in(1:end_idx)).^2);

noisefloor = mean(log_energy);



% Decay line(Line Fitting to EDC)

x = 1:end_idx;

order = 1; % order of the polynomial that is used to fit 

% p(1) is the slope, p(2) is the bias

p = polyfit(x,log_energy,order);

decay = p(1)*x + p(2); % Straight line equation



% find the intersection point

cp = find(decay < noisefloor, 1, 'first');

if isempty(cp)

    cp = end_idx;

end



end_index = cp;

end

scripts/binauralRenderer_interface/fastconv/generate_tables_for_fastconv.m

0 → 100644
+71 −0
Original line number Diff line number Diff line 
%

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

%

clear all;

close all;

clc;



%% add path for hrtf_library_loader.m

addpath('../matlab_hrir_generation_scripts/');



%% Set arguments

writeRomFileOutput = true;

writeBinaryOutput = true;

rom_file = fullfile('.', 'ivas_rom_binauralRenderer.c');

bin_file = fullfile('.', 'fastconv_rom.bin');



%% Set input files

hrir_file = fullfile('..', 'HRIRs_sofa', 'HRIR_128_Meth5_IRC_53_Q10_symL_Itrp1_48000.sofa');

brir_file = fullfile('..', 'BRIRs_sofa', 'IIS_BRIR_officialMPEG_Combined.sofa');



%% Generate C-code tables for RENDERER_BINAURAL_FASTCONV (SHD)

disp('Processing HRIRs (FOA) for FastConv renderer...');

FastConv_SHD_IR_FOA = SHD_2_ROM(hrir_file, 1, 128);



disp('Processing HRIRs (HOA2) for FastConv renderer...');

FastConv_SHD_IR_HOA2 = SHD_2_ROM(hrir_file, 2, 128);



disp('Processing HRIRs (HOA3) for FastConv renderer...');

FastConv_SHD_IR_HOA3 = SHD_2_ROM(hrir_file, 3, 128);



%% Generate C-code tables for RENDERER_BINAURAL_FASTCONV (SD)

disp('Processing HRIRs (SD) for FastConv renderer...');

FastConv_SD_IR = SD_2_ROM(hrir_file);



%% Generate C-code tables for RENDERER_BINAURAL_FASTCONV_ROOM (SD)

disp('Processing BRIRs (SD) for FastConv renderer...');

FastConv_SD_BRIR = generate_BRIR_CLDFB_FASTCONV(brir_file);



if writeRomFileOutput

    write_fastconv_rom_table(rom_file, FastConv_SHD_IR_FOA, FastConv_SHD_IR_HOA2, FastConv_SHD_IR_HOA3, FastConv_SD_IR, FastConv_SD_BRIR);

end



if writeBinaryOutput

    write_fastconv_binary_data(bin_file, FastConv_SHD_IR_FOA, FastConv_SHD_IR_HOA2, FastConv_SHD_IR_HOA3, FastConv_SD_IR, FastConv_SD_BRIR);

end

scripts/binauralRenderer_interface/fastconv/td2cldfb.m

0 → 100644
+60 −0
Original line number Diff line number Diff line 
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

%

%   (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 [ out ] = td2cldfb( in, prototype, L )

%   TD2CLDFB 

%   This script takes the td signal along with prototype and length and

%   outputs the frequemcy domain signal



N_in = length(in);

N_prototype = length(prototype);



K_prototype = ceil(N_prototype/L);

K_in = ceil(N_in/L);



N = K_prototype + K_in - 1;



hext = zeros(L*(2*K_prototype + K_in - 2),1);

hext((K_prototype-1)*L+(1:N_in)) = in;   % extend to make all shifts possible



out=zeros(L,N);



Nmid = N_prototype/2;    % midpoint

ls = -Nmid+1:Nmid;

ns = (0:L-1)';

E  = exp(-1i*(pi/L)*(ns+.5)*ls);



for k=1:N

    shift=(k-1)*L;

    out(:,k)=E*(hext(shift+(1:N_prototype)).*prototype(:));

end

end