TEAP (Toolbox for Emotion Analysis using Physiological Signals) doc
detectmuscle.m
Go to the documentation of this file.
1 function [INI,S,E] = detectmuscle(S, iter, Mode)
2 % Muscle detection with inverse filtering
3 % Artifacts are indicated with NaN's.
4 %
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)
8 %
9 % iter number of iterations [default:1]
10 % INI.MU mean of S
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]
18 %
19 %
20 % References:
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.
24 %
25 
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/
29 %
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.
34 %
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.
39 %
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/>.
42 
43 
44 if nargin<2,
45  iter=[];
46 end;
47 INI = [];
48 
49 % Muscle Detection
50 if Mode==1,
51  %% TODO: Validation
52  if isempty(iter) iter = 1; end;
53  %% inverse filter
54  TH = 5;
55  [se,INI.MU] = sem(S);
56  if TH*se<abs(INI.MU),
57  [S] = center(S);
58  end;
59 
60  while iter>0,
61  [AR,RC,PE]= lattice(S',10);
62  INI.InvFilter = ar2poly(AR);
63  E = zeros(size(S));
64  for k = 1:size(S,2),
65  E(:,k) = filter(INI.InvFilter(k,:),1,S(:,k));
66  end;
67  INI.V = std(E);
68  INI.TH = INI.V * TH;
69  iter = iter-1;
70 
71  for k = 1:size(S,2),
72  S(E(:,k)>(INI.TH(k)) | E(:,k)<(-INI.TH(k)),k) = NaN;
73  end;
74  end;
75  % the following part demonstrates a possible correction
76  for k = 1:size(S,2),
77  E(:,k) = filter(INI.InvFilter(k,:),1,S(:,k));
78  end;
79  % isnan(E), returns Artifactselse
80 
81 elseif Mode==2,
82  Fs = iter;
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))];
88 
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
94 
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))];
97 
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
103  ix3 = abs(S)>75;
104 
105  ix = filtfilt(ones(Fs/10,1),1,double(ix1|ix2|ix3))>0;
106  S(ix) = NaN;
107  INI=[];
108 
109 
110 elseif Mode==3,
111  Fs = iter;
112  E = filter(ones(1,Fs),Fs,abs(diff(S))/Fs*250);
113  Threshold = 11; % uV/sample [1]
114  S(E>Threshold) = NaN;
115 
116 
117 elseif Mode==4,
118  %% TODO: free butter
119  Fs = iter;
120  if Fs>= 250,
121  [A,B]=butter(5,25*2/Fs,'high');
122  else
123  [A,B]=butter(5,[25,125]*2/Fs);
124  end;
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;
128 
129 end;
130 
131 if nargout<3, return; end;
132 
detectmuscle
function detectmuscle(in S, in iter, in Mode)