# Efficient implementation of rolling MAD/Hampel filter

I have implemented the moving median absolute deviation (moving MAD) and it seems like bit-exact to Matlab's implementation. Nevertheless, I am sure that it is not efficient.

The usual median filter should be implemented with a double heap. The moving MAD can not be implemented this way since the absolute deviation of each element in the window can be extremely modified with any modification of the median. This forces using the quick-select for every sample, which makes the calculation very long...

Somehow, Matlab managed to implement moving MAD with computation time only twice the median filter. This suggests that they somehow managed to use a double median filter. Any ideas on how to implement it?

Edit: After a couple of years I think I have found the answer and implemented it in the attached code. Please find it below

• While I have written a sliding $\min()$ or sliding $max()$ algorithm that has complexity proportional to $\log_2(B)$ where $B$ is the buffer length, and I understand how to split that into two simultaneous sliding min and max for the top and bottom halves of the data (and the sliding median would be the average of the max of the bottom half buffer and the min of the top half buffer), i have never written the program that is an efficient sliding median. Sep 16, 2020 at 0:59
• @robertbristow-johnson most popular implementation of the sliding-median is the median-heap. I think that there are many implementations that you can find on the web... Sep 17, 2020 at 6:50
• You could also review the source code for the median filter implemented in python scipy.ndimage.median_filter Oct 26, 2020 at 14:09
• @Dan Boschen , I was looking for an efficient implementation of the median filter - which I have. I was looking for an efficient implementation of the rolling MAD, which is closely related to the Hampel filter. Unlike the median filter, the Hampel filter modifies the signal in only at the locations which are marked as spike Oct 26, 2020 at 17:37

The median-heap algorithm as you probably know keeps two heaps.

One possible efficient solution would be to implement a heap based on a self-balancing binary search tree, that on any operation keeps track of sum and size of each subtree.

Then you can implement an algorithm that will find the median(abs(X - median(X))), as well as the node that splits the nodes in median(X) - median(X - median(X)) and median(X) + median(X - median(X))

Let

• nm be the first node with value greater than median(X)
• nl the first node with value not less than median(X) - median(abs(X - median(X)))
• nr the last node with value not greater than median(X) + median(abs(X + median(X)))

For simplicity let's consider only the odd sized case.

To find nl and nr we use the fact that half of the nodes are between nl and nr and that there is no node n with abs(median(X) - n.value) between median(X) - nl.value and nr.value - median(X).

def find_mid_section(root, medianX, N):
while ri - li != N / 2
if ri - li < N / 2:
nl, nr, li, ri = expand_mid_section(nl, nr, li, ri, medianX, N)
else ri - li > N / 2:
nl, nr, li, ri = contract_mid_section(nl, nr, li, ri, medianX, N)

def expand_mid_section(nl, nr, li, ri, medianX, N):
'''
Move left to the left node or right node to the right
whichever gives the least expansion
'''
if nl.left and nr.right:
if medianX - nl.left.value < nr.right.value - medianX:
nl, li = move_left(nl, li)
else:
nr, ri = move_right(nr, ri)
if nl.left:
nl, li = move_left(nl, li)
else:
nr, ri = move_right(nr, ri)
return nl, nr, li, ri

def contract_mid_section(nl, nr, li, ri, medianX, N):
'''
Move left node to the right node or right to the left
whichever gives the least contraction
'''
if nl.right and nr.left:
if medianX - nl.right.value > nr.left.value - medianX:
nl, li = move_right(nl, li)
else:
nr, ri = move_left(nr, ri)
if nl.left:
nl, li = move_left(nl, li)
else:
nr, ri = move_right(nr, ri)
return nl, nr, li, ri

def move_left(n, i):
n = n.left
i -= 1
if n.right:
i -= n.right.size
return n, i

def move_right(n, i):
n = n.right
i += 1
if n.left:
i += n.left.size
return n, i


The median can be obtained by an index search.

The sum can be computed by a value search adding the values of each node, for instance, sum all the values strictly to the left of a value

