function [Ldn,E,Cam,CF] = chen2011(inputF,inputLdB,varargin)
%CHEN2011  Fast excitation pattern estimation
%   Usage: [Ldn,E,Cam,CF] = chen2011(inputF,inputLdB)
%
%
%
%   Input parameters:
%     inputF           : vector of input frequency values
%     inputLdB         : vector of input level
%
%   Output parameters:
%     Ldn              : loudness
%     E                : excitation at the output of the passive filter
%     Cam              : ERB-number in Cams (Glasberg and Moore, 1990; Moore, 2003)
%     CF               : center frequency of auditory filters
%
%   Optional key value pairs:
%
%     'HLcf',HLcf             parameters to simulate hearing loss - audiogram frequencies
%
%     'HLohcdB0',HLohcdB0     parameters to simulate hearing loss - OHC loss at audiogram frequencies
%
%     'HLihcdB0',HLihcdB0     parameters to simulate hearing loss - IHC loss at audiogram frequencies
%
%     'cambin',cambin         spacing [ERB] between successive auditory filter CFs
%
%     'flow',flow             lowest center frequency of an auditory filter
%
%     'fhigh',fhigh           highest center frequency of an auditory filter
%
%   Optional flags:
%
%     'outerearcorrection'     outer ear correction ('FreeField', 'PDR10', 'Eardrum')
%
%   See also: demo_chen2011 exp_chen2011
%
%   References:
%     Z. Chen, G. Hu, B. R. Glasberg, and C. Moore, Brian. A new model for
%     calculating auditory excitation patterns and loudness for cases of
%     cochlear hearing loss. J. Acoust. Soc. Am., 282(1), 2011.
%     
%     B. R. Glasberg and B. C. J. Moore. A Model of Loudness Applicable to
%     Time-Varying Sounds. J. Audio Eng. Soc, 50(5):331--342, 2002.
%     
%     B. C. J. Moore, B. R. Glasberg, and T. Baer. A Model for the Prediction
%     of Thresholds, Loudness, and Partial Loudness. J. Audio Eng. Soc,
%     45(4):224--240, 1997.
%     
%
%   Url: http://amtoolbox.org/amt-1.1.0/doc/models/chen2011.php

% Copyright (C) 2009-2021 Piotr Majdak, Clara Hollomey, and the AMT team.
% This file is part of Auditory Modeling Toolbox (AMT) version 1.1.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: Perfect
%   #Verification: Verified
%   #Author : Zhangli Chen: original code
%   #Author : Clara Hollomey (2020): integration in the AMT


if length(inputF) ~= length(inputLdB)
    error('inputF and inputLdB should be dB/Hz and have same length');
end

definput.import = {'chen2011'}; % load defaults from arg_chen2011

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

%calculation of auditory filter CF
Cam = (fc2erb(kv.flow):kv.cambin:fc2erb(kv.fhigh))';     % 40-17000Hz
CF = 1000*(10.^(Cam/21.4)-1)/4.37;   % Hz

if length(kv.HLcf) ~= length(kv.HLohcdB0) || length(kv.HLcf) ~= length(kv.HLihcdB0)
    error('HLcf, HLohcdB and HLihcdB should have same length');
elseif length(kv.HLcf) == 1
    HLohcdB = kv.HLohcdB0 * ones(length(CF),1);
    HLihcdB = kv.HLihcdB0 * ones(length(CF),1);
else
    HLohcdB = interp1(kv.HLcf,kv.HLohcdB0,CF,'linear','extrap');
    HLihcdB = interp1(kv.HLcf,kv.HLihcdB0,CF,'linear','extrap');
end
HLohcdB = max(HLohcdB,0);
HLihcdB = max(HLihcdB,0);


%%%%%%%%%%%%%%%%%%%%% step1: outer ear correction

if strcmp(flags.outerearcorrection,'FreeField')
    freefield_F = [0 20 25 31.5 40 50 63 80 100 ...
    125 160 200 250 315 400 500 630 ...
    750 800 1000 1250 1500 1600 2000 2500 ...
    3000 3150 4000 5000 6000 6300 8000 9000 ...
    10000 11200 12500 14000 15000 16000 20000];
    freefield_dB = [0 0 0 0 0 0 0 0 0 ...
    0.1 0.3 0.5 0.9 1.4 1.6 1.7 2.5 ...
    2.7 2.6 2.6 3.2 5.2 6.6 12 16.8 ...
    15.3 15.2 14.2 10.7 7.1 6.4 1.8 -0.9 ...
    -1.6 1.9 4.9 2 -2 2.5 2.5];

    inputLdB = inputLdB + interp1(freefield_F,freefield_dB,inputF,'linear','extrap');
    
