The Auditory Modeling Toolbox

Applies to version: 0.10.0

View the help

Go to function

ITDESTIMATOR - Estimate ITD from a binaural signal

Program code:

function [toa_diff,toa,IACC] = itdestimator(Obj,varargin)
%ITDESTIMATOR Estimate ITD from a binaural signal
%   Usage: itd = itdestimator(data,mode,threshlvl,lowpass,butterpoly,upper_cutfreq) 
%
%   Input parameters:
% 
%       data:       SOFA object or IR matrix with dimensions: 
%                     emitter x receiver x time
% 
%       fs:         sampling rate, used only if data provided as matrix
% 
%       mode:       (optional) Select one estimation methods
%                   ('Threshold' (default),'Cen_e2','MaxIACCr', 'MaxIACCe',
%                   'CenIACCr','CenIACCe', 'CenIACC2e', 'PhminXcor','IRGD')
%             
%                   For details concerning estimation methods see:
%                   'http://asa.scitation.org/doi/10.1121/1.4996457'
% 
%       lowpass:    (optional) Bandwidth considered. 'lp' for lowpass (default), 'bb' for broadband
%
%       peak:       (optional) Method to find the max, used in Threshold mode only. 
%                   'hp' for max (default), 'fb' for findpeak
% 
%       threshlvl:  (optional) Set threshold level for 'Threshold' mode in dB.        
%                   Default is -10 dB. 
%
%       butterpoly: (optional) Select the order of the polynom
%                   applied in the butterworth filter. ( 2 =< i =< 10 )
%                   Default is 10. 
% 
%       upper_cutfreq: (optional) Set frequency of lowpass cutoff in Hz.
%                      Default is 3000 Hz. 
% 
%       lower_cutfreq: (optional) Set frequency of highpass cutoff in Hz, 
%                      only used in IRGD mode. Default is 1000 Hz.    
% 
% 
%   Output parameters:
% 
%       itd:        interaural time difference in seconds
%       toa:        detected activation onsets for left and right channels
%       IACC:       interaural cross-correlation coefficient
%                   Available on when xcorr is used (modes: 'MaxIACCr', 'MaxIACCe',
%                   'CenIACCr','CenIACCe', 'CenIACC2e')
% 
% 
%   Purpose:
%   Estimates the ITD based on biaural impulse responses.
%   Several different estimaton methods can be chosen.
%   MaxIAACe is recommended.
% 
%   Requirements: 
%   -------------
%
%   1) SOFA API from http://sourceforge.net/projects/sofacoustics for Matlab (in e.g. thirdparty/SOFA)
% 
%
%   Examples:
%   ---------
% 
%   Obj = SOFAload(fullfile(SOFAdbPath,'baumgartner2017','hrtf b_nh15.sofa'));
%   toa_diff = itdestimator(Obj,'MaxIACCe','lp','upper_cutfreq',3000)
%   
%   With these settings the estimator uses the MaxIAAce method and applies
%   a lowpass with a cut off frequency of 3 kHz.
%   
%   The output array is structured as the SOFA Data.IR
%   If you would like to select for example only data on the horizontal
%   plane you could:
%
%   plane_idx = find( Obj.SourcePosition(:,2) == 0 );
%   plane_angle = Obj.SourcePosition(plane_idx,1);
%
%   Url: http://amtoolbox.sourceforge.net/data/amt-test/htdocs/amt-0.10.0/doc/general/itdestimator.php

% Copyright (C) 2009-2020 Piotr Majdak and the AMT team.
% This file is part of Auditory Modeling Toolbox (AMT) version 0.10.0
%
% This program is free software: you can redistribute it and/or modify
% it under the terms of the GNU General Public License as published by
% the Free Software Foundation, either version 3 of the License, or
% (at your option) any later version.
%
% This program is distributed in the hope that it will be useful,
% but WITHOUT ANY WARRANTY; without even the implied warranty of
% MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
% GNU General Public License for more details.
%
% You should have received a copy of the GNU General Public License
% along with this program.  If not, see <http://www.gnu.org/licenses/>.

%
% Copyright (C) 2009-2015 Piotr Majdak and the AMT team.
% This file is part of Auditory Modeling Toolbox (AMT) version 0.9.9
%
% This program is free software: you can redistribute it and/or modify
% it under the terms of the GNU General Public License as published by
% the Free Software Foundation, either version 3 of the License, or
% (at your option) any later version.
%
% This program is distributed in the hope that it will be useful,
% but WITHOUT ANY WARRANTY; without even the implied warranty of
% MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
% GNU General Public License for more details.
%
% You should have received a copy of the GNU General Public License
% along with this program.  If not, see <http://www.gnu.org/licenses/>.

