function [ despiked ] = despike( raw, theshold, debug )
% function [ despiked ] = despike( raw, theshold )
%
% Function despike takes an array and converts spikes to NaNs.
%   raw:  signal to be despiked
%   threshold:  size of spike to be removed (10 a good starting point)
%   debug:  1 if you want a plot, 0 otherwise

n = length(raw);
threshold = 10;

dif_r  = raw(1:n-1) - raw(2:n);  % Compute the difference of each value and the following value
dif_2r = raw(1:n-2) - raw(3:n); % Compute the different of each value from the second following value

% Find a reasonable expectation of difference between adjacent points.
cond_dif_r = abs(dif_r(find(~isnan(dif_r))));
r_med = median(cond_dif_r);

% Test to see if the one point difference is much greater than the two
% point difference (a sign of a spike) AND the one point difference is also
% greater than some multiple (threshold) of the median expectation.  This
% second test guards against the case when the two point difference is
% really small by chance, and the one point difference looks really large
% in comparison, but isn't really a spike.
spike = find(abs(dif_r(1:n-2)./dif_2r) > threshold & abs(dif_r(1:n-2)/r_med) > threshold);

if debug
    plot(raw, 'b-'); hold on; plot(spike+1, raw(spike+1), 'r*');hold off
end

% cond_dif_r = abs(dif_r(find(~isnan(dif_r))));
%
% r_med = median(cond_dif_r);
%
% length(find(cond_dif_r > 100*r_med))

raw(spike+1) = NaN;

despiked = raw;

