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)
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:
• padNaN Creates a padded data matrix from input structural array of spike times
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
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
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