Compute average sensor array

The following Matlab code will produce the average sensor array from a list of input files, and select the input file which is closest to the average (which I would recommend to use for interpolation). As input files, you can for example specify the on-line averages for your subjects.

% From a list of fiff-files (e.g. on-line average files), compute the one
% with average sensor transformation, and find the input file which is
% closest to the average (according to both translation and rotation parameters)
% file list should be specified in cell array "files{}"
% output file should be specified in string "file_out"
% both average file "_avg.fif" and file closest to it will be written
% if not output file specified, output will be written to current working directory
% OH, March 2009


% Output file for average sensor array
file_out = '/fullpath/OutputFile.fif';


% you can specify a list of files here (e.g. on-line averages for all subjects in your study)
fiff_files = {'/fullpath/file4subj1.fif', ...
              '/fullpath/file4subj2.fif', ...
              '/fullpath/file4subj3.fif'};


fid = fopen(file_out, 'a');
if fid==-1,
    fprintf(1, 'Cannot access output file %s\n', file_out);
    return;
end;
fclose(fid);

nr_files = length(fiff_files);
epoch_gm = [0];
trans_gm = [0];
for f = 1:nr_files, % read all files and average all that's relevant
    fprintf(1, '%s\n', fiff_files{f});
    data = fiff_read_evoked(fiff_files{f});
    epoch_gm = epoch_gm + data.evoked.epochs;   % average MEG data epoch
    trans_gm = trans_gm +  data.info.dev_head_t.trans;  % average sensor transformation matrix
    all_trans(:,:,f) = data.info.dev_head_t.trans;  % remember individual transformations
end;
epoch_gm = epoch_gm/nr_files;
trans_gm = trans_gm/nr_files;

data_new = data;
data_new.evoked.epoch = epoch_gm;
data_new.info.dev_head_t.trans = trans_gm;

[thispath, thisfile,thisext,thisversn] = fileparts(file_out);
file_out_avg = fullfile(thispath, [thisfile '_avg' thisext]);
fiff_write_evoked(file_out_avg, data_new);  % Write average coordinate transformation

% find input file that is closest to average...
for f = 1:nr_files,
    diff_trans_t = all_trans(1:3,4,f) - trans_gm(1:3,4);        % translation difference
    diff_trans_r = all_trans(1:3,1:3,f) - trans_gm(1:3,1:3);    % rotation difference
    norm_t(f) = norm( diff_trans_t );
    norm_r(f) = norm( diff_trans_r );
end;

[y_t, i_t] = sort( norm_t );    % sort according to translation difference
[y_r, i_r] = sort( norm_r );    % sort according to rotation difference

for f = 1:nr_files, % combine the two sorted lists
    rank_t(f) = find( (i_t-f)==0 ); % rank for translation
    rank_r(f) = find( (i_r-f)==0 ); % rank for rotation
    avg_rank(f) = (rank_t(f) + rank_r(f))/2;
    [min_val, min_file] = min( avg_rank );  % best average rank
end;

min_file = fiff_files{min_file}; % filename for file with best average rank

fprintf('\nThe winner is... %s\n', min_file);

copyfile(min_file, file_out);   % copy file with best average rank to output file

% The end of this