function data = exp_laback2021(varargin)
%EXP_LABACK2021 Experiments from Laback 2021
%   Usage: exp_laback2021;
%
%   EXP_LABACK2021 reproduce Figures from the Laback 2021 paper.
%
%   exp_laback2021 predicts the effect of different types of precursurs on 
%   ILD-based perceived lateralization. It is based on AN model by Smalt, Heinz 
%   and Strickland (2013), which is a binaural version of the Auditory Nerve 
%   Model by Zilany and Bruce (JASA 2006, 2007) that involves ipsi- and contralateral 
%   cochlear compression control via efferent feedback.   
%
%   Examples:
%   ---------
%
%   To display Fig.3 of Laback(2021) use :
%
%     exp_laback2021('fig3');
%
%   AUTHOR: Bernhard Laback, Oct. 2020                                                         
%
%   Url: http://amtoolbox.sourceforge.net/amt-0.10.0/doc/experiments/exp_laback2021.php

% Copyright (C) 2009-2020 Piotr Majdak and the AMT team.
% This file is part of Auditory Modeling Toolbox (AMT) version 1.0.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/>.

bp = which('amt_start');
basepath = bp(1:end-12);
smaltpath=fullfile(basepath,'thirdparty', 'smalt');
addpath(smaltpath);
%addpath('./general_testing');
%Appreviations
%   T Target stimulus
%   P Precursor stimulus

%General parameters
mexzbcatmodel;
Experiment = 1;     %Index of Experiment [1, 2, or 3], as in Laback, B. (2020). Peripheral and Central Factors in Contextual Lateralization based on Interaural Level Difference, Proceedings of Forum Acusticum, 2020, Lyon.
                    %currently, only experiment 1 is fully implemented
modelreps = 4;      %number of repetitions of each random noise stimlus
frozennoise = 0;    %if 1: use frozen noise; if 0: newly generated noises for different P condition
F0cue = 0;          %if 0: T and P have same envelope rate, if 1: T has double envelope rate as P
PrecILD = 10;       %ILD (in dB) of P (has to be positive; a condition with negative ILD is automatically added) 
PrecITD = 0;        %ITD (in �s) of P (has to be positive; a condition with negative ILD is automatically added) [default: 2000]
step=2;             %Stepsize of ILDs of T [default: 2]
NrILD=6;            %Number of ILDs of T, starting from zero) [default: 6]
Modeldelay = 240;   %Delay between stimulus and model response (in samples) [default: 240]
DurT  = 301e-3;     %T duration (in sec) [default: 301e-3]
DurN  = 600e-3;     %P duration (in sec) [default: 600e-3]
DurG  = 500e-3;     %Duration of gap between P and T (in sec) [default: 10e-3]
RampDur = 50;       %Duration of raised-cosine ramps for stimuli (in ms) [default: 50]
stimdb = 60;        %Stimulus intensity (in dB SPL): [default: 60] (MidLev= 60, LowLev=45, HighLev=75)
factor = 1;         %amplification factor for transposed stimuli => not needed any more

if PrecILD < 0
    error('Precursor ILD has to be >= 0');
end
if PrecITD < 0
    error('Precursor ITD has to be >= 0!');
end

%AN model parameters
Fs = 100e3;         % sampling rate in Hz (must be 100, 200 or 500 kHz)
nrep = 20;          % number of stimulus repetitions (e.g., 50);
shocks = [0 0];     % set to 1 to enable shocking MOCR ([ipsi contra])
psthbinwidth = 0.5e-3; % binwidth in seconds; %default: 0.5e-3
gapvec=zeros(1,Fs*DurG);
CF    = 4e3;        %Charecteristic frequency (in Hz) of AN neurons to be predicted [default: 4e3]  
F = CF;             %Center frequency of P and T stimuli (in Hz)
%spontvec = [0.1 5 100];     % spontaneous rate (in spikes/s) of the fiber BEFORE refractory effects stimulus parameters (default: 30)
%spontvec = [0.3 10 130];     % spontaneous rate (in spikes/s) of the fiber BEFORE refractory effects stimulus parameters (default: [0.8 10 130])
spontvec = [0.8 10 130];     % spontaneous rate (in spikes/s) of the fiber BEFORE refractory effects stimulus parameters (default: [0.8 10 130])
cohc  = 1.0;        %OHC gain factor; if 1: normal OHC function; if 0: completely lost OHC function 
cihc  = 1.0;        %IHC gain factor; if 1: normal IHC function; if 0: completely lost IHC function 
BinMOC = 1;         %if 1: binaural MOC, with the amount determined by mocr_max_def; if 0, contralateral MOC, but no ipsilateral MOC
mocr_max_def = 1.0; %Factor for MOC-induced OHC gain reduction; if 1: full reduction; if 0: no gain reduction (MOC activation)

