Home > chronux > spectral_analysis > pointtimes > psth.m

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

Generated on Fri 28-Sep-2012 12:34:30 by m2html © 2005