Spike Averaging in Matlab

% This example uses the xippmex('spike') command to retrieve spikes from
% the Ripple Grapevine processor and average them based on the units
% sorted via the Trellis Spike Scope application.
 
% Clean the world
close all; fclose('all'); clc; clear all;
 
% Initialize xippmex
status = xippmex;
if status ~= 1; error('Xippmex Did Not Initialize');  end
 
% Spike Setup
spk_ch = 1;                                 % channel to read spikes from
sort_colors = {'k','m','y','g','b'};        % unsorted (0), 1, 2, 3, 4
sorted_units = [3 4];                       % sorted units to show in plot
n_sorted = length(sorted_units);            % number of sorted units to plot
spk_samples = 52;                           % number of samples in spike segment
spk_t = 1:spk_samples;                      % spike sample vector
spk_buffer_sz = 25;                         % spikes to plot and average
spk_buffer = cell(n_sorted,1);
for i = 1:n_sorted
    spk_buffer{i} = NaN(spk_buffer_sz,spk_samples);
end
spk_average = NaN(n_sorted,spk_samples);
 
% Plot Setup
FR = 30;                                    % Hz, plot update framerate
spk_plot = figure; hold on;
for i = 1:n_sorted
    h_buffer{i} = plot(spk_t,spk_buffer{i},'color',[0.8 0.8 0.8]);
    uistack(h_buffer{i},'bottom');
    h_mu{i} = plot(spk_t,spk_average(i,:),'color',sort_colors{sorted_units(i)+1},'linewidth',2);
end
y_range = [-220 180];                       % uV
x_range = spk_t([1 end]);                   % samples
ylim(y_range); ylabel('\muV');
xlim(x_range); xlabel('Sample Index');
title(sprintf('Channel #%d, last %d spikes shown',spk_ch,spk_buffer_sz));
hold off;
 
% Clear Spike Buffera
[c,t,w,u] = xippmex('spike',spk_ch); % count, timestamps, waveforms, sorted unit
 
% Main Collection/Update Loop while figure is open
while isa(spk_plot,'handle') && isvalid(spk_plot)
 
    % grab spikes from Grapevine processor
    [c,t,w,u] = xippmex('spike',spk_ch); % count, timestamps, waveforms, sorted unit
 
    if c > 0 % update plot
 
        % loop through sorted units
        for i = 1:n_sorted
            spk_idx = find(u{1}==sorted_units(i));
            count = length(spk_idx);
            if count >= spk_buffer_sz % fill buffer with most recent spikes
                spk_buffer{i} = w{1}(spk_idx(end-spk_buffer_sz+1:end),:);
            else % replace oldest spikes with newest spikes
                spk_buffer{i} = circshift(spk_buffer{i},-count);
                spk_buffer{i}(end-count+1:end,:) = w{1}(spk_idx,:);
            end
            spk_average(i,:) = nanmean(spk_buffer{i});
            for j = 1:numel(h_buffer{i})
                set(h_buffer{i}(j),'ydata',spk_buffer{i}(j,:));
            end
            set(h_mu{i},'ydata',spk_average(i,:));
        end
        drawnow
    end
 
    % Wait until next frame update
    pause(1/FR);
 
end
disp('Exiting main loop');
 
% Close xippmex after exiting loop
disp('Closing xippmex connection');
xippmex('close');