# psth

## PURPOSE

function to plot trial averaged rate smoothed by

## SYNOPSIS

function [R,t,E] = psth(data,sig,plt,T,err,t)

## DESCRIPTION

      function to plot trial averaged rate smoothed by
      a Gaussian kernel - visual check on stationarity
      Usage: [R,t,E] = psth(data,sig,plt,T,err,t)

      Inputs:
     Note that all times have to be consistent. If data
     is in seconds, so must be sig and t. If data is in
     samples, so must sig and t. The default is seconds.
     data            structural array of spike times
     sig             std dev of Gaussian (default 50ms)
                     (minus indicates adaptive with
                      approx width equal to mod sig)
     plt = 'n'|'r' etc      (default 'r')
     T is the time interval (default all)
     err - 0 = none
           1 = Poisson
           2 = Bootstrap over trials (default)
     (both are based on 2* std err rather than 95%)
     t   = times to evaluate psth at

     The adaptive estimate works by first estimating
     the psth with a fixed kernel width (-sig) it
     then alters the kernel width so that the number
     of spikes falling under the kernel is the same on
     average but is time dependent.  Reagions rich
     in data therefore have their kernel width reduced

     Outputs:

     R = rate
     t = times
     E = errors (standard error)

## CROSS-REFERENCE INFORMATION

This function calls:

- padNaN Creates a padded data matrix from input structural array of spike times

This function is called by:

