реклама на сайте
подробности

 
 
4 страниц V   1 2 3 > »   
Reply to this topicStart new topic
> вычисление медианы массива (или произвольной моды), stm32f4 + STL
klen
сообщение Jul 28 2013, 16:50
Сообщение #1


бессмертным стать можно тремя способами
*****

Группа: Свой
Сообщений: 1 405
Регистрация: 9-05-06
Из: Москва
Пользователь №: 16 912



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

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
Go to the top of the page
 
+Quote Post
Xenia
сообщение Jul 28 2013, 18:33
Сообщение #2


Гуру
******

Группа: Модератор FTP
Сообщений: 4 479
Регистрация: 20-02-08
Из: Москва
Пользователь №: 35 237



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

Осмелюсь предложить алгоритм собственного сочинения. Содержимое массива не трогает, копий с него не снимает. Для МК (когда мало памяти) - самое оно.
Опять же 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. Готова сразиться по скорости с мострообразным алгоритмом Клена. sm.gif
Go to the top of the page
 
+Quote Post
AlexandrY
сообщение Jul 28 2013, 19:01
Сообщение #3


Ally
******

Группа: Модераторы
Сообщений: 6 232
Регистрация: 19-01-05
Пользователь №: 2 050



Цитата(klen @ Jul 28 2013, 19:50) *
результаты по быстродействию при спецификации шаблона типом uint16_t на STM32F405RGT 168MHz :
Код
размер буфера                                    добавления элемента в фифо+ вычисление медианы(максимальное значение)
4                                                                                             4,2мкс
16                                                                                             10мкс
32                                                                                            16мкс
128                                                                                           56мкс

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


Какие-то несоразмерно большие длительности. Как измерялось время?
Go to the top of the page
 
+Quote Post
Forger
сообщение Jul 28 2013, 19:49
Сообщение #4


Профессионал
*****

Группа: Свой
Сообщений: 1 215
Регистрация: 22-02-05
Пользователь №: 2 831



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

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


--------------------
Кругозор некоторых людей - круг с нулевым радиусом. Они называют его "точкой зрения".
Go to the top of the page
 
+Quote Post
AlexandrY
сообщение Jul 28 2013, 20:29
Сообщение #5


Ally
******

Группа: Модераторы
Сообщений: 6 232
Регистрация: 19-01-05
Пользователь №: 2 050



Цитата(Forger @ Jul 28 2013, 22:49) *
Я для таких целей использую фильтр - скользящее среднее, ресурсов жрет минимум, считает оч. быстро.


А "оч. быстро" в микросекундах можно с указанием для какого ARM-а?

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


Go to the top of the page
 
+Quote Post
klen
сообщение Jul 28 2013, 21:12
Сообщение #6


бессмертным стать можно тремя способами
*****

Группа: Свой
Сообщений: 1 405
Регистрация: 9-05-06
Из: Москва
Пользователь №: 16 912



Цитата(Xenia @ Jul 28 2013, 22:33) *
Плохой алгоритм - через полную сортировку медиану считает. sm.gif



плохой или не плохой это понятно, но полную сортировку не делает!
http://www.cplusplus.com/reference/algorithm/nth_element/
"Rearranges the elements in the range [first,last), in such a way that the element at the nth position is the element that would be in that position in a sorted sequence.
The other elements are left without any specific order, except that none of the elements preceding nth are greater than it, and none of the elements following it are less.
The elements are compared using operator< for the first version, and comp for the second."

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


sm.gif я согласен, тока я ручками попишу сам sm.gif, STL это был "вброс" sm.gif

время мерил осцилографом.
Go to the top of the page
 
+Quote Post
Xenia
сообщение Jul 28 2013, 21:36
Сообщение #7


Гуру
******

Группа: Модератор FTP
Сообщений: 4 479
Регистрация: 20-02-08
Из: Москва
Пользователь №: 35 237



Цитата(klen @ Jul 29 2013, 01:12) *
sm.gif я согласен, тока я ручками попишу сам sm.gif, STL это был "вброс" sm.gif
время мерил осцилографом.


