Blogs

tl;dr: Step dection is actually a special case of convolution that can be made really fast. See the Godbolt example

Overview

 Step detection in one-dimensional signals is an interesting problem. In the ideal case you can just look when the difference between two samples is greater than 0 and done like dinner. But reality isn't so simple sadly, as almost all signals in real life are noisy. One way this is often done is via convolution. I encountered this Stack Overflow answer that introduced me to that. I have previously used convolution for filtering in my digital signal processing class, but I had not considered it for other use cases.

 But it turns out this has a downside. Convolution is extremly slow. It's complexity is $𝒪(nm)$, where $n$ and $m$ are the lengths of the input signals $g\left(x\right)$ where $x=[0,n]$ and $h\left(x\right)$ where $x=[0,m]$. For step detection $m=\mathrm{2n}$, and $h\left(x\right)=\{\begin{array}{cc}1& \text{if }x=[0,n]\\ -1& \text{if }x=(n,2n]\end{array}$. This makes it very convinent to find where steps occur as the output signal is a sawtooth signal, and where the second-derivitive is non-zero is where steps occur. An image from the Stack Overflow post shows this clearly.

Step detection of signal. Blue signal is input. Orange signal is output from convolution. Vertical line is where step occurs.

 For the general case of convolution this is a $𝒪(n^2)$ calculation which sucks for large signals, for example from an oscilliscope. It is also unfortunate that the accelerator packs for Python's Numpy pacakge vary widely in their convolution speeds. Numpy accelerated with OpenBLAS is several hundred times slower than Numpy built with Intel's MKL, which is tragic as I would prefer the open-source alternative. Shown below is a function that will calculate the convolved signal.

Python

import numpy as np

def step_detect(data: np.ndarray) -> np.ndarray:
    """Convolves input 1-D signal with step detection kernel
    """
    preprocess = data - np.average(data)
    step = np.hstack((np.ones(data.size), -1 * np.ones(data.size)))
    return np.convolve(preprocess, step, mode="valid")

Then timing it with OpenBLAS:

Python

>>> from timeit import timeit 
>>> test_data = np.random.rand(15_000_000)
>>> timeit(lambda: step_detect(test_data), number=15_000_000) / 15_000_000
# This is me a year later, I don't want to rerun these but it was around 25 minutes

and with Intel MKL:

Python

>>> from timeit import timeit
>>> test_data = np.random.rand(15_000_000)
>>> timeit(lambda: step_detect(test_data), number=15_000_000) / 15_000_000
# This is me a year later, I don't want to rerun these but it was around 120 seconds

While some people may dismiss this as a microbenchmark that is not indicative of real-world performance, this actually is a real world situation. These small numbers do scale to larger input signals sadly. But what if there was a way to make it much much faster?

How to make it much much faster

 As the title spoiled it can be made to be extremly fast. While I was sitting in the queue to Amazon Game Studio's New World on release day I had plenty of time to not play the video game. During that time I though about convolution, and this thought process brought me to realizing that all that step detection is doing is summing the signal up to a certain index, then subtracting the rest of the signal. It may be easier to imagine discrete convolution as a "sliding dot product", where the input signal slides along the convolution kernel and at each index the output vector is built from the dot product of the signals.

Note
	All below was written a year later cause I forgot about this.
				

 I sort of forgot about this blog article, and all of the above was written about 1.5 years ago. But anyways with the above knowledge, you can take the unoptimized implementation, that with the input I was using it with (noisy oscilliscope signals) you can take the runtime from about 2 minutes to somewhere in the tens of milliseconds.

C

int fast_conv(const size_t len, const float array[], float output[])
{
  if (!array || !output) {
    return -1;
  }
  size_t i = 0;
  float accum = 0.0f;

// Calculate the convolution
#pragma omp simd reduction(+ : accum)
  for (i = 0; i != len; i++) {
    accum += array[i];
  }

  for (i = 0; i != len; i++) {
    output[i] = accum;
    accum -= 2.0f * array[len - i - 1];
  }
  output[i] = accum;

  return 0;
}

And the Rust equivalent

Rust

pub fn fast_conv(input: &[f32]) -> Vec {
    let mut buffer = Vec::with_capacity(input.len() + 1);

    let mut accum = input.iter().sum();

    for x in input.iter().rev() {
        buffer.push(accum);
        accum -= 2.0 * x;
    }
    buffer.push(accum);
    
    buffer
}

I didn't want to make a real Python module at the time, so I built a module around it using ctype and NumPy ctypes. In the end I sped up a program that took a couple days to run to less than an hour.