/
spm_uomeeg_findbadchans.m
116 lines (103 loc) · 3.96 KB
/
spm_uomeeg_findbadchans.m
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
function [D,bads,corrfname] = spm_uomeeg_findbadchans(S)
% Identify and mark as bad any channels with a sub-threshold correlation
% with its N nearest neighbours.
% FORMAT: [D,bads] = spm_uomeeg_findbadchans(S)
% INPUT: Struct 'S' with fields:
% S.D - MEEG object or filename of MEEG object
% S.thresh - NN correlation z-score threshold (default: -3)
% S.nNN - Number of nearest neighbours to consider (def: 3)
%
% OUTPUT:
% D - data with bad chans marked bad (also saved)
% bads - indices of bad channels
% corrfname - filename output (correlations, means, zscores)
%
% spm_uomeeg tools ** HIGHLY EXPERIMENTAL **
% by Jason Taylor (09/Mar/2018) jason.taylor@manchester.ac.uk
% *Base it on z-score instead of r thresh?
% *Base it on spectral power rather than time-domain signal?
%--------------------------------------------------
% This function looks for correlations between each channel and its 3
% nearest neighbours, which should generally be high because of the spatial
% smoothness of EEG. It goes through iterations of removing a single bad
% channel, then choosing the 3 NNs for each channel omitting those
% previously labelled bad (because it's not fair to expect a high
% correlation with a bad channel).
% Check inputs:
try thresh = S.thresh; catch, thresh = -3; end
try nNN = S.nNN; catch, nNN = 3; end
% Load SPM-format data file:
D = spm_eeg_load(S.D);
[~,fstem] = fileparts(D.fname);
fprintf('\n\n');
fprintf('++ %s\n',datestr(now));
fprintf('++ RUNNING spm_uomeeg_findbadchans ON %s\n',D.fname);
fprintf('++ USING: thresh = %g\tnNN = %d\n',thresh,nNN);
%% Find bad channels
% Note: includes those already labelled as 'bad'
% (but they will still be marked bad in the end)
goods = indchantype(D,'EEG');
% Initialise:
zcnnm = -Inf; i=0; bads=[];
allcnnm={}; allR={}; allz={};
d = selectdata(D,chanlabels(D,goods),[],[]);
d = reshape(d,size(d,1),size(d,2)*size(d,3));
[spec,fq] = spectopo(d,0,D.fsample,'plot','on');
while min(zcnnm)<thresh
i=i+1; % iterate
fprintf('++ Iteration %d ... ',i);
% Get correlations:
%d = selectdata(D,chanlabels(D,goods),[],[]);
%d = d(goods,:);
isgood = ismember(goods,indchantype(D,'EEG'));
%d = d(isgood,:);
spec = spec(isgood,:);
%R = corrcoef(reshape(d,size(d,1),size(d,2)*size(d,3))');
R = corrcoef(spec(:,fq<48)');
allR{i} = R;
% Get NNs:
pos = D.coor2D(goods)';
if ~isempty(which('knnsearch')); % requires toolbox
[idx,dist] = knnsearch(pos,pos,'dist','euclidean','k',nNN+1);
else
[idx,dist] = spm_uomeeg_getnearestchans(pos,pos,nNN+1);
dist = abs(dist); % ??
end
% Get NN correlations for each channel:
cnn = [];
for c=1:length(goods)
jind=idx(c,[2:nNN+1]); % +1 skips self!
cnn(c,:) = R(c,jind);
end
cnnm = mean(cnn,2); % mean over NNs
allcnnm{i} = cnnm;
zcnnm = (cnnm-mean(cnnm))/std(cnnm);
allz{i} = zcnnm;
% Plot topography, histogram:
% spm_eeg_plotScalpData(cnnm,pos,D.chanlabels(goods));
% set(gca,'clim',[0 .6]);
% colormap gray
% title(sprintf('Iteration %d',i));
% figure; hist(cnnm,[.05:.05:.95]); title(sprintf('Iteration %d',i));
% drawnow
if min(zcnnm)<thresh
ind = find(zcnnm==min(zcnnm));
bads(i) = goods(ind);
goods = setdiff(goods,bads(i));
fprintf('MinZ: %.2f (r=%.2f)\tBad: %s\n',zcnnm(ind),cnnm(ind),char(D.chanlabels(bads(i))));
else
%fprintf('No corr<%.2f: DONE.\n',thresh)
fprintf('MinZ: %.2f (r=%.2f)\tNo Zcorr<%.2f: DONE.\n',min(zcnnm),min(cnnm),thresh)
end
end
fprintf('++ Found %d bad channels.\n',length(bads));
% % Add (to existing) any new bad channels:
% if any(bads)
% D = badchannels(D,bads,1);
% D.save;
% end
%
% % Save correlation matrices and mean nn correlation vectors (per iteration)
% savefname = sprintf('badchan_correlations_%s',fstem);
% save(savefname,'allR','allcnnm','allz');
return