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/';
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
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{:});
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);
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
% 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);
Define the neighbours of each electrode, so bad ones can be interpolated
cfg = [];
cfg.method = 'spline';
cfg.neighbours = neighbours_EEG;
cfg.badchannel = badchannels;
cfg.senstype = 'EEG';
interpolated_data = ft_channelrepair(cfg, preprocessed_data_EEG);
Can you spot the differences?
cfg = [];
cfg.viewmode = 'butterfly';
ft_databrowser(cfg, preprocessed_data_EEG); % not interpolated
ft_databrowser(cfg, interpolated_data);
NOT interpolated
Interpolated
This is necessary as electric potentials can only compared by reference to a common reference
cfg = [];
cfg.reref = 'yes';
cfg.refchannel = 'all';
rereferenced_interpolated_data = ft_preprocessing(cfg, interpolated_data);
Now the data can be handled together
cfg = [];
preprocessed_data = ft_appenddata(cfg, preprocessed_data_MEG, ...
rereferenced_interpolated_data);
cfg = [];
cfg.method = 'summary';
cfg.keepchannel = 'yes';
cfg.channel = 'MEGMAG';
cfg.layout = 'neuromag306eeg_1005_natmeg.lay';
cleaned_data = ft_rejectvisual(cfg, preprocessed_data);
cfg.channel = 'MEGGRAD';
cleaned_data = ft_rejectvisual(cfg, cleaned_data);
cfg.channel = 'EEG';
cleaned_data = ft_rejectvisual(cfg, cleaned_data);
This moves the offset (0 msec) 41 msec back in time. This has been measured with an accelerometer
cfg = [];
cfg.offset = -41; % ms
cleaned_data = ft_redefinetrial(cfg, cleaned_data);
Downsample from 1000 Hz to 200 Hz
cfg = [];
cfg.resamplefs = 200;
cleaned_downsampled_data = ft_resampledata(cfg, cleaned_data);
% save(fullfile(output_path, 'cleaned_downsampled_data'), ... % uncomment to save the data in your current path
% 'cleaned_downsampled_data');
Q: Why is it necessary to include noise covariance estimation if you want to combine magnetometers, gradiometers and electrodes when doing source analsis?
n_events = length(events);
timelockeds = cell(1, n_events);
cfg = [];
cfg.toilim = [-0.200 0.600];
cropped_data = ft_redefinetrial(cfg, cleaned_downsampled_data);
for event_index = 1:n_events
event = events(event_index);
cfg = [];
cfg.covariance = 'yes';
cfg.covariancewindow = 'prestim';
cfg.preproc.demean = 'yes';
cfg.preproc.baselinewindow = [-0.200 0];
cfg.preproc.lpfilter = 'yes';
cfg.preproc.lpfreq = 70;
cfg.trials = cropped_data.trialinfo == event;
timelockeds{event_index} = ft_timelockanalysis(cfg, cropped_data);
end
%save(fullfile(output_path, 'timelockeds'), 'timelockeds'); %% uncomment to
%save on output_path
Just for making a comparison
n_events = length(events);
timelockeds_non_demeaned = cell(1, n_events);
cfg = [];
cfg.toilim = [-0.200 0.600];
cropped_data = ft_redefinetrial(cfg, cleaned_downsampled_data);
for event_index = 1:n_events
event = events(event_index);
cfg = [];
cfg.covariance = 'yes';
cfg.covariancewindow = 'prestim';
cfg.preproc.demean = 'no'; % this has been changed
cfg.preproc.baselinewindow = [-0.200 0];
cfg.preproc.lpfilter = 'yes';
cfg.preproc.lpfreq = 70;
cfg.trials = cropped_data.trialinfo == event;
timelockeds_non_demeaned{event_index} = ft_timelockanalysis(cfg, ...
cropped_data);
end
See the difference that demeaning does
Q: Which time course is demeaned and which is not? Which topography is demeaned and which not?
Planar gradiometers are combined by taking the absolute sum of each pair of channels
After that EEG, gradiometers and magnetometers are plotted separately on a 2d plot
Then they are plotted on 3d surfaces. Use the rotation tools to answer the questions below
Q: Describe the three topographies in as much detail as possible answering:
a) do the topographies point towards a single dipolar source or several distrubuted sources?
b) is/are the underlying source(s) most likely to be radial or tangential dipoles?
c) what is are the position(s) of the underlying source(s)? (Indicate on an image with LibreOffice/Microsoft Office/etc.)
d) what is/are the orientation(s) of the source(s)? (Indicate on an image with LibreOffice/Microsoft Office/etc.)
e) on one of the plots it is not possible to get any information about the orientation of the source. Which one is that?
f) what are the channels showing the maximum/minimum responses?
g) what are the maximum/minimum values approximately for each of the plots?
%% combine gradiometers and subsequently plot
close all
cfg = [];
cfg.method = 'sum';
combined_planar = ft_combineplanar(cfg, timelockeds{1});
% plot eeg
figure('units', 'normalized', 'outerposition', [0 0 1 1]);
cfg = [];
cfg.layout = 'neuromag306eeg1005_natmeg.lay';
cfg.channel = 'EEG';
cfg.xlim = [0.057 0.057];
ft_topoplotER(cfg, timelockeds{1});
% plot gradiometers
figure('units', 'normalized', 'outerposition', [0 0 1 1]);
cfg = [];
cfg.layout = 'neuromag306cmb.lay';
cfg.channel = 'MEGGRAD';
cfg.xlim = [0.057 0.057];
ft_topoplotER(cfg, combined_planar);
% plot magnetometers
figure('units', 'normalized', 'outerposition', [0 0 1 1]);
cfg = [];
cfg.layout = 'neuromag306mag.lay';
cfg.channel = 'MEGMAG';
cfg.xlim = [0.057 0.057];
ft_topoplotER(cfg, timelockeds{1});
%% 3d topo plots
close all
elec_pos = combined_planar.elec.chanpos;
sensors = combined_planar.grad;
mag_indices = 103:204;
grad_indices = 1:102;
mag_pos = sensors.chanpos(mag_indices, :);
grad_pos = sensors.chanpos(grad_indices, :);
time_index = 52;
figure('units', 'normalized', 'outerposition', [0 0 1 1])
hold on
ft_plot_topo3d(elec_pos, combined_planar.avg(205:end, time_index));
ft_plot_sens(combined_planar.elec, 'elecsize', 10, 'label', 'label')
c = colorbar;
c.Label.String = 'Electric Potential (V)';
title('EEG 57 msec', 'fontsize', 40);
figure('units', 'normalized', 'outerposition', [0 0 1 1])
hold on
ft_plot_topo3d(mag_pos, combined_planar.avg(mag_indices, time_index));
ft_plot_sens(combined_planar.grad, 'label', 'label', 'chantype', 'megmag')
c = colorbar;
c.Label.String = 'Magnetic Field Strength (T)';
title('Magnetometers 57 msec', 'fontsize', 40);
figure('units', 'normalized', 'outerposition', [0 0 1 1])
hold on
ft_plot_topo3d(grad_pos, combined_planar.avg(grad_indices, time_index));
ft_plot_sens(combined_planar.grad, 'label', 'label', 'chantype', 'unknown')
c = colorbar;
c.Label.String = 'Magnetic Field Gradient (T/m)';
title('Gradiometers 57 msec', 'fontsize', 40);