Efficient way to find closest spacing of a "comb" array

Hello!

I found that the most inefficient part of some code is a nested loop that is not vectorized; I'm not sure how to vectorize it though at the innermost loop.

The problem is to find, in a "comb" array, the maximum number of "teeth" within a given interval. If you don't know what a comb array is, it is mostly zeros with occasional values of 1 (the teeth). This is then convolved with a waveform, but you don't really need to know that part. It is important to know how many teeth are close together.

The comb array is multidimensional, but really only the one dimension is relevant, the one looped over in the innermost loop. The code is:

% Check for number of delays within 8 ms for comb function
delay_interval = 7;  % One ms less than 8ms
length_of_comb = length(comb_array);
num_checks = length_of_comb - delay_interval;
%multiWaitbar('Part 2',0);
for n = 1:num_row_delays
    %multiWaitbar('Part 2',n/num_row_delays);
    for m = 1:num_hole_delays
        charge_per_8ms_delay = 0; %start out at zero
        for p = 1:num_checks

              chunk_to_check = comb_array(n,m,p:p+delay_interval);
              ctc = sum(chunk_to_check);
              if ctc > charge_per_8ms_delay
                  charge_per_8ms_delay = ctc;
              end
          end
          chg_per_delay_array(n,m) = charge_per_8ms_delay;
      end
  end
  %multiWaitbar('Part 2','Close');

Any thoughts on how to make this run faster (i.e., smarter!) are welcomed!

Thanks.

Doug Anderson

Answer 1

[Guillaume]

What you're doing is calculating a sum in the third dimension over a sliding window of fixed size. This is easily achieved with a convolution with a vector of ones. To demonstrate it with a simpler 2d matrix:

comb_array = [1 1 0 0 0 0 0 1 0 0 0 0 0 0 0 0 1
              0 0 0 0 0 0 1 1 0 1 0 1 0 0 0 0 0];
ctc = conv2(comb_array, ones(1, 7), 'valid') %sliding sum over 7 columns

So, you can replace all of your loops with:

delay_interval = 8;  for delay within 8 ms. Don't subtract one anymore
ctc = convn(comb_array, ones(1, 1, delay_interval), 'valid'); sliding sum over delay_interval pages

With that sum, you're then simply calcuting its maximum in the third dimension, so:

charge_per_delay = max(ctc, [], 3);

Finally, I'll note that in your original code, with the line

length_of_comb = length(comb_array);

you mean to get the size of the third dimension. If that size ever happens to be smaller than the number of columns or rows of your comb array, then it will fail, since length gives you the size of the largest dimension. It is always better to be explicit about which dimension you want, so it should have read:

length_of_comb = size(comb_array, 3);

[Douglas Anderson]

OK, thanks, that makes a lot of sense!

Answer 2

[Jan]

As far as I can see, this is a "moving sum" with a threshold.

B = ones(1, delay_interval+1);
...
  % Replacement of the "for p" loop:
  value = filter(B, 1, squeeze(comb_array(n, m, :)));
  chg_per_delay_array(n,m) = value(find(value, 1, 'last'));