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