% AUTHOR: Laurin Steidle
% data matrix option: Robert Baumgartner


    % ---------------------- ltfatarghelper -------------------------------
    definput.import = {'itdestimator'};
    [flags,kv]=ltfatarghelper({},definput,varargin);

    
    % ---------------------- renaming input parameter ---------------------
    
    if isstruct(Obj)
      pos = Obj.API.M;
      ear = Obj.API.R;
      Ns  = Obj.API.N;
      IR = Obj.Data.IR;
      fs  = Obj.Data.SamplingRate;
    else
      pos = size(Obj,1);
      ear = size(Obj,2);
      Ns  = size(Obj,3);
      IR = Obj;
      if isempty(kv.fs)
        error('RB: No sampling rate (fs) provided.')
      end
      fs  = kv.fs;
    end

    % ---------------------- initialising variables -----------------------

    toa = zeros(pos,ear);
    toa_diff = zeros(pos,1);
    IACC = zeros(pos,1);
    
    amt_disp('itdestimator \n',flags.disp)

    % ---------------------- Applying low-pass ----------------------------

    if flags.do_lp
        amt_disp('Applying Butterworth low pass \n',flags.disp)
        amt_disp(strcat('Polynomial order of Butterworth filter: '...
            ,num2str(kv.butterpoly),'\n'),flags.disp)
        amt_disp(strcat('Cut-off frequency is: '...
            ,num2str(kv.upper_cutfreq),' Hz\n\n'),flags.disp)
        cut_off_freq_norm = kv.upper_cutfreq/(fs/2);
        [lp_a,lp_b] = butter(kv.butterpoly,cut_off_freq_norm);
        f_ir = zeros(pos,ear,Ns);
        for ii=1:pos
            for jj=1:ear  
                sir = squeeze( IR(ii,jj,:) );
                f_sir = filter(lp_a,lp_b,sir);
                f_ir(ii,jj,:) = f_sir;
            end
        end
    
    else
        amt_disp('No low pass filter is applied \n\n',flags.disp)
        f_ir = IR;
    end

    % ---------------------- estimating itd -------------------------------
    % ---------------------------------------------------------------------

    % ---------------------- Threshold ------------------------------------
    switch(flags.mode)
        case 'Threshold'
            amt_disp('Threshold mode \n',flags.disp)
            amt_disp(strcat('Threshold level is: ',...
                num2str(kv.threshlvl),'dB \n\n'),flags.disp)
            
            if flags.do_fp
                for ii=1:pos
                     for jj=1:ear
                        indB = 0.5*mag2db(squeeze(f_ir(ii,jj,:)).^2);
                        [~,B] = findpeaks(indB);
                        th_value = indB(B(1)) + kv.threshlvl;
                        toa(ii,jj) = find(indB>th_value,1);
                    end
                    toa_diff(ii) = toa(ii,1) - toa(ii,2);     
                end

            else
                for ii=1:pos
                    for jj=1:ear
                        indB = 0.5*mag2db(squeeze(f_ir(ii,jj,:)).^2);
                        th_value = max(indB) + kv.threshlvl;
                        toa(ii,jj) = find(indB>th_value,1);
                    end
                    toa_diff(ii) = toa(ii,1) - toa(ii,2);     
                end
            end
                        

    % ---------------------- Cross-Correlation ----------------------------        
        case 'Cen_e2'
            amt_disp('Cen-e2 mode \n',flags.disp)
            for ii=1:pos
                for jj = 1:ear
                    e_sir_sq = envelope(squeeze(f_ir(ii,jj,:))).^2;
                    toa(ii,jj) = centroid(transpose(1:Ns),e_sir_sq);
                end
                toa_diff(ii) = toa(ii,1) - toa(ii,2);     
            end        
            

        case 'MaxIACCr'
            amt_disp('MaxIACCr mode \n',flags.disp)
            for ii=1:pos                
                cc = xcorr(squeeze(f_ir(ii,1,:)),squeeze(f_ir(ii,2,:)));
                [IACC(ii),idx_lag] = max(abs(cc));
                toa_diff(ii) = idx_lag - Ns;                
            end
            if flags.do_guesstoa
                toa = guesstoa(toa_diff,toa, kv.avgtoa);
            end
            
        case 'MaxIACCe'
            amt_disp('MaxIACCe mode \n',flags.disp)
            for ii=1:pos
                e_sir1 = envelope(squeeze(f_ir(ii,1,:)));
                e_sir2 = envelope(squeeze(f_ir(ii,2,:)));
                cc = xcorr(e_sir1,e_sir2);
                [IACC(ii),idx_lag] = max(abs(cc));
                toa_diff(ii) = idx_lag - Ns;
            end
            if flags.do_guesstoa
                toa = guesstoa(toa_diff,toa, kv.avgtoa);
            end
            

        case 'CenIACCr'
            amt_disp('CenIACCr mode \n',flags.disp)
            x = transpose(1:(Ns*2-1));
            for ii=1:pos                
                cc = xcorr(squeeze(f_ir(ii,1,:)),squeeze(f_ir(ii,2,:)));
                pos_cc = abs(cc);
                IACC(ii) = max(pos_cc);
                toa_diff(ii) = centroid(x,pos_cc)-Ns;
            end
            if flags.do_guesstoa
                toa = guesstoa(toa_diff,toa, kv.avgtoa);
            end
            

        case 'CenIACCe'
            amt_disp('CenIACCe mode \n',flags.disp)
            x = transpose(1:(Ns*2-1));
            for ii=1:pos
                e_sir1 = envelope(squeeze(f_ir(ii,1,:)));
                e_sir2 = envelope(squeeze(f_ir(ii,2,:)));
                cc = xcorr(e_sir1,e_sir2);
                IACC(ii) = max(abs(cc));
                toa_diff(ii) = centroid(x,abs(cc))-Ns;
            end
            if flags.do_guesstoa
                toa = guesstoa(toa_diff,toa, kv.avgtoa);
            end
            

        case 'CenIACC2e'
            amt_disp('CenIACC2e mode \n',flags.disp)
            x = transpose(1:(Ns*2-1));
            for ii=1:pos              
                e_sir1 = envelope(squeeze(f_ir(ii,1,:)));
                e_sir2 = envelope(squeeze(f_ir(ii,2,:)));           
                cc = xcorr(e_sir1,e_sir2).^2;
                IACC(ii) = max(abs(cc));
                toa_diff(ii) = centroid(x,abs(cc))-Ns;    
            end
            if flags.do_guesstoa
                toa = guesstoa(toa_diff,toa, kv.avgtoa);
            end
            

        case 'PhminXcor'
            amt_disp('PhminXcor mode \n',flags.disp)
            ir_min=ARI_MinimalPhase(Obj);
            for ii=1:pos
                for jj=1:ear                    
                    cc = xcorr(squeeze(IR(ii,jj,:)),squeeze(ir_min(ii,jj,:)));
                    [~,toa(ii,jj)] = max(abs(cc));
                end
                toa_diff(ii) = toa(ii,1) - toa(ii,2);
            end
            
            
    % ---------------------- Groupdelay -----------------------------------
        case 'IRGD'
            amt_disp('IRGD mode \n',flags.disp)
            for ii = 1:pos
                for jj = 1:ear
                    f_sir = squeeze( f_ir(ii,jj,:) );
                    [gd,w] = grpdelay(transpose(double(f_sir)),1,Ns,fs);
                    toa(ii,jj)=mean(gd(find(w>kv.lower_cutfreq): ...
                                        find(w>kv.upper_cutfreq)));
                end
                toa_diff(ii) = toa(ii,1) - toa(ii,2);
            end
    end
    toa_diff = toa_diff/fs;
    toa = toa/fs;
