Pre-Processing MEG Data in SPM8
This script will epoch and average your MEG data. Filtering is optional.
The "chanfile" is a Matlab mat-file containing an array "label" that contains labels for each channel in the data file, i.e. e.g. label{1} = 'MEG0113', label{2} = 'MEG0112', ..., label{307} = 'EEG001', ..., label{380} = 'STI101'.
% Root OUTPUT DIRECTORY (e.g. '/imaging/meg.ryan/data')
bwd = '/MyDataDir/'
% File Containing CHANNEL INFORMATION for MEG and EEG (e.g.
% '/imaging/meg.ryan/Batch/chan_select_MEG_EEG_STI101.mat')
chanfile = 'MyChanFile';
% Specify (multiple) SUBJECTs)
cnt = 0;
cnt = cnt + 1;
subjects{cnt} = {'meg10_1000', '101224'}; % IDs e.g. as in /megdata/cbu/...
blocksin{cnt} = {'block1_raw_sss', 'block2_raw_sss'}; % as named after Maxfilter
cnt = cnt+1;
subjects{cnt} = {'meg10_1001', '101231'}; % IDs e.g. as in /megdata/cbu/...
blocksin{cnt} = {'block1_raw_sss', 'block2_raw_sss'}; % as named after Maxfilter
eog_thr = [100e-6]; % EOG artefact THRESHOLD (in Volts)
epoch = [-100 300]; % EPOCH for averaging (milliseconds)
% Define EVENT INFORMATION
con_values = [1:3]; % Trigger codes of interest (integers)
con_labels = {'CondA', 'CondB', 'Control'}; % Labels of conditions corresponding to trigger codes
offset = 0; % OFFSET between trigger and stimulus presentation, e.g. projector delay (milliseconds)
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% The rest should run smoothly... %
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
cd(bwd) % change to working directory
nr_sbjs = length(subjects); % number of subjects
Ncons = length(con_values);
for c=1:Ncons
con_trigs{c} = 'STI101_up'; % look for rising signal in trigger channel for event times
end
efile = {};
for ss = 1:nr_sbjs,
nr_sess = length( blocksin{ss} );
clear S D;
swd = fullfile(bwd,subjects{ss}{1},subjects{ss}{2});
fprintf(1, 'Subject: %s\n', swd);
try
eval(sprintf('!mkdir %s',swd))
end
cd(swd)
for ses = 1:nr_sess
rawfile = fullfile(bwd,subjects{ss}{1},subjects{ss}{2},sprintf('%s.fif',blocksin{ss}{ses}));
fprintf(1, 'Processing %s\n', rawfile);
S = [];
S.dataset = rawfile;
S.outfile = sprintf('spm8_%s',spm_str_manip(rawfile,'rt')); %_%s',filestem,maxflag)
tmp = load(chanfile);
S.channels = tmp.label;
D = spm_eeg_convert(S);
% spm_eeg_convert_ui(S); % if want to use GUI
% Below are other defaults
% S.blocksize = 3276800;
% S.datatype = 'float32-le';
% S.timewindow = [1 5];
% S.inputformat = [];
% S.eventpadding = 0;
% S.saveorigheader = 0;
% S.continuous = 1;
% S.checkboundary = 1;
% S.usetrials = 1; % SPM will define trials from datafile
ch = find(strcmp(D.chanlabels,'STI101'));
D(ch,:,:) = D(ch,:,:)/10^7; % Downscale trigger channel so easier to see EOG when display "Other" channels
% If you want to filter data in SPM8:
% S = [];
% S.D = D.fname;
% S.filter.type = 'butterworth';
% S.filter.order = 5;
% S.filter.band = 'bandpass';
% S.filter.PHz = [0.1 40];
% D = spm_eeg_filter(S);
S = [];
S.dataset = D.fname;
S.pretrig = epoch(1);
S.posttrig = epoch(2);
S.save = 0; % saved anyway (if S.save=1, then prompts for new filename!)
% S.reviewtrials = 1; % enable if you want to check before converting
S.reviewtrials = 0;
for c = 1:length(con_values)
S.trialdef(c).conditionlabel = con_labels{c};
S.trialdef(c).eventtype = 'STI101_up';
S.trialdef(c).eventvalue = con_values(c);
S.trialdef(c).trlshift = offset(c);
% If used binary coding (1,2,4,8...) then could read STI001, STI002, etc
% S.trialdef(c).eventtype = sprintf('STI00%d_up',c);
% S.trialdef(c).eventvalue = 5;
end
[trl, con, S] = spm_eeg_definetrial(S);
S = [];
S.D = D.fname;
S.epochinfo.trl = trl;
S.epochinfo.conditionlabels = con;
S.bc = 1;
D = spm_eeg_epochs(S);
efile{ses} = D.fname;
end % of ses loop
%%% Concatenation of sessions
S=[];
S.D = strvcat(efile);
S.recode = 'same';
D = spm_eeg_merge(S);
%%% Artifact rejection
S = []; S.D = D.fname;
S.methods(1).fun = 'flat';
S.methods(1).channels = 'MEG';
S.methods(1).settings.threshold = 0;
S.methods(1).settings.seqlength = 4;
S.methods(end+1).fun = 'flat';
S.methods(end).channels = 'EEG';
S.methods(end).settings.threshold = 0;
S.methods(end).settings.seqlength = 4;
% S.methods(end+1).fun = 'peak2peak';
S.methods(end+1).fun = 'threshchan';
S.methods(end).channels = 'EOG';
S.methods(end).settings.threshold = eog_thr(ss);
% S.methods(end+1).fun = 'threshchan';
% S.methods(end).channels = 'MEG';
% S.methods(end).settings.threshold = meg_thr(ss);
%
% S.methods(end+1).fun = 'threshchan';
% S.methods(end).channels = 'EEG';
% S.methods(end).settings.threshold = eeg_thr(ss);
D = spm_eeg_artefact(S);
nbadchan(ss) = length(D.badchannels);
nrejects(ss) = sum(D.reject);
%%% Average
D = condlist(D,con_labels); % redefine condition order for weight epochs below
D.save;
preavg = D.fname;
S=[]; S.D = preavg;
S.robust = 0;
D = spm_eeg_average(S);
nevents(ss,:) = D.repl;
S=[]; S.D = D.fname;
S.refchan = 'average';
D = spm_eeg_reref_eeg(S);
save batch_params rawfile efile nbadchan nrejects nevents
end % of subjects loop
return