Maxfilter Script in Matlab

This script will do the following for you:

• Detect bad MEG channels from the pre-HPI period in your data (first 20 seconds)
• Apply SSS including ST and movement compensation, downsampling by a factor 3 (to 3 ms), 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
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
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
cnt=cnt+1;


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

        % (2) %%%%%%%%%%%%%%%%%%%%%%%%%

        skipint = '0 20';
        mfcmd2=[
        '/neuro/bin/util/maxfilter -f ' [rawfname] ' -o ' [outfname1],...
        ' -autobad 20 -skip ' [skipint] ' -v | tee ' [logfname1]
        ];
        fprintf(1, '\n%s\n', mfcmd2);
        eval([' ! ' mfcmd2])
        delete( outfname1 );


        % 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


        % (3) %%%%%%%%%%%%%%%%%%%%%%%%%
        % -- MAXFILTER ARGUMENTS --:


        % ORIGIN and FRAME:
        %orgcmd=sprintf('-frame head -origin %g %g %g',fit(1),fit(2),fit(3));
        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);
        hpicmd

        % 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 ',...
        '     -origin 0 0 45 ',...
        '     -frame head ',...
        '     -movecomp ',...
        '     -st',...
        '     -ds 3',...
        '     -format short ',...
        '     -hp ' posfname,...
        '  -v | tee ' [logfname2]
        ];

        fprintf(1, '\nMaxfiltering... (SSS+ST)\n');
        fprintf(1, '%s\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