def get_sum_lt(root, v):
n = root
sum = root.left.sum
while n:
if n.value < v:
sum += n.value
n = n.right
sum += n.left.sum
else:
n = n.left
sum = n.right.sum + n.value


After adding one would call add_to_branch(n, n.value, 1), before removing a value one would call add_to_branch(n, -n.value, -1) to update the ancestors of a node.

def add_to_branch(n, sum, size):
while n.parent:
n = n.parent
n.sum += sum
n.size += size


I am not giving the complete algorithm, but only the key ideas that you can use to compute the moving MAD with complexity O((N-W)*log(W)), where N is the input size and W is the window size.

• I am not quite sure that I am following you. As the first line is missing the Median expression, are you deriving the MAD or AAD? Also, count(h) and count(l) are equal by definition, no? Jan 19 at 9:10
• I am deriving median absolute deviation, $MAD(X)$ above would be written as sum(abs(X - median(X)) in Matlab. I can imagine an implementation where for an odd size sample the median is in one of the heaps, in that case $count(h) - count(l) = \pm 1$
– Bob
Jan 19 at 10:22
• So, for MAD should not it be MAD(X)=median(|X-m|)? Maybe it would be more clear if you will start from the definition and show the missing step. Sorry for the inconvenience, I am probably missing something... Jan 19 at 12:45
• Sorry, I misread the definition of MAD(X), sorry. I posted an update to the answer. Unfortunately I couldn't think of a way to explain the algorithm without going to the underlying data structure, so I wrote a sort of pseudo code.
– Bob
Jan 19 at 15:24
• I have to admit that I am yet not following all the answer. How the sums are related to the MAD? Jan 20 at 3:32

This answer is inspired by the orignial answer by @Bob, which he found not correct but had many of the elements that I will use here.

For simplicity, let us assume that we have the rolling median from an already implemented median filter with length $$w$$. From all the elements within the window, there are $$z = \frac{w-1}{2}$$ values that maintain $$|v - m| >= MAD$$, where $$m$$ is the already found median of the window. Here we emphasize that we use $$v$$ only for the elements that maintain this condition and ignore the rest of the elements.

Let us mark $$h \in \{v|v\ge m\}, l \in \{v|v\le m\}$$. Following we can define two new elements from the window, $$o_1, o_2$$ for which we define the following two group of conditions:

1. conditions on values $$h - m\ge o_1 - m$$ $$m - l\ge o_1 - m$$ $$h - m\ge m - o_2$$ $$m - l\ge m - o_2$$
2. condition on orders $$\text{ord}(o_1) - \text{ord}(o_2) = z$$

Therefore, all we have to do is to create two median heaps with the second condition and use the following actions iteratively for each new value entering the window:

1. Maintain the median heaps.
2. Increase/decrease the orders of both of the heaps to maintain condition 1.
3. The MAD is the maximal deviation produced by one of the two limit order values.

Since the maintenance of the median heaps is of logarithmic order of complexity, this solution logarithmic order of complexity. The code for this solution is attached here.

• I am not sure if I understood correctly. If you add a pseudo-code it could help. If I understood correctly. after taking into account the complexity of balancing the heaps, this algorithm has worst case O(W * log(W)). As you wrote ~absolute deviation of each element in the window can be extremely modified with any modification of the median~
– Bob
Jan 22 at 8:46
• Suppose window size 1001 and input as alternating +1, -1. The median will be alternating +1 and -1. If the heaps are exclude the values equal to the median. Every iteration ~500 elements become smaller/larger than the median. Or you would keep only two counters of ~500
– Bob
Jan 22 at 8:47
• Also, how do you compute the median(abs(X - median(X))) given the state?
– Bob
Jan 22 at 8:50
• I can't tell for sure, maybe your algorithm works as you said, only the description is incomplete. By the way I liked this problem. What are the applications for of this filter?
– Bob
Jan 22 at 8:53
• mathworks.com/help/signal/ref/hampel.html Jan 22 at 11:30