This documentation page applies to an outdated major AMT version. We show it for archival purposes only.
Click here for the documentation menu and here to download the latest AMT (1.6.0).
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/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