% ONE1  Exercise 1 from CHAPTER 1 of Dayan and Abbott

% BN168   Spring 2008


% Generate spikes for 10 s (or longer if you want better statistics) using
% a Poisson spike generator with a constant rate of 100 Hz, and record
% their times of occurrence. Compute the coefficient of variation of the
% interspike intervals, and the Fano factor for spike counts obtained
% over counting intervals ranging from 1 to 100 ms. Plot the interspike
% interval histogram.

clear
epoch = 100;                        % length of epoch in sec
bin_length = .001;                  % bin length (discretization of time axis) in sec
t_disp = 1;                        % length of time to display (in seconds)

time_axis = 0:bin_length:epoch;     % time axis
N = length(time_axis);              % total number of bins in epoch
b_disp = t_disp/bin_length;         % length of time to display (in bin numbers)

rate = 100;                         % firing rate (Hz)

p = rate*bin_length;                % probability of observing a spike in any bin
p = min(p,1);

% the Poisson spike train is generated by independent Bernoulli random variables:
spike_train = rand(size(time_axis)) < p;      % independent Bernoulli trials,  

number_spikes = sum(spike_train);

figure(1)
subplot(2,1,1)
h = stem(time_axis(1:b_disp),spike_train(1:b_disp));
set(h(1),'Marker','none')           % removes the circles from the stem plot
set(h(1),'LineWidth',.1)
set(gca,'YLim',[0 1.5])
title('HOMOGENEOUS  POISSON  SPIKE  TRAIN')
xlabel('time (sec)')
text(.25,1.3,['total length of  epoch: ' num2str(epoch) ' sec'])
text(.25,1.1,['total number of spikes: ' num2str(number_spikes)])

[c_v,ISI] = coef_var(spike_train,bin_length);       % ISI and coefficient of variation of ISI
                                                    % std over mean. = 1 for an exponential RV

subplot(2,1,2)
max_interval = max(ISI);
hist(ISI,0:bin_length:max_interval);                % displays the ISI histogram
histogram = hist(ISI,0:bin_length:max_interval);    % ISI histogram
set(gca,'XLim',[0 max_interval])
set(gca,'YLim',[0 max(histogram)])
title('ISI  HISTOGRAM')
xlabel('interval (sec)')
text(.5*max_interval,.9*max(histogram),['CV: ' num2str(c_v)])

for k = 1:3
    c_i = bin_length*10^(k-1);                          % length of counting interval in seconds
    fano = fano_factor(spike_train,bin_length,c_i);     % variance of count over mean of count
                                                        % = 1 for a Poisson process
    text(.5*max_interval,(.4+.1*k)*max(histogram),['Fano factor for ' num2str(c_i) ' s: ' num2str(fano)])
end