% SH interpolation example. This script performs HRTF interpolation using
% two different methods:
% 1. Trunc (SH order truncation): baseline SH interpolation without
% additional preprocessing.
% 2. TA (time-alignment): HRTFs are time-aligned before interpolation,
% which improves the performance at low orders.
%
% For each method, we evaluate spatial orders ranging from 1 (lowest usable
% order) to 40 (very high order).
%
% The following plots are shown:
% 1. Interpolation errors (magnitude and phase) vs the target HRTF.
% 2. Interaural time differences (ITD) and level differences (ILD) and
% their error vs the target HRTF.
% 3. Loudness maps and perceptual spectral difference (PSD).
% 4. Performance per order according to various binaural models:
% - Reijniers2014: localisation performance (lateral and polar)
% - Baumgartner2020: externalisation rates
% - Jelfs2011: speech reception in presence of maskers
%
% WARNING: this script can take several minutes to run and will generate a
% few large files (~400MB in total)
%
% EXTERNAL DEPENDENCIES:
% Auditory Modeling Toolbox (amtoolbox.sourceforge.net)
%
% AUTHOR: Isaac Engel - isaac.engel(at)imperial.ac.uk
% March 2021

%% Clear variables
clear

%% Parameters
N_vec = [1,5,10:10:40]; % spatial orders to test
Nref = 44; % reference (very high) spatial order
r = 0.0875; % head radius = 8.75 cm
c = 343; % speed of sound
itdJND = 20; % ITD JND for plots, according to Klockgether 2016
ildJND = 0.6; % ILD JND for plots, according to Klockgether 2016

% HRTF file
hrirname = '../hrtfs/FABIAN_HRIR_measured_HATO_0.sofa'; % change as needed

% Working directory
workdir = '../processed_data'; % change as needed
[~,hrtfname] = fileparts(hrirname);
workdir = [workdir,'/',hrtfname];
if ~isfolder(workdir)
fprintf('Making folder %s...\n',workdir)
mkdir(workdir);
end
if ~isfolder([workdir,'/hnm'])
fprintf('Making folder %s...\n',[workdir,'/hnm'])
mkdir([workdir,'/hnm'])
end
if ~isfolder([workdir,'/results'])
fprintf('Making folder %s...\n',[workdir,'/results'])
mkdir([workdir,'/results'])
end

%% Load HRTF
SOFA_obj = SOFAload(hrirname); % load HRTF in SOFA format
[h,fs,az,el] = sofa2hrtf(SOFA_obj); % get HRTF data

% Zero-pad the HRIRs to increase frequency resolution
taps = 2048;
h = [h;zeros(taps-size(h,1),size(h,2),size(h,3))];

H = ffth(h); % to frequency domain
nfreqs = size(H,1);
f = linspace(0,fs/2,nfreqs).';

%% Get high-order Hnm
filename = sprintf('%s/hnm/ref.mat',workdir);
if isfile(filename) % load if it already exists
fprintf('Found %s. Loading...\n',filename)
load(filename,'hnm_ref','Hnm_TA_ref')
else
fprintf('Making reference hnm (also time-aligned Hnm, for Fig. 1)...\n');
% Standard
hnm_ref = toSH(h,Nref,'az',az,'el',el,'fs',fs,'mode','trunc');
% Time-aligned
hnm_TA_ref = toSH(h,Nref,'az',az,'el',el,'fs',fs,'mode','TA');
Hnm_TA_ref = ffth(hnm_TA_ref);
% Save
fprintf('Saving %s...\n',filename)
save(filename,'hnm_ref','Hnm_TA_ref','fs')
end
Hnm_ref = ffth(hnm_ref);

%% Defining test conditions
ncond = 2;
% 1. "Trunc": baseline method, with no preprocessing such as time-alignment
test_conditions{1}.name = 'Trunc';
test_conditions{1}.preproc = 'trunc';
% 2. "TA": time-align the HRTF prior to getting the SH coefficients
test_conditions{2}.name = 'TA';
test_conditions{2}.preproc = 'TA';

%% Define a few direction subsets
% 1. Nearest neighbours to 110-point Lebedev grid
gridTestReq = sofia_lebedev(110,0);
[~,indLeb] = getGridSubset([az,el],gridTestReq,0);
indLeb = find(indLeb); % get indices
azLeb = az(indLeb);
elLeb = el(indLeb);
wLeb = gridTestReq(:,3);
% 2. Median plane (lat = 0 +/- 1 deg)
[lat,~]=sph2hor(az*180/pi,90-el*180/pi); % lat/pol coordinates
indMP = find(abs(lat)<1);
azMP = az(indMP);
elMP = el(indMP);
% 3. Horizontal plane (el = 90 +/- 1 deg)
indHP = find(abs(el-pi/2)<(pi/180));
azHP = az(indHP);
% elHP = el(indHP);
indFront = find(abs(el-pi/2)<(pi/180) & abs(az)<(pi/180)); % front

