function data = exp_osses2022(varargin)
%EXP_OSSES2022 - 
%
%   Usage: data = exp_osses2022(flags)
%          data = exp_osses2022(..., 'idxs', models)
%
%   exp_osses2021 reproduces Figs. X and XX from Osses et al. (2022), where 
%   the following models are compared: 
%     'dau1997'           : Model #1 from Dau et al., (1997)
%     'zilany2014'        : Model #2 from Zilany et al., (2014)
%     'verhulst2015'      : Model #3 from Verhulst et al., (2015)
%     'verhulst2018'      : Model #4 from Verhulst et al., (2018)
%     'bruce2018'         : Model #5 from Bruce et al., (2018)
%     'king2019'          : Model #6 from King et al., (2019)
%     'relanoiborra2019'  : Model #7 from Relano-Iborra et al., (2019)
%     'osses2021'         : Model #8 from Osses and Kohlrausch (2021)
%
%   The following flags can be specified;
%
%     'redo'    Recomputes data for specified figure
%
%     'plot'    Plot the output of the experiment. This is the default.
%
%     'no_plot' Don't plot, only return data.
%  
%     'models'  Vector selecting the models, if not all models to process. 
%               For example 'models',[1 3 6] shows the data for Dau et al.,
%               (1997), Verhulst et al., (2015), and King et al., (2019) only.
%
%
%   To display Fig. 3 of Osses et al., (2022) use :
%
%     out = exp_osses2022('fig3');
%
%   To display Fig. 4 of Osses et al., (2022) use :
%
%     out = exp_osses2022('fig4');
%
%   To display Figs. 6 and 7 of Osses et al., (2022) use :
%
%     out = exp_osses2022('fig6');
%
%   To display Fig. 8 of Osses et al., (2022) use :
%
%     out = exp_osses2022('fig8');
%
%   To display Fig. 9 of Osses et al., (2022) use :
%
%     out = exp_osses2022('fig9');
%
%   To display Fig. 11 of Osses et al., (2022) use :
%
%     out = exp_osses2022('fig11');
%
%   To display Fig. 12 of Osses et al., (2022) use :
%
%     out = exp_osses2022('fig12');
%
%   References:
%     I. C. Bruce, Y. Erfani, and M. S. R. Zilany. A phenomenological model
%     of the synapse between the inner hair cell and auditory nerve:
%     Implications of limited neurotransmitter release sites. Hearing
%     Research, 360:40--54, 2018.
%     
%     M. S. A. Zilany, I. C. Bruce, and L. H. Carney. Updated parameters and
%     expanded simulation options for a model of the auditory periphery. The
%     Journal of the Acoustical Society of America, 135(1):283--286, Jan.
%     2014.
%     
%     S. Verhulst, H. Bharadwaj, G. Mehraei, C. Shera, and
%     B. Shinn-Cunningham. Functional modeling of the human auditory
%     brainstem response to broadband stimulation. jasa, 138(3):1637--1659,
%     2015.
%     
%               verhulst2018 king2019 relanoiborra2019 osses2021
%
%   Url: http://amtoolbox.sourceforge.net/amt-0.10.0/doc/experiments/exp_osses2022.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/>.

%   #Author: Alejandro Osses (2020-2021): Integration in the AMT
%   #Author: Piotr Majdak (2021): Adaptations for the AMT 1.0

data = [];

definput.import={'amt_cache'};
definput.flags.type={'missingflag','fig3','fig4','fig5','fig6','fig7', ...
                       'fig8','fig9','fig10','fig11','fig12','fig13ab','fig13cd','fig14'};

definput.flags.plot={'plot','no_plot'};
definput.keyvals.models=[];

[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;

%%%
FontSize = 13; % FontSize for all figure axes
fs = 100e3; % 100 kHz of sampling rate, this is arbitrary
dBFS = 94; % i.e., amplitude 1 = 1 Pa = 94 dB SPL re 2x10^{-5} Pa

%%% for Figs 3, 4, 5, 6 (confirm for the first figures)
models = {'dau1997','zilany2014','verhulst2015','verhulst2018','bruce2018','king2019','relanoiborra2019','osses2021'};
Colours = {'b',il_rgb('Green'),'r',il_rgb('LightSkyBlue'),il_rgb('Maroon'),'m','k',il_rgb('mediumorchid')}; % ,'k'};
Markers = {'o','s','<','>','v','^','p','d'};
MarkersSize = [10 10 10 10 10 10 10 10];
LineStyle = {'-','-','-.','-','--','-','-','-'};
LineWidth = [2 2 2 2 2 2 2 2];

if ~isempty(keyvals.models),
	warning('Not all models selected'); 
	idxs=keyvals.models;
	Markers = Markers(idxs); 
	MarkersSize = MarkersSize(idxs); 
	Colours = Colours(idxs); 
	LineStyle = LineStyle(idxs); 
	LineWidth = LineWidth(idxs); 
	models = models(idxs);
end

N_models = length(models);

figure_handle = []; % multiple figures will be generated
figure_name   = [];
   
%% ------ FIG 3 Osses, Verhulst, and Majdak 2021 -------------------------
if flags.do_fig3
    
    % Params for FFT:
    K = fs/2; % for FFT
    % Memory allocation
    h_dB    = nan(K,N_models);
    hmax_dB = nan(1,N_models);
    h0_dB   = nan(K,N_models);
    k2remove = [];
    label2use = [];
    
    for k = 1:N_models
        switch models{k}
            case {'dau1997','king2019'}
                % No middle ear filter
                me_type = [];
            case {'relanoiborra2019'}
                me_type = 'lopezpoveda2001';
                label2use{end+1} = 'relanoiborra2019, osses2021';
            case 'osses2021' 
                % uses the same as relanoiborra (no need to calculate it again)
                me_type = [];
            case 'zilany2014'
                me_type = 'zilany2009';
                label2use{end+1} = 'zilany2014, bruce2018';
            case 'bruce2018'
                me_type = [];
            case 'verhulst2015'
                me_type = 'verhulst2015';
                label2use{end+1} = 'verhulst2015';
            case 'verhulst2018'
                me_type = 'verhulst2018';
                label2use{end+1} = 'verhulst2018';
        end
    
        if ~isempty(me_type)

            [b,a] = middleearfilter(fs,me_type);
        
            Nr_cascaded = size(b,1);
            if Nr_cascaded == 1 % no cascaded
                [h,w]=freqz(b,a,K);
            else
                h = ones(K,1);
                for j = 1:Nr_cascaded
                    [htmp,w]=freqz(b(j,:),a(j,:),K);
                    h = h.*htmp;
                end
            end
            f = w/pi*(fs/2); % freq. of the FFT bins in Hz
            h_dB(:,k)   = 20*log10(abs(h));
            hmax_dB(k)  = max(h_dB(:,k));
            h0_dB(:,k)  = h_dB(:,k) - hmax_dB(k);
            
        else
            k2remove = [k2remove k];
        end
    end
    
    % Preparing the plots with the four middle ear filters:
    models(k2remove) = []; 
    Colours(k2remove) = []; 
    Markers(k2remove) = []; 
    MarkersSize(k2remove) = [];
    LineStyle(k2remove) = []; 
    LineWidth(k2remove) = [];
    
    h_dB(:,k2remove) = [];
    hmax_dB(k2remove) = [];
    h0_dB(:,k2remove) = [];
    
    %%%
    b_oear = headphonefilter(fs); % Pralong filter
    htmp=freqz(b_oear,1,K);
    hOuter_ear = 20*log10(abs(htmp));
    
    %%%
    pl = [];
    figure;
    if strcmp(models{4},'relanoiborra2019')
        hCombined = h0_dB(:,4) + hOuter_ear;
        pl(end+1) = semilogx(f,hCombined,'--','Color',[0.5 0.5 0.5],'LineWidth',4); hold on, grid on;
    end
    
    for k = 1:length(models)
        opts_plot = {'Color',Colours{k},'LineWidth',LineWidth(k),'LineStyle',LineStyle{k}};
        pl(end+1) = semilogx(f,h0_dB(:,k),opts_plot{:}); hold on, grid on;
        
        L1 = [min(f) max(f); -3 -3];
        L2 = [f'; h0_dB(:,k)'];
        P = il_InterX(L1,L2);
        
        if ~isempty(P)
            f_min03_low(k)  = P(1,1);
            f_min03_high(k) = P(1,end);
        else
            f_min03_low(k)  = nan;
            f_min03_high(k) = nan;
        end
            
        
        L1 = [min(f) max(f); -15 -15];
        P = il_InterX(L1,L2);
        
        if ~isempty(P)
            f_min15_low(k)  = P(1,1);
            f_min15_high(k) = P(1,end);
        else
            f_min15_low(k)  = nan;
            f_min15_high(k) = nan;
        end
        
    end
    
    label2use{end+1} = '(outer+middle ear)';
    xlim([80 30000]);
    % legend(label2use,'Location','NorthEastOutside');
    hl = legend([pl(2:end) pl(1)],label2use);
    set(hl,'FontSize',8);
    
    set(gca,'FontSize',FontSize);
     
    f2tick = [63 125 250 500 1000 2000 4000 8000 16000];
    set(gca,'XTick',f2tick);
    set(gca,'YTick',-21:3:3);
    ylim([-23 6])
     
    xlabel('Frequency (Hz)');
    ylabel('Amplitude (dB re. max)');
    % ht = title('Middle-ear frequency response');
    % set(ht,'FontSize',FontSize);
     
    Pos = get(gcf,'Position');
    Pos(3:4) = [640 360]; % fix width
    set(gcf,'Position',Pos);
     
    data.figure_flag   = 'do_fig3';
    data.figure_handle = gcf;
    data.figure_name   = 'fig03-me-resp'; % middle ear response
     
    data.description = 'Frequency response of the middle ears used in the selected models';
    data.hmax_dB = hmax_dB;
    data.f_min03_low  = f_min03_low;
    data.f_min03_high = f_min03_high;
    data.f_min15_low  = f_min15_low;
    data.f_min15_high = f_min15_high;
    data.evaluated_models = models;
end
   
%% ------ FIG 4 -------------------------
if flags.do_fig4

    k2remove = [];
    N_models_here = N_models; % Bruce2018, is removed later...
    
    %%% 1. Generating the input signals: %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    % Signal parameters:
    dur=100e-3; % 50 ms
    lvls=0:10:100; % Level of the test signals
    N_lvls = length(lvls);
    
    f0s = [1000 500 4000];
    
    dt = 1/fs; % \Delta t in s
    t=0:dt:dur-dt;

    cfs = il_m2hz(401); % default characteristic frequencies from Verhulst2012
           
    Pos34 = [500 350];
    
    for i = 1:length(f0s)
        f0_target = f0s(i);
    
        if f0_target == 1000
            idx_lvls = find(lvls==100); % 0-dB reference
        else
            idx_lvls = 1:N_lvls;
        end
        
        idx = find(cfs>f0_target,1,'last');
        f0 = cfs(idx);
    
        f0_lo = audtofreq(freqtoaud(f0_target)-1); % min 1 ERB
        idx_off(1) = find(cfs>f0_lo,1,'last');
        
        f0_hi = audtofreq(freqtoaud(f0_target)+1); % plus 1 ERB
        idx_off(2) = find(cfs>f0_hi,1,'last');
        
        fprintf('actual f0=%.1f Hz using Verhulst''s mapping (bin=%.0f)\n',f0,idx);

        % Basis input signal
        insig = zeros(length(t),N_lvls); % Memory allocation
        insig_orig = sin(2*pi*f0*t(:)); % the same sinusoid will be scaled at different levels

        % Up/down cosine ramp (fixed)
        dur_ramp_ms = 10;
        dur_ramp = round((dur_ramp_ms*1e-3)*fs); % duration ramp in samples
         
        rp    = ones(size(insig_orig)); 
        rp(1:dur_ramp) = rampup(dur_ramp);
        rp(end-dur_ramp+1:end) = rampdown(dur_ramp);
        %%%
 
        % Calibration of the input signals:
        for j = idx_lvls
            insig(:,j) = scaletodbspl(insig_orig,lvls(j),dBFS); % setdbspl
            insig(:,j) = rp.*insig(:,j);
        end
    
        %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
        N = length(lvls(idx_lvls)); 
        vrms        = nan(N,N_models); % for the on-frequency values
        vrms_off_lo = nan(N,N_models); % for the off-frequency values
        vrms_off_hi = nan(N,N_models); % for the off-frequency values
     
        k_test = 1;
        
        insig_here = insig(:,idx_lvls);
         
        for k = 1:N_models_here
            
            %%% Loading flags and keyvals
            [fg,kv] = il_get_flags(models{k});
            fc_flags   = fg.fc_flags;
            afb_flags  = fg.afb_flags;
            afb_kv     = kv.afb_keyvals;
            % ihc_flags  = fg.ihc_flags;
            noihc_flags  = fg.noihc_flags;
            noan_flags = fg.noan_flags; % No auditory nerve module
            nomfb_flags = fg.nomfb_flags;
            
            fname = ['fig04_' models{k} '-f0-' num2str(f0_target) '-Hz'];
            c = amt_cache('get',fname,flags.cachemode); 

            if ~isempty(c)
                bRun = 0;
                out_all    = c.out_all;
                out_lo_all = c.out_lo_all;
                out_hi_all = c.out_hi_all;
            else
                bRun = 1;
                out_all    = [];
                out_lo_all = [];
                out_hi_all = [];
            end
            
            if bRun
                amt_disp(['Calculating ' models{k} '...']);
                switch models{k}
                    case 'dau1997'
                        for j = 1:N
                            fc_kv = {'flow',f0,'fhigh',f0,'basef',f0,'dboffset',dBFS};
                            out = dau1997(insig_here(:,j),fs,fc_kv{:},afb_flags{:},noihc_flags{:},noan_flags{:});
                            out_all(:,end+1) = out;

                            fc_kv = {'flow',f0_lo,'fhigh',f0_lo,'basef',f0_lo,'dboffset',dBFS};
                            out = dau1997(insig_here(:,j),fs,fc_kv{:},afb_flags{:},noihc_flags{:},noan_flags{:});
                            out_lo_all(:,end+1) = out;

                            fc_kv = {'flow',f0_hi,'fhigh',f0_hi,'basef',f0_hi,'dboffset',dBFS};
                            out_hi = dau1997(insig_here(:,j),fs,fc_kv{:},afb_flags{:},noihc_flags{:},noan_flags{:});
                            out_hi_all(:,end+1) = out;
                        end

                    case 'zilany2014'

                        kv={'numH',0,'numM',0,'numL',0,'fiberType',4,'nrep',1};
                        for j = 1:N
                            [~,~,~,c1,c2] = zilany2014(insig_here(:,j),fs,f0,kv{:});
                            idxs = 1:length(insig_here); % to remove the zero padding
                            out = c1(:)+c2(:);
                            out = out(idxs);
                            out_all(:,end+1) = out;

                            [~,~,~,c1,c2] = zilany2014(insig_here(:,j),fs,f0_lo,kv{:});
                            out = c1(:)+c2(:);
                            out = out(idxs);
                            out_lo_all(:,end+1) = out;

                            [~,~,~,c1,c2] = zilany2014(insig_here(:,j),fs,f0_hi,kv{:});
                            out = c1(:)+c2(:);
                            out = out(idxs);
                            out_hi_all(:,end+1) = out;
                        end

                    case 'verhulst2015'
                        out_model = verhulst2015(insig_here,fs,fc_flags{:},afb_flags{:},noihc_flags{:},noan_flags{:});
                        for j = 1:N
                            out = out_model(j).v(:,idx);
                            out_all(:,end+1) = out;

                            out = out_model(j).v(:,idx_off(1));
                            out_lo_all(:,end+1) = out;

                            out = out_model(j).v(:,idx_off(2));
                            out_hi_all(:,end+1) = out;
                        end

                    case 'verhulst2018'
                        out_model = verhulst2018(insig_here,fs,fc_flags{:},afb_flags{:},noihc_flags{:},noan_flags{:});
                        for j = 1:N
                            out = out_model(j).v(:,idx);
                            out_all(:,end+1) = out;

                            out = out_model(j).v(:,idx_off(1));
                            out_lo_all(:,end+1) = out;

                            out = out_model(j).v(:,idx_off(2));
                            out_hi_all(:,end+1) = out;
                        end

                    case 'bruce2018'
                        kv={'numH',0,'numM',0,'numL',0,'nrep',1};
                        for j = 1:N
                            output = bruce2018(insig_here(:,j),fs,f0,kv{:});
                            idxs = 1:length(insig_here); % to remove the zero padding
                            out = output.C1(:)+output.C2(:);
                            out = out(idxs);
                            out_all(:,end+1) = out;

                            output = bruce2018(insig_here(:,j),fs,f0_lo,kv{:});
                            out = output.C1(:)+output.C2(:);
                            out = out(idxs);
                            out_lo_all(:,end+1) = out;

                            output = bruce2018(insig_here(:,j),fs,f0_hi, kv{:});
                            out = output.C1(:)+output.C2(:);
                            out = out(idxs);
                            out_hi_all(:,end+1) = out;
                        end

                    case 'king2019'
                        for j = 1:N
                            fc_kv = {'flow',f0,'fhigh',f0,'basef',f0,'dboffset',dBFS};
                            out = king2019(insig_here(:,j),fs,afb_kv{:},fc_kv{:},afb_flags{:},noihc_flags{:},noan_flags{:});
                            out_all(:,end+1) = out;

                            fc_kv = {'flow',f0_lo,'fhigh',f0_lo,'basef',f0_lo,'dboffset',dBFS};
                            out = king2019(insig_here(:,j),fs,afb_kv{:},fc_kv{:},afb_flags{:},noihc_flags{:},noan_flags{:});
                            out_lo_all(:,end+1) = out;

                            fc_kv = {'flow',f0_hi,'fhigh',f0_hi,'basef',f0_hi,'dboffset',dBFS};
                            out = king2019(insig_here(:,j),fs,afb_kv{:},fc_kv{:},afb_flags{:},noihc_flags{:},noan_flags{:});
                            out_hi_all(:,end+1) = out;
                        end
                        
                    case 'relanoiborra2019'
                        for j = 1:N
                            fc_kv = {'flow',f0,'fhigh',f0,'basef',f0,'erbspacebw','no_ihc','no_an'};
                            [~,~,out] = relanoiborra2019_featureextraction(insig_here(:,j),fs,fc_kv{:},afb_flags{:});
                            out_all(:,end+1) = out;

                            fc_kv = {'flow',f0_lo,'fhigh',f0_lo,'basef',f0_lo,'erbspacebw','no_ihc','no_an'};
                            [~,~,out] = relanoiborra2019_featureextraction(insig_here(:,j),fs,fc_kv{:},afb_flags{:});
                            out_lo_all(:,end+1) = out;

                            fc_kv = {'flow',f0_hi,'fhigh',f0_hi,'basef',f0_hi,'erbspacebw','no_ihc','no_an'};
                            [~,~,out] = relanoiborra2019_featureextraction(insig_here(:,j),fs,fc_kv{:},afb_flags{:});
                            out_hi_all(:,end+1) = out;
                        end
                        
                    case 'osses2021'
                        for j = 1:N
                            fc_kv = {'flow',f0,'fhigh',f0,'basef',f0,'dboffset',dBFS};
                            out = osses2021(insig_here(:,j),fs,fc_kv{:},afb_flags{:},noihc_flags{:},noan_flags{:});
                            out_all(:,end+1) = out;

                            fc_kv = {'flow',f0_lo,'fhigh',f0_lo,'basef',f0_lo,'dboffset',dBFS};
                            out = osses2021(insig_here(:,j),fs,fc_kv{:},afb_flags{:},noihc_flags{:},noan_flags{:});
                            out_lo_all(:,end+1) = out;

                            fc_kv = {'flow',f0_hi,'fhigh',f0_hi,'basef',f0_hi,'dboffset',dBFS};
                            out = osses2021(insig_here(:,j),fs,fc_kv{:},afb_flags{:},noihc_flags{:},noan_flags{:});
                            out_hi_all(:,end+1) = out;
                        end
                end
                
                if ~isempty(out_all)
                    c = [];
                    c.out_all = out_all;
                    c.out_lo_all = out_lo_all;
                    c.out_hi_all = out_hi_all;
                    amt_cache('set',fname,c);
                end
            end
            
            if ~isempty(out_all)
                % Now calculation of RMS values:
                for j = 1:N
                    vrms(j,k)        = il_rmsdb(out_all(:,j));
                    vrms_off_lo(j,k) = il_rmsdb(out_lo_all(:,j));
                    vrms_off_hi(j,k) = il_rmsdb(out_hi_all(:,j));
                end
            end
        end
    
        if i == 1
            % Preparing the plots with the four middle ear filters:
            models(k2remove) = []; 
            Colours(k2remove) = []; 
            Markers(k2remove) = []; 
            MarkersSize(k2remove) = [];
            LineStyle(k2remove) = []; 
            LineWidth(k2remove) = [];
            
            N_models_here = length(models); % number of models after removing
            
            vrms(:,k2remove) = [];
            vrms_off_lo(:,k2remove) = [];
            vrms_off_hi(:,k2remove) = [];
        end
            
        if f0_target == 1000
            v_ref_0_dB = vrms;
            % v_ref_0_dB = [   -7.3561  -44.0567  -77.9298 -101.8625  -44.0567]; % calculated at 1 kHz
            data.v_ref = v_ref_0_dB;
        else
            if f0_target == 500
                prefix4title = {'a) ','b) ','c) '};
            else
                prefix4title = {'d) ','e) ','f) '};
            end
            v_ref = v_ref_0_dB;
            
            Format = [];
            for k = 1:N_models_here
                 Format{k} = {'Color',Colours{k},'LineStyle',LineStyle{k}, ...
                    'LineWidth',LineWidth(k),'Marker',Markers{k},'MarkerFaceColor','w', ...
                    'MarkerSize',MarkersSize(k)};
            end
            %%% End formatting options
            
            figure;
            for k = 1:N_models_here
                pl(k)=plot(lvls,vrms(:,k)-v_ref(k),Format{k}{:}); 
                grid on, hold on;
            end
            
            title(sprintf('%sCF at %.0f Hz, on-frequency bin',prefix4title{1},f0));
            xlabel('Input level (dB SPL)');
            ylabel('Output level (dB re. max)');
            ylim([-103 13])
            xlim([min(lvls)-5 max(lvls)+5])
            set(gca,'XTick', 0:10:100)
            set(gca,'YTick',-100:10:10)
            figure_handle(end+1) = gcf;
            figure_name{end+1}   = sprintf('fig04-IO-at-%.0f-Hz',f0); 
            Pos = get(gcf,'Position');
            Pos(3:4) = Pos34;
            set(gcf,'Position',Pos);
            
            %%%%
            figure;
            for k = 1:N_models_here
                plot(lvls,vrms_off_lo(:,k)-v_ref(k),Format{k}{:}); 
                grid on, hold on
            end
            title(sprintf('%sCF at %.0f Hz, off frequency bin (-1 ERB)',prefix4title{2},f0_lo));
            xlabel('Input level (dB SPL)');
            ylabel('Output level (dB re. max)');
            ylim([-103 13])
            xlim([min(lvls)-5 max(lvls)+5])
            set(gca,'XTick', 0:10:100)
            set(gca,'YTick',-100:10:10)
            figure_handle(end+1) = gcf;
            figure_name{end+1}   = sprintf('fig04-IO-at-%.0f-Hz-off-lo',f0); 
            Pos = get(gcf,'Position');
            Pos(3:4) = Pos34;
            set(gcf,'Position',Pos);
            
            figure;
            for k = 1:N_models_here
                plot(lvls,vrms_off_hi(:,k)-v_ref(k),Format{k}{:}); 
                grid on, hold on
            end
            title(sprintf('%sCF at %.0f Hz, off frequency bin (+1 ERB)',prefix4title{3},f0_hi));
            xlabel('Input level (dB SPL)');
            ylabel('Output level (dB re. max)');
            ylim([-103 13])
            xlim([min(lvls)-5 max(lvls)+5])
            set(gca,'XTick', 0:10:100)
            set(gca,'YTick',-100:10:10)
            Pos = get(gcf,'Position');
            Pos(3:4) = Pos34;
            set(gcf,'Position',Pos);
            
            if i == length(f0s)
                text4leg = models;
                if strcmp(models{2},'zilany2014')
                    text4leg{2} = [text4leg{2} ',bruce2018'];
                end
                hl = legend(text4leg,'Location','SouthEast');
                set(hl,'FontSize',8);
            end
            
            figure_handle(end+1) = gcf;
            figure_name{end+1}   = sprintf('fig04-IO-at-%.0f-Hz-off-hi',f0); 
            
            data(k_test).vrms = vrms;
            data(k_test).vrms_off_lo = vrms_off_lo;
            data(k_test).vrms_off_hi = vrms_off_hi;
            data(k_test).f0 = f0;
            
            k_test = k_test+1;
        end
        
    end
    data.figure_flag   = 'do_fig4';
    data.figure_handle = figure_handle;
    data.figure_name   = figure_name;
    data.models = models;
    
    data.v_ref_0_dB = v_ref_0_dB;
    data.insig = insig;
    data.insig_orig = insig_orig;
    data.fs    = fs;
    data.ramp  = rp;
    data.dBFS  = dBFS;
    data.lvls  = lvls;
