mpib_memoreye/gaze_shift/detectMS_EK03_jld.m

function [angles, sizes]=detectMS_EK03(actgaze, smoothing)
% BS and JLD adapted from FieldTrip function ft_detect_movement.m to work without FT
% Inputs: actgaze: the x/y gaze time series of the current trial (2 x nsamples)
% Outputs: vectors of angles (directions) and sizes, respectively, of detected microsaccades. 
% Can be modified to output other things like velocity or time of the MS

%% Settings (Engbert & Kliegl 2003)
kernel = [1 1 0 -1 -1].*(500/6); % 500 is Sampling Rate
velthres=8;%6;
mindur=12;%3;

% If Smoothing

if smoothing 
    sampling_rate=500;
    smooth_window_size = 7 / 1000 * sampling_rate; % in samples
    actgaze = smoothdata(actgaze', 'gaussian', smooth_window_size)'; %Original data smooth
end


% mean-center the trial
actgaze(1,:)=actgaze(1,:)-nanmean(actgaze(1,:));
actgaze(2,:)=actgaze(2,:)-nanmean(actgaze(2,:));

% padding, to prevent edge effects usign conv
pad=mindur;
padleft=repmat(nanmean(actgaze(:,1:pad),2),1,pad);
padright=repmat(nanmean(actgaze(:,end-pad+1:end),2),1,pad);
padgaze=[padleft actgaze padright];

% velocity time series
vel=convn(padgaze, kernel, 'same');
vel=vel(:,1+pad:end-pad); % remove pad 

% velocity exceedance (saccades)
medianstd = sqrt( median(vel.^2,2) - (median(vel,2)).^2 );
radius = velthres*medianstd;
test = sum((vel./radius(:,ones(1,length(actgaze)))).^2,1);
sacsmp = find(test>1); % microsaccade's indexing

%% determine microsaccades per trial
% first find eye movements of n-consecutive time points
j = find(diff(sacsmp)==1);
j1 = [j; j+1];
com = intersect(j,j+1);
cut = ~ismember(j1,com);
sacidx = reshape(j1(cut),2,[]);

sizes=[];
angles=[];
for k=1:size(sacidx,2)
    duration = sacidx(1,k):sacidx(2,k);
    if size(duration,2) >= mindur
        begtrl = sacsmp(duration(1,1));   % begin of MS
        endtrl = sacsmp(duration(1,end)); % end of MS
        [peakvel, smptrl] = max(sqrt(sum(vel(:,begtrl:endtrl).^2,1)));
        veltrl = sacsmp(duration(1,smptrl)); % peak velocity microsaccade sample -> important for spike conversion
        % what follows is only some ugly code by BS to average a few (here 3)
        % samples before and after the MS, hoping it makes things less noisy
        afew=3-1;
        if begtrl>afew & endtrl<length(actgaze)-afew
            s1=nanmean(actgaze(:,begtrl-afew:begtrl),2); %pre-MS position
            s2=nanmean(actgaze(:,endtrl:endtrl+afew),2); %post-MS position
            % collect saccade size
            sizes=[sizes; pdist([s1'; s2'],'euclidean')];
            % calculate saccade direction
            sc=s2-s1;
            angles=[angles; atan2d(sc(2),sc(1))];
        end
    end
end    


end