Было бы хорошо, если бы и мой алгоритм запрограммировали (благо, что он очень короткий) и с вашим сравнили на одном и том же МК, в тех же самых условиях. А то, я и в АРМах не очень шарю, и осциллографа у меня нету sm.gif.
Go to the top of the page
 
+Quote Post
klen
сообщение Jul 28 2013, 22:09
Сообщение #8


бессмертным стать можно тремя способами
*****

Группа: Свой
Сообщений: 1 405
Регистрация: 9-05-06
Из: Москва
Пользователь №: 16 912



Цитата(Xenia @ Jul 29 2013, 01:36) *
Было бы хорошо, если бы и мой алгоритм запрограммировали (благо, что он очень короткий) и с вашим сравнили на одном и том же МК, в тех же самых условиях. А то, я и в АРМах не очень шарю, и осциллографа у меня нету sm.gif.

закодим Ваш вариант померим и запостим, походу дела ктонить еще чтонибудь интесное по теме сообщит.

STL не точтобы тяжелый -его просто к месту нада примнять, например хеш свой городить в мапе.....
Go to the top of the page
 
+Quote Post
demiurg_spb
сообщение Jul 29 2013, 06:58
Сообщение #9


неотягощённый злом
******

Группа: Свой
Сообщений: 2 746
Регистрация: 31-01-08
Из: Санкт-Петербург
Пользователь №: 34 643



А вот мой боевой вариант для uint16_t данных.
Код
static void bubble_sort(uint16_t* p, size_t n)
{
    if (n<2) return;
    
    for (size_t i=0; i<n; i++)
    {
        for (size_t j=1; j<(n-i); j++)
        {
            if (p[j-1]>p[j])
            {
                uint16_t temp = p[j-1];
                p[j-1] = p[j];
                p[j] = temp;
            }
        }
    }
}

uint16_t median_filter(uint16_t* p, size_t n)  // портит исходный массив
{
    if (!n)  return (0);

    bubble_sort(p, n);

    size_t i = n/2;

    if (n&1)  return (p[i]);                      // Число элементов массива нечетное.
    else      return (((uint32_t)p[i-1]+p[i])/2); // Число элементов массива четное.
}

Посмотрел теорию и понял, что сложность алгоритма с использованием сортировки пузырём O(n^2), тогда как сложность Selection algorithm вроде как O(kn).
Есть над чем задуматься. Спасибо Клёну за идею!


--------------------
“Будьте внимательны к своим мыслям - они начало поступков” (Лао-Цзы)
Go to the top of the page
 
+Quote Post
Golikov A.
сообщение Jul 29 2013, 13:29
Сообщение #10


Гуру
******

Группа: Свой
Сообщений: 4 256
Регистрация: 17-02-06
Пользователь №: 14 454



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

x - вход
y - выход
exp (0 .. 1.0] - параметр фильтра

yn = y(n-1) + (xn - y(n-1))*exp


малые изменения сглажены, большие - резкие отфильтрованы.

потом опять же если вы в итоге все равно по всему массиву ходите, чего просто не посчитать дельту между точками и средню дельту по всему массиву, это требует 2 проходов (1 если длина известна) по массиву, зато потом вы сразу знаете все точки что имеют выбросы. Можно среднюю дельту не считать, а задать ее порогом, тогда с 1 прохода выкидываете точки.

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

И уж точно не стал бы библиотеки подтягивать, в целом шаблоны бы тоже не стал. Написал бы если очень хочется 4 функции: инт16-32 флоат, дабл, остальные можно дописать по мере необходимости, благо функции маленькие.
Go to the top of the page
 
+Quote Post
RabidRabbit
сообщение Jul 29 2013, 14:18
Сообщение #11


Местный
***

Группа: Свой
Сообщений: 397
Регистрация: 3-12-09
Из: Россия, Москва
Пользователь №: 54 040



По моему личному представлению всяческие скользящие средние и прочие экспоненциальные фильтры ни разу не избавляют от одиночных ошибочных отсчётов (выбросов), и в таких случаях, по-моему, медианной фильтрации альтернативы нет sm.gif
Go to the top of the page
 
