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

Generated on Tue 24-Aug-2004 15:55:33 by m2html © 2003