вычисление медианы массива (или произвольной моды), 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! И во многих случаях для реальных измерений с АЦП он будет давать лучшее время!
А физический смысл для рандомного массива применять медианный фильтр не очень понятный!

[Tanya]

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

"Мой" алгоритм еще быстрее. Вот только имеет ли тут смысл говорить о скорости...

А физический смысл для рандомного массива применять медианный фильтр не очень понятный!

Вот по этой самой причине. Тем более, что массив вполне себе не массив, а последовательность с вполне понятными особенностями.

[AlexandrY]

"Мой" алгоритм еще быстрее. Вот только имеет ли тут смысл говорить о скорости...

Вот по этой самой причине. Тем более, что массив вполне себе не массив, а последовательность с вполне понятными особенностями.

Говорить не надо, надо только измерить.
Если имеете исходник, то выкладывайте.

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

Медиану с успехом применял для декодеров частотно модулированных сигналов.
Вот ниже график сигналов с которыми я имею дело при проектировании силовой электроники.
Вверху исходный сигнал. посередине зашумленный, внизу сравнение 3-х фильтров: желтый - скользящее среднее (31 отсчет), голубой - ФНЧ (Fpass=80Hz, Fstop=1000Hz, Astop=60dB ), розовый - медиана (31 отсчет)
Как видно медиана ближе всех повторяет исходный вид сигнала.
[attachment=78573:Median_vs_Average.png]

[Tanya]

Говорить не надо, надо только измерить.
Если имеете исходник, то выкладывайте.

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


Как видно медиана ближе всех повторяет исходный вид сигнала.

Какой исходник? Одно сравнение и один же инкремент или декремент.
Медиана у Вас разве не задерживает сигнал?

[AlexandrY]

Какой исходник? Одно сравнение и один же инкремент или декремент.
Медиана у Вас разве не задерживает сигнал?

Мы наверно о разном.

Еще раз повторил тесты уже с нормальной библиотекой C-и и с сигналом повторяющим характер зашумленного сигнала АЦП как на графике выше.
Добавил итераций для более надежной статистики.


Алгоритм Xenia оказался самым чувствительным к виду данных.
Алгоритм Ekstrom почти не чувствует изменения характера данных, а главное имеет очень маленькую дисперсию. Однако какая сила эти списки! А с виду кажутся такими громоздкими.
Где там Klen? Пора переводить на STL!

Код

INVCPU v1.00 PCB Test Suit.
System Core Clock = 120000000 Hz
Now. System ticks = 750131l, Sysytem time = 0.006251 s
Clock check. T1= 1440087l, T2= 1560069l, Delta(tick)= 119982l, Delta(s)= 0.0010
----------------------------------------------
Algorithm 'Qsort' (standart C library ):
----------------------------------------------
Results CRC = 7EC7
Number of samplles   = 65
Number of iterations = 6000
Measured max. number of cycles = 15054, time(us)=125.45
Measured min. number of cycles = 9369, time(us)=78.08
Measured aver.number of cycles = 12036, time(us)=100.30

----------------------------------------------
Algorithm 'Shell' (http://electronix.ru/forum/index.php?s=&showtopic=114436&view=findpost&p=1181687):
----------------------------------------------
Results CRC = 7EC7
Number of samplles   = 65
Number of iterations = 6000
Measured max. number of cycles = 6370, time(us)=53.08
Measured min. number of cycles = 3640, time(us)=30.33
Measured aver.number of cycles = 5454, time(us)=45.45

----------------------------------------------
Algorithm 'Select' (Numerical Recipes in C, 8.5 Selecting the Mth Largest p.341):
----------------------------------------------
Results CRC = 7EC7
Number of samplles   = 65
Number of iterations = 6000
Measured max. number of cycles = 3129, time(us)=26.07
Measured min. number of cycles = 834, time(us)=6.95
Measured aver.number of cycles = 1621, time(us)=13.51

----------------------------------------------
Algorithm 'Selip' (Numerical Recipes in C, 8.5 Selecting the Mth Largest p.343):
----------------------------------------------
Results CRC = 7EC7
Number of samplles   = 65
Number of iterations = 6000
Measured max. number of cycles = 15714, time(us)=130.95
Measured min. number of cycles = 2224, time(us)=18.53
Measured aver.number of cycles = 14548, time(us)=121.23

----------------------------------------------
Algorithm 'Xenia' (http://electronix.ru/forum/index.php?s=&showtopic=114436&view=findpost&p=1180687):
----------------------------------------------
Results CRC = 7EC7
Number of samplles   = 65
Number of iterations = 6000
Measured max. number of cycles = 21885, time(us)=182.38
Measured min. number of cycles = 930, time(us)=7.75
Measured aver.number of cycles = 10514, time(us)=87.62

----------------------------------------------
Algorithm 'Neiver' (http://electronix.ru/forum/index.php?s=&showtopic=114436&view=findpost&p=1180944):
----------------------------------------------
Results CRC = 7EC7
Number of samplles   = 65
Number of iterations = 6000
Measured max. number of cycles = 8009, time(us)=66.74
Measured min. number of cycles = 1512, time(us)=12.60
Measured aver.number of cycles = 1535, time(us)=12.79

----------------------------------------------
Algorithm 'Ekstrom' (http://embeddedgurus.com/stack-overflow/tag/median-filter/):
----------------------------------------------
Results CRC = 7EC7
Number of samplles   = 65
Number of iterations = 6000
Measured max. number of cycles = 1325, time(us)=11.04
Measured min. number of cycles = 1256, time(us)=10.47
Measured aver.number of cycles = 1265, time(us)=10.54

----------------------------------------------
Algorithm 'Klen' (http://electronix.ru/forum/index.php?s=&showtopic=114436&view=findpost&p=1180679)
----------------------------------------------
Results CRC = 7EC7
Number of samplles   = 65
Number of iterations = 6000
Measured max. number of cycles = 4405, time(us)=36.71
Measured min. number of cycles = 1595, time(us)=13.29
Measured aver.number of cycles = 2403, time(us)=20.02

[demiurg_spb]

Ещё можно попробовать сортировать пузырём, но не весь массив, а лишь половину (ИМХО, это будет самый компактный по объёму кода алгоритм).