%% Generate and save hnm for each test condition
fprintf('Generating hnms...\n')
for N=N_vec % iterate through spatial orders
fprintf('Processing order %d...\n',N);
for i=1:ncond
name = test_conditions{i}.name;
preproc = test_conditions{i}.preproc;
filename = sprintf('%s/hnm/ord%0.2d_%s.mat',workdir,N,name);
if isfile(filename)
fprintf('\tFound %s. Skipping...\n',filename)
else
fprintf('\tGenerating hnm...\n');
[hnm,varOut] = toSH(h,N,'az',az,'el',el,'fs',fs,'mode',preproc);
fprintf('\tSaving %s...\n',filename)
save(filename,'hnm','varOut')
end
end
end

%% Interpolate HRTFs and generate results
fprintf('Generating results...\n')

%% First, run analysis for the reference (original HRTF)
filename = sprintf('%s/results/ref.mat',workdir);
if isfile(filename)
fprintf('\tFound %s. Loading...\n',filename)
load(filename,'results')
else
fprintf('\tGenerating results for reference HRTF...\n')

%% Numerical analysis
% Magnitude and phase delay 110-point Lebedev grid
results.mag = 20*log10(abs(H(:,indLeb,:)));
results.pd = (-unwrap(angle(H(:,indLeb,1)))+unwrap(angle(H(:,indLeb,2))))./(2*pi*f)*1e6;
% Loudness across directions
[L,wERB,calibrationGain] = perceptualSpectrum(H,fs);
Lavgd = sum(L.*wERB,1); % ERB-weighted avg loudness per direction
results.L = Lavgd;
results.L_perDirection = L; % save to calculate PSD
results.Lavg = sum(Lavgd(:,indLeb,:).*wLeb.',2); % weighted avg on the Lebedev grid
results.calibrationGain = calibrationGain; % save to calibrate loudness
% Interaural differences for horizontal plane
results.itd = itdestimator(permute(h(:,indHP,:),[2,3,1]),'fs',fs,...
'MaxIACCe','lp','upper_cutfreq', 3000,...
'butterpoly', 10)*1e6;
results.ild = getILD(h(:,indHP,:),fs);

%% Reijniers 2014
fprintf('\t\tRunning reijniers2014...\n'); tic
% Make DTF and SOFA object for Lebedev grid directions only
dtf = getDTF(h(:,indLeb,:),fs);
SOFA_obj = makesofa(dtf,fs,azLeb,elLeb);
% Preprocessing source information (demo_reijniers2014)
[template_loc, target] = reijniers2014_preproc(SOFA_obj);
% Run virtual experiments (demo_reijniers2014)
num_exp = 100;
[doa, params] = reijniers2014(template_loc, target, 'num_exp', num_exp);
% Calculate performance measures (demo_reijniers2014)
results.lat_acc = reijniers2014_metrics(doa, 'accL'); % mean lateral error
results.lat_prec = reijniers2014_metrics(doa, 'precL'); % lateral std
results.pol_acc = reijniers2014_metrics(doa, 'accP'); % mean polar error
results.pol_prec = reijniers2014_metrics(doa, 'precP'); % polar std
results.template_loc = template_loc; % template DTF
fprintf('\t\tFinished running reijniers2014 (took %0.2f s)...\n',toc);

%% Run baumgartner2020
fprintf('\t\tRunning baumgartner2020...\n'); tic
% Make DTF for median plane directions only
dtf = getDTF(h(:,indMP,:),fs);
ndirs = numel(elMP);
results.ext = nan(ndirs,1);
template_ext = cell(ndirs,1);
for j=1:ndirs
template_ext{j} = makesofa(dtf(:,j,:),fs,azMP(j),elMP(j));
% Get externalisation values
results.ext(j) = baumgartner2020(template_ext{j},template_ext{j});
end
results.template_ext = template_ext; % template DTFs
fprintf('\t\tFinished running baumgartner2020 (took %0.2f s)...\n',toc);

%% Run jelfs2011
fprintf('\t\tRunning jelfs2011...\n'); tic
ndirs = numel(azHP);
srm = nan(ndirs,1);
target = squeeze(h(:,indFront,:)); % target fixed at front
for j = 1:ndirs
interferer = squeeze(h(:,indHP(j),:)); % interferer moves around the HP
srm(j) = jelfs2011(target,interferer,fs);
end
results.srm = srm;
fprintf('\t\tFinished running jelfs2011 (took %0.2f s)...\n',toc);

%% Save results
fprintf('\t\tSaving results in %s...\n',filename)
save(filename,'results')
end

mag_ref = results.mag; % to calculate magnitude error
pd_ref = results.pd; % to calculate phase delay error
Lref = results.L_perDirection; % to calculate PSD
Lavg_ref = results.Lavg; % to level-match loudness
calibrationGain = results.calibrationGain; % to calibrate perceptual spectrum
template_loc = results.template_loc; % template localisation model
template_ext = results.template_ext; % template externalisation model
clear results

%% Then, for all test conditions
for N=N_vec % iterate through spatial orders
fprintf('Processing order %d...\n',N);
Y = [];
for i=1:ncond % iterate through conditions
name = test_conditions{i}.name;
filename = sprintf('%s/results/ord%0.2d_%s.mat',workdir,N,name);
if isfile(filename)
fprintf('\tFound %s. Skipping...\n',filename)
else
fprintf('\tGenerating results...\n')

if isempty(Y)
Y = AKsh(N, [], az*180/pi, el*180/pi, 'real').';
end

%% Load hnm and interpolate to test grid
hnm = load(sprintf('%s/hnm/ord%0.2d_%s.mat',workdir,N,name)).hnm;
isaligned = strcmp(name,'TA');
hInterp = fromSH(hnm,fs,az,el,isaligned,r);
HInterp = ffth(hInterp);

%% Numerical analysis
% Magnitude and phase delay error 110-point Lebedev grid
mag = 20*log10(abs(HInterp(:,indLeb,:)));
pd = (-unwrap(angle(HInterp(:,indLeb,1)))+unwrap(angle(HInterp(:,indLeb,2))))./(2*pi*f)*1e6;
results.err_mag = mean(abs(mag-mag_ref),2); % avg abs difference across directions
results.err_pd = mean(abs(pd-pd_ref),2); % avg abs difference across directions
% Loudness across directions
[L,wERB] = perceptualSpectrum(HInterp,fs,calibrationGain);
Lavgd = sum(L.*wERB,1); % ERB-weighted avg loudness per direction
Lavg = sum(Lavgd(:,indLeb,:).*wLeb.',2); % weighted avg on the Lebedev grid
L = L - Lavg + Lavg_ref; % level-match with reference
Lavgd = sum(L.*wERB,1); % recalculate avg loudness per direction
Ldif = L-Lref; % loudness difference
PSD = sum(abs(L-Lref).*wERB,1); % avg abs difference per direction (PSD)
PSDavg = sum(PSD(:,indLeb,:).*wLeb.',2); % avg PSD over the Lebedev grid
results.L = Lavgd;
results.PSD = PSD;
results.PSDavg = PSDavg;
% Interaural differences for horizontal plane
results.itd = itdestimator(permute(hInterp(:,indHP,:),[2,3,1]),'fs',fs,...
'MaxIACCe','lp','upper_cutfreq', 3000,...
'butterpoly', 10)*1e6;
results.ild = getILD(hInterp(:,indHP,:),fs);

%% Reijniers 2014
fprintf('\t\tRunning reijniers2014...\n'); tic
% Make DTF and SOFA object for Lebedev grid directions only
dtf = getDTF(hInterp(:,indLeb,:),fs);
SOFA_obj = makesofa(dtf,fs,azLeb,elLeb);
% Preprocessing source information (demo_reijniers2014)
[~, target] = reijniers2014_preproc(SOFA_obj);
% Run virtual experiments (demo_reijniers2014)
num_exp = 100;
[doa, params] = reijniers2014(template_loc, target, 'num_exp', num_exp);
% Calculate performance measures (demo_reijniers2014)
results.lat_acc = reijniers2014_metrics(doa, 'accL'); % mean lateral error
results.lat_prec = reijniers2014_metrics(doa, 'precL'); % lateral std
results.pol_acc = reijniers2014_metrics(doa, 'accP'); % mean polar error
results.pol_prec = reijniers2014_metrics(doa, 'precP'); % polar std
fprintf('\t\tFinished running reijniers2014 (took %0.2f s)...\n',toc);

%% Run baumgartner2020
fprintf('\t\tRunning baumgartner2020...\n'); tic
% Make DTF for median plane directions only
dtf = getDTF(hInterp(:,indMP,:),fs);
ndirs = numel(elMP);
results.ext = nan(ndirs,1);
for j=1:ndirs
target = makesofa(dtf(:,j,:),fs,azMP(j),elMP(j));
% Get externalisation values
results.ext(j) = baumgartner2020(target,template_ext{j});
end
fprintf('\t\tFinished running baumgartner2020 (took %0.2f s)...\n',toc);

%% Run jelfs2011
fprintf('\t\tRunning jelfs2011...\n'); tic
ndirs = numel(azHP);
srm = nan(ndirs,1);
target = squeeze(hInterp(:,indFront,:)); % target fixed at front
for j = 1:ndirs
interferer = squeeze(hInterp(:,indHP(j),:)); % interferer moves around the HP
srm(j) = jelfs2011(target,interferer,fs);
end
results.srm = srm;
fprintf('\t\tFinished running jelfs2011 (took %0.2f s)...\n',toc);

%% Save results
fprintf('\t\tSaving results in %s...\n',filename)
save(filename,'results')
end
end
end

%% Plot Fig. 1 (SH spectra)
fig1size = [17.78 7];
fig1 = figure('units','centimeters','PaperUnits','centimeters',...
'PaperSize',fig1size,'Renderer','painters',...
'pos',[2 2 fig1size(1) fig1size(2)],...
'paperposition',[0 0 fig1size(1) fig1size(2)]);
% Tight subplot
gap = [.02 .008]; % gap between subplots in norm units (height width)
marg_h = [.18 .15]; % figure height margins in norm units (lower upper)
marg_w = [.05 .08]; % figure width margins in norm units (left right)
[ha, pos] = tight_subplot(1,2,gap,marg_h,marg_w);
% Fig. 1a: original HRTF
axes(ha(1)), plotSHenergy(Hnm_ref(:,:,1),fs); % left only
xlabel('f (Hz)'), ylabel('Order (n)')
axpos = get(ha(1),'pos');
annotation(fig1,'textbox',...
[axpos(1)+axpos(3)-0.2 axpos(2)+axpos(4)-0.09 0.2 0.09],...
'String',{'a) Original'},...
'HorizontalAlignment','right',...
'FontWeight','bold',...
'FontSize',7,...
'FitBoxToText','off',...
'EdgeColor','none');
set(gca,'fontsize',7)
% Fig. 1b: time-aligned HRTF
axes(ha(2)), plotSHenergy(Hnm_TA_ref(:,:,1),fs);
xlabel('f (Hz)'), ylabel(''), yticklabels({})
axpos = get(ha(2),'pos');
annotation(fig1,'textbox',...
[axpos(1)+axpos(3)-0.2 axpos(2)+axpos(4)-0.09 0.2 0.09],...
'String',{'b) Time-aligned'},...
'HorizontalAlignment','right',...
'FontWeight','bold',...
'FontSize',7,...
'FitBoxToText','off',...
'EdgeColor','none');
set(gca,'fontsize',7)
sgtitle('HRTF energy in the SH domain, with and without alignment')

%% Plot Fig. 2 (mag/phase errors per spatial order)
% Load data from file
names = {
'ord01_Trunc'
'ord05_Trunc'
'ord10_Trunc'
'ord20_Trunc'
'ord30_Trunc'
'ord40_Trunc'
};
labels = {
'N=1'
'N=5'
'N=10'
'N=20'
'N=30'
'N=40'
};
n = numel(names);
err_mag = zeros(nfreqs,n);
err_pd = zeros(nfreqs,n);
for i=1:n
name = names{i};
filename = sprintf('%s/results/%s.mat',workdir,name);
load(filename,'results');
err_mag(:,i) = results.err_mag(:,1,1); % left ear only
err_pd(:,i) = results.err_pd;
end

% Tight subplot
fig2 = figure('pos',[56.6000 98.6000 323.2000 319.6000]);
gap = [.02 .008]; % gap between subplots in norm units (height width)
marg_h = [.09 .1]; % [.09 .1] figure height margins in norm units (lower upper)
marg_w = [.1 .03]; % figure width margins in norm units (left right)
[ha, ~] = tight_subplot(2,1,gap,marg_h,marg_w);

% Top plot: magnitude error
colors = parula(n+1);
lsvec = {'-'};%,'-.'};
lwvec = [0.5,0.5];
mvec = {'^','v','x','s','o','d','p','h'};
ms = 3; % marker size
mi = int32(logspace(log10(1),log10(1025),10));

axes(ha(1))
for i=1:n
ls = lsvec{mod(i-1,numel(lsvec))+1};
lw = lwvec(mod(i-1,numel(lwvec))+1);
m = mvec{i};
semilogx(f,err_mag(:,i),'Color',colors(i,:),'LineWidth',lw,'LineStyle',ls,'Marker',m,'MarkerSize',ms,'MarkerIndices',mi); hold on
end

semilogx([f(2) 20000], [1 1],'k:')
legend(labels,'location','west')

ylim([0 40])
xlabel(''), xticklabels({}), xlim([f(2) 20000])
ylabel('Error (dB)'), grid on
axpos = get(ha(1),'pos');
annotation(fig2,'textbox',...
[axpos(1) axpos(2)+axpos(4)-0.09 0.5 0.09],...
'String',{'a. Magnitude error'},...
'HorizontalAlignment','left',...
'FontWeight','bold',...
'FontSize',7,...
'FitBoxToText','off',...
'EdgeColor','none');

% Bottom plot: phase error
axes(ha(2))
for i=1:n
ls = lsvec{mod(i-1,numel(lsvec))+1};
lw = lwvec(mod(i-1,numel(lwvec))+1);
m = mvec{i};
semilogx(f,err_pd(:,i),'Color',colors(i,:),'LineWidth',lw,'LineStyle',ls,'Marker',m,'MarkerSize',ms,'MarkerIndices',mi); hold on
end

semilogx([f(2) 20000], [20 20],'k:')
ylim([0 300])
xlim([f(2) 20000]), grid on, ylabel('Error (\mus)')
xticks([100,1000,10000,20000])
xticklabels(["100","1k","10k","20k"])
xlabel('f (Hz)')
axpos = get(ha(2),'pos');
annotation(fig2,'textbox',...
[axpos(1) axpos(2)+axpos(4)-0.09 0.5 0.09],...
'String',{'b. Phase delay error'},...
'HorizontalAlignment','left',...
'FontWeight','bold',...
'FontSize',7,...
'FitBoxToText','off',...
'EdgeColor','none');

set(ha,'fontsize',7)
set(gcf,'units','centimeters','Renderer','painters')
figpos = get(gcf,'position');
set(gcf,'PaperSize',[figpos(3) figpos(4)],'paperposition',[0 0 figpos(3) figpos(4)])
sgtitle('Interpolation errors (non-aligned HRTF)')

%% Plot Fig. 3 (mag/phase errors per spatial order, time-aligned)
% Load data from file
names = {
'ord01_TA'
'ord05_TA'
'ord10_TA'
'ord20_TA'
'ord30_TA'
'ord40_TA'
};
labels = {
'N=1'
'N=5'
'N=10'
'N=20'
'N=30'
'N=40'
};
n = numel(names);
err_mag = zeros(nfreqs,n);
err_pd = zeros(nfreqs,n);
for i=1:n
name = names{i};
filename = sprintf('%s/results/%s.mat',workdir,name);
load(filename,'results');
err_mag(:,i) = results.err_mag(:,1,1); % left ear only
err_pd(:,i) = results.err_pd;
end

% Tight subplot
fig3 = figure('pos',[56.6000 98.6000 323.2000 319.6000]);
gap = [.02 .008]; % gap between subplots in norm units (height width)
marg_h = [.09 .1]; % [.09 .1] figure height margins in norm units (lower upper)
marg_w = [.1 .03]; % figure width margins in norm units (left right)
[ha, ~] = tight_subplot(2,1,gap,marg_h,marg_w);

% Top plot: magnitude error
colors = parula(n+1);
lsvec = {'-'};%,'-.'};
lwvec = [0.5,0.5];
mvec = {'^','v','x','s','o','d','p','h'};
ms = 3; % marker size
mi = int32(logspace(log10(1),log10(1025),10));

axes(ha(1))
for i=1:n
ls = lsvec{mod(i-1,numel(lsvec))+1};
lw = lwvec(mod(i-1,numel(lwvec))+1);
m = mvec{i};
semilogx(f,err_mag(:,i),'Color',colors(i,:),'LineWidth',lw,'LineStyle',ls,'Marker',m,'MarkerSize',ms,'MarkerIndices',mi); hold on
end

semilogx([f(2) 20000], [1 1],'k:')
legend(labels,'location','west')

ylim([0 40])
xlabel(''), xticklabels({}), xlim([f(2) 20000])
ylabel('Error (dB)'), grid on
axpos = get(ha(1),'pos');
annotation(fig2,'textbox',...
[axpos(1) axpos(2)+axpos(4)-0.09 0.5 0.09],...
'String',{'a. Magnitude error'},...
'HorizontalAlignment','left',...
'FontWeight','bold',...
'FontSize',7,...
'FitBoxToText','off',...
'EdgeColor','none');

% Bottom plot: phase error
axes(ha(2))
for i=1:n
ls = lsvec{mod(i-1,numel(lsvec))+1};
lw = lwvec(mod(i-1,numel(lwvec))+1);
m = mvec{i};
semilogx(f,err_pd(:,i),'Color',colors(i,:),'LineWidth',lw,'LineStyle',ls,'Marker',m,'MarkerSize',ms,'MarkerIndices',mi); hold on
end

semilogx([f(2) 20000], [20 20],'k:')
ylim([0 300])
xlim([f(2) 20000]), grid on, ylabel('Error (\mus)')
xticks([100,1000,10000,20000])
xticklabels(["100","1k","10k","20k"])
xlabel('f (Hz)')
axpos = get(ha(2),'pos');
annotation(fig2,'textbox',...
[axpos(1) axpos(2)+axpos(4)-0.09 0.5 0.09],...
'String',{'b. Phase delay error'},...
'HorizontalAlignment','left',...
'FontWeight','bold',...
'FontSize',7,...
'FitBoxToText','off',...
'EdgeColor','none');

set(ha,'fontsize',7)
set(gcf,'units','centimeters','Renderer','painters')
figpos = get(gcf,'position');
set(gcf,'PaperSize',[figpos(3) figpos(4)],'paperposition',[0 0 figpos(3) figpos(4)])
sgtitle('Interpolation errors (time-aligned HRTF)')

%% Plot Fig. 3 (ITD/ILD per spatial order)
% Load data from file
names = {
'ord01_Trunc'
'ord05_Trunc'
'ord10_Trunc'
'ord20_Trunc'
'ord30_Trunc'
'ord40_Trunc'
'ref'
};
labels = {
'N=1'
'N=5'
'N=10'
'N=20'
'N=30'
'N=40'
'Reference'
};
n = numel(names);
ndirs = numel(azHP);
itd = zeros(ndirs,n);
ild = zeros(ndirs,n);
for i=1:n
name = names{i};
filename = sprintf('%s/results/%s.mat',workdir,name);
load(filename,'results');
itd(:,i) = results.itd;
ild(:,i) = results.ild;
end

% Tight subplot
fig4size = [21 13];
fig4 = figure('units','centimeters','pos',[2 2 fig4size(1) fig4size(2)],...
'Renderer','painters','PaperSize',[fig4size(1) fig4size(2)],...
'paperposition',[0 0 fig4size(1) fig4size(2)]);
gap = [.04 .06]; % gap between subplots in norm units (height width)
marg_h = [0.05 .15]; % [.09 .1] figure height margins in norm units (lower upper)
marg_w = [.07 .03]; % figure width margins in norm units (left right)
[ha, ~] = tight_subplot(2,2,gap,marg_h,marg_w);

colors = parula(n);
lsvec = {'-'};%,'-.'};
lwvec = [0.5,0.5];
mvec = {'^','v','x','s','o','d','p','h'};
ms = 3; % marker size
step_big = round(numel(azHP)/4);
step_small = round(step_big/n);

