% 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 .
% 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