2 % Muscle detection with inverse filtering
3 % Artifacts are indicated with NaN
's.
5 % [INI,S,E] = detectmuscle(S [, iter [,1]])
6 % [INI,S,E] = detectmuscle(S, Fs, Mode) with Mode>1
7 % [INI,S,E] = detectmuscle(S, arg2, Mode)
9 % iter number of iterations [default:1]
11 % INI.InvFilter coefficients of inverse filter
12 % S outlier replaced by NaN
13 % E isnan(E) indicates muscle artifact
14 % Mode 1: [default] inverse filtering
15 % 2: based on gradient, range and amplitude [BrainVision method]
16 % 3: slope > 11 uV/sample [1]
17 % 4: beta2 > 0.9 uV^2/Hz [1]
21 % [1] Van de Velde, M., Van Erp, G., Cluitmans, P., 1998.
22 % Detection of muscle artefact in the normal human awake EEG.
23 % Electroencephalography and Clinical Neurophysiology 107 (2), 149-158.
26 % $Id: detectmuscle.m 2202 2009-10-27 12:06:45Z schloegl $
27 % Copyright (C) 2003,2008,2009 by Alois Schloegl <a.schloegl@ieee.org>
28 % This is part of the BIOSIG-toolbox http://biosig.sf.net/
30 % BioSig is free software: you can redistribute it and/or modify
31 % it under the terms of the GNU General Public License as published by
32 % the Free Software Foundation, either version 3 of the License, or
33 % (at your option) any later version.
35 % BioSig is distributed in the hope that it will be useful,
36 % but WITHOUT ANY WARRANTY; without even the implied warranty of
37 % MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
38 % GNU General Public License for more details.
40 % You should have received a copy of the GNU General Public License
41 % along with BioSig. If not, see <http://www.gnu.org/licenses/>.
52 if isempty(iter) iter = 1; end;
61 [AR,RC,PE]= lattice(S',10);
62 INI.InvFilter = ar2poly(AR);
65 E(:,k) = filter(INI.InvFilter(k,:),1,S(:,k));
72 S(E(:,k)>(INI.TH(k)) | E(:,k)<(-INI.TH(k)),k) = NaN;
75 % the following part demonstrates a possible correction
77 E(:,k) = filter(INI.InvFilter(k,:),1,S(:,k));
79 % isnan(E), returns Artifactselse
83 % Criterion
for bad gradient:
84 % Maximal allowed voltage step / sampling point: 100.00 µV
85 % Mark as bad before event: 100.00 ms
86 % Mark as bad after event: 100.00 ms
87 ix1 = [abs(diff(S,[],1))>100;zeros(1,size(S,2))];
89 % Criterion
for bad max-min:
90 % Maximal allowed absolute difference: 150.00 µV
91 % Interval Length: 100.00 ms
92 % Mark as bad before event: 100.00 ms
93 % Mark as bad after event: 100.00 ms
95 ix2 = abs(filter([-1;zeros(Fs/10-1,1);1],1,S))>150;
96 ix2 = [ix2(Fs/20+1:end,:);zeros(Fs/20+1,size(S,2))];
98 % Criterion
for bad amplitude:
99 % Minimal allowed amplitude: -75.00 µV
100 % Maximal allowed amplitude: 75.00 µV
101 % Mark as bad before event: 100.00 ms
102 % Mark as bad after event: 100.00 ms
105 ix = filtfilt(ones(Fs/10,1),1,
double(ix1|ix2|ix3))>0;
112 E = filter(ones(1,Fs),Fs,abs(diff(S))/Fs*250);
113 Threshold = 11; % uV/sample [1]
114 S(E>Threshold) = NaN;
121 [A,B]=butter(5,25*2/Fs,
'high');
123 [A,B]=butter(5,[25,125]*2/Fs);
125 E = filtfilt(ones(1,Fs),Fs,filtfilt(A,B,S).^2);
126 Threshold = 0.9; % uV^2/Hz ??? [1]
127 S(E>Threshold) = NaN;
131 if nargout<3,
return; end;