% ITD
axes(ha(1))
for i=1:n % multiply by 0.001 to have it in ms (cleaner plot)
ls = lsvec{mod(i-1,numel(lsvec))+1};
lw = lwvec(mod(i-1,numel(lwvec))+1);
m = mvec{i};
mi = [((i-1)*step_small+1):step_big:numel(azHP)-1]; % marker indices
color = colors(i,:);
if i==n % reference
ls = ':';
lw = 0.5;
m = 'none';
color = [0 0 0];
end
polarplot(azHP,abs(itd(:,i))*0.001,'Color',color,'LineWidth',lw,'LineStyle',ls,'Marker',m,'MarkerSize',ms,'MarkerIndices',mi), hold on

end
set(gca,'ThetaZeroLocation','top','ThetaLim',[0 360],'fontsize',7)
set(gca,'ThetaTick',[0:45:360],'ThetaTickLabel',["0","45","90","135","","-135","-90","-45"],'RAxisLocation',-90)
title('ITD (ms)')

% ILD
axes(ha(2))
for i=1:n
ls = lsvec{mod(i-1,numel(lsvec))+1};
lw = lwvec(mod(i-1,numel(lwvec))+1);
m = mvec{i};
mi = [((i-1)*step_small+1):step_big:numel(azHP)-1]; % marker indices
color = colors(i,:);
if i==n % reference
ls = ':';
lw = 0.5;
m = 'none';
color = [0 0 0];
end
polarplot(azHP,abs(ild(:,i)),'Color',color,'LineWidth',lw,'LineStyle',ls,'Marker',m,'MarkerSize',ms,'MarkerIndices',mi), hold on
end
set(gca,'ThetaZeroLocation','top','ThetaLim',[0 360],'fontsize',7)
set(gca,'ThetaTick',[0:45:360],'ThetaTickLabel',["0","45","90","135","","-135","-90","-45"],'RAxisLocation',-90)
title('ILD (dB)')

