function exp_workshop2022paris(varargin)
% function exp_workshop2022paris(varargin)
%
% This script processes two waveforms using the model specified by the 
%   corresponding flag. The structure of this script is similar to those
%   adopted in AMT to reproduce data of the literature (e.g., exp_lopezpoveda2001,
%   exp_osses2022, exp_takanen2013, etc.).
%
% To run this script you can run:
%   exp_workshop2022paris('dau1997');          
%   exp_workshop2022paris('zilany2014');       
%   exp_workshop2022paris('verhulst2015');
%   exp_workshop2022paris('verhulst2018');
%   exp_workshop2022paris('bruce2018');        
%   exp_workshop2022paris('relanoiborra2019'); 
%   exp_workshop2022paris('king2019');         
%   exp_workshop2022paris('osses2021');        
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

% We will create one flag called 'model_to_use'. If undefined, you will get
%   and error.
definput.flags.model_to_use = {'unspecified','dau1997','zilany2014', ...
    'verhulst2015','verhulst2018','bruce2018','relanoiborra2019','king2019','osses2021'};
[flags,keyvals]  = ltfatarghelper({},definput,varargin);
% keyvals: in this script keyvals is unused... I still included it as an 
%   output parameter because it is done this one in every script used in AMT.

[insig,fs] = greasy; % Loading wav file from the LTFAT toolbox

if flags.do_unspecified
    help exp_workshop2022paris;
    error('See the help of this script and specify the model you want to adopt')
end
%% For each model I will give immediately get to the last step of the workshop activity:
basef = 200; % Hz, this is the frequency at which we request the plots later...

if flags.do_dau1997
    flow = 80;    % Hz, lower  CF of the filter bank (change it to basef if you want to get one CF only)
    fhigh = 8000; % Hz, higher CF of the filter bank (change it to basef if you want to get one CF only)
    subfs = fs;
    
    flags_to_use = {'basef',basef,'flow',flow,'fhigh',fhigh,'subfs',subfs};
    [outsig,fc,mfc] = dau1997(insig,fs,flags_to_use{:});
    
    idx_cf = find(fc>=basef,1,'first'); % gets the idx of the band containing the basef
    
    t = (1:size(outsig{1},1))/subfs; % subfs is the sampling frequency of the internal representation, by default equal to fs
    N_mod_channels = size(outsig{1},2); % number of modulation channels for this band
    
    text4leg = [];
    figure;
    for i = 1:N_mod_channels
        plot(t,outsig{idx_cf}(:,i)); hold on;
        text4leg{end+1} = sprintf('mfc=%.1f Hz',mfc(i));
    end
    legend(text4leg);
    xlabel('Time (s)');
    ylabel('Amplitude (MU)');
    title(sprintf('dau1997: Mod. filter bank output\n for the audio frequency band at %.1f Hz',fc(idx_cf)));
    % Extra comments about this model:
    %   - mfc(1) = 0 Hz is actually the low-pass modulation filter at mfc = 2.5 Hz
    %   - The amplitude of the representations is expressed in model units (MU), 
    %     which is an arbitrary unit, but extensively used during the development
    %     of the dau1997 model (before and after that publication).
    
    % Now, to by-pass the modulation filter we need to identify the name of
    %   the model stage (or module): the name here is 'mfb', so we need to indicate 'no_mfb'
    
    [outsig,fc] = dau1997(insig,fs,flags_to_use{:},'no_mfb');
    %   - fc adopts the same values as before, so I'll directly use idx_cf
    %   - In this case, outsig is a regular Matrix array
    
    figure;
    plot(t,outsig(:,idx_cf)); 
    
    xlabel('Time (s)');
    ylabel('Amplitude (MU)');
    title(sprintf('dau1997: Adaptation loops output\n for the audio frequency band at %.1f Hz',fc(idx_cf)));
    
end

