Home > chronux_code > exploratory > functionestimation > hybrid > sta.m

sta

PURPOSE ^

Spike Triggered Average %

SYNOPSIS ^

function[s,t,E] = sta(data_spk,data_lfp,smp,plt,w,T,D,err)

DESCRIPTION ^

 Spike Triggered Average                            % 
     Usage: [s,t,E] = sta(data_spk,data_lfp,smp,plt,w,T,D,err)                                               %
 INPUT                                              %
                                                    %
 Note that all times have to be consistent. If data_spk
 is in seconds, so must be sig and t. If data_spk is in 
 samples, so must sig and t. The default is seconds.

 data_spk    - strucuture array of spike times data 
               or NaN padded matrix
 data_lfp    - array of lfp data(samples x trials)         %
                                                    %
 Optional...                                        %
 plt 'n'|'r' etc                                    %
 width kernel smoothing in s                        %
 T = [-0.1 0.1] - extract this range about each spk %
 D = plot spike triggered average out to [D1 D2]    %
 err = calcluate error bars (bootstrap)             %
                                                    %
 OUTPUT                                             %
                                                    %
 s  spike triggered average                         %
 t  times                                           %
 E  bootstrap standard err                          %
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

CROSS-REFERENCE INFORMATION ^

This function calls: This function is called by:

SOURCE CODE ^

0001 function[s,t,E] = sta(data_spk,data_lfp,smp,plt,w,T,D,err)
0002 % Spike Triggered Average                            %
0003 %     Usage: [s,t,E] = sta(data_spk,data_lfp,smp,plt,w,T,D,err)                                               %
0004 % INPUT                                              %
0005 %                                                    %
0006 % Note that all times have to be consistent. If data_spk
0007 % is in seconds, so must be sig and t. If data_spk is in
0008 % samples, so must sig and t. The default is seconds.
0009 %
0010 % data_spk    - strucuture array of spike times data
0011 %               or NaN padded matrix
0012 % data_lfp    - array of lfp data(samples x trials)         %
0013 %                                                    %
0014 % Optional...                                        %
0015 % plt 'n'|'r' etc                                    %
0016 % width kernel smoothing in s                        %
0017 % T = [-0.1 0.1] - extract this range about each spk %
0018 % D = plot spike triggered average out to [D1 D2]    %
0019 % err = calcluate error bars (bootstrap)             %
0020 %                                                    %
0021 % OUTPUT                                             %
0022 %                                                    %
0023 % s  spike triggered average                         %
0024 % t  times                                           %
0025 % E  bootstrap standard err                          %
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); % create a zero padded data matrix from input structural array
0032    sz=size(data_spk); 
0033    if sz(1)>sz(2); data_spk=data_spk'; end;% transpose data to get in form compatible with Murray's routine
0034 else
0035    if sz(1)>sz(2); data_spk=data_spk'; end;% transpose data to get in form compatible with Murray's routine
0036 end
0037 sz=size(data_lfp);
0038 if sz(1)>sz(2); data_lfp=data_lfp'; end;% transpose data to get in form compatible with Murray's routine
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 % count up the spikes...
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 % local smoother...
0132 
0133 N = fix(w/smp);
0134 if N > 5
0135   mlfp = locsmooth(mlfp,N,fix(N/2)); 
0136 end
0137 
0138 % bootstrap errorbars...
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 %  errorbar(0.1*xax(2)+0.9*xax(1),0.1*yax(2)+0.9*yax(1), ...
0171 %         mean(stdlfp),'k')
0172 %plot(0.1*xax(2)+0.9*xax(1),0.1*yax(2)+0.9*yax(1),'k.')
0173 end
0174      
0175 title(['spike triggered average : ' ...
0176       num2str(Nspk) ' used : '     ...
0177       ' Errorbars are two std err']);
0178 %line(get(gca,'xlim'),mean(mlfp)*[1 1],'color','k')
0179 line([0 0],get(gca,'ylim'),'color','k')
0180 hold off
0181 
0182 
0183 
0184 
0185 
0186 
0187 
0188

Generated on Tue 07-Jun-2005 12:20:32 by m2html © 2003