%% ABR analysis script: this script takes the raw recordings and extracts the results
% raw recordings are available from last author upon request

%% settings
fs=19200;                   % sampling rate after 20-fold resampling
bma200_input_gain=10000;    % value gained during calibration
abr_scale_factor=1/(7e-5*bma200_input_gain);    % in mV
nboot=500;                  % number of bootstraps for statistics

% set if you want to save the results  and bootstrap options (1=on, 0=off)
saveon=0;
booton=1;

% set abr analysis window
ana_wind_start=0.001;       % offset of rms calculation window
ana_wind_dur=0.007;         % duration of rms calculation window
ana_wind_starts=round(ana_wind_start*fs);
ana_wind_stops=round((ana_wind_start+ana_wind_dur)*fs);

% define (internal) bat ID, recording date and recording number for analyses
dates = {'20170207','20170210','20170817','20191204',...
    '20170213','20170223','20170818','20191204',...
    '20170217','20170220','20170821','20191205',...
    '20170203','20170209','20170817','20191204',...
    '20170208','20170216','20170818','20191203',...
    '20170306','20170313','20170821','20191205'};

reco = {'3','3','3','2','3','3','8','2','3','5','3','2',...
        '5','3','7','2','3','3','3','3','3','7','4','2'};

animalname = {'3','3','3','3','5','5','5','5','8','8','8','8',...
    '2','2','2','2','4','4','4','4','7','7','7','7'};

%% data extraction and bootstrap analysis
for i=1:24
    % choose current bat, date, record number
    batid=char(animalname(i));
    date=dates(i);
    recn=reco(i);
    
    % location of folder with the raw recordings
    rootpath=['D:\Documents\Manuskripte\Vocal_development\ABR_recordings\B' batid '\' char(date) '\'];
    
    % circumvent naming error of B2 files
    if batid=='2'
        recname=['abr_tonepips_ed_b' batid '_m' char(recn) '_'];
    else
        recname=['abr_tonepips_vpl_b' batid '_m' char(recn) '_'];
    end
    
    % access the folder with the recordings
    d=dir(rootpath);
    filenames=char(d.name);     % read in the names of the files in this folder
    recnum=0;                   % set recnum to zero
    for fn=3:size(filenames,1)  % go through all files in this folder
        filename=deblank(filenames(fn,:));
        if 1-isempty(strfind(filename,recname))
            % extract info from filename
            li=strfind(filename,'_l');
            fi=strfind(filename,'_fc');
            xi=strfind(filename,'_');
            xi=xi(min(find(xi>fi)));
            if xi>fi+4          % if there is frequency information in the filename
                xi=max(xi);
                recnum=recnum+1;
                lstring=filename(li+2:fi-1);
                current_level=str2num(lstring);
                fstring=filename(fi+3:xi-1);
                current_freq=str2num(fstring);
                load([rootpath filename]); % load the file to be analysed
            end
            
            mrxin=squeeze(mean(rxin,2)).*abr_scale_factor;
            levels(recnum)=current_level;
            freqs(recnum)=current_freq/1e3;
            rec(recnum,:)=mrxin;
            
            %%% BOOTSTRAPPING ANALYSIS STARTS HERE
            if booton==1
                data=rxin(ana_wind_starts:ana_wind_stops,:);
                datalong=reshape(data,numel(data),1);
                
                % do a bootstrapping analysis
                samplesize=size(data,1);
                indices=1:(length(datalong)-samplesize);
                samplenumber=size(data,2);
                originalsignal=mean(data,2);
                originalrms=std(originalsignal);
                
                for bn=1:nboot
                    for rep=1:size(data,2)
                        bsdata(:,rep)=circshift(data(:,rep),ceil(rand*size(data,1)));
                    end
                    bootrms(bn)=std(mean(bsdata,2));
                end
                
                confidence=numel(find(originalrms>bootrms))/nboot;
                confidences(recnum)=confidence;
            end     %%% BOOTSTRAP ANALYSIS ENDS HERE
        end
    end
    
    % organise data for bootstrap plots
    all_levels=sort(unique(levels));
    all_freqs=sort(unique(freqs));
    out=zeros(length(all_levels),length(all_freqs));
    for sn=1:length(levels)
        cl=levels(sn);
        cf=freqs(sn);
        li=find(all_levels==cl);
        fi=find(all_freqs==cf);
        out(li,fi)=std(rec(sn,ana_wind_starts:ana_wind_stops)); % calculate rms from sample 31 because of rec onset irregularities
    end
    in=[freqs' levels' confidences'];
    sort_in=sortrows(in,[1 2]);
    ufreqs=unique(freqs);
    ulevels=unique(levels);
    
    for fn=1:length(ufreqs)
        cfreq=ufreqs(fn);
        fconf=sort_in(find(sort_in(:,1)==cfreq),3);
        threshold_ind=max(find(fconf<0.95))+1;
        if threshold_ind<=length(ulevels) & 1-isempty(threshold_ind)
            bootstrap_thresholds(fn)=ulevels(threshold_ind);
        else
            bootstrap_thresholds(fn)=nan;
        end
    end
    
    % save results in structures
    resultTHRES(i)=struct('uniqueFrequencies',ufreqs','bootstrapThresholds',bootstrap_thresholds');
    allCONFI(i)=struct('frequencies',freqs,'levels',levels,'confidences',confidences);
    heatRES(i)=struct('all_freqs',all_freqs,'all_levels',all_levels,'out',out);
end

%% saving of the variables
if saveon==1
    bats = strcat('b', animalname);
    save('abr_deafening_results.mat','dates','bats','resultTHRES','allCONFI');
end