-
Notifications
You must be signed in to change notification settings - Fork 0
/
p_burst.m
235 lines (231 loc) · 7.31 KB
/
p_burst.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
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
233
234
235
%Function p_burst
%NOTE: Calls POISSCDF function in the STATISTICS toolbox
%CALLING:
%[BOB, EOB, SOB]=p_burst(InTrain, StartT, StopT)
% InTrain: Timestamps of the spike train.
% StartT and StopT : Time limits for analysis. Spike count between
% these time limits are used to calculate Average Spike Rate MU
% for analyzing the complete spike train.
% Note: Burst results may not be valid outside limits.
%OUTPUTS:
% BOB: A row vector containing indices for Begining Of Burst.
% EOB: A row vector containing indices for End Of Burst.
% SOB: A row vector containing Surprise Of Burst.
% Note: Spike Train(BOB)gives the actual times. The Spike Train
% however must be stripped off "zero","Nan" and "Inf"
%PROBLEMS:
% If you get an error: One or more output arguments
% not assigned during call to 'distchck'
% Change the spike Train dimension by passing
% a transpose(Spike Train)
% First Version Jan 5, 1998: S.Chenchal Rao
function [BOB, EOB, SOB]=p_burst(InTrain, StartT, StopT)
%Returns Spike Index of BOB, EOB abd Surprise Of Burst
%The spike train is a series of timestamps
%******************User Parameters************************
MaxXT=30;%Max Xtra Time
MaxXS=10;%Max Xtra Spikes
MinSPInBurst=3;%5%2;%Minimum spkes in a Burst
Anchor=50;%Anchor Time
Signif=0.025;%0.025%targ_pburst ;%0.05;0.01
UserSI=-log(Signif);
Tol=1e-300;
%****************Spike Train Properties*****************
if(size(InTrain,1))>1
InTrain=InTrain';
end
SPT=InTrain(find(InTrain > 0 & InTrain < Inf))';
StopT=max(0,StopT);%traps NaN
StartT=max(0,StartT);%traps NaN
if(~StopT),StopT=max(SPT);,end
if(~StartT),StartT=min(SPT);,end
minT=min(StartT,StopT);
maxT=max(StartT,StopT);
StartT=minT;StopT=maxT;
Duration=StopT-StartT;
%Average Spike Rate MU
if(~isempty(SPT))
%Calculate Average spike rate between start and stop time
MU=(length(find(SPT>= StartT & SPT <= StopT)))/(Duration);
%if there are no spikes in the selected interval
%find MU for the entire trial
if (MU==0)
MU=(length(SPT)-1)/(max(SPT)-min(SPT));
end
%if there are no spikes in the spike train or
%n spikes in train is less than 4
if(MU==0 | length(SPT)<=4)
BOB=[];
EOB=[];
SOB=[];
return
end
else
BOB=[];
EOB=[];
SOB=[];
return
end
MaxSpikes=length(SPT);
ISI=diff(SPT);%Inter Spike Intervals
%################Parameter Initializations##############
%******Flags, Counters, Indices and Rel. OPs************
ISBURST=0;
BNo=1;CurrBNo=0;%Burst No for indexing
MinBOB=5000;%Dummy
MaxEOB=0;%Dummy
OldBOB=0;OldEOB=0;
PrevBOB=0;PrevEOB=0;
CurrBOB=0;CurrEOB=0;
CurrEOB_XS=0;CurrXS=0;CurrXT=0;%Vars for Extra Time
MaxSI=0;
%*******************Iteration vars**********************
Iterate=1;IC=1;%Iteration Counter
FromI=0;ToI=0; %Current portion of Spike Train
FspAB=1;%First spike After Burst
Temp=0;Done=0;
%******************Output Arguments*********************
BOB=[];EOB=[];SOB=[];%MUST be set to EMPTY
%########################################################
while(FspAB <= MaxSpikes-1 | ~Done)
Iterate=1;
while(Iterate)
%***************FIND EOB****************************
if (IC==1)
FromI=FspAB;
else
FromI=CurrBOB;
end
cISI=cumsum(ISI(FromI:length(ISI)))+Anchor;
Prob=(poisscdf([1:length(cISI)]', cISI.*MU))+Tol;
Prob=(1-Prob);SI=-log(Prob);
%find index that maximizes SI
Temp=(find(diff(SI)<0));
if (length(Temp)>1)
CurrEOB=Temp(min(find(Temp >1)));
if(CurrEOB==1 & IC==1)
FspAB=Temp(min(find(diff(Temp) >MinSPInBurst)))+1;
if (isempty(FspAB) | FspAB>=MaxSpikes-MinSPInBurst)
Done=1;
ISBURST=0;
Iterate=0;
end
break %Break out of Iteration
end
CurrEOB_XS=0;Curr_XT=0;
%check for extra spikes and extra time
%Find index to next MaxSI
CurrEOB_XS=min(find(SI > SI(CurrEOB)));
if(~isempty(CurrEOB_XS))
CurrXS=CurrEOB_XS-CurrEOB-1;
CurrXT=SPT(CurrEOB_XS)-SPT(CurrEOB);
if(CurrXS <= MaxXS & CurrXT <= MaxXT)
CurrEOB=CurrEOB_XS;
CurrEOB_XS=0;
end
end
elseif (length(Temp)==1 & ~isempty(Temp))
CurrEOB=Temp;
elseif(isempty(Temp))
%Index to the max
CurrEOB=length(SI);
else
end
CurrEOB=CurrEOB+FromI;
%***************FIND BOB****************************
ToI=FspAB;
BSPT=SPT(CurrEOB:-1:ToI);
cISI=cumsum(abs(diff(BSPT)))+Anchor;
Prob=poisscdf([1:length(cISI)]',cISI.*MU)+Tol;
Prob=(1-Prob);
SI=-log(Prob);
Temp=max(find(SI==max(SI)));
%What if Temp = [];??
if (isempty(Temp))
disp('ERROR finding BOB ==> Temp is EMPTY')
[Temp]
end
CurrBOB=find(SPT==BSPT(Temp+1));
%if diff between CurrEOB and recently found BOB is < criteria
if (CurrEOB-CurrBOB) < MinSPInBurst
ISBURST=0;
Iterate=0;
FspAB=CurrEOB+1;
FromI=FspAB;
break;%Break out of Iteration loop
end
%Check for convergence among 3 values for each
if(OldBOB==CurrBOB & OldEOB==CurrEOB )
%Converged
%disp(['Burst Number ',num2str(BNo),' Converged in ',num2str(IC),' iterations'])
%disp(['BOB Index ',num2str(CurrBOB),' EOB Index 'num2str(CurrEOB)])
Iterate=0;
ISBURST=1;
break
else
IC=IC+1;
PrevBOB=OldBOB;
PrevEOB=OldEOB;
OldBOB=CurrBOB;
OldEOB=CurrEOB;
end
%find max and min
CurrBNo=max(CurrBNo,BNo);
if(CurrBNo)
MaxEOB=max(CurrEOB,MaxEOB);
MinBOB=min(CurrBOB,MinBOB);
end
if(IC==10)
MaxEOB=CurrEOB;MinBOB=CurrBOB;
elseif(IC>20)
CurrEOB=MaxEOB;CurrBOB=MinBOB;
Iterate=0;
ISBURST=1;
else
end
%disp(['End of Iteration cycle : ',num2str(IC-1)])
end%while(Iterate)
if(ISBURST)
%CHECK CRITERIA
MaxSI=-log(1-(poisscdf(CurrEOB-CurrBOB,(SPT(CurrEOB)-SPT(CurrBOB))*MU)));
if(CurrEOB-CurrBOB >= MinSPInBurst & MaxSI>UserSI)
BOB(BNo)=CurrBOB;
EOB(BNo)=CurrEOB;
SOB(BNo)=MaxSI;
BNo=BNo+1;
end
FspAB=CurrEOB+1;
CurrBOB=FspAB;
IC=1;
else
IC=1;
end
%Check limits
if(FspAB >= MaxSpikes-2)
Done=1;
Iterate=0;
break
end
end%while(FspAB <= MaxSpikes-1)
%****************************DISPLAY***********************************
% %clf
% BOB
% EOB
% SOB
% if ~isempty(BOB) | ~isempty(EOB) | ~isempty(SOB)
% figure
% colors_='rbm';
% colors_=repmat(colors_,1,ceil(length(EOB)/length(colors_)));
% %It takes less time to put '+' than draw line/ticks
% %h=plot(SPT,1,'k+');
% set(line([SPT SPT],[0.95 1.05]),'color',[0 0 0],'linewidth',0.1);
% hold on
% axis([0 max(SPT) 0 2])
% set(findobj('type','axes'),'color',[1 1 1],'ytick',[],'box','on')
% for EB=1:length(EOB)
% plot(SPT(BOB(EB):EOB(EB)),1.1,strcat(colors_(EB),'+'));
% end
% %This part depends on the user...
% % MIN INTER BURST INTERVAL
% %InterBDiff=SPT(BOB(2:length(BOB)))-SPT(EOB(1:length(EOB)-1));
% end