function [totalLoudness,specLoudness,Estim,Ethq] = calcTotalLoudness(sig)
% [totalLoudness,specLoudness,Estim] = calcTotalLoudness(sig)
% calculation of excirtation pattern, specific and total loudness
% 10.12.2002 | Marc Schönwiesner, Institut for Zoology, University of Leipzig

%%%%%%%%%%%%%%%%%%%%%%
%%%% define model %%%%
%%%%%%%%%%%%%%%%%%%%%%
disp('Defining model...');
clear model;
model.fs=44100;
model.cochlea.fb.file='gt64';% 64-channel gammatone filterbank
model.cochlea.asymmcomp=1;   % Activate asymmetric compensation filterbank
model.haircell.rcf.r='half'; % Half-wave rectification
model.haircell.rcf.c=0.7;    % Compression by c^0.7
model.haircell.rcf.f='1kHz'; % Low-pass filter
model.neural.function='mean';% Mean within channel

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%% calculate excitation pattern at threshold %%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% This is done by generating pure tones at the filter
% CentFreqs with an intensity that equals the MAF at
% these frequencies.
% output: Ethq -> Exitation at MAF 
% The calculation need only be done once for a specific model;
% the resulting Ethq could be saved as .mat file and reused.

% get MAF
disp('Interpolating MAF...');
load gt64;
Npnts = size(CentFreqs,2);
[CrctLinPwr, frqNpnts, MAF] = OutMidCrct('MAF',0,0,0);
MAFinterp = interp1(frqNpnts,MAF,CentFreqs,'linear','extrap');

% make tones
disp('Generating signals...');
for i=1:Npnts
    tones(i,:)=Pascalize(tone(CentFreqs(i),5000),MAFinterp(i));
end;

% get exitation pattern (same model as above)
disp('Calculating Exitation (Ethq). This might take several minutes (64 steps)...');
A=[];
for i=1:Npnts,
  fprintf('\t%d',i);
  A=[A; AudMod(tones(i,:),model)];
end
figure;plot(CentFreqs,A);axis([0,20000,0,10]);xlabel('Frequency [kHz]');ylabel('Excitation');title('Excitation pattern of 64 tones at minimum audiable field levels');

% combine exitation for all 64 tones into one exitation pattern 
Ethq = zeros(1,Npnts);
for i=1:Npnts
    Ethq(i) = A(i,i); % 
end
figure;plot(CentFreqs,Ethq);axis([0,20000,0,10]);xlabel('Frequency [kHz]');ylabel('Excitation');title('Excitation pattern of 64 tones at MAF level - collapsed');

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%% calculate exitation pattern of stimulus %%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

disp('Calculating Exitation (Estim)...');
Estim = [];
Estim = AudMod(sig,model);
figure;plot(CentFreqs,Estim);xlabel('Frequency [kHz]');ylabel('Excitation');title('Excitation pattern of signal');

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%% calculate specific loudness %%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

% compress excitation pattern according to Moore & Glasberg (1996):
%                       a
% N'     = C * (E      )     -> Loudness produced by Stimulus
%   Stim         |ELCC|

%                       a
% N'     = C * (E Thq  )     -> Loudness produced by internal Noise
%   intN         |ELCC|

% and 
% N' = N'     - N'
%        intN     Stim

% so
%                     a            a
% N' = C * ( (E      )  - (E Thq  )    -> specific Loudness
%              |ELCC|       |ELCC|  

% with N'       ... specific Loudness
%      N'       ... specific Loudness of the Stimulus
%        Stim   
%      N'       ... specific Loudness of internal Noise
%        intN  
%      E        ... equal level loudness contour corrected excitation
%       |ELCC|  
%      E Thq    ... ELC-corrected excitation at hearing threshold
%       |ELCC| 
%      C        ... parameter (shift of function along loudness axis)
%      a        ... parameter (a<1 -> nonlinear compression)


% set constants - these parameters are fitted manually to approx. conform
% to the sone scale definition:
% 1kHz 40dB pure tone -> 1 sone, 50dB  -> 2 sone, 60dB -> 4 sone, 70 dB -> 8 sone etc.
% these values yield:
% 40dB -> 1.0 sone, 50dB  -> 2.1 sone, 60dB -> 4.6 sone, 70 dB -> 8.1 sone
C = 0.1;
a = 0.005;
disp('Calculating specific loudness (specLoudness)...');
specLoudness = C .* (Estim.^a - Ethq.^a);
specLoudness(specLoudness<0)=0; % no specific loudness below zero - rule by Moore & Glasberg
figure;plot(CentFreqs,specLoudness);xlabel('Frequency [kHz]');ylabel('Specific Loudness');title('Specific Loudness of signal');

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%% calculate total loudness %%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% integrate specific loudness
% N = sum(N'(ferb) d(ferb))
disp('Calculating ERBwidth...');
EarQ = 1/0.107939;
minBW = 24.7;
order = 1;
b = 1.019;
ERBwidth = ((CentFreqs/EarQ).^order + minBW^order).^(1/order);
disp('Calculating total loudness...');
totalLoudness = sum(specLoudness.*ERBwidth);
disp('Done.');