The Auditory Modeling Toolbox
Go to function
EXP_BAUMGARTNER2017 - - Experiments of Baumgartner et al. (2017)
Program code:
function data = exp_baumgartner2017(varargin)
%EXP_BAUMGARTNER2017 - Experiments of Baumgartner et al. (2017)
% Usage: data = exp_baumgartner2017(flag)
%
% EXP_BAUMGARTNER2017(flag) reproduces figures of the study from
% Baumgartner et al. (2017).
%
% The following flags can be specified
%
% boyd2012 models experiments from Boyd et al. (2012; Fig.1, top).
% Average externalization ratings of 1 talker for NH participants
% against mix point as a function of microphone position (ITE/BTE)
% and frequency response (BB/LP). The reference condition (ref) is
% the same as ITE/BB. Error bars show SEM.
%
% hartmann1996 models experiments from Hartmann & Wittenberg (1996; Fig.7-8)
% 1st panel: Synthesis of zero-ILD signals. Only the harmonics
% from 1 to nprime had zero interaural level difference;
% harmonics above nprime retained the amplitudes of the baseline
% synthesis. Externalization scores as a function of the boundary
% harmonic number nprime. Fundamental frequency of 125 Hz.
% 2nd panel: Synthesis of signals to test the ISLD hypothesis.
% Harmonics at and below the boundary retained only the interaural
% spectral level differences of the baseline synthesis. Higher
% harmonics retained left and right baseline harmonic levels.
% Externalization scores as a function of the boundary
% frequency.
%
% hassager2016 models experiments from Hassager et al. (2016; Fig.6).
% The mean of the seven listeners perceived sound source
% location (black) as a function of the bandwidth factor
% and the corresponding model predictions (colored).
% The model predictions have been shifted slightly to the right
% for a better visual interpretation. The error bars are one
% standard error of the mean.
%
% baumgartner2017 models experiments from Baumgartner et al. (2017).
% Effect of HRTF spectral contrast manipulations on sound externalization.
% Externalization scores were derived from paired comparisons via
% Bradley-Terry-Luce modeling.
%
% Requirements:
% -------------
%
% 1) SOFA API v0.4.3 or higher from http://sourceforge.net/projects/sofacoustics for Matlab (in e.g. thirdparty/SOFA)
%
% 2) Data in hrtf/baumgartner2017
%
% 3) Statistics Toolbox for Matlab (for some of the figures)
%
% Examples:
% ---------
%
% To display results for Fig.1 from Boyd et al. (2012) use :
%
% exp_baumgartner2017('boyd2012');
%
% To display results for Fig.7-8 from Hartmann & Wittenberg (1996) use :
%
% exp_baumgartner2017('hartmann1996');
%
% To display results for Fig.6 from Hassager et al. (2016) use :
%
% exp_baumgartner2017('hassager2016');
%
% References:
% R. Baumgartner, P. Majdak, H. Colburn, and B. Shinn-Cunningham.
% Modeling sound externalization based on listener-specific spectral
% cues. In Acoustics ‘17 Boston: The 3rd Joint Meeting of the Acoustical
% Society of America and the European Acoustics Association, Boston, MA,
% Jun 2017.
%
% A. W. Boyd, W. M. Whitmer, J. J. Soraghan, and M. A. Akeroyd. Auditory
% externalization in hearing-impaired listeners: The effect of pinna cues
% and number of talkers. J. Acoust. Soc. Am., 131(3):EL268--EL274, 2012.
% [1]arXiv | [2]www: ]
%
% W. M. Hartmann and A. Wittenberg. On the externalization of sound
% images. J. Acoust. Soc. Am., 99(6):3678--88, June 1996.
%
% References
%
% 1. http://arxiv.org/abs/ http://dx.doi.org/10.1121/1.3687015
% 2. http://dx.doi.org/10.1121/1.3687015
%
% hassager2016 baumgartner2017
%
% Url: http://amtoolbox.sourceforge.net/data/amt-test/htdocs/amt-0.10.0/doc/experiments/exp_baumgartner2017.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/>.
% AUTHOR: Robert Baumgartner, Acoustics Research Institute, Vienna, Austria
definput.import={'amt_cache','baumgartner2017'};
definput.flags.type = {'missingflag','boyd2012','hartmann1996','hassager2016','baumgartner2017'};
definput.flags.quickCheck = {'','quickCheck'};
definput.keyvals.Sintra = 2;
definput.keyvals.Sinter = 2;
definput.keyvals.Interaction = {'MSG','ITIT'};
[flags,kv] = ltfatarghelper({},definput,varargin);
if flags.do_missingflag
flagnames=[sprintf('%s, ',definput.flags.type{2:end-2}),...
sprintf('%s or %s',definput.flags.type{end-1},definput.flags.type{end})];
error('%s: You must specify one of the following flags: %s.',upper(mfilename),flagnames);
end
Ncues = 7;
% Symbol order for plotting
symb = {'-o','-s','-d','-<','->','-^','-v','-p',':*','--x'};
colors = 0.8*[zeros(1,3);hsv(Ncues+2)];
InteractionLbl = [kv.Interaction{1},' x ',kv.Interaction{2}];
%% Hassager et al. (2016)
if flags.do_hassager2016
azi = [0,50];
flp = 6000; % lowpass filtered noise stimuli
Pext_A = data_hassager2016;
B = Pext_A.B;
fncache = ['hassager2016'];
[Pext,cues,cueLbl] = amt_cache('get',fncache,flags.cachemode);
if isempty(cues)
data = data_baumgartner2017;
if flags.do_quickCheck
data = data(1:5);
end
% Pext = nan(length(B),length(data),length(azi));
% Pext = {Pext,Pext};
cues = nan(length(B),length(azi),length(data),Ncues);
for isubj = 1:length(data)
Obj = data(isubj).Obj;
for iazi = 1:length(azi)
idazi = Obj.SourcePosition(:,1) == azi(iazi) & Obj.SourcePosition(:,2) == 0;
template = squeeze(shiftdim(Obj.Data.IR(idazi,:,:),2));
for iB = 1:length(B)
amt_disp(num2str(iB),'volatile');
if isnan(B(iB))
target = template;
else
Obj_tar = sig_hassager2016(Obj,B(iB));
target = squeeze(shiftdim(Obj_tar.Data.IR(idazi,:,:),2));
end
[~,cues(iB,iazi,isubj,:),cueLbl] = baumgartner2017(target,template,'flow',100,'fhigh',flp,'lat',azi(iazi));
end
end
amt_disp([num2str(isubj),' of ',num2str(length(data)),' subjects completed.'],'progress')
end
if not(flags.do_quickCheck)
amt_cache('set',fncache,Pext,cues,cueLbl);
end
end
%% Cue weighting optimization procedure
dimSubj = 3;
dimCues = 4;
Nsubj = size(cues,dimSubj);
cues = num2cell(cues,1:dimCues-1);
% optimize sensitivities
% predCue = reshape(cues,length(B)*length(azi)*Nsubj,Ncues);
actual = repmat(Pext_A.rating(:),Nsubj,1);
Erange = 4;
Eoffset = 1;
InteractionIds = find(ismember(cueLbl(:,1),kv.Interaction));
[S,W,Pext,dE] = cueWeightOpti(cues,actual,Erange,Eoffset,InteractionIds);
PextLbl = [cueLbl(:,1);{'Comb.'};{InteractionLbl}];
cueTab = table([S;nan(2,1)],dE,[W;nan(2,1)],[cueLbl(:,2);{'Weighted combination'};{[InteractionLbl, ' interaction']}],...
'RowNames',PextLbl,...
'VariableNames',{'Sensitivity','Error','Weight','Description'});
amt_disp(flags.type)
amt_disp(cueTab)
%% Plot
BplotTicks = logspace(log10(0.25),log10(64),9);
BplotTicks = round(BplotTicks*100)/100;
BplotStr = ['Ref.';num2str(BplotTicks(:))];
BplotTicks = [BplotTicks(1)/2,BplotTicks];
B(1) = B(2)/2;
dataLbl = [{'Actual'};PextLbl];
idleg = not(ismember(dataLbl,'ITSD'));
figure
hax = tight_subplot(1,2,0,[.15,.1],[.1,.05]);
for iazi = 1:length(azi)
% subplot(1,2,iazi)
axes(hax(iazi))
h(1) = plot(B,Pext_A.rating(:,iazi),symb{1},'Color',colors(1,:));
hold on
for m = 1:length(Pext)
dx = 1+0.02*(m - (Ncues+1)/2);
h(m+1) = errorbar(dx*B,mean(Pext{m}(:,iazi,:),dimSubj),std(Pext{m}(:,iazi,:),0,dimSubj)/sqrt(Nsubj),...
symb{m+1},'Color',colors(m+1,:));
end
set(h,'MarkerFaceColor','w')
set(gca,'XTick',BplotTicks,'XTickLabel',BplotStr,'XScale','log')
axis([BplotTicks(1)/1.5,BplotTicks(end)*1.5,0.8,5.2])
xlabel('Bandwidth Factor [ERB]')
if iazi==1
ylabel('Mean Externalization Rating')
else
set(gca,'YTickLabel',{})
end
title([num2str(azi(iazi)),'\circ'])
end
leg = legend(h(idleg),dataLbl(idleg),'Location','southwest');
set(leg,'Box','off')
%% Output
data = {cueTab,Pext};
end
%% Hartmann & Wittenberg (1996)
if flags.do_hartmann1996
cond = {'ILD','ISLD'};
condPlotLbl = {'ILD to 0','ILD const.'};
nprime = [0,1,8,14,19,22,25,38];
fncache = ['hartmann1996'];
[Pext,cues,cueLbl] = amt_cache('get',fncache,flags.cachemode);
if isempty(Pext)
azi = -37;
data = data_baumgartner2017;
if flags.do_quickCheck
data = data(1:5);
end
Pext = repmat({nan(length(data),length(nprime))},2,2);
cues = repmat({nan(length(data),length(nprime),Ncues)},2,1);
for isub = 1:length(data)
Obj = data(isub).Obj;
template = sig_hartmann1996(0,'Obj',Obj,'dur',0.1);
for ee = 1:length(cond)
for nn = 1:length(nprime)
target = sig_hartmann1996(nprime(nn),'Obj',Obj,'dur',0.1,cond{ee});
[~,cues{ee}(isub,nn,:),cueLbl] = baumgartner2017(target,template,'flow',100,'fhigh',5000,'lat',azi);
end
end
amt_disp([num2str(isub),' of ',num2str(length(data)),' subjects completed.'],'progress')
end
if not(flags.do_quickCheck)
amt_cache('set',fncache,Pext,cues,cueLbl);
end
end
%% Cue weighting optimization procedure
% align actual and modeled data format
actual = [];
predCue = [];
for ee = 1:length(cond)
act = data_hartmann1996(cond{ee});
actual = cat(1,actual,act.avg.Escore(:));
tmp = squeeze(mean(cues{ee},1));
tmp = interp1(nprime,tmp,act.avg.nprime);
predCue = cat(1,predCue,tmp);
end
predCue = num2cell(predCue,1);
Erange = 3;
Eoffset = 0;
InteractionIds = find(ismember(cueLbl(:,1),kv.Interaction));
if length(InteractionIds) ~=2
error('RB: Check interaction labels.')
end
[S,W,~,dE] = cueWeightOpti(predCue,actual,Erange,Eoffset,InteractionIds);
% calculate separate errors for Exp 1 and 2
dE1 = nan(size(dE));
dE2 = nan(size(dE));
for ss = 1:length(S)
E = cue2E(S(ss),predCue{ss},Erange,Eoffset);
dE1(ss) = nanrms(E(1:5)-actual(1:5)) / (Erange-Eoffset);
dE2(ss) = nanrms(E(6:9)-actual(6:9)) / (Erange-Eoffset);
end
PextLbl = [cueLbl(:,1);{'Comb.'};{InteractionLbl}];
cueTab = table([S;nan(2,1)],dE1,dE2,dE,[W;nan(2,1)],[cueLbl(:,2);{'Weighted combination'};{[InteractionLbl, ' interaction']}],...
'RowNames',PextLbl,...
'VariableNames',{'Sensitivity','ErrorExp1','ErrorExp2','Error','Weight','Description'});
amt_disp(flags.type)
amt_disp(cueTab)
%% Individual data
dimCues = 3;
% Ncues = size(cues{1},dimCues);
Pext = cell(Ncues+1,length(cond));
for ee = 1:length(cond)
for cc = 1:Ncues
Pext{cc,ee} = cue2E(S(cc),cues{ee}(:,:,cc),Erange,Eoffset);
end
Pext{Ncues+1,ee} = W(1)*Pext{1,ee};
for cc = 2:Ncues
Pext{Ncues+1,ee} = nansum(cat(dimCues,Pext{Ncues+1,ee},W(cc)*Pext{cc,ee}),dimCues);
end
Pext{Ncues+2,ee} = sqrt(Pext{InteractionIds(1),ee}.*Pext{InteractionIds(2),ee}); % interaction term
end
%% Plot
dataLbl = [{'Actual'};PextLbl];
idleg = not(ismember(dataLbl,'ITSD'));
Ns = size(Pext{1},1);
figure
hax = tight_subplot(1,2,0,[.15,.1],[.1,.05]);
for ee = 1:length(cond)
act = data_hartmann1996(cond{ee});
% subplot(1,2,ee)
axes(hax(ee))
h(1) = plot(act.avg.nprime,act.avg.Escore,symb{1},'Color',colors(1,:));
hold on
for cc = 1:size(Pext,1)
if idleg(cc+1)
dx = 0.1*(cc - size(Pext,1)/2);
h(cc+1) = errorbar(nprime+dx,mean(Pext{cc,ee}),std(Pext{cc,ee})/sqrt(Ns),...
symb{cc+1},'Color',colors(cc+1,:));
end
end
set(h(idleg),'MarkerFaceColor','w')
if ee == 1
xlabel('n^{\prime} (Highest harmonic with ILD = 0)')
ylabel('Externalization score')
else
xlabel('n^{\prime} (Highest harmonic with altered amplitudes)')
set(gca,'YTickLabel',{})
end
title(condPlotLbl{ee})
axis([-2.5,42.5,-0.4,3.9])
idx = false(size(idleg));
if ee == 1
idx(1:Ncues) = true;
else
idx(Ncues+1:end) = true;
end
legend(h(and(idleg,idx)),dataLbl(and(idleg,idx)),'Location','southwest','box','off');
end
% Output
data = {cueTab,Pext};
end
%% Boyd et al. (2012) -----------------------------------------------------
if flags.do_boyd2012
flags.do_noise = false;
flags.do_georganti2013 = false;
flags.do_BRIR = true;
flp = [nan,6500]; % Low-pass cut-off frequency
azi = -30;
subjects = data_boyd2012;
if flags.do_noise || flags.do_BRIR
subjects = subjects([3,6,7]); % only the ones with BRIRs
else
subjects = subjects([1,3,4,6,7]); % only the ones with stimuli
end
% fprintf(['Note that only 3 of originally 7 listeners are simulated and compared because\n',...
% 'the listener-specific BRIRs for the other 4 listeners are not available.\n'])
mix = subjects(1).Resp.mix/100;
fncache = ['boyd2012_lp'];
[E,cues,cueLbl] = amt_cache('get',fncache,flags.cachemode);
if isempty(cues)
mix(1) = 1;
if flags.do_noise || flags.do_BRIR
sig = noise(50e3,1,'pink');
fs = subjects(1).BRIR.fs;
[b,a]=butter(10,2*flp(2)/fs,'low');
end
E.meta = {'actual','interaural','monaural'};
if flags.do_georganti2013
E.meta(4) = {'BSMD'};
end
E.all = repmat({nan([length(mix),2,2,length(subjects)])},1,length(E.meta));
cues = nan(length(mix),2,2,length(subjects),Ncues);
for isub = 1:length(subjects)
E.all{1}(:,:,:,isub) = cat(3,...
[subjects(isub).Resp.ITE_BB_1T(:),subjects(isub).Resp.BTE_BB_1T(:)],...
[subjects(isub).Resp.ITE_LP_1T(:),subjects(isub).Resp.BTE_LP_1T(:)]);
% ITE = subjects(isub).Obj;
if flags.do_noise
stim{1} = lconv(sig,subjects(isub).BRIR.ITE(:,:,1));
stim{2} = lconv(sig,subjects(isub).BRIR.BTE(:,:,1));
stim{3} = lconv(sig,subjects(isub).BRIR.noHead(:,:,1));
template = stim{1};
targetSet = cell(length(mix),2,2);
elseif flags.do_BRIR
template = subjects(isub).BRIR.ITE(:,:,1);
stim{1} = subjects(isub).BRIR.ITE(:,:,1);
stim{2} = subjects(isub).BRIR.BTE(:,:,1);
stim{3} = subjects(isub).BRIR.noHead(:,:,1);
targetSet = cell(length(mix),2,2);
else % speech samples
template = subjects(isub).Reference_1T;
targetSet = cat(3,...
[subjects(isub).Target.ITE_BB_1T(:),subjects(isub).Target.BTE_BB_1T(:)],...
[subjects(isub).Target.ITE_LP_1T(:),subjects(isub).Target.BTE_LP_1T(:)]);
targetSet = [repmat({template},[1,2,2]);targetSet];
end
%% Model simulations
flow = 100;
fhigh = [16e3,flp(2)];
for c = 1:2
for lp = 1:2
for m = 1:length(mix)
if m==1 && (c~=1 || lp~=1) % same reference condition
for ee = 2:length(E.all)
E.all{ee}(m,c,lp,isub) = E.all{ee}(m,1,1,isub);
targetSet{m,c,lp} = targetSet{m,1,1};
end
else
% target{m,c,lp} = sig_boyd2012(stim{c},sig,mix(m),flp(lp),fs);
% mixing
if flags.do_noise || flags.do_BRIR
targetSet{m,c,lp} = mix(m)*stim{c} + (1-mix(m))*stim{3};
% low-pass filtering
if lp == 2
targetSet{m,c,lp} = filter(b,a,targetSet{m,c,lp});
end
end
target = targetSet{m,c,lp};
% % Description in paper is ambiguous about whether broadband ILDs were
% % adjusted to the template:
% temSPL = dbspl(template);
% for ch = 1:2
% target{m,c,lp}(:,ch) = setdbspl(target{m,c,lp}(:,ch),temSPL(ch));
% end
amt_disp(num2str(mix(m)),'volatile')
if flags.do_BRIR
flags.do_spectralCueEchoSuppression = true;
end
[~,cues(m,c,lp,isub,:),cueLbl] = baumgartner2017(target,template...
,'argimport',flags,kv,'flow',flow,'fhigh',fhigh(1),'lat',azi,'reflectionOnsetTime',5e-3); % replace fhigh(1) with fhigh(lp) for manual bandwidth adjustment
if flags.do_georganti2013
geor.fs = subjects(1).fs;
geor.timeFr = 1;
geor.fmin = flow;
geor.fmax = fhigh(lp);
E.all{4}(m,c,lp,isub) = median(georganti2013(target,geor)); % median over time
end
end
end
end
end
amt_disp([num2str(isub),' of ',num2str(length(subjects)),' subjects completed.'],'progress')
end
if flags.do_georganti2013
E.all{4} = 100*E.all{4}./repmat(E.all{4}(1,1,1,:),[length(mix),2,2,1]);
end
if not(flags.do_quickCheck)
amt_cache('set',fncache,E,cues,cueLbl);
end
end
%% Cue weighting optimization procedure
dimSubj = 4;
dimCues = 5;
Erange = 100;
Eoffset = 0;
Nsubj = size(cues,dimSubj);
if flags.do_BRIR % comparison on individual basis
actual = E.all{1};
else % average comparison
actual = repmat(mean(E.all{1},4),[ones(1,Nsubj-1),Nsubj]);
end
cuesC = squeeze(num2cell(cues,1:dimCues-1));
InteractionIds = find(ismember(cueLbl(:,1),kv.Interaction));
[S,W,Pext,dE] = cueWeightOpti(cuesC,actual,Erange,Eoffset,InteractionIds);
E.all(2:Ncues+3) = Pext;
PextLbl = [cueLbl(:,1);{'Comb.'};{InteractionLbl}];
cueTab = table([S;nan(2,1)],dE,[W;nan(2,1)],[cueLbl(:,2);{'Weighted sum'};{[InteractionLbl, ' interaction']}],...
'RowNames',PextLbl,...
'VariableNames',{'Sensitivity','Error','Weight','Description'});
amt_disp(flags.type)
amt_disp(cueTab)
%% Average data
for ee = 1:length(E.all)
E.m{ee} = mean(E.all{ee},4);
E.se{ee} = std(E.all{ee},0,4);
end
%% Plot
flags.do_plot_individual = false;
condLbl = {'ITE / BB','BTE / BB','ITE / LP','BTE / LP'};
mix = mix(2:end)*100;
leglbl = [{'Actual'};PextLbl];
idleg = not(ismember(leglbl,'ITSD'));
if not(flags.do_plot_individual)
figure
hax = tight_subplot(2,2,0,[.15,.1],[.1,.05]);
for cc = 1:length(condLbl)
axes(hax(cc))
for ee = 1:length(E.m)
if idleg(ee)
dx = 2*(ee - length(E.m)/2);
h = errorbar(mix+dx,E.m{ee}(2:end,cc),E.se{ee}(2:end,cc),...
symb{ee},'Color',colors(ee,:));
set(h,'MarkerFaceColor','w')
hold on
end
end
set(gca,'XDir','reverse')
if cc == 3 || cc == 4
xlabel('Mix (%)')
else
set(gca,'XTickLabel',[])
end
if cc == 1 || cc == 3
ylabel('Externalization score')
else
set(gca,'YTickLabel',[])
end
% title(condLbl{cc})
text(110,0,condLbl{cc})
axis([-20,120,-15,115])
if cc == 4
% leg = legend(leglbl(idleg));
% set(leg,'Box','off','Location','eastoutside')
end
end
else % flags.do_plot_individual
for isub = 1:size(E.all{1},4)
figure
hax = tight_subplot(1,4,0,[.15,.1],[.1,.05]);
for cc = 1:length(condLbl)
axes(hax(cc))
for ee = 1:Nee
Eisub = E.all{ee}(:,:,:,isub);
h = plot(mix,Eisub(2:end,cc),symb{ee},'Color',colors(ee,:));
set(h,'MarkerFaceColor','w')
hold on
end
set(gca,'XDir','reverse')
xlabel('Mix (%)')
if cc == 1
ylabel('Externalization score')
else
set(gca,'YTickLabel',[])
end
title(condLbl{cc})
axis([-20,120,-5,105])
if cc == 4
leg = legend(dataLbl);
set(leg,'Box','off','Location','north')
end
end
end
end
% set(leg,'Location','eastoutside','Position',get(leg,'Position')+[.1,.2,0,0])
%% Output
data = {cueTab,E};
end
if flags.do_baumgartner2017
% Behavioral results
fncache = 'baumgartner2017_E_BTL';
data = amt_cache('get',fncache,flags.cachemode);
if isempty(data)
exp2 = data_baumgartner2017looming('exp2');
Nsubj = length(exp2.meta.subject);
iISI = exp2.meta.ISI == 0.1;
iresp = strcmp(exp2.meta.response,'approaching');
A = {1;2;3}; % arbitrary indices to define BTL model
% Alabel = {'C = 0','C = 0.5','C = 1'};
data.C = 0:0.5:1;
data.subj = [];%exp2.meta.subject;
data.Escore = [];%nan(Nsubj,length(A));
data.BTL_chistat = [];%nan(Nsubj,1);
data.azi = [exp2.rawData.azimuth];
M = zeros(3,3);
iM = [7,4,8,3,2,6];
iss = 0;
for ss = 1:Nsubj
M(iM) = exp2.data(ss,1:6,iresp,iISI)+eps;
[E,chistat,~,lL] = OptiPt(M,A);
if chistat(1) < 10 && chistat(1) > 0.001
iss = iss+1;
data.subj{iss} = exp2.meta.subject{ss};
data.BTL_chistat(iss) = chistat(1);
data.BTL_logLikeli(iss) = lL;
data.Escore(iss,:) = E/E(3);
end
end
amt_cache('set',fncache,data)
end
Nsubj = length(data.subj);
hrtf = data_baumgartner2017looming('exp2','hrtf');
% Modeling
fncache = ['baumgartner2017'];
[Pext,cues,cueLbl] = amt_cache('get',fncache,flags.cachemode);
if isempty(cues)
Pext = nan(length(data.C),Nsubj);
cues = nan(length(data.C),Nsubj,Ncues);
for isubj = 1:Nsubj
template = hrtf(ismember({hrtf.id},data.subj(isubj))).Obj;
for iC = 1:length(data.C)
target = sig_baumgartner2017looming(template,data.C(iC));
[~,cues(iC,isubj,:),cueLbl] = baumgartner2017(target,template,...
'lat',data.azi(isubj),'flow',1000,'fhigh',18e3);
end
amt_disp([num2str(isubj),' of ',num2str(Nsubj),' subjects completed.'],'progress')
end
amt_cache('set',fncache,Pext,cues,cueLbl);
end
%% Cue weighting optimization procedure
dimSubj = 2;
dimCues = 3;
Erange = 1;
Eoffset = 0;
Nsubj = size(cues,dimSubj);
actual = permute(data.Escore,[2,1]);
cuesC = squeeze(num2cell(cues,1:dimCues-1));
InteractionIds = find(ismember(cueLbl(:,1),kv.Interaction));
[S,W,Pext,dE] = cueWeightOpti(cuesC,actual,Erange,Eoffset,InteractionIds);
PextLbl = [cueLbl(:,1);{'Comb.'};{InteractionLbl}];
cueTab = table([S;nan(2,1)],dE,[W;nan(2,1)],[cueLbl(:,2);{'Weighted sum'};{[InteractionLbl, ' interaction']}],...
'RowNames',PextLbl,...
'VariableNames',{'Sensitivity','Error','Weight','Description'});
amt_disp(flags.type)
amt_disp(cueTab)
%% Figure
flags.do_plot_individual = false;
E_all = [{data.Escore'};Pext];
dataLbl = [{'Actual'};PextLbl];
figure
if not(flags.do_plot_individual)
for iE = 1:length(E_all)
dx = 0.01*(iE - length(E_all)/2);
h(iE) = errorbar(data.C+dx,mean(E_all{iE},2),std(E_all{iE},0,2)/sqrt(Nsubj),...
symb{iE},'Color',colors(iE,:));
set(h,'MarkerFaceColor','w')
hold on
end
set(gca,'XTick',data.C)
axis([-0.2,1.2,-0.1,1.1])
ylabel('Externalization')
xlabel('Spectral contrast, C')
else % flags.do_plot_individual
for isubj = 1:Nsubj
subplot(3,ceil(Nsubj/3),isubj)
for iE = 1:length(E_all)
dx = 0.01*(iE - length(E_all)/2);
h(iE) = plot(data.C+dx,E_all{iE}(:,isubj),...
symb{iE},'Color',colors(iE,:));
set(h,'MarkerFaceColor','w')
hold on
end
set(gca,'XTick',data.C)
axis([-0.2,1.2,-0.1,1.1])
ylabel('Externalization')
xlabel('Spectral contrast, C')
end
end
idleg = not(ismember(dataLbl,'ITSD'));
legend(h(idleg),dataLbl(idleg),'Location','EastOutside','box','off')
%% Output
data = {cueTab,Pext};
end
end
%%%%%%%%%% INTERNAL FUNCTIONS FOR VISUALIZATION %%%%%%%%%%
function ha = tight_subplot(Nh, Nw, gap, marg_h, marg_w)
% tight_subplot creates "subplot" axes with adjustable gaps and margins
%
% ha = tight_subplot(Nh, Nw, gap, marg_h, marg_w)
%
% in: Nh number of axes in hight (vertical direction)
% Nw number of axes in width (horizontaldirection)
% gap gaps between the axes in normalized units (0...1)
% or [gap_h gap_w] for different gaps in height and width
% marg_h margins in height in normalized units (0...1)
% or [lower upper] for different lower and upper margins
% marg_w margins in width in normalized units (0...1)
% or [left right] for different left and right margins
%
% out: ha array of handles of the axes objects
% starting from upper left corner, going row-wise as in
% going row-wise as in
%
% Example: ha = tight_subplot(3,2,[.01 .03],[.1 .1],[.1 .1])
% for ii = 1:6; axes(ha(ii)); plot(randn(10,ii)); end
% set(ha(1:4),'XTickLabel',''); set(ha,'YTickLabel','')
% Pekka Kumpulainen 20.6.2010 @tut.fi
% Tampere University of Technology / Automation Science and Engineering
if nargin<3; gap = .02; end
if nargin<4 || isempty(marg_h); marg_h = .05; end
if nargin<5; marg_w = .05; end
if numel(gap)==1;
gap = [gap gap];
end
if numel(marg_w)==1;
marg_w = [marg_w marg_w];
end
if numel(marg_h)==1;
marg_h = [marg_h marg_h];
end
axh = (1-sum(marg_h)-(Nh-1)*gap(1))/Nh;
axw = (1-sum(marg_w)-(Nw-1)*gap(2))/Nw;
py = 1-marg_h(2)-axh;
ha = zeros(Nh*Nw,1);
ii = 0;
for ih = 1:Nh
px = marg_w(1);
for ix = 1:Nw
ii = ii+1;
ha(ii) = axes('Units','normalized', ...
'Position',[px py axw axh], ...
'XTickLabel','', ...
'YTickLabel','');
px = px+axw+gap(2);
end
py = py-axh-gap(1);
end
end
function [p,chistat,u,lL_eba,lL_sat,fit,cova] = OptiPt(M,A,s)
% OptiPt parameter estimation for BTL/Pretree/EBA models
% p = OptiPt(M,A) estimates the parameters of a model specified
% in A for the paired-comparison matrix M. M is a matrix with
% absolute frequencies. A is a cell array.
%
% [p,chistat,u] = OptiPt(M,A) estimates parameters and reports
% the chi2 statistic as a measure of goodness of fit. The vector
% of scale values is stored in u.
%
% [p,chistat,u,lL_eba,lL_sat,fit,cova] = OptiPt(M,A,s) estimates
% parameters, checks the goodness of fit, computes the scale values,
% reports the log-likelihoods of the model specified in A and of the
% saturated model, returns the fitted values and the covariance
% matrix of the parameter estimates. If defined, s is the starting
% vector for the estimation procedure. Otherwise each starting value
% is set to 1/length(p).
% The minimization algorithm used is FMINSEARCH.
%
% Examples
% Given the matrix M =
% 0 36 35 44 25
% 19 0 31 37 20
% 20 24 0 46 24
% 11 18 9 0 13
% 30 35 31 42 0
%
% A BTL model is specified by A = {[1];[2];[3];[4];[5]}
% Parameter estimates and the chi2 statistic are obtained by
% [p,chistat] = OptiPt(M,A)
%
% A Pretree model is specified by A = {[1 6];[2 6];[3 7];[4 7];[5]}
% A starting vector is defined by s = [2 2 3 4 4 .5 .5]
% Parameter estimates, the chi2 statistic, the scale values, the
% log-likelihoods of the Pretree model and of the saturated model,
% the fitted values, and the covariance matrix are obtained by
% [p,chistat,u,lL_eba,lL_sat,fit,cova] = OptiPt(M,A,s)
%
% Authors: Florian Wickelmaier (wickelmaier@web.de) and Sylvain Choisel
% Last mod: 03/JUL/2003
% For detailed information see Wickelmaier, F. & Schmid, C. (2004). A Matlab
% function to estimate choice model parameters from paired-comparison data.
% Behavior Research Methods, Instruments, and Computers, 36(1), 29-40.
I = length(M); % number of stimuli
mmm = 0;
for i = 1:I
mmm = [mmm max(A{i})];
end
J = max(mmm); % number of pt parameters
if(nargin == 2)
p = ones(1,J)*(1/J); % starting values
elseif(nargin == 3)
p = s;
end
for i = 1:I
for j = 1:I
diff{i,j} = setdiff(A{i},A{j}); % set difference
end
end
p = fminsearch(@ebalik,p,optimset('Display','iter','MaxFunEvals',10000,...
'MaxIter',10000),M,diff,I); % optimized parameters
lL_eba = -ebalik(p,M,diff,I); % likelihood of the specified model
lL_sat = 0; % likelihood of the saturated model
for i = 1:I-1
for j = i+1:I
lL_sat = lL_sat + M(i,j)*log(M(i,j)/(M(i,j)+M(j,i)))...
+ M(j,i)*log(M(j,i)/(M(i,j)+M(j,i)));
end
end
fit = zeros(I); % fitted PCM
for i = 1:I-1
for j = i+1:I
fit(i,j) = (M(i,j)+M(j,i))/(1+sum(p(diff{j,i}))/sum(p(diff{i,j})));
fit(j,i) = (M(i,j)+M(j,i))/(1+sum(p(diff{i,j}))/sum(p(diff{j,i})));
end
end
chi = 2*(lL_sat-lL_eba);
df = I*(I-1)/2 - (J-1);
chistat = [chi df]; % 1-chi2cdf(chi,df)]; % goodness-of-fit statistic
u = sum(p(A{1})); % scale values
for i = 2:I
u = [u sum(p(A{i}))];
end
H = hessian('ebalik',p',M,diff,I);
C = inv([H ones(J,1); ones(1,J) 0]);
cova = C(1:J,1:J);
end
function lL_eba = ebalik(p,M,diff,I) % computes the likelihood
if min(p)<=0 % bound search space
lL_eba = inf;
return
end
thesum = 0;
for i = 1:I-1
for j = i+1:I
thesum = thesum + M(i,j)*log(1+sum(p(diff{j,i}))/sum(p(diff{i,j})))...
+ M(j,i)*log(1+sum(p(diff{i,j}))/sum(p(diff{j,i})));
end
end
lL_eba = thesum;
end
function H = hessian(f,x,varargin) % computes numerical Hessian
% based on a solution posted on Matlab Central by Paul L. Fackler
k = size(x,1);
fx = feval(f,x,varargin{:});
h = eps.^(1/3)*max(abs(x),1e-2);
xh = x+h;
h = xh-x;
ee = sparse(1:k,1:k,h,k,k);
g = zeros(k,1);
for i = 1:k
g(i) = feval(f,x+ee(:,i),varargin{:});
end
H = h*h';
for i = 1:k
for j = i:k
H(i,j) = (feval(f,x+ee(:,i)+ee(:,j),varargin{:})-g(i)-g(j)+fx)...
/ H(i,j);
H(j,i) = H(i,j);
end
end
end
function [S,W,Pext,dE] = cueWeightOpti(predCue,actual,Erange,Eoffset,InteractionIds)
% optimization routine for cue-specific sensitivities S and optimal
% weighting W of cues. Make sure that actual matches the dimension of each
% cue-specific cell of predCue.
% Definitions of functions
f_cue2E = @(s,cue) cue2E(s,cue,Erange,Eoffset);
f_predErr = @(p) nanrms(p(:)-actual(:)) / (Erange-Eoffset);
% optimize sensitivities
Ncues = length(predCue);
S = ones(Ncues,1); % sensitivity to cue
Pext = cell(Ncues+1,1);
dE = nan(Ncues+1,1);
predE = nan(length(predCue{1}(:)),Ncues);
for cc = 1:Ncues
if not(all(isnan(predCue{cc}(:))))
S(cc) = fminsearch(@(s) f_predErr(f_cue2E(s,predCue{cc})) ,S(cc));
% splot = logspace(log10(S(cc))-1,log10(S(cc))+1,100);
% for ii = 1:length(splot)
% Eplot(ii) = f_predErr(f_cue2E(splot(ii),predCue{cc}));
% end
% figure; semilogx(splot,Eplot)
end
Pext{cc} = f_cue2E(S(cc),predCue{cc});
predE(:,cc) = Pext{cc}(:);
dE(cc) = f_predErr(predE(:,cc)); % mean(abs( predE(:,cc) - actual )) / (Erange-Eoffset);
end
% optimize weighting
W = zeros(size(S)); % weighting of cue
isan = not(all(isnan(predE)));
Wtmp = abs( fminsearch(@(w) f_predErr(predE(:,isan) * abs(w)./sum(abs(w))),W(isan)) ); %rms( predE * abs(w)./sum(abs(w)) - actual ),W) );
Wtmp = Wtmp/sum(Wtmp);
dE(cc+1) = f_predErr(predE(:,isan)*Wtmp);
W(isan) = Wtmp;
% Combined externalization score
dimCues = ndims(Pext{1})+1;
Pext{Ncues+1} = W(1)*Pext{1};
for cc = 2:Ncues
Pext{Ncues+1} = nansum(cat(dimCues,Pext{Ncues+1},W(cc)*Pext{cc}),dimCues);
end
% Interaction term
dE(cc+2) = f_predErr(sqrt(predE(:,InteractionIds(1)).*predE(:,InteractionIds(2))));
Pext{Ncues+2} = sqrt(Pext{InteractionIds(1)}.*Pext{InteractionIds(2)});
end
function E = cue2E(s,cue,Erange,Eoffset)
E = Erange * exp(-cue/s) + Eoffset;
end
function y = nanrms(x)
id = not(isnan(x));
y = rms(x(id));
end
Build with Bootstrap