%% Template for processing of startle data collected in the Threat Proability Task 

%Assumptions:
%One channel of EMG data sampled from the orbicularis oculi muscle under the eye has been loaded as a data structure into EEGlab
%Event markers for startle events are included in the EEG file

%% High pass filter
        [HP28_1,HP28_2] = butter(2,28/(EEG.srate/2),'high');    % generate filter coefficient for 4th order 28Hz High pass Butterworth filter
        EEG.data = filtfilt(HP28_1,HP28_2,EEG.data);  %apply that filter to ORB channel; filter order is divided by 2 because we use the filtfilt function to avoid phase shifts in the data

%% Rectify 
        EEG.data = abs(EEG.data); %the 'abs' function rectifies the data
        
%% Low pass filter 
        [LP30_1,LP30_2] = butter(2,30/(EEG.srate/2),'low');    % generate filter coefficient for 4th order 30Hz Low pass Butterworth filter
        EEG.data = filtfilt(LP30_1,LP30_2,EEG.data);  %%apply that filter to ORB channel; filter order is divided by 2 because we use the filtfilt function to avoid phase shifts in the data

%% Epoch and baseline correct
        EventCodes = num2cell(unique([EEG.event.type])); %Extract event codes for epoching            
        EEG = pop_epoch( EEG, EventCodes, [-0.05 0.25], 'epochinfo', 'yes');  %Epoch file for events that indicate startle probes
        EEG = pop_rmbase( EEG, [-50   0]); %Baseline correct

%% Mark and reject trials with excessive artificat 
        Artifact = squeeze(max(EEG.data(1,1:find(EEG.times==0),:)) - min( EEG.data(1,1:find(EEG.times==0),:)) > 40); %Identify trials with greater than 40 microvolt deflections in baseline
        EEGNA = pop_select(EEG,'notrial', find(Artifact)); %Reject these trials 
        
%% Score startle peaks 
        STLResponse = squeeze(max(EEG.data(1,find(EEG.times==20):find(EEG.times==100),:))); %Identify peak response of startle trials in scoring window
        STLResponse = [[EEG.event.type]' STLResponse]; %Store peak values by event type 