Медианный фильтр на службе разработчика

filter

Недавно пришлось столкнуться с необходимостью программной фильтрации данных АЦП. Гугление и курение (различной документации) привело меня к двум технологиям: Фильтр низких частот (ФНЧ) и Медианный фильтр. Про ФНЧ есть весьма подробная статья в сообществе Easyelectronics, поэтому далее речь пойдёт про медианный фильтр.


Дисклеймер (отмазка): Эта статья по большей частью является практически дословным переводом статьи с сайта embeddedgurus. Однако, переводчик (я) тоже использовал приведенные алгоритмы в работе, нашёл их полезными, и, возможно, представляющими интерес для этого сообщества.

Итак, любой линейный фильтр создан, чтобы пропускать сигналы в заданной полосе частот, и максимально ослаблять все остальные. Подобные фильтры незаменимы, если требуется устранить влияние всевозможных шумов. Однако, в реальном мире встраиваемых систем, разработчик может столкнуться с фактом, что эти классические фильтры практически бесполезны против кратковременных мощных выбросов.

Этот тип шума обычно возникает от какого-либо случайного события, такого, как электростатический разряд, сработавший рядом с прибором брелок сигнализации и прочее. При этом входной сигнал может принять заведомо невозможное значение. Например, с АЦП поступили данные: 385, 389, 388, 388, 912, 388, 387. Очевидно, что значение 912 тут ложное, и должно быть отброшено. При использовании классического фильтра, почти наверняка это большое число повлияет на выходное значение очень сильно. Очевидным решением тут будет применение медианного фильтра.

В соответствии со своим названием, медианный фильтр пропускает среднее значение из множества значений. Обычно размер этой группы нечётный, чтобы избежать двусмысленности при выборе среднего значения. Основная идея-имеется некий буфер с несколькими значениями, из которых выбирается медиана.

Отличия медианы от среднего арифметического
Предположим, что в одной комнате оказалось 19 бедняков и один миллиардер. Каждый кладёт на стол деньги — бедняки из кармана, а миллиардер — из чемодана. По $5 кладёт каждый бедняк, а миллиардер — $1 млрд (109). В сумме получается $1 000 000 095. Если мы разделим деньги равными долями на 20 человек, то получим $50 000 004,75. Это будет среднее арифметическое значение суммы наличных, которая была у всех 20 человек в этой комнате.

Медиана в этом случае будет равна $5 (полусумма десятого и одиннадцатого, срединных значений ранжированного ряда). Можно интерпретировать это следующим образом. Разделив нашу компанию на две равные группы по 10 человек, мы можем утверждать, что в первой группе каждый положил на стол не больше $5, во второй же не меньше $5. В общем случае можно сказать, что медиана это то, сколько принёс с собой средний человек. Наоборот, среднее арифметическое — неподходящая характеристика, так как оно значительно превышает сумму наличных, имеющуюся у среднего человека.
ru.wikipedia.org/wiki/Медиана_(статистика)


По размеру этого множества разделим фильтры на два типа:
•Размерность = 3
•Размерность > 3

Фильтр размерностью 3
Размерность три — наименьшая из возможных. Вычислить среднее значение возможно, использовав лишь несколько операций IF. Ниже приведён код, реализующий этот фильтр:


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;
}


Фильтр размерностью >3
Для фильтра размерностью больше трёх предлагаю воспользоваться алгоритмом, предложенным Филом Экстормом (Phil Ekstrom) в Ноябрьском номере журнала «Embedded Systems», и переписанного с Dynamic C на стандартный С Найджелом Джонсом (Nigel Jones). Алгоритм использует односвязный список, и использует тот факт, что когда массив отсортирован, удаление самого старого значения, и добавление нового не нарушает сортировку.


#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;
}

Чтобы воспользоваться этим кодом, просто вызываем функцию каждый раз, когда появляется новое значение. Она вернёт медианное из последних MEDIAN_FILTER_SIZE измерений.
Этот подход требует довольно много RAM, т.к. приходится хранить и значения, и указатели. Однако он довольно быстрый (58мкс на 40МГц PIC18).

