spike_train.m

function spikes = spike_train(drive, sr, recfun, nfibers)
% spikes=spike_train(drive,sr,recfun,N) - simulate auditory-nerve spike generation
%
%   spikes: vector of spike times (s)
%  
%   drive: (spikes/s) time-varying rate
%   sr: sampling rate (Hz)
%   recfun: refractory function (or dead time)
%   N: number of fibers to simulate
%
% Spike times are double precision floats (no sampling involved). Spike 
% generation is efficient.
%
% Spikes are generated in three steps:
% 1 - 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 (not included).
% 2 - The spike train is "thinned" probabilistically.  The probability 
% of survival of a spike is equal to the instantaneous rate (at spike time)
% normalized by its max.
% 3 - The spike train is further thinned probabilistically to simulate 
% refractory effects.


if nargin==0; test_code; return; end

if nargin<4||isempty(nfibers); nfibers=1; end
if nargin < 3||isempty(recfun); recfun = 0; end % no refractory effects
if nargin < 2; error('!'); return; end

drive=drive(:);

if min(drive) < 0
	error ('instantaneous rate must be non-negative');
end

if nfibers>1 % run repeatedly, merge spikes
    spikes=[];
    for iFiber=1:nfibers
        spikes=[spikes;spike_train(drive,sr,recfun,1)];
    end
    spikes=sort(spikes);
    return;
end

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% First step: generate enough spikes to cover the duration of the driving function
max_rate = max(drive);		% peak instantaneous rate
drive = drive/max_rate;		% normalize
nspikes=numel(drive)/sr * max_rate * 1.1;% request a bit more

isis = spike_poisson(nspikes, max_rate); 
spikes = cumsum(isis);
spikes(spikes*sr > numel(drive)-1) = []; % trim to duration of input

% Second step: thin based on instantaneous rate
draw = rand(numel(spikes),1);
idxSpikes = fix(spikes .* sr) + 1;		% spike times expressed in samples
keep = draw < drive(idxSpikes);
spikes = keep .* spikes;
spikes=nonzeros(spikes);

if numel(spikes) < 1; warning('no spikes produced (rate too small?)'); end

% Third step: thin again based on interspike interval
if isa(recfun, 'function_handle') || recfun~=0
	spikes = spike_refractory(spikes, recfun);
end
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%


if nargout==0
    disp('spike_train: no output requested, plot AC histogram');
    spike_ach(spikes); % expects spike times
    disp(['rate (spikes/s): ',num2str(spike_rate(spikes))])
    clear spikes
end

end % spike_train

% test/example code
function test_code
    disp('spike_train test code');
    disp('100 Hz HWR sine, max_rate 1000 spikes/s, 1 ms dead time');
    max_rate=1000; % spikes/s
    sr=44100; % Hz, sampling rate of driving function
    f=100; % Hz
    D=10; % s
    drive=max(0,sin(2*pi*(1:round(sr*D)')/sr*f))*max_rate;
    recfun=0.001;
    nfibers=1;
    tic;
    spikes=spike_train(drive,sr,recfun,nfibers); 
    toc;
    figure(1); clf
    subplot 221
    plot(linspace(0,D,numel(drive)),drive); xlim([0 0.02]); 
    xlabel('time (s)'); ylabel('rate (spikes/s)'); title ('instantaneous rate');
    subplot 222
    binwidth=0.0001;
    spike_psth(spikes,binwidth,1/f);
    subplot 223
    spike_isih(spikes,binwidth);
    subplot 224
    spike_ach(spikes,binwidth,0.02);
end % function