THE AUDITORY MODELING TOOLBOX
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