clear
clf
dt = 0.001;
tmax = 3;
t = 0:dt:tmax;
r0 = 20;
rf = 40;
Ntrials = 10;

r = r0+(rf-r0)/tmax*t;  % generate a rate that ramps uniformly with time

plot(t,r)

spikes = zeros(Ntrials,length(t));

for trial = 1:Ntrials           % do a bunch of trials with different random spike distributions
    for i = 1:length(t)
        if rand(1) < dt*r(i)        % using a Poisson spike generator
            spikes(trial,i) = 1.0;  % include a spike if the rate is high enough
        end
    end
end

figure()
for trial = 1:Ntrials
    subplot(Ntrials,1,trial)
    plot(t,spikes(trial,:))
end


% Now do trial averaging of spike trains 

sum1 = sum(spikes)/(dt*Ntrials);


sum(sum1)

smoothsigma = 0.050;
taumax = 5*smoothsigma;
tau = -taumax:dt:taumax;
smooth = dt*exp(-(tau.*tau)/(2.0*smoothsigma*smoothsigma))/(sqrt(2.0*pi*smoothsigma*smoothsigma));

figure()
plot(tau,smooth)
norm = sum(smooth,2)        % check total normalization is 1

sum2 = zeros(size(t));

jrange = taumax/dt;

for i = 1:length(t)
    jmin = max(i-jrange,1);
    jmax = min(i+jrange,length(t));
    for j = jmin:jmax
        sum2(j) = sum2(j) + sum1(i)*smooth(j-i+jrange+1);
    end
end

figure()
subplot(2,1,1)
plot(t,sum1)

subplot(2,1,2)
plot(t,sum2,'r')