THE AUDITORY MODELING TOOLBOX
Go to function
EXP_BARUMERLI2021 - - Experiments of Barumerli et al. (2021)
Program code:
function varargout = exp_barumerli2021(varargin)
%EXP_BARUMERLI2021 - Experiments of Barumerli et al. (2021)
%
% Usage: [] = exp_barumerli2021(flag)
%
%
% To display Fig.2 use :
%
% exp_barumerli2021('fig2');
%
% To display Fig.3 use :
%
% exp_barumerli2021('fig3');
%
% To display Fig.4 use :
%
% exp_barumerli2021('fig4');
%
% To display exp_middlebrooks1999 use :
%
% exp_barumerli2021('exp_middlebrooks1999');
%
% To display exp_macpherson2003 use :
%
% exp_barumerli2021('exp_macpherson2003');
%
% To display exp_best2005 use :
%
% exp_barumerli2021('exp_best2005');
%
% Url: http://amtoolbox.sourceforge.net/amt-0.10.0/doc/experiments/exp_barumerli2021.php
% Copyright (C) 2009-2020 Piotr Majdak and the AMT team.
% This file is part of Auditory Modeling Toolbox (AMT) version 1.0.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/>.
% See also: barumerli2021 barumerli2021_featureextraction
% AUTHOR: Roberto Barumerli
% Information Engineering Dept., University of Padova, Italy, 2021
%% ------ Check input options ---------------------------------------------
definput.import = {'amt_cache'};
definput.keyvals.MarkerSize = 6;
definput.keyvals.FontSize = 9;
definput.keyvals.repetitions = 300;
definput.keyvals.sbj = [];
definput.flags.type = {'missingflag', 'fig2', 'fig3', 'fig4', 'exp_middlebrooks1999', 'exp_macpherson2003', 'exp_best2005', 'verify_calibration'};
definput.flags.plot = {'plot', 'no_plot'};
definput.flags.plot_type = {'interp','scatter'};
definput.flags.redo = {'no_redo_fast','redo_fast'};
definput.flags.test = {'no_test','test'};
[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
%% ------ fig.2 from paper
if flags.do_fig2
data_majdak = data_majdak2010('Learn_M');
data_majdak([1:5]) = [];
calibrations = amt_cache('get', 'barumerli2021_calibration',flags.cachemode);
assert(~isempty(calibrations), 'please run calibration experiment')
if size(calibrations.sigma, 1) ~= length(data_majdak)
warning('sigma values not enough for the provided subejcts')
end
fig2 = amt_cache('get', 'barumerli2021_fig2',flags.cachemode);
% Preallocation
if isempty(fig2)
s = 1;
calibrations = calibrations.sigma(s,2); % select pge
% sofa = SOFAload(['dtf_' lower '.sofa']);
sofa = SOFAload(fullfile(SOFAdbPath,'barumerli2021',['ARI_' upper(data_majdak(s).id) '_hrtf_M_dtf 256.sofa']));
mtx = data_majdak(s).mtx;
[template, target] = barumerli2021_featureextraction(sofa, ...
'pge', 'monoaural',...
'targ_az', mtx(:, 1), ...
'targ_el', mtx(:, 2));
m_noprior = barumerli2021('template', template, ...
'target', target, ...
'num_exp', 1, ...
'sigma_lateral',calibrations.values(1), ...
'sigma_lateral_bis', calibrations.values(2), ...
'sigma_polar', calibrations.values(3), ...
'sigma_prior', [],...
'sigma_sensorimotor', []);
m_prior = barumerli2021('template', template, ...
'target', target, ...
'num_exp', 1, ...
'sigma_lateral',calibrations.values(1), ...
'sigma_lateral_bis', calibrations.values(2), ...
'sigma_polar', calibrations.values(3), ...
'sigma_prior', calibrations.values(5),...
'sigma_sensorimotor', []);
m_prior_motor = barumerli2021('template', template, ...
'target', target, ...
'num_exp', 1, ...
'sigma_lateral',calibrations.values(1), ...
'sigma_lateral_bis', calibrations.values(2), ...
'sigma_polar', calibrations.values(3), ...
'sigma_prior', calibrations.values(5),...
'sigma_sensorimotor', calibrations.values(4));
results = {m_noprior, m_prior, m_prior_motor, mtx};
amt_cache('set', 'barumerli2021_fig2', results);
end
m_noprior = fig2{1,1};
m_prior = fig2{1,2};
m_prior_motor = fig2{1,3};
mtx = fig2{1,4};
FontSize = 8;
Size = 10;
fig = figure('Units', 'points', 'Position', [1000 1000 430 150]);
subplot(1,3,1)
scatter(mtx(:, 6), mtx(:, 8), Size, 0.6*[1 1 1], 'filled');
hold on
m_noprior(~((m_noprior(:, 6) > -30) & (m_noprior(:, 6) < 210)),:)=[];
scatter(m_noprior(:, 6), m_noprior(:, 8), Size, 0.2*[1 1 1]);
grid on
ax = gca;
axis equal
set(ax, 'YTick',[-90, 0, 90, 180, 270], 'XTick',[-90, 0, 90, 180, 270], 'XLim', [-90, 270], 'YLim', [-90, 270]);
ax.XLabel.String = ' Actual Angle [deg]';
ax.TitleFontWeight = 'normal';
ax.Title.String = 'Maximum likelihood';
ax.YLabel.String = 'Estimated Angle [deg]';
subplot(1,3,2)
scatter(mtx(:, 6), mtx(:, 8), Size, 0.6*[1 1 1], 'filled');
hold on
m_prior(~((m_prior(:, 6) > -30) & (m_prior(:, 6) < 210)),:)=[];
scatter(m_prior(:, 6), m_prior(:, 8), Size, 0.2*[1 1 1]);
grid on
ax = gca;
axis equal
set(ax, 'YTick',[-90, 0, 90, 180, 270], 'XTick',[-90, 0, 90, 180, 270], 'XLim', [-90, 270], 'YLim', [-90, 270]);
ax.XLabel.String = ' Actual Angle [deg]';
ax.TitleFontWeight = 'normal';
ax.Title.String = 'w/ prior';
set(ax,'YTickLabel','')
subplot(1,3,3)
scatter(mtx(:, 6), mtx(:, 8), Size, 0.6*[1 1 1], 'filled');
hold on
m_prior_motor(~((m_prior_motor(:, 6) > -30) & (m_prior_motor(:, 6) < 210)),:)=[];
scatter(m_prior_motor(:, 6), m_prior_motor(:, 8), Size, 0.2*[1 1 1]);
grid on
ax = gca;
axis equal
set(ax, 'YTick',[-90, 0, 90, 180, 270], 'XTick',[-90, 0, 90, 180, 270], 'XLim', [-90, 270], 'YLim', [-90, 270]);
ax.XLabel.String = ' Actual Angle [deg]';
ax.TitleFontWeight = 'normal';
ax.Title.String = 'w/ sensorimotor';
set(ax,'YTickLabel','')
end
%% -------------- fig 3 paper --------------------------
if flags.do_fig3
data_majdak = data_majdak2010('Learn_M');
data_majdak([1:5]) = [];
calibrations = amt_cache('get', 'barumerli2021_calibration',flags.cachemode);
assert(~isempty(calibrations), 'please run calibration experiment')
if size(calibrations.sigma, 1) ~= length(data_majdak)
warning('sigma values not enough for the provided subejcts')
end
fig3 = amt_cache('get', 'barumerli2021_fig3',flags.cachemode);
% Preallocation
if isempty(fig3)
s = 1;
calibrations = calibrations.sigma(s,2); % select pge
% sofa = SOFAload(['dtf_' lower(data_baseline(s).id) '.sofa']);
sofa = SOFAload(fullfile(SOFAdbPath,'barumerli2021',['ARI_' upper(data_majdak(s).id) '_hrtf_M_dtf 256.sofa']));
mtx = data_majdak(s).mtx;
[template, target] = barumerli2021_featureextraction(sofa, ...
'pge', 'monoaural',...
'targ_az', mtx(:, 1), ...
'targ_el', mtx(:, 2));
m_motor = barumerli2021('template', template, ...
'target', target, ...
'num_exp', 1, ...
'sigma_lateral',calibrations.values(1), ...
'sigma_lateral_bis', calibrations.values(2), ...
'sigma_polar', calibrations.values(3), ...
'sigma_prior', calibrations.values(5),...
'sigma_sensorimotor', calibrations.values(4));
results = {m_motor, mtx};
amt_cache('set', 'barumerli2021_fig3', results);
end
m_motor = fig3{1,1};
mtx= fig3{1,2};
%% FIGURE
FontSize = 9;
Size = 10;
fig = figure('Units', 'Points', 'Position', [1e3 1e3 245 150]);
%% SCATTER
subplot(1,2,1)
scatter(mtx(:, 5), mtx(:, 7), Size, 0.6*[1 1 1],'filled'); % shadow
hold on
scatter(m_motor(:, 5), m_motor(:, 7), Size, 0.2*[1 1 1]); % shadow
grid on
ax = gca;
axis equal
set(ax, 'YTick',[-90 -45 0 45 90], 'XTick',[-90 -45 0 45 90], 'XLim', [-100, 100], 'YLim', [-100, 100], 'FontSize',FontSize)
ax.XLabel.String = 'Actual Angle [deg]';
ax.YLabel.String = 'Estimated Angle [deg]';
ax.YLabel.FontSize = FontSize;
ax.XLabel.FontSize = FontSize;
ax.Title.String = 'Lateral';
ax.TitleFontWeight = 'normal';
title('Lateral')
subplot(1,2,2)
scatter(mtx(:, 6), mtx(:, 8), Size, 0.6*[1 1 1], 'filled'); % shadow
hold on
scatter(m_motor(:, 6), m_motor(:, 8), Size, 0.2*[1 1 1]); % shadow
grid on
ax = gca;
axis equal
set(ax, 'YTick',[-90, 0, 90, 180, 270], 'XTick',[-90, 0, 90, 180, 270], 'XLim', [-90, 270], 'YLim', [-90, 270], 'FontSize',FontSize)
ax.XLabel.String = 'Actual Angle [deg]';
ax.YLabel.FontSize = FontSize;
ax.XLabel.FontSize = FontSize;
ax.Title.String = 'Polar';
ax.TitleFontWeight = 'normal';
end
%% -------------- fig 4 paper --------------------------
if flags.do_fig4
data_majdak = data_majdak2010('Learn_M');
data_majdak([1:5]) = [];
calibrations = amt_cache('get', 'barumerli2021_calibration',flags.cachemode);
assert(~isempty(calibrations), 'please run calibration experiment')
if size(calibrations.sigma, 1) ~= length(data_majdak)
warning('sigma values not enough for the provided subejcts')
end
fig4 = amt_cache('get', 'barumerli2021_fig4',flags.cachemode);
% compute results
if isempty(fig4)
for s=1:length(data_majdak)
calib = calibrations.sigma(s,2); % select pge
% sofa = SOFAload(['dtf_' lower(data_baseline(s).id) '.sofa']);
sofa = SOFAload(fullfile(SOFAdbPath,'barumerli2021',['ARI_' upper(data_majdak(s).id) '_hrtf_M_dtf 256.sofa']));
mtx = data_majdak(s).mtx;
[template, target] = barumerli2021_featureextraction(sofa, ...
'pge', 'monoaural',...
'targ_az', mtx(:, 1), ...
'targ_el', mtx(:, 2));
m = barumerli2021('template', template, ...
'target', target, ...
'num_exp', 1, ...
'sigma_lateral',calib.values(1), ...
'sigma_lateral_bis', calib.values(2), ...
'sigma_polar', calib.values(3), ...
'sigma_prior', calib.values(5),...
'sigma_sensorimotor', calib.values(4));
results{s} = m;
end
amt_cache('set', 'barumerli2021_fig4', results);
end
%% COMPUTE RMS
lat = [-90, -40, 0, 40, 90];
lat_label = [-65, -20, 20, 65];
pol = [-30 30 150 210];
pol_label = [0 90 180];
% results
rmsL = zeros(length(data_majdak), length(lat)-1);
rmsL_real = zeros(length(data_majdak), length(lat)-1);
rmsP = zeros(1, length(pol)-1);
rmsP_real = zeros(1, length(pol)-1);
for s = 1:length(data_majdak)
m = fig4{s};
mtx = data_majdak(s).mtx;
for i=1:length(lat)-1
% real
mtx_temp = mtx((mtx(:, 5) > lat(i) & (mtx(:, 5) < lat(i+1))),:);
rmsL_real(s,i) = localizationerror(mtx_temp, 'rmsL');
% simulation
m_temp = m((m(:, 5) > lat(i) & (m(:, 5) < lat(i+1))),:);
rmsL(s,i) = localizationerror(m_temp, 'rmsL');
end
mtx(abs(mtx(7,:)) <= 30,:)=[];
m(abs(m(7,:)) <= 30,:)=[];
for i=1:length(pol)-1
m_temp = mtx((mtx(:, 6) > pol(i) & (mtx(:, 6) < pol(i+1))),:);
rmsP_real(s,i) = localizationerror(m_temp, 'rmsPmedianlocal');
m_temp = m((m(:, 6) > pol(i) & (m(:, 6) < pol(i+1))),:);
rmsP(s,i) = localizationerror(m_temp, 'rmsPmedianlocal');
end
end
%% FIGURE
FontSize = 9;
Size = 36;
fig = figure('Units', 'Points', 'Position', [1e3 1e3 510 180]);
sp = 1;
for s = 1:length(data_majdak)
ax(sp) = subplot(2, length(data_majdak), sp);
ax(sp).Position(4) = ax(s).Position(4)/1.612;
% lateral
scatter(lat_label, rmsL_real(s,:), Size, 0.8*[1 1 1], 'filled')
hold on
scatter(lat_label, rmsL(s,:), Size, 0.2*[1 1 1])
set(gca, 'YLim', [5 20], 'Ytick', [5:5:20], 'Xtick', lat_label, 'Xlim', [-90, 90],'FontSize',FontSize)
title({upper(data_majdak(s).id)})
grid on
% polar
ax(length(data_majdak)+sp) = subplot(2, length(data_majdak), length(data_majdak) + sp);
scatter(pol_label, rmsP_real(s,:), Size, 0.8*[1 1 1], 'filled')
hold on
scatter(pol_label, rmsP(s,:), Size, 0.2*[1 1 1])
set(gca, 'Xtick', pol_label, 'Xlim', [-90, 270], 'Ylim', [0 60], 'Ytick', [0:20:60],'FontSize',FontSize)
grid on
sp = sp + 1;
end
set(ax([2:5, 7:10]),'YTickLabel','')
for i=6:10
ax(i).Title.String = [];
ax(i).YLabel.FontSize = FontSize;
end
ax(8).XLabel.String = ' Actual Angle [deg]';
ax(8).XLabel.FontSize = FontSize;
for i=1:5
ax(i).TitleFontWeight = 'normal';
ax(i).YLabel.FontSize = FontSize;
end
ax(1).YLabel.String = 'LE [deg]';
ax(6).YLabel.String = 'PE [deg]';
end
%% ------ verify_calibration -----------------------------------------------------------
if flags.do_verify_calibration
data_majdak = data_majdak2010('Learn_M');
data_majdak([1:5]) = [];
calibrations = amt_cache('get', 'barumerli2021_calibration',flags.cachemode);
assert(~isempty(calibrations), 'please run calibration experiment')
if size(calibrations.sigma, 1) ~= length(data_majdak)
warning('sigma values not enough for the provided subejcts')
end
verify_calibration = amt_cache('get', 'barumerli2021_verify_calibration',flags.cachemode);
% Preallocation
if isempty(verify_calibration)
verify_calibration = repmat(struct('err', struct([])),length(data_majdak), size(calibrations.combination, 1));
num_calib = size(calibrations.combination,1);
for s = 1:size(calibrations.sigma, 1) % subjects
amt_disp('\n %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% \n SUBJECT %s\n', data_majdak(s).id)
% sofa = SOFAload(['dtf_' lower(data_baseline(s).id) '.sofa']);
% sofa = SOFAload(fullfile(SOFAdbPath,'barumerli2021',['dtf_' lower(calibrations.name{s}) '.sofa']));
sofa = SOFAload(fullfile(SOFAdbPath,'barumerli2021',['ARI_' upper(calibrations.name{s}) '_hrtf_M_dtf 256.sofa']));
for c = 1:num_calib% feature space
amt_disp(sprintf('\nCOMBINATION %i', c))
amt_disp(sprintf('%s ', calibrations.combination{c,:}))
[template_par, target_par] = barumerli2021_featureextraction(sofa, ...
calibrations.combination{c,1}, ...
calibrations.combination{c,2}, ...
calibrations.combination{c,3});
sigma_l = calibrations.sigma(s,c).values(1);
sigma_l2 = calibrations.sigma(s,c).values(2);
sigma_p = calibrations.sigma(s,c).values(3);
sigma_m = calibrations.sigma(s,c).values(4);
sigma_prior = calibrations.sigma(s,c).values(5);
if calibrations.sigma(s,c).values(3) == 0
sigma_p = [];
end
m = barumerli2021('template', template_par, ...
'target', target_par, ...
'num_exp', 300, ...
'sigma_lateral', sigma_l, ...
'sigma_lateral_bis', sigma_l2, ...
'sigma_polar', sigma_p, ...
'sigma_sensorimotor', sigma_m,...
'MAP',...
'sigma_prior', sigma_prior);
verify_calibration(s,c).err = barumerli2021_metrics(m,'middle_metrics');
end
end
amt_cache('set', 'barumerli2021_verify_calibration',verify_calibration);
end
% metrics
metric_calib = zeros(size(calibrations.combination, 1),3);
for s = 1:size(calibrations.sigma, 1)
amt_disp(sprintf('\nSUBJECT %i', s))
mtx = data_majdak(s).mtx;
real = barumerli2021_metrics(mtx,'middle_metrics');
for c = 1:size(calibrations.combination, 2)
% plot
amt_disp(sprintf('\t\t\tFEATURE SPACE'))
amt_disp(sprintf('%s ', calibrations.combination{c,:}))
le_tau = abs(verify_calibration(s,c).err.rmsL - real.rmsL)/real.rmsL;
pe_tau = abs(verify_calibration(s,c).err.rmsP - real.rmsP)/real.rmsP;
qe_tau_rau = abs(rau(verify_calibration(s,c).err.querr, 1, 'PC') - rau(real.querr, 1, 'PC'))/rau(real.querr, 1, 'PC');
amt_disp(['le_tau ', num2str(le_tau),' pe_tau ', num2str(pe_tau),' querr_tau ', num2str(qe_tau_rau)]);%.2f\n', , pe_tau, qe_tau_rau)
metric_calib(c,1) = metric_calib(c,1) + verify_calibration(s,c).err.rmsL;
metric_calib(c,2) = metric_calib(c,2) + verify_calibration(s,c).err.rmsP;
metric_calib(c,3) = metric_calib(c,3) + verify_calibration(s,c).err.querr;
end
end
metric_calib = metric_calib./size(calibrations.sigma, 1);
amt_disp(sprintf('\n############\nAVERAGE OVER SUBJECTS'))
for c = 1:size(calibrations.combination, 1)
amt_disp(sprintf('\t\t\tFEATURE SPACE'))
amt_disp(sprintf('%s ', calibrations.combination{c,:}))
amt_disp(['le ', num2str(metric_calib(c,1)), ' pe ', num2str(metric_calib(c,2)), ' querr ', num2str(metric_calib(c,3))]);%.2f, pe %.2f, querr %.2f \n', metric_calib(c,1), metric_calib(c,2), metric_calib(c,3))
end
end
%% ------ middlebrooks -------------------------------------
if flags.do_exp_middlebrooks1999
data_majdak = data_majdak2010('Learn_M');
data_majdak([1:5]) = [];
exp_middlebrooks = [];
if ~flags.do_redo
exp_middlebrooks = amt_cache('get', ...
'exp_middlebrooks1999',flags.cachemode);
end
calibrations = amt_cache('get', 'barumerli2021_calibration');
if isempty(calibrations)
error('please run calibration experiment')
end
if size(calibrations.sigma, 1) ~= length(calibrations.name)
error('sigma values not enough for the provided subejcts')
end
sbj_num = length(calibrations.name);
cal_num = size(calibrations.combination, 1);
num_exp = 50;
if flags.do_redo_fast
num_exp = 5;
exp_middlebrooks = [];
end
if flags.do_test
num_exp = 1;
cal_num = 1;
sbj_num = 1;
end
if isempty(exp_middlebrooks) || flags.do_redo_fast
% preprocess templates for each user
amt_disp('Processing subjects'' templates');
for s = 1:sbj_num
amt_disp(['Pre-processing subject #' num2str(s)],'progress');
%sofa = SOFAload(['dtf_' lower(calibrations.name{s}) '.sofa']);
sofa = SOFAload(fullfile(SOFAdbPath,'barumerli2021',['ARI_' upper(calibrations.name{s}) '_hrtf_M_dtf 256.sofa']));
horpolar_coords = barumerli2021_coordinates(sofa).return_positions('horizontal-polar');
idx = find((horpolar_coords(:, 2) >= -60 ...
& horpolar_coords(:, 2) <= 79) ...
| (horpolar_coords(:, 2) >= 101 & ...
horpolar_coords(:, 2) <= 240));
for c = 1:cal_num
[template(c,s), target(c,s)] = ...
barumerli2021_featureextraction(sofa, ...
calibrations.combination{c,1}, ...
calibrations.combination{c,2}, ...
calibrations.combination{c,3});
end
end
% preallocation for results
amt_disp('Allocating memory for results');
estimations = struct('m', []);
estimations = repmat(estimations, cal_num, ...
sbj_num, sbj_num); % all vs all
for c = 1:cal_num
amt_disp(sprintf('Combination #%i', c));
for s = 1:sbj_num
amt_disp(sprintf('\tSubject #%i', s));
assert(length(calibrations.sigma(s,c).values) == 5, 'something is wrong with the calibration file')
sigma_l = calibrations.sigma(s,c).values(1);
sigma_l2 = calibrations.sigma(s,c).values(2);
sigma_p = calibrations.sigma(s,c).values(3);
sigma_m = calibrations.sigma(s,c).values(4);
sigma_prior = calibrations.sigma(s,c).values(5);
for j = 1:sbj_num
amt_disp(num2str(j));
estimations(c, s, j).m = ...
barumerli2021('template', template(c, s), 'target', target(c, j), 'num_exp', num_exp, ...
'sigma_lateral', sigma_l, ...
'sigma_lateral_bis', sigma_l2, ...
'sigma_polar', sigma_p,...
'sigma_sensorimotor', sigma_m, ...
'sigma_prior', sigma_prior);
end
end
end
% compute metrics
for c = 1:size(estimations, 1)
for i = 1:size(estimations, 2)
for j = 1:size(estimations, 3)
metrics(c, i, j) = barumerli2021_metrics(estimations(c, i, j).m, 'middle_metrics');
end
end
end
exp_middlebrooks = metrics;
if ~flags.do_redo_fast
amt_cache('set','exp_middlebrooks1999',exp_middlebrooks);
end
end
metrics_all = exp_middlebrooks;
num_calib = size(metrics_all, 1);
quants = [0,0.05,0.25,0.5,0.75,0.95,1];
% iterate over calibrations
for c=1:num_calib
metrics = squeeze(metrics_all(c,:,:));
% aggregate metrics
ns = size(metrics,1);
own = logical(eye(ns));
other = not(own);
% code similar to baumgartner2014 - fig9
le_own(c,1).quantiles = quantile([metrics(own).rmsL], quants);
lb_own(c,1).quantiles = quantile([metrics(own).accL], quants);
qe_own(c,1).quantiles = quantile([metrics(own).querr], quants);
pe_own(c,1).quantiles = quantile([metrics(own).rmsP], quants);
pb_own(c,1).quantiles = quantile([metrics(own).accP], quants);
le_own(c,1).mean = mean([metrics(own).rmsL]);
lb_own(c,1).mean = mean([metrics(own).accL]);
qe_own(c,1).mean = mean([metrics(own).querr]);
pe_own(c,1).mean = mean([metrics(own).rmsP]);
pb_own(c,1).mean = mean([metrics(own).accP]);
le_other(c,1).quantiles = quantile([metrics(other).rmsL], quants);
lb_other(c,1).quantiles = quantile([metrics(other).accL], quants);
qe_other(c,1).quantiles = quantile([metrics(other).querr], quants);
pe_other(c,1).quantiles = quantile([metrics(other).rmsP], quants);
pb_other(c,1).quantiles = quantile([metrics(other).accP], quants);
le_other(c,1).mean = mean([metrics(other).rmsL]);
lb_other(c,1).mean = mean([metrics(other).accL]);
qe_other(c,1).mean = mean([metrics(other).querr]);
pe_other(c,1).mean = mean([metrics(other).rmsP]);
pb_other(c,1).mean = mean([metrics(other).accP]);
end
% load reference data
data_middle = data_middlebrooks1999;
% load actual data and use the middlebrooks grid
for i=1:size(data_majdak, 2)
metric_majdak(i) = barumerli2021_metrics(data_majdak(i).mtx, 'middle_metrics');
end
clear data_majdak
data_majdak.le.quantiles = quantile([metric_majdak.rmsL], quants);
data_majdak.lb.quantiles = quantile([metric_majdak.accL], quants);
data_majdak.qe.quantiles = quantile([metric_majdak.querr], quants);
data_majdak.pe.quantiles = quantile([metric_majdak.rmsP], quants);
data_majdak.pb.quantiles = quantile([metric_majdak.accP], quants);
data_majdak.le.mean = mean([metric_majdak.rmsL]);
data_majdak.lb.mean = mean([metric_majdak.accL]);
data_majdak.qe.mean = mean([metric_majdak.querr]);
data_majdak.pe.mean = mean([metric_majdak.rmsP]);
data_majdak.pb.mean = mean([metric_majdak.accP]);
% baumgartner data
data_baum_temp = exp_baumgartner2014('fig9', 'no_plot');
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),quants);
data_baum.pe_own.quantiles = quantile(data_baum.pe_pool(own),quants);
data_baum.pb_own.quantiles = quantile(data_baum.pb_pool(own),quants);
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),quants);
data_baum.pe_other.quantiles = quantile(data_baum.pe_pool(other),quants);
data_baum.pb_other.quantiles = quantile(data_baum.pb_pool(other),quants);
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));
% reijniers
data_reij_temp = exp_reijniers2014('fig2_barumerli2020forum', 'no_plot');
ns = size(data_reij_temp,1);
own = logical(eye(ns));
other = not(own);
data_reij.le_own.quantiles = quantile([data_reij_temp(own).rmsL], quants);
data_reij.lb_own.quantiles = quantile([data_reij_temp(own).accL], quants);
data_reij.qe_own.quantiles = quantile([data_reij_temp(own).querr], quants);
data_reij.pe_own.quantiles = quantile([data_reij_temp(own).rmsP], quants);
data_reij.pb_own.quantiles = quantile([data_reij_temp(own).accP], quants);
data_reij.le_own.mean = mean([data_reij_temp(own).rmsL]);
data_reij.lb_own.mean = mean([data_reij_temp(own).accL]);
data_reij.qe_own.mean = mean([data_reij_temp(own).querr]);
data_reij.pe_own.mean = mean([data_reij_temp(own).rmsP]);
data_reij.pb_own.mean = mean([data_reij_temp(own).accP]);
data_reij.le_other.quantiles = quantile([data_reij_temp(other).rmsL], quants);
data_reij.lb_other.quantiles = quantile([data_reij_temp(other).accL], quants);
data_reij.qe_other.quantiles = quantile([data_reij_temp(other).querr], quants);
data_reij.pe_other.quantiles = quantile([data_reij_temp(other).rmsP], quants);
data_reij.pb_other.quantiles = quantile([data_reij_temp(other).accP], quants);
data_reij.le_other.mean = mean([data_reij_temp(other).rmsL]);
data_reij.lb_other.mean = mean([data_reij_temp(other).accL]);
data_reij.qe_other.mean = mean([data_reij_temp(other).querr]);
data_reij.pe_other.mean = mean([data_reij_temp(other).rmsP]);
data_reij.pb_other.mean = mean([data_reij_temp(other).accP]);
% plot
if flags.do_plot
calib_plot_order = [3,1,2]; %[lat, dtf, pge]
% spacing
dx = 0.11;
% multiplier for horizontal shift
middle_off = 3;
majdak_off = 2;
reij_off = 1;
baum_off = 0;
Marker = 's-';
LineColor = [[0.9290, 0.6940, 0.1250]; ...
[0.4940, 0.1840, 0.5560]; ...
[0.4660, 0.6740, 0.1880]; ...
[0.3010, 0.7450, 0.9330]; ...
[0.6350, 0.0780, 0.1840]];
data_middle.Marker = 'ko-';
data_middle.LineColor = 'k';%[1 1 1]*0.3;
data_majdak.Marker = 'b^-';
data_majdak.LineColor = 'b';%[0 0 1]*0.3;
data_baum.Marker = 'd-';
data_baum.LineColor = [0.8500 0.3250 0.0980];
data_reij.Marker = 'v-';
data_reij.LineColor = [0 0.4470 0.7410];
mFig = figure;
mFig.Units = 'points';
mFig.Position = [0, 0, 510, 300];
tile_left = 0.02;%[left bottom width height]
tile_width = 0.32;
%% SUBPLOT 1
sp = 1;
subplot(1, 3, sp);
ax=gca;
ax.OuterPosition(1) = tile_left;
ax.OuterPosition(3) = tile_width;
ax.Position(3) = 0.26;
% reference
local_middlebroxplot(ax, 1-middle_off*dx,data_middle.le_own, data_middle.Marker, kv.MarkerSize, data_middle.LineColor, data_middle.LineColor);
local_middlebroxplot(ax, 2-middle_off*dx,data_middle.le_other, data_middle.Marker, kv.MarkerSize, data_middle.LineColor, data_middle.LineColor);
% baseline
local_middlebroxplot(ax, 1-majdak_off*dx,data_majdak.le, data_majdak.Marker, kv.MarkerSize, data_majdak.LineColor, data_majdak.LineColor);
% simulation
cdist = 1;
for c=calib_plot_order
local_middlebroxplot(ax, 1+cdist*dx,le_own(c,1), Marker, kv.MarkerSize, LineColor(c,:), 'w');
local_middlebroxplot(ax, 2+cdist*dx,le_other(c,1), Marker, kv.MarkerSize, LineColor(c,:), 'w');
cdist = cdist+1;
end
% reijniers2014
local_middlebroxplot(ax, 1-reij_off*dx,data_reij.le_own, data_reij.Marker, kv.MarkerSize, data_reij.LineColor, 'w');
local_middlebroxplot(ax, 2-reij_off*dx,data_reij.le_other, data_reij.Marker, kv.MarkerSize, data_reij.LineColor, 'w');
%
ylabel('Lateral RMS Error [deg]','FontSize',kv.FontSize)
set(gca,'YLim',[0 45],'YTick', 0:10:40,'XLim',[0.5 2.5],...
'XTick',1:2,'XTickLabel',{'Own' 'Other'},'FontSize',kv.FontSize,...
'TickLength',2*get(gca,'TickLength'))
%% SUBPLOT 2
sp = sp +1;
subplot(1, 3, sp);
ax=gca;
ax.OuterPosition(1) = tile_left*2 + tile_width;
ax.OuterPosition(3) = tile_width;
ax.Position(3) = 0.26;
% reference
local_middlebroxplot(ax, 1-middle_off*dx, data_middle.pe_own, data_middle.Marker, kv.MarkerSize, data_middle.LineColor, data_middle.LineColor);
local_middlebroxplot(ax, 2-middle_off*dx, data_middle.pe_other, data_middle.Marker, kv.MarkerSize, data_middle.LineColor, data_middle.LineColor);
% baseline
local_middlebroxplot(ax, 1-majdak_off*dx, data_majdak.pe, data_majdak.Marker, kv.MarkerSize, data_majdak.LineColor, data_majdak.LineColor);
% simulations
cdist = 1;
for c=calib_plot_order
local_middlebroxplot(ax, 1+cdist*dx, pe_own(c,1), Marker, kv.MarkerSize, LineColor(c,:),'w');
local_middlebroxplot(ax, 2+cdist*dx, pe_other(c,1), Marker, kv.MarkerSize, LineColor(c,:),'w');
cdist = cdist + 1;
end
% reijniers2014
local_middlebroxplot(ax, 1-reij_off*dx, data_reij.pe_own,data_reij.Marker, kv.MarkerSize, data_reij.LineColor,'w');
local_middlebroxplot(ax, 2-reij_off*dx, data_reij.pe_other,data_reij.Marker, kv.MarkerSize, data_reij.LineColor,'w');
% baumgartner2014
local_middlebroxplot(ax, 1-baum_off*dx, data_baum.pe_own,data_baum.Marker, kv.MarkerSize, data_baum.LineColor,'w');
local_middlebroxplot(ax, 2-baum_off*dx, data_baum.pe_other,data_baum.Marker, kv.MarkerSize, data_baum.LineColor,'w');
ylabel('Local Polar RMS Error [deg]','FontSize',kv.FontSize)
set(gca,'YLim',[0 65],'YTick', 0:10:60,'XLim',[0.5 2.5],...
'XTick',1:2,'XTickLabel',{'Own' 'Other'},'FontSize',kv.FontSize,...
'TickLength',2*get(gca,'TickLength'))
%% SUBPLOT 3
sp = sp +1;
sp_ref = subplot(1, 3, sp);
ax=gca;
ax.OuterPosition(1) = tile_left*3 + tile_width*2;
ax.OuterPosition(3) = tile_width;
ax.Position(3) = 0.26;
% reference
middle = local_middlebroxplot(ax, 1-middle_off*dx,data_middle.qe_own, data_middle.Marker, kv.MarkerSize, data_middle.LineColor, data_middle.LineColor);
local_middlebroxplot(ax, 2-middle_off*dx,data_middle.qe_other, data_middle.Marker, kv.MarkerSize, data_middle.LineColor, data_middle.LineColor);
% baseline
baseline = local_middlebroxplot(ax, 1-majdak_off*dx,data_majdak.qe, data_majdak.Marker, kv.MarkerSize, data_majdak.LineColor, data_majdak.LineColor);
% simulations
cdist = 1;
for c=calib_plot_order
baru(1,c) = local_middlebroxplot(ax, 1+cdist*dx,qe_own(c,1), Marker, kv.MarkerSize, LineColor(c,:),'w');
local_middlebroxplot(ax, 2+cdist*dx,qe_other(c,1), Marker, kv.MarkerSize, LineColor(c,:),'w');
cdist = cdist + 1;
end
% reijniers2014
reij = local_middlebroxplot(ax, 1-reij_off*dx,data_reij.qe_own, data_reij.Marker,kv.MarkerSize, data_reij.LineColor,'w');
local_middlebroxplot(ax, 2-reij_off*dx,data_reij.qe_other, data_reij.Marker,kv.MarkerSize, data_reij.LineColor,'w');
% baumgartner2014
baum = local_middlebroxplot(ax, 1-baum_off*dx,data_baum.qe_own, data_baum.Marker,kv.MarkerSize, data_baum.LineColor,'w');
local_middlebroxplot(ax, 2-baum_off*dx,data_baum.qe_other, data_baum.Marker,kv.MarkerSize, data_baum.LineColor,'w');
ylabel('Quadrant Errors [%]','FontSize',kv.FontSize)
set(gca,'YLim',[-5 55],'YTick', 0:10:50,'XLim',[0.5 2.5],...
'XTick',1:2,'XTickLabel',{'Own' 'Other'},'FontSize',kv.FontSize,...
'TickLength',2*get(gca,'TickLength'))
for c=calib_plot_order
switch lower(calibrations.combination{c,2})
case 'dtf'
labels{1,c} = '$T_{DTF}^\varphi$';
case 'pge'
labels{1,c} = '$T_{PGE}^\varphi$';
end
switch lower(calibrations.combination{c,3})
% case 'reijniers'
% labels{1,c} = 'Reijniers2014';
case 'interaural_difference'
labels{1,c} = '$T_{LAT}^\varphi$';
end
end
leg = legend(sp_ref, [middle, baseline, reij, baum, baru(calib_plot_order)], horzcat({'middlebrooks1999'},{'majdak2010'},{'reijniers2014'}, {'baumgartner2014'}, labels(calib_plot_order)));
% leg = legend(sp_ref, [middle, baseline, baum, baru(calib_plot_order)], horzcat({'Middlebrooks1999'},{'Majdak2010'}, {'Baumgartner2014'}, labels(calib_plot_order)));
leg.FontSize = kv.FontSize - 2;
leg.Units = 'centimeters';
leg.Interpreter = 'latex';
leg.Units = 'normalized';
leg.Position = [0.0996373699690041 0.679833338340123 0.183774799022303 0.21349999499321];
end
end
%% ------ exp_best -------------------------------------
if flags.do_exp_best2005
exp_best = [];
if ~flags.do_redo
exp_best = amt_cache('get', ...
'exp_best2005',flags.cachemode);
end
calibrations = amt_cache('get', 'barumerli2021_calibration');
if isempty(calibrations)
error('please run calibration experiment')
end
if size(calibrations.sigma, 1) ~= length(calibrations.name)
error('sigma values not enough for the provided subejcts')
end
% Settings
num_calib = size(calibrations.combination, 1);
num_sbj = length(calibrations.name);
num_exp = 50;
num_sample = 260;
num_exp_LP = num_exp;
if flags.do_redo_fast
amt_disp('Redo fast')
exp_best = [];
num_exp = 5;
end
if flags.do_test
amt_disp('Test')
exp_best = [];
num_exp = 1;
num_sample = 1;
num_exp_LP = num_exp;
num_calib = 1;
num_sbj = 1;
end
if isempty(exp_best)
% extract directions
% retrieved from figure 3 and 4 of the original paper
dir_hor = [-160,-150,-140,-130,-120,-90,-60,-50,-40,-30,-20,0,20,30,40,50,60,90,120,130,140,149,160,180];
dir_hor(dir_hor < -90) = dir_hor(dir_hor < -90) + 360; % adjust to interval [-90, 270]
dir_lat = [-90,-70,-60,-50,-45,-30,-20,-10,0,10,20,30,45,50,60,70,90];
[x,y] = meshgrid(dir_lat,dir_hor);
%sofa_obj = SOFAload(['dtf_' lower(calibrations.name{1}) '.sofa']);
sofa = SOFAload(fullfile(SOFAdbPath,'barumerli2021',['ARI_' upper(calibrations.name{1}) '_hrtf_M_dtf 256.sofa']));
coords = barumerli2021_coordinates(sofa_obj);
coords.normalize_distance();
coords_cart = coords.return_positions('cartesian');
targ_coords = barumerli2021_coordinates([dirs, ones(size(dirs,1),1)], 'horizontal-polar');
targ_cart = targ_coords.return_positions('cartesian');
targ_cart(targ_cart(:, 3) < min(coords_cart(:,3)), :) = []; % remove lowest positions
targ_coords = Coordinates(targ_cart, 'cartesian');
targ_idx = coords.find_positions(targ_coords);
targ_idx = unique(targ_idx);
coords_sph = coords.return_positions('spherical');
targ_az = coords_sph(targ_idx,1);
targ_el = coords_sph(targ_idx,2);
targ_cart = targ_coords.return_positions('cartesian');
% figure;
% scatter3(coords.pos(:,1), coords.pos(:,2), coords.pos(:,3), 10, 'k', 'filled')
% hold on
% scatter3(targ_cart(:,1), targ_cart(:,2), targ_cart(:,3), 20, 'b')
% scatter3(coords.pos(targ_idx,1), coords.pos(targ_idx,2), coords.pos(targ_idx,3), 50, 'r')
% legend({'sofa', 'targ', 'chosen'})
%
% Load Data
% Speech Samples from Harvard Word list
speechsample = amt_cache('get','../experiments%2Fexp_baumgartner2014.m/best2005speechSamples');
speech_fs = 48e3;
speechsample(num_sample+1:end) = [];
samples_num = length(speechsample);
for s = 1:samples_num
speechsample{s} = speechsample{s} ./ max(speechsample{s});
end
% prepare directions, stimuli combinations
if flags.do_test
comb_stim = 1;
else
repetitions = 5; % see paper sect II.D
num_stim_rep = 50;
comb_stim = zeros(length(targ_idx), repetitions);
for r = 1:repetitions
idx = repmat(transpose(((r-1)*num_stim_rep+1):(r*num_stim_rep)), ceil(length(targ_idx)/num_stim_rep), 1);
comb_stim(:,r) = idx(1:size(comb_stim,1),:);
end
% randomize everything
comb_stim = comb_stim(randperm(size(comb_stim, 1)), :);
end
% 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
estimations = struct('m', []);
est_expflat = repmat(estimations, num_calib, num_sbj, 1);
est_expLP = repmat(estimations, num_calib, ...
num_sbj,length(LP));
% Simulations
amt_disp('Processing subjects HRIR');
data_majdak = data_majdak2010('Learn_M');
data_majdak([1:5]) = [];
data_actual = [];
for s = 1:num_sbj
amt_disp(['Pre-processing subject #' num2str(s)],'progress');
% load data for prior
for s_i = 1:length(data_majdak)
if strcmpi(data_majdak(s_i).id, calibrations.name{s})
data_actual(s).m = data_majdak(s_i).mtx;
end
end
%sofa = SOFAload(['dtf_' lower(calibrations.name{s}) '.sofa']);
sofa = SOFAload(fullfile(SOFAdbPath,'barumerli2021',['ARI_' upper(calibrations.name{s}) '_hrtf_M_dtf 256.sofa']));
for c = 1:num_calib
amt_disp(sprintf('\tCalibration #%i', c),'progress');
% noise stimulus
[template, target_flat] = ...
barumerli2021_featureextraction(sofa, 'targ_az', targ_az, ...
'targ_el', targ_el, ...
calibrations.combination{c,1}, ...
calibrations.combination{c,2}, ...
calibrations.combination{c,3});
assert(length(calibrations.sigma(s,c).values) == 5, 'something is wrong with the calibration file')
sigma_l = calibrations.sigma(s,c).values(1);
sigma_l2 = calibrations.sigma(s,c).values(2);
sigma_p = calibrations.sigma(s,c).values(3);
sigma_m = calibrations.sigma(s,c).values(4);
sigma_prior = calibrations.sigma(s,c).values(5);
% flat spectrum estimations
est_expflat(c, s).m = ...
barumerli2021('template', template, ...
'target', target_flat, ...
'num_exp', num_exp, ...
'sigma_lateral', sigma_l, ...
'sigma_lateral_bis', sigma_l2, ...
'sigma_polar', sigma_p, ...
'sigma_sensorimotor', sigma_m, ...
'sigma_prior', sigma_prior);
for f = 1:length(LP)
amt_disp(sprintf('\tFilter %i\n', f))
estimations = [];% zeros(size(comb_stim,1)*size(comb_stim,2)*num_exp_LP, 8);
for sa = 1:size(comb_stim, 2)
for d = 1:size(comb_stim,1)
%fprintf('%i-%i-%i\n', f, sa, d)
amt_disp([num2str(f), '-' , num2str(sa), '-', num2str(d)]);%'%i-%i-%i\n', f, sa, d)
stim = filter(LP{f},1,speechsample{comb_stim(d, sa)});
trgt = barumerli2021_featureextraction(sofa, 'target', ...
'targ_az', targ_az(d), ...
'targ_el', targ_el(d), ...
'source', 'source_ir', stim, ...
'source_fs', speech_fs, ...
calibrations.combination{c,1}, ...
calibrations.combination{c,2}, ...
calibrations.combination{c,3});
m_LP = barumerli2021('template', template, 'target', trgt, 'num_exp', num_exp_LP, ...
'sigma_lateral', sigma_l, ...
'sigma_lateral_bis', sigma_l2, ...
'sigma_polar', sigma_p, ...
'sigma_sensorimotor', sigma_m, ...
'sigma_prior', sigma_prior);
estimations = [estimations; m_LP];
end
end
est_expLP(c, s, f).m = estimations;
end
end
end
% metrics
% allocate memory for results
% aggregate over different lateral angles
ale_expflat = zeros(num_calib, num_sbj,1);
ale_expLP = zeros(num_calib, num_sbj, length(LP));
ape_expflat = zeros(num_calib, num_sbj,1);
ape_expLP = zeros(num_calib, num_sbj, length(LP));
qe_expflat = zeros(num_calib, num_sbj,1);
qe_expLP = zeros(num_calib, num_sbj, length(LP));
for c = 1:num_calib
for i = 1:num_sbj
amt_disp(['Computing metrics #' num2str(i)],'progress');
% flat spectrum estimations
m_flat = est_expflat(c, i, 1).m;
ale_expflat(c, i, 1) = localizationerror(m_flat, 'accabsL');
m_flat = m_flat((abs(m_flat(:,7))<=30),:);
ape_expflat(c, i, 1) = localizationerror(m_flat, 'accabsP');
qe_expflat(c, i, 1) = localizationerror(m_flat, 'querr');
for f = 1:length(LP)
m_LP = est_expLP(c, i, f).m;
ale_expLP(c, i, f) = localizationerror(m_LP, 'accabsL');
m_LP = m_LP((abs(m_LP(:,7))<=30),:);
ape_expLP(c, i, f) = localizationerror(m_LP, 'accabsP');
qe_expLP(c, i, f) = localizationerror(m_LP, 'querr');
end
end
end
% save cache
exp_best.ale_expflat = ale_expflat;
exp_best.ale_expLP = ale_expLP;
exp_best.ape_expflat = ape_expflat;
exp_best.ape_expLP = ape_expLP;
exp_best.qe_expflat = qe_expflat;
exp_best.qe_expLP = qe_expLP;
if ~flags.do_redo_fast
amt_cache('set','exp_best2005', exp_best);
end
end
varargout{1} = [exp_best];
for c = 1:num_calib
ale_expflat = exp_best.ale_expflat(c, :);
ale_expLP = squeeze(exp_best.ale_expLP(c, :, :, :));
ape_expflat = exp_best.ape_expflat(c, :);
ape_expLP = squeeze(exp_best.ape_expLP(c, :, :, :));
qe_expflat = exp_best.qe_expflat(c, :);
qe_expLP = squeeze(exp_best.qe_expLP(c, :, :, :));
% Pool Samples
ale_pooled = mean(ale_expLP,3);
ape_pooled = mean(ape_expLP,3);
qe_pooled = mean(qe_expLP,3);
% Confidence Intervals or standard errors
tquant_speech = 1;
tquant_noise = 1;
% barumerli model
df_speech = size(ale_expLP,1)-1;
seale_speech = std(ale_pooled,0,1)*tquant_speech/(df_speech+1);
df_noise = size(ale_expflat,1)-1;
seale_noise = std(ale_expflat)*tquant_noise/(df_noise+1);
seale = [seale_noise, seale_speech];
df_speech = size(ape_expLP,1)-1;
seape_speech = std(ape_pooled,0,1)*tquant_speech/(df_speech+1);
df_noise = size(ape_expflat,1)-1;
seape_noise = std(ape_expflat)*tquant_noise/(df_noise+1);
seape = [seape_noise, seape_speech];
% averages
% baru2021
ale(:,c) = mean([ale_expflat', ale_pooled],1,'omitnan');
ape(:,c) = mean([ape_expflat', ape_pooled],1,'omitnan');
qe(:,c) = mean([qe_expflat', qe_pooled],1, 'omitnan');
end
% load real and baum2014 data
data = data_best2005;
% reijniers2014
data_reij = exp_reijniers2014('fig3_barumerli2020forum', 'no_plot');
data_reij.ale_pooled = mean(data_reij.ale_expLP,3);
data_reij.ape_pooled = mean(data_reij.ape_expLP,3);
data_reij.qe_pooled = mean(data_reij.qe_expLP,3);
df_speech = size(data_reij.ale_expLP,1)-1;
seale_speech = std(data_reij.ale_pooled,0,1)*tquant_speech/(df_speech+1);
df_noise = size(data_reij.ale_expflat,1)-1;
seale_noise = std(data_reij.ale_expflat)*tquant_noise/(df_noise+1);
data_reij.seale = [seale_noise, seale_speech];
df_speech = size(data_reij.ape_expLP,1)-1;
seape_speech = std(data_reij.ape_pooled,0,1)*tquant_speech/(df_speech+1);
df_noise = size(data_reij.ape_expflat,1)-1;
seape_noise = std(data_reij.ape_expflat)*tquant_noise/(df_noise+1);
data_reij.seape = [seape_noise, seape_speech];
data_reij.ale = mean([data_reij.ale_expflat, data_reij.ale_pooled],1);
data_reij.ape = mean([data_reij.ape_expflat, data_reij.ape_pooled],1);
data_reij.qe = mean([data_reij.qe_expflat, data_reij.qe_pooled],1);
% baumgartner2014 model
[baum_temp, ~] = exp_baumgartner2014('fig11', 'no_plot');
% num_sbj, num_filters, num_samples
data_baum.ape = zeros([size(baum_temp{1}, 2), size(ape_expLP, 2), size(ape_expLP, 3)]);
data_baum.qe = zeros([size(baum_temp{1}, 2), size(qe_expLP, 2), size(qe_expLP, 3)]);
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};
data_baum.ape_pooled = mean(data_baum.ape,3);
data_baum.qe_pooled = mean(data_baum.qe,3);
df_speech = size(data_baum.ape,1)-1;
seape_speech = std(data_baum.ape_pooled,0,1)*tquant_speech/(df_speech+1);
df_noise = size(data_baum.ape_expflat,1)-1;
seape_noise = std(data_baum.ape_expflat)*tquant_noise/(df_noise+1);
data_baum.seape = [seape_noise, seape_speech];
% 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
num_calib = size(calibrations.combination, 1);
calib_plot_order = [3,1,2];
MarkerSize = kv.MarkerSize -1;
FontSize = kv.FontSize;
LineColor = [[0.9290, 0.6940, 0.1250]; ...
[0.4940, 0.1840, 0.5560]; ...
[0.4660, 0.6740, 0.1880]; ...
[0.3010, 0.7450, 0.9330]; ...
[0.6350, 0.0780, 0.1840]];
Marker = 's-';
data.Marker = 'ko-';
data.LineColor = [0 0 0];
data_reij.LineColor = [0 0.4470 0.7410];
data_reij.Marker = 'v--';
data_baum.Marker = 'd--';
data_baum.LineColor = [0.8500, 0.3250, 0.0980];
mFig = figure;
mFig.Units = 'points';
mFig.Position = [1e3,1e3,250,370];
tile_position = [-0.025, 0, 1.1, 0.27];
% sgtitle(sprintf('%s ', calibrations.combination{c,:}), 'Interpreter', 'none')
% % compute delta
% data.ale = data.ale - data.ale(1,1);
% ale = ale - ale(1,:);
% data.ape = data.ape - data.ape(1,1);
% ape = ape - ape(1,:);
% data_baum.ape = data_baum.ape - data_baum.ape(1,1);
% data.qe([1 2 5]) = data.qe([1 2 5]) - data.qe(1,1);
% qe = qe - qe(1,:);
% data_baum.qe = data_baum.qe-data_baum.qe(1,1);
dx = 0;
xticks = 0:size(ale_pooled,2);
sp1 = subplot(3,1,1);
plot(xticks-dx,data.ale, data.Marker,'Color', data.LineColor, 'MarkerFaceColor',data.LineColor,'MarkerSize',MarkerSize)
hold on
for c = calib_plot_order
plot(xticks-dx,ale(:,c), Marker,'Color', LineColor(c,:), 'MarkerFaceColor', 'w','MarkerSize',MarkerSize)
end
plot(xticks,data_reij.ale, data_reij.Marker,'Color', data_reij.LineColor, 'MarkerFaceColor','w','MarkerSize',MarkerSize)
% ylabel('\Delta LE [deg]','FontSize',FontSize)
% set(gca,'XLim',[-0.5 4.5],'YLim',[0 20],'YMinorTick','on')
ylabel('Lat. acc. [deg]','FontSize',FontSize)
set(gca,'XLim',[-0.5 4.5],'YLim',[0 30],'YMinorTick','on')
set(gca,'XTick',xticks,'XTickLabel',[],'FontSize',FontSize)
sp2 = subplot(3,1,2);
plot(xticks,data.ape,data.Marker,'Color', data.LineColor, 'MarkerFaceColor',data.LineColor,'MarkerSize',MarkerSize)
hold on
for c=calib_plot_order
plot(xticks-dx,ape(:,c), Marker,'Color', LineColor(c,:), 'MarkerFaceColor','w','MarkerSize',MarkerSize)
end
plot(xticks,data_reij.ape, data_reij.Marker,'Color', data_reij.LineColor, 'MarkerFaceColor', 'w','MarkerSize',MarkerSize)
plot(xticks,data_baum.ape, data_baum.Marker,'Color', data_baum.LineColor, 'MarkerFaceColor','w','MarkerSize',MarkerSize)
% ylabel('\Delta PE [deg]','FontSize',FontSize)
% set(gca,'XLim',[-0.5 4.5],'YLim',[0 50],'YMinorTick','on')
ylabel('Pol. Acc. [deg]','FontSize',FontSize)
set(gca,'XLim',[-0.5 4.5],'YLim',[0 100],'YTick',10:20:100,'YMinorTick','on')
set(gca,'XTick',xticks,'XTickLabel',[],'FontSize',FontSize)
sp3 = 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);
best = plot(xticks([1 2 5]),data.qe([1 2 5]), data.Marker,'Color', data.LineColor, 'MarkerFaceColor',data.LineColor,'MarkerSize',MarkerSize);
hold on
for c=calib_plot_order
baru(1,c) = plot(xticks-dx,qe(:,c), Marker,'Color', LineColor(c,:), 'MarkerFaceColor','w','MarkerSize',MarkerSize);
end
reij = plot(xticks-dx,data_reij.qe, data_reij.Marker,'Color', data_reij.LineColor, 'MarkerFaceColor', 'w','MarkerSize',MarkerSize);
baum = plot(xticks,data_baum.qe, data_baum.Marker,'Color', data_baum.LineColor, 'MarkerFaceColor','w','MarkerSize',MarkerSize);
% ylabel('\Delta QE [%]','FontSize',FontSize)
% set(gca,'XLim',[-0.5 4.5],'YLim',[0 50],'YMinorTick','on')
ylabel('QE [%]','FontSize',FontSize)
set(gca,'XLim',[-0.5 4.5],'YLim',[0 70],'YTick',10:20:70,'YMinorTick','on')
set(gca,'XTick',xticks,'XTickLabel',data.meta,'FontSize',FontSize)
sp3.OuterPosition = tile_position;
sp3.OuterPosition(2) = 0.01;
sp3.Position(4) = 0.24;
sp2.OuterPosition = tile_position;
sp2.OuterPosition(2) = sum(sp3.OuterPosition([2,4]));
sp2.Position(4) = 0.24;
sp1.OuterPosition = tile_position;
sp1.OuterPosition(2) = sum(sp2.OuterPosition([2,4]));
sp1.Position(4) = 0.24;
for c=calib_plot_order
switch lower(calibrations.combination{c,2})
case 'dtf'
labels{1,c} = '$T_{DTF}^\varphi$';
case 'pge'
labels{1,c} = '$T_{PGE}^\varphi$';
end
if strcmpi(calibrations.combination{c,3}, 'interaural_difference')
labels{1,c} = '$T_{LAT}^\varphi$';
end
end
leg = legend([best, reij, baum, baru(calib_plot_order)], horzcat({'best2005'}, 'reijniers2014', {'baumgartner2014'}, labels(calib_plot_order)));
leg.FontSize = FontSize - 2;
leg.Units = 'normalized';
leg.Interpreter = 'latex';
leg.Position= [0.367949860921234 0.828427722881692 0.359055268858493 0.167246848265153];
end
end
%% ------ exp_macpherson -------------------------------------
if flags.do_exp_macpherson2003
exp_macpherson = [];
if ~flags.do_redo
exp_macpherson = amt_cache('get', ...
'exp_macpherson2003',flags.cachemode);
end
calibrations = amt_cache('get', 'barumerli2021_calibration');
if isempty(calibrations)
error('please run calibration experiment')
end
if size(calibrations.sigma, 1) ~= length(calibrations.name)
error('sigma values not enough for the provided subejcts')
end
% Settings
num_exp = 50;
num_sbj = length(calibrations.name);
num_calib = size(calibrations.combination, 1);
if flags.do_redo_fast
exp_macpherson = [];
num_exp = 10;
end
if flags.do_test
exp_macpherson = [];
num_exp = 1;
num_calib = 1;
num_sbj = 1;
end
if isempty(exp_macpherson)
%sofa = SOFAload(['dtf_' lower(calibrations.name{1}) '.sofa']);
sofa = SOFAload(fullfile(SOFAdbPath,'barumerli2021',['ARI_' upper(calibrations.name{1}) '_hrtf_M_dtf 256.sofa']));
% 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 = sofa.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
sexp1 = squeeze(sexp1(:,:,1));
% 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
sexp2 = squeeze(sexp2(:,:,1));
if flags.do_test
density(2:end) = []; % ripples/oct
depth(2:end) = [];
end
% preprocess templates for each user
for i = 1:num_sbj
amt_disp(['Processing subject ', num2str(i)]);
for c = 1:num_calib
amt_disp(['Pre-processing calibration #' num2str(c)],'progress');
%sofa = SOFAload(['dtf_' lower(calibrations.name{i}) '.sofa']);
sofa = SOFAload(fullfile(SOFAdbPath,'barumerli2021',['ARI_' upper(calibrations.name{i}) '_hrtf_M_dtf 256.sofa']));
% extract directions
% filter targets' coordinates
% convert from spherical to horizontal-polar coordinates
horpolar_coords = barumerli2021_coordinates(sofa).return_positions('horizontal-polar');
% polar in [60, 120]
% lateral = 0
idx = find(((horpolar_coords(:, 2) >= -60 ...
& 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');
[template(c, i), target_flat(c, i)] = ...
barumerli2021_featureextraction(sofa, ...
'targ_az', sofa.SourcePosition(idx, 1), ...
'targ_el', sofa.SourcePosition(idx, 2), ...
calibrations.combination{c,1}, ...
calibrations.combination{c,2}, ...
calibrations.combination{c,3});
amt_disp('Densities conditions');
for j = 1:length(density)
target_exp1(c, i, j) = ...
barumerli2021_featureextraction(sofa, 'target', 'source', ...
'source_ir', squeeze(sexp1(:, j)), 'source_fs', fs, ...
'targ_az', sofa.SourcePosition(idx, 1), ...
'targ_el', sofa.SourcePosition(idx, 2), ...
calibrations.combination{c,1}, ...
calibrations.combination{c,2}, ...
calibrations.combination{c,3});
end
amt_disp('Depth conditions');
for j = 1:length(depth)
target_exp2(c, i, j) = ...
barumerli2021_featureextraction(sofa, 'target', 'source', ...
'source_ir', squeeze(sexp2(:, j)), 'source_fs', fs, ...
'targ_az', sofa.SourcePosition(idx, 1), ...
'targ_el', sofa.SourcePosition(idx, 2), ...
calibrations.combination{c,1}, ...
calibrations.combination{c,2}, ...
calibrations.combination{c,3});
end
end
end
% preallocation for results
amt_disp('Allocating memory for results');
estimations = struct('m',[]);
est_expflat = repmat(estimations, num_calib, num_sbj);
est_exp1 = repmat(estimations, num_calib, ...
num_sbj,length(density));
est_exp2 = repmat(estimations, num_calib, ...
num_sbj,length(depth));
% data for prior computation
data_majdak = data_majdak2010('Learn_M');
data_majdak([1:5]) = [];
% simulations
for i = 1:num_sbj
for c = 1:num_calib
amt_disp(sprintf('\tCalibration #%i', c),'progress');
assert(length(calibrations.sigma(i,c).values) == 5, 'something is wrong with the calibration file')
sigma_l = calibrations.sigma(i,c).values(1);
sigma_l2 = calibrations.sigma(i,c).values(2);
sigma_p = calibrations.sigma(i,c).values(3);
sigma_m = calibrations.sigma(i,c).values(4);
sigma_prior = calibrations.sigma(i,c).values(5);
amt_disp(sprintf('\tSubject #%i', i),'progress');
% flat spectrum estimations
est_expflat(c, i, 1).m = barumerli2021('template', template(c, i), 'target', target_flat(c, i), ...
'num_exp', num_exp, ...
'sigma_lateral', sigma_l, ...
'sigma_lateral_bis', sigma_l2, ...
'sigma_polar', sigma_p,...
'sigma_sensorimotor', sigma_m, ...
'sigma_prior', sigma_prior);
% rippled estimations
for j = 1:length(density)
est_exp1(c, i, j).m = barumerli2021('template', template(c, i), ...
'target', target_exp1(c, i, j), ...
'num_exp', num_exp, ...
'sigma_lateral', sigma_l, ...
'sigma_lateral_bis', sigma_l2, ...
'sigma_polar', sigma_p,...
'sigma_sensorimotor', sigma_m, ...
'sigma_prior', sigma_prior);
end
for j =1:length(depth)
est_exp2(c, i, j).m = barumerli2021('template', template(c, i), ...
'target', target_exp2(c, i, j), ...
'num_exp', num_exp, ...
'sigma_lateral', sigma_l, ...
'sigma_lateral_bis', sigma_l2, ...
'sigma_polar', sigma_p,...
'sigma_sensorimotor', sigma_m, ...
'sigma_prior', sigma_prior);
end
end
end
% metrics
% allocate memory for results
% aggregate over different lateral angles
pe_exp1 = zeros(num_calib, num_sbj, length(density));
pe_exp2 = zeros(num_calib, num_sbj, length(depth));
pe_flat = zeros(num_calib, num_sbj, 1);
for c = 1:num_calib
for i = 1:num_sbj
% compute iterative regression (see Macpherson paper and localizationerror.m)
[f,r] = localizationerror(est_expflat(c,i).m, 'sirpMacpherson2000');
pe_flat(c,i) = localizationerror(est_expflat(c,i).m, f, r, 'perMacpherson2003');
for j = 1:length(density)
pe_exp1(c, i, j) = localizationerror(est_exp1(c, i, j).m, f, r, 'perMacpherson2003');
end
for j = 1:length(depth)
pe_exp2(c, i, j) = localizationerror(est_exp2(c ,i, j).m, f, r, 'perMacpherson2003');
end
end
end
% save cache
exp_macpherson.pe_flat = pe_flat;
exp_macpherson.pe_exp1 = pe_exp1;
exp_macpherson.pe_exp2 = pe_exp2;
if ~flags.do_redo_fast
amt_cache('set','exp_macpherson2003', exp_macpherson);
end
end
% Original data:
data = data_macpherson2003;
% Reijniers2014's data
data_reij = exp_reijniers2014('fig4_barumerli2020forum', 'no_plot');
% Baumgartner2014's data
% varargout{1} = {pe_exp1,pe_exp2,pe_flat,noDCN};
data_baum_temp = exp_baumgartner2014('fig10', 'no_plot');
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
% reijniers2014
data_reij.pe_exp1 = data_reij.pe_exp1 - repmat(data_reij.pe_flat(:), 1, size(data_reij.pe_exp1, 2));
data_reij.pe_exp2 = data_reij.pe_exp2 - repmat(data_reij.pe_flat(:), 1, size(data_reij.pe_exp2, 2));
% 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
% 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]);
% reijniers data
data_reij.quart_pe_flat = quantile(data_reij.pe_flat,[.25 .50 .75]);
data_reij.quart_pe_exp1 = quantile(data_reij.pe_exp1,[.25 .50 .75]);
data_reij.quart_pe_exp2 = quantile(data_reij.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]);
for c = 1:num_calib
% simulations data
pe_flat = exp_macpherson.pe_flat(c,:);
pe_exp1 = squeeze(exp_macpherson.pe_exp1(c,:,:));
pe_exp2 = squeeze(exp_macpherson.pe_exp2(c,:,:));
% simulations
pe_exp1 = pe_exp1 - repmat(pe_flat(:), 1, size(pe_exp1, 2) );
pe_exp2 = pe_exp2 - repmat(pe_flat(:), 1, size(pe_exp2, 2) );
% simulations
quart_pe_flat(c,:) = quantile(pe_flat,[.25 .50 .75]);
quart_pe_exp1(c,:,:) = quantile(pe_exp1,[.25 .50 .75]);
quart_pe_exp2(c,:,:) = quantile(pe_exp2,[.25 .50 .75]);
end
% plot
if flags.do_plot
calib_plot_order = [3,1,2]; % lat, dtf, pge
dx = 1.05;
FontSize = kv.FontSize;
MarkerSize = kv.MarkerSize;
LineColor = [[0.9290, 0.6940, 0.1250]; ...
[0.4940, 0.1840, 0.5560]; ...
[0.4660, 0.6740, 0.1880]; ...
[0.3010, 0.7450, 0.9330]; ...
[0.6350, 0.0780, 0.1840]];
data.Marker = 'ko-';
data.LineColor = [1 1 1]*0;
data_reij.Marker = 'v-';
data_reij.LineColor = [0 0.4470 0.7410];
data_baum.Marker = 'd-';
data_baum.LineColor = [0.8500 0.3250 0.0980];
% Exp1
mFig = figure;
mFig.Units = 'points';
mFig.Position = [0 0 500 265];
% sgtitle(sprintf('%s ', calibrations.combination{c,:}), 'Interpreter', 'none')
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
for c = calib_plot_order
baru(1,c) = errorbar(data.density/dx,squeeze(quart_pe_exp1(c, 2,:)),...
squeeze(quart_pe_exp1(c,2,:) - quart_pe_exp1(c,1,:)),...
squeeze(quart_pe_exp1(c,3,:) - quart_pe_exp1(c,2,:)),...
's-','MarkerSize',MarkerSize, 'Color', LineColor(c,:),'MarkerFaceColor','w');
end
hold on
reij = errorbar(data.density*dx,data_reij.quart_pe_exp1(2,:),...
data_reij.quart_pe_exp1(2,:) - data_reij.quart_pe_exp1(1,:),...
data_reij.quart_pe_exp1(3,:) - data_reij.quart_pe_exp1(2,:),...
'v--','MarkerSize',MarkerSize, 'Color', data_reij.LineColor,'MarkerFaceColor','w');
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','w');
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
sp_ref = 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
for c = calib_plot_order
errorbar(data.depth-0.5,squeeze(quart_pe_exp2(c, 2,:)),...
squeeze(quart_pe_exp2(c,2,:) - quart_pe_exp2(c,1,:)),...
squeeze(quart_pe_exp2(c,3,:) - quart_pe_exp2(c,2,:)),...
's-','MarkerSize',MarkerSize, 'Color', LineColor(c,:),'MarkerFaceColor','w');
end
hold on
errorbar(data.depth+1,data_reij.quart_pe_exp2(2,:),...
data_reij.quart_pe_exp2(2,:) - data_reij.quart_pe_exp2(1,:),...
data_reij.quart_pe_exp2(3,:) - data_reij.quart_pe_exp2(2,:),...
'v--','MarkerSize',MarkerSize, 'Color', data_reij.LineColor,'MarkerFaceColor','w');
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','w');
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:16)
errorbar(-1,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
for c = calib_plot_order
errorbar(-0.6 + (c-1)*0.5,quart_pe_flat(c,2),...
quart_pe_flat(c,2) - quart_pe_flat(c,1),...
quart_pe_flat(c,3) - quart_pe_flat(c,2),...
's-','MarkerSize',MarkerSize, 'Color', LineColor(c,:),'MarkerFaceColor','w');
end
hold on
errorbar(1,data_reij.quart_pe_flat(2),...
data_reij.quart_pe_flat(2) - data_reij.quart_pe_flat(1),...
data_reij.quart_pe_flat(3) - data_reij.quart_pe_flat(2),...
'd--','MarkerSize',MarkerSize, 'Color', data_reij.LineColor,'MarkerFaceColor','w');
errorbar(1,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','w');
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
for c = calib_plot_order
switch lower(calibrations.combination{c,2})
case 'dtf'
labels{1,c} = '$T_{DTF}^\varphi$';
case 'pge'
labels{1,c} = '$T_{PGE}^\varphi$';
end
if strcmpi(calibrations.combination{c,3}, 'interaural_difference')
labels{1,c} = '$T_{LAT}^\varphi$';
end
end
leg = legend(sp_ref, [mach, reij, baum, baru(calib_plot_order)], horzcat({'macpherson2003'}, {'reijniers2014'}, {'baumgartner2014'}, labels(calib_plot_order)));
leg.FontSize = FontSize - 2;
leg.Units = 'normalized';
leg.Interpreter = 'latex';
leg.Orientation = 'horizontal';
leg.Position = [0.147112074200069 0.944044062817491 0.739274312186318 0.0430594892744974];
% Overall correlation between actual and predicted median values
for c=calib_plot_order
m_pe_pred = [squeeze(quart_pe_exp1(c,2,:))' squeeze(quart_pe_exp2(c,2,:))'];
m_pe_actual = [data.quart_pe_exp1(2,:) data.quart_pe_exp2(2,:)];
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(sprintf('%s: r = %0.2f', labels{1,c}, r_sqr))
end
end
end
function hg = local_middlebroxplot(ax, x, data, Marker, MarkerSize, LineColor, FaceColor)
lilen = 0.05; % length of horizontal lines
hb=[];
% Symbols
i=1; hb(i) = plot(ax, x, data.quantiles(1),'x','MarkerSize',MarkerSize, 'MarkerEdgeColor', LineColor, 'MarkerFaceColor', LineColor); % min
hold on
i=i+1; hb(i) = plot(ax, x,data.quantiles(7),'x','MarkerSize',MarkerSize, 'MarkerEdgeColor', LineColor, 'MarkerFaceColor', LineColor); % max
% Horizontal lines
i=i+1; hb(i:(i+1)) = line(ax, x+0.5*[-lilen,lilen],repmat(data.quantiles(2),2),'Color',LineColor); % lower whisker
i=i+2; hb(i:(i+1)) = line(ax, x+[-lilen,lilen],repmat(data.quantiles(3),2),'Color',LineColor); % 25% Quartile
i=i+2; hb(i:(i+1)) = line(ax, x+[-lilen,lilen],repmat(data.quantiles(4),2),'Color',LineColor); % Median
i=i+2; hb(i:(i+1)) = line(ax, x+[-lilen,lilen],repmat(data.quantiles(5),2),'Color',LineColor); % 75% Quartile
i=i+2; hb(i:(i+1)) = line(ax, x+0.5*[-lilen,lilen],repmat(data.quantiles(6),2),'Color',LineColor); % upper whisker
% Vertical lines
i=i+2; hb(i:(i+1)) = line(ax, [x,x],data.quantiles(2:3),'Color',LineColor); % connector lower whisker
i=i+2; hb(i:(i+1)) = line(ax, [x,x],data.quantiles(5:6),'Color',LineColor); % connector upper whisker
i=i+2; hb(i:(i+1)) = line(ax, [x,x]-lilen,data.quantiles([3,5]),'Color',LineColor); % left box edge
i=i+2; hb(i:(i+1)) = line(ax, [x,x]+lilen,data.quantiles([3,5]),'Color',LineColor); % left box edge
% middle value
i=i+1; hb(i) = plot(x,data.mean, Marker,'MarkerSize', MarkerSize, 'MarkerFaceColor', FaceColor, 'MarkerEdgeColor',LineColor);
% create a group to avoid issues with the legend
% https://stackoverflow.com/questions/12894652/matlab-how-to-make-a-custom-legend
hg = hggroup;
set(hb,'Parent',hg);
set(get(get(hg,'Annotation'),'LegendInformation'),...
'IconDisplayStyle','off');
function rau=rau(X,N,opt)
% RAU rationalized arcsine transform
% RAU(X,N) transforms the number of correct responses X to the
% rationalized arcsine (rau). N gives the number of repetitions.
%
% This function allows to use ANOVA statistics with percent correct scores
% because: 1) RAUs are normally distributed; 2) mean and variance of RAUs
% are not correlated with eachother; and 3) likelihood that a score will
% increase/decrease will remain constant over the range.
%
% RAU=RAU(X,N,opt) defines one of the following options:
% 'Pc' ... X is given in percent correct scores (0..100%)
% 'X' ... X is given in the number of correct responses (default)
%
% The formula are based on Sherbecoe and Studebaker,
% Int. J. of Audiology 2004; 43; 442-448
%
% See also IRAU.
% 30.8.2007, Piotr Majdak
%
if exist('opt','var')
if strcmp(upper(opt),'PC')
X=X/100*N;
elseif strcmp(upper(opt),'X')
else
error('OPT must be Pc (=percent correct) or X (=number of correct responses)');
end
end
th=asin(sqrt(X/(N+1)))+asin(sqrt((X+1)/(N+1)));
rau=146/pi*(th)-23;
Build with Bootstrap