end



% -------------------------------------------------------------------------
% ---------------------- Functions ----------------------------------------
% -------------------------------------------------------------------------

% ---------------------- Centroid -----------------------------------------
function idx_cent = centroid(x,y)
    idx_cent = sum(x.*y)/sum(y);
end

% ---------------------- guess toa ----------------------------------------
function toa = guesstoa(toa_diff,toa, avgtoa)
    toa(:,1) = toa(:,1) + avgtoa + toa_diff/2;
    toa(:,2) = toa(:,2) + avgtoa - toa_diff/2;
end
    
% ---------------------- Create minimal phase -----------------------------
% as used in ziegelwanger2014
function hMmin=ARI_MinimalPhase(Obj)
    hM=Obj.Data.IR;
    hMmin=hM;

    for jj=1:Obj.API.R
        for ii=1:Obj.API.M
            h=squeeze(hM(ii,jj,:));
            
            amp1=abs(fft(h));
            amp2=amp1;
            
            an2u=-imag(hilbert(log(amp1)));
            an2u=an2u(1:floor(length(h)/2)+1);
            
            an3u=[an2u; -flipud(an2u(2:end+mod(length(h),2)-1))];
            an3=an3u-round(an3u/2/pi)*2*pi;
            
            amp2=amp2(1:floor(length(h)/2)+1);
            amp3=[amp2; flipud(amp2(2:end+mod(length(h),2)-1))];
            
            h2=real(ifft(amp3.*exp(1i*an3)));
            hMmin(ii,jj,:)=h2(1:Obj.API.N);
        end
    end
end