= 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
}}}