[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