0001 function[s,t,E] = sta(data_spk,data_lfp,smp,plt,w,T,D,err)
0002
0003
0004
0005
0006
0007
0008
0009
0010
0011
0012
0013
0014
0015
0016
0017
0018
0019
0020
0021
0022
0023
0024
0025
0026
0027
0028
0029 if nargin < 3;error('Require spike, lfp and lfp times');end
0030 if isstruct(data_spk)
0031 [data_spk]=padNaN(data_spk);
0032 sz=size(data_spk);
0033 if sz(1)>sz(2); data_spk=data_spk'; end;
0034 else
0035 if sz(1)>sz(2); data_spk=data_spk'; end;
0036 end
0037 sz=size(data_lfp);
0038 if sz(1)>sz(2); data_lfp=data_lfp'; end;
0039 verbose = 1;
0040 t = smp;
0041 if nargin < 4; plt = 'r'; end
0042 if nargin < 5; w = 0; end
0043 if nargin < 6; T = [min(t) max(t)]; end
0044 if nargin < 7; D = 0.25*[-1 1]; end
0045 if nargin < 8; err = 1;end
0046
0047 if isempty(plt); plt = 'r'; end
0048 if isempty(w); w = 0; end
0049 if isempty(T); T = [min(t) max(t)]; end
0050 if isempty(D); D = 0.25*[-1 1]; end
0051 if isempty(err); err = 1;end
0052
0053 if w > (T(2)-T(1))/2
0054 disp('Smoothing > data segment : should be in seconds : turn off smoothing')
0055 w = 0;
0056 end
0057
0058 sz = size(data_spk);
0059 NT = sz(1);
0060 mlfp = 0;
0061 slfp = 0;
0062 Nspk = 0;
0063 smp = t(2)-t(1);
0064 if D(1) <= 0 & D(2) >= 0
0065 t1 = [D(1):smp:(-smp+eps) 0:smp:D(2)+eps];
0066 else
0067 t1 = (round(D(1)/smp)*smp):smp:D(2);
0068 end
0069
0070
0071
0072 if err
0073 for n=1:NT
0074 indx = find(t>T(1)&t<T(2));
0075 lfp = data_lfp(n,indx);
0076 spk = data_spk(n,find(data_spk(n,:)>T(1) & ...
0077 data_spk(n,:)<T(2) & ...
0078 data_spk(n,:)~=0));
0079 tt = t(indx);
0080 if ~isempty(spk) > 0
0081 ND = length(spk);
0082 for s=1:ND
0083 spktime = spk(s);
0084 t0 = tt-spktime + eps;
0085 if min(t0)>(D(1)-smp) | max(t0)<(D(2)+smp); break;end
0086 Nspk = Nspk + 1;
0087 end
0088 end
0089 end
0090 Err = zeros(Nspk,length(t1));
0091 Nspk = 0;
0092 end
0093
0094 for n=1:NT
0095 indx = find(t>T(1)&t<T(2));
0096 lfp = data_lfp(n,indx);
0097 spk = data_spk(n,find(data_spk(n,:)>T(1) & ...
0098 data_spk(n,:)<T(2) & ...
0099 data_spk(n,:)~=0));
0100 tt = t(indx);
0101 if ~isempty(spk) > 0
0102 ND = length(spk);
0103 for s=1:ND
0104 spktime = spk(s);
0105 t0 = tt-spktime + eps;
0106 if min(t0) < (D(1)-smp) & max(t0) > (D(2)+smp);
0107 indx = find(t0<D(1));
0108 indx = indx(length(indx));
0109 offset = (t0(indx)-D(1))/smp;
0110 indx = indx:(indx+length(t1)-1);
0111 lfp_t1 = lfp(indx) + (lfp(indx+1)-lfp(indx))*offset;
0112 Nspk = Nspk + 1;
0113 mlfp = mlfp + lfp_t1;
0114 slfp = slfp + lfp_t1.^2;
0115 Err(Nspk,:) = lfp_t1;
0116 end
0117 end
0118 end
0119 end
0120 if Nspk == 0
0121 if verbose;disp('No spikes in interval');end
0122 t = t1;
0123 s = zeros(length(t),1);
0124 E = zeros(length(t),1);
0125 return
0126 end
0127 mlfp = mlfp/Nspk;
0128 slfp = slfp/Nspk;
0129 stdlfp = sqrt((slfp - mlfp.^2)/Nspk);
0130
0131
0132
0133 N = fix(w/smp);
0134 if N > 5
0135 mlfp = locsmooth(mlfp,N,fix(N/2));
0136 end
0137
0138
0139
0140 if err == 1;
0141 Nboot = 20;
0142 blfp = 0;
0143 slfp = 0;
0144 for n = 1:Nboot
0145 indx = floor(Nspk*rand(1,Nspk)) + 1;
0146 lfptmp = mean(Err(indx,:));
0147 if N > 5
0148 lfptmp = locsmooth(lfptmp,N,fix(N/2));
0149 end
0150 blfp = blfp + lfptmp;
0151 slfp = slfp + lfptmp.^2;
0152 end
0153 stdlfp = sqrt((slfp/Nboot - blfp.^2/Nboot^2));
0154 end
0155
0156 s = mlfp-mean(mlfp);
0157 E = stdlfp;
0158 t = t1;
0159
0160 cols = 'krbgycm';
0161 if plt == 'n';return;end
0162 plot(t1,s,plt)
0163 xax = get(gca,'xlim');
0164 yax = get(gca,'ylim');
0165 if err == 1
0166 me = real(2*mean(stdlfp));
0167 line(xax,me*[1 1],'color','b')
0168 line(xax,-me*[1 1],'color','b')
0169 line(xax,0*[1 1],'color','k')
0170
0171
0172
0173 end
0174
0175 title(['spike triggered average : ' ...
0176 num2str(Nspk) ' used : ' ...
0177 ' Errorbars are two std err']);
0178
0179 line([0 0],get(gca,'ylim'),'color','k')
0180 hold off
0181
0182
0183
0184
0185
0186
0187
0188