[FieldTrip] marking artifacts by channel + trial

Teresa Madsen braingirl at gmail.com
Wed Mar 8 21:35:12 CET 2017


I actually do a mix of approaches:  a quick look using ft_databrowser with
all channels for irregular artifacts like disconnection events, then a
channel-by-channel search for large z-value artifacts and clipping
artifacts, then I remove all those and do one last ft_databrowser review of
all channels together.  I'll attach the function I was working on, but it's
more complex than you originally asked for and not fully tested yet, so use
at your own risk.

Do you use ft_databrowser or ft_rejectvisual for visual artifact rejection?

~Teresa


On Wed, Mar 8, 2017 at 12:04 PM, Tim Meehan <timeehan at gmail.com> wrote:

> Thanks for sharing! I'm just taking a look now. It looks like you're doing
> mostly automated rejection. Or are you also doing visual rejection along
> with the z-value thresholding?
>
> Thanks again,
> Tim
>
> On Fri, Mar 3, 2017 at 5:31 PM, Teresa Madsen <braingirl at gmail.com> wrote:
>
>> Here's a rough sketch of my approach, with one custom function attached.
>> If you or others find it useful, maybe we can think about ways to
>> incorporate it into the FieldTrip code.  I've been working mostly with
>> scripts, but you've inspired me to work on functionizing the rest of it so
>> it's more shareable.
>>
>> So, assuming raw multichannel data has been loaded into FieldTrip
>> structure 'data' with unique trial identifiers in data.trialinfo...
>>
>> for ch = 1:numel(data.label)
>>   %% pull out one channel at a time
>>   cfg = [];
>>   cfg.channel = data.label{ch};
>>
>>   datch{ch} = ft_selectdata(cfg, data);
>>
>>   %% identify large z-value artifacts and/or whatever else you might want
>>
>>   cfg = [];
>>   cfg.artfctdef.zvalue.channel        = 'all';
>>   cfg.artfctdef.zvalue.cutoff         = 15;
>>   cfg.artfctdef.zvalue.trlpadding     = 0;
>>   cfg.artfctdef.zvalue.fltpadding     = 0;
>>   cfg.artfctdef.zvalue.artpadding     = 0.1;
>>   cfg.artfctdef.zvalue.rectify        = 'yes';
>>
>>   [~, artifact.zvalue] = ft_artifact_zvalue([], datch{ch});
>>
>>   %% replace artifacts with NaNs
>>   cfg = [];
>>   cfg.artfctdef.zvalue.artifact = artifact.zvalue;
>>   cfg.artfctdef.reject          = 'nan';
>>
>>   datch{ch} = ft_rejectartifact(cfg,datch{ch});
>> end
>>
>> %% re-merge channels
>> data = ft_appenddata([],datch);
>>
>> %% mark uniform NaNs as artifacts when they occur across all channels
>> % and replace non-uniform NaNs (on some but not all channels) with
>> zeroes, saving times
>> [artifact,data,times] = artifact_nan2zero_TEM(data) % custom function,
>> see attached
>>
>> %% reject artifacts by breaking into sub-trials
>> cfg = [];
>> cfg.artfctdef.nan2zero.artifact = artifact;
>> cfg.artfctdef.reject            = 'partial';
>>
>> data = ft_rejectartifact(cfg,data);
>>
>> %% identify real trials
>> trlinfo = unique(data.trialinfo,'rows','stable');
>>
>> for tr = 1:size(trlinfo,1)
>>
>>   %% calculate trial spectrogram
>>
>>   cfg = [];
>>
>>   cfg.trials = ismember(data.trialinfo, trlinfo(tr,:), 'rows');
>>   cfg.keeptrials = 'no';             % refers to sub-trials
>>
>>   cfg.method     = 'mtmconvol';
>>
>>   cfg.output     = 'powandcsd';
>>
>>   cfg.foi = 2.^(0:0.1:log2(300));    % 83 freqs, log2 spaced, 1-300 Hz
>>   cfg.tapsmofrq  = cfg.foi/10;       % smooth by 10%
>>   cfg.t_ftimwin  = 2./cfg.tapsmofrq; % for 3 tapers (K=3), T=2/W
>>   cfg.toi        = '50%';
>>   cfg.pad        = 'nextpow2';
>>
>>
>>   freq = ft_freqanalysis(cfg,data);
>>
>>   %% replace powspctrm & crsspctrum values with NaNs
>>   % where t_ftimwin (or wavlen for wavelets) overlaps with artifact
>>   for ch = 1:numel(freq.label)
>>     badt = [times{tr,ch}];
>>     if ~isempty(badt) && any(...
>>         badt > (min(freq.time) - max(freq.cfg.t_ftimwin)) & ...
>>         badt < (max(freq.time) + max(freq.cfg.t_ftimwin)))
>>       ci = find(any(strcmp(freq.label{ch}, freq.labelcmb)));
>>       for t = 1:numel(freq.time)
>>         for f = 1:numel(freq.freq)
>>           mint = freq.time(t) - freq.cfg.t_ftimwin(f);
>>           maxt = freq.time(t) + freq.cfg.t_ftimwin(f);
>>           if any(badt > mint & badt < maxt)
>>             freq.powspctrm(ch,f,t) = NaN;
>>             freq.crsspctrm(ci,f,t) = NaN;
>>           end
>>         end
>>       end
>>     end
>>   end
>>
>>   %% save corrected output
>>
>>   save(['trial' num2str(tr) 'mtmconvolTFA.mat'], 'freq', '-v7.3');
>> end
>>
>>
>>
>> On Thu, Mar 2, 2017 at 9:55 AM, Tim Meehan <timeehan at gmail.com> wrote:
>>
>>> Hi Teresa,
>>>
>>> Thanks for the reply. I'll take a look at your example if you don't mind
>>> sharing. Thanks!
>>>
>>> Tim
>>>
>>> On Thu, Mar 2, 2017 at 9:53 AM, Teresa Madsen <braingirl at gmail.com>
>>> wrote:
>>>
>>>> No, not really.  The only way I've found to do that is to loop through
>>>> my artifact rejection process on each trial individually, then merge them
>>>> back together with NaNs filling in where there are artifacts, but then that
>>>> breaks every form of analysis I want to do.  :-P
>>>>
>>>> I wonder if it would work to fill in the artifacts with 0s instead of
>>>> NaNs....I might play with that.  Let me know if you're interested in some
>>>> example code.
>>>>
>>>> ~Teresa
>>>>
>>>>
>>>> On Wed, Mar 1, 2017 at 3:55 PM, Tim Meehan <timeehan at gmail.com> wrote:
>>>>
>>>>> Hello All,
>>>>>
>>>>> When performing visual artifact rejection, I want to be able to mark
>>>>> artifacts that occur during some specific trials and only on some specific
>>>>> channels. In the tutorials I see only ways to mark bad channels (i.e.
>>>>> across all trials) or bad trials (i.e. across all channels). Does FieldTrip
>>>>> handle marking artifacts restricted to some channel/trial combination?
>>>>>
>>>>> Thanks,
>>>>> Tim
>>>>>
>>>>> _______________________________________________
>>>>> fieldtrip mailing list
>>>>> fieldtrip at donders.ru.nl
>>>>> https://mailman.science.ru.nl/mailman/listinfo/fieldtrip
>>>>>
>>>>
>>>>
>>>>
>>>> --
>>>> Teresa E. Madsen, PhD
>>>> Research Technical Specialist:  *in vivo *electrophysiology & data
>>>> analysis
>>>> Division of Behavioral Neuroscience and Psychiatric Disorders
>>>> Yerkes National Primate Research Center
>>>> Emory University
>>>> Rainnie Lab, NSB 5233
>>>> 954 Gatewood Rd. NE
>>>> Atlanta, GA 30329
>>>> (770) 296-9119
>>>> braingirl at gmail.com
>>>> https://www.linkedin.com/in/temadsen
>>>>
>>>> _______________________________________________
>>>> fieldtrip mailing list
>>>> fieldtrip at donders.ru.nl
>>>> https://mailman.science.ru.nl/mailman/listinfo/fieldtrip
>>>>
>>>
>>>
>>> _______________________________________________
>>> fieldtrip mailing list
>>> fieldtrip at donders.ru.nl
>>> https://mailman.science.ru.nl/mailman/listinfo/fieldtrip
>>>
>>
>>
>>
>> --
>> Teresa E. Madsen, PhD
>> Research Technical Specialist:  *in vivo *electrophysiology & data
>> analysis
>> Division of Behavioral Neuroscience and Psychiatric Disorders
>> Yerkes National Primate Research Center
>> Emory University
>> Rainnie Lab, NSB 5233
>> 954 Gatewood Rd. NE
>> Atlanta, GA 30329
>> (770) 296-9119
>> braingirl at gmail.com
>> https://www.linkedin.com/in/temadsen
>>
>> _______________________________________________
>> fieldtrip mailing list
>> fieldtrip at donders.ru.nl
>> https://mailman.science.ru.nl/mailman/listinfo/fieldtrip
>>
>
>
> _______________________________________________
> fieldtrip mailing list
> fieldtrip at donders.ru.nl
> https://mailman.science.ru.nl/mailman/listinfo/fieldtrip
>



