%Generate single spiketrain

rand('state',sum(100*clock));

nsec=1 ;
T=1;
deltat=0.0001;
r=500;
p=r*deltat;
numbins=round(T/deltat);
spiketrain=round(rand(1,numbins) + (p-1/2)) ;

figure;
imagesc(spiketrain)
xlabel('time (ms)','FontSize',18)
title('spike raster plot','FontSize',18)
set(gca,'FontSize',16)



%Generate many "trials" of spiketrains
numtrials=100000;
spiketrain=round(rand(numtrials,numbins) + (p-1/2));

figure;
imagesc(spiketrain)
xlabel('time','FontSize',18)
ylabel('trial')
set(gca,'FontSize',18)


%Compute the average spike rate, and standard deviation
rate_per_trial=1/T * sum(spiketrain,2)
mean_rate_per_trial = mean(rate_per_trial)
std_dev_rate_per_trial = std(rate_per_trial)


spikesPerTrial = sum(spiketrain,2);
Fano = var(spikesPerTrial) / mean(spikesPerTrial);

