Home > chronux_0.5 > 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 for n=1:NT
0075   indx = find(~isnan(data(n,:)) & data(n,:)>=T(1) & data(n,:)<=T(2));
0076   ND(n) = length(indx);
0077   D(n,1:ND(n)) = data(n,indx);
0078   m = m + ND(n); 
0079 end
0080 N_tot = m;
0081 N_max = max(ND);
0082 D = D(:,1:N_max);
0083 
0084 % if the kernel density is such that there are on average
0085 % one or less spikes under the kernel then the units are probably wrong
0086 
0087 L = N_tot/(NT*(T(2)-T(1)));
0088 % if 2*L*NT*sig < 1 | L < 0.1
0089 %   disp('Spikes very low density: are the units right? is the kernel width sensible?')
0090 %   disp(['Total events: ' num2str(fix(100*N_tot)/100) ' sig: ' ...
0091 %         num2str(fix(1000*sig)) 'ms T: ' num2str(fix(100*T)/100) ' events/sig: ' ...
0092 %         num2str(fix(100*N_tot*sig/(T(2)-T(1)))/100)])
0093 % end
0094 
0095 %    Smear each spike out
0096 %    std dev is sqrt(rate*(integral over kernal^2)/trials)
0097 %    for Gaussian integral over Kernal^2 is 1/(2*sig*srqt(pi))
0098 
0099 if nargin < 6
0100   N_pts =  fix(5*(T(2)-T(1))/sig);
0101   t = linspace(T(1),T(2),N_pts);
0102 else
0103   N_pts = length(t);
0104 end
0105   
0106 RR = zeros(NT,N_pts);
0107 f = 1/(2*sig^2);
0108 for n=1:NT
0109   for m=1:ND(n)
0110     RR(n,:) = RR(n,:) + exp(-f*(t-D(n,m)).^2);
0111   end
0112 end
0113 RR = RR*(1/sqrt(2*pi*sig^2));
0114 if NT > 1; R = mean(RR); else R = RR;end
0115 
0116 if err == 1
0117   E = sqrt(R/(2*NT*sig*sqrt(pi)));
0118 elseif err == 2
0119   Nboot = 10;
0120   mE = 0;
0121   sE = 0;
0122   for b=1:Nboot
0123     indx = floor(NT*rand(1,NT)) + 1;
0124     mtmp = mean(RR(indx,:));
0125     mE = mE + mtmp;
0126     sE = sE + mtmp.^2;
0127   end
0128   E = sqrt((sE/Nboot - mE.^2/Nboot^2));
0129 end
0130 
0131 % if adaptive warp sig so that on average the number of spikes
0132 % under the kernel is the same but regions where there is
0133 % more data have a smaller kernel
0134 
0135 if adapt 
0136   sigt = mean(R)*sig./R;
0137   RR = zeros(NT,N_pts);
0138   f = 1./(2*sigt.^2);
0139   for n=1:NT
0140     for m=1:ND(n)
0141       RR(n,:) = RR(n,:) + exp(-f.*(t-D(n,m)).^2);
0142     end
0143     RR(n,:) = RR(n,:).*(1./sqrt(2*pi*sigt.^2));
0144   end
0145   if NT > 1; R = mean(RR); else R = RR;end
0146 
0147   if err == 1
0148     E = sqrt(R./(2*NT*sigt*sqrt(pi)));
0149   elseif err == 2
0150     Nboot = 10;
0151     mE = 0;
0152     sE = 0;
0153     for b=1:Nboot
0154       indx = floor(NT*rand(1,NT)) + 1;
0155       mtmp = mean(RR(indx,:));
0156       mE = mE + mtmp;
0157       sE = sE + mtmp.^2;
0158     end
0159     E = sqrt((sE/Nboot - mE.^2/Nboot^2)); 
0160   end
0161 end  
0162 
0163 if plt == 'n';return;end
0164 plot(t,R,plt)
0165 hold on
0166 if err > 0
0167   plot(t,R+2*E,'g')
0168   plot(t,R-2*E,'g')
0169 end
0170 %axis([T(1)+(4*sig) T(2)-(4*sig) 0 1.25*max(R)])
0171 axis([T(1)+(4*sig) T(2)-(4*sig) 0 max(R+2*E)+10])
0172 xlabel('time (s)')
0173 ylabel('rate (Hz)')
0174 title(['Trial averaged rate : Gaussian Kernel :'  ...
0175         ' sigma = ' num2str(1000*sig) 'ms'])
0176 hold off
0177 
0178 
0179 
0180 
0181 
0182 
0183 
0184 
0185 
0186 
0187

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