


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)

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