THE AUDITORY MODELING TOOLBOX
This documentation page applies to an outdated major AMT version. We show it for archival purposes only.
Click here for the documentation menu and here to download the latest AMT (1.6.0).
Go to function
EXP_LINDEMANN1986 - Figures from Lindemann (1986)
Program code:
function output = exp_lindemann1986(varargin)
%EXP_LINDEMANN1986 Figures from Lindemann (1986)
% Usage: output = exp_lindemann1986(flag)
%
% EXP_LINDEMANN1986(flag) reproduces the results for the figure given
% by flag from the Lindemann (1986) paper. It will also plot the
% results. The format of its output depends on the chosen figure. If not
% otherwise stated, the cross-correlation of a pure tone sinusoids with
% f=500 Hz and different ITDs, ILDs or a combination of both is
% calculated. Because of the stationary character of the input signals,
% T_{int} = inf is used to produce only one time step in the crosscorr
% output.
%
% The following flags can be specified;
%
% 'plot' plot the output of the experiment. This is the default.
%
% 'noplot' Don't plot, only return data.
%
% 'auto' Re-calculate the file if it does not exist. Return 1 if the
% file exist, otherwise 0. This is the default
%
% 'redo' Always recalculate the results.
%
% 'cached' Always use the cached version. Default.
%
% 'fig6' Reproduce Fig.6 from Lindemann (1986). The cross-correlation is
% calculated for different ITDs and different inhibition factors
% c_s=0,0.3,1. Afterwards for every c_s the correlation is
% plotted for every used ITD dependend on the correlation-time
% delay. The output is cross-correlation result of the figure.
% The output has dimensions: number of c_s conditions x
% nitds x delay line length
%
% 'fig7' Reproduce Fig.7 from Lindemann (1986). The cross-correlation
% is calculated for different ITDs and different inhibition
% factors c_s=0,0.2,0.4,0.6,0.8,1.0. Afterwards for every
% c_s the displacement of the centroid of the auditory image
% is calculated and plotted depending on the ITD. The output is
% the displacement of the centroid for different c_s values.
% The output has dimensions: c_s x nitds
%
% 'fig8' Reproduce Fig.8 from Lindemann (1986). The cross-correlation
% is calculated for different ILDs and different inhibition factors
% c_s = 0.3, 1. Afterwards for every c_s the ILD is plotted
% depending on the correlation time. The output is the
% cross-correlation result of the figure. The dimensions of the
% output are: number of c_s conditions x nilds x delay line length.
%
% 'fig10' Reproduce Fig.10 from Lindemann (1986). The
% cross-correlation is calculated for different ILDs with an
% inhibition factor of c_s = 0.3 and a monaural detector
% factor w_f = 0.035. Afterwards the ILD is plotted depending
% on the correlation time. The output is the cross-correlation
% result of the figure. The output has dimensions:
% number of c_s conditions x nilds x delay line length
%
% 'fig11' Reproduce Fig.11 from Lindemann (1986). The centroid
% position is calculated for different ILDs, an inhibition
% factor c_s = 0.3 and a monaural detector factor w_f =
% 0.035. Afterwards for every c_s the displacement of the
% centroid is plotted dependend on the ILD. The output is the
% cross-correlation result of the to figure. The dimensions of
% the output are: number of c_s conditions x nilds x delay
% line length
%
% 'fig12' Reproduce Fig.12 from Lindemann (1986). The centroids are
% calculated for combinations of ITDs and ILDs. After the
% calculation the values for the centroids of the stimuli are
% searched to find the nearest value to 0. The corresponding ILD
% value is stored in output and plotted depending on the
% ITD. The output is the simulated results for a trading
% experiment. The ILD value for getting a centroid near the
% center for an combined ITD, ILD stimulus with a given ITD
% value. The output has dimensions: number of ITDs x 1
%
% 'fig13' Reproduce Fig.13 from Lindemann (1986). The centroids are
% calculated for ILD only and ITD/ILD combination
% stimuli. After the calculation the values for the centroids
% of the ILD only stimuli are searched to find the nearest
% value to the centroid of a given combined stimulus. The
% resulting ILD value is stored for different combinaition
% values and plotted dependend on ITD. The output is the ILD
% value for getting the same lateralization with an ILD only
% stimulus compared to a stimulus with both ITD and ILD.
% The output has dimensions: number of ILDs x number of ITDs
%
% 'fig14a' Reproduce fig.14 (a) from Lindemann (1986). The
% cross-correlations for a combination of ILDs and a ITD of
% ~1ms are calculated. This is done for different ILDs and
% different inhibition factors c_s = 0,0.2,0.4,0.6,1.
% Afterwards for every c_s the centroid of the auditory image
% is calculated and plotted dependend on the ILD. The output is
% the displacement of the centroid for different c_s values
% and ILDs. The output has dimensions: c_s x nilds
%
% 'fig14b' Reproduce Fig.14 (b) from Lindemann (1986). The
% cross-correlations for a combination of ILDs and a ITD of ~1ms
% are calculated. This is done for different small ILDs with a
% standard deviation of 0-5. Afterwards for every standard
% deviation the mean centroid displacement is calculated and
% plotted dependend on the ITD. The output is the displacement
% of the centroid as a function of the ITD averaged over
% different small ILD with a standard deviation of 0,1,2,3,4,5.
% The output has dimensions: nilds x nitds
% THIS FIGURE DOES NOT SEEM TO MATCH THE PAPER!
%
% 'fig15' Reproduce Fig.15 from Lindemann (1986). The cross-correlation
% for an ITD of -0.5ms is calculated. This is done for
% different ILDs, an inhibition factor c_s = 0.3 and a
% monaural detector factor w_f = 0.035. Afterwards for every
% c_s the ILD is plotted depending on the correlation
% time. The output is the cross-correlation result. The output
% has dimensions: number of c_s conditions x nilds x delay line
% length
%
% 'fig16' Reproduces Fig.16 from Lindemann (1986). The
% cross-correlations for combinations of ITDs and ILDs are
% calculated. Afterwards the combinations of ILD and ITD are
% looked for the ones that have two peaks in its
% cross-correlation which have the nearly same height. The
% corresponding ILD value is than stored in output and plotted
% dependend on the ITD. The output is the ILD value for which
% the cross-correlation for a combined ITD, ILD stimulus has
% two peaks with the same (nearest) height, depending on the
% ITD. The output has dimensions: number of ITDs x 1
%
% 'fig17' Reproduce Fig.17 from Lindemann (1986). The
% cross-correlation for ITD/ILD combination and ITD only
% stimuli is calculated. Afterwards the values for the
% centroids and maxima of the ITD only stimuli are searched to
% find the nearest value to the centroid and maxima of a given
% combined stimulus. The resulting ITD value is stored for
% different combinaition values. The output is the ITD
% value for getting the same lateralization with an ITD only
% stimulus compared to a stimulus with both ITD and ILD using
% the centroid of the cross-correlation. The output has
% dimensions: number of ITDs x number of ILDs.
%
% 'fig18' Reproduce Fig.18 from Lindemann (1986). The cross-correlation
% of pink noise with different interaural coherence values is
% calculated. Afterwards for every interaural coherence value
% the correlation is plotted dependend on the correlation-time
% delay. The output is the cross-correlation result of the
% figure. The output has dimensions: number of interaural
% coherences x delay line length
%
% If no flag is given, the function will print the list of valid flags.
%
% Examples:
% ---------
%
% To display Figure 6 use :
%
% exp_lindemann1986('fig6');
%
% To display Figure 7 use :
%
% exp_lindemann1986('fig7');
%
% To display Figure 8 use :
%
% exp_lindemann1986('fig8');
%
% To display Figure 10 use :
%
% exp_lindemann1986('fig10');
%
% To display Figure 11 use :
%
% exp_lindemann1986('fig11');
%
% To display Figure 12 use :
%
% exp_lindemann1986('fig12');
%
% To display Figure 13 use :
%
% exp_lindemann1986('fig13');
%
% To display Figure 14a use :
%
% exp_lindemann1986('fig14a');
%
% To display Figure 15 use :
%
% exp_lindemann1986('fig15');
%
% To display Figure 16 use :
%
% exp_lindemann1986('fig16');
%
% To display Figure 17 use :
%
% exp_lindemann1986('fig17');
%
% To display Figure 18 use :
%
% exp_lindemann1986('fig18');
%
% References:
% W. Lindemann. Extension of a binaural cross-correlation model by
% contralateral inhibition. I. Simulation of lateralization for
% stationary signals. J. Acoust. Soc. Am., 80:1608--1622, 1986.
%
%
%
% Url: http://amtoolbox.sourceforge.net/amt-0.10.0/doc/experiments/exp_lindemann1986.php
% Copyright (C) 2009-2020 Piotr Majdak and the AMT team.
% This file is part of Auditory Modeling Toolbox (AMT) version 0.10.0
%
% This program is free software: you can redistribute it and/or modify
% it under the terms of the GNU General Public License as published by
% the Free Software Foundation, either version 3 of the License, or
% (at your option) any later version.
%
% This program is distributed in the hope that it will be useful,
% but WITHOUT ANY WARRANTY; without even the implied warranty of
% MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
% GNU General Public License for more details.
%
% You should have received a copy of the GNU General Public License
% along with this program. If not, see <http://www.gnu.org/licenses/>.
% AUTHOR: Hagen Wierstorf
definput.import={'amt_cache'};
definput.flags.type={'missingflag','fig6','fig7','fig8','fig10','fig11',...
'fig12','fig13','fig14a','fig14b','fig15',...
'fig16','fig17','fig18'};
definput.flags.plot={'plot','noplot'};
[flags,keyvals] = ltfatarghelper({},definput,varargin);
if flags.do_missingflag
flagnames=[sprintf('%s, ',definput.flags.type{2:end-2}),...
sprintf('%s or %s',definput.flags.type{end-1},definput.flags.type{end})];
error('%s: You must specify one of the following flags: %s.',upper(mfilename),flagnames);
end;
%% ------ FIG 6 -----------------------------------------------------------
if flags.do_fig6
% Sampling rate
fs = 44100;
% Frequency of the sinusoid
f = 500;
T = 1/f;
fc = round(freqtoerb(f)); % corresponding frequency channel
% Model parameter
T_int = inf;
w_f = 0;
M_f = 6; % not used, if w_f==0
c_s = [0,0.3,1];
% NOTE: the longer the signal, the more time we need for computation. On the
% other side N_1 needs to be long enough to eliminate any onset effects.
% Lindemann uses N_1 = 17640. Here I uses only N_1 = 2205 which gives the same
% results for this demo.
N_1 = ceil(25*T*fs);
siglen = ceil(30*T*fs);
% Calculate crosscorrelations for 21 ITD points between 0~ms and 1~ms
nitds = 21; % number of used ITDs
ndl = round(fs/1000)+1; % length of the delay line (see lindemann1986_bincorr.m)
itd = linspace(0,1,nitds);
output=amt_cache('get','fig6',flags.cachemode);
if isempty(output)
output = zeros(length(c_s),nitds,ndl);
for ii = 1:nitds;
% Generate ITD shifted sinusoid
sig = sig_itdsin(f,itd(ii),fs);
% Use only the beginning of the signal to generate only one time instance of
% the cross-correlation and apply onset window
sig = sig(1:siglen,:);
sig = rampsignal(sig,[round(N_1/2)-1 0],'tria');
% Calculate cross-correlation for different inhibition factor c_s
for jj = 1:length(c_s)
% Calculate cross-correlation (and squeeze due to T_int==inf)
tmp = squeeze(lindemann1986(sig,fs,c_s(jj),w_f,M_f,T_int,N_1));
% Store the needed frequency channel. NOTE: the cross-correlation
% calculation starts with channel 5, so we have to subtract 4.
output(jj,ii,:) = tmp(:,fc-4);
end
end
amt_cache('set','fig6',output);
end;
if flags.do_plot
% ------ Plotting ------
% Generate time axis
tau = linspace(-1,1,ndl);
% Plot figure for every c_s condition
for jj = 1:length(c_s)
figure;
mesh(tau,itd(end:-1:1),squeeze(output(jj,:,:)));
view(0,57);
xlabel('correlation-time tau (ms)');
ylabel('interaural time difference (ms)');
set(gca,'YTick',0:0.2:1);
set(gca,'YTickLabel',{'1','0.8','0.6','0.4','0.2','0'});
tstr = sprintf('c_s = %.1f\nw_f = 0\nf = %i Hz\n',c_s(jj),f);
title(tstr);
end
end;
end;
%% ------ FIG 7 -----------------------------------------------------------
if flags.do_fig7
% Sampling rate
fs = 44100;
% Frequency of the sinusoid
f = 500;
T = 1/f;
fc = round(freqtoerb(f)); % corresponding frequency channel
% Model parameter
T_int = inf;
w_f = 0;
M_f = 6; % not used, if w_f==0
c_s = 0:0.2:1;
% NOTE: the longer the signal, the more time we need for computation. On the
% other side N_1 needs to be long enough to eliminate any onset effects.
% Lindemann uses N_1 = 17640. Here I uses only N_1 = 2205 which gives the same
% results for this demo.
N_1 = ceil(25*T*fs);
siglen = ceil(30*T*fs);
% Calculate crosscorrelations for 21 ITD points between 0~ms and 1~ms
nitds = 21; % number of used ITDs
itd = linspace(0,1,nitds);
output=amt_cache('get','fig7',flags.cachemode);
if isempty(output)
output = zeros(length(c_s),nitds);
for ii = 1:nitds
% Generate ITD shifted sinusoid
sig = sig_itdsin(f,itd(ii),fs);
% Use only the beginning of the signal to generate only one time instance of
% the cross-correlation and apply an onset window
sig = sig(1:siglen,:);
sig = rampsignal(sig,[round(N_1/2)-1 0],'tria');
% Calculate cross-correlation for different inhibition factor c_s
for jj = 1:length(c_s)
% Calculate cross-correlation (and squeeze due to T_int==inf)
tmp = squeeze(lindemann1986(sig,fs,c_s(jj),w_f,M_f,T_int,N_1));
% Store the needed frequency channel. NOTE: the cross-correlation
% calculation starts with channel 5, so we have to subtract 4.
cc = tmp(:,fc-4);
% Calculate the position of the centroid
output(jj,ii) = lindemann1986_centroid(cc);
end
end
amt_cache('set','fig7',output);
end;
if flags.do_plot
% ------ Plotting ------
figure;
for jj = 1:length(c_s)
plot(itd,output(jj,:));
hold on;
end
xlabel('interaural time difference (ms)');
ylabel('displacement of the centroid d');
tstr = sprintf('w_f = 0\nf = %i Hz\n',f);
title(tstr);
end;
end;
%% ------ FIG 8 -----------------------------------------------------------
if flags.do_fig8
% Sampling rate
fs = 44100;
% Frequency of the sinusoid
f = 500;
T = 1/f;
fc = round(freqtoerb(f)); % corresponding frequency channel
% Model parameter
T_int = inf;
w_f = 0;
M_f = 6; % not used, if w_f==0
c_s = [0.3,1];
% NOTE: the longer the signal, the more time we need for computation. On the
% other side N_1 needs to be long enough to eliminate any onset effects.
% Lindemann uses N_1 = 17640. Here I uses only N_1 = 2205 which gives the same
% results for this demo.
N_1 = ceil(25*T*fs);
siglen = ceil(30*T*fs);
% Calculate crosscorrelations for 26 ILD points between 0~dB and 25~dB
nilds = 26; % number of used ILDs
ndl = 2*round(fs/2000)+1; % length of the delay line (see bincorr.m)
ild = linspace(0,25,nilds);
output=amt_cache('get','fig8',flags.cachemode);
if isempty(output)
output = zeros(2,nilds,ndl);
for ii = 1:nilds
% Generate sinusoid with given ILD
sig = sig_ildsin(f,ild(ii),fs);
% Use only the beginning of the signal to generate only one time instance of
% the cross-correlation and apply onset window
sig = sig(1:siglen,:);
sig = rampsignal(sig,[round(N_1/2)-1 0],'tria');
% Calculate cross-correlation for different inhibition factor c_s
for jj = 1:length(c_s)
% Calculate cross-correlation (and squeeze due to T_int==inf)
tmp = squeeze(lindemann1986(sig,fs,c_s(jj),w_f,M_f,T_int,N_1));
% Store the needed frequency channel. NOTE: the cross-correlation
% calculation starts with channel 5, so we have to subtract 4.
output(jj,ii,:) = tmp(:,fc-4);
end
end
amt_cache('set','fig8',output);
end;
if flags.do_plot
% ------ Plotting ------
% Generate time axis
tau = linspace(-1,1,ndl);
% Plot figure for every c_s condition
for jj = 1:length(c_s)
figure;
mesh(tau,ild(end:-1:1),squeeze(output(jj,:,:)));
view(0,57);
xlabel('correlation-time delay (ms)');
ylabel('interaural level difference (dB)');
set(gca,'YTick',0:5:25);
set(gca,'YTickLabel',{'25','20','15','10','5','0'});
tstr = sprintf('c_{inh} = %.1f\nw_f = 0\nf = %i Hz\n',c_s(jj),f);
title(tstr);
end
end;
end;
%% ------ FIG 10 ----------------------------------------------------------
if flags.do_fig10
% Sampling rate
fs = 44100;
% Frequency of the sinusoid
f = 500;
T = 1/f;
fc = round(freqtoerb(f)); % corresponding frequency channel
% Model parameter
T_int = inf;
w_f = 0.035;
M_f = 6; % not used, if w_f==0
c_s = 0.3;
% NOTE: the longer the signal, the more time we need for computation. On the
% other side N_1 needs to be long enough to eliminate any onset effects.
% Lindemann uses N_1 = 17640. Here I uses only N_1 = 2205 which gives the same
% results for this demo.
N_1 = ceil(25*T*fs);
siglen = ceil(30*T*fs);
% Calculate crosscorrelations for 26 ILD points between 0~dB and 25~dB
nilds = 26; % number of used ILDs
ndl = 2*round(fs/2000)+1; % length of the delay line (see bincorr.m)
ild = linspace(0,25,nilds);
output=amt_cache('get','fig10',flags.cachemode);
if isempty(output)
output = zeros(length(c_s),nilds,ndl);
for ii = 1:nilds
% Generate sinusoid with given ILD
sig = sig_ildsin(f,ild(ii),fs);
% Use only the beginning of the signal to generate only one time instance of
% the cross-correlation and apply onset window
sig = sig(1:siglen,:);
sig = rampsignal(sig,[round(N_1/2)-1 0],'tria');
% Calculate cross-correlation for different inhibition factor c_s
for jj = 1:length(c_s)
% Calculate cross-correlation (and squeeze due to T_int==inf)
tmp = squeeze(lindemann1986(sig,fs,c_s(jj),w_f,M_f,T_int,N_1));
% Store the needed frequency channel. NOTE: the cross-correlation
% calculation starts with channel 5, so we have to subtract 4.
output(jj,ii,:) = tmp(:,fc-4);
end
end
amt_cache('set','fig10',output);
end;
% ------ Plotting ------
if flags.do_plot
% Generate time axis
tau = linspace(-1,1,ndl);
% Plot figure for every c_s condition
for jj = 1:length(c_s)
figure;
mesh(tau,ild(end:-1:1),squeeze(output(jj,:,:)));
view(0,57);
xlabel('correlation-time delay (ms)');
ylabel('interaural level difference (dB)');
set(gca,'YTick',0:5:25);
set(gca,'YTickLabel',{'25','20','15','10','5','0'});
tstr = sprintf('c_{inh} = %.1f\nw_f = 0.035\nf = %i Hz\n',c_s(jj),f);
title(tstr);
end
end;
end;
%% ------ FIG 11 ----------------------------------------------------------
if flags.do_fig11
% Sampling rate
fs = 44100;
% Frequency of the sinusoid
f = 500;
T = 1/f;
fc = round(freqtoerb(f)); % corresponding frequency channel
% Model parameter
T_int = inf;
w_f = [0,0,0.035];
M_f = 6; % not used, if w_f==0
c_s = [0.3,1,0.3];
% NOTE: the longer the signal, the more time we need for computation. On the
% other side N_1 needs to be long enough to eliminate any onset effects.
% Lindemann uses N_1 = 17640. Here I uses only N_1 = 2205 which gives the same
% results for this demo.
N_1 = ceil(25*T*fs);
siglen = ceil(30*T*fs);
% Caelculate crosscorrelations for 26 ILD points between 0~dB and 25~dB
nilds = 26; % number of used ILDs
ild = linspace(0,25,nilds);
output=amt_cache('get','fig11',flags.cachemode);
if isempty(output)
output = zeros(length(c_s),nilds);
for ii = 1:nilds
% Generate sinusoid with given ILD
sig = sig_ildsin(f,ild(ii),fs);
% Use only the beginning of the signal to generate only one time instance of
% the cross-correlation and apply onset window
sig = sig(1:siglen,:);
sig = rampsignal(sig,[round(N_1/2)-1 0],'tria');
% Calculate cross-correlation for different inhibition factor c_s
for jj = 1:length(c_s)
% Calculate cross-correlation (and squeeze due to T_int==inf)
tmp = squeeze(lindemann1986(sig,fs,c_s(jj),w_f(jj),M_f,T_int,N_1));
% Store the needed frequency channel
cc = tmp(:,fc-4);
% Calculate the position of the centroid
output(jj,ii) = lindemann1986_centroid(cc);
end
end
amt_cache('set','fig11',output);
end;
if flags.do_plot
% ------ Plotting ------
figure;
% Plot data from experiments
data = data_lindemann1986('fig11_yost');
plot(data(:,1),data(:,2)*0.058,'+');
hold on;
data = data_lindemann1986('fig11_sayers');
plot(data(:,1),data(:,2)*0.058,'*');
% Plot line for every condition
for jj = 1:length(c_s)
plot(ild,output(jj,:));
end
axis([0 25 0 0.8]);
legend('500 Hz, Yost','600 Hz, Sayers','500 Hz, model');
xlabel('interaural level difference (dB)');
ylabel('displacement of the centroid d');
end;
end;
%% ------ FIG 12 ----------------------------------------------------------
if flags.do_fig12
% Sampling rate
fs = 44100;
% Frequency of the sinusoid
f = 500;
T = 1/f;
fc = round(freqtoerb(f)); % corresponding frequency channel
% Model parameter
T_int = inf;
w_f = 0.035;
M_f = 6; % not used, if w_f==0
c_s = 0.3;
% NOTE: the longer the signal, the more time we need for computation. On the
% other side N_1 needs to be long enough to eliminate any onset effects.
% Lindemann uses N_1 = 17640. Here I uses only N_1 = 2205 which gives the same
% results for this demo.
N_1 = ceil(25*T*fs);
siglen = ceil(30*T*fs);
% Calculate crosscorrelations for 26 ILD points between 0~dB and 25~dB
nilds = 21; % number of used ILDs for the ILD only stimuli
nitds = 21; % number of used ITDs for the combined stimuli
ild = linspace(0,10,nilds);
itd = linspace(-1,0,nitds);
output=amt_cache('get','fig12',flags.cachemode);
if isempty(output)
% Calculate the centroids for ILD+ITD stimuli
cen = zeros(nitds,nilds);
for ii = 1:nitds
for jj = 1:nilds
% Generate sinusoid with given ILD
sig = sig_itdildsin(f,itd(ii),ild(jj),fs);
% Use only the beginning of the signal to generate only one time
% instance of the cross-correlation and apply a linear onset window
sig = sig(1:siglen,:);
sig = rampsignal(sig,[round(N_1/2)-1 0],'tria');
% Calculate cross-correlation (and squeeze due to T_int==inf)
tmp = squeeze(lindemann1986(sig,fs,c_s,w_f,M_f,T_int,N_1));
% Store the needed frequency channel
cc = tmp(:,fc-4);
% Calculate the position of the centroid
cen(ii,jj) = lindemann1986_centroid(cc);
end
end
% ------ Fiting ------
% For every ITD find the ILD with gives a centroid near 0
output = zeros(nitds,1);
for ii = 1:nitds
% Find centroid nearest to 0
[val,idx] = min(abs(cen(ii,:)));
output(ii) = ild(idx);
end
amt_cache('set','fig12',output);
end;
if flags.do_plot
% ------ Plotting ------
figure;
data = data_lindemann1986('fig12_400');
plot(data(:,1),data(:,2),'+');
hold on;
data = data_lindemann1986('fig12_600');
plot(data(:,1),data(:,2),'*');
plot(itd,output);
axis([-1 0 0 10]);
legend('400 Hz','600 Hz','500 Hz, model');
xlabel('interaural time difference (ms)');
ylabel('interaural level difference (dB)');
end;
end;
%% ------ FIG 13 ----------------------------------------------------------
if flags.do_fig13
% Sampling rate
fs = 44100;
% Frequency of the sinusoid
f = 500;
T = 1/f;
fc = round(freqtoerb(f)); % corresponding frequency channel
% Model parameter
T_int = inf;
N_1 = 1;
w_f = 0.035;
M_f = 6; % not used, if w_f==0
c_s = 0.3;
% NOTE: the longer the signal, the more time we need for computation. On the
% other side N_1 needs to be long enough to eliminate any onset effects.
% Lindemann uses N_1 = 17640. Here I uses only N_1 = 2205 which gives the same
% results for this demo.
N_1 = ceil(25*T*fs);
siglen = ceil(30*T*fs);
% Calculate crosscorrelations for 26 ILD points between 0~dB and 25~dB
nilds_p = 21; % number of used ILDs for the ILD only stimuli
nitds_t = 21; % number of used ITDs for the combined stimuli
nilds_t = 6; % number of used ILDs for the combined stimuli
ndl = 2*round(fs/2000)+1; % length of the delay line (see bincorr.m)
ild_p = linspace(-10,40,nilds_p);
itd_t = linspace(-1,1,nitds_t);
ild_t = [-3,0,3,9,15,25];
output=amt_cache('get','fig13',flags.cachemode);
if isempty(output)
% Calculate the centroids for the ILD only stimuli
cen_p = zeros(1,nilds_p);
for ii = 1:nilds_p
% Generate sinusoid with given ILD
sig = sig_ildsin(f,ild_p(ii),fs);
% Use only the beginning of the signal to generate only one time instance of
% the cross-correlation and apply a linear onset window
sig = sig(1:siglen,:);
sig = rampsignal(sig,[round(N_1/2)-1 0],'tria');
% Calculate cross-correlation (and squeeze due to T_int==inf)
tmp = squeeze(lindemann1986(sig,fs,c_s,w_f,M_f,T_int,N_1));
% Store the needed frequency channel
cc = tmp(:,fc-4);
% Calculate the position of the centroid
cen_p(ii) = lindemann1986_centroid(cc);
end
% Calculate the centroids for the combined stimuli
cen_t = zeros(nilds_t,nitds_t);
for ii = 1:nitds_t
for jj = 1:nilds_t
sig = sig_itdildsin(f,itd_t(ii),ild_t(jj),fs);
sig = sig(1:siglen,:);
sig = rampsignal(sig,[round(N_1/2)-1 0],'tria');
tmp = squeeze(lindemann1986(sig,fs,c_s,w_f,M_f,T_int,N_1));
cc = tmp(:,fc-4);
cen_t(jj,ii) = lindemann1986_centroid(cc);
end
end
% ------ Fitting ------
% For the results for the combined centroids find the nearest centroid for the
% ILD only stimuli
output = zeros(nilds_t,nitds_t);
for ii = 1:nitds_t
for jj = 1:nilds_t
idx = findnearest(cen_p,cen_t(jj,ii));
output(jj,ii) = ild_p(idx);
end
end
amt_cache('set','fig13',output);
end;
if flags.do_plot
% ------ Plotting ------
% First plot the only data from experiments
figure;
d = data_lindemann1986('fig13');
plot(d(:,1),d(:,2),'x-r', ... % -3dB
d(:,1),d(:,3),'x-b', ... % 3dB
d(:,1),d(:,4),'x-g', ... % 9dB
d(:,1),d(:,5),'x-b', ... % 15dB
d(:,1),d(:,6),'x-r') % 25dB
legend('25dB','15dB','9dB','3dB','-3dB');
axis([-1 1 -10 40]);
set(gca,'XTick',-1:0.4:1);
xlabel('interaural time difference (ms)');
ylabel('interaural level difference (dB)');
% Then plot model results
figure;
% Plot line for every condition
for jj = 1:nilds_t
plot(itd_t,output(jj,:));
hold on;
end
axis([-1 1 -10 40]);
set(gca,'XTick',-1:0.4:1);
xlabel('interaural time difference (ms)');
ylabel('interaural level difference (dB)');
end;
end;
%% ------ FIG 14a ---------------------------------------------------------
if flags.do_fig14a
% Sampling rate
fs = 44100;
% Frequency of the sinusoid
f = 500;
T = 1/f;
fc = round(freqtoerb(f)); % corresponding frequency channel
% --- FIG 14 (a) ---
% Model parameter
T_int = inf;
w_f = 0;
M_f = 6; % not used, if w_f==0
c_s = [0.0,0.2,0.4,0.6,1.0];
% NOTE: the longer the signal, the more time we need for computation. On the
% other side N_1 needs to be long enough to eliminate any onset effects.
% Lindemann uses N_1 = 17640. Here I uses only N_1 = 2205 which gives the same
% results for this demo.
N_1 = ceil(25*T*fs);
siglen = ceil(30*T*fs);
% Calculate crosscorrelations for 21 ITD points between 0~ms and 1~ms
nilds = 21; % number of used ITDs
ild = linspace(-3,3,nilds);
itd = 2000/f;
output=amt_cache('get','fig14a',flags.cachemode);
if isempty(output)
output = zeros(length(c_s),nilds);
for ii = 1:nilds
% Generate ITD shifted sinusoid
sig = sig_itdildsin(f,itd,ild(ii),fs);
% Use only the beginning of the signal to generate only one time instance of
% the cross-correlation
sig = sig(1:siglen,:);
% Apply a linear onset window with length N_1/2 to minimize onset effects
% (see lindemann1986 p. 1614)
sig = rampsignal(sig,[round(N_1/2)-1 0],'tria');
% Calculate cross-correlation for different inhibition factor c_s
for jj = 1:length(c_s)
% Calculate cross-correlation (and squeeze due to T_int==inf)
tmp = squeeze(lindemann1986(sig,fs,c_s(jj),w_f,M_f,T_int,N_1));
% Store the needed frequency channel. NOTE: the cross-correlation
% calculation starts with channel 5, so we have to subtract 4.
cc = tmp(:,fc-4);
% Calculate the position of the centroid
output(jj,ii) = lindemann1986_centroid(cc);
end
end
amt_cache('set','fig14a',output);
end;
if flags.do_plot
% ------ Plotting ------
figure;
for jj = 1:length(c_s)
plot(ild,output(jj,:));
hold on;
end
xlabel('interaural level difference (dB)');
ylabel('displacement of the centroid d');
tstr = sprintf('w_f = 0\nf = %i Hz\n',f);
title(tstr);
end;
end;
%% ------ FIG 14b ---------------------------------------------------------
if flags.do_fig14b
% Sampling rate
fs = 44100;
% Frequency of the sinusoid
f = 500;
T = 1/f;
fc = round(freqtoerb(f)); % corresponding frequency channel
% Model parameter
T_int = inf;
w_f = 0;
M_f = 6; % not used, if w_f==0
c_s = 1.0;
% NOTE: the longer the signal, the more time we need for computation. On the
% other side N_1 needs to be long enough to eliminate any onset effects.
% Lindemann uses N_1 = 17640. Here I uses only N_1 = 2205 which gives the same
% results for this demo.
N_1 = ceil(25*T*fs);
siglen = ceil(30*T*fs);
% Calculate crosscorrelations for 21 ITD points between 0~ms and 1~ms
nitds = 21; % number of used ITDS
nilds = 201; % number of used ILDs
itd = linspace(0,1,nitds);
ild_std = [1,2,3,4,5];
output=amt_cache('get','fig14b',flags.cachemode);
if isempty(output)
amt_disp('NOTE: this test function will need a lot of time!','progress');
output = zeros(length(ild_std)+1,nitds);
centmp = zeros(nilds,nitds);
for ii = 1:nitds
% Show progress
amt_disp([num2str(ii) ' of ' num2str(nitds)],'progress');
% First generate the result for std(ILD) == 0
sig = sig_itdsin(f,itd(ii),fs);
sig = sig(1:siglen,:);
sig = rampsignal(sig,[round(N_1/2)-1 0],'tria');
tmp = squeeze(lindemann1986(sig,fs,c_s,w_f,M_f,T_int,N_1));
cc = tmp(:,fc-4);
cen(1,ii) = lindemann1986_centroid(cc);
% Generate results for std(ILD) ~= 0
for nn = 1:length(ild_std)
% Generate normal distributed ILDs with mean ~ 0 and std = 1
tmp = randn(nilds,1);
tmp = tmp/std(tmp);
% Generate desired ILD distribution
ild = tmp * ild_std(nn);
% For all distributed ILD values calculate the centroid
for jj = 1:nilds
sig = sig_itdildsin(f,itd(ii),ild(jj),fs);
sig = sig(1:siglen,:);
sig = rampsignal(sig,[round(N_1/2)-1 0],'tria');
tmp = squeeze(lindemann1986(sig,fs,c_s,w_f,M_f,T_int,N_1));
cc = tmp(:,fc-4);
centmp(jj,ii) = lindemann1986_centroid(cc);
end
% Calculate the mean centroid above the ILD distribution
output(nn+1,ii) = mean(centmp(:,ii));
end
end
amt_cache('set','fig14b',output);
end;
if flags.do_plot
figure;
for nn = 1:length(ild_std)+1
plot(itd,output(nn,:)); hold on;
end
axis([0 1 0 1]);
xlabel('interaural level difference (dB)');
ylabel('displacement of the centroid d');
tstr = sprintf('w_f = 0\nc_{inh} = 1\nf = %i Hz\n',f);
title(tstr);
end;
end;
%% ------ FIG 15 ----------------------------------------------------------
if flags.do_fig15
% Sampling rate
fs = 44100;
% Frequency of the sinusoid
f = 500;
T = 1/f;
fc = round(freqtoerb(f)); % corresponding frequency channel
% Model parameter
T_int = inf;
w_f = 0.035;
M_f = 6; % not used, if w_f==0
c_s = 0.3;
% NOTE: the longer the signal, the more time we need for computation. On the
% other side N_1 needs to be long enough to eliminate any onset effects.
% Lindemann uses N_1 = 17640 (200*T). Here I uses only N_1 = 2205 (25*T)
% which gives the same results for this demo.
N_1 = ceil(25*T*fs);
siglen = ceil(30*T*fs);
% Calculate crosscorrelations for 26 ILD points between 0~dB and 25~dB
nilds = 26; % number of used ILDs
ndl = 2*round(fs/2000)+1; % length of the delay line (see bincorr.m)
ild = linspace(0,25,nilds);
itd = -0.5;
output=amt_cache('get','fig15',flags.cachemode);
if isempty(output)
output = zeros(length(c_s),nilds,ndl);
for ii = 1:nilds
% Generate sinusoid with given ILD
sig = sig_itdildsin(f,itd,ild(ii),fs);
% Use only the beginning of the signal to generate only one time instance of
% the cross-correlation and apply onset window
sig = sig(1:siglen,:);
sig = rampsignal(sig,[round(N_1/2)-1 0],'tria');
% Calculate cross-correlation for different inhibition factor c_s
for jj = 1:length(c_s)
% Calculate cross-correlation (and squeeze due to T_int==inf)
tmp = squeeze(lindemann1986(sig,fs,c_s(jj),w_f,M_f,T_int,N_1));
% Store the needed frequency channel. NOTE: the cross-correlation
% calculation starts with channel 5, so we have to subtract 4.
output(jj,ii,:) = tmp(:,fc-4);
end
end
amt_cache('set','fig15',output);
end;
if flags.do_plot
% ------ Plottinqg --------------------------------------------------------
% Generate time axis
tau = linspace(-1,1,ndl);
% Plot figure for every c_s condition
for jj = 1:length(c_s)
figure;
mesh(tau,ild(end:-1:1),squeeze(output(jj,:,:)));
view(0,57);
xlabel('correlation-time delay (ms)');
ylabel('interaural level difference (dB)');
set(gca,'YTick',0:5:25);
set(gca,'YTickLabel',{'25','20','15','10','5','0'});
tstr = sprintf('c_{inh} = %.1f\nw_f = 0.035\nf = %i Hz\n',c_s(jj),f);
title(tstr);
end
end;
end;
%% ------ FIG 16 ---------------------------------------------------------
if flags.do_fig16
s = [mfilename('fullpath'),'_fig16.mat'];
% Sampling rate
fs = 44100;
% Frequency of the sinusoid
f = 500;
T = 1/f;
fc = round(freqtoerb(f)); % corresponding frequency channel
% Model parameter
T_int = inf;
w_f = 0.035;
M_f = 6; % not used, if w_f==0
c_s = 0.3;
% NOTE: the longer the signal, the more time we need for computation. On the
% other side N_1 needs to be long enough to eliminate any onset effects.
% Lindemann uses N_1 = 17640. Here I uses only N_1 = 2205 which gives the same
% results for this demo.
N_1 = ceil(25*T*fs);
siglen = ceil(30*T*fs);
% Calculate crosscorrelations for 26 ILD points between 0~dB and 25~dB
nilds = 21; % number of used ILDs for the ILD only stimuli
nitds = 21; % number of used ITDs for the combined stimuli
ndl = 2*round(fs/2000)+1; % length of the delay line (see bincorr.m)
ild = linspace(0,10,nilds);
itd = linspace(-1,0,nitds);
output=amt_cache('get','fig16',flags.cachemode);
if isempty(output)
% Calculate the centroids for ILD+ITD stimuli
cc = zeros(nitds,nilds,ndl);
for ii = 1:nitds
for jj = 1:nilds
% Generate sinusoid with given ILD
sig = sig_itdildsin(f,itd(ii),ild(jj),fs);
% Use only the beginning of the signal to generate only one time
% instance of the cross-correlation and apply a linear onset window
sig = sig(1:siglen,:);
sig = rampsignal(sig,[round(N_1/2)-1 0],'tria');
% Calculate cross-correlation (and squeeze due to T_int==inf)
tmp = squeeze(lindemann1986(sig,fs,c_s,w_f,M_f,T_int,N_1));
% Store the needed frequency channel
cc(ii,jj,:) = tmp(:,fc-4);
end
end
% ------ Fiting ------
% For every ITD find the ILD with gives a centroid near 0
output = zeros(nitds,1);
for ii = 1:nitds
% Find two peaks of same height
idx = findpeaks(squeeze(cc(ii,:,:)));
output(ii) = ild(idx);
end
amt_cache('set','fig16',output);
end;
if flags.do_plot
% ------ Plotting ------
figure;
data = data_lindemann1986('fig16');
plot(data(:,1),data(:,2),'+');
hold on;
plot(itd,output);
axis([-1 0 0 15]);
xlabel('interaural time difference (ms)');
ylabel('interaural level difference (dB)');
end;
end;
%% ------ FIG 17 ----------------------------------------------------------
if flags.do_fig17
% Sampling rate
fs = 44100;
% Frequency of the sinusoid
f = 500;
T = 1/f;
fc = round(freqtoerb(f)); % corresponding frequency channel
% Model parameter
T_int = inf;
w_f = 0.035;
M_f = 6; % not used, if w_f==0
c_s = 0.3;
% NOTE: the longer the signal, the more time we need for computation. On the
% other side N_1 needs to be long enough to eliminate any onset effects.
% Lindemann uses N_1 = 17640. Here I uses only N_1 = 2205 which gives the same
% results for this demo.
N_1 = ceil(25*T*fs);
siglen = ceil(30*T*fs);
% Calculate crosscorrelations for 26 ILD points between 0~dB and 25~dB
nitds_p = 21; % number of used ITDs for the ITD only stimuli
nitds_t = 4; % number of used ITDs for the combined stimuli
nilds_t = 21; % number of used ILDs for the combined stimuli
itd_p = linspace(-0.36,0.72,nitds_p);
ild_t = linspace(-9,9,nilds_t);
itd_t = [0,0.09,0.18,0.27];
output=amt_cache('get','fig17',flags.cachemode);
if isempty(output)
% Calculate the centroids for the ITD only stimuli
cen_p = zeros(1,nitds_p);
max_p = zeros(1,nitds_p);
for ii = 1:nitds_p
% Generate sinusoid with given ILD
sig = sig_itdsin(f,itd_p(ii),fs);
% Use only the beginning of the signal to generate only one time instance of
% the cross-correlation and apply a linear onset window
sig = sig(1:siglen,:);
sig = rampsignal(sig,[round(N_1/2)-1 0],'tria');
% Calculate cross-correlation (and squeeze due to T_int==inf)
tmp = squeeze(lindemann1986(sig,fs,c_s,w_f,M_f,T_int,N_1));
% Store the needed frequency channel
cc = tmp(:,fc-4);
% Find the maximum position
max_p(ii) = findmax(cc);
% Calculate the position of the centroid
cen_p(ii) = lindemann1986_centroid(cc);
end
% Calculate the centroids for the combined stimuli
cen_t = zeros(nitds_t,nilds_t);
max_t = zeros(nitds_t,nilds_t);
for ii = 1:nilds_t
for jj = 1:nitds_t
sig = sig_itdildsin(f,itd_t(jj),ild_t(ii),fs);
sig = sig(1:siglen,:);
sig = rampsignal(sig,[round(N_1/2)-1 0],'tria');
tmp = squeeze(lindemann1986(sig,fs,c_s,w_f,M_f,T_int,N_1));
cc = tmp(:,fc-4);
max_t(jj,ii) = findmax(cc);
cen_t(jj,ii) = lindemann1986_centroid(cc);
end
end
% ------ Fiting ------
% For the results for the combined centroids find the nearest centroid for the
% ITD only stimuli
output.rmax = zeros(nitds_t,nilds_t);
output.rcen = output.rmax;
for ii = 1:nilds_t
for jj = 1:nitds_t
idx = findnearest(cen_p,cen_t(jj,ii));
output.rcen(jj,ii) = itd_p(idx);
idx = findnearest(max_p,max_t(jj,ii));
output.rmax(jj,ii) = itd_p(idx);
end
end
amt_cache('set','fig17',output);
end;
if flags.do_plot
% ------ Plotting ------
% First plot the only data from experiments
figure; % fig 17 (a)
d = data_lindemann1986('fig17');
plot(d(:,1),d(:,5),'x-r', ... % 0 ms
d(:,1),d(:,4),'x-b', ... % 0.09 ms
d(:,1),d(:,3),'x-g', ... % 0.18 ms
d(:,1),d(:,2),'x-b') % 0.27 ms
legend('0.27 ms','0.18 ms','0.09 ms','0 ms');
axis([-9 9 -0.36 0.72]);
set(gca,'XTick',-9:3:9);
xlabel('interaural level difference (dB)');
ylabel('interaural time difference (ms)');
% Then plot model results
figure; % fig 17 (b)
% Plot line for every condition
for jj = 1:nitds_t
plot(ild_t,output.rcen(jj,:));
hold on;
end
axis([-9 9 -0.36 0.72]);
set(gca,'XTick',-9:3:9);
xlabel('interaural level difference (dB)');
ylabel('interaural time difference (ms)');
%
figure; % fig 17 (c)
for jj = 1:nitds_t
plot(ild_t,output.rmax(jj,:));
hold on;
end
axis([-9 9 -0.36 0.72]);
set(gca,'XTick',-9:3:9);
xlabel('interaural level difference (dB)');
ylabel('interaural time difference (ms)');
end;
end;
%% ------ FIG 18 ----------------------------------------------------------
if flags.do_fig18
% Sampling rate
fs = 44100;
% Frequency of the sinusoid
f = 500;
T = 1/f;
fc = round(freqtoerb(f)); % corresponding frequency channel
% Model parameter
T_int = inf;
w_f = 0.035;
M_f = 6; % not used, if w_f==0
c_s = 0.11;
% NOTE: the longer the signal, the more time we need for computation. On the
% other side N_1 needs to be long enough to eliminate any onset effects.
% Lindemann uses N_1 = 17640. Here I uses only N_1 = 2205 which gives the same
% results for this demo.
N_1 = ceil(25*T*fs);
% Calculate crosscorrelations for 21 ITD points between 0~ms and 1~ms
niacs = 11; % number of used interaural coherence values
ndl = 2*round(fs/2000)+1; % length of the delay line (see bincorr.m)
iac = linspace(0,1,niacs);
output=amt_cache('get','fig18',flags.cachemode);
if isempty(output)
output = zeros(niacs,ndl);
for ii = 1:niacs;
% Generate ITD shifted sinusoid
sig = sig_bincorrnoise(fs,iac(ii),'pink');
% Aplly onset window
sig = rampsignal(sig,[round(N_1/2)-1 0],'tria');
% Calculate cross-correlation (and squeeze due to T_int==inf)
tmp = squeeze(lindemann1986(sig,fs,c_s,w_f,M_f,T_int,N_1));
% Store the needed frequency channel. NOTE: the cross-correlation
% calculation starts with channel 5, so we have to subtract 4.
output(ii,:) = tmp(:,fc-4);
end
amt_cache('set','fig18',output);
end;
if flags.do_plot
% ------ Plotting ------
% Generate time axis
tau = linspace(-1,1,ndl);
% Plot figure for every c_s condition
figure;
mesh(tau,iac,output);
view(0,57);
xlabel('correlation-time tau (ms)');
ylabel('degree of interaural coherence');
end;
end;
%% ------ Subfunctions ----------------------------------------------------
function idx = findnearest(list,val)
[dif,idx] = min(abs(list-val));
function idx = findpeaks(cc)
% FINDPEAKS Finds two peaks of the same height in a given cross-correlation
h = zeros(size(cc,1),1);
for ii = 1:size(cc,1)
% Find a peak in the two halfes of the cross-correlation and subtract them
h(ii) = max(cc(ii,1:30)) - max(cc(ii,31:end));
end
% Find the index with the smallest distance between the two peaks
[val,idx] = min(abs(h));
function m = findmax(list)
[m,idx] = max(list);
d = linspace(-1,1,length(list));
m = d(idx);
Build with Bootstrap