Maxfilter Script in Matlab

This script will do the following for you:

• Detect bad MEG channels from the pre-HPI period in your data (assuming HPI measurement was switched on after 20s)
• 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
• Interpolate each data set to the first one specified in the list, for each subject separately ("trans")

At the end, you will have files ending in "sss" (before trans) and "ssst" (after trans), which you can use for the interesting part of your analysis...

% based on script by Jason Taylor

pathstem = '/YourOutputPath/'; % for output data
rawpathstem = '/megdata/cbu/YourSubDir';  % input data

% 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} = {''}; % define bad MEG (not EEG) channels here (if there are any)
cnt=cnt+1;

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} = {''}; % define bad MEG (not EEG) channels here (if there are any)
cnt=cnt+1;

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} = {''}; % define bad MEG (not EEG) channels here (if there are any)
cnt=cnt+1;

%% The rest should not need editing

nr_sbj = length(subject);

try do_subjects,    % if do_subjects not defined, do all subjects
catch
    do_subjects = [1:nr_sbj];
end;

% 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;


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'] );

        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'] );                

        posfname = fullfile( outpath, [blocksout{ss}{bb} '_raw_hpi.pos'] );     % HPI info

        badfname = fullfile( outpath, [blocksout{ss}{bb} '_raw_bad.txt'] );     % bad channel info

        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} ) );

%% (2) Convert data 

        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 );

%% 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]);


        % 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):


        % Get frequencies (number of buffers in which chan was bad):
        [frq,allbad] = hist(x,unique(x));


        % 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

        % 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);

%% (3) Maxfilter incl. ST and Movecomp
        % -- MAXFILTER ARGUMENTS --:

        % ORIGIN and FRAME:
        orgcmd=sprintf('  -frame head -origin 0 0 45');


        % BAD CHANNELS:
        if length(badstxt)>0
            badcmd=['  -bad ', badstxt];
        else
            badcmd='';
        end


        % HPI ESTIMATION/MOVEMENT COMPENSATION:
        hpistep=200;hpisubt='amp';
        hpicmd=sprintf('  -hpistep %d -hpisubt %s -movecomp -hp %s',hpistep,hpisubt,posfname);

        % SSS with ST:
        stwin=4;
        stcorr=0.980;
        stcmd=sprintf('  -st %d -corr %g',stwin,stcorr);

        % Downsampling
        dsval = 3;
        dscmd=sprintf('  -ds %d', dsval');


        % -- MAXFILTER COMMAND --

        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]
        ];


        fprintf(1, '\nMaxfiltering... (SSS+ST)\n');
        fprintf(1, '\n\n%s\n\n', mfcmd3);
        eval([' ! ' mfcmd3 ]);


        % (4) %%%%%%%%%%%%%%%%%%%%%%%%%


        % TRANSFORMATION (all but first file, block 1):
        if bb>1
        
            trcmd=sprintf(' -trans %s -frame head -origin 0 0 45',b1file);

            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        

    end;    % blocks
    
end;    % subjects