if flags.do_zilany2014
    % On purpose I am indicating arbitrary flow and fhigh values. To do the
    %   exercise we would only need one CF (basef could be used as input directly)
    %   It is more interesting to show the multichannel approach here
    flow = 125; % Hz
    fhigh = 2000; % Hz
    fc = audspacebw(flow,fhigh,'hitme',basef); % You can set any input CF
    
    flags_to_use = {'human'}; % Making sure we will use the human (and not the cat settings)
    [r_mean,psth,ihc,c1,c2,r_var,output] = zilany2014(insig,fs,fc,flags_to_use{:});
    subfs = 100e3; % Hz, zilany2014 doesn't return the simulated sampling frequency
                   % For simulation accuracy, the authors recommend to adopt a high
                   % sampling frequency of 100 kHz for human settings, 200 kHz for cat settings
    idx_cf = find(fc>=basef,1,'first');
    
    t = (1:size(r_mean,1))/subfs;
    figure;
    plot(t,r_mean(:,idx_cf));
    
    xlabel('Time (s)');
    ylabel('Mean-rate firing rate (spikes/s)');
    title(sprintf('zilany2014: Auditory adaptation stage output\n for the audio frequency band at %.1f Hz',fc(idx_cf)));
        
    % Extra comments about this model:
    %   - This model does not have a modulation processing stage by default, 
    %     so I skipped the first plot.
    %   - By default the model zero pads the input signal and simulates a signal
    %     twice as long as the input signal. This can be shortened by adopting
    %     a specific keyvalue, but the authors do not recommend that because
    %     they believe that the adaptation effects after the signal offset can
    %     introduce relevant signal cues.
    
end

if flags.do_bruce2018
    % On purpose I am indicating arbitrary flow and fhigh values. To do the
    %   exercise we would only need one CF (basef could be used as input directly)
    %   It is more interesting to show the multichannel approach here
    flow = 125; % Hz
    fhigh = 2000; % Hz
    fc = audspacebw(flow,fhigh,'hitme',basef); % You can set any input CF
    
    flags_to_use = {'human'}; % Making sure we will use the human (and not the cat settings)
    [output, info] = bruce2018(insig,fs,fc,flags_to_use{:});
    subfs = info.keyvals.fsmod; % fs of the simulations
	r_mean = output.meanrate;
    idx_cf = find(fc>=basef,1,'first');
    
    t = (1:size(r_mean,1))/subfs;
    figure;
    plot(t,r_mean(:,idx_cf));
    
    xlabel('Time (s)');
    ylabel('Mean-rate firing rate (spikes/s)');
    title(sprintf('bruce2018: Auditory adaptation stage output\n for the audio frequency band at %.1f Hz',fc(idx_cf)));
        
    % Extra comments about this model:
    %   - This model does not have a modulation processing stage by default, 
    %     so I skipped the first plot.
    %   - By default the model zero pads the input signal. This zero padding
    %     can be shortened by adopting a specific keyvalue, but the authors 
    %     do not recommend that because they believe that the adaptation 
    %     effects after the signal offset can introduce relevant signal cues.
end

if flags.do_verhulst2015
    % fc_flag = 'abr'; % if fc_flag is a string ('abr','all','half') it returns 401, 1001, or 500 bands respectively
    fc_flag = basef; % Hz, it will return only one CF but the simulation will still be performed for many bands
    [output,fc] = verhulst2015(insig,fs,fc_flag);
    
    subfs    = output.fs_abr; % Hz
    subfs_AN = output.fs_an; % Hz
    
    t = (1:size(output.ic,1))/subfs;
    figure;
    plot(t,output. output.ic);
    
    xlabel('Time (s)');
    ylabel('IC potential, unscaled amplitude, proportional to (V)');
    title(sprintf('verhulst2015: Auditory adaptation stage output\n for the audio frequency band at %.1f Hz',fc_flag));
        
    t_AN = (1:size(output.an_summed,1))/subfs_AN;
    figure;
    plot(t_AN,output.an_summed);
    
    xlabel('Time (s)');
    ylabel('Auditory nerve mean firing rate (spikes/s)');
    title(sprintf('verhulst2015: Auditory adaptation stage output\n for the audio frequency band at %.1f Hz',fc_flag));
end

