Thresholding a signal using hysteresis to reduce noise

Question

I'm using a simple thresholding for noisy data to detect zeros/edges in signals. Since the noise can be pretty strong, I'm using hysteresis to improve the results. While this helps a lot, it also slows down a lot. Is there a way to improve this, maybe even a better way to calculate this? Currently I'm using the straight-forward loop approach:

% generate a signal
t = linspace(0, 10, 1000);
y = sin(2 * pi * t); 

% constants
threshold = 0;
hyst = 0.1;

% thresholding
yth = zeros(size(y));
for i = 2:length(y)
    if yth(i - 1) > 0.5 
        yth(i) = y(i) > (threshold - hyst);
    else
        yth(i) = y(i) > (threshold + hyst);
    end 
end

In comparison to yth = y > threshold this is much slower.

Answer 1

An improvement (25% time reduction; only for large inputs; on R2017b) can be gained by precomputing both possibilities for yth(i):

function yth = idea1()
  yth = false(size(y)); % change #1 - `false` vs `zeros`
  c1 = y > (th - hyst); % change #2 - precompute c1, c2
  c2 = y > (th + hyst);

  for k = 2:numel(y)
      if yth(k - 1) 
          yth(k) = c1(k);
      else
          yth(k) = c2(k);
      end 
  end
end

Usually, a big improvement can be gained by vectorization. To vectorize this problem we must understand the switching logic, so let us put it into words:

• Start from Mode 2.
• Mode 1: Take a value from c1. As long as c1 contains true values - take those. When a false value is encountered switch to Mode 2.
• Mode 2: Take a value from c2. As long as c2 contains false values - take those. When a true value is encountered switch to Mode 1.

So if we could find the transition locations, we are practically done.

After some trial and error, while I was unable to get rid of the loop in a way that improves performance, I did reach the conclusion that it was possible (idea2). Furthermore, looking at the RLE-encoded correct yth, I came up with some pretty good approximations for it (idea3 and idea4) - though these will require tweaking for other inputs.

Perhaps somebody could use it to find a more clever implementation with fewer redundant computations. My full code is provided below. The RLE encoding algorithm was adapted from this answer, and RLE decoding from here.

function q48637952
% generate a signal
t = linspace(0, 10, 1000).';
y = sin(2 * pi * t); 

% constants
th = 0;
hyst = 0.1;

%% Comaprison:
% Correctness:
R = {originalIdea(), idea1(), idea2()};
assert(isequal(R{:}));
% Runtime:
T = [timeit(@originalIdea,1), timeit(@idea1,1), timeit(@idea2,1)];
disp(T);

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function yth = originalIdea()
  yth = zeros(size(y));
  for i = 2:length(y)
      if yth(i - 1) > 0.5 
          yth(i) = y(i) > (th - hyst);
      else
          yth(i) = y(i) > (th + hyst);
      end 
  end
end

function yth = idea1()
  yth = false(size(y)); % change #1 - `false` vs `zeros`
  c1 = y > (th - hyst); % change #2 - precompute c1, c2
  c2 = y > (th + hyst);

  for k = 2:numel(y)
      if yth(k - 1) 
          yth(k) = c1(k);
      else
          yth(k) = c2(k);
      end 
  end
end

function yth = idea2()
  c = [y > (th - hyst), y > (th + hyst)];
  % Using Run-length encoding:
  [v1,l1] = rlEncode(c(:,1));
  [v2,l2] = rlEncode(c(:,2));
  rle = cat(3,[v1 l1 cumsum(l1)],[v2 l2 cumsum(l2)]);
  % rle(:,:,2) is similar to rle(:,:,1) with small changes:
  %  - col1 is circshifted by 1 and 
  %  - col2 has the 1st and last elements switched
  newenc = reshape([rle(1:2:end,3,2), rle(1:2:end,3,1)].', [], 1);
  yth = rlDecode(rle(:,1,2), [newenc(1); diff(newenc(1:end-1))]);
end

function yth = idea3()
  % Approximation of yth, should be almost indistinguishable with high resolution
  yth = [0 0 repmat(repelem([1,0],50), 1, ceil(numel(y)/(2*50)) )].';
  % The amount of zeros at the beginning as well as the value "50" are problem-specific 
  % and need to be computed using signal-specific considerations
  yth = yth(1:1000);
end

function yth = idea4()
  % Another approximation
  yth = circshift(y > th, 1);
  % The value "1" is problem-specific
end

end % q48637952

%% RLE (de)compression:
function [vals, lens] = rlEncode(vec)
  J = find(diff([vec(1)-1; vec(:)]));
  vals = vec(J);
  lens = diff([J; numel(vec)+1]);
end

function vec = rlDecode(vals, lens)
  l = cumsum([ 1; lens(:) ]);
    z = zeros(1, l(end)-1);
    z(l(1:end-1)) = 1;
    vec = vals(cumsum(z));
end