%Bruce's EOD data
clear all

load 'c:\my documents\bigdata\bruceC\data'

dt = .001;
t = 0:dt:20;
TrainTimes1 = data(:,1);
Spikes1 = zeros(size(t));
Spikes1(fix(TrainTimes1/dt)) = 1 ;

figure(1)
clf
np = 4;

subplot(np,1,1)
plot(t,Spikes1)
set(gca,'ylim',[0 2])
xlabel('Time')
Ylabel('Whole train')

burst1index = fix(2.8/dt):fix(2.8/dt)+fix(.7/dt) ;
burst1 = Spikes1(burst1index) ;
subplot(np,1,2)
plot(t(burst1index),burst1)
set(gca,'ylim',[0 2])
xlabel('Time')
Ylabel('Ref Burst')

burst2start = 13.6 ; %13.6--18.7--4.35--5.3
for i=1:400
    burst2index=fix(burst2start/dt)+i:fix(burst2start/dt)+i+fix(.7/dt) ;
    burst2 = Spikes1(burst2index) ;
    
    for s=1:14
        TimeScale(s) = dt * 1.5^(s-1);
        %calc vanRossum dist
        tc = TimeScale(s) ; % Sec
        tk = 0:dt:3*tc;
        ekernel = exp(-tk/tc);
        cSpike1 = conv(ekernel,burst1);
        cSpike2 = conv(ekernel,burst2);
        dRossum(s) =  dt/tc * sum((cSpike1-cSpike2).^2) ;
    end
    Dtotal(i) = sum(dRossum) ;
end

subplot(np,1,3)
plot(Dtotal)
xlabel('Sample offset')
Ylabel('Summed Dist')

subplot(np,1,4)
plot(burst1)
hold on
[Dmin,DminIndex] = min(Dtotal);
burst2index=fix(burst2start/dt)+DminIndex:fix(burst2start/dt)+DminIndex+fix(.7/dt) ;
burst2 = Spikes1(burst2index) ;
plot(burst2+.9,'r')
xlabel('Sample within burst')
Ylabel('Spikes')
set(gca,'ylim',[0 2])


