0001 function[N,B,E] = isi(data,T,err,Nbins,plt)
0002
0003
0004
0005
0006
0007
0008
0009
0010
0011
0012
0013
0014
0015
0016
0017
0018
0019
0020
0021
0022 if nargin < 1; error('I need data!'); end
0023 data=padNaN(data);
0024 data=data';
0025 if nargin < 2; T = [min(data(:,1)) max(max(data))]; end
0026 if nargin < 3; err = 0;end
0027 if nargin < 4; Nbins = -1; end
0028 if nargin < 5; plt = 'r'; end
0029
0030 if isempty(T); T = [min(min(data)) max(max(data))]; end
0031 if isempty(err); err = 0;end
0032 if isempty(Nbins); Nbins = -1; end
0033 if isempty(plt); plt = 'r'; end
0034
0035
0036
0037
0038 NT = length(data(1,:));
0039 for n=1:NT
0040 indx = find(data(:,n) >= T(1) & data(:,n) <= T(2) ...
0041 & ~isnan(data(:,n)));
0042 if isempty(indx)
0043 NI(n) = 0;
0044 else
0045 NI(n) = length(indx)-1;
0046 index(n).keep=indx;
0047 end
0048 end
0049
0050
0051
0052
0053 I = zeros(NT,max(NI));
0054 IT = [];
0055 for n=1:NT
0056 I(n,1:NI(n)) = diff(data(index(n).keep,n));
0057 IT = [IT I(n,1:NI(n))];
0058 end
0059
0060 Mx = max(IT);
0061 if Nbins == -1
0062 Nbins = floor(sum(NI)/30);
0063 Med = median(IT);
0064 Nbins = max(floor(Nbins*Mx/Med),10);
0065 end
0066
0067 B = linspace(0,Mx,Nbins);
0068
0069 N = zeros(NT,Nbins);
0070 for n=1:NT
0071 N(n,:) = hist(I(n,1:NI(n)),B);
0072 end
0073
0074
0075
0076 if NT > 1;Ns = sum(N)/NT;else Ns = N;end
0077 if ~strcmp(plt,'n')
0078 bar(B,NT*Ns);
0079 end
0080
0081
0082
0083 if NT > 4 && err == 1
0084 MN = 0;
0085 SN = 0;
0086 for n=1:NT
0087 JK = (NT*Ns - N(n,:))/(NT-1);
0088 MN = MN + JK;
0089 SN = SN + JK.^2;
0090 end
0091 MN = MN/NT;
0092 SN = SN/NT;
0093 E = sqrt((NT-1)*(SN - MN.^2));
0094 if ~strcmp(plt,'n')
0095 hold on
0096 errorbar(B,NT*Ns,NT*2*E,'r-')
0097 hold off
0098 end
0099 end
0100 N = NT*Ns;
0101
0102
0103
0104
0105
0106
0107
0108