Skip to content

Commit

Permalink
implemented a test script for issue #623, use all channels in case EE…
Browse files Browse the repository at this point in the history
…G selection returns empty
  • Loading branch information
robertoostenveld committed Jan 12, 2018
1 parent 6c9c587 commit 5492bab
Show file tree
Hide file tree
Showing 2 changed files with 56 additions and 35 deletions.
75 changes: 40 additions & 35 deletions ft_qualitycheck.m
Original file line number Diff line number Diff line change
Expand Up @@ -83,18 +83,18 @@
%% ANALYSIS
if strcmp(cfg.analyze,'yes')
tic

% checks
cfg = ft_checkconfig(cfg, 'dataset2files', 'yes'); % translate into datafile+headerfile

% these will be replaced by more appropriate values
info.filename = cfg.dataset;
info.datasetname = 'unknown';
info.starttime = 'unknown';
info.startdate = 'unknown';
info.stoptime = 'unknown';
info.stopdate = 'unknown';

% the exportname is also used in the cron job
exportname = qualitycheck_exportname(cfg.dataset);
[iseeg, ismeg, isctf, fltp] = filetyper(cfg.dataset);
Expand All @@ -104,12 +104,12 @@
info = read_ctf_hist(cfg.dataset);
end
end

% add info
info.event = ft_read_event(cfg.dataset);
info.hdr = ft_read_header(cfg.dataset);
info.filetype = fltp;

% trial definition
cfgdef = [];
cfgdef.dataset = cfg.dataset;
Expand All @@ -119,7 +119,7 @@
cfgdef = ft_definetrial(cfgdef);
ntrials = size(cfgdef.trl,1)-1; % remove last trial
timeunit = cfgdef.trialdef.triallength;

% channelselection for jump detection (all) and for FFT (brain)
if ismeg
allchans = ft_channelselection({'MEG','MEGREF'}, info.hdr.label);
Expand All @@ -129,12 +129,17 @@
jumpthreshold = 1e-10;
elseif iseeg
allchans = ft_channelselection('EEG', info.hdr.label);
if isempty(allchans)
% some EEG systems and data files use non-standard channel names that are not detected automatically
ft_warning('no EEG channels detected, selecting all channels');
allchans = info.hdr.label;
end
chans = allchans; % brain
allchanindx = match_str(info.hdr.label, allchans);
chanindx = match_str(chans, allchans);
jumpthreshold = 1e4;
end

% find headcoil channels
if isctf % this fails for older CTF data sets
Nx = strmatch('HLC0011', info.hdr.label); % x nasion coil
Expand All @@ -151,28 +156,28 @@
headpos.label = {'Nx';'Ny';'Nz';'Lx';'Ly';'Lz';'Rx';'Ry';'Rz'};
headpos.avg = NaN(length(headpos.label), ntrials);
headpos.grad = info.hdr.grad;

if numel(cat(1,Nx,Ny,Nz,Lx,Ly,Lz,Rx,Ry,Rz))==9
hasheadpos = true;
else
hasheadpos = false;
end

end % if

% analysis settings
cfgredef = [];
cfgredef.length = 1;
cfgredef.overlap = 0;

cfgfreq = [];
cfgfreq.output = 'pow';
cfgfreq.channel = allchans;
cfgfreq.method = 'mtmfft';
cfgfreq.taper = 'hanning';
cfgfreq.keeptrials = 'no';
cfgfreq.foilim = [0 min(info.hdr.Fs/2, 400)];

% output variables
timelock.dimord = 'chan_time';
timelock.label = allchans;
Expand All @@ -183,35 +188,35 @@
timelock.range = NaN(length(allchans), ntrials); % updated in loop
timelock.min = NaN(length(allchans), ntrials); % updated in loop
timelock.max = NaN(length(allchans), ntrials); % updated in loop

freq.dimord = 'chan_freq_time';
freq.label = allchans;
freq.freq = (cfgfreq.foilim(1):cfgfreq.foilim(2));
freq.time = (timeunit-timeunit/2:timeunit:timeunit*ntrials-timeunit/2);
freq.powspctrm = NaN(length(allchans), length(freq.freq), ntrials); % updated in loop

summary.dimord = 'chan_time';
summary.time = (timeunit-timeunit/2:timeunit:timeunit*ntrials-timeunit/2);
summary.label = {'Mean';'Median';'Min';'Max';'Range';'HmotionN';'HmotionL';'HmotionR';'LowFreqPower';'LineFreqPower';'Jumps'};
summary.avg = NaN(length(summary.label), ntrials); % updated in loop

% try add gradiometer info
if isfield(info.hdr, 'grad')
timelock.grad = info.hdr.grad;
freq.grad = info.hdr.grad;
summary.grad = info.hdr.grad;
end


% process trial by trial
for t = 1:ntrials
fprintf('analyzing trial %s of %s \n', num2str(t), num2str(ntrials));

% preprocess
cfgpreproc = cfgdef;
cfgpreproc.trl = cfgdef.trl(t,:);
data = ft_preprocessing(cfgpreproc); clear cfgpreproc;

% determine headposition
if isctf && hasheadpos
headpos.avg(1,t) = mean(data.trial{1,1}(Nx,:) * 100); % meter to cm
Expand All @@ -224,44 +229,44 @@
headpos.avg(8,t) = mean(data.trial{1,1}(Ry,:) * 100);
headpos.avg(9,t) = mean(data.trial{1,1}(Rz,:) * 100);
end

