FieldTrip Preprocessing

%% Montreal-Tutorial
% Steps:
% 1. Take a look at the data, what's in there?
% 2. Define the trials based on triggers
% 2.1. Read in the data trials
% 3. Take a look at the data again, what does it look like?
% 4. Throw out all you don't need
% 4.1. Remove TMS-Artefact
% 4.2. Filter
% 4.3. Downsample
% 4.4. Remove Blinks etc.
% 4.5. ICA

% Julian Keil, 13.04.2012
% julian.keil@gmail.com

%% Set Path and Data

inpath = ('/Users/juliankeil/Documents/Arbeit/ACN-Exchange/Tutorial/Tutorial Data'); % path to the data

cd (inpath) % GoTo path

data{1} = 'Pilot1.cnt'; % What's the dataset called?

%% 1. Take a look

cfg=[]; % empties the universally used cfg-structure
cfg.dataset=data{1}; % which dataset? Useful if a for-loop is used
cfg.trialdef.eventtype= '?'; % We don't know what events have happened
cfg.headerformat='ns_cnt32'; % Watch Out! With the neuroscan system, you have to set the bit-rate
cfg.dataformat='ns_cnt32'; % Otherwise you'll have corrupt data
cfg.eventformat='ns_cnt32';

cfg=ft_definetrial(cfg); % definetrail looks at all events in the dataset

%% 2. Define the trials based on triggers

cfg=[]; % Empty the cfg again, so you don't mess up stuff
cfg.dataset=data{1}; % which dataset? Useful if a for-loop is used
cfg.trialdef.eventtype='trigger'; % We now know that the trigger channel is called 'trigger'
%cfg.trialdef.eventvalue=6; % Only read in the trigger 13, more complex trial definitions need a trial-function
cfg.trialdef.prestim=1.5; % Seconds before the stimulus
cfg.trialdef.poststim=1.5; % Seconds after the stimulus
cfg.trialfun='trialfun_MEP'; % the normal trial-function is fine for now, otherwise call for an own
cfg.headerformat='ns_cnt32'; % Watch Out! With the neuroscan system, you have to set the bit-rate
cfg.dataformat='ns_cnt32'; % Otherwise you'll have corrupt data
cfg.eventformat='ns_cnt32';

cfg=ft_definetrial(cfg); % define trials based on these settings

%% 2.1. Actually read in the data. This might take a while. Get some fresh air.
% Watch out! Don't clear the cfg!

cfg.detrend='yes'; % De-Trending, so that we get rid of the offset and slow drift in the EEG-data
% cfg.hpfilter='yes'; % High-Pass-Filter On
% cfg.hpfreq=2; % High-Pass Frequenc
% cfg.hpfilttype='but'; % Filter Type Butterworth (IIR)
% cfg.hpfiltord=6; % % Filter Order, watch out for Instability!
% cfg.lpfilter='yes'; % Low-Pass Filter On
% cfg.lpfreq=100; % Low-Pass Frequency
% cfg.lpfilttype='but'; % Filter Type
% cfg.lpfiltord=16; % Filter Order
cfg.dftfilter='yes'; % Line Noise Filter using Discrete Fourier Transform
cfg.dftfreq=[30 60 90 120 180]; % Line Frequency (in Europe this would be 50)
%cfg.trials=[1 2 3];
tic
dat=ft_preprocessing(cfg); % Read in the trials definded above
toc
% Now we have the dat-structure containing the actual EEG-Data with the
% fields:
% hdr: Header information
% label: labels for the data-channels
% time: time-vector for each trial
% trial: actual data for each trial in chan-by-samplepoint
% fsample: sample-frequency
% sampleinfo: Beginning and End for each trial in sample points
% trialinfo: trigger value for each trial
% cfg: configuration that was used to call definetrial
%% 3. Take a look
load ~/subversion/tinnitus/matlab/Julian/NeuroScan64.mat % this contains the Neuroscan Layout.

cfg=[]; % the usual
cfg.channel={'EEG' 'CB1' 'CB2' '-CP3'};
cfg.viewmode='vertical'; % or butterfly
cfg.selectmode='eval'; % what should happen if you select a patch of data
cfg.selcfg.layout=lay; % we use the neuroscan Layout, more Layouts can be found in the fieldtrip %directory in the "templates"-Folder or you can create your own using ft_prepare_layout

ft_databrowser(cfg,dat); % we call the databrowser with the cfg-settings and the dat-structure defined above

%% 4. Now let's start trowing out all the stuff we don't need
% First clean the TMS-Artifact as it srews up the filtering and ICA
% Then Downsample the data, as you don't really need 2000 Hz sampling
% Then maybe filter a bit
% And remove artifacts by visual rejection
% Finally do some fancy ICA cleaning

%% 4.1. clean TMS artifact as described in Current Biology by Thut et al. 2011
cfg=[];
cfg.interp='spline'; % method for interpolation of the artifact
%cfg.pretime=10;
%cfg.posttime=50;

dat_cubic=cleanTMS(cfg,dat); % You'll end up with a new dataset