+Quote Post
neiver
сообщение Jul 30 2013, 07:02
Сообщение #12


Местный
***

Группа: Участник
Сообщений: 214
Регистрация: 22-03-10
Из: Саратов
Пользователь №: 56 123



Ура мерение длинной скоростью кода!!! cool.gif
Вот мой вариант:
Код
int Median(const int *data, size_t size)
{
    size_t middle = size / 2;
    int min, max;
    min = max = data[0];
    for(size_t i = 0; i < size; i++)
    {
        if(data[i] > max)
            max = data[i];
        if(data[i] < min)
            min = data[i];
    }

    int medianCandidate;
    size_t less, more;
    do
    {
        medianCandidate = (min + max) / 2;
        more = less = 0;
        for(size_t i = 0; i < size; i++)
        {
            more += data[i] > medianCandidate;
            less += data[i] < medianCandidate;
        }
        if(more < middle)
        {            
            max = medianCandidate;
        }
        else
        {
            min = medianCandidate;
        }
    }
    while(max - min > 1);
    return medianCandidate;
}

Тестировал на ПК, старый Пентиум 4 2.4 ГГц

| N | t, ms |
--------------------------
100000 4
1000000 47
10000000 516
100000000 5250
--------------------------
Как видно асимптотика получилась почти линейная O(k*N), где k - количество бит типе данных значений.

Сообщение отредактировал neiver - Jul 30 2013, 10:17
Go to the top of the page
 
+Quote Post
Altemir
сообщение Jul 30 2013, 11:43
Сообщение #13


Местный
***

Группа: Свой
Сообщений: 249
Регистрация: 2-05-06
Из: Россия, Поволжье
Пользователь №: 16 686



Добавлю свои 5 копеек. Со времён LPC2132 использую сортировку Шелла, после чего беру центральное значение массива. Массивы обрабатываю до 256 элементов. Реализация - обычная на голых сях. Самая наглядная визуализация, показывающая преимущество алгоритма, здесь: http://www.sorting-algorithms.com/ На хабре немного обсуждения по теме: http://habrahabr.ru/post/117200/
Go to the top of the page
 
+Quote Post
neiver
сообщение Jul 30 2013, 13:15
Сообщение #14


Местный
***

Группа: Участник
Сообщений: 214
Регистрация: 22-03-10
Из: Саратов
Пользователь №: 56 123



Цитата(Xenia @ Jul 28 2013, 22:33) *
Плохой алгоритм - через полную сортировку медиану считает. sm.gif

Хороший. Полную сортировку не делает. Делает Heap sort пока все елементы справа данного (середины) не станут больше него, а слева - меньше. Сложность от O(N) до O(n*Log(N)) в зависимости от того насколько данные упорядочены.
Цитата(Xenia @ Jul 28 2013, 22:33) *
До полного перебора всех точек обычно не доходит, т.к. подходящая точка находится раньше (где-то в середине массива), т.е. тоже ДОСРОЧНО.

А полного перебора и не надо для того чтоб O(N2) получить sm.gif
Go to the top of the page
 
+Quote Post
AlexandrY
сообщение Jul 30 2013, 16:03
Сообщение #15


Ally
******

Группа: Модераторы
Сообщений: 6 232
Регистрация: 19-01-05
Пользователь №: 2 050



Цитата(Xenia @ Jul 28 2013, 21:33) *
P.P.S. Готова сразиться по скорости с мострообразным алгоритмом Клена. sm.gif

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

Да, и проверяем только максимальное время, ибо случайные быстрые проходы в DSP обработке мало кого волнуют.
Go to the top of the page
 
+Quote Post

4 страниц V   1 2 3 > » 
Reply to this topicStart new topic
2 чел. читают эту тему (гостей: 2, скрытых пользователей: 0)
Пользователей: 0

 


RSS Текстовая версия Сейчас: 13th July 2025 - 05:02
Рейтинг@Mail.ru


Страница сгенерированна за 0.01546 секунд с 7
ELECTRONIX ©2004-2016