Description of psth

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
psth

PURPOSE

SYNOPSIS

DESCRIPTION

CROSS-REFERENCE INFORMATION

SOURCE CODE
