function [SRTs conditions] = joergensen2013_sim(NSpeechsamples,varargin)
%joergensen2013_sim Simulate the experiments shown in figure 2 of Jørgensen, Ewert and Dau (2013)
%   Usage: [SRTs conditions]  = joergensen2013_sim(NSpeechsamples,varargin)
%
%   Output parameters:
%     SRT   : the SRT as a function of the condition used in the experiment
%     conditions : the conditions used for a given experiment
%
%   The flag may be one of:
% 
%     'Jetal2013'           Simulate the conditions with SSN, SAM and ISTS using the
%                           CLUE speech materual, shown in fig2 of
%                           Jørgensen, Ewert and Dau (2013)
% 
%     'JandD2011reverb'     Simulate the conditions with reverberation
%                           shown in fig2 of Jørgensen, Ewert and Dau (2013)
%
%     'JandD2011specsub'    Simulate conditions with spectral subtraction
%                           shown in fig2 of Jørgensen, Ewert and Dau (2013)
%
%     'Kjems2009'           Simulate the data from Kjems et al. (2009)
%                           shown in fig2 of Jørgensen, Ewert and Dau (2013)
% 
%     'FP1990'              Simulate the data from Festen and Plomp
%                           (1990) shown in fig2 of Jørgensen, Ewert and Dau (2013)
%
%   ---------
%
%   Please cite Joergensen et al. (2013) if you use
%   this model.
%
%   See also: joergensen2013, plot_joergensen2013
%
%   References:
%     S. Jørgensen, S. D. Ewert, and T. Dau. A multi-resolution envelope
%     power based model for speech intelligibility. J. Acoust. Soc. Am.,
%     134(1):436--446, 2013.
%     
%     U. Kjems, J. B. Boldt, M. S. Pedersen, T. Lunner, and D. Wang. Role of
%     mask pattern in intelligibility of ideal binary-masked noisy speech. J.
%     Acoust. Soc. Am., 126:1415--1426, 2009.
%     
%     J. Festen and R. Plomp. Effects of fluctuating noise and interfering
%     speech on the speech-reception threshold for impaired and normal
%     hearing. J. Acoust. Soc. Am., 88(4):1725--1736, 1990.
%     
%
%   Url: http://amtoolbox.sourceforge.net/amt-0.10.0/doc/modelstages/joergensen2013_sim.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/>.

%   WARNING: This Model or Modelstage does not currently meet the AMT requirements. 

disp('**********************************************************');
disp('AMT WARNING: incomplete file --- missing: original stimuli');
disp('**********************************************************');
 
definput.flags.type = {'JandD2011specsub','JandD2011reverb','FP1990','Kjems2009','Jetal2013'};

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