end

%% ------ FIG 5 Frequency selectivity ------------------------------------
if flags.do_fig5

    error('Not implemented yet. Sorry!'); 
    
    k2remove = [];
    
    fsig = 4000; % Hz
    basef = fsig;

    % - click response, QERBs and N values
    dt = 1/fs;

    dur_stim_total = 100e-3; % Total duration of the stimulus
    % dur_click      = 80e-6; % 80 us
    dur_click      = 100e-6; % 25 ms, 50 ms

    t_clk=0:dt:dur_stim_total-dt;
    clk = zeros(1,length(t_clk));
    idxi = round(10e-3*fs+1);
    idxf = round(idxi+dur_click*fs);
    clk(idxi:idxf)=1;

    idxi = round((50e-3+10e-3)*fs+1);
    idxf = round(idxi+dur_click*fs);
    clk(idxi:idxf)=-1;
        
    lvls=[30 70]; % 0:10:90;
    N_lvls = length(lvls);
    level_factor = 2*sqrt(2)*(2e-5.*(10.^(lvls/20))); 

    insigs = zeros(length(clk),N_lvls);

    for i=1:N_lvls
        % lvl=corr_spl(i);
        insigs(:,i) = level_factor(i).*clk(:);
    end

    % Amplitude at 30 dB:     2e-5*10^(30/20)*sqrt(2)*2;
    % Amplitude at 70 dB:     2e-5*10^(70/20)*sqrt(2)*2;
        
    N_samples = size(insigs,1);
    [cf,x] = il_m2hz(401); % Get characteristic frequencies from Verhulst models
    cf_max = 10000;
    cf_min =  125;
    %%%
    N_sigs = 31; % Approximate number of bands
    bin_cf_max = find(cf>cf_max,1,'last');
    bin_cf_min = find(cf>cf_min,1,'last');
    df_bin = floor(abs(bin_cf_max-bin_cf_min)/N_sigs);
        
    idx_cf = bin_cf_min:-df_bin:bin_cf_max;
    fc_ref = cf(idx_cf);
    N_sigs = length(fc_ref); % Exact number of bands
    
    for k = 1:N_models
        %%% Loading flags and keyvals
        [fg,kv] = il_get_flags(models{k});
        fc_flags   = fg.fc_flags;
        afb_flags  = fg.afb_flags;
        afb_kv     = kv.afb_keyvals;
        noihc_flags  = fg.noihc_flags;
        noan_flags = fg.noan_flags; % No auditory nerve module
        nomfb_flags = fg.nomfb_flags;

        for j = 1:N_lvls
            fname = ['fig05_' models{k} '-QERB-' num2str(lvls(j))];
            c = amt_cache('get',fname,flags.cachemode); 

            if ~isempty(c)
                bRun = 0;

                outsig = c.outsig;
                fc     = c.fc;
            else
                bRun = 1;
                outsig = [];
            end

            if bRun
                amt_disp(['Calculating ' models{k} '...']);
                switch models{k}
                    case 'dau1997'
                        for ii = 1:N_sigs
                            fc_kv = {'flow',fc_ref(ii),'fhigh',fc_ref(ii),'basef',fc_ref(ii),'dboffset',dBFS};
                            
                            [outsig(:,ii),fc(ii)] = dau1997(insigs(:,j),fs,fc_kv{:},afb_flags{:},noihc_flags{:},noan_flags{:});

                        end
                        
                    case 'zilany2014'
                        outsig = [];
                        kv={'numH',0,'numM',0,'numL',0,'fiberType',4,'nrep',1,'reptime',1.2};
                        for ii = 1:N_sigs 
                            [~,~,~,c1,c2] = zilany2014(insigs(:,j),fs,fc_ref(ii),kv{:});                            
                            outsig(:,ii) = c1(:)+c2(:);
                        end
                        
                    case 'verhulst2015'
                        if j == 1
                            outs = verhulst2015(insigs,fs,fc_flags{:},afb_flags{:},noihc_flags{:},noan_flags{:});
                            CF_here = outs(1).cf;
                        end
                        outsig = [];
                        for ii = 1:N_sigs
                            outsig(:,ii) = outs(j).v(:,idx_cf(ii));
                        end
                        
                    case 'verhulst2018'
                        if j == 1
                            outs = verhulst2018(insigs,fs,fc_flags{:},afb_flags{:},noihc_flags{:},noan_flags{:});
                            CF_here = outs(1).cf;
                        end
                        outsig = [];
                        for ii = 1:N_sigs
                            outsig(:,ii) = outs(j).v(:,idx_cf(ii));
                        end
                        
                    case 'bruce2018'
                        kv={'numH',0,'numM',0,'numL',0,'nrep',1};
                        outsig = [];
                        for ii = 1:N_sigs 
                            output = bruce2018(insigs(:,j)',fs,fc_ref(ii),kv{:});                            
                            outsig(:,ii) = output.C1(:)+output.C2(:);
                        end

                    case 'king2019'
                        for ii = 1:N_sigs
                            fc_kv = {'flow',fc_ref(ii),'fhigh',fc_ref(ii),'basef',fc_ref(ii),'dboffset',dBFS};
                            [outsig(:,ii),fc(ii)] = king2019(insigs(:,j),fs,fc_kv{:},afb_flags{:},noihc_flags{:},noan_flags{:});
                        end
                        
                    case 'relanoiborra2019'
                        for ii = 1:N_sigs
                            fc_kv = {'flow',fc_ref(ii),'fhigh',fc_ref(ii),'basef',fc_ref(ii),'erbspacebw','no_internalnoise','no_ihc','no_an'};
                            [~,~,outsig(:,ii),fc(ii)] = relanoiborra2019_featureextraction(insigs(:,j), fs,fc_kv{:});
                        end
                        
                    case 'osses2021'
                        for ii = 1:N_sigs
                            fc_kv = {'flow',fc_ref(ii),'fhigh',fc_ref(ii),'basef',fc_ref(ii),'dboffset',dBFS};
                            [outsig(:,ii),fc(ii)] = osses2021(insigs(:,j),fs,fc_kv{:},afb_flags{:},noihc_flags{:},noan_flags{:});                            
                        end   
                end            
                if ~isempty(outsig)
                    c = [];
                    c.outsig = outsig;
                    c.fc     = fc;
                    amt_cache('set',fname,c);
                end
            end
            
            if k == 1
                idx_CF = find(fc_ref<=basef,1,'last');
            end
            K = N_samples/2;
            fz = (1:K)/K*(fs/2);

            if ~isempty(outsig)
                switch j
                    case 1 % Low intensity
                        spec_CF_lo(:,k)         = 0.5*(abs(freqz(outsig(1:K,idx_CF  ),1,K))+abs(freqz(outsig(K+1:end,idx_CF  ),1,K))); % only first half of the samples
                        spec_CF_min_1ERB_lo(:,k)= 0.5*(abs(freqz(outsig(1:K,idx_CF-1),1,K))+abs(freqz(outsig(K+1:end,idx_CF-1),1,K)));
                        for ii = 1:length(fc_ref)
                            N_here = length(outsig(:,ii));
                            insig_here = reshape(outsig(:,ii),N_here/2,2);
                            for kk = 1:2                                         
                                [BW_low_each(ii,kk,k), QERB_low_each(ii,kk,k),extra] = il_Get_ERB_estimation(insig_here(:,kk), fc_ref(ii), fs);
                               
                                Q03_low_each(ii,kk,k) = extra.Q03;
                                Q10_low_each(ii,kk,k) = extra.Q10; 
                            end
                            % [BW_low(ii,k), QERB_low(ii,k)] = il_Get_ERB_estimation(outsig(:,ii), fc_ref(ii), fs);

                            disp('')
                        end
                    case 2 % Higher intensity
                        spec_CF_hi(:,k)         = 0.5*(abs(freqz(outsig(1:K,idx_CF),1,K))+abs(freqz(outsig(K+1:end,idx_CF),1,K)));
                        spec_CF_min_1ERB_hi(:,k)= 0.5*(abs(freqz(outsig(1:K,idx_CF-1),1,K))+abs(freqz(outsig(K+1:end,idx_CF-1),1,K)));
                        for ii = 1:length(fc_ref)
                            N_here = length(outsig(:,ii));
                            insig_here = reshape(outsig(:,ii),N_here/2,2);
                            for kk = 1:2
                                [BW_high_each(ii,kk,k), QERB_high_each(ii,kk,k),extra] = il_Get_ERB_estimation(insig_here(:,kk), fc_ref(ii), fs);

                                Q03_high_each(ii,kk,k) = extra.Q03;
                                Q10_high_each(ii,kk,k) = extra.Q10; 
                            end
                        end % end ii
                end % end switch j
            end % end ~isempty
        end % end for j     
    end % end for models

    %%%
    % Preparing the plots
    models(k2remove) = []; 
    Colours(k2remove) = []; 
    Markers(k2remove) = []; 
    MarkersSize(k2remove) = [];
    LineStyle(k2remove) = []; 
    LineWidth(k2remove) = [];
     
    N_models_here = length(models); % number of models after removing
    
    BW_low_each(:,:,k2remove) = [];
    QERB_low_each(:,:,k2remove) = [];
    Q03_low_each(:,:,k2remove) = [];
    Q10_low_each(:,:,k2remove) = [];
    BW_high_each(:,:,k2remove) = [];
    QERB_high_each(:,:,k2remove) = [];
    Q03_high_each(:,:,k2remove) = [];
    Q10_high_each(:,:,k2remove) = [];
    spec_CF_lo(:,k2remove) = [];
    spec_CF_min_1ERB_lo(:,k2remove) = [];
    spec_CF_hi(:,k2remove) = [];
    spec_CF_min_1ERB_hi(:,k2remove) = [];
    
    idx = find(fc_ref<= 8000);
    fc_ref = fc_ref(idx);
    
    QERB_low_each = QERB_low_each(idx,:,:);
    Q03_low_each = Q03_low_each(idx,:,:);
    Q10_low_each = Q10_low_each(idx,:,:);
    
    QERB_high_each = QERB_high_each(idx,:,:);
    Q03_high_each = Q03_high_each(idx,:,:);
    Q10_high_each = Q10_high_each(idx,:,:);
    
    %%% -------------------------------------------------------------------
    figure; 
    
    % Analytical expressions:
    fc_ref_here = [125 fc_ref 10000];
    alpha = 0.3;  % Shera2002, Table I, col 'human'
    beta  = 12.7; % Shera2002, Table I, col 'human'
    Qerb_Shera_extended = beta*((fc_ref_here/1000).^alpha);
    Qerb_Shera = Qerb_Shera_extended(2:end-1);
    % Q10 = Qerb*0.505+0.2085; % Ibrahim2010, Eq. 40.2
    pl = [];
    pl(end+1) = semilogx(fc_ref_here,Qerb_Shera_extended,'-','Color',il_rgb('LightGray'),'LineWidth',5); hold on
    
    BW_Glasberg = 24.7*(4.37*(fc_ref_here/1000)+1);
    Qerb_Glasberg_extended = fc_ref_here./BW_Glasberg;
    Qerb_Glasberg = Qerb_Glasberg_extended(2:end-1);
    % Q10   = Qerb*0.505+0.2085; % Ibrahim2010, Eq. 40.2
    pl(end+1) = plot(fc_ref_here,Qerb_Glasberg_extended,'-','Color',il_rgb('Gray'),'LineWidth',5);
    
    for k = 1:N_models_here % little adjustments in format for the coming figures
        switch models{k}
            case 'dau1997'
                Marker = [Markers{k} '-'];
                MSize = 7;
                LW = 1;
            case 'zilany2014'
                Marker = [Markers{k} LineStyle{k}];
                MSize = 7;
                LW = LineWidth(k);
            otherwise
                Marker = LineStyle{k};
                MSize = MarkersSize(k);
                LW = LineWidth(k);
        end
        format_kv = {Marker,'Color',Colours{k},'MarkerSize',MSize,'LineWidth',LW};
        Q_low_here(:,k)     = mean(QERB_low_each(:,:,k),2); % average of the positive and negative click
        pl(end+1) = semilogx(fc_ref,Q_low_here(:,k),format_kv{:},'MarkerFaceColor',Colours{k}); hold on
        % plot(fc_ref,QERB_higher(:,k),'s--',format_kv{:},'MarkerFaceColor','w');
    end
    grid on;
    XT = [20 125 250 500 1000 2000 4000 8000 16000];
    set(gca,'XTick',XT);
    YT = 2:2:26;
    set(gca,'YTick',YT);

    ylim([1 27])
    
    leg = models;
    k = 2;
    if strcmp(leg{k},'zilany2014')
        leg{k} = 'zilany2014,bruce2018';
    end
    leg{end+1} = 'Shera (2001)';
    leg{end+1} = 'G & M (1990)';
    
    legend([pl(3:end) pl(1:2)],leg,'Location','NorthWest');
    
    xlabel('Frequency (Hz)');
    ylabel('Q factor');

    Pos = get(gcf,'Position');
    Pos(4) = 350;
    set(gcf,'Position',Pos);
    
    figure_handle(end+1) = gcf;
    figure_name{end+1}   = sprintf('fig05-QERB-low-level');
    
    % ---------------------------------------------------------------------
    for k = 1:N_models_here
        
        X_ref_empirical = mean(QERB_low_each(:,:,k),2);
        Y_10 = Q10_low_each(:,1,k);
        Y_03 = Q03_low_each(:,1,k);
        
        % Fit a polynomial p of degree 1 to the (x,y) data:
        p10(k,:) = polyfit(X_ref_empirical,Y_10,1);
        p03(k,:) = polyfit(X_ref_empirical,Y_03,1);
        
        switch models{k}
            case {'dau1997','relanoiborra2019','king2019','osses2021'}
                X_ref(:,k) = Qerb_Glasberg(:);
            case {'zilany2014','bruce2018','verhulst2015','verhulst2018'}
                X_ref(:,k) = Qerb_Shera(:);
        end
        p10_with_formula(k,:) = polyfit(X_ref(:,k),Y_10,1); 
        p03_with_formula(k,:) = polyfit(X_ref(:,k),Y_03,1); % figure; plot(X_ref,Y_10); hold on; plot(X_ref,Y_03,'r')
        
        Q03_recons(:,k) = X_ref(:,k)*p03_with_formula(k,1) + p03_with_formula(k,2);
        
        % plot(X_ref(:,k),Q03_recons(:,k),format_kv{:},'MarkerFaceColor',Colours{k}); hold on
    end
    
    for k = 1:N_models_here
        f_current = 125;
        f_high_lim = 8000;

        count(k) = 0;
        while f_current < f_high_lim
            count(k) = count(k)+1; 
            f_low(count(k),k) = f_current;
            if count(k)>10000, break; end % break the loop if got stuck
            % Q_here = interp1(fc_ref,X_ref(:,k),f_low(count(k),k),'linear','extrap');
            Q_here = interp1(fc_ref, Q03_recons(:,k), f_low(count(k),k),'linear','extrap');
            
            BW_here = f_low(count(k),k)/Q_here;
            f_high(count(k),k) = f_low(count(k),k) + BW_here;
            f_current = f_high(count(k),k);
        end
        %%% Fittings that can be useful to report at some point:
        % idx = [2 3 4];   p10_bio = prctile(p10(idx,:),50); p03_bio = prctile(p03(idx,:),50);
        % idx = [1 5 6 7]; p10_per = prctile(p10(idx,:),50); p03_per = prctile(p03(idx,:),50);
        % 
        % idx = [2 3 4];   p10_bio_me = mean(p10(idx,:)); p03_bio_me = mean(p03(idx,:));
        % idx = [1 5 6 7]; p10_per_me = mean(p10(idx,:)); p03_per_me = mean(p03(idx,:));

        % figure; 
        % plot(X_ref,p(2)*X_ref+p(1),'k-'); hold on;
        % plot(X_ref,Y_10,'or-');
        % disp('')

        % % Evaluate the fitted polynomial p and plot:
        % f = polyval(p,X_ref);
        % figure;
        % plot(X_ref,Y_10,'o',X_ref,f,'-')
        % legend('data','linear fit')
        %%% End fittings
        
        fc_max = max(f_high);
        fc_min = f_low(1,:); 
    end
    % Values reported in Table II:
    counts   = count;
    erb_step = (freqtoaud(fc_max)-freqtoaud(fc_min))./(count-1); 
    
    
    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    figure; % Now high levels only
    
    % Analytical expressions:
    pl = [];
    pl(end+1) = semilogx(fc_ref_here,Qerb_Shera_extended   ,'-','Color',il_rgb('LightGray'),'LineWidth',5); hold on
    pl(end+1) = plot(    fc_ref_here,Qerb_Glasberg_extended,'-','Color',il_rgb('Gray')     ,'LineWidth',5);
    
    for k = 1:N_models_here
        switch models{k}
            case 'dau1997'
                Marker = [Markers{k} '-'];
                MSize = 7;
                LW = 1;
            case 'zilany2014'
                Marker = [Markers{k} LineStyle{k}];
                MSize = 7;
                LW = LineWidth(k);
            otherwise
                Marker = LineStyle{k};
                MSize = MarkersSize(k);
                LW = LineWidth(k);
        end
        format_kv = {Marker,'Color',Colours{k},'MarkerSize',MSize, ...
            'LineWidth',LW};
        Q_high_here(:,k)     = mean(QERB_high_each(:,:,k),2); % average of the positive and negative click
        pl(end+1) = semilogx(fc_ref,Q_high_here(:,k),format_kv{:},'MarkerFaceColor',Colours{k}); hold on
        % plot(fc_ref,QERB_higher(:,k),'s--',format_kv{:},'MarkerFaceColor','w');
    end
    grid on;
    set(gca,'XTick',XT);
    YT = 2:2:26;
    set(gca,'YTick',YT);

    ylim([1 27])
    
    xlabel('Frequency (Hz)');
    ylabel('Q factor');
    
    Pos = get(gcf,'Position');
    Pos(4) = 350;
    set(gcf,'Position',Pos);
    
    figure_handle(end+1) = gcf;
    figure_name{end+1}   = sprintf('fig05-QERB-high-level');
    
    % ---------------------------------------------------------------------
    figure; % Difference Low - high
    
    % Analytical expressions:
    pl = [];
    for k = 1:N_models_here
        switch models{k}
            case 'dau1997'
                Marker = [Markers{k} '-'];
                MSize = 7;
                LW = 1;
            case 'zilany2014'
                Marker = [Markers{k} LineStyle{k}];
                MSize = 7;
                LW = LineWidth(k);
            otherwise
                Marker = LineStyle{k};
                MSize = MarkersSize(k);
                LW = LineWidth(k);
        end
        format_kv = {Marker,'Color',Colours{k},'MarkerSize',MSize,'LineWidth',LW};
        pl(end+1) = semilogx(fc_ref,Q_low_here(:,k)-Q_high_here(:,k),format_kv{:},'MarkerFaceColor',Colours{k}); hold on
        % plot(fc_ref,QERB_higher(:,k),'s--',format_kv{:},'MarkerFaceColor','w');
    end
    grid on;
    
    Pos = get(gcf,'Position');
    Pos(4) = 250;
    set(gcf,'Position',Pos);
    
    set(gca,'XTick',XT);
    ylim([-4.5 12.5])
    YT = -4:2:12;
    set(gca,'YTick',YT);
    
    xlabel('Frequency (Hz)');
    ylabel('Q factor difference');
     
    figure_handle(end+1) = gcf;
    figure_name{end+1}   = sprintf('fig05-QERB-difference-high-low');
    
    X_ref_empirical = [];
    for k = 1:N_models_here
        X_ref_empirical(:,k) = Q_high_here(:,k);
        Y_10 = Q10_high_each(:,1,k);
        Y_03 = Q03_high_each(:,1,k); 
    end
    
    Q03_recons_hi = squeeze(Q03_high_each(:,1,:));
    
    f_low = [];
    f_high = [];
    count = [];
    for k = 1:N_models_here
        f_current = 125;
        f_high_lim = 8000;

        count(k) = 0;
        while f_current < f_high_lim
            count(k) = count(k)+1;
            f_low(count(k),k) = f_current;
            if count(k)>10000, break; end % break the loop if got stuck
            % Q_here = interp1(fc_ref,X_ref(:,k),f_low(count(k),k),'linear','extrap');
            Q_here = interp1(fc_ref, Q03_recons_hi(:,k), f_low(count(k),k),'linear','extrap');
            if isnan(Q_here)
                Q_here = Q03_recons_hi(1,k);
            end           
            
            BW_here = f_low(count(k),k)/Q_here;
            f_high(count(k),k) = f_low(count(k),k) + BW_here;
            
            f_current = f_high(count(k),k);
            disp('')
            % BW = 
            % f_high = 
        end
        
        fc_max = max(f_high);
        fc_min = f_low(1,:); 
    end
    counts(2,:)   = count;
    erb_step(2,:) = (freqtoaud(fc_max)-freqtoaud(fc_min))./(count-1); 
    
    data.figure_handle = figure_handle;
    data.figure_name   = figure_name;
    data.figure_flag   = 'do_fig5';
    data.models        = models;
    data.counts        = counts;
    data.counts_description = 'Number of bands required to cover a frequency range from 125 to 8000 Hz';
    data.erb_step      = erb_step;
    data.erb_step_description = 'ERB spacing that would ensure filters overlapped at their -3 dB points';
end

%% ------ FIG 6 or 7 -----------------------------------------------------
if flags.do_fig6 || flags.do_fig7
    
    k2remove = [];
    N_models_here = N_models; % Bruce2018, is removed later...
    
    t_ms_all = [];
    vihc_all = [];
    ti_all = [];
    
    %%% 1. Generating the input signals: %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%   
    % Signal parameters:
    dur = 0.08; % 80 ms
    lvl = 80; % 80 dB SPL 
    dur_i = 0.02; % starting time for AC assessment
    dur_f = 0.07; % ending time for AC assessment
    
    % Arbitrary zero padding before and after the signal:
    dur_sil = 50e-3;
    sil_bef = zeros(round(dur_sil*fs),1);

    idxi = length(sil_bef)+1; % index corresponding to the first signal sample after the zero padding
    
    % AC will be calculated between dur_i and dur_f, excluding the silence:
    idxi_ac = idxi+round(dur_i*fs)-1;
    idxf_ac = idxi+round(dur_f*fs)-1;
 
    [cf,x] = il_m2hz(401); % Get characteristic frequencies from Verhulst models
    
    cf_max = 4000;
    cf_min =  125;
    %%%
    N_sigs = 12;
    
    vac  = nan(N_models, N_sigs);
    vdc  = nan(size(vac));
    
    %%% Getting f0 ('f0_method'):
    df_bin = 25;
    idx_cf(N_sigs) = find(cf>cf_max,1,'last');

    for i = N_sigs-1:-1:1
        idx_cf(i) = idx_cf(i+1)+df_bin;
    end
    f0 = cf(idx_cf);
    %%%
    
    dt = 1/fs; % \Delta t in s
    t=0:dt:dur-dt;
    
    for i = 1:N_sigs
        insig(:,i) = sin(2*pi*f0(i)*t); % the same sinusoid will be scaled at different levels
    end
    % Up/down cosine ramp (fixed)
    dur_ramp_ms = 5; % in ms
    dur_ramp = round((dur_ramp_ms*1e-3)*fs); % duration ramp in samples

    rp    = ones(size(insig(:,1))); 
    rp(1:dur_ramp) = rampup(dur_ramp);
    rp(end-dur_ramp+1:end) = rampdown(dur_ramp);
    
    % Adjusting the signal level before the ramp and silences are added:
    insig = scaletodbspl(insig,lvl,dBFS);
    % Adding the silence and ramps to the signals:
    insig = [repmat(sil_bef,1,N_sigs); repmat(rp,1,N_sigs).*insig];
    %%%
     
    for k = 1:N_models_here

        bInclude = 1; % models with bInclude = 0 no IHC AD/DC is assessed later...
        
        %%% Loading flags and keyvals
        [fg,kv] = il_get_flags(models{k});
        fc_flags   = fg.fc_flags;
        afb_flags  = fg.afb_flags;
        afb_kv     = kv.afb_keyvals;
        ihc_flags  = fg.ihc_flags;
        noan_flags = fg.noan_flags; % No auditory nerve module
        nomfb_flags = fg.nomfb_flags;
        %%%

        fname = ['fig06_and_07_ihc-' models{k}];
        c = amt_cache('get',fname,flags.cachemode); 
        
        if ~isempty(c)
            bRun = 0;
            
            vihc = c.vihc;
            fs   = c.fs;
            fc   = c.fc;
        else
            bRun = 1;
            vihc = [];
        end

        amt_disp(['Calculating ' models{k} '...']);
        switch models{k}
            case 'dau1997'
                zoom_higher(k) = 1;
                if bRun
                    for j = 1:N_sigs
                        fc_here = f0(j);
                        fc_kv = {'flow',fc_here,'fhigh',fc_here,'basef',fc_here,'dboffset',dBFS};
                        [outsig,fc] = dau1997(insig(:,j),fs,fc_kv{:},afb_flags{:},ihc_flags{:},noan_flags{:});
                        vihc(:,j) = outsig;
                    end
                    
                    c = [];
                    c.vihc = vihc;
                    c.fs   = fs;
                    c.fc   = fc;
                    amt_cache('set',fname,c);
                end
            
            case 'zilany2014'
                zoom_higher(k) = 3;
                if bRun
                    kv={'numH',0,'numM',0,'numL',0,'fiberType',4,'nrep',1};
                    for j = 1:N_sigs
                        [~,~,vihc(:,j)] = zilany2014(insig(:,j),fs,f0(j),kv{:});
                    end
                    c = [];
                    c.vihc = vihc;
                    c.fs   = fs;
                    c.fc   = fc;
                    amt_cache('set',fname,c);
                end
            
            case 'verhulst2015'
                zoom_higher(k) = 3;
                if bRun
                    outsig = verhulst2015(insig,fs,fc_flags{:},afb_flags{:},ihc_flags{:},noan_flags{:});
                    for j = 1:N_sigs
                        vihc(:,j) = outsig(j).ihc(:,idx_cf(j));
                    end
                    c = [];
                    c.vihc = vihc;
                    c.fs   = fs;
                    c.fc   = fc;
                    amt_cache('set',fname,c);
                end
                
            case 'verhulst2018'
                zoom_higher(k) = 3;
                if bRun
                    outsig = verhulst2018(insig,fs,fc_flags{:},afb_flags{:},ihc_flags{:},noan_flags{:});

                    for j = 1:N_sigs
                        vihc(:,j) = outsig(j).ihc(:,idx_cf(j));
                    end
                    c = [];
                    c.vihc = vihc;
                    c.fs   = fs;
                    c.fc   = fc;
                    amt_cache('set',fname,c);
                end
            
            case 'bruce2018'
                bInclude = 0;
                k2remove = [k2remove k];
                  % Code to test bruce2018, use only with selected models
%                 zoom_higher(k) = 3;
%                 if bRun
%                     kv={'numH',0,'numM',0,'numL',0,'nrep',1};
%                     for j = 1:N_sigs
%                         [~,output] = bruce2018(insig(:,j),fs,f0(j),kv{:});
%                         vihc(:,j)=output.vihc; 
%                     end
%                     c = [];
%                     c.vihc = vihc;
%                     c.fs   = fs;
%                     c.fc   = fc;
%                     amt_cache('set',fname,c);
%                 end

                
            case 'king2019'
                zoom_higher(k) = 1;
                if bRun
                    for j = 1:N_sigs
                        fc_here = f0(j);
                        fc_kv = {'flow',fc_here,'fhigh',fc_here,'basef',fc_here,'dboffset',dBFS};
                        [outsig,fc] = king2019(insig(:,j),fs,fc_kv{:},afb_kv{:},afb_flags{:},ihc_flags{:},noan_flags{:});
                        vihc(:,j) = outsig;
                    end
                    
                    c = [];
                    c.vihc = vihc;
                    c.fs   = fs;
                    c.fc   = fc;
                    amt_cache('set',fname,c);
                end
                
            case 'relanoiborra2019'  
                zoom_higher(k) = 3;
                if bRun
                    for j = 1:N_sigs
                        fc_here = f0(j);
                        fc_kv = {'flow',fc_here,'fhigh',fc_here,'basef',fc_here,'erbspacebw','no_an'};
                        [~,~,outsig,fc] = relanoiborra2019_featureextraction(insig(:,j), fs,fc_kv{:});
                        vihc(:,j) = outsig;
                    end
                    c = [];
                    c.vihc = vihc;
                    c.fs   = fs;
                    c.fc   = fc;
                    amt_cache('set',fname,c);
                end
                
            case 'osses2021'
                zoom_higher(k) = 3;
                if bRun
                    for j = 1:N_sigs
                        fc_here = f0(j);
                        fc_kv = {'flow',fc_here,'fhigh',fc_here,'basef',fc_here,'dboffset',dBFS};
                        [outsig,fc] = osses2021(insig(:,j),fs,fc_kv{:},afb_flags{:},ihc_flags{:},noan_flags{:});
                        vihc(:,j) = outsig;
                    end
                    c = [];
                    c.vihc = vihc;
                    c.fs   = fs;
                    c.fc   = fc;
                    amt_cache('set',fname,c);
                end
        end

        if bInclude
            dc_reg = 1:length(sil_bef);
            ac_reg = idxi_ac:idxf_ac;
                
            [ac_target,dc_target,vr] = il_get_ac_dc_osses2022(vihc,ac_reg,dc_reg);

            % vrest(k) = median(vihc(1:length(sil_bef),1)); % median of the silent initial segment            
            vac(k,:) = ac_target;
            vdc(k,:) = dc_target; % mean(vihc(idxi_ac:idxf_ac,:))-vrest(k);
            vrest(k,1) = vr(1);
            
            % Plotting: -------------------------------------------------------
            % % Normalisation to 1.3 the maximum value:
            offy_step = 1.5*max(max(abs(vihc)));

            t_ms = 1000*(1:size(vihc,1))/fs;
            t2plot = [35 105];
            idxi_plot = find(t_ms<=t2plot(1),1,'last');
            idxf_plot = find(t_ms<=t2plot(2),1,'last');
            idxs = idxi_plot:idxf_plot-1;

            v_dec_norm = (vihc-vrest(k))/offy_step; % normalisation to maximum value
            
            idxs_N = 9:12;
            v_dec_norm(:,idxs_N) = zoom_higher(k)*v_dec_norm(:,idxs_N);
            
            idxs2store = idxs(1:5000);
            t_ms_all = [t_ms_all t_ms(idxs2store)];
            ti_all(end+1) = size(vihc_all,1)+1;
            vihc_all = [vihc_all; v_dec_norm(idxs2store,:)];
            
            vmin(k,:) = min(vihc(idxs2store,:));
            
        end
    end     
     
    models(k2remove) = []; 
    Colours(k2remove) = []; 
    Markers(k2remove) = []; 
    MarkersSize(k2remove) = [];
    LineStyle(k2remove) = []; 
    LineWidth(k2remove) = [];
    
    zoom_higher(k2remove)  = [];
    vac(k2remove,:)  = [];
    vdc(k2remove,:)  = [];
    vmin(k2remove,:) = [];
    
    N_models_here = length(models); % number of models after removing
             
    %%% All formatting options in one variable. The formatting options
    %     for the bruce2018 are overwritten for visibility reasons, 
    %     given that it provides the same outputs of zilany2014
    Format = [];
    for k = 1:N_models_here
        Format{k} = {'Color',Colours{k},'LineStyle',LineStyle{k}, ...
        'LineWidth',LineWidth(k),'Marker',Markers{k},'MarkerFaceColor','w', ...
        'MarkerSize',MarkersSize(k)};
    end
    %%% End formatting options
          
    ti_all(end+1) = length(t_ms_all);
    toffset = t_ms(idxs(1));
    tf = 0;
    idxi = 1;
    YL = [-.5 N_sigs];
    plt = [];
    
    idx_split = 5;
    
    panel_labels = {'a) ','b) ','c) ','d) ','e) ','f) ','g) '};
    for k = 1:N_models_here
        if strcmp(models{k},'zilany2014')
            leg{k} = [panel_labels{k} 'zilany2014,bruce2018'];
        else
            leg{k} = [panel_labels{k} models{k}];
        end
    end
        
    for k = 1:N_models_here
        
        if k<idx_split
            figure(1);
            if k == 1
                figure_handle(end+1) = gcf;
                figure_name{end+1} = 'fig06-ihc-waveforms';
            end
        else
            figure(2);
            if k == idx_split
                figure_handle(end+1) = gcf;
                figure_name{end+1} = 'fig06-ihc-waveforms-2';
            end
        end
        
        idxi = ti_all(k);
        idxf = ti_all(k+1)-1;
        t_here    = t_ms_all(idxi:idxf);
        vihc_here = vihc_all(idxi:idxf,:);
        
        for j = 1:N_sigs
            offy = j-1; % (j-1) lower to higher freqs from bottom to top. Use (N_sigs-j) otherwise
            if j == 1
                plt(end+1) = plot(t_here-toffset+tf,vihc_here(:,j)+offy,'Color',Colours{k}); 
            else
                plot(t_here-toffset+tf,vihc_here(:,j)+offy,'Color',Colours{k}); 
            end
            hold on
            if j == 1 && (k == 1  || k == idx_split)
                ylim(YL);
            end
            plot(tf*[1 1],YL,'k-');
            
            if (k == 1 || k == idx_split)
                text(2,offy+.3,sprintf('%.0f Hz',f0(j)),'Color','k','FontSize',8);
                
                grid on
                
                dur_plot = 1000*length(idxs2store)/fs;
                XL = [0 (idx_split-1)*dur_plot];
                xlim(XL);
                Pos = get(gcf,'Position');
                Pos(3:4) = [650 600];
                set(gcf,'Position',[0 0 650 650]);
                
                XT  = 0:10:XL(end);
                XTL = [0 repmat([10:10:dur_plot],1,(idx_split-1))];
                
                XT(dur_plot/10+1:dur_plot/10:end) = [];
                XTL(dur_plot/10+1:dur_plot/10:end) = [];
                set(gca,'XTick',XT);
                set(gca,'XTickLabel',XTL);
                YT = 0:N_sigs;
                set(gca,'YTick',YT);
                set(gca,'YTickLabel','');
                
            end
            
            if j >= 9
                if k < idx_split
                    text((k*dur_plot-8),offy-.05,sprintf('x %.0f',zoom_higher(k)),'Color','k','FontSize',8);
                else
                    text(((k-idx_split+1)*dur_plot-8),offy-.05,sprintf('x %.0f',zoom_higher(k)),'Color','k','FontSize',8);
                end
            end
        end % end for j 
        
        if mod(k,4) == 0
            x_coor = (4-1)*.25+.02;
        else
            x_coor = ((mod(k,4))-1)*.25+.02;
        end
        y_coor = .98;
        text(x_coor,y_coor,leg{k},'Units','Normalize','FontSize',8);
        
        if k == idx_split-1
            tf = 0;
        else
            tf = tf+(t_here(end)+dt*1000-toffset);
        end
        ylabel('Amplitude (a.u.)');
        xlabel('Time (ms)');
    end
    plot(tf*[1 1],YL,'k-'); % last time for the second figure
       
    ra  = vac ./vdc;
    
    hl = [];
    
    % Fig. 7 --------------------------------------------------------------
    figure;
    for k = 1:N_models_here
        pl(k) = loglog(f0,ra(k,:),Format{k}{:}); grid on, hold on
        xlabel('Frequency (Hz)');
        ylabel('IHC AC/DC ratio');         
    end
    set(gca,'XTick',[80 125 250 500 1000 2000 4000 8000]);
    YT = get(gca,'YTick'); 
    set(gca,'YTickLabel',YT)
    xlim([100 6000])
    ylim([0.008 180])
    
    leg = models;
    if strcmp(leg{2},'zilany2014')
        leg{2} = 'zilany2014,bruce2018';
    end
    legend(hl,leg,'Location','SouthWest')
     
    figure_handle(end+1) = gcf;
    figure_name{end+1} = 'fig07-IHC-AC-DC';
    
    data.figure_flag   = 'do_fig6 or do_fig7';
    data.figure_handle = figure_handle;
    data.figure_name   = figure_name;
    
    data.insig = insig;
    data.fs    = fs;
    data.ramp  = rp;
    
    data.models = models;
    data.ACDCratio = ra;
    data.AC    = vac;
    data.DC    = vdc;
    
    data.dBFS  = dBFS;
    data.lvl   = lvl;
    data.f0_approx = round(f0');
    data.f0_exact  = round(cf(idx_cf)');
    
    data.vrest = vrest;
    data.vmin  = vmin;
end

%% ------ FIG. 8 ---------------------------------------------------------
if flags.do_fig8
    
    outsig_all = [];
    t_an_all = [];
    psth_binwidth = 0.5e-3; % width of the PSTH bin, in seconds
    nrep=100; % number of repetitions in PSTH calculations
    
    % 1. Stimulus creation:
    dur = 300e-3; % 300 ms
    lvl = 70;
    dBFS = 94; % AMT default
     
    t = (1:dur*fs)/fs; t = t(:); % creates 't' as a column array
    fc = 4000;
    dur_ramp_ms = 2.5;
    dur_ramp = round((dur_ramp_ms*1e-3)*fs); % duration ramp in samples
     
    insig = sin(2*pi*fc.*t);
    insig = scaletodbspl(insig,lvl,dBFS); % calibration before applying the ramp
     
    rp    = ones(size(insig)); 
    rp(1:dur_ramp) = rampup(dur_ramp);
    rp(end-dur_ramp+1:end) = rampdown(dur_ramp);
    insig = rp.*insig;
     
    insig = [zeros(50e-3*fs,1); insig; zeros(200e-3*fs,1)]; % 50 and 200 ms 
          % of silence before and after the sine tone
    
	ti_steady = 300*1e-3; % ms
	tf_steady = 340*1e-3; % ms
    
    for k = 1:N_models
        %%% Loading flags and keyvals
        [fg,kv] = il_get_flags(models{k});
        afb_flags = fg.afb_flags;
        ihc_flags = fg.ihc_flags;
        an_flags  = fg.an_flags;
        an_kv     = kv.an_keyvals;
        nomfb_flags = fg.nomfb_flags;
        %%%
        
        fname = ['fig08_' models{k} '-an-wave'];
        c = amt_cache('get',fname,flags.cachemode); 
        
        if ~isempty(c)
            bRun  = 0;
             
            outsig = c.outsig;
            fs_an  = c.fs_an;
            fc     = c.fc;
            if isfield(c,'out_psTH')
                out_psTH = c.out_psTH;
            end
        else
            bRun = 1;
            outsig = [];
        end
        
        YTL = []; % empty tick label
        if bRun
            amt_disp(['Calculating ' models{k} '...']);
            switch models{k}
                case 'dau1997'
                    fc_kv = {'flow',fc,'fhigh',fc,'basef',fc,'dboffset',dBFS};

                    outsig = dau1997(insig,fs,fc_kv{:},afb_flags{:}, ...
                        ihc_flags{:},an_flags{:},nomfb_flags{:});
                    fs_an = fs;
                    
                    c = [];
                    c.outsig = outsig;
                    c.fs_an  = fs_an;
                    c.fc     = fc;
                    amt_cache('set',fname,c);
                    
                case 'zilany2014'                      
                    kv={'fiberType',4,'numH',12,'numM',4,'numL',4,'nrep',nrep,'psth_binwidth',psth_binwidth};
                    [outsig,psth] = zilany2014(insig,fs,fc,kv{:});
                    fs_an = fs;
                    out_psTH.psth = psth;
                    out_psTH.psth_binwidth = psth_binwidth;
                    
                    c = [];
                    c.outsig = outsig;
                    c.fs_an  = fs_an;
                    c.fc     = fc;
                    c.out_psTH = out_psTH;
                    amt_cache('set',fname,c);
                    
                case 'verhulst2015'
                    fc_flag = fc; % one frequency will be simulated only
                    out = verhulst2015(insig,fs,fc_flag,an_flags{:},an_kv{:},nomfb_flags{:});
                    outsig = out(1).an_summed/19;
                    fs_an = 20000;
                    
                    c = [];
                    c.outsig = outsig;
                    c.fs_an  = fs_an;
                    c.fc     = fc;
                    amt_cache('set',fname,c);
                    
                case 'verhulst2018'
                    fc_flag = fc; % one frequency will be simulated only
                    out = verhulst2018(insig,fs,fc_flag,an_flags{:},an_kv{:},nomfb_flags{:});
                    outsig = out(1).an_summed/19;
                    fs_an = 20000;
                    
                    c = [];
                    c.outsig = outsig;
                    c.fs_an  = fs_an;
                    c.fc     = fc;
                    amt_cache('set',fname,c);

                case 'bruce2018'
                    kv = {'numH',12,'numM',4,'numL',4,'psthbinwidth_mr',psth_binwidth,'nrep',nrep};                 
                    out = bruce2018(insig,fs,fc,kv{:});
                    outsig = out.meanrate;                    
                    out_psTH.psth = out.psth; %neurogram_ft/nrep*250; %/2/(size(out.psth_ft,3));
                    out_psTH.psth_binwidth = psth_binwidth; %out.t_ft(2);
                    
                    fs_an = fs;
                    
                    c = [];
                    c.outsig = outsig;
                    c.fs_an  = fs_an;
                    c.fc     = fc;
                    c.out_psTH = out_psTH;
                    amt_cache('set',fname,c);
                    
                case 'king2019'
                    fc_kv = {'flow',fc,'fhigh',fc,'basef',fc,'dboffset',dBFS,'compression_n',0.3};

                    outsig = king2019(insig,fs,fc_kv{:},afb_flags{:}, ...
                        ihc_flags{:},an_flags{:},nomfb_flags{:});
                    fs_an = fs;
                    
                    c = [];
                    c.outsig = outsig;
                    c.fs_an  = fs_an;
                    c.fc     = fc;
                    amt_cache('set',fname,c);

                case 'relanoiborra2019'
                    fc_kv = {'flow',fc,'fhigh',fc,'basef',fc,'erbspacebw'};
                    [~,~,outsig] = relanoiborra2019_featureextraction(insig, fs,fc_kv{:});
                    fs_an = fs;
                    
                    c = [];
                    c.outsig = outsig;
                    c.fs_an  = fs_an;
                    c.fc     = fc;
                    amt_cache('set',fname,c);

                case 'osses2021'
                    fc_kv = {'flow',fc,'fhigh',fc,'basef',fc,'dboffset',dBFS};

                    outsig = osses2021(insig,fs,fc_kv{:},afb_flags{:}, ...
                        ihc_flags{:},an_flags{:},nomfb_flags{:});
                    fs_an = fs;
                    
                    c = [];
                    c.outsig = outsig;
                    c.fs_an  = fs_an;
                    c.fc     = fc;
                    amt_cache('set',fname,c);
            end
        end
        
        switch models{k}
            case {'zilany2014','bruce2018'}
                pst_hist = out_psTH.psth;
                dt = out_psTH.psth_binwidth;
                t_psth = (1:length(pst_hist))*dt;
        end
        
        t_an = (1:length(outsig))/fs_an; 
        t_an = t_an(:); % creates 't_an' as a column array
        
        outsig_all{k} = outsig;
        t_an_all{k}   = t_an;
                
        id_steady = find(t_an>=ti_steady & t_an<=tf_steady);
        
        onset  = max(outsig);
        steady = mean(outsig(id_steady));
        
        ra(k) = onset/steady;
        % fprintf('%s: Onset = %.1f (amplitude units), steady = %.1f %s, ratio = %.3f\n',models{k},onset,units_amplitude,steady,units_amplitude,ra(k));
    
        data_sim(k).onset = onset;
        data_sim(k).onset = steady;
        fs_an_all(k) = fs_an;
    end
     
    if flags.do_plot
        if ~exist('idxs','var'), idxs=1:8; end
        for k = 1:N_models
            
            switch idxs(k)
                case {1,7,8}
                    figure(1); 
                    handle_data(k) = plot(t_an_all{k}*1000,outsig_all{k},'Color',Colours{k},'LineWidth',2); hold on, grid on;
                    
                    if k == 1
                        figure_handle(end+1) = gcf; % multiple figures will be generated
                        figure_name{end+1}   = ['fig08-tone-4-kHz-' models{k}];
                    end

                    units_amplitude = '(MU)';
                    YL = [-300 1500]; % a priori knowledge
                    stepY = 100;
                    YT = YL(1)+stepY:stepY:YL(2)-stepY;
                    
                    if k == 1
                        YTL = [];
                        for ii = 1:length(YT)
                            if mod(YT(ii),200)==0
                                YTL{ii} = num2str(YT(ii));
                            else
                                YTL{ii} = '';
                            end
                        end
                    end

                case {2,5}
                    figure;
                    
                    handle_data(k) = stairs(t_psth*1000,pst_hist,'Color',0.5*Colours{k},'LineWidth',2); hold on, grid on;
                    plot(t_an_all{k}*1000,outsig_all{k},'Color',Colours{k},'LineWidth',2); hold on, grid on;
                    % if k == 2
                    figure_handle(end+1) = gcf; % multiple figures will be generated
                    figure_name{end+1}   = ['fig08-tone-4-kHz-' models{k}];
                    % end
                    units_amplitude = '(spikes/s)';
                    YL = [0 1000]; % a priori knowledge
                    stepY = 100;
                    YT = YL(1)+stepY:stepY:YL(2)-stepY;

                case {3,4}
                    figure(3)
                    handle_data(k) = plot(t_an_all{k}*1000,outsig_all{k},'Color',Colours{k},'LineWidth',2); hold on, grid on;
                    if k == 3
                        figure_handle(end+1) = gcf; % multiple figures will be generated
                        figure_name{end+1}   = ['fig08-tone-4-kHz-' models{k}];
                    end
                    
                    units_amplitude = '(spikes/s)';
                    YL = [0 1000]; % a priori knowledge
                    stepY = 100;
                    YT = YL(1)+stepY:stepY:YL(2)-stepY;

                case 6
                    factor = 1e3;
                    figure;
                    handle_data(k) = plot(t_an_all{k}*1000,factor*outsig_all{k},'Color',Colours{k},'LineWidth',2); hold on, grid on;
                    if k == 6
                        figure_handle(end+1) = gcf; % multiple figures will be generated
                        figure_name{end+1}   = ['fig08-tone-4-kHz-' models{k}];
                    end
                    
                    units_amplitude = '(arbitrary units x 10^{-3})';
                    YL = factor*1.8e-3*[-1 1]; % a priori knowledge    
                    stepY = factor*.3e-3;
                    YT = YL(1)+stepY:stepY:YL(2)-stepY;
                    % for i = 1:length(YT)
                    %     YTL{end+1} = sprintf('%.4f',YT(i));
                    % end
            end
            xlabel('Time (ms)');
            ylabel(['Amplitude ' units_amplitude])

            set(gca,'XTick',50:50:450);
            xlim([25 475]);
            
            ylim(YL); % ([-300 1480])
            set(gca,'YTick',YT);
            switch k
                case {1,7,8}
                    set(gca,'YTickLabel',YTL);
            end
                        
            offy = 5;

            Xhere = 1000*[t_an(min(id_steady)) t_an(max(id_steady))];
            Yhere = (YL(1)+offy)*[1 1];
            plot(Xhere,Yhere,'--','LineWidth',2,'Color',[0.5 0.5 0.5]);
            
            Yhere = (YL(2)-offy)*[1 1];
            plot(Xhere,Yhere,'--','LineWidth',2,'Color',[0.5 0.5 0.5]);
            
            Yhere = YL;
            Xhere = 1000*t_an(min(id_steady))*[1 1];
            plot(Xhere,Yhere,'--','LineWidth',2,'Color',[0.5 0.5 0.5]);
            
            Xhere = 1000*t_an(max(id_steady))*[1 1];
            plot(Xhere,Yhere,'--','LineWidth',2,'Color',[0.5 0.5 0.5]);

            Pos = get(gcf,'Position');
            Pos(3:4) = [420 250]; % [320 360];
            set(gcf,'Position',Pos);
        end
    end
    
    data.figure_flag = 'do_fig8';
    data.ra = ra;
    
    data.insig = insig;
    data.fs    = fs;
    data.data_sim = data_sim;
    data.models = models;
    data.figure_handle = figure_handle;
    data.handle_data   = handle_data;
    data.figure_name   = figure_name;  
    data.outsig_all = outsig_all;
    data.fs_an_all  = fs_an_all;
end
%% FIG 9
if flags.do_fig9 
    % Code adapted from testPopRateLevel_BEZ2018_mine.m
    
    %%% 1. Generating the input signals: %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%   
    % Signal parameters:
    CF    = 4000; %8e3;  % CF in Hz;
    lvls = 0:10:100;
    N_lvls = length(lvls);
    p0 = 2e-5; % reference pressure

    % dur     = 50e-3; % stimulus duration in seconds
    dur     = 300e-3; % stimulus duration in seconds
    % dur4sim = 2*dur; % duration to simulate, double as long as the input stimulus length
    dur_ramp_ms = 2.5; % rise/fall time in seconds
    dur_ramp = round((dur_ramp_ms*1e-3)*fs); % duration ramp in samples
    
    t = 0:1/fs:dur-1/fs; % time vector

    insig_orig = sin(2*pi*CF*t'); % column array
    rp    = ones(size(insig_orig)); 
    rp(1:dur_ramp) = rampup(dur_ramp);
    rp(end-dur_ramp+1:end) = rampdown(dur_ramp);
        
    for kk = 1:N_lvls
        lvl = lvls(kk);
        insig(:,kk) = scaletodbspl(insig_orig,lvl,dBFS); % calibrated, unramped stimulus
    end 
    insig = repmat(rp,1,N_lvls).*insig;
    %%% End generating the signal (the level adjustment is done later)

    numH = 12; % Automate this
    numM =  4;
    numL =  4;
    numTot = numH+numM+numL;
    numCum = cumsum([numL numM numH]);
    rates_avg = nan(numTot,N_lvls,N_models);
    rates     = nan(numTot,N_lvls,N_models);
    
    for k = 1:N_models
        
        psTH_H = [];
        psTH_M = [];
        psTH_L = [];
        
        %%% Loading flags and keyvals
        [fg,kv] = il_get_flags(models{k});
        afb_flags = fg.afb_flags;
        ihc_flags = fg.ihc_flags;
        an_flags = fg.an_flags;
        nomfb_flags = fg.nomfb_flags;
        an_kv    = kv.an_keyvals;
        %%%
        
        fname = ['fig09_rate-level-' models{k}];
        c = amt_cache('get',fname,flags.cachemode); 
         
        if ~isempty(c)
            bRun  = 0;
        
            psTH_H = c.psTH_H;
            psTH_M = c.psTH_M;
            psTH_L = c.psTH_L;
            rates  = c.rates;
            fs_an  = c.fs_an;
            idxi   = c.idxi;
            idxf   = c.idxf;
        else
            bRun = 1;
            % outsig = [];
        end
        
        if bRun
            amt_disp(['Calculating ' models{k} '...']);
            switch models{k}
                case 'dau1997'
                    fc_kv = {'flow',CF,'fhigh',CF,'basef',CF,'dboffset',dBFS};

                    for kk = 1:N_lvls
                        psTH_H(:,kk) = dau1997(insig(:,kk),fs,fc_kv{:}, ...
                            afb_flags{:},ihc_flags{:},an_flags{:},nomfb_flags{:});
                    end
                    rates = [];
                    psTH_M = nan(size(psTH_H));
                    psTH_L = nan(size(psTH_H));
                    fs_an = fs;
                    
                    ronset  = round(15e-3*fs_an)+1; % start after 15 ms (strong onset)
                    roffset = round(dur*fs_an); % end 
                    idxi = ronset;
                    idxf = roffset;
        
                case 'zilany2014'
                    for kk = 1:N_lvls
                        kv={'fiberType',4,'numH',0,'numM',0,'numL',numL,'nrep',10};
                        psTH_L(:,kk) = zilany2014(insig(:,kk),fs,CF,kv{:});
                        kv={'fiberType',4,'numH',0,'numM',numM,'numL',0,'nrep',10};
                        psTH_M(:,kk) = zilany2014(insig(:,kk),fs,CF,kv{:});
                        kv={'fiberType',4,'numH',numH,'numM',0,'numL',0,'nrep',10};
                        psTH_H(:,kk) = zilany2014(insig(:,kk),fs,CF,kv{:});
                        kv={'fiberType',4,'numH',numH,'numM',numM,'numL',numL,'nrep',10};
                        rates(:,kk,k) = zilany2014(insig(:,kk),fs,CF,kv{:});
                    end
                    fs_an = fs;

                    ronset  = round(15e-3*fs_an)+1; % start after 15 ms (strong onset)
                    roffset = round(dur*fs_an); % end 
                    idxi = ronset;
                    idxf = roffset;
                    
                case 'verhulst2015' 
                    fc_flag = CF;
                    out = verhulst2015(insig,fs,fc_flag,an_flags{:},an_kv{:},nomfb_flags{:});

                    for kk = 1:N_lvls
                      psTH_L(:,kk) = out(kk).anfL;
                      psTH_M(:,kk) = out(kk).anfM;
                      psTH_H(:,kk) = out(kk).anfH;
                    end
                    rates = [];
                    fs_an = out(1).fs_an;

                    ronset  = round(15e-3*fs_an)+1; % start after 15 ms (strong onset)
                    roffset = round(dur*fs_an); % end 
                    idxi = ronset;
                    idxf = roffset;

                case 'verhulst2018' 
                    fc_flag = CF;
                    out = verhulst2018(insig,fs,fc_flag,an_flags{:},an_kv{:},nomfb_flags{:});

                    for kk = 1:N_lvls
                        psTH_L(:,kk) = out(kk).anfL;
                        psTH_M(:,kk) = out(kk).anfM;
                        psTH_H(:,kk) = out(kk).anfH;
                    end
                    rates = [];
                    fs_an = out(1).fs_an;
                    
                    ronset  = round(15e-3*fs_an)+1; % start after 15 ms (strong onset)
                    roffset = round(dur*fs_an); % end 
                    idxi = ronset;
                    idxf = roffset;
                    
                case 'bruce2018'
                    kv = {'numH',numH,'numM',numM,'numL',numL,'nrep',10,'outputPerSynapse'}; 
                    rates = [];
                    for kk = 1:N_lvls
                      out = bruce2018(insig(:,kk),fs,CF,kv{:});
                      psTH_L(:,kk)=squeeze(mean(out.meanrate(:,1,1:numL),3))';
                      psTH_M(:,kk)=squeeze(mean(out.meanrate(:,1,numL+1:numL+numM),3))';
                      psTH_H(:,kk)=squeeze(mean(out.meanrate(:,1,numL+numM+1:end),3))';
                      rates(:,kk,k)= squeeze(mean(out.meanrate(idxi:idxf,1,:),1));
                    end
                    fs_an = fs;

                    ronset  = round(15e-3*fs_an)+1; % start after 15 ms (strong onset)
                    roffset = round(dur*fs_an); % end 
                    idxi = ronset;
                    idxf = roffset;
                
                case {'king2019'}
                    fc_kv = {'flow',CF,'fhigh',CF,'basef',CF,'dboffset',dBFS};

                    for kk = 1:N_lvls
                        psTH_H(:,kk) = king2019(insig(:,kk),fs,fc_kv{:}, ...
                            afb_flags{:},ihc_flags{:},an_flags{:},nomfb_flags{:});
                    end
                    rates = [];
                    psTH_M = nan(size(psTH_H));
                    psTH_L = nan(size(psTH_H));
                    fs_an = fs;
                    
                    ronset  = round(15e-3*fs_an)+1; % start after 15 ms (strong onset)
                    roffset = round(dur*fs_an); % end 
                    idxi = ronset;
                    idxf = roffset;
                    
                case {'relanoiborra2019'}
                    fc_kv = {'flow',CF,'fhigh',CF,'basef',CF,'erbspacebw','no_internalnoise'};

                    for kk = 1:N_lvls
                        [~,~,psTH_H(:,kk)] = relanoiborra2019_featureextraction(insig(:,kk), fs,fc_kv{:});
                    end
                    rates = [];
                    psTH_M = nan(size(psTH_H));
                    psTH_L = nan(size(psTH_H));
                    fs_an = fs;
                    
                    ronset  = round(15e-3*fs_an)+1; % start after 15 ms (strong onset)
                    roffset = round(dur*fs_an); % end 
                    idxi = ronset;
                    idxf = roffset;
                    
                case 'osses2021'
                    fc_kv = {'flow',CF,'fhigh',CF,'basef',CF,'dboffset',dBFS};

                    for kk = 1:N_lvls
                        psTH_H(:,kk) = osses2021(insig(:,kk),fs,fc_kv{:}, ...
                            afb_flags{:},ihc_flags{:},an_flags{:},nomfb_flags{:});
                    end
                    rates = [];
                    psTH_M = nan(size(psTH_H));
                    psTH_L = nan(size(psTH_H));
                    fs_an = fs;
                    
                    ronset  = round(15e-3*fs_an)+1; % start after 15 ms (strong onset)
                    roffset = round(dur*fs_an); % end 
                    idxi = ronset;
                    idxf = roffset;
        
            end
            
            c = [];
            c.psTH_H = psTH_H; c.psTH_M = psTH_M; c.psTH_L = psTH_L;
            c.rates = rates; 
            c.fs_an = fs_an;
            c.idxi = idxi;
            c.idxf = idxf;
            amt_cache('set',fname,c);
        end
        
        for kk = 1:N_lvls
            rates_avg(          1:numCum(1),kk,k)= mean(psTH_L(idxi:idxf,kk));
            rates_avg(numCum(1)+1:numCum(2),kk,k)= mean(psTH_M(idxi:idxf,kk));
            rates_avg(numCum(2)+1:end      ,kk,k)= mean(psTH_H(idxi:idxf,kk));
        end
                
        if flags.do_plot
            
            Pos34 = [300 300]; % width and height of the resulting plot
            factor = 1;
            switch models{k}
                case {'dau1997','relanoiborra2019','osses2021'}
                    figure(1);
                    if k == 1
                        figure_handle(end+1) = gcf;
                        figure_name{end+1} = ['fig09-AN-firing-rates-model-' models{k}];
                    end
                    YL = 'Amplitude (MU)';
                    YLim = [-10 190];
                    YT = 0:20:180;
                    
                case {'zilany2014','bruce2018'}
                    figure;
                    YL = 'Firing rate (Spikes/s)';
                    figure_handle(end+1) = gcf;
                    figure_name{end+1} = ['fig09-AN-firing-rates-model-' models{k}];
                    if strcmp(models{k},'zilany2014')
                        YLim = [-20 380];
                        YT = 0:30:360;
                    else
                        YLim = [-20 320];
                        YT = 0:30:300;
                    end
                    
                case {'verhulst2015','verhulst2018'}
                    figure
                    
                    figure_handle(end+1) = gcf;
                    figure_name{end+1} = ['fig09-AN-firing-rates-model-' models{k}];
                    
                    YL = 'Firing rate (Spikes/s)';
                    YLim = [-20 320];
                    YT = 0:30:300;
                    
                case 'king2019'
                    factor = 1e4;
                    
                    figure;
                    
                    if factor == 1e4
                        YL = 'Amplitude (a.u. x 10^{-4})';
                    else
                        YL = 'Amplitude (a.u.)';
                    end
                    figure_handle(end+1) = gcf;
                    figure_name{end+1} = ['fig09-AN-firing-rates-model-' models{k}];
                    YLim = 1e-4*[-0.2 5.2]*factor;
                    YT = (0:.5:5)*1e-4*factor;
                    
            end
            
            switch models{k}
                case {'zilany2014','bruce2018'}
                plot(lvls,factor*squeeze(rates(1: 3,:,k)),'-' ,'Color',il_rgb('Gray')); hold on, grid on
                plot(lvls,factor*squeeze(rates(4: 6,:,k)),'--','Color',il_rgb('Gray'));
                plot(lvls,factor*squeeze(rates(7:19,:,k)),'-' ,'Color',il_rgb('Gray'));
            end
            
            hl1 = plot(lvls,factor*squeeze(rates_avg(1:numCum(1),:,k)),'o-','Color',Colours{k},'MarkerFaceColor','w','LineWidth',3,'MarkerSize',10,'MarkerFaceColor','w'); hold on, grid on
            hl2 = plot(lvls,factor*squeeze(rates_avg(numCum(1)+1:numCum(2),:,k)),'s--','Color',Colours{k},'LineWidth',3,'MarkerSize',10,'MarkerFaceColor','w');
            hl3 = plot(lvls,factor*squeeze(rates_avg(numCum(2)+1:numCum(3),:,k)),'d-','Color',Colours{k},'LineWidth',3,'MarkerSize',10,'MarkerFaceColor',Colours{k});
            
            handle_dataL(k) = hl1(1);
            handle_dataM(k) = hl2(1);
            handle_data(k) = hl3(1);
            
            rates_avg_all{k} = squeeze(rates_avg(7:19,:,k));
            if k == 1
                rates_avg_description = 'rate levels for ''HSR''';
            end
        end
        
        xlim([min(lvls)-2 max(lvls)+2])
        ylim(YLim)
        
        xlabel('Stimulus Level (dB SPL)')
        ylabel(YL)
        
        set(gca,'XTick',0:10:100);
        set(gca,'YTick',YT);
        
        Pos = get(gcf,'Position');
        Pos(3:4) = Pos34;
        set(gcf,'Position',Pos);
    end
    
    data.figure_flag   = 'do_fig9';
    data.figure_handle = figure_handle;
    data.figure_name   = figure_name;    
    data.handle_data   = handle_data;
    data.handle_dataL  = handle_dataL; % handle for LSR data
    data.handle_dataM  = handle_dataM; % handle for MSR data
    data.models        = models;
    data.lvls = lvls;
    data.rate_avg_all  = rates_avg_all;
    data.rate_avg_description = rates_avg_description;
end

%% ------ Fig 10 ---------------------------------------------------------
if flags.do_fig10

    error('Not implemented yet. Sorry!'); 
    
    
    % 1. Stimulus creation:
    dur = 50e-3; % 300 ms
    lvls = 0:10:100;
    N_lvls = length(lvls);
    psth_binwidth = 0.5e-3; % width of the PSTH bin, in seconds
    nrep=100;
    numH = 12; 
    numM =  4;
    numL =  4;

    t = (1:dur*fs)/fs; t = t(:); % creates 't' as a column array
    fc = 4000;
    dur_ramp_ms = 2.5;
    dur_ramp = round((dur_ramp_ms*1e-3)*fs); % duration ramp in samples
     
    insig_orig = sin(2*pi*fc.*t);
    insig = zeros(length(t),N_lvls); % memory allocation
    
    for j = 1:N_lvls
        lvl = lvls(j);
        insig(:,j) = scaletodbspl(insig_orig,lvl,dBFS); % calibration before applying the ramp
    end
    
    rp    = ones(size(insig,1),1); 
    rp(1:dur_ramp) = rampup(dur_ramp);
    rp(end-dur_ramp+1:end) = rampdown(dur_ramp);
    
    insig = repmat(rp,1,N_lvls).*insig;
    
    insig = [zeros(50e-3*fs,N_lvls); insig]; % 50 ms of silence before and after the sine tone
     
    ti_steady = 300*1e-3; % ms
    tf_steady = 340*1e-3; % ms
     
    for k = 1:N_models
        %%% Loading flags and keyvals
        [fg,kv]     = il_get_flags(models{k});
        afb_flags   = fg.afb_flags;
        ihc_flags   = fg.ihc_flags;
        an_flags    = fg.an_flags;
        an_kv       = kv.an_keyvals;
        nomfb_flags = fg.nomfb_flags;
        %%%
     
        fname = ['fig10_' models{k} '-adaptation-onset-steady'];
        c = amt_cache('get',fname,flags.cachemode); 
         
        if ~isempty(c)
            bRun  = 0;
             
            outsig = c.outsig;
            fs_an  = c.fs_an;
            units_amplitude = c.units_amplitude;
            if isfield(c,'out_psTH')
                out_psTH = c.out_psTH;
            end
        else
            bRun = 1;
            outsig = [];
            amt_disp(['Calculating ' models{k} '...']);
        end
        
        factor = 1; % later used for the plots
        switch models{k}
            
            case {'dau1997','osses2021'}
                if bRun
                    fc_kv = {'flow',fc,'fhigh',fc,'basef',fc,'dboffset',dBFS};

                    for j = 1:N_lvls
                        outsig(:,j) = dau1997(insig(:,j),fs,fc_kv{:}, ...
                            afb_flags{:},ihc_flags{:},an_flags{:},nomfb_flags{:});
                    end
                    fs_an = fs;
%                     adapt_min = -230.9; % MU, Osses2018, Appendix C, page 182
                    % outsig = outsig - adapt_min;

                    units_amplitude = '(Model Units)';
                    
                    c = [];
                    c.fs_an   = fs_an;
                    c.outsig  = outsig;
                    c.units_amplitude = units_amplitude;
                    amt_cache('set',fname,c);
                end
     
            case 'zilany2014'
                if bRun
                    kv={'fiberType',4,'numH',numH,'numM',numM,'numL',numL,'nrep',nrep,'psth_binwidth',psth_binwidth};
                    for j = 1:N_lvls
                        [outsig(:,j),out_psTH.psth(:,j)] = zilany2014(insig(:,j),fs,fc,kv{:});
                        if j == 1
                            out_psTH.psth_binwidth = psth_binwidth;
                        end
                    end

                    fs_an = fs;
                    units_amplitude = '(spikes/s)';
                    
                    c = [];
                    c.fs_an   = fs_an;
                    c.outsig  = outsig;
                    c.units_amplitude = units_amplitude;
                    c.out_psTH = out_psTH;
                    amt_cache('set',fname,c);
                end    
    
            case 'verhulst2015'
                if bRun
                    fc_flag = fc;
                    out = verhulst2015(insig,fs,fc_flag,an_flags{:},an_kv{:},nomfb_flags{:});

                    N_tot = out(1).keyvals.numH+out(1).keyvals.numM+out(1).keyvals.numL;
                    for j = 1:N_lvls
                        outsig(:,j) = out(j).an_summed/N_tot;
                    end
                    fs_an = 20000;
                    units_amplitude = '(spikes/s)';
                    
                    c = [];
                    c.fs_an   = fs_an;
                    c.outsig  = outsig;
                    c.units_amplitude = units_amplitude;
                    amt_cache('set',fname,c);
                end
                
            case 'verhulst2018'
                if bRun
                    fc_flag = fc;
                    out = verhulst2018(insig,fs,fc_flag,an_flags{:},an_kv{:},nomfb_flags{:});
                    N_tot = out(1).keyvals.numH+out(1).keyvals.numM+out(1).keyvals.numL;
                    for j = 1:N_lvls
                        outsig(:,j) = out(j).an_summed/N_tot;
                    end
                    fs_an = 20000;
                    units_amplitude = '(spikes/s)';
                    
                    c = [];
                    c.fs_an   = fs_an;
                    c.outsig  = outsig;
                    c.units_amplitude = units_amplitude;
                    amt_cache('set',fname,c);
                end
                
            case 'bruce2018'
                if bRun                    
                    kv = {'numH',numH,'numM',numM,'numL',numL,'psthbinwidth_mr',psth_binwidth,'nrep',nrep};
                    outsig=[]; out_psTH=[];
                    for j = 1:N_lvls
                        out = bruce2018(insig(:,j),fs,fc,kv{:});
                        outsig(:,j) = out.meanrate;
                        out_psTH.psth(:,j) = out.psth;
                    end
                    out_psTH.psth_binwidth = psth_binwidth;                  
                    fs_an = fs;
                    units_amplitude = '(spikes/s)';
                    
                    c = [];
                    c.fs_an   = fs_an;
                    c.outsig  = outsig;
                    c.units_amplitude = units_amplitude;
                    c.out_psTH = out_psTH;
                    amt_cache('set',fname,c);
                end    
                
            case 'king2019'
                if bRun
                    fc_kv = {'flow',fc,'fhigh',fc,'basef',fc,'dboffset',dBFS};

                    for j = 1:N_lvls
                        outsig(:,j) = king2019(insig(:,j),fs,fc_kv{:},'compression_n',0.3, ...
                            afb_flags{:},ihc_flags{:},an_flags{:},nomfb_flags{:});
                    end
                    fs_an = fs;
                    
                    c = [];
                    c.fs_an   = fs_an;
                    c.outsig  = outsig;
                    c.units_amplitude = units_amplitude;
                    amt_cache('set',fname,c);
                end 
                factor = 1e3;
                if factor == 1e3
                    units_amplitude = '(a.u. x 10^{-3})';
                else
                    units_amplitude = '(a.u.)';
                end
               
            case 'relanoiborra2019'
                if bRun
                    fc_kv = {'flow',fc,'fhigh',fc,'basef',fc,'erbspacebw'};

                    for j = 1:N_lvls
                      [~,~,outsig(:,j)] = relanoiborra2019_featureextraction(insig(:,j), fs,fc_kv{:});
                    end
                    fs_an = fs;
                    
                    units_amplitude = '(a.u.)';
                    
                    c = [];
                    c.fs_an   = fs_an;
                    c.outsig  = outsig;
                    c.units_amplitude = units_amplitude;
                    amt_cache('set',fname,c);
                end  
    
        end
        t_an = (1:length(outsig))/fs_an; 
        t_an = t_an(:); % creates 't_an' as a column array
        id_steady = find(t_an>=ti_steady & t_an<=tf_steady);
     
        onset_here  = max(outsig);
        steady_here = mean(outsig(id_steady,:),1);
        steady_min  = min(outsig(id_steady,:));
        steady_max  = max(outsig(id_steady,:));
        
        onset(k,1:N_lvls)  = onset_here;
        steady(k,1:N_lvls) = steady_here;
    
        ra(k,1:N_lvls) = onset_here./steady_here;
        
        % fprintf('%s: Onset = %.1f %s, steady = %.1f %s, ratio = %.3f\n',models{k},onset,units_amplitude,steady,units_amplitude,ra(k));
    
        if flags.do_plot
     
            XT = 0:10:100;
            Pos34 = [300 300];
                        
            switch models{k}
                case {'dau1997','relanoiborra2019','osses2021'} % dau1997
                    figure(1);
                    
                    if k == 1
                        figure_handle(end+1) = gcf; % multiple figures will be generated
                        figure_name{end+1}   = sprintf('fig10-tone-4-kHz-IO-onset-%s',models{k});
                    end
                    YL = [-50 1650];
                    stepY = 100;
                    YT = 0:stepY:1600;
                    
                    if k == 1
                        YTL = [];
                        for ii = 1:length(YT)
                            if mod(YT(ii),200)==0
                                YTL{ii} = num2str(YT(ii));
                            else
                                YTL{ii} = '';
                            end
                        end
                    end
                    
                case {'zilany2014','bruce2018'} % zilany and bruce
                    figure;
                    
                    figure_handle(end+1) = gcf; % multiple figures will be generated
                    figure_name{end+1}   = sprintf('fig10-tone-4-kHz-IO-onset-%s',models{k});
                    
                    YL = [-50 1650];
                    stepY = 100;
                    YT = 0:stepY:1600;
                    
                case {'verhulst2015','verhulst2018'}
                    figure;
                    
                    figure_handle(end+1) = gcf; % multiple figures will be generated
                    figure_name{end+1}   = sprintf('fig10-tone-4-kHz-IO-onset-%s',models{k});
                    
                    YL = [-50 1650];
                    stepY = 100;
                    YT = 0:stepY:1600;
                    
                case 'king2019'
                    figure;
                    figure_handle(end+1) = gcf; % multiple figures will be generated
                    figure_name{end+1}   = sprintf('fig10-tone-4-kHz-IO-onset-%s',models{k});
                    
                    YL = [-0.2 5.2]*1e-3*factor;
                    YT = (0:.5:5)*1e-3*factor;
                    
                    
            end
            plot(lvls,factor*onset_here,'Color',Colours{k},'Marker',Markers{k},'LineWidth',2,'MarkerFaceColor',Colours{k}); 
            hold on, grid on;
    
            switch models{k}
                case {'zilany2014','bruce2018'}
                    onset_psTH  = max(out_psTH.psth);
                    
                    plot(lvls,factor*onset_psTH,'--','Color',0.5*Colours{k},'Marker',Markers{k},'LineWidth',2,'MarkerFaceColor','w'); 
                    hold on, grid on;
            end
            
            xlabel('Input level (dB)');
            ylabel(sprintf('Onset amplitude %s',units_amplitude));
            
            Pos = get(gcf,'Position');
            Pos(3:4) = Pos34;
            set(gcf,'Position',Pos);
            
            xlim([-3 103]);
            set(gca,'XTick',XT);
            
            ylim(YL);
            set(gca,'YTick',YT);
            
            if ~strcmp(models{k},'king2019')
                % For all models except King that has another scaling
                set(gca,'YTickLabel',YTL);
            end
            
        end
    end
     
    data.figure_flag = 'do_fig10';
    data.ra = ra;
    data.onset = onset;
    data.steady = steady;
    data.figure_handle = figure_handle;
    data.figure_name   = figure_name; 
    data.models        = models;
end

%% ------ FIG 11 ------------------------
if flags.do_fig11
    
    figure_handle_ax = []; 
    % Changes on 12/01/2021:
    %       - Verhulst2015 was not using the same SR, now it is...
    %       - The models now use fc_model=4012 Hz instead of 4000 Hz
    
    % This is similar to Verhulst et al. 2018, Fig. 3C:
    %%% 1. Generating the input signals: %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%   
    % Signal parameters:
    fc_target = 4000; % Hz, frequency of the carrier
    cf = il_m2hz(401); % Get characteristic frequencies from Verhulst models
    idx_cf = find(cf>fc_target,1,'last');
    fc = cf(idx_cf); 
        
    fmod     =  100; % Hz, frequency of the modulator
    mdepth   = 1; % value between 0 and 1
    dur      = 500e-3; % stimulus duration in seconds
    dur_ramp = 2.5e-3; % s, duration of the ramp
    psth_binwidth = 0.5e-3; % width of the PSTH bin, in seconds
    nrep=100;
    numH = 12; 
    numM =  4;
    numL =  4;

    lvl = 60; % level, dB
    
    N_samples = round(dur*fs);
    dur_ramp_samples = round((dur_ramp)*fs);

    % Creating a cosine ramp:
    ramp = ones(N_samples,1);
    ramp(1:dur_ramp_samples)         = rampup(dur_ramp_samples);
    ramp(end-dur_ramp_samples+1:end) = rampdown(dur_ramp_samples);

    % AM stimulus and calibration:
    t = (0:N_samples-1)/fs; t=t(:);
    carrier = sin(2*pi*fc*t); % starts at phase = 0
    env = (1 + mdepth * sin(2*pi*fmod*t-pi/2) ); % modulator starts at minimum (phase=-pi/2)
    insig = env .* carrier; % Amplitude-modulated signal
    insig = scaletodbspl(insig,lvl,dBFS);
    insig = ramp.*insig;

    %%%
    psthbinwidth = 0.5e-3; % 0.75 ms
    percent = 90; % percentage overlap
    %%%
    outs_anf = [];
    
    for k = 1:N_models
        %%% Loading flags and keyvals:
        [fg,kv] = il_get_flags(models{k});
        fc_flags   = fg.fc_flags;
        afb_flags  = fg.afb_flags;
        ihc_flags  = fg.ihc_flags;
        an_flags   = fg.an_flags; % No auditory nerve module
        nomfb_flags= fg.nomfb_flags;
        an_kv = kv.an_keyvals;
        
        fname = ['fig11_' models{k} '-an'];
        c = amt_cache('get',fname,flags.cachemode); 
        
        if ~isempty(c)
            bRun  = 0;
            
            an_summed = c.an_summed;
            fs_an     = c.fs_an;
            anfH      = c.anfH;
            anfM      = c.anfM;
            anfL      = c.anfL;
            if isfield(c,'out_psTH')
                out_psTH = c.out_psTH;
            end
        else
            bRun = 1;
            an_summed = [];
            amt_disp(['Calculating ' models{k} '...']);
        end
        
        %%% Figure settings (overloaded later for dau1997):
        text4ylabel='(spikes/s)';
        YL = [-48 330];
        factor = 1;
        %%%
        switch models{k}
            case {'dau1997','osses2021'}
                if bRun
                    fc_kv = {'flow',fc,'fhigh',fc,'basef',fc,'dboffset',dBFS};
                    out = dau1997(insig,fs,fc_kv{:},afb_flags{:}, ...
                        ihc_flags{:},an_flags{:},nomfb_flags{:});
                    fs_an = fs;
                    an_summed = out;
                    anfH = out;
                    anfM = 0*out;
                    anfL = 0*out;
                    
                    c = [];
                    c.an_summed = an_summed;
                    c.fs_an = fs_an;
                    c.anfH  = anfH;
                    c.anfM  = anfM;
                    c.anfL  = anfL;
                    amt_cache('set',fname,c);
                end
                    
                text4ylabel='(MU)';
                YL = [-230 330];
                
            case {'zilany2014'}
                if bRun
                    kv={'fiberType',4,'numH',numH,'numM',numM,'numL',numL,'nrep',nrep,'psth_binwidth',psth_binwidth};

                    [an_summed,out_psTH.psth] = zilany2014(insig,fs,fc,kv{:});
                    fs_an = fs;
                    anfH = an_summed; 
                    anfM = an_summed;
                    anfL = an_summed;
                    out_psTH.psth_binwidth = psth_binwidth;

                    c = [];
                    c.an_summed = an_summed;
                    c.fs_an = fs_an;
                    c.anfH  = anfH;
                    c.anfM  = anfM;
                    c.anfL  = anfL;
                    c.out_psTH = out_psTH;
                    amt_cache('set',fname,c);
                end
                
            case {'verhulst2015'}
                if bRun
                    hear_profile = 'Flat00'; % NH audiogram, default
                    fc_kv = {'hearing_profile',hear_profile};
                    out = verhulst2015(insig,fs,fc_flags{:},fc_kv{:},an_kv{:},an_flags{:},nomfb_flags{:});
                    fs_an = out.fs_abr;

                    if strcmp(an_kv{1},'numH')
                        numH = an_kv{2};
                    end
                    if strcmp(an_kv{3},'numM')
                        numM = an_kv{4};
                    end
                    if strcmp(an_kv{5},'numL')
                        numL = an_kv{6};
                    end
                    
                    num_tot = numH+numM+numL;
                    an_summed = out.an_summed(:,idx_cf)/num_tot;
                    anfH = out.anfH;
                    anfM = out.anfM;
                    anfL = out.anfL;
                    
                    c = [];
                    c.an_summed = an_summed;
                    c.fs_an = fs_an;
                    c.anfH  = anfH;
                    c.anfM  = anfM;
                    c.anfL  = anfL;
                    amt_cache('set',fname,c);
                end
                
            case {'verhulst2018'}
                if bRun
                    %%% Model parameters (only one hearing profile, Flat00 and 13-3-3):
                    hear_profile = 'Flat00';
                    fc_kv = {'hearing_profile',hear_profile};
                    out = verhulst2018(insig,fs,fc_flags{:},fc_kv{:},an_kv{:},an_flags{:},nomfb_flags{:});
                    fs_an = out.fs_abr;

                    if strcmp(an_kv{1},'numH')
                        numH = an_kv{2};
                    end
                    if strcmp(an_kv{3},'numM')
                        numM = an_kv{4};
                    end
                    if strcmp(an_kv{5},'numL')
                        numL = an_kv{6};
                    end
                    
                    num_tot = numH+numM+numL;
                    an_summed = out.an_summed(:,idx_cf)/num_tot;
                    anfH = out.anfH;
                    anfM = out.anfM;
                    anfL = out.anfL;
                    
                    c = [];
                    c.an_summed = an_summed;
                    c.fs_an = fs_an;
                    c.anfH  = anfH;
                    c.anfM  = anfM;
                    c.anfL  = anfL;
                    amt_cache('set',fname,c);
                end
                
            case {'bruce2018'}
                if bRun
                    kv = {'numH',numH,'numM',numM,'numL',numL,'psthbinwidth_mr',psth_binwidth,'nrep',nrep}; 
                    out = bruce2018(insig,fs,fc,kv{:});
                    an_summed = out.meanrate;
                    anfH=[]; anfM = []; anfL = [];
                    out_psTH.psth = out.psth;
                    out_psTH.psth_binwidth = psth_binwidth; 
                    fs_an = fs;
                
                    c = [];
                    c.fs_an = fs_an;
                    c.anfH  = anfH;
                    c.anfM  = anfM;
                    c.anfL  = anfL;
                    c.out_psTH = out_psTH;
                    amt_cache('set',fname,c);
                end
                
            case 'king2019'
                if bRun
                    fc_kv = {'flow',fc,'fhigh',fc,'basef',fc,'dboffset',dBFS};
                    out = king2019(insig,fs,fc_kv{:},afb_flags{:}, ...
                        ihc_flags{:},an_flags{:},nomfb_flags{:});
                    fs_an = fs;
                    an_summed = out;
                    anfH = out;
                    anfM = 0*out;
                    anfL = 0*out;
                    
                    c = [];
                    c.fs_an = fs_an;
                    c.an_summed = an_summed;
                    c.anfH  = anfH;
                    c.anfM  = anfM;
                    c.anfL  = anfL;
                    amt_cache('set',fname,c);
                end
                    
                factor = 1000;
                if factor == 1000
                    text4ylabel='(a.u. x 10^-3)';
                else
                    text4ylabel='(a.u.)';
                end
                YL = 1e-3*[-1 1];
                
            case 'relanoiborra2019'
                if bRun
                    fc_kv = {'flow',fc,'fhigh',fc,'basef',fc,'erbspacebw'};
                    [~,~,out] = relanoiborra2019_featureextraction(insig, fs,fc_kv{:});                  
                    fs_an = fs;
                    an_summed = out;
                    anfH = out;
                    anfM = 0*out;
                    anfL = 0*out;
                    
                    c = [];
                    c.fs_an = fs_an;
                    c.anfH  = anfH;
                    c.anfM  = anfM;
                    c.anfL  = anfL;
                    c.out_psTH = out;
                    amt_cache('set',fname,c);
                end
                    
                text4ylabel='(MU)';
                YL = [-230 330];    
        end
    
        t_anf = (1:length(an_summed))'/fs_an;

        L_bin = round(psthbinwidth*fs_an); % samples
        L_overlap = round((percent/100)*L_bin); 
        t_anf = buffer(t_anf, L_bin, L_overlap,'nodelay');
        t_anf = t_anf(1,:);
        
        anf = buffer(an_summed, L_bin, L_overlap,'nodelay'); 
        anf = mean(anf); % one value per bin
        
        if flags.do_plot
            offx = 2.5;
            XL = [350-offx 400+offx];
            
            switch models{k}
                case {'dau1997','relanoiborra2019','osses2021'}
                    if k == 1
                        fig_functional = figure;
                        figure_handle_ax(end+1) = gca; 
                        figure_handle(end+1) = gcf; 
                        figure_name{end+1}   = ['fig11-Adaptation-' models{k}];
                        % title(sprintf('Model: %s -- Adaptation output\n(%.0f-%.0f-%.0f neurons; bin size=%.2f ms, %.1f percent overlap)',models{k},numH,numM,numL,psthbinwidth*1000,percent));
                        hold on, grid on;
                    else
                        figure(fig_functional);
                    end
                    if strcmp(models{k},'osses2021')
                        Style = '--';
                    else
                        Style = '-';
                    end
                    
                case {'zilany2014','bruce2018'}
                    figure;
                    
                    figure_handle_ax(end+1) = gca; 
                    figure_handle(end+1) = gcf; 
                    figure_name{end+1}   = ['fig11-Adaptation-' models{k}];
                    % title(sprintf('Model: %s -- Adaptation output\n(%.0f-%.0f-%.0f neurons; bin size=%.2f ms, %.1f percent overlap)',models{k},numH,numM,numL,psthbinwidth*1000,percent));
                    hold on, grid on;
                    Style = '-';
                    
                case {'verhulst2015','verhulst2018'}
                    if exist('fig_verhulst','var'), 
                      figure(fig_verhulst); 
                    else
                      fig_verhulst = figure; 
                    end
                    figure_handle_ax(end+1) = gca; 
                    figure_handle(end+1) = gcf; 
                    figure_name{end+1}   = ['fig11-Adaptation-' models{k}];
                    % title(sprintf('Model: %s -- Adaptation output\n(%.0f-%.0f-%.0f neurons; bin size=%.2f ms, %.1f percent overlap)',models{k},numH,numM,numL,psthbinwidth*1000,percent));
                    hold on, grid on;

                    Style = '-';
                    
                case 'king2019'
                    figure;
                    if k == 6
                        figure_handle_ax(end+1) = gca; 
                        figure_handle(end+1) = gcf; 
                        figure_name{end+1}   = ['fig11-Adaptation-' models{k}];
                        % title(sprintf('Model: %s -- Adaptation output\n(%.0f-%.0f-%.0f neurons; bin size=%.2f ms, %.1f percent overlap)',models{k},numH,numM,numL,psthbinwidth*1000,percent));
                        hold on, grid on;
                    end
                    Style = '-';
                    
            end
            plot(t_anf*1000,factor*anf,Style,'Color',Colours{k},'LineWidth',2); grid on, hold on
            
            switch models{k}
                case {'zilany2014','bruce2018'}
                    psth_count = out_psTH.psth;
                    t_psth = (1:length(psth_count))*out_psTH.psth_binwidth;
                    stairs(t_psth*1000,psth_count,Style,'Color',0.5*Colours{k},'LineWidth',2); grid on, hold on
                    
                    outs_anf(k).psth_count = psth_count;
                    outs_anf(k).t_psth = t_psth;
            end
            
            xlim(XL);
            ylim(factor*YL);
            xlabel('Time (ms)');
            ylabel(text4ylabel);
            
            Pos = get(gcf,'Position');
            Pos(3:4) = [360 400];
            set(gcf,'Position',Pos);
        end
        
        outs_anf(k).t_anf = t_anf;
        outs_anf(k).fs_an = fs_an;
        outs_anf(k).anf  = anf;
        outs_anf(k).anfH = anfH;
        outs_anf(k).anfM = anfM;
        outs_anf(k).anfL = anfL;
    end
    % hl = legend(text4leg,'Location','SouthEast');
    % set(hl,'FontSize',8);
    
    data.models = models;
    data.insig = insig;
    data.fs    = fs;
    data.outs_anf = outs_anf;
    data.figure_flag   = 'do_fig11';
    data.figure_handle = figure_handle;
    data.figure_name   = figure_name;    
    data.figure_handle_ax = figure_handle_ax;
end

%%% ------ FIG 12 Osses, Verhulst, and Majdak 2021 ---------------------
if flags.do_fig12 
    
    dur = 300e-3;
    rampdur = 10e-3;
    ncomponents = 10;
    dB_incr = 0;
    lvl = 50;
    nrep=100;
    numH = 12; 
    numM =  4;
    numL =  4;
    psth_binwidth=0.0005;
    
    basef = 1000; 
    baseaud = freqtoaud(basef);

    N_freqs     = 7;
    erb_step    = 3;
    freqs       =  audtofreq(baseaud-N_freqs+1:1:baseaud);
    freqs4tones = freqs(1:erb_step:end);
    
    [insig,f2test] =  il_Profile_Analysis(dur, rampdur, ncomponents, dB_incr, lvl, fs,freqs4tones);
    insig = insig(:);
    
    ti = 220; % ms
    tf = 270;
    ti_e = ti+10; % 210;
    %%%
        
    for k = 1:length(models)
        out = [];
        suff_str = [];
    
        %%% Loading flags and keyvals
        [fg,kv] = il_get_flags(models{k});
        afb_flags = fg.afb_flags;
        ihc_flags = fg.ihc_flags;
        an_flags = fg.an_flags;
        nomfb_flags = fg.nomfb_flags;
        an_kv    = kv.an_keyvals;
        %%%
    
        fname = ['fig12_profile_' models{k}];
        c = amt_cache('get',fname,flags.cachemode); 
        
        if ~isempty(c)
            bRun  = 0;
            
            out   = c.out;
            fc    = c.fc;
            fs_an = c.fs_an;
        else
            bRun = 1;
            out = [];
            amt_disp(['Calculating ' models{k} '...']);
        end
    
        fc_green = il_m2hz(401);
        
        for j = 1:N_freqs
            bin_nrs(j) = find(fc_green>freqs(j),1,'last');
        end
        % bin_nrs = bin_nrs(1:erb_step:end);
        CFs = fc_green(bin_nrs);
        
        switch models{k}
            case 'dau1997'
                if bRun
                    for j = 1:length(CFs)
                        fc_kv = {'basef',CFs(j),'flow',CFs(j),'fhigh',CFs(j),'dboffset',dBFS};
                        out_tmp = dau1997(insig,fs,afb_flags{:},ihc_flags{:}, ...
                            an_flags{:},nomfb_flags{:},fc_kv{:});
                        out(:,j) = out_tmp;
                        fc(j) = freqs(j);
                    end
                    
                    fs_an = fs;
                    
                    c = [];
                    c.out   = out;
                    c.fc    = fc;
                    c.fs_an = fs_an;
                    amt_cache('set',fname,c);
                end
                fc_ref = fc;
                ZL = [-240 1000];
                ymax = 1600;
                ymax_str = sprintf('%.0f MU',ymax/2);
                extra_str = 'a) ';
                
            case 'zilany2014'
                if bRun
                    kv={'fiberType',4,'numH',numH,'numM',numM,'numL',numL,'nrep',nrep};
                    out = zilany2014(insig,fs,fc_ref,kv{:});
                    out = out(1:length(insig),:);

                    fs_an = fs;
                    
                    c = [];
                    c.out   = out;
                    c.fc    = fc_ref;
                    c.fs_an = fs_an;
                    amt_cache('set',fname,c);
                end
                
                ymax = 1000;
                ymax_str = sprintf('%.0f spikes/s',ymax/2);
                extra_str = 'b) ';
                
            case 'verhulst2015'
                if bRun
                    out_tmp = verhulst2015(insig,fs,'abr',an_flags{:},an_kv{:},'no_cn','no_ic');

                    fc = fc_green(bin_nrs);
                    fs_an = out_tmp.fs_an;
                    
                    if strcmp(kv.an_keyvals{1},'numH')
                        numH = kv.an_keyvals{2};
                    end
                    if strcmp(kv.an_keyvals{3},'numM')
                        numM = kv.an_keyvals{4};
                    end
                    if strcmp(kv.an_keyvals{5},'numL')
                        numL = kv.an_keyvals{6};
                    end
                    out = out_tmp.an_summed(:,bin_nrs)/(numL+numM+numH);
                    
                    c = [];
                    c.out   = out;
                    c.fc    = fc;
                    c.fs_an = fs_an;
                    amt_cache('set',fname,c);
                end
                
                extra_str = 'c) ';
                ymax = 1000;
                % ymax_str = sprintf('%.0f\n spikes/s',ymax/2);
                ymax_str = sprintf('%.0f spikes/s',ymax/2);
                
            case 'verhulst2018'
                if bRun
                    out_tmp = verhulst2018(insig,fs,'abr',an_flags{:},an_kv{:},'no_cn','no_ic');

                    fc = fc_green(bin_nrs);
                    fs_an = out_tmp.fs_an;
                    if strcmp(kv.an_keyvals{1},'numH')
                        numH = kv.an_keyvals{2};
                    end
                    if strcmp(kv.an_keyvals{3},'numM')
                        numM = kv.an_keyvals{4};
                    end
                    if strcmp(kv.an_keyvals{5},'numL')
                        numL = kv.an_keyvals{6};
                    end
                    out = out_tmp.an_summed(:,bin_nrs)/(numL+numM+numH);
                    
                    c = [];
                    c.out   = out;
                    c.fc    = fc;
                    c.fs_an = fs_an;
                    amt_cache('set',fname,c);
                end
                
                extra_str = 'd) ';
                ymax = 1000;
                ymax_str = sprintf('%.0f spikes/s',ymax/2);
                
            case 'bruce2018'
                
                bUse_PSTH = 1;
                if bRun
                    fc = fc_ref;
                    kv = {'numH',numH,'numM',numM,'numL',numL,'psthbinwidth_mr',psth_binwidth,'nrep',nrep}; 
                    tmp = bruce2018(insig,fs,fc_ref,kv{:});
                    out=tmp.psth; 

                    fs_an = 1/psth_binwidth; 
                    Nf = round(dur*fs_an);
                    out = out(1:Nf,:);
                    
                    c = [];
                    c.out   = out;
                    c.fc    = fc;
                    c.fs_an = fs_an;
                    amt_cache('set',fname,c);
                end
                
                if bUse_PSTH == 0
                    % ZL = [0 1000];
                    ymax = 300;
                else
                    ymax = 1000;
                end
                ymax_str = sprintf('%.0f spikes/s',ymax/2);
                extra_str = 'e) ';
                
            case 'king2019'
                
                compression_n = 0.3; ymax = 2.5e-3; % all channels compressed
                % compression_n = 1; ymax = .125;
                if bRun
                    fc = [];
                    for j = 1:length(fc_ref)
                        fc_kv = {'basef',fc_ref(j),'flow',fc_ref(j),'fhigh',fc_ref(j),'dboffset',dBFS};

                        [out_tmp,fc(j)] = king2019(insig,fs,afb_flags{:},ihc_flags{:}, ...
                            an_flags{:},nomfb_flags{:},fc_kv{:},'compression_n',compression_n);
                        out(:,j) = out_tmp;
                    end
                    fs_an = fs;
                    
                    c = [];
                    c.out   = out;
                    c.fc    = fc;
                    c.fs_an = fs_an;
                    amt_cache('set',fname,c);
                end
                 
                extra_str = 'f) ';
                % ymax_str = sprintf('%.4f\n a.u',ymax/2);
                ymax_str = sprintf('%.4f a.u',ymax/2);
                
            case 'relanoiborra2019'
                if bRun
                    fc = [];
                    for j = 1:length(fc_ref)
                        fc_kv = {'basef',fc_ref(j),'flow',fc_ref(j),'fhigh',fc_ref(j),'erbspacebw'};
                        [~,~,out_tmp,fc(j)] = relanoiborra2019_featureextraction(insig, fs,fc_kv{:});                  
                          
                        out(:,j) = out_tmp;
                    end
                    % out = out(:,1:2*erb_step:end);
                    % fc  = fc(1:2*erb_step:end);
                    fs_an = fs;
                    
                    c = [];
                    c.out   = out;
                    c.fc    = fc;
                    c.fs_an = fs_an;
                    amt_cache('set',fname,c);
                end
                
                ZL = [-240 1000];
                ymax = 1600;
                ymax_str = sprintf('%.0f MU',ymax/2);   
                extra_str = 'g) ';
                
            case 'osses2021'
                if bRun
                    for j = 1:length(fc_ref)
                        fc_kv = {'basef',fc_ref(j),'flow',fc_ref(j),'fhigh',fc_ref(j),'dboffset',dBFS};
                        out_tmp = osses2021(insig,fs,afb_flags{:},ihc_flags{:}, ...
                            an_flags{:},nomfb_flags{:},fc_kv{:});
                        out(:,j) = out_tmp;
                        fc(j) = fc_ref(j);
                    end
                    
                    fs_an = fs;
                    
                    c = [];
                    c.out   = out;
                    c.fc    = fc;
                    c.fs_an = fs_an;
                    amt_cache('set',fname,c);
                end
                
                extra_str = 'h) ';
                ZL = [-240 1000];
                ymax = 1600;
                ymax_str = sprintf('%.0f MU',ymax/2);
        end
        
        idxi_=round(ti*1e-3*fs_an);
        idxf_=round(tf*1e-3*fs_an);
        
        offy = [];
        
        t_ms   = 1000*(1:size(insig,1))/fs;
        tan_ms = 1000*(1:size(out,1))/fs_an;
         
        bw_erb = audfiltbw(fc);
        figure;
        for i = 1:length(fc)
            % [env_here,idx_here] = il_Get_envelope(out(:,i),fs_an,fc(i)-bw_erb(i)/2); % envelope extraction assuming 'positive neurone firing'
         
            me = mean(out(idxi_:idxf_,i)); % looks for peaks above the avg in the steady section
            [env_here,idx_here] = il_Get_envelope2(out(:,i),me);
            
            offy(i) = (i-1);
         
            tan_ms_here = tan_ms(idx_here);
            
            idxs = find(tan_ms >= ti-10 & tan_ms <= tf-7.5); % 10 ms before ti
            plot(tan_ms(idxs),out(idxs,i)/ymax + offy(i),'Color',Colours{k}); hold on, grid on
            
            idxs = find(tan_ms_here >= ti-10 & tan_ms_here <= tf-7.5);
            plot(tan_ms_here(idxs),env_here(idxs)/ymax + offy(i),'Color',[0.6 0.6 0.6],'LineWidth',2);
        
            env_mean(k,i) = mean(env_here(idxs));
            env_std(k,i)  = std(env_here(idxs));
            % text(min(t_ms(idxs_pre))+2,offy+.6,sprintf('%.0f Hz',fc(i)),'Color','k');
            env_scale(k) = ymax;
            
        end
        
        %%% YAxis: scale
        tf_here=tf-7.5-.5;
        plot(tf_here*[1 1],[0 0.5],'k-','LineWidth',2);
        plot([tf_here-2 tf_here],0.5*[1 1],'k-','LineWidth',2);
        plot([tf_here-2 tf_here],    [0 0],'k-','LineWidth',2);
        text(tf_here-6,-0.35,ymax_str,'FontWeight','Bold');
        %%% XAxis: scale
        if k == 7 || k == 8
            plot([ti_e+10 ti_e+20],-0.45*[1 1],'k-','LineWidth',2);
            plot([ti_e+10 ti_e+10],[-0.45 -0.35],'k-','LineWidth',2);
            plot([ti_e+20 ti_e+20],[-0.45 -0.35],'k-','LineWidth',2);
            text(ti_e+10,-0.85,'10 ms','FontWeight','Bold');
        end
        %%%
        ylabel('Simulated CF_n (Hz)')
        xlim([ti tf]);
        ylim([-.5 7])
        
        YL = get(gca,'YLim');
        
        factor_std = 40;
        env_std_here = env_std(k,:)*(factor_std/env_scale(k));
        plot(tf+2.5-10+env_std_here, env_mean(k,:)/ymax + offy,'o--','Color',il_rgb('Maroon'),'LineWidth',2);
        plot((tf+2.5-10)*[1 1],YL,'k-');
        
        if k == 1
            x_fc   = (1:length(fc));
            
            YTL = [];
            for j = 1:length(CFs)
                YTL{j} = [num2str(round(CFs(j))) ' Hz'];
            end
        end
        set(gca,'YTick',x_fc-1);
        set(gca,'YTickLabel',YTL);
        
        set(gca,'XTick',ti+10:10:tf-10);
        set(gca,'XTickLabel',[])
        
        title([extra_str models{k}])
        
        Pos = get(gcf,'Position');
        Pos(3:4) = [380 280]; % [500 300];
        set(gcf,'Position',Pos);
        
        figure_handle(end+1) = gcf;
        figure_name{end+1} = ['fig12-profile-' models{k} suff_str];
        
        disp('')
    end
    data.figure_flag = 'do_fig12';
    
    data.env_mean = env_mean;
    data.env_std  = env_std;
    data.env_scale = env_scale;
    data.env_scale_std = env_scale/factor_std;
    
    % data.env_std(8,:) ./ (env_scale(8)/factor_std)
    
    data.insig = insig;
    data.fs    = fs;
    data.models = models;
    data.figure_handle = figure_handle;
    data.figure_name   = figure_name;  
    
end

%% ------ FIG 13 ------------------------
if flags.do_fig13ab || flags.do_fig13cd   
  
  
    error('Not implemented yet. Sorry!'); 
    

    % Adapted from g20190618_investigating_IC_CN.m
    %%% Stimulus generation:
    if flags.do_fig13ab, L=30; else L=70; end
    dBFS = 94;
    fs = 100e3;
    dur= 300e-3; % ms
    t  = 0:1/fs:dur-1/fs;
    t  = t(:); % column vector
    fc = 1000; % 4000 % Hz
    dur_samples = length(t);
    numH = 12; 
    numM =  4;
    numL =  4;
    ic_delay_inh=0.0011; % 1.1 ms as in the manuscript

    fmod_step = 5;
    fmods = 10:fmod_step:130; % 250;
    N_signals = length(fmods);
    
    insig = zeros(dur_samples,N_signals); % Memory allocation

    % Up/down cosine ramp (fixed)
    dur_ramp_ms = 5;
    dur_ramp = round((dur_ramp_ms*1e-3)*fs); % duration ramp in samples

    rp    = ones(dur_samples,1); 
    rp(1:dur_ramp) = rampup(dur_ramp);
    rp(end-dur_ramp+1:end) = rampdown(dur_ramp);

    %Stimulus
    for i = 1:length(fmods)
        fmod = fmods(i);
        mod_index=1;
        carrier=sin(2*pi*fc.*t);
        modulator=mod_index*cos(2*pi*fmod.*t+pi);
        insig_tmp=(1+modulator).*carrier;
        
        insig_tmp = scaletodbspl(insig_tmp,L,dBFS);

        insig(:,i) = rp.*insig_tmp;
    end
        
    ti = 190e-3;
    tf = 290e-3;
    for k = 1:N_models
        
        %%% Loading flags and keyvals
        [fg,kv] = il_get_flags(models{k});
        afb_flags = fg.afb_flags;
        ihc_flags = fg.ihc_flags;
        an_flags = fg.an_flags;
        an_keyvals = kv.an_keyvals;
        mfb_flags = fg.mfb_flags;
        mfb_keyvals = kv.mfb_keyvals;
        %%%
        
        switch L
            case 70
                suff_lvl='';
            otherwise
                suff_lvl = ['-' num2str(L) '-dB'];
        end
        fname = ['fig13_MTF-' models{k} suff_lvl];
        c = amt_cache('get',fname,flags.cachemode); 
  
        if ~isempty(c)
            bRun  = 0;
            
            outsig = c.outsig;
            subfs  = c.subfs;
            mfc    = c.mfc;
        else
            bRun = 1;
            outsig = [];
        end
    
        mfc_target = 100;
        switch models{k}
            case 'dau1997'
                if bRun
                    afb_kv = {'basef',fc,'flow',fc,'fhigh',fc,'dboffset',dBFS};
                    for j = 1:N_signals
                        [out_tmp,fc,mfc_here] = dau1997(insig(:,j),fs,afb_kv{:},afb_flags{:},ihc_flags{:},an_flags{:},mfb_flags{:});
                        mfc_idx = find(mfc_here<mfc_target,1,'last');
                        outsig(:,j) = out_tmp{1}(:,mfc_idx);
                        subfs = fs;
                        if j == 1
                            mfc = mfc_here(mfc_idx);
                        end
                    end
                    
                    c = [];
                    c.outsig = outsig;
                    c.subfs  = subfs;
                    c.mfc    = mfc;
                    amt_cache('set',fname,c);
                end

            case 'zilany2014'
                if bRun
                    kv={'fiberType',4,'numH',numH,'numM',numM,'numL',numL,'nrep',10};                    
                    for j = 1:length(fmods)
                        mean_rate = zilany2014(insig(:,j),fs,fc,kv{:});
                        [ic_out,~,~,par] = carney2015(mean_rate,90,fs,'ic_delay_inh',ic_delay_inh);
                        N_samples = length(insig);
                        outsig(:,j) = ic_out(1:N_samples);
                        if j == 1
                            subfs = fs;
                            mfc = 90;
                            amt_disp(sprintf('Zilany2014: CN: tau_exc=%f, tan_inh=%f, D=%f, S=%f',par.tau_ex_cn,par.tau_inh_cn,par.cn_delay,par.Sinh_cn));
                            amt_disp(sprintf('Zilany2014: IC: tau_exc=%f, tan_inh=%f, D=%f, S=%f',par.tau_ex_ic,par.tau_inh_ic,par.ic_delay_inh,par.Sinh_ic));
                        end                       
                    end

                    c = [];
                    c.outsig = outsig;
                    c.subfs  = subfs;
                    c.mfc    = mfc;
                    amt_cache('set',fname,c);
                end
                
            case 'verhulst2015'
                if bRun
                    fc_green = il_m2hz(401);
                    idx = find(fc_green>fc,1,'last');
                    
                    out = verhulst2015(insig,fs,'abr',afb_flags{:},ihc_flags{:},an_flags{:}, ...
                        kv.an_keyvals{:},mfb_flags{:},'no_ihc'); % no_ihc just reduce the data load here
                    subfs = out(1).fs_abr;
                    
                    %%%
                    if strcmp(kv.an_keyvals{1},'numH')
                        numH = kv.an_keyvals{2};
                    end
                    if strcmp(kv.an_keyvals{3},'numM')
                        numM = kv.an_keyvals{4};
                    end
                    if strcmp(kv.an_keyvals{5},'numL')
                        numL = kv.an_keyvals{6};
                    end
                    numTot = numH + numM + numL;
                    %%%
                    
                    for j = 1:length(fmods)
                        outsig(:,j) = out(j).ic(:,idx)/numTot;
                    end
                    mfc = [];
                    
                    c = [];
                    c.outsig = outsig;
                    c.subfs  = subfs;
                    c.mfc    = mfc;
                    amt_cache('set',fname,c);
                end
            
            case 'verhulst2018'
                if bRun
                    fc_green = il_m2hz(401);
                    idx = find(fc_green>fc,1,'last');
                    
                    out = verhulst2018(insig,fs,'abr',afb_flags{:},ihc_flags{:},an_flags{:}, ...
                        kv.an_keyvals{:},mfb_flags{:},'no_ihc','no_anfH','no_anfM','no_anfL');
                    subfs = out(1).fs_abr;
                    
                    %%%
                    if strcmp(kv.an_keyvals{1},'numH')
                        numH = kv.an_keyvals{2};
                    end
                    if strcmp(kv.an_keyvals{3},'numM')
                        numM = kv.an_keyvals{4};
                    end
                    if strcmp(kv.an_keyvals{5},'numL')
                        numL = kv.an_keyvals{6};
                    end
                    numTot = numH + numM + numL;
                    %%%
                    
                    for j = 1:length(fmods)
                        outsig(:,j) = out(j).ic(:,idx)/numTot;
                    end
                    mfc = [];
                    
                    c = [];
                    c.outsig = outsig;
                    c.subfs  = subfs;
                    c.mfc    = mfc;
                    amt_cache('set',fname,c);
                end
                
            case 'bruce2018'
                if bRun
                    psth_binwidth=0.0005;
                    kv = {'numH',numH,'numM',numM,'numL',numL,'nrep',10,'psthbinwidth_mr',psth_binwidth}; 
                    N_samples = size(insig,1);
                    subfs = fs;
                    mfc = 90;
                    for j = 1:length(fmods)                      
                      tmp = bruce2018(insig(:,j),fs,fc,kv{:});
                        % This is how Carney calls her model: PSTH_FT/Number_of_fibers*Model_Sampling_rate
                      ic_out = carney2015(tmp.psth_ft/(numH+numM+numL)*100000,mfc,fs,'ic_delay_inh',ic_delay_inh);
                      outsig(:,j) = ic_out(1:N_samples);
                    end
                    
                    c = [];
                    c.outsig = outsig;
                    c.subfs  = subfs;
                    c.mfc    = mfc;
                    amt_cache('set',fname,c);
                end
                
            case 'relanoiborra2019'
                if bRun
                    afb_kv = {'basef',fc,'flow',fc,'fhigh',fc,'erbspacebw'};
                    for j = 1:N_signals
                        %[out_tmp,mfc_here,~,fc(j)] = relanoiborra2019_featureextraction(insig(:,j), fs, afb_kv{:}); 
                        [out_tmp,mfc_here] = relanoiborra2019_featureextraction(insig(:,j), fs, afb_kv{:}); 
                        
                        mfc_idx = find(mfc_here<mfc_target,1,'last');
                        outsig(:,j) = squeeze(out_tmp(:,1,mfc_idx));
                        subfs = fs;
                        if j == 1
                            mfc = mfc_here(mfc_idx);
                        end
                    end
                    c = [];
                    c.outsig = outsig;
                    c.subfs  = subfs;
                    c.mfc    = mfc;
                    amt_cache('set',fname,c);
                end
                
            case 'king2019'
                if bRun
                    for j = 1:N_signals
                        afb_kv = {'basef',fc,'flow',fc,'fhigh',fc,'compression_n',.3,'dboffset',dBFS,'subfs',20000};
                        [out_tmp,fc,mfc_here,extras] = king2019(insig(:,j),fs,afb_kv{:});
                        mfc_idx = find(mfc_here<mfc_target,1,'last');

                        out_tmp = squeeze(out_tmp);
                        mfc = mfc_here(mfc_idx);
                        outsig(:,j) = out_tmp(:,mfc_idx);
                        if j == 1
                            subfs = extras.subfs;
                        end
                    end
                    c = [];
                    c.outsig = outsig;
                    c.subfs  = subfs;
                    c.mfc    = mfc;
                    amt_cache('set',fname,c);
                end
                
            case 'osses2021'
                if bRun
                    afb_kv = {'basef',fc,'flow',fc,'fhigh',fc,'dboffset',dBFS};
                    for j = 1:N_signals
                        [out_tmp,fc,mfc_here] = osses2021(insig(:,j),fs,afb_kv{:},afb_flags{:},ihc_flags{:},an_flags{:},mfb_flags{:});
                        mfc_idx = find(mfc_here<mfc_target,1,'last');
                        outsig(:,j) = out_tmp{1}(:,mfc_idx);
                        subfs = fs;
                        if j == 1
                            mfc = mfc_here(mfc_idx);
                        end
                    end
                    
                    c = [];
                    c.outsig = outsig;
                    c.subfs  = subfs;
                    c.mfc    = mfc;
                    amt_cache('set',fname,c);
                end

        end
            
        if ~isempty(outsig)
            
            switch models{k}
                case {'dau1997','relanoiborra2019','osses2021'}
                    YL = [-100 1000];
                case {'zilany2014','verhulst2015','verhulst2018','bruce2018'}
                    YL = [-200 500];
                case 'king2019'
                    YL = [-.1 1]*1e-4;
            end
                        
            idxi = round(ti*subfs)+1;
            idxf = round(tf*subfs);
            
            % Maximum:
            [vals_ma,idx_ma] = max(outsig(idxi:idxf,:));
            
            % min to max:
            Mi = min( outsig(idxi:idxf,:) );
            Ma = max( outsig(idxi:idxf,:) );
            vals_ma_mi = Ma-Mi;
            
            % this was the default:
            vals_mean_hwr  = mean( abs(outsig(idxi:idxf,:)) );
            vals_med_fwr   = median( abs(outsig(idxi:idxf,:)) );
            vals_hwr       = mean(max(outsig(idxi:idxf,:),0));
            
            vals_med  = (median( outsig(idxi:idxf,:) ));
            vals_mean = (mean( outsig(idxi:idxf,:) ));
                    
            valsL = prctile( (outsig(idxi:idxf,:)),5);
            valsU = prctile( (outsig(idxi:idxf,:)),95);
                    
            kk = 1;
            MFic(k,1:N_signals,kk)      = vals_ma; % figure; plot(fmods,vals); hold on
            MFic_norm(k,1:N_signals,kk) = vals_ma / max(abs(vals_ma));
            
            kk = 2;
            MFic(k,1:N_signals,kk)      = vals_ma_mi; % vals_mean_hwr; % figure; plot(fmods,vals); hold on
            MFic_norm(k,1:N_signals,kk) = vals_ma_mi / max(abs(vals_ma_mi));
            
            bDo = 0; % bDo = 1;
            if bDo
                close all
                figure(100); 
                for iii = 1:25
                    subplot(1,2,1)
                    plot(outsig(:,iii)); ylim(YL); grid on

                    hold on;
                    plot([idxi idxi],YL,'r--');
                    plot([idxf idxf],YL,'r--');

                    title(['model: ' models{k} '-' num2str(fmods(iii))]); 

                    hold off

                    subplot(1,2,2)
                    plot(fmods(iii),vals_ma(iii),'bo'); hold on, grid on 
                    plot(fmods(iii),vals_mean(iii),'r>');
                    plot(fmods(iii),vals_mean_hwr(iii),'ks');
                    plot(fmods(iii),vals_ma_mi(iii),'md');
                    hl = legend('max','mean','mean-hwr','ma min mi'); set(hl,'box','off');
                    xlim([0 130])
                    var_here = [vals_ma vals_mean vals_mean_hwr vals_ma_mi];
                    ylim([min(var_here) max(var_here)])
                    pause(0.5);
                end
            end
            
            switch L
                case 30
                    MF_label{1} = 'a) 30 dB, max';
                    MF_label{2} = 'b) 30 dB, max-min';
                    
                case 70
                    MF_label{1} = 'c) 70 dB, max';
                    MF_label{2} = 'd) 70 dB, max-min';
            end
        end
    end
    
    for kk = 1:length(MF_label)
        figure;
        for k = 1:N_models
            opts = {'LineStyle',LineStyle{k},'Color',Colours{k},'LineWidth',LineWidth(k)}; % ,'MarkerFaceColor','w','Markers',Markers{k},'MarkerSize',MarkersSize(k)}
            plot(fmods,MFic_norm(k,:,kk),opts{:}); grid on, hold on
        end

        % set(gca,'Fontsize',16)
        xlabel('Modulation frequency (Hz)')
        ylabel('On-CF response (Normalised)')
        xlim([min(fmods)-5 max(fmods)+5]);
        ylim([-0.05 1.15])

        Pos = get(gcf,'Position');
        Pos(3:4) = [400 360]; % [550 360];
        set(gcf,'Position',Pos);

        if kk == 1
            XTL = [];
            for i = 1:length(fmods)
                if mod(i,4) == 1
                    XTL{i} = num2str(fmods(i));
                else
                    XTL{i} = '';
                end
            end
        end
        set(gca,'XTick',fmods(1:2:end));
        set(gca,'XTickLabel',XTL(1:2:end));
        % set(gca,'XTick',fmods(1:2:end));
        set(gca,'YTick',-1:.1:1);
        figure_handle(end+1) = gcf;
        switch kk
            case 1
                suff_la = '';
            otherwise
                suff_la = ['-' num2str(kk)];
        end
        figure_name{end+1} = ['fig13-modulation-strength' suff_lvl suff_la];
        
        % text(0.05,0.92,MF_label{kk},'Units','Normalized','FontSize',14,'FontWeight','Bold');
        title(MF_label{kk});
    end
    
    data.figure_flag = 'do_fig13';
    data.fmods = fmods;
    data.MFic = MFic;
    data.MFic_norm = MFic_norm;
    data.models = models;
    data.figure_handle = figure_handle;
    data.figure_name   = figure_name; 
    
end
%% FIG 14
if flags.do_fig14
    % Local script: g20191107_rerun_Verhulst2018a.m
    
    error('Not implemented yet. Sorry!'); 
        
    % Generating click
    
    lvl       = 70; % lvl_dBnHL + 30;
    p0=2e-5;
    dur_click = 1; % 6; 
    N_samples = dur_click*fs;
    numH = 12; 
    numM =  4;
    numL =  4;
    ic_delay_inh=0.0011; % 1.1 ms as in the manuscript
    
    % 
    t=(0:1/fs:dur_click);
    click_duration = (100e-6*fs); % 100 us click (Osses and Verhulst 2019 used 80 us)
    % stim=zeros(length(t),1); %the simulation runs as long as the stimulus
     
    N_length = fs/10;
    sil_skip = N_length-click_duration; % 90.1 ms
     
    samples_click      = [];
    samples_click_even = [];
    
    Nr_clicks = floor(fs*dur_click/(sil_skip+click_duration));
     
    idx_offset = round(10e-3*fs); % click starts 10 ms after time 0
    for i = 1:Nr_clicks
        start_sample = (i-1)*N_length+idx_offset;
        if mod(i,2) == 1
            samples_click      = [samples_click start_sample+click_duration:start_sample+click_duration+click_duration];
        else
            samples_click_even = [samples_click_even start_sample+click_duration:start_sample+click_duration+click_duration];
        end 
    end
     
    A = 2*sqrt(2)*p0*10^(lvl/20);
    insig = zeros(N_samples,1);
    insig(samples_click)      =  A;
    insig(samples_click_even) = -A;
    
    for k = 1:N_models
        %%% Loading flags and keyvals
        [fg,kv] = il_get_flags(models{k});
        afb_flags = fg.afb_flags;
        ihc_flags = fg.ihc_flags;
        an_flags = fg.an_flags;
        an_keyvals = kv.an_keyvals;
        mfb_flags = fg.mfb_flags;
        nomfb_flags = fg.nomfb_flags;
        mfb_keyvals = kv.mfb_keyvals;
        %%%
    
        fname = ['fig14_click-' models{k}];
        c = amt_cache('get',fname,flags.cachemode); 
  
        if ~isempty(c)
            bRun  = 0;
            
            out   = c.out;
            subfs = c.subfs;
            mfc   = c.mfc;
        else
            bRun = 1;
            out = [];
        end
    
        mfc_target = 100;
        switch models{k}
            case 'dau1997'
                if bRun
                    subfs = 20000;
                    afb_kv = {'subfs',subfs,'dboffset',dBFS};
                 
                    t_start = tic; 
                    [out_tmp,fc_here,mfc_here] = dau1997(insig,fs,afb_kv{:},afb_flags{:},ihc_flags{:},an_flags{:},mfb_flags{:});
                    out.t_elapsed = toc(t_start);
                 
                    % % Adaptation loops: no time count
                    % %    dau1997 cannot return (at this point) two output signals
                    % out_tmp_an = dau1997(insig,fs,afb_kv{:},afb_flags{:},ihc_flags{:},an_flags{:},nomfb_flags{:});
                    % out_tmp_an = fftresample(out_tmp_an,round(length(out_tmp_an)/fs*subfs)); % Resampling (as in dau1997):
                    
                    mfc_idx = find(mfc_here<mfc_target,1,'last');
                    
                    Nr_fc = length(fc_here);
                    outsig = []; 
                    % outsig_adapt = [];
                    fcs = [];
                    for j = 1:Nr_fc
                        if size(out_tmp{j},2) >= mfc_idx
                            fcs(end+1) = fc_here(j); % only adding the fcs of filters containing mfc_idx
                            outsig(:,end+1) = out_tmp{j}(:,mfc_idx);
                            
                            % outsig_adapt(:,end+1) = out_tmp_an(:,j);
                        end
                    end
                    
                    out.fcs    = fcs;
                    out.oustig = outsig;
                    out.outsig_description = 'Modulation filter bank output';
                    
                    % out.outsig_adapt = outsig_adapt;
                    % out.outsig_adapt_description = 'Adaptation loops output';
                
                    mfc = mfc_here(mfc_idx);
                 
                    c = [];
                    c.out   = out;
                    c.subfs = subfs;
                    c.mfc   = mfc;
                    amt_cache('set',fname,c);
                end
    
                factor = 1;
                unit = 'MU';
                YL  = [-20 560];
                YT = 0:50:500;
                
            case 'zilany2014'
                if bRun
                    fcs = il_m2hz(401);
                    fcs = fcs(end:-6:1);
                    fcs = fcs(fcs>125);
%                              fcs=fcs(1:10:end); % for quick debugging only
                    psth_binwidth=0.0005;
                    kv={'fiberType',4,'numH',numH,'numM',numM,'numL',numL,'nrep',10,'psth_binwidth',psth_binwidth}; 
                    subfs = fs;
                    
                    t_start = tic;
                    [mean_rate,psth] = zilany2014(insig,fs,fcs,kv{:});
                    [ic_out,~,cn_out,par] = carney2015(mean_rate,90,fs,'ic_delay_inh',ic_delay_inh);
                    t_elapsed = toc(t_start);
                    out.an(1:N_samples,:) = mean_rate(1:N_samples,:);
                    out.cn(1:N_samples,:) = cn_out(1:N_samples,:);
                    out.ic(1:N_samples,:) = ic_out(1:N_samples,:);
                    out.psth = psth;
                    out.t_elapsed = t_elapsed*ones(size(fcs))/length(fcs);
                    out.psth_binwidth = psth_binwidth;
                    out.fcs = fcs;
                    out.fcs_description = 'CFs that were tested';
                    mfc = [];
                    
                    amt_disp(sprintf('Zilany2014: CN: tau_exc=%f, tan_inh=%f, D=%f, S=%f',par.tau_ex_cn,par.tau_inh_cn,par.cn_delay,par.Sinh_cn));
                    amt_disp(sprintf('Zilany2014: IC: tau_exc=%f, tan_inh=%f, D=%f, S=%f',par.tau_ex_ic,par.tau_inh_ic,par.ic_delay_inh,par.Sinh_ic));                            
                    
                    out.an_all = sum(out.an,2);
                    out.cn_all = sum(out.cn,2);
                    out.ic_all = sum(out.ic,2);
                    out = rmfield(out,{'an','cn','ic'}); % to reduce size
                    
                    c = [];
                    c.out   = out;
                    c.subfs = subfs;
                    c.mfc   = mfc;
                    amt_cache('set',fname,c);
                end
                
                factor = 1;
                unit = 'spikes/s';
                YL  = [-280 280];
                YT = -240:40:240;
                
            case 'verhulst2015'
                if bRun
                    t_start = tic; 
                    out = verhulst2015(insig,fs,'abr',afb_flags{:},ihc_flags{:},an_flags{:}, ...
                        kv.an_keyvals{:},mfb_flags{:});
                    t_elapsed = toc(t_start);
                    
                    subfs = out(1).fs_abr;
                    out.t_elapsed = t_elapsed;
                    out = rmfield(out,{'v','fs_bm','ihc','anfH','anfM','anfL','an_summed','ic'});
                    mfc = [];
                    
                    c = [];
                    c.out   = out;
                    c.subfs = subfs;
                    c.mfc   = mfc;
                    amt_cache('set',fname,c);  
                end
                
                factor = 1e6; unit = '\muV';
                YL  = [-.32 .32];
                YT = -.3:.05:.3;
                                
            case 'verhulst2018'
                if bRun
                    t_start = tic; 
                    out = verhulst2018(insig,fs,'abr',afb_flags{:},ihc_flags{:},an_flags{:}, ...
                        kv.an_keyvals{:},mfb_flags{:});
                    t_elapsed = toc(t_start);
                    
                    subfs = out(1).fs_abr;
                    out.t_elapsed = t_elapsed;
                    out = rmfield(out,{'v','fs_bm','ihc','anfH','anfM','anfL','an_summed','ic'});
                    mfc = [];
                    
                    c = [];
                    c.out   = out;
                    c.subfs = subfs;
                    c.mfc   = mfc;
                    amt_cache('set',fname,c);
                end
                
                factor = 1e6; unit = '\muV';
                YL  = [-.32 .32];
                YT = -.3:.05:.3;
    
            case 'bruce2018'
                if bRun
                    fcs = il_m2hz(401);
                    fcs = fcs(end:-6:1);
                    fcs = fcs(fcs>125);
                    psth_binwidth=0.0005;
                    kv = {'numH',numH,'numM',numM,'numL',numL,'nrep',10,'psthbinwidth_mr',psth_binwidth}; 
%                            fcs=fcs(1:10:end); % for quick debugging only
                    tmp = bruce2018(insig,fs,fcs,kv{:});
                      % This is how Carney calls her model: PSTH_FT/Number_of_fibers*Model_Sampling_rate
                      %   [ic_out,cn_out] = carney2015(tmp.psth_ft/(numH+numM+numL)*100000,90,fs,'ic_delay_inh',ic_delay_inh);
                      % But we also can use PSTH directly: 
                    [ic_out,cn_out] = carney2015(tmp.psth,90,1/psth_binwidth,'ic_delay_inh',ic_delay_inh);

                    t_elapsed = toc(t_start);
                    out.t_elapsed = t_elapsed*ones(size(fcs))/length(fcs);
                    subfs = 1/psth_binwidth;
                    out.cn_all = sum(cn_out(1:round(N_samples/fs/psth_binwidth),:),2);
                    out.ic_all = sum(ic_out(1:round(N_samples/fs/psth_binwidth),:),2);
                    out.fcs = fcs;
                    out.fcs_description = 'CFs that were tested';
                    mfc = [];
                    
                    c = [];
                    c.out   = out;
                    c.subfs = subfs;
                    c.mfc   = mfc;
                    amt_cache('set',fname,c);
                end
                
                factor = 1;
                unit = 'spikes/s';
                YL  = [-280 280];
                YT = -240:40:240;
            
            case 'king2019'
                
                if bRun
                    subfs = 20000;
                    afb_kv = {'basef',4000,'flow',80,'fhigh',8000, ... % arbitrary
                        'compression_n',.3,'dboffset',dBFS,'subfs',subfs, ...
                        'modbank_Nmod',10}; % 10 modulation filters instead of 5
                     
                    t_start = tic; 
                    [outsig,fcs,mfc_here,step_here] = king2019(insig,fs,afb_kv{:},afb_flags{:},ihc_flags{:},an_flags{:},mfb_flags{:});
                    out.t_elapsed = toc(t_start);
                 
                    mfc_idx = find(mfc_here<mfc_target,1,'last');
                    outsig = squeeze(outsig(:,:,mfc_idx));
                                        
                    out.fcs    = fcs;
                    out.oustig = outsig;
                    out.outsig_description = 'Modulation filter bank output';
                    
                    % outsig_adapt = step_here.adapted_response;
                    % outsig_adapt = fftresample(outsig_adapt,round(length(outsig_adapt)/fs*subfs));
                    % 
                    % out.outsig_adapt = outsig_adapt;
                    % out.outsig_adapt_description = 'Adaptation by filtering';
                
                    mfc = mfc_here(mfc_idx);
                 
                    c = [];
                    c.out   = out;
                    c.subfs = subfs;
                    c.mfc   = mfc;
                    amt_cache('set',fname,c);
                end
    
                factor = 1e4;
                unit = 'a.u. x 10^{-4}';
                YL  = [-0.1 2.1];
                YT = [0:.2:2];
                
            case 'relanoiborra2019'
                
                if bRun
                    afb_kv = {'erbspacebw'};
                    t_start = tic; 
%                     [out_tmp,fc_here,mfc_here] = relanoiborra2019_preproc(insig,fs,afb_kv{:},afb_flags{:},ihc_flags{:},an_flags{:},mfb_flags{:});
                    [out_tmp,mfc_here,~,fc_here] = relanoiborra2019_featureextraction(insig, fs, afb_kv{:});                  
                    out.t_elapsed = toc(t_start);
                    % out_tmp = [length fc mfc], i.e., [100000 60 12]
                    
                    mfc_idx = find(mfc_here<mfc_target,1,'last');
                    
                    Nr_fc = length(fc_here);
                    outsig = []; % zeros(size(out_tmp{1}(:,1)));
                    % outsig_adapt = [];
                    fcs = [];
                    for j = 1:Nr_fc
                        if size(out_tmp,3) >= mfc_idx
                            fcs(end+1) = fc_here(j); % only adding the fcs of filters containing mfc_idx
                            outsig(:,end+1) = out_tmp(:,j,mfc_idx);
                            
                            % outsig_adapt(:,end+1) = out_tmp_an(:,j);
                        end
                    end
                    
                    out.fcs    = fcs;
                    out.oustig = outsig;
                    out.outsig_description = 'Modulation filter bank output';
                    
                    % out.outsig_adapt = outsig_adapt;
                    % out.outsig_adapt_description = 'Adaptation loops output';
                
                    mfc = mfc_here(mfc_idx);
                 
                    c = [];
                    c.out   = out;
                    c.subfs = fs;
                    c.mfc   = mfc;
                    amt_cache('set',fname,c);
                end
    
                factor = 1;
                unit = 'MU';
                YL  = [-20 560];
                YT = 0:50:500;
                
            case 'osses2021'
                if bRun
                    subfs = 20000;
                    afb_kv = {'subfs',subfs,'dboffset',dBFS};
                 
                    t_start = tic; 
                    [out_tmp,fc_here,mfc_here] = osses2021(insig,fs,afb_kv{:},afb_flags{:},ihc_flags{:},an_flags{:},mfb_flags{:});
                    out.t_elapsed = toc(t_start);
                 
                    mfc_idx = find(mfc_here<mfc_target,1,'last');
                    
                    Nr_fc = length(fc_here);
                    outsig = []; 
                    
                    fcs = [];
                    for j = 1:Nr_fc
                        if size(out_tmp{j},2) >= mfc_idx
                            fcs(end+1) = fc_here(j); % only adding the fcs of filters containing mfc_idx
                            outsig(:,end+1) = out_tmp{j}(:,mfc_idx);
                        end
                    end
                    
                    out.fcs    = fcs;
                    out.oustig = outsig;
                    out.outsig_description = 'Modulation filter bank output';
                    
                    mfc = mfc_here(mfc_idx);
                 
                    c = [];
                    c.out   = out;
                    c.subfs = subfs;
                    c.mfc   = mfc;
                    amt_cache('set',fname,c);
                end
    
                factor = 1;
                unit = 'MU';
                YL  = [-20 560];
                YT = 0:50:500;
        end
    
        if ~isempty(out)
            
            w5 = [];
            switch models{k}
                case {'dau1997','king2019','relanoiborra2019','osses2021'}
                    num_cfs = length(out.fcs);
                    w5 = sum(out.oustig,2)/num_cfs;
                    % w3 = zeros(size(w5)); % no comparable to w3
                    % w1 = sum(out.outsig_adapt,2)/num_cfs;
                    
                case {'zilany2014','bruce2018'}
                    num_cfs = length(out.fcs);
                    % w1 = out.an_all/num_cfs; w3 = out.cn_all/num_cfs;
                    w5 = out.ic_all/num_cfs;
                case {'verhulst2015','verhulst2018'}
                    % w1 = out.w1; w3 = out.w3;
                    w5 = out.w5;
            end
            
            % -------------------------------------------------------------
            % Time elapsed:
            if length(out.t_elapsed) ~= 1
                t_elapsed_here = sum(out.t_elapsed);
            else
                t_elapsed_here = out.t_elapsed;
            end
            t_elapsed(1,k) = t_elapsed_here;
            t_elapsed(2,k) = num_cfs;
            t_elapsed(3,k) = t_elapsed(1,k)/num_cfs;
            if k == 1
                t_elapsed_description = 'Row 1: Total time elapsed (s); Row 2: Total number of simulated channels; Row 3 = time per channel (Row 1/Row 2)';
            end
            fprintf('Model: %s; time elapsed=%.3f s (%.4f s per channel, %.0f channels)\n',models{k},t_elapsed(1,k),t_elapsed(3,k),t_elapsed(2,k));
            % END Time elapsed --------------------------------------------
           
            t_ms = 1000*(1:length(w5))/subfs;
            
            disp('Using Wave-V estimate only...')
            efr_here = w5;
            
            figure; % For each model there is one new figure
            figure_handle(end+1) = gcf;
            figure_name{end+1} = ['fig14-click-' models{k}];
                        
            opts = {'LineStyle','-','Color',Colours{k},'LineWidth',LineWidth(k)};
            
            % if k == 1
            %     warning('bReshape is temporal...')
            % end
            
            t_ms_orig = t_ms;
               
            %%% Getting the clicks to be plotted only
            idxs = [1 2 9 10]; % Click numbers to be plotted
            
            N = length(efr_here)/Nr_clicks;
            M = Nr_clicks;

            efr_here = reshape(efr_here,N,M);
            t_ms_orig = t_ms;
            t_ms     = reshape(t_ms    ,N,M);

            efr_here = efr_here(:,idxs);
            Nr_clicks_here = length(idxs);

            t_ms = t_ms(:,idxs);

            efr_here = efr_here(:); % truncated clicks back to a one-column array
            %%%
            
            pl(k) = plot(t_ms_orig(1:length(efr_here)),factor*efr_here,opts{:}); hold on, grid on
            title(models{k})
            
            xlabel('Time (ms)');

            step_time = 40; % 30
            idxi = 210;
            idxf = idxi + 180;
            XT = [10:step_time:170 idxi:step_time:idxf];

            idxi = floor(t_ms(1,3));
            idxf = idxi + 180;
            XT_for_tick = [10:step_time:170 idxi:step_time:idxf];
            XTL = [];
            for ii = 1:length(XT)
                if mod(ii,2) == 1
                    XTL{ii} = num2str(XT_for_tick(ii));
                else
                    XTL{ii} = '';
                end
            end
            set(gca,'XTick'     ,XT);
            set(gca,'XTickLabel',XTL);

            xlim([-5 390]);
            
            ylim(YL);
            
            YTL = [];
            for ii = 1:length(YT)
                if mod(ii,2) == 1
                    YTL{ii} = num2str(YT(ii));
                else
                    YTL{ii} = '';
                end
            end
            set(gca,'YTick',YT);
            set(gca,'YTickLabel',YTL);
            
            % ylabel(['Amplitude (' unit ')']);
            ylabel(['(' unit ')']);
            
            Pos = get(gcf,'Position');
            Pos(3:4) = [380 300];
            set(gcf,'Position',Pos);
            
            
            YL = get(gca,'YLim');
            plot(200*[1 1],YL,'k-');

            text(0.05,.92,sprintf('Clicks #%.0f-%.0f',idxs(1:2)),'Units','Normalized','FontSize',14,'FontWeight','Bold');
            text(0.55,.92,sprintf('Clicks #%.0f-%.0f',idxs(3:4)),'Units','Normalized','FontSize',14,'FontWeight','Bold');
        end % end ~ifempty(out)
        
        Nr_samples = length(efr_here); 
        output_buf = reshape(efr_here,round(Nr_samples/Nr_clicks_here),Nr_clicks_here);
        
        Amp_max(k,:) = max(output_buf);
        Amp_min(k,:) = min(output_buf);
    end % end for k
    
    data.Amp_max = Amp_max;
    data.Amp_min = Amp_min;
    data.Amp_pp  = Amp_max - Amp_min;
    data.models  = models;
    
    data.figure_flag = 'do_fig14';
    data.figure_handle = figure_handle;
    data.figure_name   = figure_name;
end


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Start of subfunctions:
function [flags,keyvals] = il_get_flags(model_str)

kSR_H  = 68.5; % spikes/s, spontaneous rate: only Verhulst 2018
kSR_M  = 10; % spikes/s, spontaneous rate: Verhulst 2015 and 2018
kSR_L  = 1; % spikes/s, spontaneous rate: Verhulst 2015 and 2018

numH = 12;
numM =  4;
numL =  4;

bAMT = 1; % Is AMT v1.0

switch model_str
    case 'dau1997'
        fc_flags   = [];
        afb_flags = {};
        afb_keyvals = [];
        ihc_flags  = {'ihc','ihc_dau1996'}; 
        noihc_flags  = {'no_ihc'}; 
        an_flags   = {'adt','adt_dau1997'};
        noan_flags = {'no_adt','no_mfb'};
        an_keyvals = [];
        mfb_flags  = {'mfb','mfb_dau1997'};
        nomfb_flags  = {'no_mfb'};
        mfb_keyvals = [];
        
    case 'osses2021'
        fc_flags   = [];
        afb_flags  = {'afb_osses2021'};
        afb_keyvals = [];
        ihc_flags  = {'ihc','ihc_breebaart2001'}; 
        noihc_flags  = {'no_ihc'}; 
        an_flags   = {'adt', 'adt_osses2021'};
        noan_flags = {'no_adt','no_mfb'};
        an_keyvals = [];
        mfb_flags  = {'mfb','mfb_jepsen2008'}; 
        nomfb_flags  = {'no_mfb'};
        mfb_keyvals = [];
        
    case 'relanoiborra2019'
        fc_flags   = [];
        afb_flags  = {'no_internalnoise','erbspacebw'};
        afb_keyvals = [];
        ihc_flags  = {'ihc','ihc_relanoiborra2019'}; 
        noihc_flags  = {'no_ihc'}; 
        an_flags   = {'adt', 'adt_relanoiborra2019'};
        noan_flags = {'no_adt','no_mfb'};
        an_keyvals = [];
        mfb_flags  = {'mfb','mfb_jepsen2008'};
        nomfb_flags  = {'no_mfb'};
        mfb_keyvals = [];
        
    case 'king2019'
        fc_flags   = [];
        afb_flags  = {'afb','compression_brokenstick'};
        afb_keyvals = {'compression_n',0.3};
        ihc_flags  = {'ihc'}; 
        noihc_flags  = {'no_ihc'}; 
        an_flags   = {'adt'};
        noan_flags = {'no_adt','no_mfb'};
        an_keyvals = [];
        mfb_flags  = {'mfb'};
        nomfb_flags  = {'no_mfb'};
        mfb_keyvals = [];
        
    case 'verhulst2015'
        fc_flags   = {'abr'}; % abr = 401 frequency channels
        afb_flags  = {'no_outerear','middleear'}; % these are the defaults
        afb_keyvals = [];
        ihc_flags  = {'ihc'}; % this is the default 
        noihc_flags  = {'no_ihc'};
        an_flags   = {'an'};
        noan_flags = {'no_an','no_cn','no_ic'}; % if no_an, then no_cn and no_ic
        an_keyvals = {'numH',numH,'numM',numM,'numL',numL,'kSR_H',kSR_H,'kSR_M',kSR_M,'kSR_L',kSR_L};
        mfb_flags  = {'no_cn','ic'};
        nomfb_flags  = {'no_cn','no_ic'};
        mfb_keyvals = [];
        
    case 'verhulst2018'
        fc_flags   = {'abr'}; % abr = 401 frequency channels
        afb_flags  = {'no_outerear','middleear'}; % these are the defaults
        afb_keyvals = [];
        ihc_flags  = {'ihc'}; % this is the default 
        noihc_flags  = {'no_ihc'};
        an_flags   = {'an'};
        noan_flags = {'no_an','no_cn','no_ic'};
        an_keyvals = {'numH',numH,'numM',numM,'numL',numL,'kSR_H',kSR_H,'kSR_M',kSR_M,'kSR_L',kSR_L};
        mfb_flags  = {'no_cn','ic'};
        nomfb_flags  = {'no_cn','no_ic'};
        mfb_keyvals = [];
        
    case 'zilany2014'
        species   = 'human'; 
        noiseType = 'fixedFGn';
        % noiseType = 'varFGn';
        fc_flags   = [];
        if bAMT
            afb_flags = {species};
        else 
            afb_flags = {'middleear',species}; % this is the default, making sure that 'species' is always loaded
        end
        afb_keyvals = [];
        ihc_flags  = {'ihc'}; % this is the default 
        noihc_flags  = {'no_ihc'};
        an_flags   = {'an',species,noiseType};
        if numH ~= 0
            an_flags(end+1) = {'anfH'};
        end
        if numM ~= 0
            an_flags(end+1) = {'anfM'};
        end
        if numL ~= 0
            an_flags(end+1) = {'anfL'};
        end
        noan_flags = {'no_an',species,noiseType};
        an_keyvals = {'numH',numH,'numM',numM,'numL',numL,'kSR_H',kSR_H,'kSR_M',kSR_M,'kSR_L',kSR_L, ...
            'psth_binwidth',0.5e-3};
        mfb_flags  = {'cn','ic'};
        nomfb_flags  = {'no_cn','no_ic'};
        mfb_keyvals = {'BMF',90,'apply_ic_hwr',0};
        
    case 'bruce2018'
        species   = 'human'; 
        noiseType = 'fixedFGn';
        % noiseType = 'varFGn';

        fc_flags   = [];
        if bAMT
            afb_flags = {species};
        else 
            afb_flags = {'middleear',species}; % this is the default, making sure that 'species' is always loaded
        end
        afb_keyvals = [];
        ihc_flags  = {'ihc'}; % this is the default 
        noihc_flags  = {'no_ihc'};
        an_flags   = {'an',species,noiseType};
        if numH ~= 0
            an_flags(end+1) = {'anfH'};
        end
        if numM ~= 0
            an_flags(end+1) = {'anfM'};
        end
        if numL ~= 0
            an_flags(end+1) = {'anfL'};
        end
        noan_flags = {'no_an',species,noiseType};
        an_keyvals = {'numH',numH,'numM',numM,'numL',numL,'kSR_H',kSR_H,'kSR_M',kSR_M,'kSR_L',kSR_L, ...
            'psth_binwidth',0.5e-3};
        mfb_flags  = {'cn','ic'};
        nomfb_flags  = {'no_cn','no_ic'};
        mfb_keyvals = {'BMF',90,'apply_ic_hwr',0};
end

flags = [];
flags.fc_flags   = fc_flags;
flags.afb_flags  = afb_flags;
flags.ihc_flags  = ihc_flags;
flags.noihc_flags  = noihc_flags;
flags.an_flags   = an_flags;
flags.noan_flags = noan_flags;
flags.mfb_flags  = mfb_flags;
flags.nomfb_flags= nomfb_flags;

keyvals = [];
keyvals.afb_keyvals = afb_keyvals;
keyvals.an_keyvals  = an_keyvals;
keyvals.mfb_keyvals = mfb_keyvals;

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function rms_val = il_rmsdb(insig)

rms_val = 20*log10(rms(insig));

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function [ac_target,dc_target,vrest] = il_get_ac_dc_osses2022(ihc_target,ac_reg,dc_reg)
% function [ac_target,dc_target,vrest] = il_get_ac_dc_osses2022(ihc_target,ac_reg,dc_reg)
%
% DC: is computed from the signal (mean value using the ac_reg samples), and 
%     then this value is referenced to the resting potential
% AC: I compute it here as the peak-to-peak voltage, in line with Russel and
%     Palmer 1986

vrest     = median(ihc_target(dc_reg,:)); % median of the silent initial segment            
dc_target = mean(ihc_target(ac_reg,:))-vrest;
        
v_min   = min(ihc_target(ac_reg,:));
v_max   = max(ihc_target(ac_reg,:));

dist_ac_pp = v_max-v_min; % peak to peak
ac_target  = dist_ac_pp; % approximation to the RMS of the AC  

if nargout == 0
    t_ds = 1:size(ihc_target,1);
    
    for j = 1:size(ihc_target,2);
        figure;
        plot(t_ds,ihc_target(:,j)); hold on;

        plot(t_ds(dc_reg),ihc_target(dc_reg,j),'r');
        dc2plot = [vrest(j) vrest(j)+dc_target(j)];
        plot(t_ds(dc_reg(end))*[1 1],dc2plot,'r');
        plot([min(t_ds) max(t_ds)],vrest(j)*[1 1],'r--','LineWidth',2);

        plot(t_ds(ac_reg),ihc_target(ac_reg,j),'m-');
        ac2plot = (v_min(j)+dist_ac_pp(j)/2)*[1 1];
        plot(t_ds([ac_reg(1) ac_reg(end)]),ac2plot,'k--','LineWidth',2);
        plot(t_ds([ac_reg(1000) ac_reg(1000)]),dc2plot,'k-');
        
        xlabel('Time (samples)');
        ylabel('IHC receptor potential (V)');
        title(sprintf('Signal nr. %.0f Hz \n(AC=%.4f; DC=%.4f; ratio=%.4f)',j, ...
            dc_target(j),ac_target(j),ac_target(j)/dc_target(j)));
        legend('IHC response','DC component','AC component');
    end
end

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function [f,x] = il_m2hz(Nr_sections,k)
% function [f,x] = il_m2hz(Nr_sections,k)
%
% 1. Description:
%       See Greenwood 1990, Equation 1 with A = 165.4, k = 0.85.
%       bApex_to_base is referenced to Helicotrema.
%
% Programmed by Alejandro Osses, WAVES, UGent, Belgium, 2018-2020

if nargin < 2
    k = 0.85; % As used in the models by Verhulst et al
end

if nargin == 0
    Nr_sections = 500;
end

cochleaLength = 0.035; % m
helicotremaLength = 0.001; % m
A = 165.4118; % For human specie
a = 61.765; % 
bm_length = cochleaLength-helicotremaLength;

switch Nr_sections
    case {201,401}
        % Nr_sections = 500;
        f_range = [12000 112];
        step_bm = bm_length/500;
    otherwise
        f_range = [];
        step_bm = bm_length/Nr_sections;
end

A0 = 20682; % cochlear_model2018.py, L114
switch Nr_sections
    case 201
        x = step_bm  :2*step_bm:bm_length;

    case 401
        x = step_bm  :step_bm:bm_length;

    case 500
        x = step_bm/2:step_bm:bm_length;

    case 1000
        x = step_bm  :step_bm:bm_length;

    otherwise
        error('Check this function: %s',mfilename);
end
    
f = A0*(10.^(-a*x)) - A*k;

if ~isempty(f_range)
    idxi = find(f>= f_range(1),1,'last'); 
    idxf = find(f>= f_range(2),1,'last'); 
            
    f = f(idxi:idxf);
    x = x(idxi:idxf);
end

if nargout == 0
    fprintf('The output CF derived from %.0f cochlear sections with frequencies between %.1f and %.1f [Hz]\n',length(f),f(1),f(end));
end

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function col_rgb = il_rgb(colour)

switch colour
    case 'Gray'
        col_rgb = [0.5 0.5 0.5];
    case 'LightGray'
        col_rgb = [0.8242 0.8242 0.8242];
    case 'Green'
        col_rgb = [0 0.5 0];
    case 'LightSkyBlue'
        col_rgb = [0.5273 0.8047 0.9792];
    case 'Maroon'
        col_rgb = [0.5 0 0];
    case 'mediumorchid'
        col_rgb = [0.7266 0.3320 0.8242];
end

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function [env,idx_env] = il_Get_envelope2(insig,me)

method = 2;

switch method
    case 1
        error('Not used, tested on 7/4/2021')
    case 2
        
        % bDebug = 0;
        % if bDebug
        %     plot(insig); hold on
        % end
        env = [];
        idx_env = [];
        if nargin < 4
            me = mean(insig); % prctile(insig,75);
        end
        
        L1 = [1:length(insig); insig'];
        L2 = [1 length(insig); me me];
        
        P = il_InterX(L1,L2);
        N_P = size(P,2);
        
        % if bDebug
        %     YL = get(gca,'YLim');
        %     for j = 1:N_P
        %         plot(round(P(1,j))*[1 1],YL,'k-');
        %     end
        % end
        
        for j = 1:2:N_P-2
            idx1 = round(P(1,j));
            idx2 = round(P(1,j+1));
            [env(end+1),idx_here] = max(insig(idx1:idx2));
            
            idx_env(end+1) = idx_here+idx1-1;
        end
end

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function [BW_ERB, Q_ERB, extra]= il_Get_ERB_estimation(impulse_ch, CF, fs)

N = length(impulse_ch); % N-point FFT
K = floor(N/2);         % K-non redundant points

freq_here = (1:K)*(fs/2)/K;

Amp     = (2*abs(fft(impulse_ch,N))/N).^2; Amp = Amp(:);
[M,Mn]  = max(Amp(1:K));
E       = sum(Amp(1:K));
BW_ERB  = (E/M)*fs/N;
Q_ERB   = CF/BW_ERB;

Amp_dB = 10*log10(Amp(1:K));
Amp_dB = Amp_dB-max(Amp_dB);

L1 = [freq_here; transpose(Amp_dB)];
L2 = [min(freq_here) max(freq_here); [-10 -10]];
P = il_InterX(L1, L2);
if ~isempty(P) && size(P,1)==2
    flow_10 = P(1,1);
    fhigh_10 = P(1,2);
    BW10 = P(1,2)-P(1,1);
    Q10  = CF/BW10;
else
    flow_10 = nan;
    fhigh_10 = nan;
    BW10 = nan;
    Q10  = nan;
end

L2 = [min(freq_here) max(freq_here); [-3 -3]];
P = il_InterX(L1, L2);
if ~isempty(P) && size(P,1)==2
    flow_03 = P(1,1);
    fhigh_03 = P(1,2);
    BW03 = P(1,2)-P(1,1);
    Q03  = CF/BW03;
else
    flow_03 = nan;
    fhigh_03 = nan;
    BW03 = nan;
    Q03  = nan;
end
extra.BW_ERB = BW_ERB;
extra.BW10 = BW10;
extra.BW03 = BW03;
extra.Q_ERB = Q_ERB;
extra.Q10 = Q10;
extra.Q03 = Q03;
extra.flow_03 = flow_03;
extra.flow_10 = flow_10;
extra.fhigh_03 = fhigh_03;
extra.fhigh_10 = fhigh_10;

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function P = il_InterX(L1,L2)
%IL_INTERX Intersection of curves
%   P = IL_INTERX(L1,L2) returns the intersection points of two curves L1 
%   and L2. The curves L1,L2 can be either closed or open and are described
%   by two-row-matrices, where each row contains its x- and y- coordinates.
%
%   P has the same structure as L1 and L2, and its rows correspond to the
%   x- and y- coordinates of the intersection points of L1 and L2. If no
%   intersections are found, the returned P is empty.
%
%   Two words about the algorithm: Most of the code is self-explanatory.
%   The only trick lies in the calculation of C1 and C2. To be brief, this
%   is essentially the two-dimensional analog of the condition that needs
%   to be satisfied by a function F(x) that has a zero in the interval
%   [a,b], namely
%           F(a)*F(b) <= 0
%   C1 and C2 exactly do this for each segment of curves 1 and 2
%   respectively. If this condition is satisfied simultaneously for two
%   segments then we know that they will cross at some point. 
%   Each factor of the 'C' arrays is essentially a matrix containing 
%   the numerators of the signed distances between points of one curve
%   and line segments of the other.
%
% This function is based on InterX.m (v3.0, 21 Sept. 2010) by author 'NS'
% that is available within the Mathworks FileExchange.

% Reordering the data before intersecting L1 and L2:
x1  = L1(1,:)';  x2 = L2(1,:);
y1  = L1(2,:)';  y2 = L2(2,:);
dx1 = diff(x1); dy1 = diff(y1);
dx2 = diff(x2); dy2 = diff(y2);

%...Determine 'signed distances'   
S1 = dx1.*y1(1:end-1) - dy1.*x1(1:end-1);
S2 = dx2.*y2(1:end-1) - dy2.*x2(1:end-1);

C1 = le(il_D(bsxfun(@times,dx1,y2)-bsxfun(@times,dy1,x2),S1),0);
C2 = le(il_D((bsxfun(@times,y1,dx2)-bsxfun(@times,x1,dy2))',S2'),0)';

%...Obtain the segments where an intersection is expected
[i,j] = find(C1 & C2); 
if isempty(i),P = zeros(2,0);return; end;

%...Transpose and prepare for output
i=i'; dx2=dx2'; dy2=dy2'; S2 = S2';
L = dy2(j).*dx1(i) - dy1(i).*dx2(j);
i = i(L~=0); j=j(L~=0); L=L(L~=0);  %...Avoid divisions by 0

%...Solve system of eqs to get the common points
P = unique([dx2(j).*S1(i) - dx1(i).*S2(j), ...
            dy2(j).*S1(i) - dy1(i).*S2(j)]./[L L],'rows')';

function u = il_D(x,y)
u = bsxfun(@minus,x(:,1:end-1),y).*bsxfun(@minus,x(:,2:end),y);

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function [waveform,freq_array] = il_Profile_Analysis(dur, rampdur, ncomponents, dB_incr, stim_dB, Fs, freq_array)
% function [waveform,freq_array] = il_Profile_Analysis(dur, rampdur, ncomponents, dB_incr, stim_dB, Fs, freq_array)
% 
% This function produces a vector of log-spaced components based on Lentz (JASA 2005)
%Variables:
% dur = duration in sec  (e.g. 0.5 sec)
% ncomponents is # of components in stimulus (spaced between 200 & 5000 Hz)
% dB_incr = 20*log(Delta_A/A), where A is the amplitude of the standard component at 1000 Hz
% ramp_time: onset/offset in seconds
% stim_dB = overall stimulus level in dB SPL (both inetervals are scaled to same dB SPL)
% Fs = sampling rate (Hz)

npts = dur * Fs; % # pts in stimulus
t = (0:(npts-1))/Fs; % time vector

if nargin < 7
    freq_array = logspace(log10(200),log10(5000),ncomponents); %column vector of n equally spaced compenents in the frequency range 200 to 5000Hz
end
waveform = 0;  % initialize variable

for f = freq_array  % step through each frequency component in the complex
    phase = 2*pi * rand(1,1); %Randomly vary the starting phase of each component
    if mod(f,4) ~= 0     %rounding each component to the nearest 4Hz (see Lentz, 2005)
        freq =  4*(ceil(f/4));
    else
        freq = f;
    end
    % Note that the increment is a tone added in phase to the component at 1000 Hz
    incr_size = 10.^(dB_incr/20); % convert increment from dB into a linear scalar  (if dB_incr = 0, incr_size = 1 and the amplitude at 1000 Hz is doubled)
    
    %Adding components with or without incr
    if freq == 1000
        waveform = waveform + (1 + incr_size) * cos(2 * pi * freq * t + phase);
    else
        waveform = waveform + cos(2 * pi * freq * t + phase);
    end
end

waveform = tukeywin(npts, 2*rampdur/dur)' .* waveform; % gate
waveform = 20e-6 * 10.^(stim_dB/20) * waveform/rms(waveform);% convert signal into Pascals