%% 4.2. Downsample
cfg=[];
cfg.resamplefs=300; % New Samplerate in Hz
cfg.detrend='no'; % We already detrended our data.

dat_res=ft_resampledata(cfg,dat);

%% 4.3. Filter a little (or a lot ;-))

cfg=[];
cfg.channel={'EEG' 'CB1' 'CB2'};
cfg.hpfilter='yes'; % High-Pass-Filter On
cfg.hpfreq=3; % High-Pass Frequenc
cfg.hpfilttype='but'; % Filter Type Butterworth (IIR)
cfg.hpfiltord=2; % % Filter Order, watch out for Instability!
cfg.lpfilter='yes'; % Low-Pass Filter On
cfg.lpfreq=25; % Low-Pass Frequency
cfg.lpfilttype='but'; % Filter Type
cfg.lpfiltord=16; % Filter Order

dat_filt=ft_preprocessing(cfg,dat_res); % Here we call again for preprocessing,
%but this time with an actual dataset at and. Don't get confused, we're
%just using the same function again. Also, you might want to check the
%filter settings so this doesn't screw up the data!

%% 4.4. Visual Artifact Rejection
% Watch out!
% I fyou get an error when plotting the channels, you might have to set
% cfg.mp.interactive='yes'; in rejectvisual_summary. It's a but, not a
% feature.

cfg=[];
cfg.channel={'EEG' 'CB1' 'CB2'}; % We'll only look at the EEG-Channels, the rest at the moment doesn't matter to us
cfg.method='summary';% We'll look at all the info we can get
cfg.latency=[-.7 .7];% Set the latenc
cfg.keepchannel='yes';% If you want to interpolate channels, you have to keep the bad ones!
cfg.layout=lay; % Also, you have to actually set CB1 and CB2 as EEG

dat_clear=ft_rejectvisual(cfg,dat_res);

%% 4.4.1. Interpolate Channels
% First you have to define, which channels are neighbors to which channels
% Then you can interpolate from the surrounding channels
% However, we need to do this in 3D-space, so we need the actual
% 3-Dimensional Senesor positions:

elec.pnt=Scan64_pos; % the 3D-Positions
elec.label=Scan64_lab; % and the names
%% 4.4.1.1. Neighbors
cfg=[];
cfg.method='distance'; % how should the neighbors be selected?
cfg.neighbourdist=50; % I have no Idea what range this has, just make sure, that you get meaningful neighbors
cfg.elec=elec; % Here we need the 3d-positions!
dummy=ft_prepare_neighbours(cfg); % between 5 and 10 neighbors is a good value, always good to check!

%% 4.4.1.2. Check!
cfg=[];
cfg.neighbours=dummy; % what neighbor-structure
cfg.elec=elec; % what layout
ft_neighbourplot(cfg)
% Again, ist's best to check the actual neighbors for each channel. On
% average you'll want to end up with ahout 5 to 10 neighbors for each
% channel, at least have 2, otherwise you'll just end up copying one
% channel.

%% 4.4.1.3. Repair
cfg=[];
cfg.badchannel={'FT8' 'CP3' 'FP1'}; % Who's bad?
cfg.neighbours=dummy; % What channels should be used to fix?

dat_clear.elec=elec; % and we need the 3d positions

dat_fix=ft_channelrepair(cfg,dat_clear);

%% 4.5. ICA
% in order to get rid of blink artefacts and external noise, the
% independent component analysis can be useful. But Beware, don't exclude
% to much! Whereas teh ICA can split the signal into independent
% components, this does not mean that one component does not contain
% information from more than one source!

%% 4.5.1. Compute ICA-Components
cfg=[];
cfg.channel={'EEG' 'CB1' 'CB2'}; % Super important to only use the EEG-Channels, otherwise the ICA won't work
cfg.method='fastica'; % Whcih method should be used?
cfg.numcomponent=50; % if the ICA does not converge, limit the number of components to compute
cfg.demean='no';
comp=ft_componentanalysis(cfg,dat_fix);

%% 4.5.2. Take a look at the components

cfg=[];
cfg.layout=lay;
cfg.viewmode='component'; % same as above, just with a different mode
ft_databrowser(cfg,comp);

%% 4.5.3. And take out the ones that are clearly artefacts

cfg=[];
cfg.component=[9]; % I focus on blinks and muscle artefacts
dat_ica=ft_rejectcomponent(cfg,comp,dat_fix);


%% 4.5.4 Compare before & after
figure;plot(dat_ica.time{1},dat_ica.trial{1}(22,:),'r');hold;plot(dat_clear.time{1},dat_clear.trial{1}(22,:),'b')

%% 6 Save
%mkdir 2ndtry
cd 2ndtry
save dat_orig.mat dat -V7.3
%save dat_cubic.mat dat_cubic -V7.3
save dat_clear.mat dat_clear -V7.3
%save dat_filt.mat dat_filt -V7.3
save dat_res.mat dat_res -V7.3
save dat_ica.mat dat_ica comp -V7.3
save dat_fix.mat dat_fix -V7.3