-- 
Teresa E. Madsen, PhD
Research Technical Specialist:  *in vivo *electrophysiology & data analysis
Division of Behavioral Neuroscience and Psychiatric Disorders
Yerkes National Primate Research Center
Emory University
Rainnie Lab, NSB 5233
954 Gatewood Rd. NE
Atlanta, GA 30329
(770) 296-9119
braingirl at gmail.com
https://www.linkedin.com/in/temadsen
-------------- next part --------------
An HTML attachment was scrubbed...
URL: <http://mailman.science.ru.nl/pipermail/fieldtrip/attachments/20170308/22203c07/attachment-0001.html>
-------------- next part --------------
function [data] = AutoArtReject_TEM(cfg,data)
% AutoArtReject_TEM performs automated artifact rejection, processing each
%   channel independently, removing clipping & large zvalue artifacts based
%   on automated thresholds (best to run on all data from a given subject
%   in one data structure, so the same threshold is applied consistently
%   across conditions), and returning the data structure re-merged across
%   channels, with NaNs in place of artifacts.
%
% Input cfg structure should contain:
%   interactsubj   = true or false, whether to select visual artifacts per
%     subject (i.e., per call to this function) & review all channels after
%     automated detection
%   interactch     = true or false, whether to preview detected artifacts
%     & select visual artifacts per channel
%   artfctdef.clip = struct as defined in ft_artifact_clip, but absdiff is
%     applied before the data is passed to that function, so it's actually
%     comparing these thresholds to the 2nd derivative of the original data
%   artfctdef.zvalue        = struct as defined in ft_artifact_zvalue
%   artfctdef.minaccepttim  = scalar as defined in ft_rejectartifact
%
% To facilitate data-handling and distributed computing you can use
%   cfg.inputfile   =  ...
%   cfg.outputfile  =  ...
% If you specify one of these (or both) the input data will be read from a
% *.mat file on disk and/or the output data will be written to a *.mat
% file. These *.mat files should contain only a single variable 'data',
% corresponding with ft_datatype_raw.
%
% written 3/2/17 by Teresa E. Madsen

