Despikes neural data using linear interpolation

despike.m

function dataDespiked = despike(data,Fs,spikeTs,spikeWidth)
% data: raw neural recording data
% Fs: sampling frequency
% spikeTs: spike timestamps in seconds
% spikeWith: width of average spike in ms to replace

% display a plot of the replaced spikes at the end
doDebug = true;

% convert spikeTs (in s) to samples
spikeSamples = round(spikeTs * Fs);
% convert spikeWidth (in ms) to samples
rmSamples = round(Fs * (spikeWidth/1000) / 2);

f = waitbar(0,'setting up inerpolation...');
% fill in all data points around spikes with NaN
dataNaN = double(data);
for iSpike = 1:numel(spikeTs)
    waitbar(iSpike / numel(spikeTs),f);
    nanRange = spikeSamples(iSpike) - rmSamples:spikeSamples(iSpike) + rmSamples - 1;
    if nanRange(1) > 0 && nanRange(end) < numel(dataNaN)
        dataNaN(nanRange) = NaN;
    end
end
waitbar(1,f,'interpolating...');
% interpolate NaN data using linear method
% see help > interp1 to change method
v = 1:numel(dataNaN);
xq = v;
v = v(~isnan(dataNaN));
dataNaN = dataNaN(~isnan(dataNaN));
dataDespiked = interp1(v,dataNaN,xq,'linear');
close(f);

% plot original and despiked data and show spike locations
if doDebug
    showSec = 0.5;
    dataSnip = double(data(1:round(Fs * showSec)));
    dataNaN_interpSnip = double(dataDespiked(1:round(Fs * showSec)));

    figuree(1000,300);
    plot(dataSnip);
    hold on;
    plot(dataNaN_interpSnip,'-k');

    tsSnip = spikeTs(spikeTs < showSec);
    tsSnip_samp = round(tsSnip * Fs);
    plot(tsSnip_samp,dataSnip(tsSnip_samp),'rx');
    xlim([1 numel(dataSnip)]);
    title([num2str(showSec),'s snip of data']);
    legend('original','despiked');
end