


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) %
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

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';return;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