%% Extract and plot eABR data
% written by Charlotte Navntoft (chamnav@elektro.dtu.dk)
% date: 30.04.2018

% This script provides a simple code to extract and plot eABR data recorded 
% with a TDT system. The trigger from the stimulation platform and the 
% recorded eABR are stored as variables. The trigger is used to find and 
% define the start of the epoch window.

% Further signal processing, including removal of the CI artifact, 
% removal of breathing trials, baseline correction, filtering, and 
% averaging techniqies, can be applied afterwards.


% The data need to be saved in a structure called "data"
% with the following fields:
% data.SR: one scalar representing the sampling rate set by the TDT.
% data.ABR: one column vector with the recording of the ABR probe
% data.Trigger: One column vector with a probe connected to the Trigger 
%
% It is assumed that the sampling rate is the same for the Trigger and the
% ABR, and length(data.ABR) == length(data.Trigger)

% An example of an eABR recording is provided: filename 180110CN_Block-17.mat
% The stimuli is a biphasic pulse with 50us/phase and a 10us interphase gap
% presented at a current level of 300 uA with a repitition rate of 23.3 Hz. 
% The stimuli is presented in monopolar configuration on e1 (second most 
% apical electrode) and with 400 repetitions. 

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

% PARAMETERS

dur_epoch_ms = 40;  % set duration of epochs depending on the pulse presentation rate
offset_dur_ms = -2; % set epoch offset relative to the trigger
trigger_threshold = 0.5; % trigger threshold to detect the trigger
Epoch_Buffer_Size = 500; % Buffer size for the epoch. Increase if expecting more epoch



%% Load the data

data_file = uigetfile;
load(data_file);   % load the data

%% Extract eABR epochs


size_epoch = floor(dur_epoch_ms*data.SR/1000); %Nb of sample for 1 epochs
epochs = zeros(Epoch_Buffer_Size,size_epoch); % Create a 2D buffer
offset_sample = floor(offset_dur_ms*data.SR/1000); % Nb of Sample before trigger

ii=1; % index that will scan all the sample
i_epoch = 0; % Index for the epoch

while ii<(length(data.Trigger)-size_epoch)
    if data.Trigger(ii)>trigger_threshold % If the trigger is higger than a threshold 
                                          % and Epoch is detected  
        i_epoch = i_epoch+1; % increment the epoch counter
        epochs(i_epoch,:) = data.ABR((1:size_epoch)+ii+offset_sample); % fill the epoch buffer
        ii = ii+size_epoch; % increase tne sample counter to the end of the epoch
    else
        ii=ii+1; % if no trigger is found go to the next sample
    end
end


disp([' Number of Epoch detected : ' num2str(i_epoch)]); % print number of epochs detected 


epochs_mean = mean(epochs(1:i_epoch,:));

%% Plot epochs
time_axis = 1e3*(1:size_epoch)/data.SR+offset_dur_ms; % in ms

figure
    plot(time_axis,1e3*epochs_mean);
    title(['data file:  ' data_file ])
    xlabel('Time relative to trigger (ms)')
    ylabel('Amplitude (mV)')
    xlim([-1 9])
    ylim([-1 1])
    