pupil_preprocessing_matlab/pupilTutorial.m at master · DonnerLab/pupil_preprocessing_matlab · GitHub

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
%% PUPIL PREPROCESSING
%
% 1. convert edf file to asc
% 2. create a FieldTrip-style data structure
% 3. interpolate Eyelink-defined and additionally detected blinks
% (optional): regress out blink- and saccade-linked pupil response
% 4. epoch using a custom trial-definition function
% 5. plot event-related pupil responses
%
% Anne Urai, 2016

% setup path
clear; clc; close all;
thispath = '~/Dropbox/code/PupilPreprocessing/';
addpath(thispath);

% first, define the asc filename
edfFile = 'EL_P22_s5_b1_2015-06-28_11-14-27.edf';
ascFile = regexprep(edfFile, 'edf', 'asc');

% set path to FieldTrip - get this from http://www.fieldtriptoolbox.org/download
addpath('~/Documents/fieldtrip/');
ft_defaults;

% ============================================== %
% 1. convert edf file to asc
% ============================================== %

if ~exist(ascFile, 'file'),
    if ismac,
        edf2ascPath = [thispath '/edf2asc-mac'];
    elseif isunix,
        edf2ascPath = [thispath '/edf2asc-linux'];
    else
        error('Sorry, I don''t have an edf2asc converter for Windows')
    end
    % use this converter to create the asc file
    % failsafe mode avoids errors when some samples are missing
    system(sprintf('%s %s -failsafe -input', edf2ascPath, edfFile));
end
assert(exist(ascFile, 'file') > 1, 'Edf not properly converted');

% ============================================== %
% 2. create a FieldTrip-style data structure
% ============================================== %

% read in the asc EyeLink file
asc = read_eyelink_ascNK_AU(ascFile);

% create events and data structure, parse asc
[data, event, blinksmp, saccsmp] = asc2dat(asc);

% ============================================== %
% 3. interpolate Eyelink-defined and additionally detected blinks
% ============================================== %

plotMe = true;
newpupil = blink_interpolate(data, blinksmp, plotMe);
data.trial{1}(find(strcmp(data.label, 'EyePupil')==1),:) = newpupil;

% ============================================== %
% (optional): regress out blink- and saccade-linked pupil response
% http://journals.plos.org/plosone/article?id=10.1371/journal.pone.0155574
% ============================================== %

pupildata = data.trial{1}(~cellfun(@isempty, strfind(lower(data.label), 'eyepupil')),:);
newpupil = blink_regressout(pupildata, data.fsample, blinksmp, saccsmp, 1, 1);
% put back in fieldtrip format
data.trial{1}(~cellfun(@isempty, strfind(lower(data.label), 'eyepupil')),:) = newpupil;

% zscore since we work with the bandpassed signal
data.trial{1}(find(strcmp(data.label, 'EyePupil')==1),:) = ...
    zscore(data.trial{1}(find(strcmp(data.label, 'EyePupil')==1),:));

% ==================================================================
% define trials
% ==================================================================

cfg                         = [];
cfg.dataset                 = ascFile;
cfg.event                   = event;
cfg.trialfun                = 'my_trialfun';
cfg.trialdef.pre            = 0;
cfg.trialdef.post           = 2;
cfg.fsample                 = asc.fsample;
cfg.sj                      = 22;
cfg.session                 = 5;
cfg.block                   = 1;
cfg                         = ft_definetrial(cfg);
% epoch
data                        = ft_redefinetrial(cfg, data);

% note that this trialfun defines four events per trial: two stimuli (that
% have to be compared), a response, and feedback). We have the trial start
% at fixation and save the samples of each of those events (as well as more
% information about stimulus identity, choice and accuracy) in
% data.trialinfo:
% [fixoffset refoffset stimtype stimoffset resptype respcorrect respoffset ...
% feedbacktype feedbackoffset trlcnt blockcnt session]

% ==================================================================
% downsample before saving
% ==================================================================

cfg             = [];
cfg.resamplefs  = 100;
cfg.fsample     = data.fsample;

% see Niels' message on the FT mailing list
samplerows = find(data.trialinfo(1,:)>100); % indices of the rows with sample values (and not event codes)
data.trialinfo(:,samplerows) = round(data.trialinfo(:,samplerows) * (cfg.resamplefs/cfg.fsample));

% use fieldtrip to resample
data    = ft_resampledata(cfg, data);

% ==================================================================
% visualize the pupil timecourse
% ==================================================================

cfg                 = [];
cfg.channel         = 'EyePupil';
cfg.trials(1).name  = 'correct';
cfg.trials(1).idx   = find(data.trialinfo(:, 8) == 1);
cfg.trials(2).name  = 'error';
cfg.trials(2).idx   = find(data.trialinfo(:, 8) == 0);

clf; subplot(2,2,[1 2]);
plotEventRelated(cfg, data);
title(ascFile, 'interpreter', 'none');
ylabel('Pupil response (z)');

% ==================================================================
% save file
% ==================================================================

filename = regexprep(ascFile, 'asc', 'mat');
save(filename, 'data');