вычисление медианы массива (или произвольной моды), 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

[AlexandrY]

Ну что, пришло время подвести некоторые итоги.

Ниже результаты тестирования алгоритмов на платформе на основе Kinetis MK60 (Cortex-M4).
Программа выполнялась из внутренней Flash, данные во внутренней RAM. Частота ядра 120 МГц, частота шины 60 МГц, частота Flash 24 МГц.
Компилятор Keil ARMCC v5.03. Ключи компиляции: -c --cpu Cortex-M4.fp -D__MICROLIB -g -O3 -Otime --apcs=interwork

Я проводил тестирование на интересующем меня объеме выборки. Кому надо дальше может тестировать сам. Проект выложу.
Как видно абсолютный победитель - алгоритм Ekstrom .
Но надо сказать, что в нем есть нюанс. В данных не должно быть числа равного символу STOPPER.
В алгоритме STOPPER равен 0, значит в данных не должно быть 0, что я сделал формирую случайные данные.

Если такая кривизна не устраивает, то лучший алгоритм select из рецептов на C.

Код

INVCPU v1.00 PCB Test Suit.
System Core Clock = 120000000 Hz
Now. System ticks = 750055l, Sysytem time = 0.006250 s
Clock check. T1= 1440076l, T2= 1560048l, Delta(tick)= 119972l, Delta(s)= 0.0010
----------------------------------------------
Algorithm 'Qsort' (standart C library ):
----------------------------------------------
Results CRC = D6EF
Number of samplles   = 65
Number of iterations = 1000
Measured max. number of cycles = 67287, time(us)=560.73
Measured min. number of cycles = 23532, time(us)=196.10
Measured aver.number of cycles = 26937, time(us)=224.48

----------------------------------------------
Algorithm 'Select' (Numerical Recipes in C, 8.5 Selecting the Mth Largest p.341):
----------------------------------------------
Results CRC = D6EF
Number of samplles   = 65
Number of iterations = 1000
Measured max. number of cycles = 3142, time(us)=26.18
Measured min. number of cycles = 1237, time(us)=10.31
Measured aver.number of cycles = 2087, time(us)=17.39

----------------------------------------------
Algorithm 'Selip' (Numerical Recipes in C, 8.5 Selecting the Mth Largest p.343):
----------------------------------------------
Results CRC = D6EF
Number of samplles   = 65
Number of iterations = 1000
Measured max. number of cycles = 18478, time(us)=153.98
Measured min. number of cycles = 2198, time(us)=18.32
Measured aver.number of cycles = 16979, time(us)=141.49

----------------------------------------------
Algorithm 'Xenia' (http://electronix.ru/forum/index.php?s=&showtopic=114436&view=findpost&p=1180687):
----------------------------------------------
Results CRC = D6EF
Number of samplles   = 65
Number of iterations = 1000
Measured max. number of cycles = 40263, time(us)=335.53
Measured min. number of cycles = 920, time(us)=7.67
Measured aver.number of cycles = 19902, time(us)=165.85

----------------------------------------------
Algorithm 'Ekstrom' (http://embeddedgurus.com/stack-overflow/tag/median-filter/):
----------------------------------------------
Results CRC = D6EF
Number of samplles   = 65
Number of iterations = 1000
Measured max. number of cycles = 1243, time(us)=10.36
Measured min. number of cycles = 1182, time(us)=9.85
Measured aver.number of cycles = 1187, time(us)=9.89

----------------------------------------------
Algorithm 'Klen' (http://electronix.ru/forum/index.php?s=&showtopic=114436&view=findpost&p=1180679)
----------------------------------------------
Results CRC = D6EF
Number of samplles   = 65
Number of iterations = 1000
Measured max. number of cycles = 3517, time(us)=29.31
Measured min. number of cycles = 1732, time(us)=14.43
Measured aver.number of cycles = 2470, time(us)=20.58

Да, а Klen-у выкинуть свой компилятор.
Скомпилированное Keil-ом 127 сэмплов по его алгоритму вычисляет за 53 мкс на 120 МГц ядре против 56 мкс на 168 МГц ядре скомпилированное GCC.

Если отключить опцию MICROLIB в компиляторе то алгоритм qsotr ускорится где-то в 9 раз.

[KRS]

Как видно абсолютный победитель - алгоритм Ekstrom .

Это можно было бы понять и без тестов, как я и писал выше он использует то, что при добавление очередного элемента и удалении предыдущего массив уже отсортирован!

Но как видно минимально возможное время все таки у алгоритма Xenia! И во многих случаях для реальных измерений с АЦП он будет давать лучшее время!
А физический смысл для рандомного массива применять медианный фильтр не очень понятный!

[demiurg_spb]

Но как видно минимально возможное время все таки у алгоритма Xenia! И во многих случаях для реальных измерений с АЦП он будет давать лучшее время!

Ещё огромным плюсом алгоритма Xenia является то, что он не портит исходный массив и не требует его копирования, что тоже и время и расход ОЗУ!