function [t,s] = fastgammatrain(duration,meanrate,order)
% function [t,s] = fastgammatrain(duration,meanrate,order)
%
% returns spike train output of an integrate and fire neuron with
% a random spiking threshold which is distributed according to
% a gamma distribution of order n. n=1 gives a poisson train.
% increasing n gives increasingly regular spike trains.
%
% where
% duration = duration of spiketrain (s)
% meanrate = mean firing rate of the spike train (spikes/s)
% order = order of the gamma distribution for resetting threshold
%
% return parameters are:
% t = vector containing the time index for spike train (s)
% s = spike train (1’s represent spikes, 0’s represent no spikes)
delta_t = 0.001; % resolution of spike train = 1 ms
t = 0:delta_t:duration;
idur=length(t);
Vrandth = gamrnd(order,1/order,1,idur);
itonextspike =round(Vrandth/(meanrate*delta_t));
itonextspike = itonextspike(itonextspike > 2); % refractory period = 2 ms
ispikes = cumsum(itonextspike);
ispikes = ispikes(ispikes < (duration/delta_t));
nspikes = length(ispikes);
s = sparse(1,ispikes,1,1,idur,nspikes);
if (nargout == 0)
plot(t,s)
end
s = full(s);