Microstim_FF-FB/SCRIPTS/FUNCTIONS/SACCADE_Analysis.m

%% SACCADE_Analysis.m
EYE_DataFile = fullfile(fld.data,'EYE.mat');
load(EYE_DataFile)

%% Words meaning ----------------------------------------------------------
%1-4 Homogenous Texture [noMS SOA1 SOA2 SOA3]
%5-8 Figure position #1 [noMS SOA1 SOA2 SOA3]
%9-12 Figure position #2 [noMS SOA1 SOA2 SOA3]

%% Process the data -------------------------------------------------------
wNoMs = [1 5 9]; wWithMS= [4 8 12];

fixations.noMS.X = []; fixations.noMS.Y = [];
fixations.withMS.X = []; fixations.withMS.Y = [];
fixations.noMS.D = []; fixations.withMS.D = [];
fixations.noMS.V = []; fixations.withMS.V = [];
fixations.noMS.microsacc = []; fixations.withMS.microsacc = [];
fixations.noMS.microsaccn = []; fixations.withMS.microsaccn = [];

TANKS_DATA = NEW_PERFORMANCES.FEEDBACK.TANKS;
for tank = 1:size(TANKS_DATA,2)
    performance = TANKS_DATA(tank).newPerformance;
    wordData    = TANKS_DATA(tank).oldTrials(:,7);
    
    % ----- NO MICROSTIMULATION ------------
    indNoMS = ismember(wordData,wNoMs) & performance == 1;
    fixations.noMS.X = [fixations.noMS.X; ...
        nanmean(TANKS_DATA(tank).gazeXYCoords(indNoMS,1))];
    fixations.noMS.Y = [fixations.noMS.Y; ...
        nanmean(TANKS_DATA(tank).gazeXYCoords(indNoMS,2))];
    fixations.noMS.D = [fixations.noMS.D; ...
        nanmean(TANKS_DATA(tank).gazeAmplitude(indNoMS,1))];
    
    mEyeShift=[]; cEyeShift=[]; sdEyePos=[]; 
    sz_micsacc = []; num_microsacc = [];
    for tr=find(indNoMS==1)'
        ev = TANKS_DATA(tank).gazePos{tr};
        if ~isempty(ev)
            mEyeShift = [mEyeShift; ...
                nanmean(sqrt(diff(ev(:,1)).^2 + diff(ev(:,2)).^2))];
            cEyeShift = [cEyeShift; ...
                sum(sqrt(diff(ev(:,1)).^2 + diff(ev(:,2)).^2))];
            sdEyePos = [sdEyePos; std(ev(:,1)) std(ev(:,2)) tank tr];
            
            % get microsaccades
            micsacc = micsaccdeg([zeros(length(ev),1) ev], Sampf);
            n_micsacc = size(micsacc,1);
            if ~isempty(micsacc)
                sz_micsacc = [sz_micsacc; ...
                    micsacc(:,4) tank*ones(size(micsacc,1),1) ...
                    tr*ones(size(micsacc,1),1)...
                    ]; 
                num_microsacc=[num_microsacc; n_micsacc];
            else
                num_microsacc=[num_microsacc; 0];
            end
        end
    end
    fixations.noMS.V = [fixations.noMS.V; mEyeShift cEyeShift sdEyePos];
    fixations.noMS.microsacc = [fixations.noMS.microsacc; sz_micsacc];
    fixations.noMS.microsaccn = [fixations.noMS.microsaccn; num_microsacc];

    % ----- WITH MICROSTIMULATION ------------
    indWithMS = ismember(wordData,wWithMS) & performance == 1;
    fixations.withMS.X = [fixations.withMS.X;  ...
        nanmean(TANKS_DATA(tank).gazeXYCoords(indWithMS,1))];
    fixations.withMS.Y = [fixations.withMS.Y;  ...
        nanmean(TANKS_DATA(tank).gazeXYCoords(indWithMS,2))];
    fixations.withMS.D = [fixations.withMS.D; ...
        nanmean(TANKS_DATA(tank).gazeAmplitude(indWithMS,1))];
    
    mEyeShift=[]; cEyeShift=[]; sdEyePos=[]; 
    sz_micsacc = []; num_microsacc = [];
    for tr=find(indWithMS==1)'
        ev = TANKS_DATA(tank).gazePos{tr};
        if ~isempty(ev)
            mEyeShift = [mEyeShift; ...
                nanmean(sqrt(diff(ev(:,1)).^2 + diff(ev(:,2)).^2))];
            cEyeShift = [cEyeShift; ...
                sum(sqrt(diff(ev(:,1)).^2 + diff(ev(:,2)).^2))];
            sdEyePos = [sdEyePos; ...
                std(ev(:,1)) std(ev(:,2)) tank tr];
            
            % get microsaccades
            micsacc = micsaccdeg([zeros(length(ev),1) ev], Sampf);
            n_micsacc = size(micsacc,1);
            if ~isempty(micsacc)
                sz_micsacc = [sz_micsacc; ...
                    micsacc(:,4) tank*ones(size(micsacc,1),1) ...
                    tr*ones(size(micsacc,1),1)...
                    ]; 
                num_microsacc=[num_microsacc; n_micsacc];
            else
                num_microsacc=[num_microsacc; 0];
            end
        end
    end
    fixations.withMS.V = [fixations.withMS.V; mEyeShift cEyeShift sdEyePos];
    fixations.withMS.microsacc = [fixations.withMS.microsacc; sz_micsacc];
    fixations.withMS.microsaccn = [fixations.withMS.microsaccn; num_microsacc];