Выводы
Как и большинство других вещей в мире встраиваемых систем, Медианный фильтр имеет свою цену. Например, он вносит задержку в одно чтение при непрерывно растущих входных значениях. К тому же, этот фильтр сильно искажает информацию о частоте сигнала. Конечно, если нам интересна только постоянная составляющая, это не создает особых проблем.

11 комментариев

avatar
А для чего ты его использовал, если не секрет?
avatar
И еще вопросик, можно ли одновременно прицепить оверсэмплинг?
avatar
Есть прибор, в котором надо контролировать выходную мощность. Мощность измеряется детектором, и затем проверяется, нет ли превышения порога. Но ложные срабатывания крайне нежелательны. Для защиты от них в ТЗ было заложено среднее арифметическое 3-х измерений с АЦП, т.е. натурально 3 раза измеряем, суммируем, делим на 3. Но по факту, это очень плохо работало, т.к. большой выброс всё равно сильно влияет на получившееся значение (скажем, нормальное значение было 300, порог 320, измерили 299, 305, 600, среднее — 401, т.е. защита сработает). Плюс деление на 8 бит МК тоже не айс. Медианный фильтр обработает такую ситуацию корректно, и не требует деления.
avatar
Хорошая статья, спасибо.

По поводу реализации фильтра размерностью 3 — по моему так проще и менее ресурсоемко:
#define MAX(a, b) ((a) > (b)) ? (a) : (b)
#define MIN(a, b) ((a) < (b)) ? (a) : (b)
#define MIDDLE(a, b, c) MAX(MIN(a, b), MIN(a, c))

или, если без дефайнов:
min1 = (a < b) ? a : b;
min2 = (a < c) ? a : c;
middle = (min1 > min2) ? min1 : min2

Не, так не работает.
Сейчас еще подумаю :)
avatar
Подумал.
Вот так вроде правильно и все еще менее ресурсоемко:
#define MAX(a, b) ((a) > (b)) ? (a) : (b)
#define MIN(a, b) ((a) < (b)) ? (a) : (b)
#define MIDDLE(a, b, c) MAX(MAX(MIN(a, b), MIN(a, c)), MIN(b, c))
avatar
Использовать медианный фильтр надо очень осторожно. Он относится к нелинейным фильтрам, а Это приводит к нескольким критичным моментам:
1. Происходит появление новых гармоник в спектре выходного сигнала (линейный фильтр такого не делает)
2. По выходному сигналу невозможно восстановить исходный.
Фильтр со скользящим средним=с арифметическим усреднением, тоже не хорош, поскольку его частотная характеристика — периодическая функция, которая при увеличении частоты сигнала не падает в ноль, т.е. некоторые частоты выше частоты отсечки такой фильтр совсем не сгладит.
На мой взгляд, одним из наиболее удобным в реализации цифровым фильтром является экспоненциальный фильтр, преобразование занимает мало-мало места-времени, частотная характеристика имеет спад 6dB/oct, аналогично как и у обычного электрического RC-фильтра.
avatar
Дык расскажи нам про экспоненциальный фильтр в виде статьи с кодом. Будет полезно.
avatar
Ну хорошо, но не раньше выходных
avatar
На сколько я понимаю, с подобной задачей справляется ФНЧ (LPF).
И с точки зрения куда ФНЧ намного компактнее.
avatar
Дак в статье рассказали, чем медианный фильтр отличается от обычных БИХ/КИХ.
avatar
«ФНЧ хорошо подходит для аддитивного белого гауссовского шума, а медианный-для одиночных мощных выбросов» ©Автор оригинальной статьи в комментах. Я прогнал на симуляторе оба этих фильтра с одним и тем же набором входных данных, и получил, что для того, чтобы отфильтровать импульс, равный троекратному значению нормального сигнала, и при этом не выйти за требуемый порог (нормальный уровень+10%), альфа и бэта ФНЧ фильтра должны быть таковы, что время установления Ntau составит несколько секунд, что для моей задачи неприемлемо. В случае же медианного фильтра же время время реакции составляет 1-5 измерений (в зависимости от выбранной размерности).
Только зарегистрированные и авторизованные пользователи могут оставлять комментарии.