%Define time vectors for P and T
nvec = 0: 1/Fs : DurN-1/Fs; 
tvec = 0: 1/Fs : DurT-1/Fs;

%Generate noises for P and T
PrecNoise=randn(length(nvec),1);
TargNoise=randn(length(tvec),1);

if Experiment == 1
    %P filter
    fl = CF/sqrt(2); 
    fu = fl*(2^1.2); 
    [B_P,A_P]=butter(4,[ fl/(Fs/2) fu/(Fs/2) ]);

    %T filter
    fu = CF*sqrt(2);
    [B_T,A_T]=butter(4,[ fl/(Fs/2) fu/(Fs/2) ]);

elseif Experiment >= 2
    %P filter
    %F0: 125 Hz
    %%%%%%%%%%%
    fu = 138; %CF: 125 Hz; BW: 25 Hz;  
    fl = 113; 
    [B_P,A_P]=butter(3,[ fl/(Fs/2) fu/(Fs/2) ]);
    [D,C]=butter(12,2000/(Fs/2) , 'low'); %low-pass filter for filtering envelope signal 
        
    %T filter
    if F0cue == 0
        B_T = B_P;
        A_T = A_P;
    elseif F0cue == 1
        fu = 263; %CF: 250 Hz; BW: 25 Hz; 
        fl = 238; 
        [B_T,A_T]=butter(3,[ fl/(Fs/2) fu/(Fs/2) ]);
    end  
end


if BinMOC % binaural efferent model fit 2:1 dB ipsi/contra ratio
    mocr_threshold=   -130.5076;  % 2:1
    mocr_slope =    .0163840; % 2:1
    mocr_binauralratio = [  0.4654    0.2332] ;% 2:1 => differs slightly from Tab. 1
else %contralateral MOC only
    mocr_threshold = -130.5076;  % 2:1
    mocr_slope =    .0163840; % 2:1
    mocr_binauralratio = [  0 0.2332] ;% 2:1 => differs slightly from Tab. 1
end

%experimental loop=========================================================
NrFiberTypes = length(spontvec);
for iiiii=1:NrFiberTypes %three different fiber types (based on SR)
    
    spont = spontvec(iiiii);  % spontaneous rate (in spikes/s) of the fiber BEFORE refractory effects stimulus parameters (default: 30)
for iiii=1:modelreps %loop for repetition
for iii=1:4 %loop for four different stimulus conditions (1: no T, 2:diotic P, 3: ipsi P, 4: contra P)
       
%Define precursor conditions
if iii==1 %no precursor
    precmode = 0; %precursor mode: 0 = without precursor, 1 = with precursor 
    ILDNoise = 0; %ILD of precursor noise, positive value favoring left ear
    ILDNoiseLin = 10^(ILDNoise/20);
elseif iii==2 %diotic precursor
    precmode = 1;  
    ILDNoise = 0; 
    ILDNoiseLin = 10^(ILDNoise/20);
    PrecITD_side = 0;
    PrecITD_t=0;
elseif iii==3 %precursor favoring LEFT ear (ipsilateral, because target is always left leading)
    precmode = 1; 
    ILDNoise = PrecILD;
    ILDNoiseLin = 10^(ILDNoise/20);
    PrecITD_side = -1;
    PrecITD_t=PrecITD;
elseif iii==4 %precursor favoring right ear (contralater, because target is always left leading)
    precmode = 1; 
    ILDNoise = PrecILD*-1; 
    ILDNoiseLin = 10^(ILDNoise/20);
    PrecITD_side = 1;
    PrecITD_t=PrecITD;
