function data = exp_dau1997(varargin)
%EXP_DAU1997 - 
%
%   Usage: data = exp_dau1997(flags)
%
%   exp_dau1997 reproduces Figs. 4 and 14 from Osses et al. (2020), where the 
%   dau1997 model was used. The figures are similar to Figs. 4.14, C.9B, 
%   and C.11B from Osses (2018).
%
%
%   The following flags can be specified;
%
%     'redo'    Recomputes data for specified figure
%
%     'plot'    Plot the output of the experiment. This is the default.
%
%     'noplot'  Don't plot, only return data.
%
%     'fig4_osses2020'    Reproduce Fig. 4 of Osses et al. (2020).
%
%     'fig14_osses2020'    Reproduce Fig. 14 of  Osses et al. (2020).
%
%   Fig. 4 - Two internal representations of a piano sound ('P1') using the  
%   PEMO model with two configurations of the adaptation loops are shown:
%   Overshoot limitation with a factor of 5, as suggested in Osses et al. (2020), and 
%   with a factor of 10 (see, Dau et al., 1997).
%   To display Fig. 4 of Osses et al. (2020) use :
%
%     out = exp_dau1997('fig4_osses2020');
%
%   Fig. 14 - The effect of the overshoot limitation with factors of 5 and 10
%   are shown for a 4-kHz pure tone of 70 dB SPL that includes 2.5-ms up/down 
%   ramps. For these plots the outer and middle ear stages are skipped. One
%   gammatone filter at 4 kHz is used, followed by the ihc stage (ihc_breebaart),
%   and the adaptation loops (adt_osses2020 for lim=5, adt_dau for lim=10).
%   To display Fig. 14 of Osses et al. (2020) use :
%
%     out = exp_dau1997('fig14_osses2020');
%
%   References:
%     T. Dau, B. Kollmeier, and A. Kohlrausch. Modeling auditory processing
%     of amplitude modulation. I. Detection and masking with narrow-band
%     carriers. J. Acoust. Soc. Am., 102:2892--2905, 1997a.
%     
%
%   Url: http://amtoolbox.sourceforge.net/amt-0.10.0/doc/experiments/exp_dau1997.php

% Copyright (C) 2009-2020 Piotr Majdak and the AMT team.
% This file is part of Auditory Modeling Toolbox (AMT) version 0.10.0
%
% This program is free software: you can redistribute it and/or modify
% it under the terms of the GNU General Public License as published by
% the Free Software Foundation, either version 3 of the License, or
% (at your option) any later version.
%
% This program is distributed in the hope that it will be useful,
% but WITHOUT ANY WARRANTY; without even the implied warranty of
% MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
% GNU General Public License for more details.
%
% You should have received a copy of the GNU General Public License
% along with this program.  If not, see <http://www.gnu.org/licenses/>.

%   AUTHOR: Alejandro Osses

data = [];

definput.import={'amtredofile'};
definput.flags.type={'missingflag','fig4_osses2020','fig14_osses2020'};

definput.flags.plot={'plot','noplot'};

[flags,keyvals]  = ltfatarghelper({},definput,varargin);


if flags.do_missingflag
  flagnames=[sprintf('%s, ',definput.flags.type{2:end-2}),...
             sprintf('%s or %s',definput.flags.type{end-1},definput.flags.type{end})];
  error('%s: You must specify one of the following flags: %s.',upper(mfilename),flagnames);
end;

