function results = moore1997(inSig,fs,varargin)
%MOORE1997  Loudness model for stationary signals
%   Usage: results = moore1997(inSig, fs);
%          results = moore1997(inSig, fs, ...);
%
%   Input parameters:
%     insig : Audio signal (single audio channel).
%     fs    : Sampling frequency (in Hz) of the input signal. 
%
%   Output parameters:
%     results: Structure containing the model results: 
%
%              - eLdB*: Excitation pattern, i.e., vector with the intensities
%                (in dB) according to Eq. 3 from Moore et al. (1997). 
%
%              - erbN*: Vector with the ERB rates (in Cam) of the 
%                considered auditory filters. 
% 
%              - fc*: Vector with the center frequencies (in Hz) of the
%                considered auditory filters. 
%
%              - specLoud*: Specific loudness as calculated in Sec. 1.6 of
%                Moore et al. (1997). 
%
%              - monauralLoudness*: the specific loudness specLoud integrated
%                over all auditory filters. The meaning of this output is unclear 
%                and remains here for backward compatibility only.
%
%   
%   Optional key-value pairs:
%
%     'fs',intfs  Model-internal sampling frequency (in Hz). The input signal
%                 insig will be resampled to intfs. Default: 32000 Hz.
%                 It is strongly recommended to not change this default. 
%
%     'flow',flow   Lowest frequency (in Hz) for the evaluation of the outer-
%                   and middle-ear transfer function. Default: 20 Hz.
%
%     'fhigh',fhigh   Highest frequency (in Hz) for the evaluation of the outer-
%                     and middle-ear transfer function. Default: 16000 Hz.
%
%     'order',order   Order of the FIR filter (in samples) used to derive the 
%                     outer- and middle-ear transfer function. Default: 4096 samples.
%
%     'erbStep',erbStep   Spacing (in Cam) between successive excitation patterns 
%                         on the ERB scale. 
%
%     'erbFcMin',erbFcMin   Center frequency (in Hz) of the lowest excitation pattern.
%
%     'erbFcMax',erbFcMax   Center frequency (in Hz) of the highest excitation pattern.
%
%   This function calculates the excitation patterns (see Sec. 1.4 in Moore et al., 1997) 
%   and the specific loudness in quiet (Sec. 1.6) as described in Moore et al. (1997). Note that
%   the calculation of the partial loudness (Secs. 1.7ff) is not implemented. 
%
%
%   See also: data_glasberg2002 exp_moore1997 glasberg2002
%
%   References:
%     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.6.0/doc/models/moore1997.php


%   #StatusDoc: OK
%   #StatusCode: OK
%   #Verification: Unknown
%   #Requirements: M-Signal
%   #Author: Thomas Deppisch (2017): Original implementation for the AMT.
%   #Author: Clara Hollomey (2021): Documentation added and tntegration in the AMT. 
%   #Authro: Piotr Majdak (2024): Issue #208 and major documentation upgrade. 

% This file is licensed unter the GNU General Public License (GPL) either 
% version 3 of the license, or any later version as published by the Free Software 
% Foundation. Details of the GPLv3 can be found in the AMT directory "licences" and 
% at <https://www.gnu.org/licenses/gpl-3.0.html>. 
% You can redistribute this file and/or modify it under the terms of the GPLv3. 
% This file is distributed without any warranty; without even the implied warranty 
% of merchantability or fitness for a particular purpose. 

%% model

definput.import = {'moore1997'}; % load from arg_moore1997

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

if fs ~= kv.fs
    inSig = resample(inSig, kv.fs, fs);
    fs = kv.fs;
end

fVec = kv.flow:kv.fhigh;   
data = data_glasberg2002('tfOuterMiddle1997','fVec',fVec);

% create FIR filter
tfLinear = 10.^(data.tfOuterMiddle/10);
outerMiddleFilter = fir2(kv.order, linspace(0, 1, length(fVec)), tfLinear);
earSig = filtfilt(outerMiddleFilter,1,inSig);   % why does filter(..) not work?

