TEAP (Toolbox for Emotion Analysis using Physiological Signals) doc
detect_spikes_bursts.m
Go to the documentation of this file.
1 function [HDR, s] = detect_spikes_bursts(fn, chan, varargin)
2 % DETECT_SPIKES_BURSTS detects spikes and bursts of spikes in
3 % neural recordings.
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.
7 %
8 %
9 % HDR = detect_spikes_bursts(filename, chan)
10 % ... = detect_spikes_bursts(HDR, data)
11 % ... = detect_spikes_bursts(... ,'OptimumJK')
12 % ... = detect_spikes_bursts(... ,'-o',outputFilename)
13 % ... = detect_spikes_bursts(... ,'-e',eventFilename)
14 % ... = detect_spikes_bursts(... ,'-b',burstFilename)
15 % ... = detect_spikes_bursts(... ,'slopeThreshold',slopeThreshold)
16 % ... = detect_spikes_bursts(... ,'winlen',winlen)
17 % ... = detect_spikes_bursts(... ,'dT_Burst',dT_Burst)
18 % ... = detect_spikes_bursts(... ,'dT_Exclude',dT_Exclude)
19 % ... = detect_spikes_bursts(... ,[slopeThreshold [, winlen [, dT_Burst [,dT_Exclude ]]] ])
20 % ... = detect_spikes_bursts(... ,'OptimumJK')
21 % [HDR, data] = detect_spikes_bursts(...)
22 %
23 % Input:
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.
41 % outputFilename
42 % name of file for storing the resulting data with the
43 % detected spikes and bursts in GDF format.
44 % eventFilename
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
48 % burstFilename
49 % filename for the "burst table", containing basic properties of each burst,
50 % (it is an ASCII file in <tab>-delimited format)
51 %
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.
54 %
55 %
56 % Output:
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
64 %
65 % References:
66 %
67 
68 % $Id: detect_spikes_bursts.m 2818 2011-11-07 12:23:25Z schloegl $
69 % Copyright (C) 2011 by Alois Schloegl <alois.schloegl@gmail.com>
70 % This is part of the BIOSIG-toolbox http://biosig.sf.net/
71 %
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.
76 %
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.
81 %
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/>.
84 
85 if nargin < 1,
87 end;
88 
89 if nargin < 2,
90  chan = 0;
91 end;
92 
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]
98 outFile = [];
99 evtFile = [];
100 burstFile = [];
101 
102 
103 %%%%% analyze input arguments %%%%%
104 k = 1;
105 K = 0;
106 while k <= length(varargin)
107  if ischar(varargin{k})
108  if (strcmp(varargin{k},'-o'))
109  k = k + 1;
110  outFile = varargin{k};
111  elseif (strcmp(varargin{k},'-e'))
112  k = k + 1;
113  evtFile = varargin{k};
114  elseif (strcmp(varargin{k},'-b'))
115  k = k + 1;
116  burstFile = varargin{k};
117  elseif (strcmpi(varargin{k},'slopeThreshold'))
118  k = k + 1;
119  slopeThreshold = varargin{k};
120  elseif (strcmpi(varargin{k},'winlen'))
121  k = k + 1;
122  dT = varargin{k};
123  elseif (strcmpi(varargin{k},'dT_Burst'))
124  k = k + 1;
125  dT_Burst = varargin{k};
126  elseif (strcmpi(varargin{k},'dT_Exclude'))
127  k = k + 1;
128  dT_Exclude = varargin{k};
129  elseif (strcmpi(varargin{k},'OptimumJK'))
130  slopeThreshold = 2.5; %
131  dT = .450e-3;
132  dT_Exclude = 2e-3;
133  else
134  warning(sprintf('unknown input argument <%s>- ignored',varargin{k}));
135  end;
136  elseif isnumeric(varargin{k})
137  K = K + 1
138  switch (K)
139  case {1}
140  slopeThreshold = varargin{k}; %% [V/s]
141  case {2}
142  dT = varargin{k}; %%% smoothing window length [s]
143  case {3}
144  dT_Burst = varargin{k}; %%% smaller ISI is a burst [s]
145  case {4}
146  dT_Exclude = varargin{k}; %%% smaller ISI is a burst [s]
147  otherwise
148  warning(sprintf('unknown input argument <%f> - ignored',varargin{k}));
149  end;
150  end;
151  k = k+1;
152 end;
153 
154 
155 Fs = 20000; % assumed samplerate
156 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
157 % load data
158 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
159  if ischar(fn) && exist(fn,'file') && any(size(chan)==1)
160  winlen = Fs*.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);
165  end;
166  if chan==0, chan = 1:HDR.NS; end;
167  elseif isstruct(fn)
168  HDR = fn;
169  s = chan;
170  HDR.NS = size(s,2);
171  chan = 1:HDR.NS;
172  else
173  help(mfilename);
174  end
175 
176 
177  EVENT = HDR.EVENT;
178  if ~isfield(EVENT,'DUR');
179  EVENT.DUR = zeros(size(EVENT.POS));
180  end;
181  if ~isfield(EVENT,'CHN');
182  EVENT.CHN = zeros(size(EVENT.TYP));
183  end;
184 
185 
186 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
187 % Set Parameters for Spike Detection
188 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
189  B = [.5; zeros(HDR.SampleRate*dT - 1, 1); -.5];
190  HDR.BurstTable = [];
191 
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 %%%%%%%%%%%%%%%%%%%
196  [unit, scale] = physicalunits(HDR.PhysDimCode(ch));
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)]);
202  end;
203 
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) ];
208  end;
209  end;
210 
211 
212 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
213 % Set Parameters for Burst Detection
214 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
215  HDR.EVENT = EVENT;
216  vararg = {'dT_Burst', dT_Burst, 'dT_Exclude', dT_Exclude};
217  if ~isempty(evtFile)
218  vararg{end+1}='-e';
219  vararg{end+1}=evtFile;
220  end;
221  if ~isempty(burstFile)
222  vararg{end+1}='-b';
223  vararg{end+1}=burstFile;
224  end;
225  HDR = spikes2bursts(HDR,vararg{:});
226 
227 
228 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
229 % Output
230 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
231 
232  if ~isempty(outFile)
233  %%% write data to output
234  HDR.TYPE = 'GDF';
235  HDR.VERSION = 3;
236  %[p,f,e]=fileparts(fn);
237  HDR.FILE = [];
238  HDR.FileName = outFile;
239  HDR.FILE.Path = '';
240  HDR.PhysMax = max(s);
241  HDR.PhysMin = min(s);
242  HDR.DigMax(:) = 2^15-1;
243  HDR.DigMin(:) = -2^15;
244  HDR.GDFTYP(:) = 3;
245  HDR.FLAG.UCAL = 0;
246  HDR.NRec = size(s,1);
247  HDR.SPR = 1;
248  HDR.NS = size(s,2);
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);
254  else
255  HDR = swrite(HDR,s);
256  HDR = sclose(HDR);
257  end;
258  end;
259 
sclose
function sclose(in HDR)
physicalunits
function physicalunits(in arg1)
sload
function sload(in FILENAME, in varargin)
swrite
function swrite(in HDR, in data)
detect_spikes_bursts
function detect_spikes_bursts(in fn, in chan, in varargin)
sopen
function sopen(in arg1, in PERMISSION, in CHAN, in MODE, in arg5, in arg6)
spikes2bursts
function spikes2bursts(in fn, in varargin)