function [pm , ng, nn] = glitch(p,factor, pm)
% function [pm , ng, nn] = glitch(p,factor, pm)
%
% Locates and counts glitches in a vector, using a criterion
% based on the absolute values of 3 consecutive second differences
% exceeding a threshold.
%
% inputs:
%   p       column or row vector, typically a time or space series
%   factor  scalar; threshold is factor times the mean absolute value
%           of the second differences. A good starting value is 2.
%   pm      (optional) an initial mask with zeros at bad points in p.
%           If not supplied, it will be generated as ~isnan(p).
% outputs:
%   pm      input mask pm, updated with zeros at newly located glitches.
%   ng      count of new glitches found.
%   nn      count of original glitches.
%

% factor < 5; typical value factor = 2;


[nr, nc] = size(p);
np = max(nr, nc);
if np == nr,
   is_col = 1;
   p = p(:)';
else
   is_col = 0;
end

pi   = 1:np;
if nargin == 2,
   pm   = ~isnan(p);
end
nn   = np - sum(pm);  % Number of NaNs in p.
done = 0;

dd_m = mean(abs(diff(diff(p(pm)))));
if (dd_m < 1.0e-6), dd_m = 1.0e-6;  end
tol = factor*dd_m;

while (done == 0),

  pddb = [0 abs(diff(diff(p(pm)))) > tol 0];

  pi_m = pi(pm);

  i_b0 = find(diff([0 pddb]) == 1);

  if (length(i_b0) > 0)

     i_b0(i_b0+2 > sum(pm)) = [];

     i_bb = find(pddb(i_b0+1) == 1 & pddb(i_b0+2) == 1);

     i_del = i_b0(i_bb) + 1;

     if (length(i_del) > 0),
         pm(pi_m(i_del)) = zeros(1,length(i_del));
     else
        done = 1;
     end
   else
     done = 1;
   end

   %tol = 2*tol;

end

ng = sum(~pm) - nn;  % Number of new glitches found.

if is_col == 1,
   pm = pm(:);  % I am not sure if there is any point in ensuring
end             % that the mask has the same orientation as p.

return

