function [r_mean,psth,ihc,c1,c2,r_var] = zilany2014(stim,fsstim,fc,varargin)
%ZILANY2014 Auditory-nerve filterbank (improved)
%
%   Usage: [ANresp,fc] = zilany2014(lvl,stim,fsstim);
%          [ANresp,fc,psth,ihc,t_var] = zilany2014(lvl,stim,fsstim);
%          [ANresp,fc,psth,ihc,t_var] = zilany2014(lvl,stim,fsstim); 
%
%   Input parameters:
%     stim        : Pressure waveform of stimulus (timeseries)
%     fsstim      : Sampling frequency of stimulus  
%     fc         : Frequency vector containing the CFs. 
%                  Use fc=audspace(lo,hi,numCF,'erb'); to space equally on the 
%                  ERB frequency scale. 
%
%   Output parameters:
%     r_mean     : Instantaneous mean spiking rate (incl. refractoriness) 
%                  of different AN fibers at corresponding CFs. Size: [time CFs]
%     psth       : Spike histogram
%     ihc        : Output from inner hair cells (IHCs) in Volts
%     c1         : Output from the chirping filter C1
%     c2         : Output from the chirping filter C2
%     r_var      : Instananeous variance in the discharge rate of the ANs
%
%
%   This function takes the following optional key/value pairs:
%
%     'fsmod',fsmod   Model sampling rate. It is possible to run the model 
%                     at a range of fsmod between 100 kHz and 500 kHz.
%                     Default value is 200kHz for cats and 100kHz for humans.
%
%     'fiberType',fT  Type of the fiber based on spontaneous rate (SR)
%                     1: Low SR, SR fixed to 0.1 spikes/s
%                     2: Medium SR, SR fixed to 4 spikes/s
%                     3: High SR, SR fixed to 100 spikes/s
%                     4: Custom, defined by the fibre numbers in numH, numM and numL
%
%     'numH'          Number of high SR fibres. Only if fiberType is 4. 
%     'numM'          Number of medium SR fibres. Only if fiberType is 4. 
%     'numL'          Number of low SR fibres. Only if fiberType is 4.
%
%     'cohc',cohc     OHC scaling factor: 1 denotes normal OHC function (default); 
%                     0 denotes complete OHC dysfunction.
%
%     'cihc',cihc     IHC scaling factor: 1 denotes normal IHC function (default); 
%                     0 denotes complete IHC dysfunction.
%
%     'nrep',nrep     Number of repetitions for the mean rate, 
%                     rate variance & psth calculation. 
%
%   ZILANY2014 accepts the following flag:
%
%     'human'         Use model parameters for humans. This is the default.
%
%     'cat'           Use model parameters for cats.
%
%     'fixedFGn'      Fractional Gaussian noise will be the same in every 
%                     simulation. This is the default.
%
%     'varFGn'        Fractional Gaussian noise will be different in every 
%                     simulation.
%
%     'approxPL'      Use approxiate implementation of the power-law
%                     functions. This is the default.
%
%     'actualPL'      Use actual implementation of the power-law functions.
%
%     'shera2002'     Selects the BM tuning from Shera et al. (2002) (default).
% 
%     'glasberg1990'  Selects the BM tuning from Glasberg & Moore (1990)
%
%   ZILANY2014(...) returns modeled responses of multiple AN fibers tuned to 
%   various characteristic frequencies (CFs). Middle-ear filtering is 
%   included and corresponds to middleearfilter(...,'zilany2009');
%   Please cite the references below if you use this model.
%
%   See also: demo_zilany2014 zilany2014_synapse zilany2014_innerhaircells
%             zilany2014_ffgn zilany2014 zilany2007
%             data_baumgartner2016 plot_roenne2012 plot_roenne2012_chirp
%             plot_roenne2012_tonebursts demo_carney2015 baumgartner2016_spectralanalysis
%             roenne2012_click roenne2012_chirp carney2015_generateneurogram
%             roenne2012_tonebursts exp_osses2022
%             bruce2018 roenne2012 baumgartner2013
%             baumgartner2016
%
%   Demos: demo_zilany2014
%
%   References:
%     C. A. Shera, J. J. J. Guinan, and O. A. J. Revised estimates of human
%     cochlear tuning from otoacoustic and behavioral measurements.
%     Proceedings of the National Academy of Sciences of the United States of
%     America, 99(5):3318--3323, 2002.
%     
%     B. R. Glasberg and B. Moore. Derivation of auditory filter shapes from
%     notched-noise data. Hearing Research, 47(1-2):103--138, 1990.
%     
%     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.
%     
%     M. Zilany, I. Bruce, P. Nelson, and L. Carney. A phenomenological model
%     of the synapse between the inner hair cell and auditory nerve:
%     Long-term adaptation with power-law dynamics. J. Acoust. Soc. Am.,
%     126(5):2390 -- 2412, 2009.
%     
%
%   Url: http://amtoolbox.sourceforge.net/amt-0.10.0/doc/models/zilany2014.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/>.