% ------- simulate experiment with reverberation
if flags.do_Jetal2013
    %in the absence of the original stimuli - using those from joergensen2011 instead
    data = amt_load('relanoiborra2019', 'single_150_SentArray22kHz_varLength.mat');
    sentenceArray = data.sentenceArray;
    
    amt_disp(['start Jetal2013: ' datestr(now, 'dd-mm-yyyy HH:MM:SS')]);
    sentenceFileLevel = -26;
    %noise_names = {'SSN_CLUE_22kHz.wav','SSN_MOD_CLUE','ISTS_eq'};stimuli
    %not available
    noise_names = 'SSN_CLUE_22kHz.wav';
    speechSPL = 65;
    fixedNoiseLevel = 0;
    SNRs = -27:3:3;
    conditions = {'SSN'  'SAM' 'ISTS'};
    
    IOparameters = [0.61 0.5 8000 0.6]; %
    
    for q = 1:NSpeechsamples
        
        x = sentenceArray{q}';
        fs = 22050;
        x = x*10^((speechSPL-sentenceFileLevel)/20);
        
        % the level of the sentence is set such that the long-term RMS of all
        % sentences are presented at a level of 65 dB SPL.
        
        Ts = 1/fs;
        T = length(x)/fs;
        t = 0:Ts:T;
        t = t(1:end-1);
        N = length(t);
        
        % load the noise files
        for k = 1:length(noise_names)
            
            clear tmp tmp2

            %[tmp, fs] = audioread (noise_names{k});
            [tmp, fs] = amt_load('relanoiborra2019', noise_names);
            %if fs_tmp ~= fs
            %    noise_scaled{k} = resample(tmp,fs,fs_tmp);
            %else
                noise_scaled{k} = tmp;
            %end
            
            Nsegments = floor(length(noise_scaled{k})/N);
            % pick a random segment
            startIdx = randi(Nsegments-2 ,1)*N;
            
            noise_glob{k} = noise_scaled{k}(startIdx:startIdx+N -1);
        end
        
        for n = 1:length(conditions)
            
            for k = 1:length(SNRs)
                noise = noise_glob{n};
                noise = noise/rms(noise)*10^((speechSPL-SNRs(k))/20);
                if size(noise) ~= size(x)
                    noise = noise';
                end
                test = noise + x;
                
                tmp = joergensen2013(test,noise,fs,IOparameters);
                SNRenvs(k,n,q) = tmp.SNRenv;
                Pcorrect(k,n,q) = tmp.P_correct;
                
            end
            
        end
        amt_disp(['sentence nr: ' num2str(q) ' ' datestr(now, 'dd-mm-yyyy HH:MM:SS')]);
        
    end
  
    result.Pcorrect = Pcorrect;
    result.SNRenvs = SNRenvs;
    result.conditions =conditions;
    result.SNRs =SNRs;

    
    
    
    %% Average across speech samples
    Pc_est_mean = mean(result.Pcorrect,3);
    
    
    % ------------------  Estimating changes in SRTs based on the mean Pcorrect -------------
    %   The first column of Pc_est_mean should always be the reference
    selection = 1:length(conditions);
    [dSRT SRTs] = joergensen2011_pctodsrt(Pc_est_mean,result.SNRs,selection);
end

if flags.do_JandD2011specsub
    %in the absence of the original stimuli - using those from joergensen2011 instead
    data = amt_load('relanoiborra2019', 'single_150_SentArray22kHz_varLength.mat');
    sentenceArray = data.sentenceArray;
    
    amt_disp(['start JandD2011specsub: ' datestr(now, 'dd-mm-yyyy HH:MM:SS')]);
    sentenceFileLevel = -26;
    noise_names = 'SSN_CLUE_22kHz.wav';
    speechSPL = 65;
    SNRs = -9:3:9;
    conditions = [0 0.5 1 2 4 8];
    
    IOparameters = [0.61 0.5 8000 0.6]; %
    
    for q = 1:NSpeechsamples
        k = 1;
        x = sentenceArray{q}';
        fs = 22050;
        x = x*10^((speechSPL-sentenceFileLevel)/20);
        
        % the level of the sentence is set such that the long-term RMS of all
        % sentences are presented at a level of 65 dB SPL.
        
        Ts = 1/fs;
        T = length(x)/fs;
        t = 0:Ts:T;
        t = t(1:end-1);
        N = length(t);
        
        % load the noise file
        
        clear tmp
        %[tmp, fs] = audioread (noise_names);
        [tmp, fs] = amt_load('relanoiborra2019', noise_names);
        %if fs_tmp ~= fs
        %    noise_scaled = resample(tmp,fs,fs_tmp);
        %else
            noise_scaled = tmp;
        %end
        
        Nsegments = floor(length(noise_scaled)/N);
        % pick a random segment
        startIdx = randi(Nsegments-2 ,1)*N;
        noise_glob = noise_scaled(startIdx:startIdx+N -1);
        
        
        for n = 1:length(conditions)
            
            for k = 1:length(SNRs)
                noise = noise_glob;
                noise = noise/rms(noise)*10^((speechSPL-SNRs(k))/20);
                if size(noise) ~= size(x)
                    noise = noise';
                end
                test = noise + x;
                
                %                     --------- spec sub -----------------
                W=1024/2; % frame length
                padz=1024/2; %zero padding (pad with padz/2 from the left and padz/2 from the right )
                % % Note that (W+padz) is the final frame window and hence the fft length (it is normally chose as a power of 2)
                SP=0.5; %Shift percentage is 50%
                
                factor = conditions(n);
                ProcMix  = joergensen2011_specsub(test,noise,W,padz,SP,factor);
                
                ProcMixInv = joergensen2011_specsub(x-noise,-noise,W,padz,SP,factor);
