% Forum Acusticum 2023
% EEA summer school, Hot-topic 8
% Toolboxes for auditory modelling and psychoacoustic research. 
% 
% This script prototypes a sound localization model based on ITD while playing 
% around with the AMT functionalities. The code works and runs smoothly but it
% does NOT follow the AMT guidelines. Moreover, it does not leverage many of the 
% AMT functionalities like default parameters or caching. 
% As a policy, we ask authors interested in publishing their work in the AMT 
% to follow our folder's structure to increase code reproducibility. 
% 
% In this hands-on, the task is to re-organize this script's code to comply
% with the AMT guidelines. 
% 
% 
% Here, it follows the todo list: 
% - Divide this script into a model script (models/surname2023.m) and an
%       experiment file (experiements/exp_surname2023.m).
% - Move the code dedicated to plot into the plot folder
%       (plots/plot_surname2023.m).
% - Create an experiment flag in exp_surname2023 by using ltfatarghelper. 
%       After this addition you should call exp_surname2023('simple_test') to
%       reproduce the result of this script.
% - By using ltfatarghelper, define default parameters in you model. After
%       this the model interface should be 
%       varargout = surname2023(binaural_sound, fs, varargin).
% - Use amt_disp instead of disp or fprintf.
% - By using the caching feature in the AMT, check if the stimulus is
%       locally available by calling amt_cache('get', 'exp_surname2023')
%       otherwise generate and save it with amt_cache('set', 'exp_surname2023',
%       binaural_stimulus) for future use. 
% - Advanced: create a new experiment to be called with the code
% 		exp_surname2023('response_distribution') to visualize 
% 		the response distribution. Hint: given the same binaural sound, 
% 		estimate the angle with the model many times, 
% 		collect all responses and plot an histogram.
%
% We recommend to call exp_yourname2023.m to check if your code can reproduce
% the same results of this script.
%
% If you are lost in the process, you can check the code of
% workshop2022vienna to get an idea on how to implement each point. Another
% interesting pointer is the template model that provides an empty
% structure. You could also check other models. But be careful!
% The AMT has many models and not all are following correctly our guidelines. 
% That is why we provide the following webpage amtoolbox.org/models.php to
% report the quality of the current implementions. We recommend to
% take a look to the models with all three green ticks. 
%

%   #Author: Roberto Barumerli and Piotr Majdak (2023): initial implementation

% This file is licensed unter the GNU General Public License (GPL) either 
% version 3 of the license, or any later version as published by the Free Software 
% Foundation. Details of the GPLv3 can be found in the AMT directory "licences" and 
% at <https://www.gnu.org/licenses/gpl-3.0.html>. 
% You can redistribute this file and/or modify it under the terms of the GPLv3. 
% This file is distributed without any warranty; without even the implied warranty 
% of merchantability or fitness for a particular purpose. 

clearvars
close all

%% generate a binaural sound for our experiment
% set lateral direction at which the stimulus should be presented
lateral_direction_deg = -40; % deg

% load HRTF dataset
hrtfs = amt_load('amt_workshop', 'ARI_NH12_hrtf_M_dtf 256.sofa');
fs = hrtfs.Data.SamplingRate;

% generate sound of a source
sound = sig_whitenoiseburst(fs);

% spatialize sound 
% generate binaural sound by convolving HRIR with the sound source
% (we use SOFAspat from the SOFA toolbox)
binaural_sound = SOFAspat(sound, hrtfs, lateral_direction_deg, 0);

% Uncomment the following line to listen to it! 
% soundsc(.5*binaural_sound, fs);

%% Here, we call the model to estimate the lateral angle
% The model is wrapped into the function at the end of this script
itd_method = 'MaxIACCr'; % we want to use a specific method to compute the itd

% model start
% Some default values for our parameters
sound_speed = 340.29; % m/s
head_rad = 0.0875; % m
human_noise_std = 5; % degrees
spline_points = 100;

% We estimate the ITD using common/itdestimator.m

% Let's prepare the matrix because itdestimator is a bit picky.
% It requires the input matrix to have a matrix with (sounds x ear x time)
sound_temp(:,1,:) = binaural_sound(:,1)';
sound_temp(:,2,:) = binaural_sound(:,2)';

% Call itdestimator with few fixed parameters
perceived_itd = itdestimator(sound_temp, itd_method, 'lp', ...
    'upper_cutfreq', 3000, 'fs', fs);

% Then, we convert the itd value into an angle using 
% the Woodworth formula (Eq. 1.9 in Bosun Xie's book). 
% The formula is ITD(theta) = head_rad / sound_speed*(sin(theta) + theta) 
% with 0<=theta<=pi/2.
% We need to invert such equation to get to the angle theta and we solve our
% problem by using a simple interpolation method. 
az_support_rad = linspace(-pi/2, pi/2, spline_points);
response_rad = -spline(head_rad/sound_speed *(az_support_rad + sin(az_support_rad)), ...
                        az_support_rad, ...
                        perceived_itd);
% since the formula returns the angle in radiants we convert it degrees
% which is the common unit of measurement in the AMT for angles.
response_deg = rad2deg(response_rad);

% Finally, we want to humanize our model and listeners responses are 
% affected by uncertainty. We model this by adding some gaussian noise with
% zero mean. 
response_deg = response_deg + randn(1,1)*human_noise_std; 
% model end

%% let print the results of our experiment
fprintf("Perceived ITD: %.2f ms\n", perceived_itd/1e-3)
fprintf("True lateral direction: %.2f deg\n", lateral_direction_deg)
fprintf("Estimated lateral direction: %.2f deg\n", response_deg)

%% prepare the result
% here we do something fancy to visualize the actual and
% responded angles on an arc 

% we need some prep before plotting to generate the arc
% rotate to improve visualization 
targets = lateral_direction_deg + 90;
estimates = response_deg + 90; 

% generate arc for visualization
N = 100;
r_angl = linspace(0, 180, N);                             
radius = 1; 
% generate the cartesian coordinate matrix
[arc_x,arc_y] = sph2cart(deg2rad(r_angl), 0, 1);

% let's plot
figure
plot(arc_x, arc_y, 'k') 
hold on
plot([0 0], [0 1], '--', 'Color', [1 1 1]*0.2) 
scatter(0, 0, 1e3, [0 0 0])
plot([0 cosd(targets(1))], [0 sind(targets)+.1], 'r')
plot([0 cosd(estimates)], [0 sind(estimates)+.1], 'b')
axis([-1.25*radius  1.25*radius    0  1.25*radius]) 
set(gca, ...
    'XTick', [], ...
    'YTick', []);
legend({'', 'Front position', 'Subject head', ...
    'True location', 'Estimated location'}, 'Location', 'southoutside')
axis equal