Maxfilter Script in Matlab
% 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/lanaction/... (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/lanaction/... (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/lanaction/... (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