MATLAB: How to rapidly filter a sorted vector for elements separated by a minimum value, processing one-way from the first value

filtering separation

Using MATLAB R2012b, I have a (large) sorted vector of elements. Starting from the initial value, I need to discard all elements that are within some particular range of any previous undiscarded element. Simply looking for any set with the given minimum separation will not give the unique solution I require.
(Slow, but obvious) MWE:
function vec2 = testfun(vec,threshold)
  idx = 0.1;
  vec2 = vec;
  while not(isempty(idx));
      testdiff = diff(vec2);
      idx = find(testdiff<threshold,1);
      vec2(idx+1) = [];
  end
end
When I call this on
vec = [1,13,27,30,33,39,53,54,100];
threshold = 25;
I should only retain:
vec2 == [1,27,53,100]
but calling
testfun(vec,10)
should give me
vec2 = [1, 13, 27, 39, 53, 100]
When I call it with a script to randomly generate larger vectors:
numval = 1e5;
threshold = 25;
vec = sort(rand(numval,1).*numval*1e2);
vec2 = testfun(vec,threshold);
the profiler indicates the following breakdown within the function:
vec2(idx+1) = [];             66%

dx = find(vecdiff<thresho...  20%
vecdiff = diff(vec2);         12%
Clearly, the re-indexing and allocation of vec2 after deleting elements is slow. Yet, if I don't take those elements out (e.g., I could set them = 0 instead of []), I have to find a way to ignore negative and zero elements in the diff output, but continue to search for results below the threshold value from the previous retained element. The total runtime in the function appears to behave quadratically, although the idx = find() line begins to take up more of the total time as the length of the vector increases (517s, 46/40/13% for numval = 4e5).

Best Answer

function vec = testfun2(vec, threshold)
idx  = 1;
n    = length(vec);
keep = false(1, n);
while ~isempty(idx)
  keep(idx) = true;
  iThresh   = vec(idx) + threshold;
  idx       = find(vec > iThresh, 1);  % [EDITED]
end
vec = vec(keep);
[EDITED, 2017-02-17 08:24 UTC]:
A version which runs in linear time:
function vec = testfun3(vec, threshold)
n    = length(vec);
keep = false(1, n);
lim  = vec(1) + threshold;
keep(1) = true;
for k = 2:n
   if vec(k) > lim
      keep(k) = true;
      lim     = vec(k) + threshold;
   end
end
vec = vec(keep);
Timings under R2016b/64:
numval = 1e5;
threshold = 25;
vec = sort(rand(numval,1).*numval*1e2);
tic; r1 = testfun(vec,  threshold); toc   % Original
tic; r2 = testfun2(vec, threshold); toc   % FIND(vec > iThresh)
tic; r3 = testfun3(vec, threshold); toc   % Dull loop
isequal(r1, r2, r3)
 Elapsed time is 44.117896 seconds.
 Elapsed time is 10.090938 seconds.
 Elapsed time is 0.003627 seconds.
 1  % Results are equal
Whoooops.
Could you please cross-check this? A speedup of factor 12000 usually means, that something important has been overseen.
Related Question