Commit 987ffcde authored by kinuthia's avatar kinuthia
Browse files

Matlab script to generate left, right and iac for late reverb tables

- only 48kHz generated so far
- TODO:
  - generate 32kHz and 16kHz
  - write generated tables to file
parent 8bef1c48
Loading
Loading
Loading
Loading
+98 −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 [left_avg_power, right_avg_power, ia_coherence] = compute_lr_energies_and_iac(HRTF_mdfts, in_freq_count, out_freq_count)
    % Compute left/right and coherence of entire HRTF dataset
    %
    % HRTF_mdfts - HRTF dataset in frequency domain.
    % in_freq_count - number of HRTF bins
    % out_freq_count - number of bins in target for energies and coherence
    % values

    hrtf_count = size(HRTF_mdfts,3);

    left_avg_power = zeros(1, out_freq_count);
    right_avg_power = zeros(1, out_freq_count);
    ia_coherence = zeros(1, out_freq_count);

    inp_freq_step = 0.5 / in_freq_count;
    inp_freq_offset = 0.5 * inp_freq_step;
    out_freq_step = 0.5 / out_freq_count;

    norm_freqs = out_freq_step * (1:out_freq_count);
    tbl_indexs = ( norm_freqs - inp_freq_offset ) / inp_freq_step;

    base_idx = 1;
    relative_pos = 0;

    % loop over output frequency bins
    for out_bin_idx=1:out_freq_count
        if (tbl_indexs(out_bin_idx) > 1)
            base_idx = floor(tbl_indexs(out_bin_idx));
            relative_pos = tbl_indexs(out_bin_idx) - base_idx;
            % extrapolation (above last bin), choose nearest
            if(base_idx > in_freq_count-1)
                base_idx = in_freq_count-1;
                relative_pos = 1;
            end
        end

        IA_coherence = zeros(2,1);
        avg_pwr_lr = zeros(2,2);

        for hrtf_idx=1:hrtf_count-1
            lr_pair_0 = HRTF_mdfts(base_idx,:,hrtf_idx);
            lr_pair_1 = HRTF_mdfts(base_idx + 1,:,hrtf_idx);
            avg_pwr_lr(1,:) = avg_pwr_lr(1,:) + real(lr_pair_0).^2 + imag(lr_pair_0).^2;
            avg_pwr_lr(2,:) = avg_pwr_lr(2,:) + real(lr_pair_1).^2 + imag(lr_pair_1).^2;
            IA_coherence(1,1) = IA_coherence(1,1) + real(lr_pair_0(1)) * real(lr_pair_0(2)) + imag(lr_pair_0(1)) * imag(lr_pair_0(2));
            IA_coherence(2,1) = IA_coherence(2,1) + real(lr_pair_1(1)) * real(lr_pair_1(2)) + imag(lr_pair_1(1)) * imag(lr_pair_1(2));
        end

        % compute averages and IA coherence
        avg_pwr_lr = avg_pwr_lr/hrtf_count;
        IA_coherence = IA_coherence/hrtf_count;
        IA_coherence(1,1) = IA_coherence(1,1) / sqrt((avg_pwr_lr(1,1) * avg_pwr_lr(1,2)));
        IA_coherence(2,1) = IA_coherence(2,1) / sqrt((avg_pwr_lr(2,1) * avg_pwr_lr(2,2)));
    
        % Limiting to (0...1) range in case of small numerical errors or negative values
        for i=1:2
            IA_coherence(i,1) = min(IA_coherence(i,1),1);
            IA_coherence(i,1) = max(IA_coherence(i,1),0);
        end

        % Computing weighted average of 2 nearest values (1 below + 1 above) for linear interpolation
        weight_1st = 1 - relative_pos;
        left_avg_power(out_bin_idx) = weight_1st * avg_pwr_lr(1,1) + relative_pos * avg_pwr_lr(1,2);
        right_avg_power(out_bin_idx) = weight_1st * avg_pwr_lr(2,1) + relative_pos * avg_pwr_lr(2,2);
        ia_coherence(out_bin_idx) = weight_1st * IA_coherence(1,1) + relative_pos * IA_coherence(2,1);
    end
end
 No newline at end of file
+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.
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function data_struct = generate_lr_energies_and_iac(sofa_path, sofa_file_name)
    % Generate left and right energies and coherence for late reverb from
    % 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
    %
    % Typical usage:
    %   generate_lr_energies_and_iac('../HRIRs_sofa', 'HRIR_128_Meth5_IRC_53_Q10_symL_Itrp1_48000.sofa')
    %
    
    data_struct = struct.empty(1,0);
    sr           = [48000, 32000, 16000];
    sr_short     = [48, 32, 16];
    sr_dft_size  = [240, 160, 80];

    H = hrtf_library_loader();
    H.readSOFA(char(fullfile(sofa_path, sofa_file_name)));
    
    % 48kHz
    IRs = permute(H.get_Discrete_HRTFs(), [1,3,2]);
    HRTF_mdfts = m_dft(IRs, sr_dft_size(1));
    [left_avg_power, right_avg_power, ia_coherence] = compute_lr_energies_and_iac(HRTF_mdfts, sr_dft_size(1), 257);
    
    % 32 kHz
    % TODO
    
    % 16 kHz
    % TODO
    
    %
    data_struct(1).lr_energies_iac_48kHz = [left_avg_power; right_avg_power; ia_coherence];
    data_struct(1).lr_energies_iac_32kHz = [];
    data_struct(1).lr_energies_iac_16kHz = [];
    data_struct(1).sr = sr;
    data_struct(1).sr_short = sr_short;
    data_struct(1).sr_dft_size = sr_dft_size;
    
end
 No newline at end of file
+4 −0
Original line number Diff line number Diff line
@@ -161,6 +161,10 @@ classdef hrtf_library_loader < handle
            IR = permute(cat(3,this.getHRTF_L(XYZ),this.getHRTF_R(XYZ)),[1,3,2]);
        end

        function discrete_HRTFs = get_Discrete_HRTFs(this)
            discrete_HRTFs = this.Discrete_HRTFs;
        end
    
    end % methods
  
    methods (Access=private)