[FieldTrip] marking artifacts by channel + trial

Teresa Madsen braingirl at gmail.com
Fri Mar 3 23:31:04 CET 2017


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
-------------- next part --------------
An HTML attachment was scrubbed...
URL: <http://mailman.science.ru.nl/pipermail/fieldtrip/attachments/20170303/ecd1e7f3/attachment-0002.html>
-------------- next part --------------
function [artifact,data,times] = artifact_nan2zero_TEM(data)
% ARTIFACT_NAN2ZERO_TEM marks NaNs that occur uniformly across all channels
%   as artifacts, taking FT format data structure & returning the same
%   format as ft_artifact_xxx, for input to ft_rejectartifact.  Non-uniform
%   NaNs (those not present on all channels) are replaced with 0s to avoid
%   breaking analysis functions.  Also returns times of replaced NaNs by
%   trial & channel, so they can be changed back to NaNs in freq output.
%
% written 3/2/17 by Teresa E. Madsen

artifact = [];
times = cell(numel(data.trial), numel(data.label));

for tr = 1:numel(data.trial)
  % find NaNs to mark as artifacts (present uniformly across all channels)
  trlnan = isnan(data.trial{tr});   % identify NaNs by channel & timepoint
  allnan = all(trlnan,1);   % need to specify dim in case of single channel
  
  % find, save timepoints, & replace non-uniform NaNs (not on all chs) w/ 0
  replacenan = trlnan & repmat(~allnan,size(trlnan,1),1);
  for ch = 1:numel(data.label)
    times{tr,ch} = data.time{tr}(replacenan(ch,:));  % ID before replacing
  end
  data.trial{tr}(replacenan) = 0;   % replace these w/ 0s
  
  if any(allnan)
    % determine the file sample #s for this trial
    trsamp = data.sampleinfo(tr,1):data.sampleinfo(tr,2);
    
    while any(allnan)
      % start from the end so sample #s don't shift
      endnan = find(allnan,1,'last');
      allnan = allnan(1:endnan);  % remove any non-NaNs after this
      
      % find last non-NaN before the NaNs
      beforenan = find(~allnan,1,'last');
      
      if isempty(beforenan)   % if no more non-NaNs
        begnan = 1;
        allnan = false;   % while loop ends
      else  % still more to remove - while loop continues
        begnan = beforenan + 1;
        allnan = allnan(1:beforenan);  % remove the identified NaNs
      end
      
      % identify file sample #s that correspond to beginning and end of
      % this chunk of NaNs and append to artifact
      artifact = [artifact; trsamp(begnan) trsamp(endnan)]; %#ok<AGROW>
    end   % while any(tnan)
  end   % if any(tnan)
end   % for tr = 1:numel(data.trial)

end


More information about the fieldtrip mailing list