% load data into time, voltage and spiketimes
x = load('extdata');
time = x.time;
voltage = x.voltage;
spiketimes = x.spiketimes;
waveformt = x.waveformt;
waveform1 = x.waveform1;
waveform2 = x.waveform2;

% indices into voltage trace of spike times:
dt = time(2) - time(1);
tinx = round(spiketimes/dt) + 1;

% plot voltage trace with detected spikes:
figure(1);
plot(time, voltage, '-b')
hold on
scatter(time(tinx), voltage(tinx), 'r', 'filled');
xlabel('time [s]');
ylabel('voltage');
xlim([0.1, 0.4])  % zoom in
hold off

% spike waveform snippets:
snippetwindow = 0.005;  % milliseconds
w = ceil(snippetwindow/dt);
vs = zeros(length(tinx), 2*w);
for k=1:length(tinx)
    vs(k,:) = voltage(tinx(k)-w:tinx(k)+w-1);
end
% time axis for snippets:
ts = time(1:size(vs,2));
ts = ts - ts(floor(length(ts)/2));

% plot all snippets:
figure(2);
hold on
plot(1000.0*ts, vs, '-b');
title('spike snippets')
xlabel('time [ms]');
ylabel('voltage');
hold off
savefigpdf(gcf, 'spikesorting1.pdf', 12, 6);

% pca:
cv = cov(vs);
[ev , ed] = eig(cv);
[d, dinx] = sort(diag(ed), 'descend');

% features:
figure(2)
subplot(1, 2, 1);
plot(1000.0*ts, ev(:,dinx(1)), 'r', 'LineWidth', 2);
xlabel('time [ms]');
ylabel('eigenvector 1');
subplot(1, 2, 2);
plot(1000.0*ts, ev(:,dinx(2)), 'g', 'LineWidth', 2);
xlabel('time [ms]');
ylabel('eigenvector 2');
savefigpdf(gcf, 'spikesorting2.pdf', 12, 6);

% plot covariance matrix:
figure(3);
subplot(1, 2, 1);
imagesc(cv);
xlabel('time bin');
ylabel('time bin');
title('covariance matrix');
caxis([-0.1 0.1])

% spectrum of eigenvalues:
subplot(1, 2, 2);
scatter(1:length(d), d, 'b', 'filled');
title('spectrum');
xlabel('index');
ylabel('eigenvalue');

savefigpdf(gcf, 'spikesorting3.pdf', 12, 6);

% project onto eigenvectors:
nx = vs * ev(:,dinx(1));
ny = vs * ev(:,dinx(2));

% clustering (two clusters):
%kx = kmeans([nx, ny], 2);
% nx smaller or greater a threshold:
kthresh = 1.6;
kx = ones(size(nx));
kx(nx<kthresh) = 2;

% plot pca coordinates:
figure(4)
scatter(nx(kx==1), ny(kx==1), 'r', 'filled');
hold on;
scatter(nx(kx==2), ny(kx==2), 'g', 'filled');
hold off;
xlabel('projection onto eigenvector 1');
ylabel('projection onto eigenvector 2');

savefigpdf(gcf, 'spikesorting4.pdf', 12, 10);

% show sorted spike waveforms:
figure(5)
subplot(1, 2, 1);
hold on
plot(1000.0*ts, vs(kx==1,:), '-r');
plot(1000.0*waveformt, waveform2, '-k', 'LineWidth', 2);
xlim([1000.0*ts(1) 1000.0*ts(end)])
xlabel('time [ms]');
ylabel('waveform 1');
hold off

subplot(1, 2, 2);
hold on
plot(1000.0*ts, vs(kx==2,:), '-g');
plot(1000.0*waveformt, waveform1, '-k', 'LineWidth', 2);
xlim([1000.0*ts(1) 1000.0*ts(end)])
xlabel('time [ms]');
ylabel('waveform 2');
hold off
savefigpdf(gcf, 'spikesorting5.pdf', 12, 6);

% mark neurons in recording:
figure(1);
hold on;
scatter(time(tinx(kinx1)), voltage(tinx(kinx1)), 'r', 'filled');
scatter(time(tinx(kinx2)), voltage(tinx(kinx2)), 'g', 'filled');
hold off;

savefigpdf(gcf, 'spikesorting6.pdf', 12, 6);

% ISIs:
figure(7);
subplot(1, 3, 1)
allisis = diff(spiketimes); 
bins = [0:0.005:0.2];
[h, b] = hist(allisis, bins);
bar(1000.0*b, h/sum(h)/mean(diff(b)))
title('all spikes')
xlabel('ISI [ms]')
xlim([0, 200.0])

subplot(1, 3, 2)
spikes1 = time(tinx(kx==1));
isis1 = diff(spikes1); 
[h, b] = hist(isis1, bins);
bar(1000.0*b, h/sum(h)/mean(diff(b)))
title('neuron 1')
xlabel('ISI [ms]')
xlim([0, 200.0])

subplot(1, 3, 3)
spikes2 = time(tinx(kx==2));
isis2 = diff(spikes2); 
[h, b] = hist(isis2, bins);
bar(1000.0*b, h/sum(h)/mean(diff(b)))
title('neuron 2')
xlabel('ISI [ms]')
xlim([0, 200.0])

savefigpdf(gcf, 'spikesorting7.pdf', 12, 6);