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

Generated on Tue 16-Aug-2005 21:33:45 by m2html © 2003