legend(labels,'position',[0.4646 0.6219 0.1135 0.1883]);

itderr = abs(itd-itd(:,end));
ilderr = abs(ild-ild(:,end));

axes(ha(3)), violinplot(itderr(:,1:end-1),[],'ShowData',false,'BoxWidth',0.03); grid on
set(gca,'YTickLabelMode','auto')
hold on, plot([0 n],[itdJND itdJND],'k:')
xticklabels(labels(1:end-1))
ylabel('Abs. ITD error (ms)')
set(gca,'fontsize',6)
ylim([0 700])

axes(ha(4)), violinplot(ilderr(:,1:end-1),[],'ShowData',false,'BoxWidth',0.03); grid on
set(gca,'YTickLabelMode','auto')
hold on, plot([0 n],[ildJND ildJND],'k:')
xticklabels(labels(1:end-1))
ylabel('Abs. ILD error (dB)')
set(gca,'fontsize',6)
ylim([0 10])

sgtitle('Interaural differences (non-aligned HRTF)')

%% Plot Fig. 5 (ITD/ILD per spatial order)
% Load data from file
names = {
'ord01_TA'
'ord05_TA'
'ord10_TA'
'ord20_TA'
'ord30_TA'
'ord40_TA'
'ref'
};
labels = {
'N=1'
'N=5'
'N=10'
'N=20'
'N=30'
'N=40'
'Reference'
};
n = numel(names);
ndirs = numel(azHP);
itd = zeros(ndirs,n);
ild = zeros(ndirs,n);
for i=1:n
name = names{i};
filename = sprintf('%s/results/%s.mat',workdir,name);
load(filename,'results');
itd(:,i) = results.itd;
ild(:,i) = results.ild;
end

