2 % DETECT_SPIKES_BURSTS detects spikes and bursts of spikes in
4 % Spikes are detected when voltages increase is larger than slopeThreshold
5 % (
default 20 V/s) within a window of length winlen (
default 0.0002 s). An interspike interval
6 % large than dT_Burst (
default 75 ms) define the start of the next burst.
24 % filename: name of source file
25 % chan list of channels that should be analyzed (
default is 0: all channels)
26 % HDR header structure obtained by SOPEN, SLOAD, or meXSLOAD
27 % data signal data that should be analyzed
28 % slopeThreshold [
default: 5V/s] Spike is detected when
29 % slope (over time winlen) exceeds
this value
30 % winlen [
default: .2e-3 s] windowlength in seconds
for computing slope
31 % dT_Burst [
default: 50e-3 s] am inter-spike-interval (ISI) exceeding
this value,
32 % marks the beginning of a
new burst
33 % dT_Exclude an interspike interval smaller than
this value, indicates a
34 %
double detection, and the second detection is deleted.
35 % in
case of several consecutive ISI
's smaller than this value,
36 % all except the first spikes are deleted.
37 % 'OptimumJK
' overrides previously defined settings, and uses the optimal setting based on
38 % an JK's data set. It sets slopeThreshold = 2.5, dT=450e-6, dT_Exclude = 2e-3;
39 % These settings were found based on 18 different recordings containing
40 % 23999 spikes, 3967 bursts, and over 9700 seconds of recorded data.
42 % name of file
for storing the resulting data with the
43 % detected spikes and bursts in GDF format.
45 % filename to store
event inforamation in GDF format.
this is similar to
46 % the outputFile, except that the signal data is not included and is, therefore,
47 % much smaller than the outputFile
49 % filename
for the
"burst table", containing basic properties of each burst,
50 % (it is an ASCII file in <tab>-delimited format)
52 % Arguments can appear in any order and multiple times (except
for filename, chan, HDR and data),
53 % In
case of conflicting definitions, the latest definition has highest precedence and is used.
57 % HDR header structure as defined in biosig
58 % HDR.EVENT includes the detected spikes and bursts.
59 % HDR.BurstTable contains
for each burst (each in a row) the following 5 numbers:
60 % channel number, sweep number, OnsetTime within sweep [s],
61 % number of spikes within burst, and average inter-spike interval (ISI) [ms]
62 % data signal data, one channel per column
63 % between segments, NaN values
for 0.1s are introduced
69 % Copyright (C) 2011 by Alois Schloegl <alois.schloegl@gmail.com>
70 % This is part of the BIOSIG-toolbox http:
72 % BioSig is free software: you can redistribute it and/or modify
73 % it under the terms of the GNU General Public License as published by
74 % the Free Software Foundation, either version 3 of the License, or
75 % (at your option) any later version.
77 % BioSig is distributed in the hope that it will be useful,
78 % but WITHOUT ANY WARRANTY; without even the implied warranty of
79 % MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
80 % GNU General Public License
for more details.
82 % You should have received a copy of the GNU General Public License
83 % along with BioSig. If not, see <http://www.gnu.org/licenses/>.
93 %%%%%
default settings %%%%%
94 dT = .2e-3; %%% smoothing window length [s]
95 dT_Burst = 50e-3; %%% smaller ISI is a burst [s]
96 dT_Exclude = []; %%%
for identifying
double detections, spikes with smaller ISI are excluded
97 slopeThreshold = 5; %% [V/s]
103 %%%%% analyze input arguments %%%%%
106 while k <= length(varargin)
107 if ischar(varargin{k})
108 if (strcmp(varargin{k},
'-o'))
110 outFile = varargin{k};
111 elseif (strcmp(varargin{k},
'-e'))
113 evtFile = varargin{k};
114 elseif (strcmp(varargin{k},
'-b'))
116 burstFile = varargin{k};
117 elseif (strcmpi(varargin{k},
'slopeThreshold'))
119 slopeThreshold = varargin{k};
120 elseif (strcmpi(varargin{k},
'winlen'))
123 elseif (strcmpi(varargin{k},
'dT_Burst'))
125 dT_Burst = varargin{k};
126 elseif (strcmpi(varargin{k},
'dT_Exclude'))
128 dT_Exclude = varargin{k};
129 elseif (strcmpi(varargin{k},
'OptimumJK'))
130 slopeThreshold = 2.5; %
134 warning(sprintf('unknown input argument <%s>- ignored',varargin{k}));
136 elseif isnumeric(varargin{k})
140 slopeThreshold = varargin{k}; %% [V/s]
142 dT = varargin{k}; %%% smoothing window length [s]
144 dT_Burst = varargin{k}; %%% smaller ISI is a burst [s]
146 dT_Exclude = varargin{k}; %%% smaller ISI is a burst [s]
148 warning(sprintf(
'unknown input argument <%f> - ignored',varargin{k}));
155 Fs = 20000; % assumed samplerate
156 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
158 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
159 if ischar(fn) && exist(fn,
'file') && any(size(chan)==1)
161 [s, HDR] =
sload(fn, 0, 'NUMBER_OF_NAN_IN_BREAK', winlen);
162 if Fs < HDR.SampleRate,
163 winlen = HDR.SampleRate * .1;
164 [s, HDR] =
sload(fn, 0, 'NUMBER_OF_NAN_IN_BREAK', winlen);
166 if chan==0, chan = 1:HDR.NS; end;
178 if ~isfield(EVENT,'DUR');
179 EVENT.DUR = zeros(size(EVENT.POS));
181 if ~isfield(EVENT,'CHN');
182 EVENT.CHN = zeros(size(EVENT.TYP));
186 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
187 % Set Parameters for Spike Detection
188 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
189 B = [.5; zeros(HDR.SampleRate*dT - 1, 1); -.5];
192 for ch = chan(:)'; % look for each channel
193 % only voltage channels are considered
194 if (bitand(HDR.PhysDimCode(ch), hex2dec('ffe0')) == 4256), %%
physicalunits('V'),
195 %%%%%%% Spike Detection %%%%%%%%%%%%%%%%%%%
197 tmp = scale * filter(B, dT, s(:,ch));
198 OnsetSpike = find( diff (tmp > slopeThreshold) > 0); %% spike onset time [samples]
199 % --- remove
double detections < 1 ms
200 if ~isempty(dT_Exclude)
201 OnsetSpike = OnsetSpike([1; 1+find(diff(OnsetSpike) > Fs * dT_Exclude)]);
204 EVENT.TYP = [EVENT.TYP; repmat(hex2dec('0201'), size(OnsetSpike))];
205 EVENT.POS = [EVENT.POS; OnsetSpike];
206 EVENT.DUR = [EVENT.DUR; repmat(1, size(OnsetSpike))];
207 EVENT.CHN = [EVENT.CHN; repmat(ch, size(OnsetSpike,1), 1) ];
212 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
213 % Set Parameters for Burst Detection
214 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
216 vararg = {
'dT_Burst', dT_Burst,
'dT_Exclude', dT_Exclude};
219 vararg{end+1}=evtFile;
221 if ~isempty(burstFile)
223 vararg{end+1}=burstFile;
228 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
230 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
233 %%% write data to output
236 %[p,f,e]=fileparts(fn);
238 HDR.FileName = outFile;
240 HDR.PhysMax = max(s);
241 HDR.PhysMin = min(s);
242 HDR.DigMax(:) = 2^15-1;
243 HDR.DigMin(:) = -2^15;
246 HDR.NRec = size(s,1);
249 HDR.Dur = 1/HDR.SampleRate;
250 HDR = rmfield(HDR,
'AS');
251 HDR =
sopen(HDR,
'w');
252 if (HDR.FILE.FID < 0)
253 fprintf(2,
'Warning can not open file <%s> - GDF file can not be written\n',HDR.FileName);