end

TANKS_DATA = NEW_PERFORMANCES.FEEDFORWARD.TANKS;
for tank = 1:size(TANKS_DATA,2)
    performance = TANKS_DATA(tank).newPerformance;
    wordData    = TANKS_DATA(tank).oldTrials(:,7);
    
    % ----- NO MICROSTIMULATION ------------
    indNoMS = ismember(wordData,wNoMs) & performance == 1;
    fixations.noMS.X = [fixations.noMS.X; ...
        nanmean(TANKS_DATA(tank).gazeXYCoords(indNoMS,1))];
    fixations.noMS.Y = [fixations.noMS.Y; ...
        nanmean(TANKS_DATA(tank).gazeXYCoords(indNoMS,2))];
    fixations.noMS.D = [fixations.noMS.D; ...
        nanmean(TANKS_DATA(tank).gazeAmplitude(indNoMS,1))];
    
    mEyeShift=[]; cEyeShift=[]; sdEyePos=[]; 
    sz_micsacc = []; num_microsacc = [];
    for tr=find(indNoMS==1)'
        ev = TANKS_DATA(tank).gazePos{tr};
        if ~isempty(ev)
            mEyeShift = [mEyeShift; ...
                nanmean(sqrt(diff(ev(:,1)).^2 + diff(ev(:,2)).^2))];
            cEyeShift = [cEyeShift; ...
                sum(sqrt(diff(ev(:,1)).^2 + diff(ev(:,2)).^2))];
            sdEyePos = [sdEyePos; std(ev(:,1)) std(ev(:,2)) tank tr];
            
            % get microsaccades
            micsacc = micsaccdeg([zeros(length(ev),1) ev], Sampf);
            n_micsacc = size(micsacc,1);
            if ~isempty(micsacc)
                sz_micsacc = [sz_micsacc; ...
                    micsacc(:,4) tank*ones(size(micsacc,1),1) ...
                    tr*ones(size(micsacc,1),1)...
                    ]; 
                num_microsacc=[num_microsacc; n_micsacc];
            else
                num_microsacc=[num_microsacc; 0];
            end
        end
    end
    fixations.noMS.V = [fixations.noMS.V; mEyeShift cEyeShift sdEyePos];
    fixations.noMS.microsacc = [fixations.noMS.microsacc; sz_micsacc];
    fixations.noMS.microsaccn = [fixations.noMS.microsaccn; num_microsacc];

    % ----- WITH MICROSTIMULATION ------------
    indWithMS = ismember(wordData,wWithMS) & performance == 1;
    fixations.withMS.X = [fixations.withMS.X;  ...
        nanmean(TANKS_DATA(tank).gazeXYCoords(indWithMS,1))];
    fixations.withMS.Y = [fixations.withMS.Y;  ...
        nanmean(TANKS_DATA(tank).gazeXYCoords(indWithMS,2))];
    fixations.withMS.D = [fixations.withMS.D; ...
        nanmean(TANKS_DATA(tank).gazeAmplitude(indWithMS,1))];
    
    mEyeShift=[]; cEyeShift=[]; sdEyePos=[]; 
    sz_micsacc = []; num_microsacc = [];
    for tr=find(indWithMS==1)'
        ev = TANKS_DATA(tank).gazePos{tr};
        if ~isempty(ev)
            mEyeShift = [mEyeShift; ...
                nanmean(sqrt(diff(ev(:,1)).^2 + diff(ev(:,2)).^2))];
            cEyeShift = [cEyeShift; ...
                sum(sqrt(diff(ev(:,1)).^2 + diff(ev(:,2)).^2))];
            sdEyePos = [sdEyePos; std(ev(:,1)) std(ev(:,2)) tank tr];
            
            % get microsaccades
            micsacc = micsaccdeg([zeros(length(ev),1) ev], Sampf);
            n_micsacc = size(micsacc,1);
            if ~isempty(micsacc)
                sz_micsacc = [sz_micsacc; ...
                    micsacc(:,4) tank*ones(size(micsacc,1),1) ...
                    tr*ones(size(micsacc,1),1)...
                    ]; 
                num_microsacc=[num_microsacc; n_micsacc];
            else
                num_microsacc=[num_microsacc; 0];
            end
        end
    end
    fixations.withMS.V = [fixations.withMS.V; mEyeShift cEyeShift sdEyePos];
    fixations.withMS.microsacc = [fixations.withMS.microsacc; sz_micsacc];
    fixations.withMS.microsaccn = [fixations.withMS.microsaccn; num_microsacc];
