function [InstantaneousSpecificLoudnessLeft, InstantaneousSpecificLoudnessRight] = moore2016_monauralinstspecloudness(signal, Fs, dBMax)
%MOORE2016_MONAURALINSTSPECLOUDNESS calculates instantaneous specific loudness over time
%
%   Input parameters:
%     signal : input signal
%     fs : sampling frequency
%     dbMax : the RMS SPL of a sinusoid with a peak amplitude of 1.
%
%   Output parameters:
%     InstantaneousSpecificLoudnessLeft  : instantaneous specific loudness left ear
%     InstantaneousSpecificLoudnessRight : instantaneous specific loudness right ear
%
%   (left, right) signal and sampling rate Fs must be passed. Fs must be 32 kHz, 
%   reading of the wav file and resampling must be done before calling this function. 
%   The signal must contain two columns.
%
%   version for TVL 2016 based on ANSI S3.4-2007 and Moore & Glasberg (2007)
%
%   Url: http://amtoolbox.org/amt-1.1.0/doc/modelstages/moore2016_monauralinstspecloudness.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/>.

npts               = Fs / 1000 * 64;                     % points for FFT, 2048
nSegmentDuration   = 1;                                  % duration of segment in ms
nSamplesPerSegment = Fs / 1000 * nSegmentDuration;       % 32
nSegmentsInSignal  = floor( (length(signal) - npts) / Fs * 1000 / nSegmentDuration );

% Hann windows for 6 FFTs; 1st column 64 ms, 6th column 2 ms
wHann = zeros(npts,6);
for i = 1:6
    wHann(:,i) = [ zeros( (1 - 1/2^(i-1) ) / 2 * npts, 1 ); hann( npts / 2^(i-1) ); zeros( (1 - 1/2^(i-1) ) / 2 * npts, 1 ) ];
end

% indices which shall be used from the ffts, 20-80 Hz from the first...
vLimitingF = [20 80 500 1250 2540 4050 15000]; 
vLimitingIndices = zeros(1,7);
for i=1:7
    vLimitingIndices(i) = floor( vLimitingF(i) / (Fs/npts) ) + 1;
end

InstantaneousSpecificLoudnessLeft  = zeros(nSegmentsInSignal + 1, length(1.75:0.25:39) );
InstantaneousSpecificLoudnessRight = zeros(nSegmentsInSignal + 1, length(1.75:0.25:39) );

%% loop: window segment, apply stationary loudness model, move window
%% by 1 ms

% tic
for iSegment = 0:nSegmentsInSignal
%     
%     if ( mod(iSegment,50) == 0 )
%         disp([num2str(iSegment) ' ms of ' num2str(nSegmentsInSignal) ' ms analyzed'])
%     end
    
    [fLeftRelevant, LLeftRelevant, fRightRelevant, LRightRelevant] = moore2016_spectrum(signal( (iSegment*nSamplesPerSegment+1) : (iSegment*nSamplesPerSegment+npts), : ), Fs, dBMax, wHann, vLimitingIndices);
    
    if ( isempty( LLeftRelevant ) )
        SpecificLoudnessLeft = zeros(1,150);
    else
        ExcitationLevelsLeft = moore2016_excitationpattern(fLeftRelevant', LLeftRelevant');
        SpecificLoudnessLeft = moore2016_specloudnessbinaural( ExcitationLevelsLeft );
    end
    if ( isempty( LRightRelevant ) )
        SpecificLoudnessRight = zeros(1,150);
    else
        ExcitationLevelsRight = moore2016_excitationpattern(fRightRelevant', LRightRelevant');
        SpecificLoudnessRight = moore2016_specloudnessbinaural( ExcitationLevelsRight );
    end
%     [ ThisSpecificLoudness, ThisLoudness ] = MonauralSpecificLoudness2BinauralSpecificLoudness025( SpecificLoudnessLeft, SpecificLoudnessRight );
    InstantaneousSpecificLoudnessLeft(iSegment+1,:) = SpecificLoudnessLeft;
    InstantaneousSpecificLoudnessRight(iSegment+1,:) = SpecificLoudnessRight;
%     InstantaneousLoudness(iSegment+1)         = ThisLoudness;
end
% toc
end 

function out = moore2016_specloudnessbinaural( ExcitationLevels )

% version for TVL 2016 based on ANSI S3.4-2007 and Moore & Glasberg (2007)
% Calculate the specific loudness out of excitation patterns at 0.25-ERB
% steps, using the constant C of the model that accounts for binaural
% loudness

% ERB-scale
ERBc = 1.75:0.25:39;
fc = erb2fc(ERBc);

% Parameters
C = 0.06267;
C = 0.063026;
C = 0.0631;
G = get_G(fc);
Alpha = get_Alpha(fc);
A = get_A(fc);

% Excitation and threshold in linear power units
Excitation = 10 .^ ( ExcitationLevels ./ 10 );
Threshold = excitationthreshold(fc);
Threshold = 10 .^ ( Threshold ./ 10 );

