Advertisement
Not a member of Pastebin yet?
Sign Up,
it unlocks many cool features!
- %% 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
Advertisement
Add Comment
Please, Sign In to add comment
Advertisement