function y = spike_train(x, sr)
% y=spike_train(x,sr) - simulate spike generation process
% as inhomogenous poisson process

if min(x) < 0
	error ('driving function must be non-negative');
end

% Spikes are generated by an inhomogenous Poisson process with 
% refractory effects.
%
% Spikes are generated in three steps. 
% First, a spike train is generated by a _homogenous_ Poisson process
% with a rate equal to the peak instantaneous rate of the driving function.
% Spike times are derived from interspike intervals that follow an 
% exponential distribution.  The first spike is supposed to occur at 0.
% In a second step, spikes are weeded out probabilistically.  The probability 
% of "survival" of a spike is equal to the normalized driving function.
% In a third step, spikes are weeded out with a probabiltiy that depends
% on the interval that precedes them, to simulate refractory effects.

% First step:
max_rate = max(x);		% peak instantaneous rate
x = x ./ max_rate;		% normalize
N = size(x, 2);
D = N / sr;				% duration of driving function

% generate enough intervals to cover duration of driving function
ns = D * max_rate;
z = poisson(ns, max_rate);

% transform to times, remove any out of range of driving function
t = [0,cumsum(z)];
mask = t * sr < N-1;
t = t .* mask;
t = nonzeros(t)';

% Second step:
% For each spike, draw a random number, and compare it to the
% instantaneous rate at the time of the spike.  Remove
% unlucky spikes.
ns = size(t, 2);
zz = rand(1, ns);
tt = fix(t .* sr) + 1;		% spike times expressed in samples
xx = x(tt);				
mask = zz < xx;
t = mask .* t;
t = nonzeros(t)';

% Third step:
% Refractory effects: remove spikes preceded by interval shorter than 
% the dead time.  This is a crude model of refractory effects.
dead = 0.001;		% ms
y = spike_dead(t, dead);