%% http://journals.plos.org/ploscompbiol/article?id=10.1371/journal.pcbi.1003258

clear;

dt = 0.001;
s = zeros(1000,1);
s(300:900) = 100*sind((0:600)/3.4);

o = zeros(1000,1);
V(1) = 0;
gam = 0.5;
lam = 0.01;

for i = 1:length(s)-1,
    V(i+1) = V(i) - 1/lam.*dt.*V(i) + gam.*dt.*s(i) - gam.^2.*o(i);
    if V(i+1)>gam.^2./2, o(i+1) = 1; end
end

ker = exp(-(0:1000)./lam.*dt); ker = ker./sum(ker)./dt;
shat = conv(ker,o).*gam;

figure;
hold on;
plot(s,'k:')
plot(100*V,'r-')
o(find(o==0)) = NaN;
plot(o,'b.')
plot(shat(1:end-1001),'g-')
xlabel('Time (ms)')
legend('stimulus S','potential V','spikes','reconstructed S')