# 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:

- locsmooth Running line fit (using local linear regression) - 1d only, continuous
- padNaN Creates a padded data matrix from input structural array of spike times

This function is called by:

- staogram

## 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    sz=size(data_spk);
    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,data_spk(n,:)>T(1) && data_spk(n,:)<T(2) && data_spk(n,:)~=0);
    0077     tt = t(indx);
    0078     if ~isempty(spk) > 0
    0079       ND = length(spk);
    0080       for s=1:ND
    0081         spktime = spk(s);
    0082         t0 = tt-spktime + eps;
    0083         if min(t0)>(D(1)-smp) || max(t0)<(D(2)+smp); break;end
    0084         Nspk = Nspk + 1;
    0085       end
    0086     end
    0087   end
    0088   Err = zeros(Nspk,length(t1));
    0089   Nspk = 0;
    0090 end
    0091
    0092 for n=1:NT
    0093   indx = find(t>T(1)&t<T(2));
    0094   lfp = data_lfp(n,indx);
    0095   spk = data_spk(n,data_spk(n,:)>T(1) && data_spk(n,:)<T(2) && data_spk(n,:)~=0);
    0096   tt = t(indx);
    0097   if ~isempty(spk) > 0
    0098     ND = length(spk);
    0099     for s=1:ND
    0100       spktime = spk(s);
    0101       t0 = tt-spktime + eps;
    0102       if min(t0) < (D(1)-smp) && max(t0) > (D(2)+smp);
    0103         indx = find(t0<D(1));
    0104         indx = indx(length(indx));
    0105         offset = (t0(indx)-D(1))/smp;
    0106         indx = indx:(indx+length(t1)-1);
    0107         lfp_t1 = lfp(indx) + (lfp(indx+1)-lfp(indx))*offset;
    0108         Nspk = Nspk + 1;
    0109         mlfp = mlfp + lfp_t1;
    0110         slfp = slfp + lfp_t1.^2;
    0111         Err(Nspk,:) = lfp_t1;
    0112       end
    0113     end
    0114   end
    0115 end
    0116 if Nspk == 0
    0117   if verbose;disp('No spikes in interval');end
    0118   t = t1;
    0119   s = zeros(length(t),1);
    0120   E = zeros(length(t),1);
    0121   return
    0122 end
    0123 mlfp = mlfp/Nspk;
    0124 slfp = slfp/Nspk;
    0125 stdlfp = sqrt((slfp - mlfp.^2)/Nspk);
    0126
    0127 % local smoother...
    0128
    0129 N = fix(w/smp);
    0130 if N > 5
    0131   mlfp = locsmooth(mlfp,N,fix(N/2));
    0132 end
    0133
    0134 % bootstrap errorbars...
    0135
    0136 if err == 1;
    0137   Nboot = 20;
    0138   blfp = 0;
    0139   slfp = 0;
    0140   for n = 1:Nboot
    0141     indx = floor(Nspk*rand(1,Nspk)) + 1;
    0142     lfptmp = mean(Err(indx,:));
    0143     if N > 5
    0144       lfptmp = locsmooth(lfptmp,N,fix(N/2));
    0145     end
    0146     blfp = blfp + lfptmp;
    0147     slfp = slfp + lfptmp.^2;
    0148   end
    0149   stdlfp = sqrt((slfp/Nboot - blfp.^2/Nboot^2));
    0150 end
    0151
    0152 s = mlfp-mean(mlfp);
    0153 E = stdlfp;
    0154 t = t1;
    0155
    0156 %cols = 'krbgycm';
    0157 if plt == 'n';return;end
    0158 plot(t1,s,plt)
    0159 xax = get(gca,'xlim');
    0160 %yax = get(gca,'ylim');
    0161 if err == 1
    0162   me = real(2*mean(stdlfp));
    0163   line(xax,me*[1 1],'color','b')
    0164   line(xax,-me*[1 1],'color','b')
    0165   line(xax,0*[1 1],'color','k')
    0166 %  errorbar(0.1*xax(2)+0.9*xax(1),0.1*yax(2)+0.9*yax(1), ...
    0167 %         mean(stdlfp),'k')
    0168 %plot(0.1*xax(2)+0.9*xax(1),0.1*yax(2)+0.9*yax(1),'k.')
    0169 end
    0170
    0171 title(['spike triggered average : ' ...
    0172       num2str(Nspk) ' used : '     ...
    0173       ' Errorbars are two std err']);
    0174 %line(get(gca,'xlim'),mean(mlfp)*[1 1],'color','k')
    0175 line([0 0],get(gca,'ylim'),'color','k')
    0176 hold off
    0177
    0178
    0179
    0180
    0181
    0182
    0183
    0184