if flags.do_verhulst2018
    % fc_flag = 'abr'; % if fc_flag is a string ('abr','all','half') it returns 401, 1001, or 500 bands respectively
    fc_flag = basef; % Hz, it will return only one CF but the simulation will still be performed for many bands
    output = verhulst2018(insig,fs,fc_flag);
    
    subfs    = output.fs_abr; % Hz
    subfs_AN = output.fs_an; % Hz
    
    num_tot = output.keyvals.numH + output.keyvals.numM + output.keyvals.numL;
    
    t = (1:size(output.ic,1))/subfs;
    figure;
    plot(t,10e6*output.keyvals.M5*output.ic);
    
    xlabel('Time (s)');
    ylabel('IC potential (\muV)');
    title(sprintf('verhulst2018: Auditory adaptation stage output\n for the audio frequency band at %.1f Hz',fc_flag));
        
    t_AN = (1:size(output.an_summed,1))/subfs_AN;
    figure;
    plot(t_AN,output.an_summed/num_tot);
    
    xlabel('Time (s)');
    ylabel('Auditory nerve mean firing rate (spikes/s)');
    title(sprintf('verhulst2018: Auditory adaptation stage output\n for the audio frequency band at %.1f Hz',fc_flag));  
end

if flags.do_relanoiborra2019
     
    % This is basically the same code as for dau1997
    
    flow = 80;    % Hz, lower  CF of the filter bank (change it to basef if you want to get one CF only)
    fhigh = 8000; % Hz, higher CF of the filter bank (change it to basef if you want to get one CF only)
    
    flags_to_use = {'basef',basef,'flow',flow,'fhigh',fhigh, ...
        'erbspacebw', ...    % ERB space is a flag that allows to define fcs 'at will'
        'no_internalnoise'}; % but because CFs are 'at will' we need to disable the internal noise generator that
                             %   expects hard-coded frequencies
    [outsig_mfb, mfc, outsig_afb, fc] = relanoiborra2019_featureextraction(insig,fs,flags_to_use{:});
    
    idx_cf = find(fc>=basef,1,'first'); % gets the idx of the band containing the basef
    
    t = (1:size(outsig_mfb,1))/fs; % subfs is the sampling frequency of the internal representation, by default equal to fs
    
    idx_fc_mod = find(mfc<=fc(idx_cf)/4); % In this model modulation channels with mfc < fc / 4 should be discarded.
                    % For implementation reasons, this removal of information is performed in the decision (back-end) stage
                    % of this model (relanoiborra2019_decision.m) and NOT at the preprocessing stage, where it should
                    % be conceptually performed. For this reason, I am applying this information removal here.
    outsig_mfb = outsig_mfb(:,:,idx_fc_mod);
    
    N_mod_channels = size(outsig_mfb,3); % number of modulation channels for this band
    
    
    text4leg = [];
    figure;
    for i = 1:N_mod_channels
        plot(t,outsig_mfb(:,idx_cf,i)); hold on;
        text4leg{end+1} = sprintf('mfc=%.1f Hz',mfc(i));
    end
    legend(text4leg);
    xlabel('Time (s)');
    ylabel('Amplitude (MU)');
    title(sprintf('relanoiborra2019: Mod. filter bank output\n for the audio frequency band at %.1f Hz',fc(idx_cf)));
    % Extra comments about this model:
    %   - mfc(1) = 0 Hz is actually the low-pass modulation filter at mfc = 2.5 Hz
    %   - The amplitude of the representations is expressed in model units (MU), 
    %     which is an arbitrary unit, but extensively used during the development
    %     of the dau1997 model (before and after that publication).
    
    % Now, the modulation filter CANNOT be by-passed, but the output outsig_afb
    %   contains the output of the adaptation loops already.
    
    figure;
    plot(t,outsig_afb(:,idx_cf)); 
    
    xlabel('Time (s)');
    ylabel('Amplitude (MU)');
    title(sprintf('relanoiborra2019: Adaptation loops output\n for the audio frequency band at %.1f Hz',fc(idx_cf)));
        
end