% Tight subplot
fig5 = figure('units','centimeters','pos',[2 2 fig4size(1) fig4size(2)],...
'Renderer','painters','PaperSize',[fig4size(1) fig4size(2)],...
'paperposition',[0 0 fig4size(1) fig4size(2)]);
gap = [.04 .06]; % gap between subplots in norm units (height width)
marg_h = [0.05 .15]; % [.09 .1] figure height margins in norm units (lower upper)
marg_w = [.07 .03]; % figure width margins in norm units (left right)
[ha, ~] = tight_subplot(2,2,gap,marg_h,marg_w);

colors = parula(n);
lsvec = {'-'};%,'-.'};
lwvec = [0.5,0.5];
mvec = {'^','v','x','s','o','d','p','h'};
ms = 3; % marker size
step_big = round(numel(azHP)/4);
step_small = round(step_big/n);

% ITD
axes(ha(1))
for i=1:n % multiply by 0.001 to have it in ms (cleaner plot)
ls = lsvec{mod(i-1,numel(lsvec))+1};
lw = lwvec(mod(i-1,numel(lwvec))+1);
m = mvec{i};
mi = [((i-1)*step_small+1):step_big:numel(azHP)-1]; % marker indices
color = colors(i,:);
if i==n % reference
ls = ':';
lw = 0.5;
m = 'none';
color = [0 0 0];
end
polarplot(azHP,abs(itd(:,i))*0.001,'Color',color,'LineWidth',lw,'LineStyle',ls,'Marker',m,'MarkerSize',ms,'MarkerIndices',mi), hold on

