Set up general paths

Change these to appropriate paths for your operating system and setup

clear variables
restoredefaultpath
addpath /home/lau/matlab/fieldtrip-20180129_workshop_teaching/
ft_defaults

meg_path = '/home/share/workshop_source_reconstruction/data/MEG/';
mri_path = '/home/share/workshop_source_reconstruction/data/MRI/';

Set up subject specific paths

Make subject and recording specific paths (the cell array “subjects_and_dates” can be expanded if you have more subjects than one)

subjects_and_dates = ...
                    {
                        'NatMEG_0177/170424/'
                    };
                
filenames = {...
        'tactile_stim_raw_tsss_mc.fif'
        'tactile_stim_raw_tsss_mc-1.fif'
        'tactile_stim_raw_tsss_mc-2.fif'
            };
        
n_filenames = length(filenames);  

output_path = fullfile(meg_path, subjects_and_dates{1});
events = 8; % right index finger
% to do all five fingers set events as below:
% events = [1 2 4 8 16]; %% right little, ring, middle, index, thumb

Set preprocessing options

  1. Elekta data files (.fif) are split into several files due to a limitation of 2 GB for each files
  2. Event code is read in (8), but note that all can be read in (1, 2, 4, 8, 16)
  3. We anticipate the measured delay of 41 ms between the trigger and the actual stimulation
  4. MEG data and EEG data are handled separately since EEG data need to be referenced to a reference “channel” (here the average of all channels)
  5. Not much preprocessing is done to the data (we’ll do that later for the respective analyses (timelocked and time frequency representation (TFR)))
  6. The split files are appended to one another to ease data handling
split_files_MEG = cell(1, n_filenames);
split_files_EEG = split_files_MEG;

for filename_index = 1:n_filenames

    filename = filenames{filename_index};
    full_path = fullfile(meg_path, subjects_and_dates{1}, filename);
 
    % define trials
    
    cfg = [];
    cfg.dataset = full_path;
    cfg.trialdef.prestim = 1.609; % seconds % adjusted with 41 ms
    cfg.trialdef.poststim = 1.691; % seconds
    cfg.trialdef.eventvalue = events;
    cfg.trialdef.eventtype = 'STI101';
    
    cfg = ft_definetrial(cfg);
    
    % preprocess
 
    cfg.demean = 'no';
    cfg.lpfilter = 'no';
    cfg.hpfilter = 'no';
    cfg.dftfilter = 'no';
    cfg.channel = 'MEG';

    split_files_MEG{filename_index} = ft_preprocessing(cfg);
    
    cfg.channel = 'EEG';
    cfg.reref = 'no';
    cfg.refchannel = 'all';
    
    split_files_EEG{filename_index} = ft_preprocessing(cfg);
         
end

% append split files

cfg = [];

preprocessed_data_MEG = ft_appenddata(cfg, split_files_MEG{:});
preprocessed_data_EEG = ft_appenddata(cfg, split_files_EEG{:});

Inspect raw data

Try to inspect the data and see if you can identify some bad electrodes (and MEG sensors)


cfg = [];
cfg.continuous = 'yes';
cfg.viewmode = 'vertical'; % also try 'butterfly'

ft_databrowser(cfg, preprocessed_data_MEG);

ft_databrowser(cfg, preprocessed_data_EEG);

Remove bad electrodes

These would mess up the source reconstructions if left in
Please note that you need not find them now, as they are specified in the next step (interpolation)


% Find bad channels (#27+34+40+79+78+96) % numbers of bad channels found
cfg = [];
cfg.method      = 'summary';
cfg.keepchannel = 'no';
cfg.channel     = 'EEG';
cfg.layout      = 'neuromag306eeg1005_natmeg.lay';

temp = ft_rejectvisual(cfg, preprocessed_data_EEG);

Interpolate bad channels and re-reference data

  1. The data for the removed channels is interpolated back in based on their neighbours
  2. The data is rereferenced to the average of all the channels
%% Interpolate bad channels
% Get the names of the selected channels
% badchannels = setdiff(preprocessed_data_EEG.label, temp.label);
badchannels = union(preprocessed_data_MEG.label, {'EEG027' 'EEG034' 'EEG040' 'EEG079' 'EEG078' 'EEG096'});

% Interpolate bad channels: Find neighbours
cfg             = [];
cfg.method      = 'triangulation';
cfg.senstype   = 'EEG'; % 
neighbours_EEG = ft_prepare_neighbours(cfg, preprocessed_data_EEG);
% plotting neighbours for inspection
cfg            = [];
cfg.neighbours = neighbours_EEG;
cfg.senstype   = 'EEG';
ft_neighbourplot(cfg, preprocessed_data_EEG);

Plot neighbours

Define the neighbours of each electrode, so bad ones can be interpolated

Do the actual interpolation

cfg = [];
cfg.method                    = 'spline';
cfg.neighbours                = neighbours_EEG;
cfg.badchannel                = badchannels;
cfg.senstype                  = 'EEG';
interpolated_data = ft_channelrepair(cfg, preprocessed_data_EEG);
Plot the interpolated and non-interplolated

Can you spot the differences?

cfg = [];
cfg.viewmode = 'butterfly';

ft_databrowser(cfg, preprocessed_data_EEG); % not interpolated
ft_databrowser(cfg, interpolated_data);

NOT interpolated