juliankeil

FieldTrip Preprocessing

Aug 8th, 2012
1,834
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

Adblocker detected! Please consider disabling it...

We've detected AdBlock Plus or some other adblocking software preventing Pastebin.com from fully loading.

We don't have any obnoxious sound, or popup ads, we actively block these annoying types of ads!

Please add Pastebin.com to your ad blocker whitelist or disable your adblocking software.

×