%   #StatusDoc: Good
%   #StatusCode: Good
%   #Verification: Unknown
%   #Requirements: MATLAB MEX M-Signal
%   #Author: Muhammad Zilany 
%   #Author: Robert Baumgartner: adapted to the AMT
%   #Author: Piotr Majdak (2021): C1 and C2 outputs

if nargin<3
  error('%s: Too few input parameters.',upper(mfilename));
end;


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Loading default values, reading input arguments:

% Define input flags and values
definput.import = {'zilany2014'}; % load defaults from arg_zilany2014

[flags,kv]  = ltfatarghelper({'fsmod','fiberType','cohc','cihc','nrep'},definput,varargin);

% Species settings
if flags.do_cat, kv.fsmod = max(kv.fsmod,200e3); end

% Stimulus settings
stim	= resample(stim,kv.fsmod,fsstim); % stim fs = mod fs
fs = kv.fsmod;
tdres   = 1/fs; % time interval of processing (in s), i.e., reciprocal of the sampling rate

stim = stim(:)'; % stim will be a row vector

numCF=length(fc);

% reptime is the time between stimulus repetitions in seconds; -> set 
% twice the duration of stim
reptime = kv.reptime*length(stim)/fs;

% species must be 1 for cat, 2 for human (shera2002) or 3 for human (glasberg1990)
if flags.do_cat, species=1; end
if flags.do_human
  if flags.do_shera2002, species=2; end
  if flags.do_glasberg1990, species=3; end
end

% noiseType is for fixed fGn (noise will be same in every simulation) or variable fGn: "0" for fixed fGn and "1" for variable fGn
noiseType = flags.do_varFGn + 0;

% implnt is for "approxiate" or "actual" implementation of the power-law functions: "0" for approx. and "1" for actual implementation
implnt = flags.do_actualPL + 0;

% Call AN model and loop for all fibers tuned to different CFs
ihc = zeros(round(reptime*fs)*kv.nrep,numCF);
r_mean = zeros(round(reptime*fs),numCF);
r_var = zeros(round(reptime*fs),numCF);
psth = zeros(round(reptime*fs),numCF);
if nargout>=5,
  c1 = zeros(round(reptime*fs*kv.nrep),numCF);
  c2 = zeros(round(reptime*fs*kv.nrep),numCF);
end

if kv.fiberType==4,
  fiberTypes=[ones(kv.numL,1);ones(kv.numM,1)*2;ones(kv.numH,1)*3]; % custom fiber types
else
  fiberTypes=kv.fiberType;
end
for jj = 1:numCF
  if nargout>=5,
    [ihc(:,jj),c1(:,jj),c2(:,jj)] = zilany2014_innerhaircells(stim,fc(jj),kv.nrep,tdres,reptime,kv.cohc,kv.cihc,species);
  else
    ihc(:,jj) = zilany2014_innerhaircells(stim,fc(jj),kv.nrep,tdres,reptime,kv.cohc,kv.cihc,species);    
  end

  for ii=1:length(fiberTypes)
	amt_disp(['CF = ' int2str(jj) '/' int2str(numCF) '; spont = ' int2str(ii) '/' int2str(length(fiberTypes))],'volatile');
    [r_mean1,r_var1,psth1] = zilany2014_synapse(ihc(:,jj)',fc(jj),kv.nrep,tdres,fiberTypes(ii),noiseType,implnt);
    r_mean(:,jj)=r_mean(:,jj)+r_mean1;
    r_var(:,jj)=r_var(:,jj)+r_var1;
    psth(:,jj)=psth(:,jj)+psth1;
  end
  amt_disp();
  if length(fiberTypes)>1
    r_mean(:,jj)=r_mean(:,jj)/length(fiberTypes);
    r_var=r_var/length(fiberTypes);
    psth(:,jj)=psth(:,jj)/length(fiberTypes);
  end
end

if ~isempty(kv.psth_binwidth),
  psth_bin = zeros(size(buffer(psth(:,1),round(kv.psth_binwidth*fs)),2),numCF);
  for ii=1:numCF
    psth_bin(:,ii)=sum(buffer(psth(:,ii),round(kv.psth_binwidth*fs)))/kv.nrep/kv.psth_binwidth;
  end
  psth=psth_bin;
end