## SOURCE CODE

    0001 function [R,t,E] = psth(data,sig,plt,T,err,t)
    0002 %  function to plot trial averaged rate smoothed by
    0003 %  a Gaussian kernel - visual check on stationarity
    0004 %  Usage: [R,t,E] = psth(data,sig,plt,T,err,t)
    0005 %
    0006 %  Inputs:
    0007 % Note that all times have to be consistent. If data
    0008 % is in seconds, so must be sig and t. If data is in
    0009 % samples, so must sig and t. The default is seconds.
    0010 % data            structural array of spike times
    0011 % sig             std dev of Gaussian (default 50ms)
    0012 %                 (minus indicates adaptive with
    0013 %                  approx width equal to mod sig)
    0014 % plt = 'n'|'r' etc      (default 'r')
    0015 % T is the time interval (default all)
    0016 % err - 0 = none
    0017 %       1 = Poisson
    0018 %       2 = Bootstrap over trials (default)
    0019 % (both are based on 2* std err rather than 95%)
    0020 % t   = times to evaluate psth at
    0021 %
    0022 % The adaptive estimate works by first estimating
    0023 % the psth with a fixed kernel width (-sig) it
    0024 % then alters the kernel width so that the number
    0025 % of spikes falling under the kernel is the same on
    0026 % average but is time dependent.  Reagions rich
    0027 % in data therefore have their kernel width reduced
    0028 %
    0029 % Outputs:
    0030 %
    0031 % R = rate
    0032 % t = times
    0033 % E = errors (standard error)
    0034
    0035 if nargin <1 ; error('I need data!');end
    0036 [data]=padNaN(data); % create a zero padded data matrix from input structural array
    0037 sz=size(data);
    0038 % transposes data so that the longer dimension is assumed to be time
    0039 % Procedure used to bring input data into form compatible with that required
    0040 % for Murray's routine
    0041 if sz(1)>sz(2); data=data';end;
    0042 if nargin < 2; sig = 0.05;end
    0043 if nargin < 3; plt = 'r';end
    0044 if nargin < 5; err = 2; end
    0045
    0046 if isempty(sig); sig = 0.05;end
    0047 if isempty(plt); plt = 'r';end
    0048 if isempty(err); err = 2; end
    0049
    0050 adapt = 0;
    0051 if sig < 0; adapt = 1; sig = -sig; end
    0052
    0053 %  to avoid end effects increase time interval by the width
    0054 %  of the kernel (otherwise rate is too low near edges)
    0055
    0056 if nargin < 4;
    0057   T(1) = min(data(:,1));
    0058   T(2) = max(max(data));
    0059 else
    0060   T(1) = T(1)-4*sig;
    0061   T(2) = T(2)+4*sig;
    0062 end
    0063
    0064 % calculate NT and ND and filter out times of interest
    0065
    0066 NT = length(data(:,1));
    0067 if NT < 4 && err == 2
    0068   disp('Using Poisson errorbars as number of trials is too small for bootstrap')
    0069   err = 1;
    0070 end
    0071
    0072 m = 1;
    0073 D = zeros(size(data));
    0074 ND=zeros(1,NT);
    0075 for n=1:NT
    0076   indx = find(~isnan(data(n,:)) & data(n,:)>=T(1) & data(n,:)<=T(2));
    0077   ND(n) = length(indx);
    0078   D(n,1:ND(n)) = data(n,indx);
    0079   m = m + ND(n);
    0080 end
    0081 N_tot = m;
    0082 N_max = max(ND);
    0083 D = D(:,1:N_max);
    0084
    0085 % if the kernel density is such that there are on average
    0086 % one or less spikes under the kernel then the units are probably wrong
    0087
    0088 L = N_tot/(NT*(T(2)-T(1)));
    0089 if 2*L*NT*sig < 1 || L < 0.1
    0090   disp('Spikes very low density: are the units right? is the kernel width sensible?')
    0091   disp(['Total events: ' num2str(fix(100*N_tot)/100) ' sig: ' ...
    0092         num2str(fix(1000*sig)) 'ms T: ' num2str(fix(100*T)/100) ' events/sig: ' ...
    0093         num2str(fix(100*N_tot*sig/(T(2)-T(1)))/100)])
    0094 end
    0095
    0096 %    Smear each spike out
    0097 %    std dev is sqrt(rate*(integral over kernal^2)/trials)
    0098 %    for Gaussian integral over Kernal^2 is 1/(2*sig*srqt(pi))
    0099
    0100 if nargin < 6
    0101   N_pts =  fix(5*(T(2)-T(1))/sig);
    0102   t = linspace(T(1),T(2),N_pts);
    0103 else
    0104   N_pts = length(t);
    0105 end
    0106
    0107 RR = zeros(NT,N_pts);
    0108 f = 1/(2*sig^2);
    0109 for n=1:NT
    0110   for m=1:ND(n)
    0111     RR(n,:) = RR(n,:) + exp(-f*(t-D(n,m)).^2);
    0112   end
    0113 end
    0114 RR = RR*(1/sqrt(2*pi*sig^2));
    0115 if NT > 1; R = mean(RR); else R = RR;end
    0116
    0117 if err == 1
    0118   E = sqrt(R/(2*NT*sig*sqrt(pi)));
    0119 elseif err == 2
    0120   Nboot = 10;
    0121   mE = 0;
    0122   sE = 0;
    0123   for b=1:Nboot
    0124     indx = floor(NT*rand(1,NT)) + 1;
    0125     mtmp = mean(RR(indx,:));
    0126     mE = mE + mtmp;
    0127     sE = sE + mtmp.^2;
    0128   end
    0129   E = sqrt((sE/Nboot - mE.^2/Nboot^2));
    0130 end
    0131
    0132 % if adaptive warp sig so that on average the number of spikes
    0133 % under the kernel is the same but regions where there is
    0134 % more data have a smaller kernel
    0135
    0136 if adapt
    0137   sigt = mean(R)*sig./R;
    0138   RR = zeros(NT,N_pts);
    0139   f = 1./(2*sigt.^2);
    0140   for n=1:NT
    0141     for m=1:ND(n)
    0142       RR(n,:) = RR(n,:) + exp(-f.*(t-D(n,m)).^2);
    0143     end
    0144     RR(n,:) = RR(n,:).*(1./sqrt(2*pi*sigt.^2));
    0145   end
    0146   if NT > 1; R = mean(RR); else R = RR;end
    0147
    0148   if err == 1
    0149     E = sqrt(R./(2*NT*sigt*sqrt(pi)));
    0150   elseif err == 2
    0151     Nboot = 10;
    0152     mE = 0;
    0153     sE = 0;
    0154     for b=1:Nboot
    0155       indx = floor(NT*rand(1,NT)) + 1;
    0156       mtmp = mean(RR(indx,:));
    0157       mE = mE + mtmp;
    0158       sE = sE + mtmp.^2;
    0159     end
    0160     E = sqrt((sE/Nboot - mE.^2/Nboot^2));
    0161   end
    0162 end
    0163
    0164 if plt == 'n';return;end
    0165 plot(t,R,plt)
    0166 hold on
    0167 if err > 0
    0168   plot(t,R+2*E,'g')
    0169   plot(t,R-2*E,'g')
    0170 end
    0171 %axis([T(1)+(4*sig) T(2)-(4*sig) 0 1.25*max(R)])
    0172 axis([T(1)+(4*sig) T(2)-(4*sig) 0 max(R+2*E)+10])
    0173 xlabel('time (s)')
    0174 ylabel('rate (Hz)')
    0175 title(['Trial averaged rate : Gaussian Kernel :'  ...
    0176         ' sigma = ' num2str(1000*sig) 'ms'])
    0177 hold off
    0178
    0179
    0180
    0181
    0182
    0183
    0184
    0185
    0186
    0187
    0188