end
%clear stim diffpsth RespL RespR diffvec timeout meout mocr c1filterout c2filterout c1vihc c2vihc vihc synout psth500k psthbins psthtime pr psth wstart wend;% pr of spike in each bin

if precmode == 0
               
    if frozennoise == 0 %random noise sample for each trial
        TargNoise=randn(length(tvec),1);
    end
    if Experiment == 1
        stimTarg = filter(B_T, A_T, TargNoise);
    elseif Experiment >= 2                        
        LFband=filter(B_T, A_T, TargNoise);
        index=find(LFband >= 0);
        LFband_env=zeros(1, length(tvec));
        LFband_env(index)=LFband(index);
        LFband_filt=filter(D,C,LFband_env);
        TargNoiseTrans = LFband_filt .* sin(2*pi*tvec*CF);
        stimTarg = TargNoiseTrans'*factor;
    end
    stimTarg = local_raisedcosine(stimTarg,Fs,RampDur);
    stimTarg= stimTarg * sqrt(2)* 20e-6*10^(stimdb/20) / rms(stimTarg); 
    
    TotDur=DurT; %Duration of stimulus
elseif  precmode == 1 
    if frozennoise == 0 %random noise sample for each trial
        %Generate noise vectors
        TargNoise=randn(length(tvec),1);
        PrecNoise=randn(length(nvec),1);
    end
    
    if Experiment == 1
        stimTarg = filter(B_T, A_T, TargNoise);
        stimPrec = filter(B_P, A_P, PrecNoise);
    elseif Experiment >= 2   
        %Target
        LFband=filter(B_T, A_T, TargNoise);
        index=find(LFband >= 0);
        LFband_env=zeros(1, length(tvec));
        LFband_env(index)=LFband(index);
        LFband_filt=filter(D,C,LFband_env);
        TargNoiseTrans = LFband_filt .* sin(2*pi*tvec*CF);
        stimTarg = TargNoiseTrans'*factor;
        
        %Precursor
        LFband=filter(B_P, A_P, PrecNoise);
        index=find(LFband >= 0);
        LFband_env=zeros(1, length(nvec));
        LFband_env(index)=LFband(index);
        LFband_filt=filter(D,C,LFband_env);
        PrecNoiseTrans = LFband_filt .* sin(2*pi*nvec*CF);
        stimPrec = PrecNoiseTrans'*factor;
    end
    stimTarg = local_raisedcosine(stimTarg,Fs,RampDur);
    stimTarg= stimTarg * sqrt(2)* 20e-6*10^(stimdb/20) / rms(stimTarg);
    targetlevel=rms(stimTarg);    
    
    stimPrec = local_raisedcosine(stimPrec,Fs,RampDur);
    preclevel=rms(stimPrec);
    
    %Add ILD to precursor
    stimPrec = stimPrec * targetlevel/preclevel;
    
    noiseILDvec_(1,:) = stimPrec * sqrt(ILDNoiseLin);
    noiseILDvec_(2,:) = stimPrec / sqrt(ILDNoiseLin);
    
    %Apply ITD to precursor
    noiseILDvec = local_ITD(noiseILDvec_,PrecITD_side,PrecITD_t, Fs);
    noiseILDvec = noiseILDvec_;
    TotDur=DurT+DurN+DurG; %Duration of stimulus
