%%EYE-BLINK ARTIFACT REJECTION
%%Ellen Lau - November 2006

%%This function removes eye-blink artifacts from MEG data by thresholding
%%amplitude of RMS of 30+ selected eye-blink sensitive channels. User
%%chooses threshold based on testing with particular subject's data.

%%usage: c162.clean = blinkreject(c162.rawEp,2E-13);

%%currently the channel file is hardcoded in. if it turns out that the
%%channels are subject-specific, this will be input; otherwise, seems
%%easier to just alter this line as needed


function [cleandata] = blinkreject(data,thresh)


cleandata = [];
p = 0;
for x = 1:size(data,3)                %% epoch loop
    maxdiff(x) = 0;
    signal = rms(data(:,:,x),'R0593 ALL EYEBLINK CHANNELS.txt');

    for y = 51:(size(signal)-50)        %% time loop: start at point 51, end at end-50, leave room for window
        ampa = mean(signal(y-50:y));    %% mean amplitude A, of 50 points before current point
        ampb = mean(signal(y:y+50));    %% mean amplitude B, of 50 points following current point
        diff = abs(ampb-ampa);          %% absolute value difference between the two
        if (diff > maxdiff(x)) maxdiff(x) = diff;     %%find the maximum difference value in the epoch
        end;
    end
    if (maxdiff(x) < thresh) cleandata(:,:,p+1) = data(:,:,x);   %% create new array with only epochs below max value threshold
    end;
    [m,n,p] = size(cleandata);               %% update position in new clean array for next epoch
end

bar(maxdiff)                      %% plot difference values to assess algorithm