% update values
timelock.avg(:,t) = mean(data.trial{1}(allchanindx,:),2);
timelock.median(:,t) = median(data.trial{1}(allchanindx,:),2);
timelock.range(:,t) = max(data.trial{1}(allchanindx,:),[],2) - min(data.trial{1}(allchanindx,:),[],2);
timelock.min(:,t) = min(data.trial{1}(allchanindx,:),[],2);
timelock.max(:,t) = max(data.trial{1}(allchanindx,:),[],2);

% detect jumps
for c = 1:size(data.trial{1}(allchanindx,:),1)
timelock.jumps(c,t) = length(find(diff(data.trial{1,1}(allchanindx(c),:)) > jumpthreshold));
end

% FFT and noise estimation
redef = ft_redefinetrial(cfgredef, data); clear data;
FFT = ft_freqanalysis(cfgfreq, redef); clear redef;
freq.powspctrm(:,:,t) = FFT.powspctrm;
summary.avg(9,t) = mean(mean(findpower(0, 2, FFT, chanindx))); % Low Freq Power
summary.avg(10,t) = mean(mean(findpower(cfg.linefreq-1, cfg.linefreq+1, FFT, chanindx))); clear FFT; % Line Freq Power

toc
end % end of trial loop

% determine headmotion: distance from initial trial (in cm)
if isctf && hasheadpos
summary.avg(6,:) = sqrt(sum((headpos.avg(1:3,:)-repmat(headpos.avg(1:3,1),1,size(headpos.avg,2))).^2,1)); % N
summary.avg(7,:) = sqrt(sum((headpos.avg(4:6,:)-repmat(headpos.avg(4:6,1),1,size(headpos.avg,2))).^2,1)); % L
summary.avg(8,:) = sqrt(sum((headpos.avg(7:9,:)-repmat(headpos.avg(7:9,1),1,size(headpos.avg,2))).^2,1)); % R
end

% summarize/mean and store variables of brain info only
summary.avg(1,:) = mean(timelock.avg(chanindx,:),1);
summary.avg(2,:) = mean(timelock.median(chanindx,:),1);
summary.avg(3,:) = mean(timelock.min(chanindx,:),1);
summary.avg(4,:) = mean(timelock.max(chanindx,:),1);
summary.avg(5,:) = mean(timelock.range(chanindx,:),1);
summary.avg(11,:) = mean(timelock.jumps(chanindx,:),1);

% save to .mat
if strcmp(cfg.savemat, 'yes')
if isctf && hasheadpos
Expand All @@ -271,12 +276,12 @@
save(exportname, 'info','timelock','freq','summary');
end
end

end % end of analysis

%% VISUALIZATION
if strcmp(cfg.visualize, 'yes')

% load data
if strcmp(cfg.analyze, 'no')
if ~isempty(cfg.matfile)
Expand All @@ -287,15 +292,15 @@
fprintf('loading %s \n', exportname);
load(exportname);
end

% determine number of 1-hour plots to be made
nplots = ceil(length(freq.time)/(cfg.plotunit/10));

% create GUI-like figure(s)
for p = 1:nplots
fprintf('visualizing %s of %s \n', num2str(p), num2str(nplots));
toi = [p*cfg.plotunit-(cfg.plotunit-5) p*cfg.plotunit-5]; % select 1-hour chunks

tmpcfg.latency = toi;
temp_timelock = ft_selectdata(tmpcfg, timelock);
temp_freq = ft_selectdata(tmpcfg, freq);
Expand All @@ -308,7 +313,7 @@
draw_figure(info, temp_timelock, temp_freq, temp_summary, toi);
clear temp_timelock; clear temp_freq; clear temp_summary; clear toi;
end

% export to .PNG and .PDF
if strcmp(cfg.saveplot, 'yes')
[pathstr,name,extr] = fileparts(exportname);
Expand Down Expand Up @@ -560,7 +565,7 @@ function draw_figure(varargin)
'Units','normalized',...
'color','white',...
'Position',[.05 .08 .9 .52]);

hmotions = ([summary.avg(8,:)' summary.avg(7,:)' summary.avg(6,:)'])*10;
boxplot(h.HmotionAxes, hmotions, 'orientation', 'horizontal', 'notch', 'on');
set(h.HmotionAxes,'YTick',1:3);
Expand Down Expand Up @@ -601,7 +606,7 @@ function draw_figure(varargin)
'Units','normalized',...
'color','white',...
'Position',[.08 .73 .89 .22]);

plot(h.HmotionTimecourseAxes, summary.time, clipat(summary.avg(6,:)*10, 0, 10), ...
summary.time, clipat(summary.avg(7,:)*10, 0, 10), ...
summary.time, clipat(summary.avg(8,:)*10, 0, 10), 'LineWidth',2);
Expand Down Expand Up @@ -752,7 +757,7 @@ function draw_figure(varargin)
'color','white',...
'Units','normalized',...
'Position',[0.4 0.05 0.55 0.4]);

MEGchans = ft_channelselection('MEG', timelock.label);
MEGchanindx = match_str(timelock.label, MEGchans);
cfgtopo = [];
Expand Down
16 changes: 16 additions & 0 deletions test/test_issue623.m
Original file line number Diff line number Diff line change
@@ -0,0 +1,16 @@
function test_issue623

% MEM 2gb
% WALLTIME 00:20:00

%%

filename = dccnpath('/home/common/matlab/fieldtrip/data/test/original/eeg/brainvision/Mischa.vhdr');

%%

cfg = [];
cfg.dataset = filename;
cfg.saveplot = 'no';
cfg.savemat = 'no';
[info, timelock, freq, summary] = ft_qualitycheck(cfg);

0 comments on commit 5492bab

Please sign in to comment.