removeSpikes.m

function [spike,seq,Data,GPS,tBox,tHead,stat,spikeData] = removeSpikes(seq,Data,GPS,tBox,tHead,stat,HZ)

% [spike,seq,Data,GPS,tBox,tHead,stat,spikeData] =
% removeSpikes(seq,Data,GPS,tBox,tHead,stat)
%
% optionally, run spike = removeSpikes(seq) to find spikes due to a glitch
% in the NIMS systems that causes a spike in the data, accompanied by the
% following conditions:
% 1) two consecutive breaks in the sequence number continuity such that the
%    sum total of the gaps is 255:
%    (gap.ind(j) - gap.ind(j-1) == 1) && (gap.len(j) + gap.len(j-1) == 255) 
% 2) status indicator shows an error: (stat(gap.ind(j)) ~= 129
% The spike is then located at index j.
% 
% Here, we remove the spike from all the data streams, replacing it in the
% sequence string by a regular 1 sec gap.
%
% Optional output spikeData are the data that have been removed on output.
% E.g., for 8 Hz data, Data(:,8*(spike(4)-1)+1:8*spike(4)) will contain the
% block of data that remains in the data stream while spikeData(:,8*3+1:8*4)
% will contain the corresponding deleted data block.
%
% @ Anna Kelbert, 11 Aug 2020.

% find all possible gaps based on sequence numbers
ind = find(diff(seq)~=1);
spikeData = [];

spike = []; j = 1;
for i = 2:length(ind)
    if (ind(i) - ind(i-1) == 1) && (seq(ind(i)+1) - seq(ind(i)-1) - 1 == 256) && stat(ind(i) ~= 129)
        spike(j) = ind(i); j = j+1;
    end
end

if nargin == 1 || isempty(spike)
    return
end

% by default, assume 1 Hz data
if nargin < 7
    HZ = 1;
end

% count true duplicates
nspike = length(spike);

% find the spike blocks in the data stream
data_spike = zeros(HZ,length(spike));
for i = 1:length(spike)
    data_spike(:,i) = (HZ*(spike(i)-1)+1:HZ*spike(i));
end
data_spike = reshape(data_spike,1,HZ*length(spike));

% an optional output used for keeping track of spiky Data values
if nargout >= 8
    spikeData = Data(:,data_spike);
end

% remove them from the data
if nspike > 0

    ind = 1:length(seq);
    ind = setdiff(ind,spike);
    
    data_ind = 1:(HZ*length(seq));
    data_ind = setdiff(data_ind,data_spike);

    Data = Data(:,data_ind);
    GPS = GPS(ind);
    tBox = tBox(ind);
    tHead = tHead(ind);
    stat = stat(ind);
    seq = seq(ind);

    disp(['Deleted ' num2str(nspike) ' spiky data blocks from ' num2str(HZ) ' Hz data']);
end