% Plot estimated loudness maps of two HRTFs, calculated with the 
% perceptually-weighted magnitude spectrum proposed in [1].
%
% EXTERNAL DEPENDENCIES:
%   SOFA API for Matlab (github.com/sofacoustics/API_MO)
%
% REFERENCES:
%   [1] Armstrong, Calum, et al. "A perceptual spectral difference model
%       for binaural signals." Audio Engineering Society Convention 145.
%       Audio Engineering Society, 2018.
%
% AUTHOR: Isaac Engel - isaac.engel(at)imperial.ac.uk
% February 2021

%% Clear variables
clear

%% Parameters
filenames = { % put here the paths to the HRTFs' SOFA files
    '../hrtfs/FABIAN_HRIR_measured_HATO_0.sofa'
    '../hrtfs/D2_44K_16bit_256tap_FIR_SOFA.sofa'
};
labels = { % put here the names for the legend
    'FABIAN'
    'KEMAR'
};
n = numel(filenames);
nplotsH = 2; % number of plots in horizontal
nplotsV = 1; % number of plots in vertical
assert(n<=nplotsH*nplotsV,'Too few plots!')
clims = [];% colorbar limits, [] for automatic

%% Prepare figure
figsize = [16.4253 6.9003]; % change as appropriate
fig = figure('units','centimeters','pos',[2 2 figsize(1) figsize(2)],...
    'Renderer','painters','PaperSize',[figsize(1) figsize(2)],...
    'paperposition',[0 0 figsize(1) figsize(2)]);
gap = [.07 .008]; % gap between subplots in norm units (height width)
marg_h = [.13 .06]; % figure height margins in norm units (lower upper)
marg_w = [.06 .11]; % figure width margins in norm units (left right)
[ha, ~] = tight_subplot(nplotsV,nplotsH,gap,marg_h,marg_w);
colormap parula

%% Load HRTFs and plot loudness for all directions
for i=1:n
    SOFA_obj = SOFAload(filenames{i});
    [h,fs,az,el,r] = sofa2hrtf(SOFA_obj);
    H = ffth(h); % to frequency domain
    
    % Calculate loudness per direction
    [L,wERB] = perceptualSpectrum(H,fs);
    Lavgd = sum(L.*wERB,1); % ERB-weighted avg loudness per direction

    axes(ha(i))
    plotSph(az,el,Lavgd(:,:,1)) % plot left ear only
    
    % Some plot settings
    title(labels{i}) % reference
    set(gca,'fontsize',7,'XDir','reverse')
    grid(gca,'off')
    if isempty(clims)
        clims = caxis;
    end
    caxis(clims)
    if i==n % colorbar only for the last plot
        c = colorbar; c.Label.String = 'Loudness (sones)';
        c.Position = [0.9 0.13 0.02 0.81]; % change as appropriate
    end
    if mod(i-1,nplotsH)==0 % show Y axis only for plots on the left
        yticks(-60:30:60)
        ylabel('Elevation (deg)')
    else
        yticklabels({})
    end
    if i>(nplotsH*(nplotsV-1)) % show X axis only for plots on last row
        xlabel('Azimuth (deg)')
    else
        xticklabels({})
    end
end
for i=(n+1):(nplotsH*nplotsV)
    axis(ha(i),'off') % hide unused axes
end