if flags.do_king2019
    % This is very similar to the code for dau1997
    
    flow = 80;    % Hz, lower  CF of the filter bank (change it to basef if you want to get one CF only)
    fhigh = 8000; % Hz, higher CF of the filter bank (change it to basef if you want to get one CF only)
    subfs = fs;
    
    flags_to_use = {'basef',basef,'flow',flow,'fhigh',fhigh,'subfs',subfs};
    [outsig, fc, mfc, step] = king2019(insig,fs,flags_to_use{:});
    
    idx_cf = find(fc>=basef,1,'first'); % gets the idx of the band containing the basef
    
    t = (1:size(outsig,1))/subfs; % subfs is the sampling frequency of the internal representation, by default equal to fs
    N_mod_channels = size(outsig,3); % number of modulation channels for this band
    
    text4leg = [];
    figure;
    for i = 1:N_mod_channels
        plot(t,outsig(:,idx_cf,i)); hold on;
        text4leg{end+1} = sprintf('mfc=%.1f Hz',mfc(i));
    end
    legend(text4leg);
    xlabel('Time (s)');
    ylabel('Amplitude (a.u.)');
    title(sprintf('king2019: Mod. filter bank output\n for the audio frequency band at %.1f Hz',fc(idx_cf)));
    % Extra comments about this model:
    %   - The amplitude of the representations is expressed in arbitrary units (a.u.)
    
    % Now, to by-pass the modulation filter we need to identify the name of
    %   the model stage (or module): the name here is 'mfb', so we need to indicate 'no_mfb'
    outsig = king2019(insig,fs,flags_to_use{:},'no_mfb');
   
    figure;
    plot(t,outsig(:,idx_cf)); 
    
    xlabel('Time (s)');
    ylabel('Amplitude (a.u.)');
    title(sprintf('king2019: Adaptation loops output\n for the audio frequency band at %.1f Hz',fc(idx_cf)));
    
end

if flags.do_osses2021 
    % This is basically the same code as for dau1997
    
    flow = 80;    % Hz, lower  CF of the filter bank (change it to basef if you want to get one CF only)
    fhigh = 8000; % Hz, higher CF of the filter bank (change it to basef if you want to get one CF only)
    subfs = fs;
    
    flags_to_use = {'basef',basef,'flow',flow,'fhigh',fhigh,'subfs',subfs};
    [outsig,fc,mfc] = osses2021(insig,fs,flags_to_use{:});
    
    idx_cf = find(fc>=basef,1,'first'); % gets the idx of the band containing the basef
    
    t = (1:size(outsig{1},1))/subfs; % subfs is the sampling frequency of the internal representation, by default equal to fs
    N_mod_channels = size(outsig{1},2); % number of modulation channels for this band
    
    text4leg = [];
    figure;
    for i = 1:N_mod_channels
        plot(t,outsig{idx_cf}(:,i)); hold on;
        text4leg{end+1} = sprintf('mfc=%.1f Hz',mfc(i));
    end
    legend(text4leg);
    xlabel('Time (s)');
    ylabel('Amplitude (MU)');
    title(sprintf('osses2021: Mod. filter bank output\n for the audio frequency band at %.1f Hz',fc(idx_cf)));
    % Extra comments about this model:
    %   - mfc(1) = 0 Hz is actually the low-pass modulation filter at mfc = 2.5 Hz
    %   - The amplitude of the representations is expressed in model units (MU), 
    %     which is an arbitrary unit, but extensively used during the development
    %     of the dau1997 model (before and after that publication).
    
    % Now, to by-pass the modulation filter we need to identify the name of
    %   the model stage (or module): the name here is 'mfb', so we need to indicate 'no_mfb'
    
    [outsig,fc] = osses2021(insig,fs,flags_to_use{:},'no_mfb');
    %   - fc adopts the same values as before, so I'll directly use idx_cf
    %   - In this case, outsig is a regular Matrix array
    
    figure;
    plot(t,outsig(:,idx_cf)); 
    
    xlabel('Time (s)');
    ylabel('Amplitude (MU)');
    title(sprintf('osses2021: Adaptation loops output\n for the audio frequency band at %.1f Hz',fc(idx_cf)));
    
end