%% load data, if needed

if nargin < 2 && isfield(cfg,'inputfile')
  load(cfg.inputfile);
end

%% preview data & mark unusual cross-channel artifacts

if cfg.interactsubj
  cfgtmp  = [];
  
  cfgtmp = ft_databrowser(cfgtmp,data);
  visual = cfgtmp.artfctdef.visual.artifact;  % not a field of artifact
  % because this will be reused for all channels, while the rest of
  % artifact is cleared when starting each new channel
else
  visual = [];
end

%% perform artifact detection on each channel separately

excludech = false(size(data.label));
datch = cell(size(data.label));

for ch = 1:numel(data.label)
  artifact = [];
  
  %% divide data into channels
  
  cfgtmp           = [];
  cfgtmp.channel   = data.label{ch};
  
  datch{ch} = ft_selectdata(cfgtmp,data);
  
  %% identify large zvalue artifacts
  
  cfgtmp                    = [];
  cfgtmp.artfctdef.zvalue   = cfg.artfctdef.zvalue;
  if ~isfield(cfgtmp.artfctdef.zvalue,'interactive')
    if interactch
      cfgtmp.artfctdef.zvalue.interactive = 'yes';
    else
      cfgtmp.artfctdef.zvalue.interactive = 'no';
    end
  end
  
  [~, artifact.zvalue] = ft_artifact_zvalue(cfgtmp,datch{ch});
  
  %% take 1st derivative of signal
  
  cfgtmp          = [];
  cfgtmp.absdiff  = 'yes';
  
  datd1 = ft_preprocessing(cfgtmp,datch{ch});
  
  %% define clipping artifacts
  % applies absdiff again, so it's actually working on 2nd derivative data
  
  cfgtmp                  = [];
  cfgtmp.artfctdef.clip   = cfg.artfctdef.clip;
  
  [~, artifact.clip] = ft_artifact_clip(cfgtmp,datd1);
  
  %% review artifacts if needed
  
  cfgtmp                               = [];
  cfgtmp.artfctdef.clip.artifact       = artifact.clip;
  cfgtmp.artfctdef.zvalue.artifact     = artifact.zvalue;
  cfgtmp.artfctdef.visual.artifact     = visual;
  
  if cfg.interactch
    % any new visual artifacts will be automatically added to cfgtmp
    cfgtmp = ft_databrowser(cfgtmp,datch{ch});
    keyboard  % dbcont when satisfied
    % excludech(ch) = true;   % exclude this channel if desired
  end
  
  clearvars d1dat
  
  %% replace artifactual data with NaNs
  
  cfgtmp.artfctdef.reject           = 'nan';
  cfgtmp.artfctdef.minaccepttim     = cfg.artfctdef.minaccepttim;
  
  datch{ch} = ft_rejectartifact(cfgtmp,datch{ch});
  
  % if any trials were rejected completely, exclude this channel, or it
  % won't merge properly
  if numel(datch{ch}.trial) ~= numel(data.trial)
    excludech(ch) = true;
  end
end   % for ch = 1:numel(data.label)

%% remerge each channel file into one cleaned data file

cfgtmp  = [];

if isfield(cfg,'outputfile')
  cfgtmp.outputfile   = cfg.outputfile;
end

data = ft_appenddata(cfgtmp,datch(~excludech));

%% visualize result

if interactsubj
  cfgtmp  = [];
  
  cfgtmp = ft_databrowser(cfgtmp,data); %#ok<NASGU> just for debugging
  keyboard % dbcont when satisfied
end
end


More information about the fieldtrip mailing list