THE AUDITORY MODELING TOOLBOX
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).
Go to function
EXP_REIJNIERS2014 - - Experiments of Reijniers et al. (2014)
Program code:
function varargout = exp_reijniers2014(varargin)
%EXP_REIJNIERS2014 - Experiments of Reijniers et al. (2014)
% Usage: [] = exp_reijniers2014(flag)
%
% EXP_REIJNIERS2014(flag) reproduces figures of the study from
% Reijniers et al. (2014). Note that HRTFs from the ARI database
% are used for Fig. 4-6, which is most probably different to the HRTFs
% used for the calculations in the paper. For Fig. 2 and 3 HRTFs from the
% Symare database are used, which are also most probably non original.
%
%
% The following flags can be specified
%
% 'fig2a' Reproduce Fig. 2 (a):
% The results for simulated trial of a locali-
% zation experiment for subject 1 [see also Fig. 3(a)]
% in front Theta = (-34,45)deg (a) is shown. The
% figures show the simulated posterior angular
% probabilities and the templates corresponding to the
% original direction Theta (solid black line) and
% estimated direction Theta_est (dotted black line).
% The simulated auditory input is shown as the solid
% red line.
%
% 'fig2b' Reproduce Fig. 2 (b):
% The results for simulated trial of a locali-
% zation experiment for subject 1 [see also Fig. 3(a)]
% in back Theta = (14,171)deg (b) is shown. The figures
% show the simulated posterior angular probabilities
% and the templates corresponding to the original
% direction Theta (solid black line) and estimated
% direction Theta_est (dotted black line). The
% simulated auditory input is shown as the solid red
% line.
%
% 'fig3' Reproduce Fig. 3:
% The simulated mean spherical error as function
% of the source position for 3 different subjects. The
% average was taken over 500 localization trials for
% each source position.
%
% 'fig4' Reproduce Fig. 4(a):
% The mean localization performance when the
% simulation results are averaged over 100 subjects.
% Superimposed arrows indicate the size and direction
% of local population response biases for different
% source positions.
%
% 'fig5' Reproduce Fig. 5:
% Sensitivity analysis of the Bayesian localization
% model. The simulated mean spherical error is shown
% as function of the source position, when each of the
% model parameters is varied separately. The average
% was taken over 500 localization trials for each
% source position and 100 subjects were pooled.
% (a) The model with input parameters as described in
% the main text. The standard deviation is doubled,
% respectively for (b) the noise on the ITD, (c) the
% internal noise and the variation on (d) the source
% spectrum. In (e), the bandwidth of the source is
% reduced to [300Hz-8kHz].
%
% 'fig6' Reproduce Fig. 6:
% The mean spherical error for different values of the
% SNR. SNR=75dB corresponds to the control situation,
% see Fig.5(a), as the magnitude in all frequency
% channels is above the system noise level.
%
% 'tab1_barumerli2020aes' Reproduce Tab. 1:
% Comparison between actual (majdak2010) and simulated,
% performances by relying on different perceptual
% metrics (middlebrooks1999b). The variable 'multiplier'
% allows to tune the internal noise.
%
% 'fig2_barumerli2020forum' Reproduce Fig.2 of Barumerli et al. (2020):
% comparison between virtual estimations
% and real data (see
% data_middlebrooks1999()) for individual
% and non-individual DTFs
%
% 'fig3_barumerli2020forum' Reproduce Fig.3 of Barumerli et al. (2020):
% estimations' evaluation with band limited spectra.
% See Best et al. 2005 and Fig 11 baumgartner2014
%
% 'fig4_barumerli2020forum' Reproduce Fig.4 of Barumerli et al. (2020):
% estimations' evaluation with rippled
% spectra. See Macpherson and Middlebrooks 2003
% and Fig 10 baumgartner2014
%
% Further, cache flags (see amt_cache) and plot flags can be specified
% (Warning: Enforcing recalculation of cached data might take several
% hours).
% 'no_plot' Do not compute plots. Flag for cluster execution.
%
% 'interp' Plot scattered interpolated data (default).
%
% 'scatter' Plot only discrete scattered data instead of inter-
% polated scattered data.
%
% 'redo_fast' Quickly recalculate data for figures without using
% cache. To speed up computation the amount of locali-
% zation trials for Fig. 3 is reduced to 20 and the
% amount of subjects for Fig. 4-6 is reduced to 12.
%
% Requirements:
% -------------
%
% 1) SOFA API v1.1 or higher from
% http://sourceforge.net/projects/sofacoustics for Matlab (e.g. in
% thirdparty/SOFA)
%
% Examples:
% ---------
%
% To display right part of Fig. 2 (a) use :
%
% exp_reijniers2014('fig2a');
%
% To display Fig. 3 and do a quick recalcultaion use :
%
% exp_reijniers2014('fig3','redo_fast');
%
% To display Fig. 6 and show a scatter plot :
%
% exp_reijniers2014('fig6','scatter','redo_fast');
%
%
% See also: reijniers2014 plot_reijniers2014 reijniers2014_preproc
%
% References:
% R. Barumerli, P. Majdak, R. Baumgartner, J. Reijniers, M. Geronazzo,
% and F. Avanzini. Predicting directional sound-localization of human
% listeners in both horizontal and vertical dimensions. In Audio
% Engineering Society Convention 148. Audio Engineering Society, 2020.
%
% R. Barumerli, P. Majdak, R. Baumgartner, M. Geronazzo, and F. Avanzini.
% Evaluation of a human sound localization model based on bayesian
% inference. In Forum Acusticum, 2020.
%
% J. Reijniers, D. Vanderleist, C. Jin, C. S., and H. Peremans. An
% ideal-observer model of human sound localization. Biological
% Cybernetics, 108:169--181, 2014.
%
%
% Url: http://amtoolbox.sourceforge.net/amt-0.10.0/doc/experiments/exp_reijniers2014.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: Jonas Reijniers
% Modified and adapted for amtoolbox by
% Roberto Barumerli and Michael Sattler,
% Acoustics Research Institute, Vienna, Austria, 2019
%% ------ Check input options ---------------------------------------------
definput.import = {'amt_cache'};
definput.keyvals.MarkerSize = 6;
definput.keyvals.FontSize = 12;
definput.flags.type = {'missingflag', 'fig2a','fig2b', ...
'fig3','fig4','fig5','fig6', ...
'tab1_barumerli2020aes', ...
'fig2_barumerli2020forum', ...
'fig3_barumerli2020forum',...
'fig4_barumerli2020forum'};
definput.flags.plot = {'plot', 'no_plot'};
definput.flags.plot_type = {'interp','scatter'};
definput.flags.redo = {'no_redo_fast','redo_fast'};
[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
%% ------ fig2a -----------------------------------------------------------
if flags.do_fig2a
if flags.do_redo_fast || flags.do_redo
num_exp = 1;
fig2a = [];
else
num_exp = 1; % definedin case of flag 'redo'
fig2a = amt_cache('get','fig2a',flags.cachemode);
end
if isempty(fig2a)
% Load SOFA file
filename='HA01.sofa';
path = ['db://reijniers2014/',filename];
SOFA_obj=SOFAload(path);
% Get index of closest directions to Theta shown in Fig. 2(a)
az = -34;
el = 45;
% Preprocessing source information for both directions
[template, target] = reijniers2014_preproc(SOFA_obj, ...
'targ_az', az, 'targ_el', el);
[fig2a.doa, fig2a.params] = ...
reijniers2014(template, target, 'num_exp', num_exp);
amt_cache('set','fig2a',fig2a);
end
if flags.do_plot
% plot left side of Fig. 2 (a)
plot_reijniers2014(fig2a.params.template_coords, ...
squeeze(fig2a.params.post_prob), ...
'target', fig2a.doa.real, ...
flags.plot_type,flags.type);
title('Fig. 2(a), posterior probability density');
% plot right side of fig 2(a)
figure('NumberTitle', 'off', 'Name', 'Fig. 2 (a) templates');
set(gcf, 'Position', [100,150, 1100, 550]);
s1 = subplot(1,3,1);
set(s1, 'position', [0.07, 0.11, 0.05, 0.8] );
plot(0, fig2a.params.T_target(fig2a.params.Tidx.itd),'k.', ...
0, fig2a.params.T_template(fig2a.params.Tidx.itd),'b.', ...
0, fig2a.params.X(fig2a.params.Tidx.itd),'r.','MarkerSize',20);
set(gca,'Color',[0.99,0.95,0.93]);
ylim([-10,10]);
xlim([-0.2,0.2]);
ax1 = gca;
title('ITD','FontSize', 14);
ylabel('itd (jnd)','FontSize', 13);
set(ax1,'xticklabel',[]);
s2 = subplot(1,3,2);
set(s2, 'position', [0.19, 0.11, 0.35, 0.8] );
plot(fig2a.params.freq_channels,...
fig2a.params.T_template(fig2a.params.est_idx, fig2a.params.Tidx.Hm),'k', ...
fig2a.params.freq_channels,fig2a.params.T_target(fig2a.params.Tidx.Hm),'k:', ...
fig2a.params.freq_channels,squeeze(fig2a.params.X(1,1,fig2a.params.Tidx.Hm)),'r', ...
fig2a.params.freq_channels,squeeze(fig2a.params.X(1,1,fig2a.params.Tidx.Hm)),'r.', ...
'MarkerSize',20,'LineWidth',1.5);
legend({'$T_{-}(\Theta)$','$T_{-}(\widehat{\Theta})$', ...
'$X_{-}(\Theta)$'},'Interpreter','latex','FontSize',13);
title('spectral difference','FontSize', 14);
ylabel('logmagnitude (dB)','FontSize',13);
xlabel('log_{10}(f) (Hz)','Interpreter','tex','FontSize',13);
ylim([-20,30]);
xlim([min(fig2a.params.freq_channels) max(fig2a.params.freq_channels)]);
set(gca,'Xscale','log')
s3 = subplot(1,3,3);
set(s3, 'position', [0.61, 0.11, 0.35, 0.8] );
plot(fig2a.params.freq_channels,...
fig2a.params.T_template(fig2a.params.est_idx, fig2a.params.Tidx.Hp),'k', ...
fig2a.params.freq_channels,fig2a.params.T_target(fig2a.params.Tidx.Hp),'k:', ...
fig2a.params.freq_channels,squeeze(fig2a.params.X(1,1,fig2a.params.Tidx.Hp)),'r', ...
fig2a.params.freq_channels,squeeze(fig2a.params.X(1,1,fig2a.params.Tidx.Hp)),'r.', ...
'MarkerSize',20,'LineWidth',1.5);
legend({'$T_{+}(\Theta)$','$T_{+}(\widehat{\Theta})$', ...
'$X_{+}(\Theta)$'},'Interpreter','latex','FontSize',13);
title('spectral sum','FontSize', 14);
ylabel('logmagnitude (dB)','FontSize',13);
xlabel('log_{10}(f) (Hz)','Interpreter','tex','FontSize',13);
ylim([-20,30]);
xlim([min(fig2a.params.freq_channels) max(fig2a.params.freq_channels)]);
set(gca,'Xscale','log');
end
end
%% ------ fig2b -----------------------------------------------------------
if flags.do_fig2b
if flags.do_redo_fast || flags.do_redo
num_exp = 1;
fig2b = [];
else
num_exp = 1;
fig2b = amt_cache('get','fig2b',flags.cachemode);
end
if isempty(fig2b)
% Load SOFA file
filename='HA01.sofa';
path = ['db://reijniers2014/',filename];
SOFA_obj=SOFAload(path);
% Get index of closest directions to Theta shown in Fig. 2(a)
az = 14;
el = 171;
% Preprocessing source information for both directions
[template, target] = reijniers2014_preproc(SOFA_obj, ...
'targ_az', az, 'targ_el', el);
% Model
[fig2b.doa, fig2b.params] = ...
reijniers2014(template, target, 'num_exp', num_exp);
amt_cache('set','fig2b',fig2b);
end
if ~flags.do_no_plot
% plot left side of Fig. 2b
plot_reijniers2014(fig2b.params.template_coords, ...
squeeze(fig2b.params.post_prob), ...
'target', fig2b.doa.real, ...
flags.plot_type,flags.type);
title('Fig. 2(b), posterior probability density');
% plot right side of fig 2(b)
figure('NumberTitle', 'off', 'Name', 'Fig. 2 (b) templates');
set(gcf, 'Position', [100,150, 1100, 550]);
s1 = subplot(1,3,1);
set(s1, 'position', [0.07, 0.11, 0.05, 0.8] );
plot(0, fig2b.params.T_target(fig2b.params.Tidx.itd),'k.', ...
0, fig2b.params.T_template(fig2b.params.Tidx.itd),'b.', ...
0, fig2b.params.X(fig2b.params.Tidx.itd),'r.','MarkerSize',20);
set(gca,'Color',[0.99,0.95,0.93]);
ylim([-10,10]);
xlim([-0.2,0.2]);
ax1 = gca;
title('ITD','FontSize', 14);
ylabel('itd (jnd)','FontSize', 13);
set(ax1,'xticklabel',[]);
s2 = subplot(1,3,2);
set(s2, 'position', [0.19, 0.11, 0.35, 0.8] );
plot(fig2b.params.freq_channels,...
fig2b.params.T_template(fig2b.params.est_idx, fig2b.params.Tidx.Hm),'k', ...
fig2b.params.freq_channels,fig2b.params.T_target(fig2b.params.Tidx.Hm),'k:', ...
fig2b.params.freq_channels,squeeze(fig2b.params.X(1,1,fig2b.params.Tidx.Hm)),'r', ...
fig2b.params.freq_channels,squeeze(fig2b.params.X(1,1,fig2b.params.Tidx.Hm)),'r.', ...
'MarkerSize',20,'LineWidth',1.5);
legend({'$T_{-}(\Theta)$','$T_{-}(\widehat{\Theta})$',...
'$X_{-}(\Theta)$'},'Interpreter','latex','FontSize',13);
title('spectral difference','FontSize', 14);
ylabel('logmagnitude (dB)','FontSize',13);
xlabel('log_{10}(f) (Hz)','Interpreter','tex','FontSize',13);
ylim([-20,30]);
xlim([min(fig2b.params.freq_channels) max(fig2b.params.freq_channels)]);
set(gca,'Xscale','log');
s3 = subplot(1,3,3);
set(s3, 'position', [0.61, 0.11, 0.35, 0.8] );
plot(fig2b.params.freq_channels,...
fig2b.params.T_template(fig2b.params.est_idx, fig2b.params.Tidx.Hp),'k', ...
fig2b.params.freq_channels,fig2b.params.T_target(fig2b.params.Tidx.Hp),'k:', ...
fig2b.params.freq_channels,squeeze(fig2b.params.X(1,1,fig2b.params.Tidx.Hp)),'r', ...
fig2b.params.freq_channels,squeeze(fig2b.params.X(1,1,fig2b.params.Tidx.Hp)),'r.', ...
'MarkerSize',20,'LineWidth',1.5);
legend({'$T_{+}(\Theta)$','$T_{+}(\widehat{\Theta})$', ...
'$X_{+}(\Theta)$'},'Interpreter','latex','FontSize',13);
title('spectral sum','FontSize', 14);
ylabel('logmagnitude (dB)','FontSize',13);
xlabel('log_{10}(f) (Hz)','Interpreter','tex','FontSize',13);
ylim([-20,30]);
xlim([min(fig2b.params.freq_channels) max(fig2b.params.freq_channels)]);
set(gca,'Xscale','log');
end
end
%% ------ fig3 ------------------------------------------------------------
if flags.do_fig3
if flags.do_redo_fast
num_exp = 20;
fig3 = [];
else
num_exp = 500;
fig3 = amt_cache('get','fig3',flags.cachemode);
end
if isempty(fig3)
filenames={'HA03.sofa'; 'HA01.sofa';'HA06.sofa'};
for i = 1:length(filenames)
path = ['db://reijniers2014/',filenames{i}];
SOFA_obj(i)=SOFAload(path);
end
% Preallocation
fig3 = struct('metrics', struct([]), 'err',struct([]), 'doa_real',struct([]));
fig3 = repmat(fig3,length(filenames),1);
amt_disp('Processing three subjects in parallel...','progress');
parfor i = 1:length(filenames)
amt_disp(['Processing subject #' num2str(i)],'progress');
% Get directions from SOFA files
[template, target] = reijniers2014_preproc(SOFA_obj(i));
doa = reijniers2014(template, target,'num_exp',num_exp);
fig3(i).err = reijniers2014_metrics(doa);
fig3(i).doa_real = doa.real;
fig3(i).metrics = reijniers2014_metrics(doa, 'middle_metrics');
end
if ~flags.do_redo_fast
amt_cache('set','fig3',fig3);
end
end
if flags.do_plot
% plot for each subject
for i = 1:length(fig3)
[~, cbar] = plot_reijniers2014(fig3(i).doa_real, ...
fig3(i).err, ...
flags.plot_type, flags.type);
title(sprintf('Fig. 3(a), Subject %i', i));
caxis([0,20]);
cbar.Label.String = 'Mean spherical error [^\circ]';
cbar.Label.FontSize = 12;
ct1=get(cbar,'TickLabels');
for k=1:numel(ct1)
ct1{k}=sprintf('%s',ct1{k});
end
set(cbar,'xticklabel',ct1);
end
end
for i = 1:length(fig3)
amt_disp(sprintf('\nSUBJECT %i', i))
amt_disp(sprintf('\t\t\tSIM'))
amt_disp(sprintf('lateral_bias [deg]:\t%0.2f', ...
mean([fig3(i).metrics.accL])))
amt_disp(sprintf('lateral_rms_error [deg]:\t%0.2f', ...
mean([fig3(i).metrics.rmsL])))
amt_disp(sprintf('elevation_bias [deg]:\t%0.2f', ...
mean([fig3(i).metrics.accE])))
amt_disp(sprintf('local_rms_polar [deg]:\t%0.2f', ...
mean([fig3(i).metrics.rmsP])))
amt_disp(sprintf('quadrant_err [%%]:\t%0.2f', ...
mean([fig3(i).metrics.querr])))
end
end
%% ------ fig4 ------------------------------------------------------------
if flags.do_fig4
if flags.do_redo_fast
num_exp = 20;
num_sub = 12;
fig4 = [];
else
num_exp = 500;
num_sub = 100;
fig4 = amt_cache('get','fig4',flags.cachemode);
end
if isempty(fig4)
offset = 13; % start at 14th file in folder
% Get names of all used hrtfs
x=amt_load('reijniers2014','hrtfnames.mat');
% Load SOFA files
for i=1:num_sub
path = ['db://reijniers2014/',x.hrtfnames{i+offset}];
SOFA_obj(i)=SOFAload(path);
end
% Preallocation
fig4 = struct('exp',struct(), ...
'err',struct([]), ...
'bias',struct([]), ...
'doa_real',struct([]));
fig4 = repmat(fig4,num_sub,1);
% Compute mean spherical error for all subjects
amt_disp(['Processing ' num2str(num_sub) ...
' subjects in parallel...'],'progress');
parfor i = 1:num_sub
amt_disp(['Processing subject #' num2str(i)],'progress');
% Get directions from SOFA files
[template, target] = reijniers2014_preproc(SOFA_obj(i));
doa = reijniers2014(template, target,'num_exp',num_exp);
[fig4(i).err, fig4(i).bias] = reijniers2014_metrics(doa);
fig4(i).doa_real = doa.real;
fig4(i).exp = reijniers2014_metrics(doa, 'middle_metrics');
end
if ~flags.do_redo_fast
amt_cache('set','fig4',fig4);
end
end
% remove ill formed hrtf
if(length(fig4(51).err) ~= length(fig4(1).err))
fig4(51) = [];
end
if flags.do_plot
% Calculate averages
mean_error = zeros(size(fig4(1).doa_real, 1), 1);
mean_bias = zeros(size(fig4(1).doa_real, 1), 3);
for k = 1:length(fig4)
mean_error = mean_error + fig4(k).err;
mean_bias = mean_bias + fig4(k).bias;
end
mean_error = abs(mean_error/length(fig4));
mean_bias = mean_bias/length(fig4);
[~, cbar] = plot_reijniers2014(fig4(1).doa_real, ... % assuming same dirs
mean_error,'bias', ...
mean_bias,flags.plot_type,flags.type);
title('Fig. 4(a), simulation');
caxis([0,35]);
cbar.Label.String = 'Mean spherical error [^\circ]';
cbar.Label.FontSize = 12;
ct1=get(cbar,'TickLabels');
for k=1:numel(ct1)
ct1{k}=sprintf('%s',ct1{k});
end
set(cbar,'xticklabel',ct1);
end
metrics = struct2cell(fig4);
metrics = cell2mat(metrics(1,:));
amt_disp(sprintf('\t\t\tSIM'))
amt_disp(sprintf('lateral_bias [dg]:\t%0.2f', ...
mean([metrics.accL], 'all')))
amt_disp(sprintf('lateral_rms_error [deg]:\t%0.2f', ...
mean([metrics.rmsL], 'all')))
amt_disp(sprintf('elevation_bias [deg]:\t%0.2f', ...
mean([metrics.accE], 'all')))
amt_disp(sprintf('local_rms_polar [deg]:\t%0.2f', ...
mean([metrics.rmsP], 'all')))
amt_disp(sprintf('quadrant_err [%%]:\t%0.2f', ...
mean([metrics.querr], 'all')))
end
%% ------ fig5 ------------------------------------------------------------
if flags.do_fig5
if flags.do_redo_fast
num_exp = 20;
num_sub = 12;
fig5 = [];
else
num_exp = 500;
num_sub = 100;
fig5 = amt_cache('get','fig5',flags.cachemode);
end
if isempty(fig5)
offset = 13; % start at 14th file in folder
% Get names of all used hrtfs
x=amt_load('reijniers2014','hrtfnames.mat');
% Load SOFA files
for i=1:num_sub
path = ['db://reijniers2014/',x.hrtfnames{i+offset}];
SOFA_obj(i)=SOFAload(path);
end
% preallocation
fig5 = struct('a',struct('err',[], 'exp', struct([])), ...
'b',struct('err',[], 'exp', struct([])), ...
'c',struct('err',[], 'exp', struct([])), ...
'd',struct('err',[], 'exp', struct([])), ...
'e',struct('err',[], 'exp', struct([])), ...
'doa_real', []);
fig5 = repmat(fig5,num_sub,1);
amt_disp(['Processing ' num2str(num_sub) ...
' subjects in parallel...'],'progress');
% Compute mean sphercial error for all subjects
parfor i = 1:num_sub
%% Preprocessing source and template information
amt_disp(['Processing subject #' num2str(i)],'progress');
[template, target] = reijniers2014_preproc(SOFA_obj(i));
% Model
% Fig. 5 (a)
amt_disp('Doing control','progress');
doa = reijniers2014(template, target, ...
'num_exp',num_exp);
fig5(i).a.err = reijniers2014_metrics(doa);
fig5(i).a.exp = reijniers2014_metrics(doa, 'middle_metrics');
% Fig. 5 (b)
amt_disp('Doing sigma_itd doubled','progress');
doa = reijniers2014(template, target, ...
'sig_itd',0.569*2, 'num_exp',num_exp);
fig5(i).b.err = reijniers2014_metrics(doa);
fig5(i).b.exp = reijniers2014_metrics(doa, 'middle_metrics');
% Fig. 5 (c)
amt_disp('Doing sigma_I doubled','progress');
doa = reijniers2014(template, target, ...
'sig_I',3.5*2,'num_exp',num_exp);
fig5(i).c.err = reijniers2014_metrics(doa);
fig5(i).c.exp = reijniers2014_metrics(doa, 'middle_metrics');
% Fig. 5 (d)
amt_disp('Doing sigma_S doubled','progress');
doa = reijniers2014(template, target, ...
'sig_S',3.5*2,'num_exp',num_exp);
fig5(i).d.err = reijniers2014_metrics(doa);
fig5(i).d.exp = reijniers2014_metrics(doa, 'middle_metrics');
% Fig. 5 (e)
amt_disp('Doing LPF','progress');
[templpf, targlpf] = reijniers2014_preproc(SOFA_obj(i), ...
'fb_low', 3e2, 'fb_high', 8e3);
doa = reijniers2014(templpf, targlpf, ...
'num_exp',num_exp);
fig5(i).e.err = reijniers2014_metrics(doa);
fig5(i).e.exp = reijniers2014_metrics(doa, 'middle_metrics');
fig5(i).doa_real = doa.real;
end
if ~flags.do_redo_fast
amt_cache('set','fig5',fig5);
end
end
graphs_lb = {'a', 'b', 'c', 'd', 'e'};
if flags.do_plot
% remove anomalous subject
fig5(51) = [];
num_sub = length(fig5);
% Calculate averages and plot
mean_error = zeros(length(fig5(1).a.err), length(graphs_lb));
titles = {'Fig. 5 (a), control', ...
'Fig. 5 (b), $2\sigma_{itd}$', ...
'Fig. 5 (c), $2\sigma_{I}$',...
'Fig. 5 (d), $2\sigma_{S}$',...
'Fig. 5 (e), LPF'};
for j =1:length(graphs_lb)
for k = 1:num_sub
mean_error(:,j) = mean_error(:,j) + fig5(k).(graphs_lb{j}).err;
end
mean_error(:,j) = abs((mean_error(:,j)/num_sub));
[~, cbar] = plot_reijniers2014(fig5(1).doa_real, ...
mean_error(:,j),flags.plot_type,flags.type);
title(titles{j},'Interpreter','Latex');
caxis([0,20]);
cbar.Label.String = 'Mean spherical error [^\circ]';
cbar.Label.FontSize = 12;
ct1=get(cbar,'TickLabels');
for k=1:numel(ct1)
ct1{k}=sprintf('%s',ct1{k});
end
set(cbar,'xticklabel',ct1);
end
end
accL = zeros(length(fig5(1).a.err), length(graphs_lb));
rmsL = zeros(length(fig5(1).a.err), length(graphs_lb));
accE = zeros(length(fig5(1).a.err), length(graphs_lb));
rmsP = zeros(length(fig5(1).a.err), length(graphs_lb));
querr = zeros(length(fig5(1).a.err), length(graphs_lb));
for j =1:length(graphs_lb)
for k = 1:num_sub
accL(k,j) = accL(k,j) + fig5(k).(graphs_lb{j}).exp.accL;
rmsL(k,j) = rmsL(k,j) + fig5(k).(graphs_lb{j}).exp.rmsL;
accE(k,j) = accE(k,j) + fig5(k).(graphs_lb{j}).exp.accE;
rmsP(k,j) = rmsP(k,j) + fig5(k).(graphs_lb{j}).exp.rmsP;
querr(k,j) = querr(k,j) + fig5(k).(graphs_lb{j}).exp.querr;
end
amt_disp(sprintf('\n EXP %s', graphs_lb{j}))
amt_disp(sprintf('\t\t\tSIM'))
amt_disp(sprintf('lateral_bias [deg]:\t%0.3f', ...
mean(accL(:,j), 'all')))
amt_disp(sprintf('lateral_rms_error [deg]:\t%0.3f', ...
mean(rmsL(:,j), 'all')))
amt_disp(sprintf('elevation_bias [deg]:\t%0.3f', ...
mean(accE(:,j), 'all')))
amt_disp(sprintf('local_rms_polar [deg]:\t%0.3f', ...
mean(rmsP(:,j), 'all')))
amt_disp(sprintf('quadrant_err [%%]:\t%0.3f', ...
mean(querr(:,j), 'all')))
end
end
%% ------ fig6 ------------------------------------------------------------
if flags.do_fig6
if flags.do_redo_fast
num_exp = 20;
num_sub = 12;
fig6 = [];
else
num_exp = 500;
num_sub = 100;
fig6 = amt_cache('get','fig6',flags.cachemode);
end
SNRs = [75, 50, 40, 30];
graphs_lb = {'a', 'b', 'c', 'd'};
if isempty(fig6)
offset = 13; % start at 14th file in folder
% Get names of all used hrtfs
x=amt_load('reijniers2014','hrtfnames.mat');
% Load SOFA files
for i=1:num_sub
path = ['db://reijniers2014/',x.hrtfnames{i+offset}];
SOFA_obj(i)=SOFAload(path);
end
% preallocation
fig6 = struct('a',struct('err',[], 'exp', struct([])), ...
'b',struct('err',[], 'exp', struct([])), ...
'c',struct('err',[], 'exp', struct([])), ...
'd',struct('err',[], 'exp', struct([])), ...
'doa_real', []);
fig6 = repmat(fig6,num_sub,1);
% Compute mean sphercial error for all subjects
amt_disp(['Processing ' num2str(num_sub) ...
' subjects in parallel...'],'progress');
parfor i = 1:num_sub
amt_disp(['Processing subject #' num2str(i)],'progress');
[template, target] = reijniers2014_preproc(SOFA_obj(i));
for j = 1:length(SNRs)
amt_disp(sprintf('Doing SNR=%idB', SNRs(j)),'progress');
doa = reijniers2014(template, target, ...
'SNR',SNRs(j),'num_exp',num_exp);
fig6(i).(graphs_lb{j}).err = reijniers2014_metrics(doa);
fig6(i).(graphs_lb{j}).exp = reijniers2014_metrics(doa, 'middle_metrics');
fig6(i).doa_real = doa.real;
end
end
if ~flags.do_redo_fast
amt_cache('set','fig6',fig6);
end
end
if flags.do_plot
num_sub = length(fig6);
% Calculate averages and plot
mean_error = zeros(length(fig6(1).doa_real), length(graphs_lb));
for j =1:length(graphs_lb)
for k = 1:num_sub
mean_error(:,j) = mean_error(:,j) + fig6(k).(graphs_lb{j}).err;
end
mean_error(:,j) = abs((mean_error(:,j)/num_sub));
[~, cbar] = plot_reijniers2014(fig6(1).doa_real, ...
mean_error(:,j),flags.plot_type,flags.type);
title(sprintf('Fig.6 (%c), SNR=%idB',graphs_lb{j}, SNRs(j)),...
'Interpreter','Latex');
caxis([0,20]);
cbar.Label.String = 'Mean spherical error [^circ]';
cbar.Label.FontSize = 12;
ct1=get(cbar,'TickLabels');
for k=1:numel(ct1)
ct1{k}=sprintf('%s',ct1{k});
end
set(cbar,'xticklabel',ct1);
end
end
accL = zeros(length(fig6(1).a.err), length(graphs_lb));
rmsL = zeros(length(fig6(1).a.err), length(graphs_lb));
accE = zeros(length(fig6(1).a.err), length(graphs_lb));
rmsP = zeros(length(fig6(1).a.err), length(graphs_lb));
querr = zeros(length(fig6(1).a.err), length(graphs_lb));
for j =1:length(graphs_lb)
for k = 1:num_sub
accL(:,j) = accL(:,j) + fig6(k).(graphs_lb{j}).exp.accL;
rmsL(:,j) = rmsL(:,j) + fig6(k).(graphs_lb{j}).exp.rmsL;
accE(:,j) = accE(:,j) + fig6(k).(graphs_lb{j}).exp.accE;
rmsP(:,j) = rmsP(:,j) + fig6(k).(graphs_lb{j}).exp.rmsP;
querr(:,j) = querr(:,j) + fig6(k).(graphs_lb{j}).exp.querr;
end
amt_disp(sprintf('\n EXP %s', graphs_lb{j}))
amt_disp(sprintf('\t\t\tSIM'))
amt_disp(sprintf('lateral_bias [deg]:\t%0.2f', ...
mean([accL], 'all')))
amt_disp(sprintf('lateral_rms_error [deg]:\t%0.2f', ...
mean([rmsL], 'all')))
amt_disp(sprintf('elevation_bias [deg]:\t%0.2f', ...
mean([accE], 'all')))
amt_disp(sprintf('local_rms_polar [deg]:\t%0.2f', ...
mean([rmsP], 'all')))
amt_disp(sprintf('quadrant_err [%%]:\t%0.2f', ...
mean([querr], 'all')))
end
end
%% ------ tab1_barumerli2020aes ----------------------------------------------
if flags.do_tab1_barumerli2020aes
data_baseline = data_majdak2010('Learn_M');
% remove subjects with no data
data_baseline(1:5) = [];
% uncertainties tuner
multiplier = 3;
amt_disp(sprintf('Uncertanties multiplied by a factor: %i\n', multiplier));
tab1_barumerli2020aes = [];
if ~flags.do_redo
tab1_barumerli2020aes = amt_cache('get','tab1_barumerli2020aes',flags.cachemode);
end
if isempty(tab1_barumerli2020aes)
amt_disp('Loading SOFA files');
parfor i = 1:length(data_baseline)
data_baseline(i).sofa = ...
SOFAload(['https://sofacoustics.org/data/database/ari/',...
'dtf_' lower(data_baseline(i).id) '.sofa']);
end
% Preallocation
tab1_barumerli2020aes = struct('doa', struct([]));
tab1_barumerli2020aes = repmat(tab1_barumerli2020aes,length(data_baseline),1);
parfor i = 1:length(data_baseline)
amt_disp(['Processing subject #' num2str(i)],'progress');
% Get directions from SOFA files
targ_az = data_baseline(i).mtx(:,1);
targ_el = data_baseline(i).mtx(:,2);
% preprocessing
[template, target] = ...
reijniers2014_preproc(data_baseline(i).sofa, ...
'targ_az', targ_az, 'targ_el', targ_el);
% model esecution
[tab1_barumerli2020aes(i).doa] = ...
reijniers2014(template, target, 'num_exp', 1, ...
'sig_itd', 0.569*multiplier, ...
'sig_I', 3.5*multiplier, ...
'sig_S', 3.5*multiplier, ...
'sig', 5*multiplier);
end
amt_cache('set','tab1_barumerli2020aes',tab1_barumerli2020aes);
end
for i = 1:length(tab1_barumerli2020aes)
% lateral_bias
exp(i).accL = reijniers2014_metrics(tab1_barumerli2020aes(i).doa, 'accL');
% lateral_rms_error
exp(i).rmsL = reijniers2014_metrics(tab1_barumerli2020aes(i).doa, 'rmsL');
% elevation_bias
exp(i).accE = reijniers2014_metrics(tab1_barumerli2020aes(i).doa, 'accE');
% local_rms_polar
exp(i).rmsP = ...
reijniers2014_metrics(tab1_barumerli2020aes(i).doa, 'rmsPmedianlocal');
% quadrant_err
exp(i).querr = ...
reijniers2014_metrics(tab1_barumerli2020aes(i).doa, 'querrMiddlebrooks');
real(i).accL = localizationerror(data_baseline(i).mtx, 'accL');
real(i).rmsL = localizationerror(data_baseline(i).mtx, 'rmsL');
real(i).accE = localizationerror(data_baseline(i).mtx, 'accE');
real(i).rmsP = ...
localizationerror(data_baseline(i).mtx, 'rmsPmedianlocal');
real(i).querr = ...
localizationerror(data_baseline(i).mtx, 'querrMiddlebrooks');
end
amt_disp(sprintf('\t\t\tSIM\t\tREAL'))
amt_disp(sprintf('lateral_bias [deg]:\t%0.2f\t\t%0.2f', ...
mean([exp.accL]), mean([real.accL])))
amt_disp(sprintf('lateral_rms_error [deg]:\t%0.2f\t\t%0.2f', ...
mean([exp.rmsL]), mean([real.rmsL])))
amt_disp(sprintf('elevation_bias [deg]:\t%0.2f\t\t%0.2f', ...
mean([exp.accE]), mean([real.accE])))
amt_disp(sprintf('local_rms_polar [deg]:\t%0.2f\t\t%0.2f', ...
mean([exp.rmsP]), mean([real.rmsP])))
amt_disp(sprintf('quadrant_err [%%]:\t%0.2f\t\t%0.2f', ...
mean([exp.querr]), mean([real.querr])))
end
%% ------ fig2_barumerli2020forum -------------------------------------
if flags.do_fig2_barumerli2020forum
fig2_barumerli2020forum = [];
if ~flags.do_redo
fig2_barumerli2020forum = amt_cache('get', ...
'fig2_barumerli2020forum',flags.cachemode);
end
if isempty(fig2_barumerli2020forum)
% load 23 DTFs from baumgartner2014
amt_disp('Loading SOFA files');
sbj_dtf = data_baumgartner2014('pool','cached');
% preprocess templates for each user
amt_disp('Processing subjects'' templates');
for i = 1:length(sbj_dtf)
amt_disp(['Pre-processing subject #' num2str(i)],'progress');
[sbj_template(i), sbj_target(i)] = reijniers2014_preproc(sbj_dtf(i).Obj);
end
% preallocation for results
amt_disp('Allocating memory for results');
estimations = struct('doa', struct([]));
estimations = repmat(estimations, ...
length(sbj_dtf),length(sbj_dtf)); % all vs all
for i = 1:length(sbj_dtf)
amt_disp(['Processing subject #' num2str(i)],'progress');
parfor j = 1:length(sbj_dtf)
amt_disp(num2str(j));
%TODO: select points in the sphere grid
estimations(i, j).doa = ...
reijniers2014(sbj_template(i), sbj_target(j), 'num_exp', 1);
end
end
% compute metrics
for i = 1:size(estimations, 1)
for j = 1:size(estimations, 2)
% lateral bias
metrics(i, j).accL = reijniers2014_metrics(estimations(i, j).doa, 'accL');
% lateral_rms_error
metrics(i, j).rmsL = reijniers2014_metrics(estimations(i, j).doa, 'rmsL');
% elevation_bias
metrics(i, j).accE = reijniers2014_metrics(estimations(i, j).doa, 'accE');
% polar bias
metrics(i, j).accP = ...
reijniers2014_metrics(estimations(i, j).doa, 'accP');
% local rms polar
metrics(i, j).rmsP = ...
reijniers2014_metrics(estimations(i, j).doa, 'rmsPmedianlocal');
% quadrant_err
metrics(i, j).querr = ...
reijniers2014_metrics(estimations(i, j).doa, 'querrMiddlebrooks');
end
end
fig2_barumerli2020forum = metrics;
amt_cache('set','fig2_barumerli2020forum',fig2_barumerli2020forum);
end
metrics = fig2_barumerli2020forum;
% aggregate metrics
ns = size(metrics,1);
own = logical(eye(ns));
other = not(own);
quants = [0,0.05,0.25,0.5,0.75,0.95,1];
% code similar to baumgartner2014 - fig9
le_own.quantiles = quantile([metrics(own).rmsL], quants);
lb_own.quantiles = quantile([metrics(own).accL], quants);
qe_own.quantiles = quantile([metrics(own).querr], quants);
pe_own.quantiles = quantile([metrics(own).rmsP], quants);
pb_own.quantiles = quantile([metrics(own).accP], quants);
le_own.mean = mean([metrics(own).rmsL]);
lb_own.mean = mean([metrics(own).accL]);
qe_own.mean = mean([metrics(own).querr]);
pe_own.mean = mean([metrics(own).rmsP]);
pb_own.mean = mean([metrics(own).accP]);
le_other.quantiles = quantile([metrics(other).rmsL], quants);
lb_other.quantiles = quantile([metrics(other).accL], quants);
qe_other.quantiles = quantile([metrics(other).querr], quants);
pe_other.quantiles = quantile([metrics(other).rmsP], quants);
pb_other.quantiles = quantile([metrics(other).accP], quants);
le_other.mean = mean([metrics(other).rmsL]);
lb_other.mean = mean([metrics(other).accL]);
qe_other.mean = mean([metrics(other).querr]);
pe_other.mean = mean([metrics(other).rmsP]);
pb_other.mean = mean([metrics(other).accP]);
data = data_middlebrooks1999;
% baumgartner data
data_baum_temp = exp_baumgartner2014('fig9', 'noplot');
data_baum.qe_pool = data_baum_temp(1).qe;
data_baum.pe_pool = data_baum_temp(1).pe;
data_baum.pb_pool = data_baum_temp(1).pb;
ns = size(data_baum.pe_pool,1);
own = eye(ns) == 1;
other = not(own);
data_baum.pb_pool = abs(data_baum.pb_pool);
data_baum.qe_own.quantiles = quantile(data_baum.qe_pool(own),[0,0.05,0.25,0.5,0.75,0.95,1]);
data_baum.pe_own.quantiles = quantile(data_baum.pe_pool(own),[0,0.05,0.25,0.5,0.75,0.95,1]);
data_baum.pb_own.quantiles = quantile(data_baum.pb_pool(own),[0,0.05,0.25,0.5,0.75,0.95,1]);
data_baum.qe_own.mean = mean(data_baum.qe_pool(own));
data_baum.pe_own.mean = mean(data_baum.pe_pool(own));
data_baum.pb_own.mean = mean(data_baum.pb_pool(own));
data_baum.qe_other.quantiles = quantile(data_baum.qe_pool(other),[0,0.05,0.25,0.5,0.75,0.95,1]);
data_baum.pe_other.quantiles = quantile(data_baum.pe_pool(other),[0,0.05,0.25,0.5,0.75,0.95,1]);
data_baum.pb_other.quantiles = quantile(data_baum.pb_pool(other),[0,0.05,0.25,0.5,0.75,0.95,1]);
data_baum.qe_other.mean = mean(data_baum.qe_pool(other));
data_baum.pe_other.mean = mean(data_baum.pe_pool(other));
data_baum.pb_other.mean = mean(data_baum.pb_pool(other));
% plot
if flags.do_plot
dx = 0.22;
dx_lat = 0.15;
Marker = 's-';
LineColor = [0 0.4470 0.7410];
data.Marker = 'ko-';
data.LineColor = [1 1 1]*0.3;
data_baum.Marker = 'd-';
data_baum.LineColor = [0.8500 0.3250 0.0980];
mFig = figure;
mFig.Units = 'centimeters';
mFig.Position = [5,5,35,10];
subplot(1, 5, 1)
middlebroxplot(1-dx_lat,data.le_own.quantiles,kv.MarkerSize, data.LineColor)
plot(1-dx_lat,data.le_own.mean,data.Marker,'MarkerSize',kv.MarkerSize,'MarkerFaceColor',data.LineColor)
middlebroxplot(1+dx_lat,le_own.quantiles,kv.MarkerSize, LineColor)
plot(1+dx_lat,le_own.mean,Marker,'MarkerSize', kv.MarkerSize, 'MarkerFaceColor', LineColor, 'MarkerEdgeColor', LineColor)
middlebroxplot(2-dx_lat,data.le_other.quantiles,kv.MarkerSize, data.LineColor)
plot(2-dx_lat,data.le_other.mean,data.Marker,'MarkerSize',kv.MarkerSize,'MarkerFaceColor',data.LineColor)
middlebroxplot(2+dx_lat,le_other.quantiles,kv.MarkerSize, LineColor)
plot(2+dx_lat,le_other.mean,Marker,'MarkerSize',kv.MarkerSize, 'MarkerFaceColor', LineColor, 'MarkerEdgeColor', LineColor)
ylabel('RMS Lateral Error [deg]','FontSize',kv.FontSize)
set(gca,'YLim',[-10 60],'XLim',[0.5 2.5],...
'XTick',1:2,'XTickLabel',{'Own' 'Other'},'FontSize',kv.FontSize,...
'TickLength',2*get(gca,'TickLength'))
subplot(1, 5, 2)
middlebroxplot(1-dx_lat,data.lb_own.quantiles,kv.MarkerSize, data.LineColor)
plot(1-dx_lat,data.lb_own.mean,data.Marker,'MarkerSize',kv.MarkerSize,'MarkerFaceColor',data.LineColor)
middlebroxplot(1+dx_lat,lb_own.quantiles,kv.MarkerSize, LineColor)
plot(1+dx_lat,lb_own.mean,Marker,'MarkerSize',kv.MarkerSize, 'MarkerFaceColor', LineColor, 'MarkerEdgeColor', LineColor)
middlebroxplot(2-dx_lat,data.lb_other.quantiles,kv.MarkerSize, data.LineColor)
plot(2-dx_lat,data.lb_other.mean,data.Marker,'MarkerSize',kv.MarkerSize,'MarkerFaceColor',data.LineColor)
middlebroxplot(2+dx_lat,lb_other.quantiles,kv.MarkerSize, LineColor)
plot(2+dx_lat,lb_other.mean,Marker,'MarkerSize',kv.MarkerSize,'MarkerFaceColor', LineColor, 'MarkerEdgeColor', LineColor)
ylabel('Magnitude Lateral Bias [deg]','FontSize',kv.FontSize)
set(gca,'YLim',[-10 60],'XLim',[0.5 2.5],...
'XTick',1:2,'XTickLabel',{'Own' 'Other'},'FontSize',kv.FontSize,...
'TickLength',2*get(gca,'TickLength'))
subplot(1, 5, 3)
middlebroxplot(1-dx,data.qe_own.quantiles,kv.MarkerSize, data.LineColor)
midd = plot(1-dx,data.qe_own.mean,data.Marker,'MarkerSize',kv.MarkerSize,'MarkerFaceColor',data.LineColor);
middlebroxplot(1,qe_own.quantiles,kv.MarkerSize, LineColor)
reij = plot(1,qe_own.mean,Marker,'MarkerSize',kv.MarkerSize,'MarkerFaceColor', LineColor, 'MarkerEdgeColor', LineColor);
middlebroxplot(1+dx,data_baum.qe_own.quantiles,kv.MarkerSize, data_baum.LineColor)
baum = plot(1+dx,data_baum.qe_own.mean,data_baum.Marker,'MarkerSize',kv.MarkerSize,'MarkerFaceColor', data_baum.LineColor, 'MarkerEdgeColor', data_baum.LineColor);
middlebroxplot(2-dx,data.qe_other.quantiles,kv.MarkerSize, data.LineColor)
plot(2-dx,data.qe_other.mean,data.Marker,'MarkerSize',kv.MarkerSize,'MarkerFaceColor',data.LineColor)
middlebroxplot(2,qe_other.quantiles,kv.MarkerSize, LineColor)
plot(2,qe_other.mean,Marker,'MarkerSize',kv.MarkerSize,'MarkerFaceColor', LineColor, 'MarkerEdgeColor', LineColor)
middlebroxplot(2+dx,data_baum.qe_other.quantiles,kv.MarkerSize, data_baum.LineColor)
plot(2+dx,data_baum.qe_other.mean,data_baum.Marker,'MarkerSize',kv.MarkerSize,'MarkerFaceColor', data_baum.LineColor, 'MarkerEdgeColor', data_baum.LineColor)
ylabel('Quadrant Errors (%)','FontSize',kv.FontSize)
set(gca,'YLim',[0 50],'XLim',[0.5 2.5],...
'XTick',1:2,'XTickLabel',{'Own' 'Other'},'FontSize',kv.FontSize,...
'TickLength',2*get(gca,'TickLength'))
leg = legend([midd, reij, baum], {'Actual', 'SA', 'SP'}, 'Location', 'none');
leg.Units = 'normalized';
leg.Position = [0.4745,0.831536390309064,0.097524379329044,0.159029645257883];
subplot(1, 5, 4)
middlebroxplot(1-dx,data.pe_own.quantiles,kv.MarkerSize, data.LineColor)
plot(1-dx,data.pe_own.mean,data.Marker,'MarkerSize',kv.MarkerSize,'MarkerFaceColor',data.LineColor)
middlebroxplot(1,pe_own.quantiles,kv.MarkerSize, LineColor)
plot(1,pe_own.mean,Marker,'MarkerSize',kv.MarkerSize,'MarkerFaceColor',LineColor, 'MarkerEdgeColor', LineColor)
middlebroxplot(1+dx,data_baum.pe_own.quantiles,kv.MarkerSize, data_baum.LineColor)
plot(1+dx,data_baum.pe_own.mean,data_baum.Marker,'MarkerSize',kv.MarkerSize,'MarkerFaceColor',data_baum.LineColor, 'MarkerEdgeColor', data_baum.LineColor)
middlebroxplot(2-dx,data.pe_other.quantiles,kv.MarkerSize, data.LineColor)
plot(2-dx,data.pe_other.mean,data.Marker,'MarkerSize',kv.MarkerSize,'MarkerFaceColor',data.LineColor)
middlebroxplot(2,pe_other.quantiles,kv.MarkerSize, LineColor)
plot(2,pe_other.mean,Marker,'MarkerSize',kv.MarkerSize,'MarkerFaceColor',LineColor, 'MarkerEdgeColor', LineColor)
middlebroxplot(2+dx,data_baum.pe_other.quantiles,kv.MarkerSize, data_baum.LineColor)
plot(2+dx,data_baum.pe_other.mean,data_baum.Marker,'MarkerSize',kv.MarkerSize,'MarkerFaceColor',data_baum.LineColor, 'MarkerEdgeColor', data_baum.LineColor)
ylabel('Local Polar RMS Error (deg)','FontSize',kv.FontSize)
set(gca,'YLim',[-10 60],'XLim',[0.5 2.5],...
'XTick',1:2,'XTickLabel',{'Own' 'Other'},'FontSize',kv.FontSize,...
'TickLength',2*get(gca,'TickLength'))
subplot(1, 5, 5)
middlebroxplot(1-dx,data.pb_own.quantiles,kv.MarkerSize, data.LineColor)
plot(1-dx,data.pb_own.mean,data.Marker,'MarkerSize',kv.MarkerSize,'MarkerFaceColor',data.LineColor)
middlebroxplot(1,pb_own.quantiles,kv.MarkerSize, LineColor)
plot(1,pb_own.mean,Marker,'MarkerSize',kv.MarkerSize,'MarkerFaceColor',LineColor, 'MarkerEdgeColor', LineColor)
middlebroxplot(1+dx,data_baum.pb_own.quantiles,kv.MarkerSize, data_baum.LineColor)
plot(1+dx,data_baum.pb_own.mean,data_baum.Marker,'MarkerSize',kv.MarkerSize,'MarkerFaceColor',data_baum.LineColor, 'MarkerEdgeColor', data_baum.LineColor)
middlebroxplot(2-dx,data.pb_other.quantiles,kv.MarkerSize, data.LineColor)
plot(2-dx,data.pb_other.mean,data.Marker,'MarkerSize',kv.MarkerSize,'MarkerFaceColor',data.LineColor)
middlebroxplot(2,pb_other.quantiles,kv.MarkerSize, LineColor)
plot(2,pb_other.mean,Marker,'MarkerSize',kv.MarkerSize,'MarkerFaceColor',LineColor, 'MarkerEdgeColor', LineColor)
middlebroxplot(2+dx,data_baum.pb_other.quantiles,kv.MarkerSize, data_baum.LineColor)
plot(2+dx,data_baum.pb_other.mean,data_baum.Marker,'MarkerSize',kv.MarkerSize,'MarkerFaceColor', data_baum.LineColor, 'MarkerEdgeColor', data_baum.LineColor)
ylabel('Magnitude of Elevation Bias (deg)','FontSize',kv.FontSize)
set(gca,'YLim',[-10 60],'XLim',[0.5 2.5],...
'XTick',1:2,'XTickLabel',{'Own' 'Other'},'FontSize',kv.FontSize,...
'TickLength',2*get(gca,'TickLength'))
end
end
%% ------ fig3_barumerli2020forum -------------------------------------
if flags.do_fig3_barumerli2020forum
fig3_barumerli2020forum = [];
if ~flags.do_redo
fig3_barumerli2020forum = amt_cache('get', ...
'fig3_barumerli2020forum',flags.cachemode);
end
if isempty(fig3_barumerli2020forum)
% Settings
num_exp = 1;
% Load Data
% Speech Samples from Harvard Word list
speechsample = amt_cache('get','../experiments%2Fexp_baumgartner2014.m/best2005speechSamples');
samples_num = length(speechsample);
% FIR Low-pass filters at 8kHz
% Brick-wall (aka sinc-filter): fir1(200,1/3) -> -60 dB
x=amt_load('baumgartner2014','highfreqatten_filters.mat');
LP{1} = [1 zeros(1,100)];
LP{2} = x.fir20db;
LP{3} = x.fir40db;
LP{4} = x.fir60db;
% Model Data
sbj_dtf = data_baumgartner2014('pool', 'cached');
estimations = struct('doa', struct([]));
est_expflat = repmat(estimations, length(sbj_dtf), 1);
est_expLP = repmat(estimations, ...
length(sbj_dtf),length(LP),length(speechsample));
% Simulations
amt_disp('Processing subjects HRIR');
for i = 1:length(sbj_dtf)
amt_disp(['Pre-processing subject #' num2str(i)],'progress');
% noise stimulus
[sbj_template(i), sbj_target_flat(i)] = ...
reijniers2014_preproc(sbj_dtf(i).Obj);
% flat spectrum estimations
est_expflat(i, 1).doa = ...
reijniers2014(sbj_template(i), sbj_target_flat(i), 'num_exp', num_exp);
amt_disp('Computing localization');
for f = 1:length(LP)
amt_disp(sprintf('Filter %i\n', f))
parfor s = 1:samples_num
stim = filter(LP{f},1,speechsample{s});
[~, trgt] = ...
reijniers2014_preproc(sbj_dtf(i).Obj, ...
'source_ir', stim);
est_expLP(i, f, s).doa = ...
reijniers2014(sbj_template(i), trgt, 'num_exp', num_exp);
end
end
end
% metrics
% allocate memory for results
% aggregate over different lateral angles
ale_expflat = zeros(length(sbj_dtf),1);
ale_expLP = zeros(length(sbj_dtf),length(LP), samples_num);
ape_expflat = zeros(length(sbj_dtf),1);
ape_expLP = zeros(length(sbj_dtf),length(LP), samples_num);
qe_expflat = zeros(length(sbj_dtf),1);
qe_expLP = zeros(length(sbj_dtf),length(LP), samples_num);
for i = 1:length(sbj_dtf)
amt_disp(['Computing metrics #' num2str(i)],'progress');
% flat spectrum estimations
ale_expflat(i,1) = reijniers2014_metrics(est_expflat(i, 1).doa, 'accabsL');
ape_expflat(i,1) = reijniers2014_metrics(est_expflat(i, 1).doa, 'accabsP');
qe_expflat(i,1) = reijniers2014_metrics(est_expflat(i, 1).doa, 'querr');
for f = 1:length(LP)
for s = 1:samples_num
ale_expLP(i, f, s) = reijniers2014_metrics(est_expLP(i, f, s).doa, 'accabsL');
ape_expLP(i, f, s) = reijniers2014_metrics(est_expLP(i, f, s).doa, 'accabsP');
qe_expLP(i, f, s) = reijniers2014_metrics(est_expLP(i, f, s).doa, 'querr');
end
end
end
% save cache
fig3_barumerli2020forum.ale_expflat = ale_expflat;
fig3_barumerli2020forum.ale_expLP = ale_expLP;
fig3_barumerli2020forum.ape_expflat = ape_expflat;
fig3_barumerli2020forum.ape_expLP = ape_expLP;
fig3_barumerli2020forum.qe_expflat = qe_expflat;
fig3_barumerli2020forum.qe_expLP = qe_expLP;
amt_cache('set','fig3_barumerli2020forum', fig3_barumerli2020forum);
end
varargout{1} = [fig3_barumerli2020forum];
ale_expflat = fig3_barumerli2020forum.ale_expflat;
ale_expLP = fig3_barumerli2020forum.ale_expLP;
ape_expflat = fig3_barumerli2020forum.ape_expflat;
ape_expLP = fig3_barumerli2020forum.ape_expLP;
qe_expflat = fig3_barumerli2020forum.qe_expflat;
qe_expLP = fig3_barumerli2020forum.qe_expLP;
% load real and baum2014 data
data = data_best2005;
% DCN enabled
[baum_temp, ~] = exp_baumgartner2014('fig11', 'noplot');
data_baum.ape = zeros(size(ape_expLP));
data_baum.qe = zeros(size(qe_expLP));
for i = 1:size(ape_expLP,3)
data_baum.ape(:,:,i) = transpose(baum_temp{1}(:,:,i));
data_baum.qe(:,:,i) = transpose(baum_temp{2}(:,:,i));
end
data_baum.ape_expflat = baum_temp{3};
data_baum.qe_expflat = baum_temp{4};
% Pool Samples
ale_pooled = mean(ale_expLP,3);
ape_pooled = mean(ape_expLP,3);
qe_pooled = mean(qe_expLP,3);
data_baum.ape_pooled = mean(data_baum.ape,3);
data_baum.qe_pooled = mean(data_baum.qe,3);
% Confidence Intervals or standard errors
% reijniers model
df_speech = size(ale_expLP,1)-1;
tquant_speech = 1;
seale_speech = std(ale_pooled,0,1)*tquant_speech/(df_speech+1);
df_noise = size(ale_expflat,1)-1;
tquant_noise = 1;
seale_noise = std(ale_expflat)*tquant_noise/(df_noise+1);
seale = [seale_noise, seale_speech];
df_speech = size(ape_expLP,1)-1;
tquant_speech = 1;
seape_speech = std(ape_pooled,0,1)*tquant_speech/(df_speech+1);
df_noise = size(ape_expflat,1)-1;
tquant_noise = 1;
seape_noise = std(ape_expflat)*tquant_noise/(df_noise+1);
seape = [seape_noise, seape_speech];
% baumgartner2014 model
df_speech = size(data_baum.ape,1)-1;
tquant_speech = 1;
seape_speech = std(data_baum.ape_pooled,0,1)*tquant_speech/(df_speech+1);
df_noise = size(data_baum.ape_expflat,1)-1;
tquant_noise = 1;
seape_noise = std(data_baum.ape_expflat)*tquant_noise/(df_noise+1);
data_baum.seape = [seape_noise, seape_speech];
% averages
% reij2014
ale = mean([ale_expflat, ale_pooled],1);
ape = mean([ape_expflat, ape_pooled],1);
qe = mean([qe_expflat, qe_pooled],1);
% baum2014
data_baum.ape = mean([data_baum.ape_expflat', data_baum.ape_pooled],1);
data_baum.qe = mean([data_baum.qe_expflat', data_baum.qe_pooled],1);
if flags.do_plot
MarkerSize = kv.MarkerSize;
FontSize = kv.FontSize;
LineColor = [0 0.4470 0.7410];
Marker = 's-';
data.Marker = 'o-';
data.LineColor = [1 1 1]*0.3;
data_baum.Marker = 'd--';
data_baum.LineColor = [0.8500 0.3250 0.0980];
mFig = figure;
mFig.Units = 'centimeters';
mFig.Position = [5,5,13.5,15];
dx = 0;
xticks = 0:size(ale_pooled,2);
subplot(3,1,1)
plot(xticks-dx,data.ale, data.Marker,'Color', data.LineColor, 'MarkerFaceColor',data.LineColor,'MarkerSize',MarkerSize)
hold on
plot(xticks-dx,ale, Marker,'Color', LineColor, 'MarkerFaceColor',LineColor,'MarkerSize',MarkerSize)
ylabel('Lateral Error (deg)','FontSize',FontSize)
set(gca,'XTick',xticks,'XTickLabel',[],'FontSize',FontSize)
set(gca,'XLim',[-0.5 4.5],'YLim',[0 75],'YMinorTick','on')
subplot(3,1,2)
plot(xticks,data.ape,data.Marker,'Color', data.LineColor, 'MarkerFaceColor',data.LineColor,'MarkerSize',MarkerSize)
hold on
plot(xticks-dx,ape, Marker,'Color', LineColor, 'MarkerFaceColor',LineColor,'MarkerSize',MarkerSize)
plot(xticks,data_baum.ape, data_baum.Marker,'Color', data_baum.LineColor, 'MarkerFaceColor',data_baum.LineColor,'MarkerSize',MarkerSize)
ylabel('Polar Error (deg)','FontSize',FontSize)
set(gca,'XTick',xticks,'XTickLabel',[],'FontSize',FontSize)
set(gca,'XLim',[-0.5 4.5],'YLim',[0 75],'YMinorTick','on')
subplot(3,1,3)
best = plot(xticks([1 2 5]),data.qe([1 2 5]), data.Marker,'Color', data.LineColor, 'MarkerFaceColor',data.LineColor,'MarkerSize',MarkerSize);
hold on
reij = plot(xticks-dx,qe, Marker,'Color', LineColor, 'MarkerFaceColor',LineColor,'MarkerSize',MarkerSize);
baum = plot(xticks,data_baum.qe, data_baum.Marker,'Color', data_baum.LineColor, 'MarkerFaceColor',data_baum.LineColor,'MarkerSize',MarkerSize);
ylabel('Quadrant Err. (%)','FontSize',FontSize)
set(gca,'XTick',xticks,'XTickLabel',data.meta,'FontSize',FontSize,...
'XLim',[-0.5 4.5],'YLim',[-3 54],'YMinorTick','on')
leg = legend([best, reij, baum], {'Actual', 'SA', 'SP'});
leg.FontSize = FontSize - 1;
leg.Units = 'centimeters';
leg.Position = [9.281,12.762,3.466,1.561];
end
end
if flags.do_fig4_barumerli2020forum
fig4_barumerli2020forum = [];
if ~flags.do_redo
fig4_barumerli2020forum = amt_cache('get', ...
'fig4_barumerli2020forum',flags.cachemode);
end
if isempty(fig4_barumerli2020forum)
% load 23 DTFs from baumgartner2014
amt_disp('Loading SOFA files');
sbj_dtf = data_baumgartner2014('pool','global');
num_exp = 5;
% generate stimulus
% copyed from exp_baumgartner2014/do_fig10
density = [0.25, 0.5, 0.75, 1, 1.5, 2, 3, 4, 6, 8]; % ripples/oct
depth = 10:10:40; % ripple depth (peak-to-trough) in dB
% 250-ms bursts, 20-ms raised-cosine fade in/out, flat from 0.6-16kHz
fs = sbj_dtf(1).Obj.Data.SamplingRate;
flow = 1e3; % lower corner frequency of ripple modification in Hz
fhigh = 16e3; % upper corner frequency of ripple modification in Hz
Nf = 2^10; % # Frequency bins
f = 0:fs/2/Nf:fs/2; % frequency bins
id600 = find(f<=600,1,'last'); % index of 600 Hz (lower corner frequency of stimulus energy)
idlow = find(f<=flow,1,'last'); % index of flow (ripples)
idhigh = find(f>=fhigh,1,'first'); % index of fhigh (ripples)
N600low = idlow - id600 +1; % # bins without ripple modification
Nlowhigh = idhigh - idlow +1; % # bins with ripple modification %
O = log2(f(idlow:idhigh)/1e3); % freq. trafo. to achieve equal ripple density in log. freq. scale
% Raised-cosine "(i.e., cos^2)" ramp 1/8 octave wide
fup = f(idlow)*2^(1/8); % upper corner frequency of ramp upwards
idup = find(f<=fup,1,'last');
Nup = idup-idlow+1;
rampup = cos(-pi/2:pi/2/(Nup-1):0).^2;
fdown = f(idhigh)*2^(-1/8); % lower corner frequency of ramp downwards
iddown = find(f>=fdown,1,'first');
Ndown = idhigh-iddown+1;
rampdown = cos(0:pi/2/(Ndown-1):pi/2).^2;
ramp = [rampup ones(1,Nlowhigh-Nup-Ndown) rampdown];
ramp = [-inf*ones(1,id600-1) zeros(1,N600low) ramp -inf*ones(1,Nf - idhigh)];
% Ripples of Experiment I
Sexp1 = zeros(Nf+1,length(density),2); % 3rd dim: 1:0-phase 2:pi-phase
Sexp1(idlow:idhigh,:,1) = (40/2* sin(2*pi*density'*O+ 0))'; % depth: 40dB, 0-phase
Sexp1(idlow:idhigh,:,2) = (40/2* sin(2*pi*density'*O+pi))'; % depth: 40dB, pi-phase
Sexp1 = repmat(ramp',[1,length(density),2]) .* Sexp1;
Sexp1 = [Sexp1;Sexp1(Nf:-1:2,:,:)];
Sexp1(isnan(Sexp1)) = -100;
sexp1 = ifftreal(10.^(Sexp1/20),2*Nf);
sexp1 = circshift(sexp1,Nf); % IR corresponding to ripple modification
% Ripples of Experiment II
Sexp2 = zeros(Nf+1,length(depth),2); % 3rd dim: 1:0-phase 2:pi-phase
Sexp2(idlow:idhigh,:,1) = (depth(:)/2*sin(2*pi*1*O+ 0))'; % density: 1 ripple/oct, 0-phase
Sexp2(idlow:idhigh,:,2) = (depth(:)/2*sin(2*pi*1*O+pi))'; % density: 1 ripple/oct, pi-phase
Sexp2 = repmat(ramp',[1,length(depth),2]) .* Sexp2;
Sexp2 = [Sexp2;Sexp2(Nf-1:-1:2,:,:)];
Sexp2(isnan(Sexp2)) = -100;
sexp2 = ifftreal(10.^(Sexp2/20),2*Nf);
sexp2 = circshift(sexp2,Nf); % IR corresponding to ripple modification
% preprocess templates for each user
amt_disp('Processing subjects HRIR');
for i = 1:length(sbj_dtf)
% extract directions
% filter targets' coordinates
% convert from spherical to horizontal-polar coordinates
horpolar_coords = SOFAconvertCoordinates(...
sbj_dtf(i).Obj.SourcePosition, ...
'spherical', 'horizontal-polar');
% polar in [60, 120]
% lateral = 0
idx = find(((horpolar_coords(:, 2) >= 300 ...
| horpolar_coords(:, 2) <= 60) ...
| (horpolar_coords(:, 2) >= 120 & ...
horpolar_coords(:, 2) <= 240)) ...
& (horpolar_coords(:, 1) <= 30 & horpolar_coords(:, 1) >= -30));
amt_disp(['Pre-processing subject #' num2str(i)],'progress');
[sbj_template(i), sbj_target_flat(i)] = ...
reijniers2014_preproc(sbj_dtf(i).Obj, ...
'targ_az', sbj_dtf(i).Obj.SourcePosition(idx, 1), ...
'targ_el', sbj_dtf(i).Obj.SourcePosition(idx, 2));
amt_disp('Densities conditions');
parfor j = 1:length(density)
[~, sbj_target_exp1(i, j)] = ...
reijniers2014_preproc(sbj_dtf(i).Obj, ...
'source_ir', squeeze(sexp1(:, j, :)), ...
'targ_az', sbj_dtf(i).Obj.SourcePosition(idx, 1), ...
'targ_el', sbj_dtf(i).Obj.SourcePosition(idx, 2));
end
amt_disp('Depth conditions');
parfor j = 1:length(depth)
[~, sbj_target_exp2(i, j)] = ...
reijniers2014_preproc(sbj_dtf(i).Obj, ...
'source_ir', squeeze(sexp2(:, j, :)), ...
'targ_az', sbj_dtf(i).Obj.SourcePosition(idx, 1), ...
'targ_el', sbj_dtf(i).Obj.SourcePosition(idx, 2));
end
end
% preallocation for results
amt_disp('Allocating memory for results');
estimations = struct('doa', struct([]));
est_expflat = repmat(estimations, length(sbj_dtf));
est_exp1 = repmat(estimations, ...
length(sbj_dtf),length(density));
est_exp2 = repmat(estimations, ...
length(sbj_dtf),length(depth));
% simulations
for i = 1:length(sbj_dtf)
amt_disp(['Processing subject #' num2str(i)],'progress');
% flat spectrum estimations
est_expflat(i, 1).doa = ...
reijniers2014(sbj_template(i), sbj_target_flat(i), 'num_exp', num_exp);
% rippled estimations
parfor j = 1:length(density)
est_exp1(i, j).doa = ...
reijniers2014(sbj_template(i), sbj_target_exp1(i, j), 'num_exp', num_exp);
end
parfor j =1:length(depth)
est_exp2(i, j).doa = ...
reijniers2014(sbj_template(i), sbj_target_exp2(i, j), 'num_exp', num_exp);
end
end
% metrics
% allocate memory for results
% aggregate over different lateral angles
pe_exp1 = zeros(length(sbj_dtf),length(density));
pe_exp2 = zeros(length(sbj_dtf),length(depth));
pe_flat = zeros(length(sbj_dtf), 1);
for i = 1:length(sbj_dtf)
% compute regression (see paper)
[f,r] = reijniers2014_metrics(est_expflat(i).doa,'sirpMacpherson2000');
pe_flat(i) = reijniers2014_metrics(est_expflat(i).doa,f,r,'perMacpherson2003');
for j = 1:size(est_exp1, 2)
pe_exp1(i, j) = reijniers2014_metrics(est_exp1(i, j).doa, f, r, 'perMacpherson2003');
end
for j = 1:size(est_exp2, 2)
pe_exp2(i, j) = reijniers2014_metrics(est_exp2(i, j).doa, f, r, 'perMacpherson2003');
end
end
% save cache
fig4_barumerli2020forum.pe_flat = pe_flat;
fig4_barumerli2020forum.pe_exp1 = pe_exp1;
fig4_barumerli2020forum.pe_exp2 = pe_exp2;
amt_cache('set','fig4_barumerli2020forum', fig4_barumerli2020forum);
end
% simulations data
pe_flat = fig4_barumerli2020forum.pe_flat;
pe_exp1 = fig4_barumerli2020forum.pe_exp1;
pe_exp2 = fig4_barumerli2020forum.pe_exp2;
% Original data:
data = data_macpherson2003;
% Baumgartner2014's data
% varargout{1} = {pe_exp1,pe_exp2,pe_flat,noDCN};
data_baum_temp = exp_baumgartner2014('fig10', 'noplot');
data_baum.pe_exp1 = data_baum_temp{1,1};
data_baum.pe_exp2 = data_baum_temp{1,2};
data_baum.pe_flat = data_baum_temp{1,3};
% Phase condition handling
% average across the phase condition
% real data
data.pe_exp1 = mean(data.pe_exp1,3);
data.pe_exp2 = mean(data.pe_exp2,3);
% baumgartner data
data_baum.pe_exp1 = mean(data_baum.pe_exp1,3);
data_baum.pe_exp2 = mean(data_baum.pe_exp2,3);
idphase = 1;
% Increase
% simulations
pe_exp1 = pe_exp1 - pe_flat(:);
pe_exp2 = pe_exp2 - pe_flat(:);
% baumgartner data
data_baum.pe_exp1 = data_baum.pe_exp1 - repmat(data_baum.pe_flat(:),1,size(data_baum.pe_exp1,2));
data_baum.pe_exp2 = data_baum.pe_exp2 - repmat(data_baum.pe_flat(:),1,size(data_baum.pe_exp2,2));
% Statistics
% simulations
quart_pe_flat = quantile(pe_flat,[.25 .50 .75]);
quart_pe_exp1 = quantile(pe_exp1,[.25 .50 .75]);
quart_pe_exp2 = quantile(pe_exp2,[.25 .50 .75]);
% real data
data.quart_pe_flat = quantile(data.pe_flat,[.25 .50 .75]);
data.quart_pe_exp1 = quantile(data.pe_exp1,[.25 .50 .75]);
data.quart_pe_exp2 = quantile(data.pe_exp2,[.25 .50 .75]);
% baumgartner data
data_baum.quart_pe_flat = quantile(data_baum.pe_flat,[.25 .50 .75]);
data_baum.quart_pe_exp1 = quantile(data_baum.pe_exp1,[.25 .50 .75]);
data_baum.quart_pe_exp2 = quantile(data_baum.pe_exp2,[.25 .50 .75]);
% plot
if flags.do_plot
dx = 1.05;
FontSize = kv.FontSize;
MarkerSize = kv.MarkerSize;
LineColor = [0 0.4470 0.7410];
data.Marker = 'ko-';
data.LineColor = [1 1 1]*0.3;
data_baum.Marker = 'd-';
data_baum.LineColor = [0.8500 0.3250 0.0980];
% Exp1
mFig = figure;
mFig.Units = 'centimeters';
mFig.Position = [5,5,13.5,12];
subplot(2,8,1:8)
mach = errorbar(data.density,data.quart_pe_exp1(2,:,idphase),...
data.quart_pe_exp1(2,:,idphase) - data.quart_pe_exp1(1,:,idphase),...
data.quart_pe_exp1(3,:,idphase) - data.quart_pe_exp1(2,:,idphase),...
'o-','MarkerSize',MarkerSize, 'Color', data.LineColor, ...
'MarkerFaceColor', data.LineColor);
hold on
reij = errorbar(data.density/dx,quart_pe_exp1(2,:,idphase),...
quart_pe_exp1(2,:,idphase) - quart_pe_exp1(1,:,idphase),...
quart_pe_exp1(3,:,idphase) - quart_pe_exp1(2,:,idphase),...
's-','MarkerSize',MarkerSize, 'Color', LineColor,'MarkerFaceColor',LineColor);
hold on
baum = errorbar(data.density*dx,data_baum.quart_pe_exp1(2,:,idphase),...
data_baum.quart_pe_exp1(2,:,idphase) - data_baum.quart_pe_exp1(1,:,idphase),...
data_baum.quart_pe_exp1(3,:,idphase) - data_baum.quart_pe_exp1(2,:,idphase),...
'd--','MarkerSize',MarkerSize, 'Color', data_baum.LineColor,'MarkerFaceColor',data_baum.LineColor);
set(gca,'XScale','log','YMinorTick','on')
set(gca,'XLim',[0.25/1.2 8*1.2],'XTick',data.density,'YLim',[-16 59],'FontSize',FontSize)
xlabel('Ripple Density (ripples/octave)','FontSize',FontSize)
ylabel({'Increase in';'Polar Error Rate (%)'},'FontSize',FontSize)
% Exp2
subplot(2,8,9:13)
errorbar(data.depth,data.quart_pe_exp2(2,:,idphase),...
data.quart_pe_exp2(2,:,idphase) - data.quart_pe_exp2(1,:,idphase),...
data.quart_pe_exp2(3,:,idphase) - data.quart_pe_exp2(2,:,idphase),...
'o-','MarkerSize',MarkerSize, 'Color', data.LineColor, ...
'MarkerFaceColor', data.LineColor);
hold on
errorbar(data.depth-1,quart_pe_exp2(2,:,idphase),...
quart_pe_exp2(2,:,idphase) - quart_pe_exp2(1,:,idphase),...
quart_pe_exp2(3,:,idphase) - quart_pe_exp2(2,:,idphase),...
's-','MarkerSize',MarkerSize, 'Color', LineColor,'MarkerFaceColor',LineColor);
hold on
errorbar(data.depth+1,data_baum.quart_pe_exp2(2,:,idphase),...
data_baum.quart_pe_exp2(2,:,idphase) - data_baum.quart_pe_exp2(1,:,idphase),...
data_baum.quart_pe_exp2(3,:,idphase) - data_baum.quart_pe_exp2(2,:,idphase),...
'd--','MarkerSize',MarkerSize, 'Color', data_baum.LineColor,'MarkerFaceColor',data_baum.LineColor);
set(gca,'XLim',[data.depth(1)-5 data.depth(end)+5],'XTick',data.depth,...
'YLim',[-16 59],'YMinorTick','on','FontSize',FontSize)
xlabel('Ripple Depth (dB)','FontSize',FontSize)
ylabel({'Increase in';'Polar Error Rate (%)'},'FontSize',FontSize)
ytick = get(gca,'YTick');
ticklength = get(gca,'TickLength');
% Baseline
subplot(2,8,14:15)
errorbar(-0.5,data.quart_pe_flat(2),...
data.quart_pe_flat(2) - data.quart_pe_flat(1),...
data.quart_pe_flat(3) - data.quart_pe_flat(2),...
'o-','MarkerSize',MarkerSize, 'Color', data.LineColor, ...
'MarkerFaceColor', data.LineColor);
hold on
errorbar(0,quart_pe_flat(2),...
quart_pe_flat(2) - quart_pe_flat(1),...
quart_pe_flat(3) - quart_pe_flat(2),...
's-','MarkerSize',MarkerSize, 'Color', LineColor,'MarkerFaceColor',LineColor);
hold on
errorbar(0.5,data_baum.quart_pe_flat(2),...
data_baum.quart_pe_flat(2) - data_baum.quart_pe_flat(1),...
data_baum.quart_pe_flat(3) - data_baum.quart_pe_flat(2),...
'd--','MarkerSize',MarkerSize, 'Color', data_baum.LineColor,'MarkerFaceColor',data_baum.LineColor);
set(gca,'XLim',[-3 3],'XTick',0,'XTickLabel',{'Baseline'},...
'YLim',[-15 59],'YTick',ytick,'TickLength',3*ticklength,...
'FontSize',FontSize,'YAxisLocation','right')
xlabel(' ','FontSize',FontSize)
ylabel({'Polar Error Rate (%)'},'FontSize',FontSize)
%legend
leg = legend([mach, reij, baum], {'Actual', 'SA', 'SP'});
leg.FontSize = FontSize - 1;
leg.Units = 'centimeters';
leg.Position = [9.281,10,3.466,1.561];
% Overall correlation between actual and predicted median values
m_pe_pred = [quart_pe_exp1(2,:,idphase) quart_pe_exp2(2,:,idphase)];
m_pe_actual = [data.quart_pe_exp1(2,:,idphase) data.quart_pe_exp2(2,:,idphase)];
r = corrcoef(m_pe_pred,m_pe_actual);
r_sqr = r(2);
amt_disp('Correlation between actual and predicted median values (15 conditions):')
amt_disp(['w/ PSGE: r = ' num2str(r_sqr,'%0.2f')])
end
end
function middlebroxplot(x,quantiles,MarkerSize,LineColor)
lilen = 0.1; % length of horizontal lines
% Symbols
plot(x,quantiles(1),'x','MarkerSize',MarkerSize, 'MarkerEdgeColor', LineColor) % min
hold on
plot(x,quantiles(7),'x','MarkerSize',MarkerSize, 'MarkerEdgeColor', LineColor) % max
% Horizontal lines
line(x+0.5*[-lilen,lilen],repmat(quantiles(2),2),'Color',LineColor) % lower whisker
line(x+[-lilen,lilen],repmat(quantiles(3),2),'Color',LineColor) % 25% Quartile
line(x+[-lilen,lilen],repmat(quantiles(4),2),'Color',LineColor) % Median
line(x+[-lilen,lilen],repmat(quantiles(5),2),'Color',LineColor) % 75% Quartile
line(x+0.5*[-lilen,lilen],repmat(quantiles(6),2),'Color',LineColor) % upper whisker
% Vertical lines
line([x,x],quantiles(2:3),'Color',LineColor) % connector lower whisker
line([x,x],quantiles(5:6),'Color',LineColor) % connector upper whisker
line([x,x]-lilen,quantiles([3,5]),'Color',LineColor) % left box edge
line([x,x]+lilen,quantiles([3,5]),'Color',LineColor) % left box edge
Build with Bootstrap