elseif strcmp(flags.outerearcorrection,'PDR10')
    PDR10_F =  [100  200  300  500 700  1000 1500 2000 2500 3000 3500 4000 5000 6000 7000];
    PDR10_dB = [-19.0 -9.5 -4.5 0 0.5 0 1.5  3.0  3.0  3.0  4.0  3.0  1.5  -3.0  -28.0];
    
    inputLdB = inputLdB + interp1(PDR10_F,PDR10_dB,inputF,'linear','extrap');
    
elseif strcmp(flags.outerearcorrection,'Eardrum')
else
    error('No such correction');
end

%%%%%%%%%%%%%%%%%%%%% step2: middle ear correction
MidEar_F = [20 25 31.5 40 50 63 80 100 125 ...
            160 200 250 315 400 500 630 750 800 1000 ...
            1250 1500 1600 2000 2500 3000 3150 4000 5000 6000 ...
            6300 8000 9000 10000 11200 12500 14000 15000 16000 20000];
MidEar_dB = [-33.2 -28.2 -23.2 -19.4 -16.3 -13.3 -10.2 -8.0 -6.1 ...
            -4.7 -3.5 -2.8 -2.4 -1.9 -1.8 -2.1 -2.5 -2.3 -2.6 ...
            -3.7 -5.5 -6.7 -11.4 -14.5 -11.5 -11.0 -10.5 -10.8 -12.8 ...
            -13.6 -16.5 -15.8 -15.0 -16.9 -18.8 -20.7 -21.9 -22.3 -24.1];

inputLdB = inputLdB + interp1(MidEar_F,MidEar_dB,inputF,'linear','extrap');

inputL = 10.^(inputLdB/10);                 % Intensity in linear

%%%%%%%%%%%%%%%%%%%%% step3: passive filter
tl = CF./(0.1084*CF+2.3301);
tu = 15.0 * ones(length(CF),1);

E_pf = zeros(length(CF),1);
for i = 1:length(CF)
    g = inputF/CF(i)-1;
    indexl = find(g<0);
    gl = abs(g(indexl));
    E_pf(i) = sum((1+gl*tl(i)).*exp(-gl*tl(i)).*inputL(indexl));
    indexu = find(g>=0);
    gu = g(indexu);
    E_pf(i) = E_pf(i) + sum((1+gu*tu(i)).*exp(-gu*tu(i)).*inputL(indexu));
end
E_pf = max(E_pf,10^(-10));
EdB_pf = 10*log10(E_pf);

EdB_pf = max(EdB_pf,0);    

%%%%%%%%%%%%%%%%%%%% step4: gain decided by passive input
GdBmax = CF./(0.0191*CF+1.1) - HLohcdB;

GdB = GdBmax.*( 1 - 1./(1+exp(-0.05*(EdB_pf-(100-GdBmax)))) + 1./(1+exp(-0.05*(0-(100-GdBmax)))));   
index = find(EdB_pf>30);
GdB(index) = GdB(index) - 0.003*(EdB_pf(index)-30).^2;

GdB = min(max(GdB,-20),GdBmax);
G = 10.^(GdB/10); 

%%%%%%%%%%%%%%%%%%%% step5: active tip filter (af)
pl = CF./(0.0272*CF+5.4365);
pu = 27.9 * ones(length(CF),1);

E_af = zeros(length(CF),1);
for i = 1:length(CF)
    g = inputF/CF(i)-1;
    indexl = find(g<0);
    gl = abs(g(indexl));
    E_af(i) = G(i) * sum((1+gl*pl(i)).*exp(-gl*pl(i)).*inputL(indexl));
    indexu = find(g>=0);
    gu = g(indexu);
    E_af(i) = E_af(i) + G(i) * sum((1+gu*pu(i)).*exp(-gu*pu(i)).*inputL(indexu));
end

E = E_pf + E_af;
E= max(E,10^(-10));
EdB = 10*log10(E);

EdB = EdB- HLihcdB.*(1-0.5./(1+exp(-0.2*((EdB-52)-(HLihcdB+20))))); 
E = 10.^(EdB/10);

Ldn = sum(E) * kv.cambin * 1.525*1e-8;