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

[Forger]

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

Код

template <typename Sample, unsigned int WIDTH>
class FilterAverage
{
public:
    FilterAverage(void)
    {
        for (unsigned index = 0; index < WIDTH; ++index) samples[index] = 0;
        totalSum = 0;
        currentSampleIndex = 0;
    }

    Sample evaluate(Sample sample)
    {
        totalSum = totalSum - samples[currentSampleIndex];
        samples[currentSampleIndex] = sample;
        totalSum = totalSum + sample;
        Sample averageSample = totalSum / WIDTH;
        if (++currentSampleIndex >= WIDTH) currentSampleIndex = 0;
        return(averageSample);
    }

    unsigned int getWidth(void) const { return(WIDTH); }

private:
    Sample        samples[WIDTH];
    Sample        totalSum; // TODO: нужен соотв. тип, иначе переполнение при счете !!!!!
    unsigned int    currentSampleIndex;
};

Пользовать проще простого:
объявляем экзэмпляр:

Код

FilterAverage<signed short, TEMPERATURE_FILTER_WIDTH> temperatureFilter;

пользуем:

Код

signed short temperature = temperatureFilter.evaluate(temperature);

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