Homec4science
Diffusion iCAPs public (master) Browse

iCAPs/packages/spm12/external/mne/mne_ex_evoked_grad_amp.m
2cb18aae5832master

Last Change
Raw File

mne_ex_evoked_grad_amp.m
View Options

function [res] = mne_ex_evoked_grad_amp(inname,bmin,bmax,outname)
%
% function [res] = mne_ex_evoked_grad_amp(inname,bmin,bmax,outname)
%
% Compute the magnitude of the tangential gradient at each
% sensor location using the planar gradiometer data and
% optionally output the result to a fif file.
%
% inname The input file name. All average data sets are
% read and processed
% bmin,bmax Baseline limits in seconds
% outname Optional output file name
%
%
% Function returns the data which was or would have been written
% to the file
%
%
% Author : Matti Hamalainen, MGH Martinos Center
% License : BSD 3-clause
%
me='MNE:mne_ex_evoked_grad_amp';
global FIFF;
if isempty(FIFF)
FIFF = fiff_define_constants();
end
if nargin == 1
do_baseline = false;
elseif (nargin == 3 || nargin == 4)
do_baseline = true;
else
error(me,'Wrong number of arguments');
end
%
% Read the data
%
data = fiff_read_evoked_all(inname);
%
% Figure out the planar gradiometer pairs
%
pairs = zeros(data.info.nchan,2);
npair = 0;
k = 1;
while k < data.info.nchan
%
% First check the coil types
%
coil1 = data.info.chs(k).coil_type;
coil2 = data.info.chs(k+1).coil_type;
if (coil1 == coil2 && ...
(coil1 == 2 || coil1 == 3012 || coil1 == 3013))
one = data.info.ch_names{k};
two = data.info.ch_names{k+1};
lastone = one(length(one));
lasttwo = two(length(two));
%
% Then the channel names
%
if (strcmp(one(1:3),'MEG') && strcmp(one(1:3),'MEG'))
if (strcmp(one(1:(length(one)-1)),two(1:(length(two)-1))) && ...
((lastone == '2' && lasttwo == '3') || ...
(lastone == '3' && lasttwo == '2')))
npair = npair + 1;
pairs(npair,1) = k;
pairs(npair,2) = k+1;
k = k + 1;
end
end
end
k = k + 1;
end
if npair == 0
error(me,'No planar gradiometers in these data');
end
%
% Compute the amplitudes
%
fprintf(1,'Computing the amplitudes');
if do_baseline
fprintf(1,' (Baseline = %7.1f ... %7.1f ms)',1000*bmin,1000*bmax);
end
fprintf(1,'...');
for k = 1:length(data.evoked)
epochs = data.evoked(k).epochs;
%
% Setup baseline limits
%
if do_baseline
b1 = double(data.info.sfreq*bmin - data.evoked(k).first);
b2 = double(data.info.sfreq*bmax - data.evoked(k).first);
if b1 < 1
b1 = 1;
end
if b2 > size(epochs,2)
b2 = size(epochs,2)
end
else
b1 = 1;
b2 = 1;
end
%
% Go through all pairs
%
for p = 1:npair
one = pairs(p,1);
two = pairs(p,2);
if b2 > b1
base1 = sum(epochs(one,b1:b2))/(b2-b1);
base2 = sum(epochs(two,b1:b2))/(b2-b1);
epochs(one,:) = sqrt((epochs(one,:)-base1).*(epochs(one, ...
:)-base1)+(epochs(two,:)-base2).*(epochs(two,:)-base2));
else
epochs(one,:) = sqrt(epochs(one,:).*epochs(one, ...
:)+epochs(two,:).*epochs(two,:));
end
end
data.evoked(k).epochs = epochs;
fprintf(1,'.');
end
fprintf(1,'[done]\n');
%
% Compose the selection name list
%
pairs = pairs(1:npair,:);
for k = 1:npair
ch_sel_names{k} = data.info.ch_names{pairs(k,1)};
end
%
% Omit MEG channels but include others
%
for p = 1:data.info.nchan
if (data.info.chs(p).kind ~= FIFF.FIFFV_MEG_CH)
k = k + 1;
ch_sel_names{k} = data.info.ch_names{p};
end
end
%
% Modify the bad channel list
%
if ~isempty(data.info.bads)
nbad = length(data.info.bads);
for k = 1:npair
one = data.info.ch_names{pairs(k,1)};
two = data.info.ch_names{pairs(k,2)};
%
% If one channel of the planar gradiometer is marked bad,
% add the other to the bad channel list
%
if (~isempty(strmatch(one,data.info.bads)) && ...
isempty(strmatch(two,data.info.bads)))
nbad = nbad + 1;
data.info.bads{nbad} = two;
elseif (isempty(strmatch(one,data.info.bads)) && ...
~isempty(strmatch(two,data.info.bads)))
nbad = nbad + 1;
data.info.bads{nbad} = one;
end
end
end
%
% Do the picking
%
res = fiff_pick_channels_evoked(data,ch_sel_names);
%
% Optionally write an output file
%
if nargin == 4
fiff_write_evoked(outname,res);
fprintf(1,'Wrote %s\n',outname);
end
c4science · Help