end

%% Plot fixations X distribution ------------------------------------------

bins = [-0.3:0.01:0.3];
[cNoMS,xNoMS] = hist(fixations.noMS.X,bins);fNoMS = cNoMS/sum(cNoMS);
[cWithMS,xWithMS] = hist(fixations.withMS.X,bins);
fWithMS = cWithMS/sum(cWithMS);
figure(1);clf;
subplot(2,2,1);hold all;
set(gca,'FontSize',12)
plot(xNoMS,fNoMS,'Color',[0.85 0.18 0.17])
plot(xWithMS,fWithMS,'Color',[0.35 0.35 0.35])
xlabel('Average Fixations X')
ylabel('Frequency')
legend({'noMS' 'withMS'})

[cNoMS,xNoMS] = hist(fixations.noMS.Y,bins);fNoMS = cNoMS/sum(cNoMS);
[cWithMS,xWithMS] = hist(fixations.withMS.Y,bins);
fWithMS = cWithMS/sum(cWithMS);
subplot(2,2,2);hold all;
set(gca,'FontSize',12)
plot(xNoMS,fNoMS,'Color',[0.85 0.18 0.17])
plot(xWithMS,fWithMS,'Color',[0.35 0.35 0.35])
xlabel('Average Fixations Y')
legend({'noMS' 'withMS'})

%% Plot one to one correspondance -----------------------------------------
subplot(2,2,3)
hold all;

for ses = 1:length(fixations.noMS.X)
    plot([1 2],[fixations.noMS.X(ses) fixations.withMS.X(ses)],'-k')
    scatter(1,fixations.noMS.X(ses),15,[1 0 0])
    scatter(2,fixations.withMS.X(ses),15,[1 0 0])
end
xlim([0.7 2.3])
set(gca,'XTick',[1 2])
set(gca,'XTickLabel',{'noMS' 'witMS'})

subplot(2,2,4)
hold all;

for ses = 1:length(fixations.noMS.Y)
    plot([1 2],[fixations.noMS.Y(ses) fixations.withMS.Y(ses)],'-k')
    scatter(1,fixations.noMS.Y(ses),15,[1 0 0])
    scatter(2,fixations.withMS.Y(ses),15,[1 0 0])
