[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