%% ------ FIG 4 Osses and Kohlrausch 2020 ---------------------------------
if flags.do_fig4_osses2020

    % Goes to 'http://sofacoustics.org/data/amt-0.10.0/auxdata/' and loads
    %   the data under the folder 'dau1997', with name 'P1-GH05-Cd5_1-dur-1300-ms.wav'
    [insig, fs] = amt_load('dau1997','P1-GH05-Cd5_1-dur-1300-ms.wav');
    
    tobs = 0.25; % s, only the first 0.25 s of the waveform
    insig = insig(1:fs*tobs);
    
    subfs = 16000; % sampling frequency for the internal representation
    [outsig05,fc,mfc] = dau1997_preproc(insig,fs, ...
        'subfs',subfs,'gtf_osses2020','ihc_breebaart','adt_osses2020','mfb_jepsen2008');
                                
    outsig10 = dau1997_preproc(insig,fs, ...
        'subfs',subfs,'gtf_osses2020','ihc_breebaart','adt_dau','mfb_jepsen2008');
     
    N = length(fc); % number of audio frequencies
    K = length(mfc);
    %%% Memory allocation:
    Ik_05 = zeros(1,K); % one value for each modulation frequency
    Ik_10 = zeros(1,K); % one value for each modulation frequency
    
    Im_05 = zeros(1,N);  % one value for each audio frequency
    Im_10 = zeros(1,N);  % one value for each audio frequency
    %%%
    
    for j = 1:N
        Im_05(j) = Im_05(j)+1/subfs*sum(outsig05{j}(:).^2); % all mod. filters together
        Im_10(j) = Im_10(j)+1/subfs*sum(outsig10{j}(:).^2);    
    end
    
    for i = 1:K
        for j = 1:N
            Nr_mod_filters = size(outsig05{j},2);
            if i <= Nr_mod_filters 
                % summed up only if the mod filter i is present.
                Ik_05(i) = Ik_05(i)+1/subfs*sum(outsig05{j}(:,i).^2);  
                Ik_10(i) = Ik_10(i)+1/subfs*sum(outsig10{j}(:,i).^2);
            else
                % Nothing to do, in this case mfc(i) > 1/4*fc(j)
            end  
        end
    end
    
    Itot_05 = sum(Ik_05);
    Itot_10 = sum(Ik_10);
    
    fc_erb = freqtoaud(fc);
    mfc_nr = 1:K;
    
    if flags.do_plot
        
        % Panel A
        figure;
        plot(fc_erb,100*Im_05/Itot_05,'ro-'); hold on; grid on;
        plot(fc_erb,100*Im_10/Itot_10,'bs--');

        set(gca,'XTick',3:33);
        set(gca,'YTick',0:1:10);
        ylim([-1 11])
        xlim([2 34])
        legend('lim=5 (Osses2020)','lim=10');
        Pos = get(gcf,'Position');
        Pos(3) = 800;
        set(gcf,'Position',Pos); % setting width of the figure

        xlabel('Audio centre frequency f_c (ERB_N)');
        ylabel('Percentage (%)');

        title('A. Information-based audio-frequency analysis')
        
        % Panel B
        figure;
        plot(mfc_nr,100*Ik_05/Itot_05,'ro-'); hold on; grid on;
        plot(mfc_nr,100*Ik_10/Itot_10,'bs--'); 

        set(gca,'XTick',1:12);
        set(gca,'YTick',0:3:24);
        xlim([0 13])
        ylim([-1 27])
        legend('lim=5 (Osses2020)','lim=10');

        xlabel('Modulation centre frequency mf_c (Nr.)');
        ylabel('Percentage (%)');
    
        title('B. Information-based modulation-frequency analysis')
    end
    
    data.figure_flag = 'do_fig4_osses2020';
    data.fc = fc;
    data.mfc = mfc;
    data.Ik_05 = Ik_05;
    data.Ik_10 = Ik_10;
    data.Ik = 'Model Units (MU)';
    data.Im_05 = Im_05;
    data.Im_10 = Im_10;
    data.Im_unit = 'Model Units (MU)';
    data.Itot_05 = Itot_05;
    data.Itot_10 = Itot_10;
    data.Itot_unit = 'Model Units (MU)';
    data.description = 'Energy content (MU) for each audio frequency band n, and each modulation frequency band k';
end
  
%% ------ FIG 14 Osses and Kohlrausch 2020 --------------------------------
if flags.do_fig14_osses2020
    
    % 1. Stimulus creation:
    fs = 44100;
    dur = 300e-3; % 300 ms
    lvl = 70;
    dBFS = 100; % 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 = setdbspl(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
    t = (1:length(insig))/fs; t = t(:); % creates 't' as a column array
          
    % No outer, no middle ear:
    outsig = auditoryfilterbank(insig,fs,'basef',fc,'flow',fc,'fhigh',fc);
    % 'haircell' envelope extraction
    outsig = ihcenvelope(outsig,fs,'ihc_breebaart');

    % non-linear adaptation loops
    outsig05 = adaptloop(outsig,fs,'adt_osses2020');
    outsig10 = adaptloop(outsig,fs,'adt_dau');

    onset05 = max(outsig05);
    onset10 = max(outsig10);
    
    
    id_steady = find(t>=.330-dur_ramp_ms*1e-3 & t<=.350-dur_ramp_ms*1e-3);

    steady05 = mean(outsig05(id_steady));
    steady10 = mean(outsig10(id_steady));
    
    ra05 = onset05/steady05;
    ra10 = onset10/steady10;
    
    fprintf('Lim  5: Onset = %.1f MU, steady = %.1f MU\n',onset05,steady05);
    fprintf('Lim 10: Onset = %.1f MU, steady = %.1f MU\n',onset10,steady10);
    
    leg4plot{1} = sprintf('lim =  5: ratio onset/steady=%.1f',ra05);
    leg4plot{2} = sprintf('lim = 10: ratio onset/steady=%.1f',ra10);
    
    if flags.do_plot
    
        figure; 
        plot(t,outsig05,'r-','LineWidth',2); hold on, grid on;
        plot(t,outsig10,'b-')
       
        xlabel('Time (s)');
        ylabel('Amplitude \Psi (Model Units)')
        legend(leg4plot);
        
        set(gca,'XTick',[50:50:450]*1e-3);
        xlim([0.025 0.475]);
        ylim([-300 1480])
        set(gca,'YTick',-200:100:1400)
    end
    
    data.figure_flag = 'do_fig14_osses2020';
    data.ra05 = ra05;
    data.onset05 = onset05;
    data.steady05 = steady05;
    data.ra10 = ra10;
    data.onset10 = onset10;
    data.steady10 = steady10;    
    data.onset_unit = 'Model Units (MU)';
    data.steady_unit = 'Model Units (MU)';
end