7

I'm working on a system with hard real-time constraints in c++ and I need a very fast way to calculate the rolling/moving/streaming median of a set of numbers of size N=100 to 300. Normally this size would be trivial but in this case the algorithm will run about 1000/2000 times per 0.1ms

Every calculation a single value in range (0-1) will be added so (depending on the container) the previous values will already be (weakly) sorted.

Requirements:

No dynamic memory allocation as the median window will be fixed size and any allocation could take too much time

FIFO behaviour in order to remove the last added and insert a new value

Current Approach

Currently I'm considering a min-max heap rolling median approach where you keep a reference to the oldest element and every element has a reference to the next-oldest element. However, I'm unsure how this would work as you would need to remove that element from any position in the heap and it could also be in either of the two heaps.

  • Apologies, I have deleted the other post as people suggested it would fit better on software engineering – Bluefarmer Apr 13 '18 at 12:58
  • @Bluefarmer I don't see why this isn't a good question for this forum. Just delete the other post on stack/flow. Can't say I really know a good answer to you're question though! Tough requirements. May I ask why you must track the median so often? Is there something that is actually polling this value as fast as you update it? – Christopher Apr 13 '18 at 12:58
  • @Christopher Yeah, I'm using the median to determine how to modify data in the frequency spectrum of an audio signal and re-synthesise it back to the time domain after modifications so it's being polled every time a new value is given. – Bluefarmer Apr 13 '18 at 13:05
  • Does a histogram approach work? Quantize sample values into a histogram bin; for each new sample, add one to the corresponding histogram bin; for each sample exiting the window, subtract one. This approach requires memory for both the histogram and the sample window; the latter to allow constant-time processing of each exiting sample without depending on a time-shifted source. The acceptance of this approach depends on whether quantization (or a scheme of it) exists for your application needs. I have a feeling that this question may be suitable for dsp.stackexchange.com – rwong Apr 13 '18 at 13:40
  • 3
    Did you actually profile the more or less naive approach of keeping a sorted array by and insert/remove the fifo values by shifting (with a binary search for the lookup)? By which order of magnitude is that too slow? "N=100 to 300" looks to my like a range where this might still be quick enough for many real-time applications. – Doc Brown Apr 13 '18 at 16:04
1

If I were going to do this, I'd probably use a balanced tree (e.g., AVL or red-black) where each node also keeps track of the size of its left sub-tree (and you keep track of the overall size). This would act as an index into a circular buffer storing the data itself.

So, when a new item arrives, you find the oldest value in the circular buffer. Search for that node in the tree and delete it [O(log N) complexity]. As you're descending the tree to delete it, you update the count in each node, to signify which sub-tree is getting smaller.

Overwrite that node in the circular buffer with the new value. Insert a node into the tree for the new value [also O(log N) complexity]. Again, as you descend the tree to insert the new value, you update the counts in the nodes to signify the new sub-tree size.

You find the median via the tree--start from the root node. Check the size of the left sub-tree to determine whether the median is in the left or right sub-tree. Travel down the tree until you reach the median. Since the tree depth is proportional to log N, this operation is also O(log N).

That leaves one minor detail: the tree will typically use dynamic allocation. You can, however, keep this from being a problem pretty easily: pre-allocate a block of space for all the tree nodes you need. Build them into (for example) a linked list of free nodes. When you need to allocate a node, you always just grab the one from the head of the list. When you need to free a node, you just add it to the front of the list.

Once your tree is "full", you'll basically just be adding a node to the list, and immediate re-using it to create the new node. The "linked list" will basically just act be a pointer to that one node while its free.

1

This is a refinement of Jerry Coffin's idea.

  • Use a nearly balanced tree, where all nodes reside directly in the circular buffer.
  • Initialize it with dummy values, so that the size stays constant all the time.
  • Use left, right and parent pointers (or indexes) in each node.
  • Store and manage rankin each node.
  • Instead of removing an inserting an element, do an inplace replace by simply changing the value.
  • After the replacement, the tree is usually no more ordered.
  • Use rotations to restore the order, such that
    • There's always at most one node violating the order.
    • The tree keeps its RB or AVL property.

Doing this is surely faster than removing and insering nodes. When you're lucky, the node moves only a few times. With a probablity of 50%, the old and the new element belong to the same half-tree, etc.

The rotations are to be devised and this may be a bit complicated, as you need to restore the order and also to keep the tree balanced at the same time.

I guess, you can gain a factor of four this way, but that's just a guess.

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.