% 3 ranges - calculate all for each case (computationally cheap)
N1 = C .* ( ( G .* Excitation + A ) .^ Alpha - A .^Alpha );
N2 = C .* ( 2 .* Excitation ./ ( Excitation + Threshold ) ) .^ 1.5 .* ( ( G .* Excitation + A ) .^ Alpha - A .^Alpha );
N3 = C .* ( Excitation ./ 1.0707 ) .^ 0.2;

out = ( Excitation > Threshold & Excitation < 10^10 ) .* N1 + ( Excitation <= Threshold ) .* N2 + (Excitation >= 10^10 ) .* N3;

% figure;
% hold on;
% semilogy(ExcitationLevels, out);
% xlim([0 110]);
% ylim([0.005 10]);
end

%end

function out = get_Alpha(f)

TableG2Alpha = [ -25.0   -20.0   -15.0   -10.0   -5.0    0.0;
                 0.26692 0.25016 0.23679 0.22228 0.21055 0.20000];

G = get_G(f);
G = 10 .* log(G) ./ log(10);

out = interpolation(TableG2Alpha(1,:), TableG2Alpha(2,:), G,'linear');

end

function out = get_A(f)

TableG2A = [ -24.54531 -23.78397 -22.78169 -21.76854 -20.74442 -19.78305 -18.90431 -18.01605 -17.11816 -16.21055 -15.32375 -14.59341 ...
             -13.91727 -13.29726 -12.73537 -12.23364 -11.75255 -11.23866 -10.75136 -10.29164  -9.86051  -9.45902  -9.08823  -8.72191 ...
              -8.35715  -8.01199  -7.68715  -7.38338  -7.10145  -6.84213  -6.60623  -6.39458  -6.14589  -5.89392  -5.65071  -5.41661 ...
              -5.19198  -4.97718  -4.77258  -4.57857  -4.39555  -4.20148  -4.00538  -3.81442  -3.62882  -3.27454  -3.10633  -2.94438 ...
              -2.78894  -2.64027  -2.50042  -2.37015  -2.24820  -2.13487  -2.03046  -1.93531  -1.84973  -1.77405  -1.70863  -1.65382 ...
              -1.60997  -1.57745  -1.51786  -1.44522  -1.37466  -1.30624  -1.24006  -1.17621  -1.11478  -1.05587  -0.99956  -0.94596 ...
              -0.89518  -0.84731  -0.80246  -0.75663  -0.69834  -0.64029  -0.58251  -0.52501  -0.46783  -0.41099  -0.35451  -0.29842 ...
              -0.24274  -0.18750  -0.13273  -0.07845  -0.02470   0.00000;
               8.85200   8.63150   8.35840   8.10120   7.85850   7.65258   7.49124   7.32853   7.16458   6.99954   6.83957   6.70559 ...
               6.58115   6.46869   6.36813   6.27970   6.19583   6.10719   6.02404   5.94640   5.87876   5.82490   5.77551   5.72719 ...
               5.67939   5.63443   5.59236   5.55322   5.51708   5.48399   5.45402   5.42723   5.39588   5.36425   5.33386   5.30472 ...
               5.27688   5.25064   5.22800   5.20667   5.18660   5.16539   5.14401   5.12326   5.10314   5.06490   5.04681   5.02944 ...
               5.01280   4.99693   4.98203   4.96818   4.95524   4.94324   4.93219   4.92215   4.91312   4.90515   4.89827   4.89251 ...
               4.88790   4.88449   4.87824   4.87063   4.86324   4.85609   4.84918   4.84252   4.83611   4.82998   4.82412   4.81854 ...
               4.81327   4.80830   4.80365   4.79890   4.79286   4.78685   4.78088   4.77494   4.76904   4.76318   4.75736   4.75159 ...
               4.74586   4.74019   4.73456   4.72899   4.72349   4.72096 ];

G = get_G(f);
G = 10 .* log(G) ./ log(10);

out = interpolation(TableG2A(1,:), TableG2A(2,:), G);

end

function out = get_G(f)
%addpath(fullfile (pwd, 'functions'));
% low level gain of the cochlear amplifier relative to the gain at 500 Hz

LinearThreshold = 10.^( excitationthreshold(f) ./ 10 );
out = 10^(3.63/10) ./ LinearThreshold;
end

function out = excitationthreshold( f )
% version for TVL 2016 based on ANSI S3.4-2007 and Moore & Glasberg (2007)         
Above500 = 3.63;
Threshold = [ 50    63    80    100   125   160   200  250  315  400  500  630  750  800  1000;
              28.18 23.90 19.20 15.68 12.67 10.09 8.08 6.30 5.30 4.50 3.63 3.63 3.63 3.63 3.63 ];

Below500 = interpolation(Threshold(1,:), Threshold(2,:), f, 'linear');
          
out = ( f < 500 ) .* Below500 + ( f >= 500 ) .* Above500;
end