%                 Estimating the noise alone using the approach by Hagerman and Olofsson (2004)
                NoiseEst = (ProcMix  - ProcMixInv)/2;
                
                test = ProcMix(1.5*W:N); %cutting off first and last bit
                
                noise = NoiseEst(1.5*W:N);
                %                      -----------------------------
                
                tmp = joergensen2013(test,noise,fs,IOparameters);
                SNRenvs(k,n,q) = tmp.SNRenv;
                Pcorrect(k,n,q) = tmp.P_correct;
                
            end
            
        end
        amt_disp(['sentence nr: ' num2str(q) ' ' datestr(now, 'dd-mm-yyyy HH:MM:SS')]);
        
    end
    
    result.Pcorrect = Pcorrect;
    result.SNRenvs = SNRenvs;
    result.conditions =conditions;
    result.SNRs =SNRs;

    
    
    
    %% Average across speech samples
    Pc_est_mean = mean(result.Pcorrect,3);
    
    
    % ------------------  Estimating changes in SRTs based on the mean Pcorrect -------------
    %   The first column of Pc_est_mean should always be the reference
    selection = 1:length(conditions);
    [dSRT SRTs] = joergensen2011_pctodsrt(Pc_est_mean,result.SNRs,selection);
end

if flags.do_JandD2011reverb
    %in the absence of the original stimuli - using those from joergensen2011 instead
    data = amt_load('relanoiborra2019', 'single_150_SentArray22kHz_varLength.mat');
    sentenceArray = data.sentenceArray;    
    
    amt_disp(['start JandD2011reverb: ' datestr(now, 'dd-mm-yyyy HH:MM:SS')]);
    sentenceFileLevel = -26;
    noise_name = 'SSN_CLUE_22kHz';
    speechSPL = 65;
    SNRs = -9:3:9;
    conditions = [0 0.4 0.7 1.3 2.3];
    
    IOparameters = [0.61 0.5 8000 0.6]; %
    
    for q = 1:NSpeechsamples
        
        x = sentenceArray{q}';
        fs = 22050;
        x = x*10^((speechSPL-sentenceFileLevel)/20);
        
        % the level of the sentence is set such that the long-term RMS of all
        % sentences are presented at a level of 65 dB SPL.
        
        Ts = 1/fs;
        T = length(x)/fs;
        t = 0:Ts:T;
        t = t(1:end-1);
        N = length(t);
        
        % load the noise file
        noise_names = 'SSN_CLUE_22kHz.wav';
        clear tmp
        %[tmp, fs] = audioread (noise_names{k});
        [tmp, fs] = amt_load('relanoiborra2019', noise_names);
%        if fs_tmp ~= fs
%            noise_scaled = resample(tmp,fs,fs_tmp);
%        else
            noise_scaled = tmp;