% compute fft
spect = fft(earSig);
fftLen = length(spect);
oneHz = (fftLen+2)/kv.fs;  % number of frequency bins representing 1Hz

numBins = round(fftLen/2+1);
compInt = 2*abs(spect(1:numBins)).^2/(numBins*fs);  % psd
compFq = linspace(0,fs/2,numBins);
nPoints = length(compFq);

% calculate ERB numbers corresponding to ERB mid frequencies    
erbNMin = f2erbrate(kv.erbFcMin);
erbNMax = f2erbrate(kv.erbFcMax);

erbN = erbNMin:kv.erbStep:erbNMax;    % numbers of erb bands
erbFc = erbrate2f(erbN);               % center frequency of erb bands

erbLoFreq = erbrate2f(erbN-0.5); % lower limit of each ERB filter
erbHiFreq = erbrate2f(erbN+0.5); % upper limit of each ERB filter

%calculate intensity for each ERB (dB/ERB)
erbInt = zeros(size(erbFc));
for ii=1:length(erbFc)
    range = round(erbLoFreq(ii)*oneHz):round(erbHiFreq(ii)*oneHz);
    erbInt(ii) = sum(compInt(range));   % intensity sum in each erb
end

erbdB = 10*log10(erbInt./(20e-6)^2);   % intensity level in each erb using reference SPL of 20 uPa
p511 = 4*1000/f2erb(1000);    % p for fc=1kHz and a level of 51dB (at 1kHz filters are symmetrical)
erbdB2F = interp1([0 erbFc fs/2], [min(erbdB) erbdB min(erbdB)], compFq);   % map erbFc to compFq

eL = zeros(size(erbN));
for e = 1:length(erbN)
    erb = f2erb(erbFc(e));
    p51 = 4*erbFc(e)/erb;
    intensity = 0;
    for comp = 1:nPoints
        g = (compFq(comp)-erbFc(e))/erbFc(e);
        if g<0
            p = p51 - 0.35*(p51/p511) * (erbdB2F(comp)-51);
        else
            p = p51;
        end
        g = abs(g);
        w = (1+p*g)*exp(-p*g);
        intensity = intensity+w*compInt(comp);  %intensity per erb
    end
    eL(e) = intensity;
end
results.eLdB = 10*log10(eL./(20e-6)^2); % get dB SPL (20uPa reference)
results.erbN = erbN;
results.fc = erbFc;

%% calculating specific loudness 
dataSL = data_glasberg2002('specLoud','fVec',erbFc);
tQdB = dataSL.tQ;
tQ = 10.^(tQdB./10);
tQdB500 = dataSL.tQ500;
%gdB = dataSL.g;    % low level gain in cochlea amplifier
g = 10.^((tQdB500-tQdB)/10);
a = dataSL.a;    % parameter for linearization around absolute threshold
alpha = dataSL.alpha;    % compressive exponent
c = dataSL.c; % constant to get loudness scale to sone

specLoud = zeros(size(eL));
specLoud1 = c*(2*eL./(eL+tQ)).^1.5 .*((g.* eL + a).^alpha-a.^alpha); % Eq. 6?
specLoud2 = c * ((g .*eL + a).^alpha - a.^alpha); % Eq. 8?
specLoud3 = c*(eL./1.04e6).^0.5; % Eq. 9?
specLoud(eL<tQ) = specLoud1(eL<tQ);
specLoud(eL<=10^10 & eL>tQ) = specLoud2(eL<=10^10 & eL>tQ);
specLoud(eL>10^10) = specLoud3(eL>10^10); % end of Sec. 1.6 in the paper

results.monauralLoudness = sum(specLoud,2) * kv.erbStep;     % integrate over the erbs
results.specLoud = specLoud; % specific loudness