Медианный фильтр

03/12/2013 - 00:42 Pavel Bobkov

Введение

Если ваше инженерное образование похоже на мое, тогда вы наверняка много знаете о различных типах линейных фильтров, основная задача которых, пропустить сигнал в одном диапазоне частот и задержать сигналы в остальных диапазонах. Эти фильтры, конечно, незаменимы для многих типов шумов. Однако в реальном мире встраиваемых систем требуется немного времени, чтобы понять, что классические линейные фильтры бесполезны против импульсного шума (burst noise, popcorn noise). 

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

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

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

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

До недавнего времени я выделял три класса медианных фильтров, отличающихся количеством используемых значений:

- фильтр использующий 3 значения (наименьший возможный фильтр),
- фильтр использующий 5, 7 или 9 значений (наиболее используемые),
- фильтр использующий 11 или больше значений.

Сейчас я придерживаюсь более простой классификации:

- фильтр использующий 3 значения,
- фильтр использующий больше 3 значений.

Медианный фильтр на 3

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


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

Для фильтра размером больше 3, я предлагаю вам использовать алгоритм, описанный Филом Экстромом в ноябрьском номере журнала Embedded Systems Programming за 2000 год. Экстром использует связный список. Этот подход хорош тем, что когда массив отсортирован, удаление старого значения и добавление нового не вносит в массив существенный беспорядок. Поэтому этот подход хорошо работает с фильтрами больших размеров.

Имейте ввиду, в оригинальном опубликованном коде были некоторые баги, которые Экстром потом исправил. Учитывая, что на embedded.com сейчас сложно что-то найти, я решил опубликовать свою реализацию его кода. Изначально код был написан на Dynamic C, но для этого поста был портирован на стандартный Си. Код предположительно рабочий, но его полная проверка остается на вашей совести.


#define NULL         0
#define STOPPER   0 /* Smaller than any datum */
#define MEDIAN_FILTER_SIZE 5

uint16_t MedianFilter(uint16_t datum)
{

struct pair{
  struct pair *point; /* Pointers forming list linked in sorted order */
  uint16_t value; /* Values to sort */
};

/* Buffer of nwidth pairs */
static struct pair buffer[MEDIAN_FILTER_SIZE] = {0};
/* Pointer into circular buffer of data */
static struct pair *datpoint = buffer; 
/* Chain stopper */
static struct pair small = {NULL, STOPPER};
/* Pointer to head (largest) of linked list.*/
static struct pair big = {&small, 0}; 

/* Pointer to successor of replaced data item */
struct pair *successor; 
/* Pointer used to scan down the sorted list */
struct pair *scan; 
/* Previous value of scan */
struct pair *scanold;
/* Pointer to median */
struct pair *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.

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

Медианная фильтрация на базе сортировки

В старой версии этой статьи для медианных фильтров размером 5, 7 или 9, я поддерживал подход на основе алгоритмов сортировки. Сейчас я изменил свою мнение. Однако, если вы хотите использовать их, я предоставляю вам базовый код:


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 */
   ...

}

Заключение

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

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

Nigel Jones "Median filtering". Вольный перевод ChipEnable.Ru

Comments   

# Saxel 2013-12-03 06:42
В принципе, как вариант. Ни разу до этого медианные фильтры не использовал, все как-то плавающее среднее больше по душе. Надо будет попробовать.
# Neiver 2013-12-03 09:43
На Electronix-е эта тема обсуждалась. С тестами производидельно сти http://electronix.ru/forum/index.php?showtopic=114436&hl=median&st=0
# Peter 2013-12-03 19:47
Так что ли?
#define MIDDLE_3(a,b,c) ((a)>(b)? \
((b)>(c)?(b):(( a)
# Peter 2013-12-03 19:50
Гм. Тэги забыл)
Code:
#define MIDDLE_3(a,b,c) ((a)>(b)? \
((b)>(c)?(b):((a)<(c)?(a):(c))): \
((b)<(c)?(b):((a)<(c)?(c):(a))))
# Pashgan 2013-12-03 21:47
Да.

У вас недостаточно прав для комментирования.