function [phi,phi_std,itd,ild,cfreqs] = wierstorf2013_estimateazimuth(insig,lookup,varargin)
%WIERSTORF2013_ESTIMATEAZIMUTH Estimate the perceived azimuth using a binaural model
% Usage: [phi,phi_std,itd,ild,cfreqs] = wierstorf2013_estimateazimuth(sig,lookup,'fs',44100,'dietz2011')
% [phi,phi_std,itd,ild,cfreqs] = wierstorf2013_estimateazimuth(sig,lookup)
%
% Input parameters:
% sig : binaural singal
% lookup : lookup table to map ITDs to angles (struct)
%
% Output parameters:
% phi : estimated azimuth / deg
% phi_std : standard deviation of the estimated azimuth / deg
% itd : calculated ITD (s)
% ild : calculated ILD (dB)
% cfreqs : center frequencies of used auditory filters (Hz)
%
% WIERSTORF2013_ESTIMATEAZIMUTH(sig,lookup) uses a binaural model to
% estimate the perceived direction for a given binaural signal. Therefore,
% it needs the struct lookup, which maps ITD values to the corresponding
% angles. This can be created with the ITD2ANGLE_LOOKUPTABLE function.
% The azimuth values are first calculated for every frequency channel and
% after that their median is calculated. In this process the different
% frequency channels could be weighted and outlier could be removed, see the
% options below. The default setting does not apply any weighting of the
% frequency channels and removes outlier that deviate more than 30 degree from
% the median.
%
% WIERSTORF2013_ESTIMATEAZIMUTH accepts the following optional parameters:
%
% 'fs',fs Sampling rate
%
% 'dietz2011' Use the dietz2011 binaural model to estimate the
% azimuth value. This is the default.
%
% 'lindemann1986' Use the lindemann1986 binaural model to estimate
% the azimuth value.
%
% 'no_spectral_weighting' Apply equal weighting of all frequency channels.
% This is the default behavior.
%
% 'rms_weighting' Weight the frequency channels according their rms
% value of the signal.
%
% 'raatgever_weighting' Weight the frequency channels after the empirical
% curve from Raatgever that has a maximum around
% 600 Hz. Note, that this works well only for
% special stimuli.
%
% 'remove_outlier' Remove frequency channels from azimuth
% calculation that deviate more than 30 degree from
% the median azimuth.
%
% 'include_outlier' Use all azimuth values to calculate the median.
% Note, this can lead to NaN if one of the
% frequency channels has a NaN as direction.
%
% See also: wierstorf2013, dietz2011, lindemann1986, itd2angle_lookuptable
%
% References:
% M. Dietz, S. D. Ewert, and V. Hohmann. Auditory model based direction
% estimation of concurrent speakers from binaural signals. Speech
% Communication, 53(5):592--605, 2011.
%
% W. Lindemann. Extension of a binaural cross-correlation model by
% contralateral inhibition. I. Simulation of lateralization for
% stationary signals. The Journal of the Acoustical Society of America,
% 80:1608--1622, 1986.
%
% J. Raatgever. On the binaural processing of stimuli with different
% interaural phase relations. PhD thesis, TU Delft, 1980.
%
% R. Stern, A. Zeiberg, and C. Trahiotis. Lateralization of complex
% binaural stimuli: A weighted-image model. The Journal of the Acoustical
% Society of America, 84:156--165, 1988.
%
% H. Wierstorf, A. Raake, and S. Spors. Binaural assessment of
% multi-channel reproduction. In J. Blauert, editor, The technology of
% binaural listening, chapter 10. Springer, Berlin--Heidelberg--New York
% NY, 2013.
%
%
% Url: http://amtoolbox.org/amt-1.6.0/doc/modelstages/wierstorf2013_estimateazimuth.php
% #StatusDoc: Perfect
% #StatusCode: Perfect
% #Verification: Verified
% #Requirements: SOFA SFS
% #Author: Hagen Wierstorf
% This file is licensed unter the GNU General Public License (GPL) either
% version 3 of the license, or any later version as published by the Free Software
% Foundation. Details of the GPLv3 can be found in the AMT directory "licences" and
% at <https://www.gnu.org/licenses/gpl-3.0.html>.
% You can redistribute this file and/or modify it under the terms of the GPLv3.
% This file is distributed without any warranty; without even the implied warranty
% of merchantability or fitness for a particular purpose.
% AUTHOR: Hagen Wierstorf
%% ===== Checking of input parameters ==================================
if nargin<2
error('%s: Too few input parameters.',upper(mfilename));
end;
if ~isnumeric(insig) || min(size(insig))~=2
error('%s: insig has to be a numeric two channel signal!',upper(mfilename));
end
if ~isstruct(lookup)
error('%s: lookup has to be a struct!',upper(mfilename));
end
definput.keyvals.fs = 44100;
definput.flags.binaural_model = {'dietz2011','lindemann1986'};
definput.flags.spectral_weighting = {'no_spectral_weighting','rms_weighting','raatgever_weighting'};
definput.flags.outlier = {'remove_outlier','include_outlier'};
[flags,kv] = ltfatarghelper({},definput,varargin);
%% ===== Computation ====================================================
%
% === Calculate azimuth values for every frequency channel ===
if flags.do_dietz2011
ic_threshold=0.98;
% Run the Dietz model on signal
[fine,cfreqs,ild] = dietz2011(insig,kv.fs,'nolowpass','fhigh',1400);
% Unwrap ITDs and get the azimuth values
itd = dietz2011_unwrapitd(fine.itd,ild,fine.f_inst,2.5);
phi = itd2angle(itd,lookup);
% Calculate the median over time for every frequency channel of the azimuth
for n = 1:size(phi,2)
idx = fine.ic(:,n)>ic_threshold&[diff(fine.ic(:,n))>0; 0]; % compare eq. 9 in Dietz (2011)
angle = phi(idx,n);
idx = ~isnan(angle);
if size(angle(idx),1)==0
azimuth(n) = NaN;
azimuth_std(n) = NaN;
else
azimuth(n) = median(angle(idx));
azimuth_std(n) = std(angle(idx));
end
end
% Calculate ITD and ILD values
itd = median(itd,1);
% weights for rms-weighting
rms_weights = fine.rms;
elseif flags.do_lindemann1986
% run Lindemann model on signal
c_s = 0.3; % stationary inhibition
w_f = 0; % monaural sensitivity
M_f = 6; % decrease of monaural sensitivity
T_int = 6; % integration time
N_1 = 1764; % sample at which first cross-correlation is calculated
[cc_tmp,~,ild,cfreqs] = lindemann1986(insig,kv.fs,c_s,w_f,M_f,T_int,N_1);
cc_tmp = squeeze(cc_tmp);
% Calculate tau (delay line time) axes
tau = linspace(-1,1,size(cc_tmp,2));
% find max in cc
itd = zeros(size(cc_tmp,1),size(cc_tmp,3));
for ii=1:size(cc_tmp,1)
for jj=1:size(cc_tmp,3)
[~,idx] = max(cc_tmp(ii,:,jj));
itd(ii,jj) = tau(idx)/1000;
end
end
azimuth = itd2angle(itd,lookup);
azimuth_std = std(azimuth);
azimuth = median(azimuth);
% TODO: rms weights
end
% === Weights for cross-frequency integration ===
if flags.do_no_spectral_weighting
w = ones(size(azimuth));
end
if flags.do_rms_weighting
w = rms_weights;
end
if flags.do_raatgever_weighting
% Calculate a spectral weighting after Stern1988, after the data of
% Raatgever1980
b1 = -9.383e-2;
b2 = 1.126e-4;
b3 = -3.992e-8;
w = 10.^( -(b1*f+b2*(f).^2+b3*(f).^3)/10 );
end
if flags.do_remove_outlier
% Remove outliers
[azimuth,azimuth_std,itd,cfreqs,w] = remove_outlier(azimuth,azimuth_std,itd,cfreqs,w);
end
% Calculate mean about frequency channels
if length(azimuth)==0
phi = NaN;
phi_std = NaN;
else
if flags.do_no_spectral_weighting
phi = median(azimuth);
phi_std = median(azimuth_std);
else
phi = sum(azimuth.*w)/sum(w);
phi_std = sum(azimuth_std.*w)/sum(w);
end
end
end % of main function
%% ===== Subfunctions ====================================================
function [azimuth,azimuth_std,itd,cfreqs,w] = remove_outlier(azimuth,azimuth_std,itd,cfreqs,w)
% NOTE: the following was enabled for the original paper
% remove unvalid ITDs
%azimuth = azimuth(abs(itd)<0.001);
%azimuth_std = azimuth_std(abs(itd)<0.001);
%cfreqs = cfreqs(abs(itd)<0.001);
%w = w(abs(itd)<0.001);
%itd = itd(abs(itd)<0.001);
% remove NaN
w = w(~isnan(azimuth));
itd = itd(~isnan(azimuth));
cfreqs = cfreqs(~isnan(azimuth));
azimuth = azimuth(~isnan(azimuth));
azimuth_std = azimuth_std(~isnan(azimuth_std));
% remove outliers more than 30deg away from median
if length(azimuth)>0
itd = itd(abs(azimuth-median(azimuth))<30);
cfreqs = cfreqs(abs(azimuth-median(azimuth))<30);
w = w(abs(azimuth-median(azimuth))<30);
azimuth_std = azimuth_std(abs(azimuth-median(azimuth))<30);
azimuth = azimuth(abs(azimuth-median(azimuth))<30);
end
end