juliankeil

FieldTrip Preprocessing

Aug 8th, 2012
1,680
Never
Not a member of Pastebin yet? Sign Up, it unlocks many cool features!
  1. %% Montreal-Tutorial
  2. % Steps:
  3. % 1. Take a look at the data, what's in there?
  4. % 2. Define the trials based on triggers
  5. % 2.1. Read in the data trials
  6. % 3. Take a look at the data again, what does it look like?
  7. % 4. Throw out all you don't need
  8. % 4.1. Remove TMS-Artefact
  9. % 4.2. Filter
  10. % 4.3. Downsample
  11. % 4.4. Remove Blinks etc.
  12. % 4.5. ICA
  13.  
  14. % Julian Keil, 13.04.2012
  15. % julian.keil@gmail.com
  16.  
  17. %% Set Path and Data
  18.  
  19. inpath = ('/Users/juliankeil/Documents/Arbeit/ACN-Exchange/Tutorial/Tutorial Data'); % path to the data
  20.  
  21. cd (inpath) % GoTo path
  22.  
  23. data{1} = 'Pilot1.cnt'; % What's the dataset called?
  24.  
  25. %% 1. Take a look
  26.  
  27. cfg=[]; % empties the universally used cfg-structure
  28. cfg.dataset=data{1}; % which dataset? Useful if a for-loop is used
  29. cfg.trialdef.eventtype= '?'; % We don't know what events have happened
  30. cfg.headerformat='ns_cnt32'; % Watch Out! With the neuroscan system, you have to set the bit-rate
  31. cfg.dataformat='ns_cnt32'; % Otherwise you'll have corrupt data
  32. cfg.eventformat='ns_cnt32';
  33.  
  34. cfg=ft_definetrial(cfg); % definetrail looks at all events in the dataset
  35.  
  36. %% 2. Define the trials based on triggers
  37.  
  38. cfg=[]; % Empty the cfg again, so you don't mess up stuff
  39. cfg.dataset=data{1}; % which dataset? Useful if a for-loop is used
  40. cfg.trialdef.eventtype='trigger'; % We now know that the trigger channel is called 'trigger'
  41. %cfg.trialdef.eventvalue=6; % Only read in the trigger 13, more complex trial definitions need a trial-function
  42. cfg.trialdef.prestim=1.5; % Seconds before the stimulus
  43. cfg.trialdef.poststim=1.5; % Seconds after the stimulus
  44. cfg.trialfun='trialfun_MEP'; % the normal trial-function is fine for now, otherwise call for an own
  45. cfg.headerformat='ns_cnt32'; % Watch Out! With the neuroscan system, you have to set the bit-rate
  46. cfg.dataformat='ns_cnt32'; % Otherwise you'll have corrupt data
  47. cfg.eventformat='ns_cnt32';
  48.  
  49. cfg=ft_definetrial(cfg); % define trials based on these settings
  50.  
  51. %% 2.1. Actually read in the data. This might take a while. Get some fresh air.
  52. % Watch out! Don't clear the cfg!
  53.  
  54. cfg.detrend='yes'; % De-Trending, so that we get rid of the offset and slow drift in the EEG-data
  55. % cfg.hpfilter='yes'; % High-Pass-Filter On
  56. % cfg.hpfreq=2; % High-Pass Frequenc
  57. % cfg.hpfilttype='but'; % Filter Type Butterworth (IIR)
  58. % cfg.hpfiltord=6; % % Filter Order, watch out for Instability!
  59. % cfg.lpfilter='yes'; % Low-Pass Filter On
  60. % cfg.lpfreq=100; % Low-Pass Frequency
  61. % cfg.lpfilttype='but'; % Filter Type
  62. % cfg.lpfiltord=16; % Filter Order
  63. cfg.dftfilter='yes'; % Line Noise Filter using Discrete Fourier Transform
  64. cfg.dftfreq=[30 60 90 120 180]; % Line Frequency (in Europe this would be 50)
  65. %cfg.trials=[1 2 3];
  66. tic
  67. dat=ft_preprocessing(cfg); % Read in the trials definded above
  68. toc
  69. % Now we have the dat-structure containing the actual EEG-Data with the
  70. % fields:
  71. % hdr: Header information
  72. % label: labels for the data-channels
  73. % time: time-vector for each trial
  74. % trial: actual data for each trial in chan-by-samplepoint
  75. % fsample: sample-frequency
  76. % sampleinfo: Beginning and End for each trial in sample points
  77. % trialinfo: trigger value for each trial
  78. % cfg: configuration that was used to call definetrial
  79. %% 3. Take a look
  80. load ~/subversion/tinnitus/matlab/Julian/NeuroScan64.mat % this contains the Neuroscan Layout.
  81.  
  82. cfg=[]; % the usual
  83. cfg.channel={'EEG' 'CB1' 'CB2' '-CP3'};
  84. cfg.viewmode='vertical'; % or butterfly
  85. cfg.selectmode='eval'; % what should happen if you select a patch of data
  86. 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
  87.  
  88. ft_databrowser(cfg,dat); % we call the databrowser with the cfg-settings and the dat-structure defined above
  89.  
  90. %% 4. Now let's start trowing out all the stuff we don't need
  91. % First clean the TMS-Artifact as it srews up the filtering and ICA
  92. % Then Downsample the data, as you don't really need 2000 Hz sampling
  93. % Then maybe filter a bit
  94. % And remove artifacts by visual rejection
  95. % Finally do some fancy ICA cleaning
  96.  
  97. %% 4.1. clean TMS artifact as described in Current Biology by Thut et al. 2011
  98. cfg=[];
  99. cfg.interp='spline'; % method for interpolation of the artifact
  100. %cfg.pretime=10;
  101. %cfg.posttime=50;
  102.  
  103. dat_cubic=cleanTMS(cfg,dat); % You'll end up with a new dataset
  104.            
  105. %% 4.2. Downsample
  106. cfg=[];
  107. cfg.resamplefs=300; % New Samplerate in Hz
  108. cfg.detrend='no'; % We already detrended our data.
  109.  
  110. dat_res=ft_resampledata(cfg,dat);
  111.  
  112. %% 4.3. Filter a little (or a lot ;-))
  113.  
  114. cfg=[];
  115. cfg.channel={'EEG' 'CB1' 'CB2'};
  116. cfg.hpfilter='yes'; % High-Pass-Filter On
  117. cfg.hpfreq=3; % High-Pass Frequenc
  118. cfg.hpfilttype='but'; % Filter Type Butterworth (IIR)
  119. cfg.hpfiltord=2; % % Filter Order, watch out for Instability!
  120. cfg.lpfilter='yes'; % Low-Pass Filter On
  121. cfg.lpfreq=25; % Low-Pass Frequency
  122. cfg.lpfilttype='but'; % Filter Type
  123. cfg.lpfiltord=16; % Filter Order
  124.  
  125. dat_filt=ft_preprocessing(cfg,dat_res); % Here we call again for preprocessing,
  126. %but this time with an actual dataset at and. Don't get confused, we're
  127. %just using the same function again. Also, you might want to check the
  128. %filter settings so this doesn't screw up the data!
  129.  
  130. %% 4.4. Visual Artifact Rejection
  131. % Watch out!
  132. % I fyou get an error when plotting the channels, you might have to set
  133. % cfg.mp.interactive='yes'; in rejectvisual_summary. It's a but, not a
  134. % feature.
  135.  
  136. cfg=[];
  137. cfg.channel={'EEG' 'CB1' 'CB2'}; % We'll only look at the EEG-Channels, the rest at the moment doesn't matter to us
  138. cfg.method='summary';% We'll look at all the info we can get
  139. cfg.latency=[-.7 .7];% Set the latenc
  140. cfg.keepchannel='yes';% If you want to interpolate channels, you have to keep the bad ones!
  141. cfg.layout=lay; % Also, you have to actually set CB1 and CB2 as EEG
  142.  
  143. dat_clear=ft_rejectvisual(cfg,dat_res);
  144.  
  145. %% 4.4.1. Interpolate Channels
  146. % First you have to define, which channels are neighbors to which channels
  147. % Then you can interpolate from the surrounding channels
  148. % However, we need to do this in 3D-space, so we need the actual
  149. % 3-Dimensional Senesor positions:
  150.  
  151. elec.pnt=Scan64_pos; % the 3D-Positions
  152. elec.label=Scan64_lab; % and the names
  153. %% 4.4.1.1. Neighbors
  154. cfg=[];
  155. cfg.method='distance'; % how should the neighbors be selected?
  156. cfg.neighbourdist=50; % I have no Idea what range this has, just make sure, that you get meaningful neighbors
  157. cfg.elec=elec; % Here we need the 3d-positions!
  158. dummy=ft_prepare_neighbours(cfg); % between 5 and 10 neighbors is a good value, always good to check!
  159.    
  160. %% 4.4.1.2. Check!
  161. cfg=[];
  162. cfg.neighbours=dummy; % what neighbor-structure
  163. cfg.elec=elec; % what layout
  164. ft_neighbourplot(cfg)
  165. % Again, ist's best to check the actual neighbors for each channel. On
  166. % average you'll want to end up with ahout 5 to 10 neighbors for each
  167. % channel, at least have 2, otherwise you'll just end up copying one
  168. % channel.
  169.  
  170. %% 4.4.1.3. Repair
  171. cfg=[];
  172. cfg.badchannel={'FT8' 'CP3' 'FP1'}; % Who's bad?
  173. cfg.neighbours=dummy; % What channels should be used to fix?
  174.  
  175. dat_clear.elec=elec; % and we need the 3d positions
  176.  
  177. dat_fix=ft_channelrepair(cfg,dat_clear);
  178.  
  179. %% 4.5. ICA
  180. % in order to get rid of blink artefacts and external noise, the
  181. % independent component analysis can be useful. But Beware, don't exclude
  182. % to much! Whereas teh ICA can split the signal into independent
  183. % components, this does not mean that one component does not contain
  184. % information from more than one source!
  185.  
  186. %% 4.5.1. Compute ICA-Components
  187. cfg=[];
  188. cfg.channel={'EEG' 'CB1' 'CB2'}; % Super important to only use the EEG-Channels, otherwise the ICA won't work
  189. cfg.method='fastica'; % Whcih method should be used?
  190. cfg.numcomponent=50; % if the ICA does not converge, limit the number of components to compute
  191. cfg.demean='no';
  192. comp=ft_componentanalysis(cfg,dat_fix);
  193.  
  194. %% 4.5.2. Take a look at the components
  195.  
  196. cfg=[];
  197. cfg.layout=lay;
  198. cfg.viewmode='component'; % same as above, just with a different mode
  199. ft_databrowser(cfg,comp);
  200.  
  201. %% 4.5.3. And take out the ones that are clearly artefacts
  202.  
  203. cfg=[];
  204. cfg.component=[9]; % I focus on blinks and muscle artefacts
  205. dat_ica=ft_rejectcomponent(cfg,comp,dat_fix);
  206.  
  207.  
  208. %% 4.5.4 Compare before & after
  209. figure;plot(dat_ica.time{1},dat_ica.trial{1}(22,:),'r');hold;plot(dat_clear.time{1},dat_clear.trial{1}(22,:),'b')
  210.  
  211. %% 6 Save
  212. %mkdir 2ndtry
  213. cd 2ndtry
  214. save dat_orig.mat dat -V7.3
  215. %save dat_cubic.mat dat_cubic -V7.3
  216. save dat_clear.mat dat_clear -V7.3
  217. %save dat_filt.mat dat_filt -V7.3
  218. save dat_res.mat dat_res -V7.3
  219. save dat_ica.mat dat_ica comp -V7.3
  220. save dat_fix.mat dat_fix -V7.3
RAW Paste Data