end
xlim([0.7 2.3])
set(gca,'XTick',[1 2])
set(gca,'XTickLabel',{'noMS' 'witMS'})

%% Plot average amplitude -------------------------------------------------
f=figure;
subplot(1,4,1);hold on;
errorbar([1 2],[nanmean(fixations.noMS.D) nanmean(fixations.withMS.D)],...
    [nanstd(fixations.noMS.D) nanstd(fixations.withMS.D)],'ks',...
    'linestyle','none','linewidth',2,'markersize',12,'markerfacecolor','k');
set(gca, 'xlim',[0 3],'ylim',[0 0.3],'xtick',1:2,'xticklabels',{'no MS','MS'})
ylabel('Gaze amplitude');
title('Mean +/- stdev')

subplot(1,4,2:4);hold on;
[y1,x1]=hist(fixations.noMS.D,100); y1b=y1/sum(y1);
[y2,x2]=hist(fixations.withMS.D,100); y2b=y2/sum(y2);
bar(x1,y1b,'r','Edgecolor','none','barwidth',1);
bar(x2,-y2b,'k','Edgecolor','none','barwidth',1);
xlabel('Gaze Amplitude')
ylabel('Frequency')
legend({'noMS' 'withMS'})
title('Gaze amplitude distributions')

%% stats on amplitude -----------------------------------------------------
gazemat = [fixations.noMS.D;fixations.withMS.D];
group = {[zeros(size(fixations.noMS.D)); ones(size(fixations.withMS.D))]};
gn={'MS'};
[P,T,STATS,TERMS]=anovan(gazemat,group)

%% Plot average eyeshifts [averaged per trial] ----------------------------
f2=figure;
subplot(2,4,1);hold on;
errorbar([1 2],[nanmean(fixations.noMS.V(:,1)) nanmean(fixations.withMS.V(:,1))],...
    [nanstd(fixations.noMS.V(:,1)) nanstd(fixations.withMS.V(:,1))],'ks',...
    'linestyle','none','linewidth',2,'markersize',12,'markerfacecolor','k');
set(gca,'xlim',[0 3],'ylim',[0 0.02],'xtick',1:2,'xticklabels',{'no MS','MS'})
ylabel('Aver. Gaze Shift');
title('Mean +/- stdev')

subplot(2,4,2:4);hold on;
[y1,x1]=hist(fixations.noMS.V(:,1),100); y1b=y1/sum(y1);
[y2,x2]=hist(fixations.withMS.V(:,1),100); y2b=y2/sum(y2);
bar(x1,y1b,'r','Edgecolor','none','barwidth',1);
bar(x2,-y2b,'k','Edgecolor','none','barwidth',1);
xlabel('Aver. Gaze Shifts')
ylabel('Frequency')
legend({'noMS' 'withMS'})
title('Average gaze shift distributions')

%% stats on amplitude -----------------------------------------------------
gazemat = [fixations.noMS.V(:,1);fixations.withMS.V(:,1)];
group = {[zeros(size(fixations.noMS.V(:,1))); ...
    ones(size(fixations.withMS.V(:,1)))]};
gn={'MS'};
[P,T,STATS,TERMS]=anovan(gazemat,group)

%% Plot cumulative eyeshifts per trial ------------------------------------
subplot(2,4,5);hold on;
errorbar([1 2],[nanmean(fixations.noMS.V(:,2)) ...
    nanmean(fixations.withMS.V(:,2))],...
    [nanstd(fixations.noMS.V(:,2)) ...
    nanstd(fixations.withMS.V(:,2))],'ks',...
    'linestyle','none','linewidth',2,'markersize',12,...
    'markerfacecolor','k');
set(gca,'xlim',[0 3],'ylim',[0 6],'xtick',1:2,'xticklabels',{'no MS','MS'})
ylabel('Cum. Gaze shift');
title('Mean +/- stdev')

