%% Read file in SOFA format

% César D. Salvador
% cesardsalvador@gmail.com
% http://www.ais.riec.tohoku.ac.jp/~salvador/

% This script requires the SOFA Matlab/Octave API version 1.0.2
% https://sourceforge.net/projects/sofacoustics/

close all
clear all
clc

SOFAstart;

%listener_model = 'gen01_ht';
%listener_model = 'gen01_h';
%listener_model = 'gen02_h';
listener_model = 'ind01_h';
%listener_model = 'ind02_h';

source_distribution = 'sph_equivert_dir614_dist91';


%% Structure array Obj from SOFA file
sofa_folder = 'D:\prog\matlab\data\hrir\sofa\';			% set your local folder
sofa_file_name = [listener_model, '_', source_distribution,'.sofa'];
Obj = SOFAload([sofa_folder, sofa_file_name]);
Obj = orderfields(Obj);


%% Dimensions of data
Fs = Obj.Data.SamplingRate;                             % sampling frequency in Hertz
Ns = Obj.API.N;                                         % number of samples along time describing a measurement for one position
sample_circshift = Obj.Data.Delay;                      % number of samples used in circular shift of HRIRs (circular delay)
Nazim = Obj.SourcePositionNumberAzimuths;               % number of source azimuths
Nelev = Obj.SourcePositionNumberElevations;             % number of source elevations (without the north and south poles)
Ndir = Obj.SourcePositionNumberDirections;              % number of source directions (Ndir = Nazim * Nelev + 2)
Ndist = Obj.SourcePositionNumberDistances;              % number of source distances
Npos = Obj.API.M;                                       % number of source positions around the listener (Npos = Ndir * Ndist)
c = Obj.SpeedOfSoundInMPS;                              % speed of sound in air (m/s)


%% One-dimensional domains
t = (0:Ns-1)*1/Fs;                                      % time (seconds)
f = (0:Ns/2)*Fs/Ns;                                     % frequency (Hertz)
k = 2*pi*f/c;                                           % wave number (1/meter)


%% Source positions in spherical and cartesian coordinates
azim = Obj.SourcePosition(:, 1);
elev = Obj.SourcePosition(:, 2);
dist = Obj.SourcePosition(:, 3);
[x(:, 1), x(:, 2), x(:, 3)] = sph2cart(azim*pi/180, elev*pi/180, dist);
azim = reshape(azim, [Ndir Ndist]);                     % azimuthal angle in degree (size: Ndir * Ndist)
elev = reshape(elev, [Ndir Ndist]);                     % elevation angle in degree (size: Ndir * Ndist)
dist = reshape(dist, [Ndir Ndist]);                     % radial distance in metre (size: Ndir * Ndist)
x = reshape(x, [Ndir Ndist 3]);                         % cartesian coordinates in metre (size: Ndir * Ndist * 3)


%% Head-related impulse responses (size: Ns * Ndir * Ndist)
HRIR_left = permute(reshape(squeeze(Obj.Data.IR(:, 1, :)), [Ndir Ndist Ns]), [3 1 2]);
HRIR_right = permute(reshape(squeeze(Obj.Data.IR(:, 2, :)), [Ndir Ndist Ns]), [3 1 2]);


%% Head-related transfer functions (size: Ns/2+1 * Ndir * Ndist)
HRTF_left = fft(HRIR_left); HRTF_left = HRTF_left(1:Ns/2+1, :, :);
HRTF_right = fft(HRIR_right); HRTF_right = HRTF_right(1:Ns/2+1, :, :);
HRTF_interaural = HRTF_left./HRTF_right;


%% Parameters for visualization
font_size = 14;
set(0, 'defaultfigurecolor', 'w', ...
    'defaultaxesfontname',  'times', ...
    'defaulttextfontSize', font_size, ...
    'defaultaxesfontsize', font_size, ...
    'defaultaxeslayer', 'top')

clim_mono = [-18 18];
ctick_mono = clim_mono(1):6:clim_mono(2);
clim_inter = [0 39];
ctick_inter = clim_inter(1):6:clim_inter(2);
view_angle = [142 25];
fmax = 20;
fmin = 0.1;
flim = [fmin fmax];
ftick_log = 2.^(-2:4);
distmin = 10;
distmax = 100;
distlim = [distmin distmax];
disttick = [12.5 25 50 100];


%% Visualization of HRTFs for a single distance and a single frequency
fview = 10000;                                           % frequency in Hertz to visualize
dview = 0.25;                                            % distance in metre to visualize

[~, indfreq] = min(abs(f - fview));
[~, inddist] = min(abs(dist(1, :) - dview));

a = squeeze(x(:, inddist, :));
msh = convhull(a);
V = squeeze(20*log10(abs(HRTF_left(indfreq, :, inddist))));

figure
trisurf(msh, a(:, 1), a(:, 2), a(:, 3), V)
%shading interp
title(['Left-ear HRTFs for ', num2str(dist(1, inddist)*1e2), ' cm and ', sprintf('%.1f', f(indfreq)*1e-3), ' kHz'], 'interpreter', 'latex')
xlabel('$x$ (m)', 'interpreter', 'latex')
ylabel('$y$ (m)', 'interpreter', 'latex')
zlabel('$z$ (m)', 'interpreter', 'latex')
axis tight
axis equal
axis vis3d
box on
set(gca, 'clim', clim_mono)
titstr = '$20\log_{10}|H_{\mbox{left}}|$';
hcb = colorbar('eastoutside', 'ytick', ctick_mono, 'xticklabel', ctick_mono, 'ylim', clim_mono);
set(get(hcb, 'ylabel'), 'String', titstr, 'interpreter', 'latex', 'fontsize', font_size);
view(view_angle)


%% Visualization of HRTFs for a single direction
azimview = 100;                                           % azimuth in degree to visualize
elevview = 20;                                            % elevation in degree to visualize

[~, inddir] = min( sqrt( (azim(:, 1) - azimview).^2 + (elev(:, 1) - elevview).^2 ) );

V = 20*log10(abs(squeeze(HRTF_interaural(:, inddir, :))));

figure
surface(dist(1, :)*1e2, f*1e-3, V, 'edgecolor', 'none');
shading interp
title(['Interaural HRTFs for azim. ', sprintf('%.1f', azim(inddir, 1)), '$^\circ$ and elev. ', sprintf('%.1f', elev(inddir, 1)), '$^\circ$'], 'interpreter', 'latex')
xlabel('Distance (cm)')
ylabel('Frequency (kHz)')
set(gca, 'xtick', disttick, 'xticklabel', disttick, 'fontsize', font_size)
set(gca, 'ytick', ftick_log, 'yticklabel', ftick_log, 'fontsize', font_size)
set(gca, 'xscale', 'log', 'yscale', 'log')
xlim(distlim)
ylim(flim)
%grid on
box on
set(gca, 'clim', clim_inter)
titstr = '$20\log_{10}|H_{\mbox{interaural}}|$';
hcb = colorbar('eastoutside', 'ytick', ctick_inter, 'xticklabel', ctick_inter, 'ylim', clim_inter);
set(get(hcb, 'ylabel'), 'String', titstr, 'interpreter', 'latex', 'fontsize', font_size);

