<?xml version="1.0" encoding="utf-8"?><!DOCTYPE article  PUBLIC '-//OASIS//DTD DocBook XML V4.4//EN'  'http://www.docbook.org/xml/4.4/docbookx.dtd'><article><articleinfo><title>MaxfilterMatlabScript</title><revhistory><revision><revnumber>10</revnumber><date>2013-03-08 10:02:31</date><authorinitials>localhost</authorinitials><revremark>converted to 1.6 markup</revremark></revision><revision><revnumber>9</revnumber><date>2013-02-04 14:48:22</date><authorinitials>OlafHauk</authorinitials></revision><revision><revnumber>8</revnumber><date>2011-04-28 15:47:02</date><authorinitials>YaaraErez</authorinitials></revision><revision><revnumber>7</revnumber><date>2011-03-23 16:58:59</date><authorinitials>YaaraErez</authorinitials></revision><revision><revnumber>6</revnumber><date>2011-02-24 12:48:14</date><authorinitials>YaaraErez</authorinitials></revision><revision><revnumber>5</revnumber><date>2011-02-24 12:40:40</date><authorinitials>YaaraErez</authorinitials></revision><revision><revnumber>4</revnumber><date>2010-12-09 12:08:09</date><authorinitials>YaaraErez</authorinitials></revision><revision><revnumber>3</revnumber><date>2010-12-09 12:06:31</date><authorinitials>YaaraErez</authorinitials></revision><revision><revnumber>2</revnumber><date>2010-07-06 17:10:57</date><authorinitials>YaaraErez</authorinitials></revision><revision><revnumber>1</revnumber><date>2010-07-06 17:03:43</date><authorinitials>YaaraErez</authorinitials></revision></revhistory></articleinfo><para><inlinemediaobject><imageobject><imagedata fileref="https://imaging.mrc-cbu.cam.ac.uk/meg/MaxfilterMatlabScript?action=AttachFile&amp;do=get&amp;target=mrclogo.gif"/></imageobject><textobject><phrase>mrclogo.gif</phrase></textobject></inlinemediaobject> </para><section><title>Maxfilter Script in Matlab</title><para>This script will do the following for you: </para><itemizedlist><listitem><para>Detect bad MEG channels from the pre-HPI period in your data (assuming HPI measurement was switched on after 20s) </para></listitem><listitem><para>Apply SSS including ST and movement compensation, downsampling by a factor 3 (to 3 ms, if sampling frequency is 1000 Hz), assuming head origin [0 0 45] for all data sets </para></listitem><listitem><para>Interpolate each data set to the first one specified in the list, for each subject separately (&quot;trans&quot;) </para></listitem></itemizedlist><para>The script automatically detects bad channels from the pre-HPI period, but also allows you to specify further bad MEG channels at the top. </para><para>At the end, you will have files ending in &quot;sss&quot; (before trans) and &quot;ssst&quot; (after trans), which you can use for the interesting part of your analysis... </para><screen><![CDATA[% based on script by Jason Taylor
clear subject blocksin blocksout badchannels;
]]><![CDATA[
pathstem = '/YourOutputPath/'; % for output data
rawpathstem = '/megdata/cbu/YourSubDir';  % input data
]]><![CDATA[
% Define data for individual subjects as follows:
cnt = 1;
subject{cnt} = {'meg01_0001', '012345'};
blocksin{cnt} = {'block1', 'block2', 'block3', 'block4'};   % as named during recording, in /megdata/cbu/YourSubDir/... (may differ across subjects)
blocksout{cnt} = {'block1', 'block2', 'block3', 'block4'};      % should be consistent for all subjects
badchannels{cnt, 1} = {'0741', '1533'}; badchannels{cnt, 2} = {'1533', '0713'}; badchannels{cnt, 3} = {''}; badchannels{cnt, 4} = {''}; % define bad MEG (not EEG) channels here (if there are any)
cnt=cnt+1;
]]><![CDATA[
subject{cnt} = {'meg02_0002', '123456'};
blocksin{cnt} = {'block1', 'block2', 'block3', 'block4'};   % as named during recording, in /megdata/cbu/YourSubDir/... (may differ across subjects)
blocksout{cnt} = {'block1', 'block2', 'block3', 'block4'};      % should be consistent for all subjects
badchannels{cnt, 1} = {'0741', '1533'}; badchannels{cnt, 2} = {'1533', '0713'}; badchannels{cnt, 3} = {''}; badchannels{cnt, 4} = {''}; % define bad MEG (not EEG) channels here (if there are any)
cnt=cnt+1;
]]><![CDATA[
subject{cnt} = {'meg03_0003', '234557'};
blocksin{cnt} = {'block1', 'block2', 'block3', 'block4'};   % as named during recording, in /megdata/cbu/YourSubDir/... (may differ across subjects)
blocksout{cnt} = {'block1', 'block2', 'block3', 'block4'};      % should be consistent for all subjects
badchannels{cnt, 1} = {'0741', '1533'}; badchannels{cnt, 2} = {'1533', '0713'}; badchannels{cnt, 3} = {''}; badchannels{cnt, 4} = {''}; % define bad MEG (not EEG) channels here (if there are any)
cnt=cnt+1;
]]><![CDATA[
%% The rest should not need editing
]]><![CDATA[
nr_sbj = length(subject);
]]><![CDATA[
try do_subjects,    % if do_subjects not defined, do all subjects
catch
    do_subjects = [1:nr_sbj];
end;
]]><![CDATA[
% Check file names and paths
checkflag = 0;
for ss = do_subjects,
    nr_bls = length( blocksin{ss} );
    if length(blocksin{ss}) ~= length(blocksout{ss}),
        checkflag = 1;
        fprintf(1, 'Different number of input and output names for subject %d (%s, %s)\n', ss, subject{ss}{1}, subject{ss}{2});
    end;
    for bb = 1:nr_bls,        
        rawpath = fullfile( rawpathstem, subject{ss}{1}, subject{ss}{2} );
        rawfname = fullfile( rawpath, [blocksin{ss}{bb} '_raw.fif'] );
        outpath = fullfile( pathstem, subject{ss}{1}, subject{ss}{2} );
        if ~exist( outpath, 'dir' ),
            success = mkdir( outpath );
            if ~success,
                checkflag = 1;
                fprintf(1, 'Could not create directory %s\n', outpath);
            end;
        end;
        if ~exist( rawfname, 'file' ),
            checkflag = 1;
            fprintf(1, '%s does not exist\n', rawfname);
        end;        
    end;
end;
if checkflag,
    fprintf(1, 'You''ve got some explaining to do.\n');
    return;
end;
]]><![CDATA[
]]><![CDATA[
for ss = do_subjects,
    nr_bls = length( blocksin{ss} );
    
    for bb = 1:nr_bls,
 
        rawpath = fullfile( rawpathstem, subject{ss}{1}, subject{ss}{2} );
        rawfname = fullfile( rawpath, [blocksin{ss}{bb} '_raw.fif'] );
        
        outpath = fullfile( pathstem, subject{ss}{1}, subject{ss}{2} );
        
        outfname1 = fullfile( outpath, [blocksout{ss}{bb} '_raw_tmp.fif'] );    % files after bad channel check
        logfname1 = fullfile( outpath, [blocksout{ss}{bb} '_raw_tmp.log'] );
]]><![CDATA[
        outfname2 = fullfile( outpath, [blocksout{ss}{bb} '_raw_sss.fif'] );    % files after SSS+ST
        logfname2 = fullfile( outpath, [blocksout{ss}{bb} '_raw_sss.log'] );
        
        outfname3 = fullfile( outpath, [blocksout{ss}{bb} '_raw_ssst.fif'] );   % files after interpolation to first specified session
        logfname3 = fullfile( outpath, [blocksout{ss}{bb} '_raw_ssst.log'] );                
]]><![CDATA[
        posfname = fullfile( outpath, [blocksout{ss}{bb} '_raw_hpi.pos'] );     % HPI info
]]><![CDATA[
        badfname = fullfile( outpath, [blocksout{ss}{bb} '_raw_bad.txt'] );     % bad channel info
]]><![CDATA[
        markbadfname = fullfile( outpath, [blocksout{ss}{bb} '_raw_markbad.fif'] );
        
        fprintf(1, '\n Now processing %s with %d pre-specified bad channels.\n', rawfname, length( badchannels{ss, bb} ) );
]]><![CDATA[
%% (2) Convert data 
]]><![CDATA[
        skipint = '0 20';
        mfcmd2=[
        '/neuro/bin/util/maxfilter -f ' [rawfname] ' -o ' [outfname1],...
        ' -autobad 20 -skip ' [skipint] ' -v | tee ' [logfname1]
        ];
        fprintf(1, '\n\n%s\n\n', mfcmd2);
        eval([' ! ' mfcmd2])
        delete( outfname1 );
]]><![CDATA[
%% Get bad channels
        % Get bad channels from log file, store in file:
        badcmd=[
        'cat ' [logfname1] ' | sed -n ''/Static/p'' | cut -f 5- -d '' '' > ' [badfname]
        ];
        fprintf(1, 'Looking for bad channels\n');
        fprintf(1, '\n%s\n', badcmd);
        eval([' ! ' badcmd]);
]]><![CDATA[
]]><![CDATA[
        % Read bad channels in to matlab variable:
        fprintf(1, '\nReading bad channel information\n');
        x=dlmread([badfname],' ');
        x=reshape(x,1,prod(size(x)));
        x=x(x>0); % Omit zeros (padded by dlmread):
]]><![CDATA[
]]><![CDATA[
        % Get frequencies (number of buffers in which chan was bad):
        [frq,allbad] = hist(x,unique(x));
]]><![CDATA[
]]><![CDATA[
        % Mark bad based on threshold (currently 5 buffers):
        bads=allbad(frq>5);
        badstxt = sprintf('%s%s%s',num2str(bads))
        if sum(badstxt)>0
            dlmwrite([markbadfname],badstxt,'delimiter',' ');
        else
            eval(['! touch ' [markbadfname] ])
        end
]]><![CDATA[
        % If extra bad channels defined, append them here
        if ~isempty( badchannels{ss,bb} ),
            for i=1:length(badchannels{ss,bb}),
                badstxt = [badstxt ' ' badchannels{ss,bb}{i}];
            end;
        end;
        fprintf(1, '\nThe following channels are marked as bad: %s\n\n', badstxt);
]]><![CDATA[
%% (3) Maxfilter incl. ST and Movecomp
        % -- MAXFILTER ARGUMENTS --:
]]><![CDATA[
        % ORIGIN and FRAME:
        orgcmd=sprintf('  -frame head -origin 0 0 45');
]]><![CDATA[
]]><![CDATA[
        % BAD CHANNELS:
        if length(badstxt)>0
            badcmd=['  -bad ', badstxt];
        else
            badcmd='';
        end
]]><![CDATA[
]]><![CDATA[
        % HPI ESTIMATION/MOVEMENT COMPENSATION:
        hpistep=200;hpisubt='amp';
        hpicmd=sprintf('  -hpistep %d -hpisubt %s -movecomp -hp %s',hpistep,hpisubt,posfname);
]]><![CDATA[
        % SSS with ST:
        stwin=4;
        stcorr=0.980;
        stcmd=sprintf('  -st %d -corr %g',stwin,stcorr);
]]><![CDATA[
        % Downsampling
        dsval = 3;
        dscmd=sprintf('  -ds %d', dsval');
]]><![CDATA[
]]><![CDATA[
        % -- MAXFILTER COMMAND --
]]><![CDATA[
        if exist(outfname2),
            fprintf(1, 'Deleting %s\n', outfname2);
            delete( outfname2 );
        end;
    
        mfcmd3=[
        ' /neuro/bin/util/maxfilter -f ' [rawfname] ' -o ' [outfname2],...
        '  -ctc /neuro/databases/ctc/ct_sparse.fif' ' ',...
        '  -cal /neuro/databases/sss/sss_cal.dat' ' ',...
        '  -autobad off ',...
        '  -skip 0 20 ',...
        stcmd,...       % temporal SSS
        dscmd,...       % downsampling
        badcmd,...      % bad channels
        orgcmd,...      % head frame and origin
        hpicmd,...      % movement compensation
        '  -format short ',...
        '  -v | tee ' [logfname2]
        ];
]]><![CDATA[
]]><![CDATA[
        fprintf(1, '\nMaxfiltering... (SSS+ST)\n');
        fprintf(1, '\n\n%s\n\n', mfcmd3);
        eval([' ! ' mfcmd3 ]);
]]><![CDATA[
]]><![CDATA[
        % (4) %%%%%%%%%%%%%%%%%%%%%%%%%
]]><![CDATA[
]]><![CDATA[
        % TRANSFORMATION (all but first file, block 1):
        if bb>1
        
            trcmd=sprintf(' -trans %s -frame head -origin 0 0 45',b1file);
]]><![CDATA[
            mfcmd4=[
            '/neuro/bin/util/maxfilter -f ' [outfname2] ' -o ' [outfname3],...
            '  -autobad off ', trcmd, ' -force -v | tee ' logfname3
            ];
            fprintf(1, '\nMaxfiltering... -trans\n');
            fprintf(1, '%s\n', mfcmd4);
            eval([' ! ' mfcmd4 ])
            
        else,
            
            b1file = outfname2;     % file used for future "trans"
            copyfile( outfname2, outfname3 );
            
        end;    % if bb>1        
]]><![CDATA[
    end;    % blocks
    
end;    % subjects]]></screen></section></article>