subplot(2,4,6:8);hold on;
[y1,x1]=hist(fixations.noMS.V(:,2),100); y1b=y1/sum(y1);
[y2,x2]=hist(fixations.withMS.V(:,2),100); y2b=y2/sum(y2);
bar(x1,y1b,'r','Edgecolor','none','barwidth',1);
bar(x2,-y2b,'k','Edgecolor','none','barwidth',1);
xlabel('Cum. Gaze Shift')
ylabel('Frequency')
legend({'noMS' 'withMS'})
title('Cumulative gaze shift distributions')

%% stats on amplitude -----------------------------------------------------
gazemat = [fixations.noMS.V(:,2);fixations.withMS.V(:,2)];
group = {[zeros(size(fixations.noMS.V(:,2))); ...
    ones(size(fixations.withMS.V(:,2)))]};
gn={'MS'};
[P,T,STATS,TERMS]=anovan(gazemat,group)

%% Plot microsaccades amplitudes and number per trial ---------------------
figure;
subplot(2,4,1);hold on;
errorbar([1 2],[nanmean(fixations.noMS.microsacc(:,1)) ...
    nanmean(fixations.withMS.microsacc(:,1))],...
    [nanstd(fixations.noMS.microsacc(:,1)) ...
    nanstd(fixations.withMS.microsacc(:,1))],'ks',...
    'linestyle','none','linewidth',2,'markersize',12,'markerfacecolor','k');
set(gca,'xlim',[0 3],'ylim',[0 0.4],'xtick',1:2,'xticklabel',{'no MS','MS'})
ylabel('Microsacc amplitude');
title('Mean +/- stdev')

subplot(2,4,2:4);hold on;
[y1,x1]=hist(fixations.noMS.microsacc(:,1),100); y1b=y1/sum(y1);
[y2,x2]=hist(fixations.withMS.microsacc(:,1),100); y2b=y2/sum(y2);
bar(x1,y1b,'r','Edgecolor','none','barwidth',1);
bar(x2,-y2b,'k','Edgecolor','none','barwidth',1);
set(gca,'xlim',[0 0.6],'ylim',[-0.12 0.12]);
xlabel('Microsaccade amplitude')
ylabel('Frequency')
legend({'no MS' 'MS'})
title('Microsaccade amplitude distributions')

subplot(2,4,5);hold on;
errorbar([1 2],[nanmean(fixations.noMS.microsaccn(:,1)) ...
    nanmean(fixations.withMS.microsaccn(:,1))],...
    [nanstd(fixations.noMS.microsaccn(:,1)) ...
    nanstd(fixations.withMS.microsaccn(:,1))],'ks',...
    'linestyle','none','linewidth',2,'markersize',12,'markerfacecolor','k');
set(gca,'xlim',[0 3],'ylim',[-3 8],'xtick',1:2,'xticklabel',{'no MS','MS'})
ylabel('# microsacc / trial');
title('Mean +/- stdev')

subplot(2,4,6:8);hold on;
[y1,x1]=hist(fixations.noMS.microsaccn(:,1),0:30); y1b=y1/sum(y1);
[y2,x2]=hist(fixations.withMS.microsaccn(:,1),0:30); y2b=y2/sum(y2);
bar(x1,y1b,'r','Edgecolor','none','barwidth',1);
bar(x2,-y2b,'k','Edgecolor','none','barwidth',1);
set(gca,'xlim',[-1 15]);
xlabel('# microsaccades per trial')
ylabel('Frequency')
legend({'no MS' 'MS'})
title('Number of microsaccades')

%% stats on microsccades --------------------------------------------------
% amplitudes
microsacmat = [fixations.noMS.microsacc(:,1);...
    fixations.withMS.microsacc(:,1)];
group = {[zeros(size(fixations.noMS.microsacc(:,1))); ...
    ones(size(fixations.withMS.microsacc(:,1)))]};
gn={'MS'}; [P,T,STATS,TERMS]=anovan(microsacmat,group)
% number
microsacmat2 = [fixations.noMS.microsaccn(:,1);...
    fixations.withMS.microsaccn(:,1)];
group = {[zeros(size(fixations.noMS.microsaccn(:,1))); ...
    ones(size(fixations.withMS.microsaccn(:,1)))]};
gn={'MS'}; [P2,T2,STATS2,TERMS2]=anovan(microsacmat2,group)