clear all;
close all;
clc;

epochMarkers = {'S202','S203'};     % markers in the analysis
currentEpoch = [-0.2 0.6];          % epoch size
baselinePoints = 100;               % number of data points for the baseline
artifactCriteria = 100;             % maximum artifact difference criteria

% load the data
EEG = pop_loadbv;

% ensure you have channel locations
EEG=pop_chanedit(EEG,'lookup','Standard-10-20-Cap81.ced');

% reference the data using TP9 and TP10, channels 63 and 60
referenceChannel = (EEG.data(63,:) + EEG.data(60,:))/2;
EEG.data = EEG.data - referenceChannel;

% filter the data
[b,a] = butter(2,[0.1 30]/(EEG.srate/2)); % define Butterworth filter
for counter = 1:size(EEG.data,1) % filter data by channel and sample
    EEG.data(counter,:) = filtfilt(b,a,double(EEG.data(counter,:)));  
end
% do the notch filter
Qfactor = 35; 
wo = 60/(EEG.srate/2);  bw = wo/Qfactor;
[b,a] = iirnotch(wo,bw);
for counter = 1:size(EEG.data,1) % filter data by channel and sample
    EEG.data(counter,:) = filtfilt(b,a,double(EEG.data(counter,:)));  
end

% epoch the data
EEG = pop_epoch(EEG,epochMarkers,currentEpoch);

% baseline correction
for trialCounter = 1:size(EEG.data,3)
    for channelCounter = 1:size(EEG.data,1)
        EEG.data(channelCounter,:,trialCounter) = EEG.data(channelCounter,:,trialCounter) - mean(EEG.data(channelCounter,1:baselinePoints,trialCounter));
    end
end   

% mark artifacts using a difference criteria
artifactDifference = zeros(size(EEG.data,1),size(EEG.data,3));
artifactSize = zeros(size(EEG.data,1),size(EEG.data,3));
for channelCounter = 1:size(EEG.data,1)
    for segmentCounter = 1:size(EEG.data,3)
        voltage_max = [];
        voltage_max = max(EEG.data(channelCounter,:,segmentCounter));
        voltage_min = [];
        voltage_min = min(EEG.data(channelCounter,:,segmentCounter));
        difference = voltage_max - voltage_min;
        if difference > artifactCriteria
            artifactDifference(channelCounter,segmentCounter) = 1;
        end
        artifactSize(channelCounter,segmentCounter) = difference; 
    end
end
artifactTrials = mean(artifactDifference,1);

% sort the data into a matrix for ERP analysis
for conditionCounter = 1:size(epochMarkers,2)
    tempData = [];
    tempCounter = 1;
    for trialCounter = 1:size(EEG.data,3)
        if strcmp(EEG.event(trialCounter).type,epochMarkers(conditionCounter)) & artifactTrials(trialCounter) == 0
            tempData(:,:,tempCounter) = EEG.data(:,:,trialCounter);
            tempCounter = tempCounter + 1;
        end
    end
    [fftResults(:,:,conditionCounter),frequencies(:,:,conditionCounter)] = doFFT(tempData,EEG.srate);
end

% plot some stuff
subplot(1,2,1);
bar(fftResults(50,1:28,1));
ylim([0 50]);
title('Oddball');

subplot(1,2,2);
bar(fftResults(50,1:28,2));
ylim([0 50]);
title('Control');