%        end
        
        Nsegments = floor(length(noise_scaled)/N);
        % pick a random segment
        startIdx = randi(Nsegments-2 ,1)*N;
        noise_glob = noise_scaled(startIdx:startIdx+N -1);
        
        
        for n = 1:length(conditions)
            
            for k = 1:length(SNRs)
                noise = noise_glob;
                noise = noise/rms(noise)*10^((speechSPL-SNRs(k))/20);
                if size(noise) ~= size(x)
                    noise = noise';
                end
                test = noise + x;
                
                %              %     ------------ applying reverberation
                if n>1
                  [tmp,Fs] = sig_joergensen2011(conditions(n));
                                    
                  tmp = resample(tmp,fs,Fs); % downsampling to 22.05 kHz
                  tmp_test = fconv(tmp',test); %
                  tmp_noise = fconv(tmp',noise);
                  test_env = abs(hilbert(tmp_test));
                  
                  [bb, aa] = butter(4, 20*2/fs);
                  test_env = filter(bb,aa,test_env);
                  cut =  0.05;
                  threshold = floor(max(test_env)*cut);
                  %         finding the index in the env vector corresponding to the
                  %         threshold:
                  idx = find(test_env(floor(length(test_env)/4):end) > threshold, 1,'last' );
                  
                  idx_end = idx + floor(length(tmp_test)/4);
                  test_env(idx_end);
                  idx_start = floor(0.047*fs);
                  
                  test = tmp_test(idx_start:idx_end);
                  noise = tmp_noise(idx_start:idx_end);
                end
                % ---------------------------------
                
                tmp = joergensen2013(test,noise,fs,IOparameters);
                SNRenvs(k,n,q) = tmp.SNRenv;
                Pcorrect(k,n,q) = tmp.P_correct;
                
            end
            
        end
        amt_disp(['sentence nr: ' num2str(q) ' ' datestr(now, 'dd-mm-yyyy HH:MM:SS')]);
        
    end
    
    result.Pcorrect = Pcorrect;
    result.SNRenvs = SNRenvs;
    result.conditions =conditions;
    result.SNRs =SNRs;

    
    
    
    %% Average across speech samples
    Pc_est_mean = mean(result.Pcorrect,3);
    
    
    % ------------------  Estimating changes in SRTs based on the mean Pcorrect -------------
    %   The first column of Pc_est_mean should always be the reference
    selection = 1:length(conditions);
    [dSRT SRTs] = joergensen2011_pctodsrt(Pc_est_mean,result.SNRs,selection);
end
%AMT: Stimuli for these last two not available
 if flags.do_FP1990
     amt_disp('amt warning: not enough data to calculate this.');
          SRTs = []; 
     conditions = [];
%     load PlompMimpen_130_SentArray22kHz
%     amt_disp(['start FP1990: ' datestr(now, 'dd-mm-yyyy HH:MM:SS')]);
%     sentenceFileLevel = 10^((-17.93)/20);
%     noise_names = {'MaleS00.wav','MaleF00.wav','P&MRunningSpeech_FemaleTR.wav'};
%     conditions = {'SSN','SMN','RT'};
%     SNRs = [  -21 -18 -12 -9 -6  -3 0 3];
%     noisedBA = 80;
%     fs = 22050;
%     IOparameters = [0.36 0.5 8000 0.6]; %
%     
%     clear noise_scaled
%     for k = 1:length(noise_names)
%         clear tmp tmp2
%         [tmp, fs] = audioread (noise_names{k});
%         tmp = resample(tmp,fs,fs_tmp);
%         tmp = tmp/rms(tmp);
%         noise_scaled{k} = Leq2dBA(tmp,fs,noisedBA); %here, all noises have same level dBA
%         noiseSPL_dBA(k) = LeqdBA(noise_scaled{k},fs);
%     end
%     
%    
%     
%     for q = 1:NSpeechsamples
%         
%         x_unscaled = sentenceArray{q}';
%         
%         
%         Ts = 1/fs;
%         T = length(x_unscaled)/fs;
%         t = 0:Ts:T;
%         t = t(1:end-1);
%         N = length(t);
%         
%         % load the noise files
%         for k = 1:length(noise_names)
%             
%             Nsegments = floor(length(noise_scaled{k})/N);
%             % pick a random segment
%             startIdx = randi(Nsegments-2 ,1)*N;
%             
%             noise_glob{k} = noise_scaled{k}(startIdx:startIdx+N -1);
%         end
%         
%         
%         
%         for n = 1:length(conditions)
%             
%             for k = 1:length(SNRs)
%                 noise = noise_glob{n};
%                 speech_dBA = (noisedBA+SNRs(k));
%                 SPL_sent(k) = 20*log10(rms(Leq2dBA(noise_scaled{1},fs,speech_dBA))); %#check!
%                 SpeechGain = 10^(SPL_sent(k)/20);
%                 x = x_unscaled ./ sentenceFileLevel *SpeechGain;
%                 SpeechSPL_dBA(k) = LeqdBA(x,fs);
%                 
%                 if size(noise) ~= size(x)
%                     noise = noise';
%                 end
%                 test = noise + x;
%                 
%                 
%                 tmp = joergensen2013(test,noise,fs,IOparameters);
%                 SNRenvs(k,n,q) = tmp.SNRenv;
%                 Pcorrect(k,n,q) = tmp.P_correct;
%                 
%             end
%             
%         end
%         amt_disp(['sentence nr: ' num2str(q) ' ' datestr(now, 'dd-mm-yyyy HH:MM:SS')]);
%         
%     end
%     
%     result.Pcorrect = Pcorrect;
%     result.SNRenvs = SNRenvs;
%     result.conditions =conditions;
%     result.SNRs =SNRs;
% 
%     
%     
%     
%     %% Average across speech samples
%     Pc_est_mean = mean(result.Pcorrect,3);
%     
%     
%     % ------------------  Estimating changes in SRTs based on the mean Pcorrect -------------
%     %   The first column of Pc_est_mean should always be the reference
%     selection = 1:length(conditions);
%     [dSRT SRTs] = joergensen2011_pctodsrt(Pc_est_mean,result.SNRs,selection);
%     
 end

 if flags.do_Kjems2009
     amt_disp('amt warning: not enough data to calculate this.');
     SRTs = []; 
     conditions = [];
%      load DANTALE2_144single_SentArray44kHz_varLength
%      amt_disp(['start Kjems2009: ' datestr(now, 'dd-mm-yyyy HH:MM:SS')]);
%             sentenceFileLevel = 10^((-19.79)/20);
%             noise_names = {'ssn_noise_20k','cafe_noise_20k','car_noise_20k','Bottle_noise_20k'};
%             fs = 22050;
%             noisedBA = 65;
%             SNRs = -30:3:10;
%             conditions = {'SSN','Cafe','Car','Bottle'};
%             
%             IOparameters = [0.42 0.5 50 0.9]; %
%             clear noise_scaled
%             
%             for k = 1:length(noise_names)
%                 [tmp, fs] = audioread (noise_names{k});
%                 tmp = resample(tmp,fs,fs_tmp);
%                 tmp = scaletodbspl(tmp,0);
%                 if k == 1
%                     [tmp GaindB]= Leq2dBA(tmp,fs,noisedBA); % the level in dBA is set for the SSN to 65
%                 end
%                 noise_scaled{k} = tmp/rms(tmp).*10^(GaindB/20); %here, all noises have same level in dB, thus different levels dBA!
%                 noiseSPL_dBA(k) = LeqdBA(noise_scaled{k},fs);
%             end
%       
%     for q = 1:NSpeechsamples
%         
%         x_unscaled = sentenceArray{q}';
%         x_unscaled = resample(x_unscaled,fs,44100);
%         
%         Ts = 1/fs;
%         T = length(x_unscaled)/fs;
%         t = 0:Ts:T;
%         t = t(1:end-1);
%         N = length(t);
%         
%         % load the noise files
%         for k = 1:length(noise_names)
%             
%             Nsegments = floor(length(noise_scaled{k})/N);
%             % pick a random segment
%             startIdx = randi(Nsegments-2 ,1)*N;
%             
%             noise_glob{k} = noise_scaled{k}(startIdx:startIdx+N -1);
%         end
%    
%         for n = 1:length(conditions)
%             
%             for k = 1:length(SNRs)
%                 noise = noise_glob{n};
%                 SpeechGain = 10^((GaindB + SNRs(k))/20);
%                  
%                 x = x_unscaled ./ sentenceFileLevel *SpeechGain;
%                 SpeechSPL_dBA(k) = LeqdBA(x,fs);
%                 
%                 if size(noise) ~= size(x)
%                     noise = noise';
%                 end
%                 test = noise + x;
%              
%                 tmp = joergensen2013(test,noise,fs,IOparameters);
%                 SNRenvs(k,n,q) = tmp.SNRenv;
%                 Pcorrect(k,n,q) = tmp.P_correct;
%                 
%             end
%             
%         end
%         amt_disp(['sentence nr: ' num2str(q) ' ' datestr(now, 'dd-mm-yyyy HH:MM:SS')]);
%         
%     end
%     
%     result.Pcorrect = Pcorrect;
%     result.SNRenvs = SNRenvs;
%     result.conditions =conditions;
%     result.SNRs =SNRs;
% 
%         
%     
%     %% Average across speech samples
%     Pc_est_mean = mean(result.Pcorrect,3);
%     
%     % ------------------  Estimating changes in SRTs based on the mean Pcorrect -------------
%     %   The first column of Pc_est_mean should always be the reference
%     selection = 1:length(conditions);
%     [dSRT SRTs] = joergensen2011_pctodsrt(Pc_est_mean,result.SNRs,selection);
%     
 end

function [y]=fconv(h,x)
%FCONV Fast Convolution
%   [y] = FCONV(h, x) convolves x and h.  The output of this
%         function is scaled.
%
%      x = input vector
%      h = input vector
%
%      See also CONV
%
%Version 2.0
%2003-2004, Stephen G. McGovern

Ly=length(x)+length(h)-1;   %
Ly2=pow2(nextpow2(Ly));     % Find smallest power of 2 that is > Ly
m= max(abs(x));
X=fft(x, Ly2);		    % Fast Fourier transform
H=fft(h, Ly2);		    % Fast Fourier transform
Y=X.*H;
y=real(ifft(Y, Ly2));       % Inverse fast Fourier transform
y=y(1:1:Ly);                % Take just the first N elements

m=m/max(abs(y));
y=m*y;


function [ret Signalscale]=Leq2dBA(sig,Fs,dBA)


InSig=sig/rms(sig).*10^(dBA/20);           % scale to dB SPL
x2=InSig;
[B,A] = adsgn(Fs);                    %
r=filter(B,A,InSig);


Laeq=20*log10(rms(r));       %dB(A)
Lleq=20*log10(rms(x2));      %dB
cor=Lleq-Laeq;
Signalscale=dBA+cor;
ret=InSig/rms(InSig).*10^(Signalscale/20);

function ret=LeqdBA(sig,Fs)

InSig=sig;           % scale to dB SPL

[B,A] = adsgn(Fs);
r=filter(B,A,InSig);
ret=20*log10(rms(r));       %dB(A)


function [B,A] = adsgn(Fs)
% ADSGN  Design of a A-weighting filter.
%    [B,A] = ADSGN(Fs) designs a digital A-weighting filter for
%    sampling frequency Fs. Usage: Y = FILTER(B,A,X).
%    Warning: Fs should normally be higher than 20 kHz. For example,
%    Fs = 48000 yields a class 1-compliant filter.
%
%    Requires the Signal Processing Toolbox.
%
%    See also ASPEC, CDSGN, CSPEC.

% Author: Christophe Couvreur, Faculte Polytechnique de Mons (Belgium)
%         couvreur@thor.fpms.ac.be
% Last modification: Aug. 20, 1997, 10:00am.

% References:
%    [1] IEC/CD 1672: Electroacoustics-Sound Level Meters, Nov. 1996.


% Definition of analog A-weighting filter according to IEC/CD 1672.
f1 = 20.598997;
f2 = 107.65265;
f3 = 737.86223;
f4 = 12194.217;
A1000 = 1.9997;
pi = 3.14159265358979;
NUMs = [ (2*pi*f4)^2*(10^(A1000/20)) 0 0 0 0 ];
DENs = conv([1 +4*pi*f4 (2*pi*f4)^2],[1 +4*pi*f1 (2*pi*f1)^2]);
DENs = conv(conv(DENs,[1 2*pi*f3]),[1 2*pi*f2]);

% Use the bilinear transformation to get the digital filter.
[B,A] = bilinear(NUMs,DENs,Fs);