end
set(gca,'ThetaZeroLocation','top','ThetaLim',[0 360],'fontsize',7)
set(gca,'ThetaTick',[0:45:360],'ThetaTickLabel',["0","45","90","135","","-135","-90","-45"],'RAxisLocation',-90)
title('ITD (ms)')

% ILD
axes(ha(2))
for i=1:n
ls = lsvec{mod(i-1,numel(lsvec))+1};
lw = lwvec(mod(i-1,numel(lwvec))+1);
m = mvec{i};
mi = [((i-1)*step_small+1):step_big:numel(azHP)-1]; % marker indices
color = colors(i,:);
if i==n % reference
ls = ':';
lw = 0.5;
m = 'none';
color = [0 0 0];
end
polarplot(azHP,abs(ild(:,i)),'Color',color,'LineWidth',lw,'LineStyle',ls,'Marker',m,'MarkerSize',ms,'MarkerIndices',mi), hold on
end
set(gca,'ThetaZeroLocation','top','ThetaLim',[0 360],'fontsize',7)
set(gca,'ThetaTick',[0:45:360],'ThetaTickLabel',["0","45","90","135","","-135","-90","-45"],'RAxisLocation',-90)
title('ILD (dB)')

legend(labels,'position',[0.4646 0.6219 0.1135 0.1883]);