end
for MOCstat=1:2 %Moc on vs. off modes
    if MOCstat==1
        mocr_max=0;   % minimum gain reduction (relative to COHC value)
    elseif MOCstat==2
        mocr_max=mocr_max_def;   % maximum gain reduction (relative to COHC value)
    end
    
    ILDstart=0;
    %step=5; %4
    %NrILD=3; %4
      for ii=1:NrILD
            %ILD = ILDstart-1+(ii*step); %in dB
            ILD = (ii*step)-step; %in dB %for calculating all ILDs around zero: add -12
            ILDLin = 10^(ILD/20);
            ILDvec(ii)=ILD;            
        if  precmode == 0    
            %Apply ILD
            stim(1,:) = stimTarg*sqrt(ILDLin); %positive ILD value enhances left ear
            stim(2,:) = stimTarg/sqrt(ILDLin);
        elseif  precmode == 1 
            %Apply ILD
            targL = stimTarg*sqrt(ILDLin); %positive ILD value enhances left ear
            targR = stimTarg/sqrt(ILDLin);
            %Total sigal
            %stim(1, length(stimnoise)+1 : length(stimnoise)+length(stimtarg)) = targL;
            %stim(2, length(stimnoise)+1 : length(stimnoise)+length(stimtarg)) = targR;
            stim(1,:) = [noiseILDvec(1,:) gapvec targL']; %stimnoise
            stim(2,:) = [noiseILDvec(2,:) gapvec targR'];
            %figure(1); plot(stim(1,:));
            dfjh=1;
        end

        %pin(2,:) = pin_contra; % comment this line out to run monaural version of model
        [timeout,meout,mocr,c1filterout,c2filterout,c1vihc,c2vihc,vihc,synout,psth500k] ...
            =  zbcatmodel_MOCR(stim,CF,nrep,1/Fs,TotDur*1.5,cohc,cihc,spont,mocr_max,mocr_threshold,mocr_slope,mocr_binauralratio,shocks);

        numears = size(stim,1);   % 1=monaural, 2=binarual
        psthbins = round(psthbinwidth*Fs);  % number of psth500k bins per psth bin
        psthtime = timeout(1:psthbins:end); % time vector for psth
        pr = squeeze(sum(reshape(psth500k,numears,psthbins,length(psth500k)/psthbins),2))/nrep; % pr of spike in each bin
        psth = pr/psthbinwidth; % psth in units of spikes/s        
        
        if 0 %use psth for determining internal ILD
            if precmode == 0   
                wstart=Modeldelay;
            elseif precmode == 1
                wstart=round(length(psth)*0.2/0.9);
            end
            wend=length(psth);
            diffpsth(ii,:,MOCstat) = psth(1,wstart:wend) - psth(2,wstart:wend);
        else %use synapse for determining internal ILD
            if precmode == 0   
                wstart=Modeldelay;   
                wend = Modeldelay + length(stimTarg);
            elseif precmode == 1
                %wstart=round(length(synout)*0.205/0.9); %optimum value: round(length(synout)*0.205/0.9)
                %Optimal
                wstart = Modeldelay + length(noiseILDvec) + length(gapvec); %240 corresponding to processing delay => %adding 15000 samples reverts contra/ipsi effect (as does subtractin) => how possible.
                wend = Modeldelay + length(noiseILDvec) + length(gapvec) + length(stimTarg); 
                %Short
                %wstart = Modeldelay + length(noiseILDvec) + length(gapvec) +  length(stimTarg)*0.3; %240 corresponding to processing delay => %adding 15000 samples reverts contra/ipsi effect (as does subtractin) => how possible.
                %wend = Modeldelay + length(noiseILDvec) + length(gapvec) + length(stimTarg); 
                %Long
                %wstart = Modeldelay + length(noiseILDvec) + length(gapvec) -  length(stimTarg)*0.3; %240 corresponding to processing delay => %adding 15000 samples reverts contra/ipsi effect (as does subtractin) => how possible.
                %wend = Modeldelay + length(noiseILDvec) + length(gapvec) + length(stimTarg); 
                %figure(1); plot(synout(1,:))
                sdjh=1;
            end
            %wend = Modeldelay + length(gapvec) + length(noiseILDvec) + length(stimtarg); 
            %figure(1); plot(synout(1,1:wend));
            MedWin=400; %window length for running average (smooting of response)
            diffpsth(ii,:,MOCstat) = medfilt1(synout(1,wstart:wend),MedWin) - medfilt1(synout(2,wstart:wend), MedWin);
        end

        if 0 %use psth
            disp('Left Ear Spike Rate:'); mean(psth(1,wstart:wend))
            disp('Right Ear Spike Rate:'); mean(psth(2,wstart:wend))
            disp('Spike Rate Difference:'); 
            diff=mean(psth(1,wstart:wend)) - mean(psth(2,wstart:wend))
            RespL(MOCstat, ii) = mean(psth(1,wstart:wend));
            RespR(MOCstat, ii) = mean(psth(2,wstart:wend));
            diffvec(MOCstat, ii) = diff;
        else %use synapse output
            disp('ILD condition:'); 
            ii
            RespL(MOCstat, ii) = mean(medfilt1(synout(1,wstart:wend),MedWin));
            RespR(MOCstat, ii) = mean(medfilt1(synout(2,wstart:wend),MedWin));
            diffvec(MOCstat, ii) = RespL(MOCstat, ii) - RespR(MOCstat, ii);           
        end
    end
end


    if iii==1
        ResultMatrix(1:2,:) = diffvec(1:2,:);
        RespLMatrix(1:2,:) = RespL(1:2,:);
        RespRMatrix(1:2,:) = RespR(1:2,:);
    elseif iii==2
        ResultMatrix(3:4,:) = diffvec(1:2,:);
        RespLMatrix(3:4,:) = RespL(1:2,:);
        RespRMatrix(3:4,:) = RespR(1:2,:);
    elseif iii==3
        ResultMatrix(5:6,:) = diffvec(1:2,:);
        RespLMatrix(5:6,:) = RespL(1:2,:);
        RespRMatrix(5:6,:) = RespR(1:2,:);
    elseif iii==4
        ResultMatrix(7:8,:) = diffvec(1:2,:);
        RespLMatrix(7:8,:) = RespL(1:2,:);
        RespRMatrix(7:8,:) = RespR(1:2,:);
    end   
    
end% %different P conditions
    if iiiii == 1
        ResultMatrixvecLOW(iiii,:,:)=ResultMatrix;
        RespLMatrixLOW(iiii,:,:)=RespLMatrix;
        RespRMatrixLOW(iiii,:,:)=RespRMatrix;
    elseif iiiii == 2
        ResultMatrixvecMID(iiii,:,:)=ResultMatrix;
        RespLMatrixMID(iiii,:,:)=RespLMatrix;
        RespRMatrixMID(iiii,:,:)=RespRMatrix;
    elseif iiiii == 3
        ResultMatrixvecHIGH(iiii,:,:)=ResultMatrix;
        RespLMatrixHIGH(iiii,:,:)=RespLMatrix;
        RespRMatrixHIGH(iiii,:,:)=RespRMatrix;
    end
    sdhg=1;
end %repetition
end %different fiber types

%Calculate Mean and SD across model repetitions
for ILDvec=1:ii
    for condvec=1:length(ResultMatrix(:,1))
        Results_Ave_LOW(ILDvec,condvec) = mean(ResultMatrixvecLOW(:,condvec,ILDvec));
        Results_STD_LOW(ILDvec,condvec) = std(ResultMatrixvecLOW(:,condvec,ILDvec));
        Results_Ave_MID(ILDvec,condvec) = mean(ResultMatrixvecMID(:,condvec,ILDvec));
        Results_STD_MID(ILDvec,condvec) = std(ResultMatrixvecMID(:,condvec,ILDvec));
        Results_Ave_HIGH(ILDvec,condvec) = mean(ResultMatrixvecHIGH(:,condvec,ILDvec));
        Results_STD_HIGH(ILDvec,condvec) = std(ResultMatrixvecHIGH(:,condvec,ILDvec));
        
        Results_Ave_LOW_L(ILDvec,condvec) = mean(RespLMatrixLOW(:,condvec,ILDvec));
        Results_STD_LOW_L(ILDvec,condvec) = std(RespLMatrixLOW(:,condvec,ILDvec));
        Results_Ave_MID_L(ILDvec,condvec) = mean(RespLMatrixMID(:,condvec,ILDvec));
        Results_STD_MID_L(ILDvec,condvec) = std(RespLMatrixMID(:,condvec,ILDvec));
        Results_Ave_HIGH_L(ILDvec,condvec) = mean(RespLMatrixHIGH(:,condvec,ILDvec));
        Results_STD_HIGH_L(ILDvec,condvec) = std(RespLMatrixHIGH(:,condvec,ILDvec));
 
        Results_Ave_LOW_R(ILDvec,condvec) = mean(RespRMatrixLOW(:,condvec,ILDvec));
        Results_STD_LOW_R(ILDvec,condvec) = std(RespRMatrixLOW(:,condvec,ILDvec));
        Results_Ave_MID_R(ILDvec,condvec) = mean(RespRMatrixMID(:,condvec,ILDvec));
        Results_STD_MID_R(ILDvec,condvec) = std(RespRMatrixMID(:,condvec,ILDvec));
        Results_Ave_HIGH_R(ILDvec,condvec) = mean(RespRMatrixHIGH(:,condvec,ILDvec));
        Results_STD_HIGH_R(ILDvec,condvec) = std(RespRMatrixHIGH(:,condvec,ILDvec));  
    end
end

ResultMatrix= ResultMatrix';

%Plotting (as in Laback, 2020).
if Experiment == 1
    %Calculate Weighted Mean
    W=[70 20 10]; %Weighting Coefficients for Low, Mid, and High-SR fibers
    Results_Ave_WMean = ( Results_Ave_LOW*W(1) + Results_Ave_MID*W(2) + Results_Ave_HIGH*W(3) ) / sum(W); 
    Results_STD_WMean = ( Results_STD_LOW*W(1) + Results_STD_MID*W(2) + Results_STD_HIGH*W(3) ) / sum(W); 
   
    ILDvec=0:2:10;
    figure(1); %Low-SR fibers
    subplot(2, 2, 1)
    hold on;
    errorbar(ILDvec,Results_Ave_LOW(:,2)', Results_STD_LOW(:,2)', 'b'); %NoPrec
    errorbar(ILDvec,Results_Ave_LOW(:,4)', Results_STD_LOW(:,4)', 'r'); %Diotic
    errorbar(ILDvec,Results_Ave_LOW(:,6)', Results_STD_LOW(:,6)', 'g'); %Ipsi
    errorbar(ILDvec,Results_Ave_LOW(:,8)', Results_STD_LOW(:,8)', 'c'); %Contra    
    axis([0 10 -10 70]); xlabel('Target ILD (dB)'); ylabel('Response Azimuth (�)'); text(6,60,'LOW-SR fibers')
    hold off;

    %figure(2); %Mid-SR fibers
        subplot(2, 2, 2)
    hold on;
    errorbar(ILDvec,Results_Ave_MID(:,2)', Results_STD_MID(:,2)', 'b'); %NoPrec
    errorbar(ILDvec,Results_Ave_MID(:,4)', Results_STD_MID(:,4)', 'r'); %Diotic
    errorbar(ILDvec,Results_Ave_MID(:,6)', Results_STD_MID(:,6)', 'g'); %Ipsi
    errorbar(ILDvec,Results_Ave_MID(:,8)', Results_STD_MID(:,8)', 'c'); %Contra    
    axis([0 10 -10 160]); xlabel('Target ILD (dB)'); ylabel('Response Azimuth (�)'); text(6,0,'MID-SR fibers')
    hold off;
    
    %figure(3); %High-SR fibers
        subplot(2, 2, 3)
    hold on;
    errorbar(ILDvec,Results_Ave_HIGH(:,2)', Results_STD_HIGH(:,2)', 'b'); %NoPrec
    errorbar(ILDvec,Results_Ave_HIGH(:,4)', Results_STD_HIGH(:,4)', 'r'); %Diotic
    errorbar(ILDvec,Results_Ave_HIGH(:,6)', Results_STD_HIGH(:,6)', 'g'); %Ipsi
    errorbar(ILDvec,Results_Ave_HIGH(:,8)', Results_STD_HIGH(:,8)', 'c'); %Contra    
    axis([0 10 0 160]); xlabel('Target ILD (dB)'); ylabel('Response Azimuth (�)'); text(6,140,'HIGH-SR fibers')
    hold off;

    %figure(4); %Weighted Mean
        subplot(2, 2, 4)
    hold on;
    errorbar(ILDvec,Results_Ave_WMean(:,2)', Results_STD_WMean(:,2)', 'b'); %NoPrec
    errorbar(ILDvec,Results_Ave_WMean(:,4)', Results_STD_WMean(:,4)', 'r'); %Diotic
    errorbar(ILDvec,Results_Ave_WMean(:,6)', Results_STD_WMean(:,6)', 'g'); %Ipsi
    errorbar(ILDvec,Results_Ave_WMean(:,8)', Results_STD_WMean(:,8)', 'c'); %Contra    
    axis([0 10 0 70]); xlabel('Target ILD (dB)'); ylabel('Response Azimuth (�)'); text(1,60,'Weighted Mean')
    hold off;

elseif Experiment == 2

end

if 0 %plotting
    figure(60); plot(ILDvec, diffvec(1,:), 'g--', ILDvec, diffvec(2,:),'r-'); axis([0, 15, 0 100]); %axis([0, 15, 0 35]);
     xlabel('External ILD (dB)')
     ylabel('Interaural spike rate difference (spikes/s)'); legend('Without MOCR','With MOCR')
     %save 'C:\Users\Bernhard Laback\Documents\Events\2015\20151026 Closing Auditory Loop Hannover Laback\Simulations\ILD\Matlabvectors\6k 100Hz_SAM SR=5 1_2 ratio sym' diffvec diffpsth

    figure(61); plot(ILDvec, RespL(1,:), 'g--', ILDvec, RespL(2,:),'r--' , ILDvec, RespR(1,:), 'g-',ILDvec, RespR(2,:),'r-'); axis([0, 15, 0 300]); %axis([0, 15, 80 130]);
    xlabel('External ILD (dB)')
    ylabel('Interaural spike rate difference (spikes/s)'); legend('Left, Without MOCR','Left, With MOCR','Right, Without MOCR','Right, With MOCR' )

    if precmode == 0
        sdp = size(diffpsth);
        for j= 1:2
            figure(7+j); surf(1:sdp(2),1:sdp(1),diffpsth(:,:,j),'EdgeColor', 'None');
            xlabel('PSTH bins'); ylabel ('External ILD (dB)'); zlabel('Interaural spike rate difference (spikes/s)'); 
            colorbar; %view(2);  
            AZ = 0; EL = 90;  view(AZ,EL);
            %axis([0, 15, 0, 6, -100, 300]);
        end
        figure(10); surf(1:sdp(2),1:sdp(1),diffpsth(:,:,2)-diffpsth(:,:,1),'EdgeColor', 'None'); colorbar; AZ = 0; EL = 90;  view(AZ,EL); 
        xlabel('PSTH bins'); ylabel ('External ILD (dB)'); zlabel('Interaural spike rate difference (spikes/s)');
    end
    figure(11); plot(1:length(diffpsth(ii,:,1)), diffpsth(ii,:,1), 1:length(diffpsth(ii,:,1)), diffpsth(ii,:,2)); 
    xlabel('PSTH bins'); 
    ylabel ('Interaural spike rate difference (spikes/s)'); 
    %figure(12); plot(1:length(diffsynout(ii,:,1)), diffsynout(ii,:,1), 1:length(diffsynout(ii,:,1)), diffsynout(ii,:,2)); xlabel('PSTH bins'); ylabel ('Interaural spike rate difference (spikes/s)');
end
end

function [out] = local_ITD(in,side,ITD, fs)
%Function to apply ITD to input signal 
%B. Laback, 03.06.2016  
%Call: [out] = ITD(in,side,ITD, fs)
%Paramters
%in         binaural signal (2-D vector) to which ITD shall be applied
%side       -1: leading left; 0:zero ITD (zeros added at end of signal at both ears; +1: leading right 
%ITD        ITD in us
%fs         sampling rate (in Hz)
%out        output signal containing ITD


inputL = in(1,:);
inputR = in(2,:);

ITDsamples = round(ITD/(1000000/fs));
Lag = zeros(1,ITDsamples);

%Introduce ITD
if side < 0
    out(1,:) = [inputL Lag];
    out(2,:) = [Lag inputR];
elseif side == 0
    out(1,:)=[inputL Lag];
    out(2,:)=[inputR Lag];
elseif side > 0
    out(1,:) = [Lag inputL];
    out(2,:) = [inputR Lag];
end
end

function signal = local_raisedcosine(signal,fs,timeS)
%signal = local_raisedcosine(signal,fs,timeS)
%applies raised-cosine Onset and Offset ramp
%Inputs:
%signal...input signal mono
%fs...sampling rate [Hz]
%timeS is ramp time [s]
%AUTHOR: Bernhard Laback, 2021

size = length(signal);

time = timeS/1000;
rs = round(time*fs); 
for i=1:(rs-1)
   ampfac = 0.5 - 0.5*cos(pi*i/rs);
   signal(i) = ampfac*signal(i);
end

for i=(size-rs+1):size 
   ampfac = 0.5 - 0.5*cos(pi*(size-i)/rs);
	signal(i) = ampfac*signal(i);
end
end 
