вычисление медианы массива (или произвольной моды), stm32f4 + STL Опции

[klen]

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

1. Итак, была проблема со съемом сигнала с потенциометра "крутилка рукой"- заказчик не захотел энкодеры ;( . в силу независящих от меня причин оказалось на резистеоре определенное количество шума в виде импульсных выбросов. попросили их "вырезать". проц - STM32F405RGT 168MHz
2. После анализа хар-к помехи пришла мысль выдернуть ее медианным фильтром.
3. дока по быстрому вычислению мод масивом прилагается http://klen.org/Files/Dosc/math/median.pdf
4. Как всегда захотелось решить задачу в общем виде. поэтому был пременен код STL в реализации GCC 4.9.0.
был нарисован шаблонный класс который реализует фифо и вычисляет текущую медиану, оформлено в виде заголовка С++:

CODE

/*
* moda.h
*
* Created on: 24 июля 2013 г.
* Author: klen
* URL: klen.org
*/

#ifndef __MODA_H__
#define __MODA_H__


#include <vector>
#include <algorithm>

//---------------------------------------------------------------------------------------------
template <typename T> class TModaFilterContainer
{
private:
protected:
public:

std::vector<T> vals; // входная последовательность
std::vector<T> gist; // копия входной последовательности для сортировки и поиска моды

TModaFilterContainer () {};
TModaFilterContainer (size_t size) { vals.resize(size); };
~TModaFilterContainer() {};
void Add(T val) {vals.push_front(val); }
void Remove () {vals.pop_back(); }
void Clear() { vals.clear(); }
void Resize(size_t size) {vals.resize(size); }
void Update(T val) { vals.pop_back(); vals.insert(vals.begin(), val); }
};

//---------------------------------------------------------------------------
template <typename T> class TMedianFilterSort : public TModaFilterContainer<T>
{
public:
TMedianFilterSort():TModaFilterContainer<T>() {};
TMedianFilterSort(size_t size):TModaFilterContainer<T>(size) {};
~TMedianFilterSort() {};

T Find()
{
this->gist = this->vals;
div_t divresult;
divresult = div( this->gist.size() , 2 );
size_t n = divresult.quot;
sort(this->gist.begin(), this->gist.end());
if ( divresult.rem )
return this->vals[n];
else
return (this->vals[n-1] + this->vals[n]) / static_cast<T>(2);
}

};
//--------------------------------------------------------------------------
template <typename T> class TMedianFilterNth : public TModaFilterContainer<T>
{
public:
TMedianFilterNth():TModaFilterContainer<T>() {};
TMedianFilterNth(size_t size):TModaFilterContainer<T>(size) {};
~TMedianFilterNth() {};

T Find()
{
this->gist = this->vals;
size_t n = this->gist.size() / 2;
nth_element(this->gist.begin(), this->gist.begin() + n, this->gist.end());
return this->gist[n];
}
};

// EXAMPLE
#if 0

#include "median.h"

typedef TMedianFilterNth<int16_t> median_t;
volatile int16_t mediana;


median_t m(8);

mediana = m.Find();

m.Update(8);
mediana = m.Find();
m.Update(9);
mediana = m.Find();
m.Update(10);
mediana = m.Find();
m.Update(8);
mediana = m.Find();
m.Update(10);
mediana = m.Find();
m.Update(11);
mediana = m.Find();
m.Update(0);
mediana = m.Find();
m.Update(1);
mediana = m.Find();
m.Update(2);
mediana = m.Find();
m.Clear();

// результат работы:

0 0 0 0 0 0 0 0
mediana:0
8 0 0 0 0 0 0 0
mediana:0
9 8 0 0 0 0 0 0
mediana:0
10 9 8 0 0 0 0 0
mediana:0
8 10 9 8 0 0 0 0
mediana:8
10 8 10 9 8 0 0 0
mediana:8
11 10 8 10 9 8 0 0
mediana:9
0 11 10 8 10 9 8 0
mediana:9
1 0 11 10 8 10 9 8
mediana:9
2 1 0 11 10 8 10 9
mediana:9

#endif


#endif /* __MODA_H__ */

результаты по быстродействию при спецификации шаблона типом uint16_t на STM32F405RGT 168MHz :

Код

размер буфера                                    добавления элемента в фифо+ вычисление медианы(максимальное значение)
4                                                                                             4,2мкс
16                                                                                             10мкс
32                                                                                            16мкс
128                                                                                           56мкс

размер кода потянутого шаблонами и алгоритмами ( класс vector + алгоритм thh ) ~1100 байт

5. можно вектор заменить на двухсвязанную очередь deque - по идее должно работать быстрее при добавлении елемента. можно ручками чтоб ваще быстро было - но это время разработки!
6. компилял компиллером из своей сборки http://electronix.ru/forum/index.php?showt...0&start=870
7. это должно работать на любом вменяемом С++ компиллере и наличии стандартной библиотеки шаблонов STL

надюсь с пользой запостил.

Сообщение отредактировал IgorKossak - Aug 1 2013, 06:36

[Xenia]

Плохой алгоритм - через полную сортировку медиану считает.

Осмелюсь предложить алгоритм собственного сочинения. Содержимое массива не трогает, копий с него не снимает. Для МК (когда мало памяти) - самое оно.
Опять же STL-библитека ему не нужна, ибо никаких функций он не вызывает и никаких FIFO-буферов не использует.

Код

datatype median( datatype array, int length)  // массив и его длина
{
  int  slit = length/2;
  for( int i=0; i < length; i++)
  { int s1=0, s2=0;
    datatype val = array[i];
    for( int j=0; j < length; j++)
    { if( array[j] < val)
      { if( ++s1 > slit) goto nexti;
      }
      else if( array[j] > val)
      { if( ++s2 > slit) goto nexti;
      }
    }
    return val;
nexti:
  }
  return 0;  // чистая формальность, досюда исполнение никогда не доходит
}

А если хранить переменные s1, s2, slit и val в регистрах, то скорость совсем хороша.

Алгоритм основан на последовательной проверке точек массива на их "медианность".
Проверка состоит в сравнении точки со всеми остальными - подсчетом двух сумм (s1 и s2) - количества точек, которые оказались больше данной, и количества тех, которые оказались меньше нее.
Превышение показания хотя бы одного из счетчиков свыше половины длины массива (slit) бракует данную точку ДОСРОЧНО и переходит к проверке следующей (переход на nexti).
До полного перебора всех точек обычно не доходит, т.к. подходящая точка находится раньше (где-то в середине массива), т.е. тоже ДОСРОЧНО.

P.S. datatype - произвольный тип данных, для которого определены операции > и <. При желании можно хоть через template все это обобщить.
P.P.S. Готова сразиться по скорости с мострообразным алгоритмом Клена.

[AlexandrY]

P.P.S. Готова сразиться по скорости с мострообразным алгоритмом Клена.

 
Вы будете смеяться Xenia, но ваш алгоритм в 4 раза! медленней по сравнению с применением библиотечной функции qsort на unsigned short.
А на unsigned int в 4.5 раза.
Алгоритм klen пока самый быстрый, если ему верить.
Хотя в Рецептах на C кажется нашел алгоритм в два раза быстрее, но пока не отладил.

Да, и проверяем только максимальное время, ибо случайные быстрые проходы в DSP обработке мало кого волнуют.

[Xenia]

Вы будете смеяться Xenia, но ваш алгоритм в 4 раза! медленней по сравнению с применением библиотечной функции qsort на unsigned short.
А на unsigned int в 4.5 раза.

А какой длины массив использовался при проверке?
Ведь qsort там не пузырьковый, у него логарифмическая скорость. А следовательно на больших массивах эта фора обязательно скажется. Однако в сфере МК память ОЗУ обычно невелика, а потому невелики и массивы. Обычно редко кто ищет медиану в массивах, более длинных чем 1024 элемента. А при такой длине поиск медианы на основе функции qsort врядли даст заметное преимущество.

[demiurg_spb]

Обычно редко кто ищет медиану в массивах, более длинных чем 1024 элемента.

Тут я поддержу Ксению, т.к. в моих задачах (не связанных с применением DSP) средняя глубина медианной фильтрации варьируется от 5 до 31 отсчёта.
Давайте определимся с правилами игры. Предлагаю предоставлять результаты своих замеров для следующих глубин фильтра: 10, 100, 1K и 10K

[Tanya]

варьируется от 5 до 31 отсчёта.
Давайте определимся с правилами игры.

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

[demiurg_spb]

Кроме того сортировка проводится не случайным образом перемешанной последовательности, а априорно известной. Таким образом, возможно, что самым лучшим для данного случая будет не самый быстрый для общего случая алгоритм.

Не понял эту глубокую мысль...
У нас каждый раз новая неизвестная последовательность (не важно все 10 отсчётов новых или только один из десяти (скользящее окно)).
Т.е. всегда что-то новое и неизвестное. А то что мы отсортировали на прошлом шаге для нас не представляет никакой ценности на текущем шаге, т.к. наши отсчёты после сортировки теряют привязку к своей очерёдности (то бишь ко времени).

[Tanya]

Не понял эту глубокую мысль...

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

[demiurg_spb]

Не всегда. У меня есть разработка - тахометр, который используют в весьма разболтанных системах (люфты страшные) + метки на колесе неоднородно размещены бывают, так для нормального функционирования приходится и медианный фильтр и фнч 2 порядка ставить. И никакие конденсаторы тут не спасают.

[Tanya]

. И никакие конденсаторы тут не спасают.

Верю. Но я же написала "возможно"... Медианный фильтр нелинейный и удаляет "неправильные" выбросы.
А откуда они берутся? Разрывается контакт (повисает в воздухе), а входное сопротивление АЦП не очень...
Есть еще вариант - сделать ограничение скорости нарастания (тоже нелинейная штука) - при рассогласовании запомненного значения с АЦПшным оно медленно увеличивается или уменьшается (на фиксированную величину). Это будет самый быстрый способ. Я так думаю.
Недостаток есть - если шаловливые ручонки непрерывно крутят туда-сюда ручки, то будет заниженное значение, если вход АЦП немного подтянут к земле.

[Axel]

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

Иногда (часто) не в ручонках дело. Типичный пример - потенциометры на педали газа автомобилей. Пример из программы ECU привести не могу - не сохранился, но когда однажды пришлось разбираться, то обнаружился вполне серьезный фильтр (насколько удалось понять после дизассемблера).