itderr = abs(itd-itd(:,end));
ilderr = abs(ild-ild(:,end));

axes(ha(3)), violinplot(itderr(:,1:end-1),[],'ShowData',false,'BoxWidth',0.03); grid on
set(gca,'YTickLabelMode','auto')
hold on, plot([0 n],[itdJND itdJND],'k:')
xticklabels(labels(1:end-1))
ylabel('Abs. ITD error (ms)')
set(gca,'fontsize',6)
ylim([0 700])

axes(ha(4)), violinplot(ilderr(:,1:end-1),[],'ShowData',false,'BoxWidth',0.03); grid on
set(gca,'YTickLabelMode','auto')
hold on, plot([0 n],[ildJND ildJND],'k:')
xticklabels(labels(1:end-1))
ylabel('Abs. ILD error (dB)')
set(gca,'fontsize',6)
ylim([0 10])

sgtitle('Interaural differences (time-aligned HRTF)')

%% Plot Fig. 6 (loudness map)
% Load data from file
names = {
'ref'
'ord05_Trunc'
'ord05_TA'
};
labels = {
'Reference'
'Trunc'
'TA'
};
n = numel(names);
ndirs = numel(az);
L = cell(n,1);
PSD = zeros(ndirs,n);
PSDavg = zeros(n,1);
for i=1:n
name = names{i};
filename = sprintf('%s/results/%s.mat',workdir,name);
load(filename,'results');
L{i} = results.L;
if isfield(results,'PSD')
PSD(:,i) = results.PSD(:,:,1).'; % left ear only
PSDavg(i) = results.PSDavg(:,:,1); % left ear only
end
end

% Tight subplot
fig6size = [16.4253 8];
fig6 = figure('units','centimeters','pos',[2 2 fig6size(1) fig6size(2)],...
'Renderer','painters','PaperSize',[fig6size(1) fig6size(2)],...
'paperposition',[0 0 fig6size(1) fig6size(2)]);
gap = [.07 .008]; % gap between subplots in norm units (height width)
marg_h = [.5 .15]; % [.09 .1] figure height margins in norm units (lower upper)
marg_w = [.06 .1]; % figure width margins in norm units (left right)
[ha, ~] = tight_subplot(1,3,gap,marg_h,marg_w);
clims = [1 5]; % change as needed

colormap parula
for i=1:n
axes(ha(i))
plotSph(az,el,L{i}(:,:,1)) % left ear only
if i==1 || i==6
yticks(-60:30:60)
ylabel('Elevation (deg)')
else
yticklabels({})
end
xlabel('Azimuth (deg)')
if PSDavg(i)>0
title(sprintf('%s [%0.2f]',labels{i},PSDavg(i)));
else
title(labels{i}) % reference
end

set(gca,'fontsize',7,'XDir','reverse')
grid(gca,'off')
caxis(clims)
if i==n
c = colorbar; c.Label.String = 'Loudness (sones)';
c.Position = [0.9174 0.4904 0.0111 0.3804];
end
end

