NOTE: I have heavily edited this blog post since it was originally published, based on some recent testing
If your engineering education was anything like mine then I’m sure that you learned all about different types of linear filters whose essential objective was to pass signals within a certain frequency band and to reject as far as possible all others. These filters are of course indispensable for many types of ‘noise’. However in the real world of embedded systems it doesn’t take one too long to realize that these classical linear filters are useless against burst noise. This kind of noise typically arises from a quasi-random event. For example a 2-way radio may be keyed next to your product or an ESD event may occur close to your signal. Whenever this happens your input signal may transiently go to a ridiculous value. For example I have often seen A2D readings that look something like this: 385, 389, 388, 388, 912, 388, 387. The 912 value is presumably anomalous and as such should be rejected. If you try and use a classical linear filter then you will almost certainly find that the 912 reading actually ends up having a significant impact on the output. The ‘obvious’ answer in this case is to use a median filter. Despite the supposed obviousness of this, it’s my experience that median filters are used remarkably infrequently in embedded systems. I don’t know why this is, but my guess is that it is a combination of a lack of knowledge of their existence, coupled with difficulty of implementation. Hopefully this post will go some way to rectifying both issues.
As its name suggests, a median filter is one which takes the middle of a group of readings. It’s normal for the group to have an odd number of members such that there is no ambiguity about the middle value. Thus the general idea is that one buffers a certain number of readings and takes the middle reading.
Now Until recently I recognized three classes of median filter, based purely on the size of the filter. They were:
- Filter size of 3 (i.e. the smallest possible).
- Filter size of 5, 7 or 9 (the most common).
- Filter size of 11 or more.
However, I now espouse a simple dichotomy
- Filter size of 3
- Filter size > 3
Filter size of 3
The filter size of three is of course the smallest possible filter. It’s possible to find the middle value simply via a few if statements. The code below is based on an algorithm described here. Clearly this is small and fast code.
uint16_t middle_of_3(uint16_t a, uint16_t b, uint16_t c)
{
uint16_t middle;
if ((a <= b) && (a <= c))
{
middle = (b <= c) ? b : c;
}
else if ((b <= a) && (b <= c))
{
middle = (a <= c) ? a : c;
}
else
{
middle = (a <= b) ? a : b;
}
return middle;
}
Filter size > 3
For filter sizes greater than 3 I suggest you turn to an algorithm described by Phil Ekstrom in the November 2000 edition of Embedded Systems Programming magazine. With the recent hatchet job on embedded.com I can’t find the original article. However there is a copy here. Ekstrom’s approach is to use a linked list. The approach works essentially by observing that once an array is sorted, the act of removing the oldest value and inserting the newest value doesn’t result in the array being significantly unsorted. As a result his approach works well – particularly for large filter sizes.
Be warned that there are some bugs in the originally published code (which Ekstrom corrected). However given the difficulty of finding anything on embedded.com nowadays I have opted to publish my implementation of his code. Be warned that the code below was originally written in Dynamic C and has been ported to standard C for this blog posting. It is believed to work. However it would behoove you to check it thoroughly before use!
#define STOPPER 0 /* Smaller than any datum */
#define MEDIAN_FILTER_SIZE (13)
uint16_t median_filter(uint16_t datum)
{
struct pair
{
struct pair *point; /* Pointers forming list linked in sorted order */
uint16_t value; /* Values to sort */
};
static struct pair buffer[MEDIAN_FILTER_SIZE] = {0}; /* Buffer of nwidth pairs */
static struct pair *datpoint = buffer; /* Pointer into circular buffer of data */
static struct pair small = {NULL, STOPPER}; /* Chain stopper */
static struct pair big = {&small, 0}; /* Pointer to head (largest) of linked list.*/
struct pair *successor; /* Pointer to successor of replaced data item */
struct pair *scan; /* Pointer used to scan down the sorted list */
struct pair *scanold; /* Previous value of scan */
struct pair *median; /* Pointer to median */
uint16_t i;
if (datum == STOPPER)
{
datum = STOPPER + 1; /* No stoppers allowed. */
}
if ( (++datpoint - buffer) >= MEDIAN_FILTER_SIZE)
{
datpoint = buffer; /* Increment and wrap data in pointer.*/
}
datpoint->value = datum; /* Copy in new datum */
successor = datpoint->point; /* Save pointer to old value's successor */
median = &big; /* Median initially to first in chain */
scanold = NULL; /* Scanold initially null. */
scan = &big; /* Points to pointer to first (largest) datum in chain */
/* Handle chain-out of first item in chain as special case */
if (scan->point == datpoint)
{
scan->point = successor;
}
scanold = scan; /* Save this pointer and */
scan = scan->point ; /* step down chain */
/* Loop through the chain, normal loop exit via break. */
for (i = 0 ; i < MEDIAN_FILTER_SIZE; ++i)
{
/* Handle odd-numbered item in chain */
if (scan->point == datpoint)
{
scan->point = successor; /* Chain out the old datum.*/
}
if (scan->value < datum) /* If datum is larger than scanned value,*/
{
datpoint->point = scanold->point; /* Chain it in here. */
scanold->point = datpoint; /* Mark it chained in. */
datum = STOPPER;
};
/* Step median pointer down chain after doing odd-numbered element */
median = median->point; /* Step median pointer. */
if (scan == &small)
{
break; /* Break at end of chain */
}
scanold = scan; /* Save this pointer and */
scan = scan->point; /* step down chain */
/* Handle even-numbered item in chain. */
if (scan->point == datpoint)
{
scan->point = successor;
}
if (scan->value < datum)
{
datpoint->point = scanold->point;
scanold->point = datpoint;
datum = STOPPER;
}
if (scan == &small)
{
break;
}
scanold = scan;
scan = scan->point;
}
return median->value;
}
To use this code, simply call the function every time you have a new input value. It will return the median of the last MEDIAN_FILTER_SIZE readings. This approach can consume a fair amount of RAM as one has to store both the values and the pointers. However if this isn’t a problem for you then it really is a nice algorithm that deserves to be in your tool box as it is dramatically faster than algorithms based upon sorting.
Median filtering based on sorting
In the original version of this article I espoused using a sorting based approach to median filtering when the filter size was 5, 7 or 9. I no longer subscribe to this belief. However for those of you that want to do it, here’s the basic outline:
if (ADC_Buffer_Full)
{
uint_fast16_t adc_copy[MEDIAN_FILTER_SIZE];
uint_fast16_t filtered_cnts;
/* Copy the data */
memcpy(adc_copy, ADC_Counts, sizeof(adc_copy));
/* Sort it */
shell_sort(adc_copy, MEDIAN_FILTER_SIZE);
/* Take the middle value */
filtered_cnts = adc_copy[(MEDIAN_FILTER_SIZE - 1U) / 2U];
/* Convert to engineering units */
...
}
Final Thoughts
Like most things in embedded systems, median filters have certain costs associated with them. Clearly median filters introduce a delay to a step change in value which can be problematic at times. In addition median filters can completely clobber frequency information in the signal. Of course if you are only interested in DC values then this is not a problem. With these caveats I strongly recommend that you consider incorporating median filters in your next embedded design.
