% load data into time, voltage and spiketimes load( 'extdata' ); % 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 ); clf; plot( time, voltage, '-b' ) hold on scatter( time(tinx), voltage(tinx), 'r', 'filled' ); xlabel( 'time [s]' ); ylabel( 'voltage' ); xlim([0.1, 0.4]) hold off % spike waveform snippets: w = ceil( 0.005/dt ); vs = zeros(length(tinx), 2*w); for k=1:length(tinx) vs(k,:) = voltage(tinx(k)-w:tinx(k)+w-1); end ts = time(1:size(vs,2)); ts = ts - ts(floor(length(ts)/2)); % % plot snippets: figure( 2 ); clf; hold on for k=1:size(vs,1) plot( ts, vs(k,:), '-b' ); end xlabel( 'time [ms]' ); ylabel( 'voltage' ); hold off % pca: cv = cov( vs ); [ ev , ed ] = eig( cv ); [d,dinx] = sort( diag(ed), 'descend' ); figure( 3 ); clf; subplot( 4, 2, 1 ); imagesc( cv ); xlabel( 'time bin' ); ylabel( 'time bin' ); title( 'covariance matrix' ); caxis([-0.1 0.1]) % spectrum of eigenvalues: subplot( 4, 2, 2 ); scatter( 1:length(d), d, 'b', 'filled' ); xlabel( 'index' ); ylabel( 'eigenvalue' ); % features: subplot( 4, 2, 5 ); plot( 1000.0*ts, ev(:,dinx(1)), 'r', 'LineWidth', 2 ); xlabel( 'time [ms]' ); ylabel( 'eigenvector 1' ); subplot( 4, 2, 6 ); plot( 1000.0*ts, ev(:,dinx(2)), 'g', 'LineWidth', 2 ); xlabel( 'time [ms]' ); ylabel( 'eigenvector 2' ); % project onto eigenvectors: nx = vs * ev(:,dinx(1)); ny = vs * ev(:,dinx(2)); %scatter( nx, ny, 'b', 'filled', 'MarkerEdgeColor', 'white' ); % 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; subplot( 4, 1, 2 ); scatter( nx(kx==1), ny(kx==1), 10.0, 'r', 'filled' ); %, 'MarkerEdgeColor', 'white' ); hold on; scatter( nx(kx==2), ny(kx==2), 10.0, 'g', 'filled' ); %, 'MarkerEdgeColor', 'white' ); hold off; xlabel( 'projection onto eigenvector 1' ); ylabel( 'projection onto eigenvector 2' ); % show sorted spike waveforms: subplot( 4, 2, 7 ); hold on kinx1 = find(kx==1); for k=1:length(kinx1) plot( 1000.0*ts, vs(kinx1(k),:), '-r' ); end 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( 4, 2, 8 ); hold on kinx2 = find(kx==2); for k=1:length(kinx2) plot( 1000.0*ts, vs(kinx2(k),:), '-g' ); end 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 % spike trains: figure( 1 ); hold on; scatter( time(tinx(kinx1)), voltage(tinx(kinx1)), 10.0, 'r', 'filled' ); scatter( time(tinx(kinx2)), voltage(tinx(kinx2)), 10.0, 'g', 'filled' ); hold off; % ISIs: figure(4); clf; subplot(1, 3, 1) allisis = diff(spiketimes); hist(allisis, 20); title('all spikes') xlabel('ISI [ms]') subplot(1, 3, 2) spikes1 = time(tinx(kinx1)); isis1 = diff(spikes1); hist(isis1, 20); title('neuron 1') xlabel('ISI [ms]') subplot(1, 3, 3) spikes2 = time(tinx(kinx2)); isis2 = diff(spikes2); hist(isis2, 20); title('neuron 2') xlabel('ISI [ms]') pause