ax = axes('innerposition',[0.4 0.05 0.2 0.3]);
violinplot(PSD(:,2:end),[],'ShowData',false,'BoxWidth',0.03); grid on
% boxplot(PSD(:,2:end)), grid on
% hold on, plot(mean(PSD(:,2:end)),'dg')
xticklabels(labels(2:end))
ylabel('PSD (sones)')
set(gca,'fontsize',7)

sgtitle('Loudness maps and perceptual spectral difference (PSD)')

%% Plot Fig. 7 (models' outputs per spatial order)
% Load data from file
names = {
'Trunc'
'TA'
'ref'
};
labels = {
'Trunc'
'TA'
'Reference'
};
n = numel(names);
m = numel(N_vec);
PSD = nan(n,m);
lat_prec = nan(n,m);
pol_prec = nan(n,m);
ext = nan(n,m);
srm = nan(n,m);
for i=1:n
name = names{i};
if strcmp(name,'ref')
filename = sprintf('%s/results/ref.mat',workdir);
load(filename,'results');
PSD(i,:) = 0;
lat_prec(i,:) = results.lat_prec;
pol_prec(i,:) = results.pol_prec;
ext(i,:) = mean(results.ext);
srm(i,:) = mean(results.srm);
else
for j=1:m
N=N_vec(j);
filename = sprintf('%s/results/ord%0.2d_%s.mat',workdir,N,name);
load(filename,'results');
PSD(i,j) = results.PSDavg(:,:,1); % left ear only
lat_prec(i,j) = results.lat_prec;
pol_prec(i,j) = results.pol_prec;
ext(i,j) = mean(results.ext);
srm(i,j) = mean(results.srm);
end
end
end
% Tight subplot
fig7size = [18.3938 9.7790];
fig7 = figure('units','centimeters','pos',[2 2 fig7size(1) fig7size(2)],...
'Renderer','painters','PaperSize',[fig7size(1) fig7size(2)],...
'paperposition',[0 0 fig7size(1) fig7size(2)]);
gap = [.06 .05]; % gap between subplots in norm units (height width)
marg_h = [.1 .1]; % [.09 .1] figure height margins in norm units (lower upper)
marg_w = [.07 .02]; % figure width margins in norm units (left right)
[ha, ~] = tight_subplot(2,3,gap,marg_h,marg_w);

colors = parula(n);
lsvec = {'-'};%,'-.'};
mvec = {'^','v','x','s','o','d','p','h'};
lwvec = [0.5,0.5];
ms = 3; % marker size
mi = 1:n; % [1:5:44]; % marker indices

colors(end,:) = [0 0 0]; % last one is reference

for i=1:n
if i<n
ls = lsvec{mod(i-1,numel(lsvec))+1};
lw = lwvec(mod(i-1,numel(lwvec))+1);
m = mvec{i};
else
ls = ':'; % reference
lw = 0.5;
m = 'None';
end
plot(ha(1),N_vec,PSD(i,:),'Color',colors(i,:),'LineWidth',lw,'LineStyle',ls,'Marker',m,'MarkerSize',ms,'MarkerIndices',mi),hold(ha(1),'on')
plot(ha(2),N_vec,lat_prec(i,:),'Color',colors(i,:),'LineWidth',lw,'LineStyle',ls,'Marker',m,'MarkerSize',ms,'MarkerIndices',mi),hold(ha(2),'on')
plot(ha(3),N_vec,pol_prec(i,:),'Color',colors(i,:),'LineWidth',lw,'LineStyle',ls,'Marker',m,'MarkerSize',ms,'MarkerIndices',mi),hold(ha(3),'on')
plot(ha(4),N_vec,ext(i,:),'Color',colors(i,:),'LineWidth',lw,'LineStyle',ls,'Marker',m,'MarkerSize',ms,'MarkerIndices',mi),hold(ha(4),'on')
plot(ha(5),N_vec,srm(i,:),'Color',colors(i,:),'LineWidth',lw,'LineStyle',ls,'Marker',m,'MarkerSize',ms,'MarkerIndices',mi),hold(ha(5),'on')
end
legh = legend(ha(1),labels,'location','best');

ylabel(ha(1),'PSD (sones)')%, ylim(ha(1),[min(PSD(:),max(PSD(:)])
ylabel(ha(2),'Lateral precision (deg)')
ylabel(ha(3),'Polar precision (deg)')
ylabel(ha(4),'Externalisation')%, ylim(ha(4),[0.3 1.01])
ylabel(ha(5),'SRM (dB)') % , ylim(ha(4),[0.3 1.01])
for i=1:5
grid(ha(i),'on')
xlim(ha(i),[1 44])
xticks(ha(i),[1,5:5:44])
if i<=3
% xticklabels(ha(i),{})
else
xlabel(ha(i),'Spatial order (N)')
end
set(ha(i),'fontsize',7)
end
set(ha(6),'visible','off')

sgtitle('Perceptual models'' output (per order)')