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

 
 
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
Xenia
сообщение Jul 30 2013, 16:35
Сообщение #16


Гуру
******

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



Цитата(AlexandrY @ Jul 30 2013, 20:03) *
Вы будете смеяться Xenia, но ваш алгоритм в 4 раза! медленней по сравнению с применением библиотечной функции qsort на unsigned short.
А на unsigned int в 4.5 раза.


А какой длины массив использовался при проверке?
Ведь qsort там не пузырьковый, у него логарифмическая скорость. А следовательно на больших массивах эта фора обязательно скажется. Однако в сфере МК память ОЗУ обычно невелика, а потому невелики и массивы. Обычно редко кто ищет медиану в массивах, более длинных чем 1024 элемента. А при такой длине поиск медианы на основе функции qsort врядли даст заметное преимущество.
Go to the top of the page
 
+Quote Post
Axel
сообщение Jul 31 2013, 05:02
Сообщение #17


Местный
***

Группа: Свой
Сообщений: 480
Регистрация: 21-11-04
Пользователь №: 1 188



Мне тоже стало интересно... Попробовал фильтр Ekstrom'а для обработки сигналов touch screen (5 - 15 отсчетов). Понравилось (время не измерял).
Go to the top of the page
 
+Quote Post
demiurg_spb
сообщение Jul 31 2013, 07:25
Сообщение #18


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

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



Цитата(Xenia @ Jul 30 2013, 20:35) *
Обычно редко кто ищет медиану в массивах, более длинных чем 1024 элемента.

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


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


Гуру
******

Группа: Модераторы
Сообщений: 8 752
Регистрация: 6-01-06
Пользователь №: 12 883



Цитата(demiurg_spb @ Jul 31 2013, 11:25) *
варьируется от 5 до 31 отсчёта.
Давайте определимся с правилами игры.

Вот тоже скажу... Нет никакого смысла в оцифровке скрипящего потенциометра последовательностью в 100 точек. Не более 10, если не слишком часто...
Кроме того сортировка проводится не случайным образом перемешанной последовательности, а априорно известной. Таким образом, возможно, что самым лучшим для данного случая будет не самый быстрый для общего случая алгоритм.
Go to the top of the page
 
+Quote Post
demiurg_spb
сообщение Jul 31 2013, 10:07
Сообщение #20


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

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



Цитата(Tanya @ Jul 31 2013, 12:41) *
Кроме того сортировка проводится не случайным образом перемешанной последовательности, а априорно известной. Таким образом, возможно, что самым лучшим для данного случая будет не самый быстрый для общего случая алгоритм.
Не понял эту глубокую мысль...
У нас каждый раз новая неизвестная последовательность (не важно все 10 отсчётов новых или только один из десяти (скользящее окно)).
Т.е. всегда что-то новое и неизвестное. А то что мы отсортировали на прошлом шаге для нас не представляет никакой ценности на текущем шаге, т.к. наши отсчёты после сортировки теряют привязку к своей очерёдности (то бишь ко времени).


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


Гуру
******

Группа: Модераторы
Сообщений: 8 752
Регистрация: 6-01-06
Пользователь №: 12 883



Цитата(demiurg_spb @ Jul 31 2013, 14:07) *
Не понял эту глубокую мысль...

Поясню. Последовательность не произвольная, а реальная - две полочки, а между ними выбросы. Или одна полочка без выбросов. Возможно, что со всем этим можно легко справиться с помощью одного резистора и одного конденсатора.
Go to the top of the page
 
+Quote Post
demiurg_spb
сообщение Jul 31 2013, 11:06
Сообщение #22


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

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



Цитата(Tanya @ Jul 31 2013, 14:59) *
Не всегда. У меня есть разработка - тахометр, который используют в весьма разболтанных системах (люфты страшные) + метки на колесе неоднородно размещены бывают, так для нормального функционирования приходится и медианный фильтр и фнч 2 порядка ставить. И никакие конденсаторы тут не спасают.


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


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

Группа: Модераторы
Сообщений: 1 951
Регистрация: 27-08-04
Из: Санкт-Петербург
Пользователь №: 555



А ведь если фильтр скользящий, на этапе добавление очередного измерения, уже есть отсортированный массив.
Но из него удаляется один элемент, при этом он остается отсортированным, надо просто вставить новый на нужное место. Ту сортировка вставками будет самой быстрой.
Go to the top of the page
 
+Quote Post
Xenia
сообщение Jul 31 2013, 11:50
Сообщение #24


Гуру
******

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



Цитата(demiurg_spb @ Jul 31 2013, 11:25) *
Давайте определимся с правилами игры. Предлагаю предоставлять результаты своих замеров для следующих глубин фильтра: 10, 100, 1K и 10K


Да-да! И тогда я победю Клёна! sm.gif
Go to the top of the page
 
+Quote Post
Tanya
сообщение Jul 31 2013, 12:06
Сообщение #25


Гуру
******

Группа: Модераторы
Сообщений: 8 752
Регистрация: 6-01-06
Пользователь №: 12 883



Цитата(demiurg_spb @ Jul 31 2013, 15:06) *
. И никакие конденсаторы тут не спасают.

Верю. Но я же написала "возможно"... Медианный фильтр нелинейный и удаляет "неправильные" выбросы.
А откуда они берутся? Разрывается контакт (повисает в воздухе), а входное сопротивление АЦП не очень...
Есть еще вариант - сделать ограничение скорости нарастания (тоже нелинейная штука) - при рассогласовании запомненного значения с АЦПшным оно медленно увеличивается или уменьшается (на фиксированную величину). Это будет самый быстрый способ. Я так думаю.
Недостаток есть - если шаловливые ручонки непрерывно крутят туда-сюда ручки, то будет заниженное значение, если вход АЦП немного подтянут к земле.
Go to the top of the page
 
+Quote Post
Axel
сообщение Jul 31 2013, 12:30
Сообщение #26


Местный
***

Группа: Свой
Сообщений: 480
Регистрация: 21-11-04
Пользователь №: 1 188



Цитата(Tanya @ Jul 31 2013, 15:06) *
Недостаток есть - если шаловливые ручонки непрерывно крутят туда-сюда ручки, то будет заниженное значение, если вход АЦП немного подтянут к земле.

Иногда (часто) не в ручонках дело. Типичный пример - потенциометры на педали газа автомобилей. Пример из программы ECU привести не могу - не сохранился, но когда однажды пришлось разбираться, то обнаружился вполне серьезный фильтр (насколько удалось понять после дизассемблера).
Go to the top of the page
 
+Quote Post
AlexandrY
сообщение Jul 31 2013, 13:14
Сообщение #27


Ally
******

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



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

Ниже результаты тестирования алгоритмов на платформе на основе 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-у выкинуть свой компилятор. wink.gif
Скомпилированное Keil-ом 127 сэмплов по его алгоритму вычисляет за 53 мкс на 120 МГц ядре против 56 мкс на 168 МГц ядре скомпилированное GCC.

Если отключить опцию MICROLIB в компиляторе то алгоритм qsotr ускорится где-то в 9 раз.
Go to the top of the page
 
+Quote Post
Altemir
сообщение Jul 31 2013, 14:03
Сообщение #28


Местный
***

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



AlexandrY
"Number of samplles" - это длина массива? Какой тип чисел был во входном массиве? signed long не пробовали с длиной 128 элементов?
Go to the top of the page
 
+Quote Post
AlexandrY
сообщение Jul 31 2013, 14:10
Сообщение #29


Ally
******

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



Цитата(Altemir @ Jul 31 2013, 17:03) *
AlexandrY
"Number of samplles" - это длина массива? Какой тип чисел был во входном массиве? signed long не пробовали с длиной 128 элементов?

Был использован тип unsigned int
Number of samplles - это длинна массива. Выбрана нечетной чтобы был элемент посередине.
Измерялось время между подачей очередного случайного сэмпла (от 1 до 10000)в массив и получением медианы как отклика.
Всего пропускалось по 1000 сэмплов.
Вначале массивы инициализировались единицами (из-за особенности алгоритма Ekstrom )
Go to the top of the page
 
+Quote Post
demiurg_spb
сообщение Jul 31 2013, 14:10
Сообщение #30


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

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



Цитата(KRS @ Jul 31 2013, 15:22) *
А ведь если фильтр скользящий, на этапе добавление очередного измерения, уже есть отсортированный массив.
Но из него удаляется один элемент, при этом он остается отсортированным, надо просто вставить новый на нужное место. Ту сортировка вставками будет самой быстрой.
Это понятно, но вы посмотрите на размер кода такого алгоритма 'Ekstrom'.
Я это тоже изобретал, но отверг по причинам невпихуемости результирующей прошивки.


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


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

Группа: Модераторы
Сообщений: 1 951
Регистрация: 27-08-04
Из: Санкт-Петербург
Пользователь №: 555



Цитата(AlexandrY @ Jul 31 2013, 17:14) *
Как видно абсолютный победитель - алгоритм Ekstrom .

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

Но как видно минимально возможное время все таки у алгоритма Xenia! И во многих случаях для реальных измерений с АЦП он будет давать лучшее время!
А физический смысл для рандомного массива применять медианный фильтр не очень понятный!
Go to the top of the page
 
+Quote Post
demiurg_spb
сообщение Jul 31 2013, 14:28
Сообщение #32


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

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



Цитата(KRS @ Jul 31 2013, 18:19) *
Но как видно минимально возможное время все таки у алгоритма Xenia! И во многих случаях для реальных измерений с АЦП он будет давать лучшее время!
Ещё огромным плюсом алгоритма Xenia является то, что он не портит исходный массив и не требует его копирования, что тоже и время и расход ОЗУ!


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


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

Группа: Модераторы
Сообщений: 1 951
Регистрация: 27-08-04
Из: Санкт-Петербург
Пользователь №: 555



Для определенных измерений его еще можно улучшить, например изначально беря за индекс медианы предыдущий индекс и двигаться по кольцу...
Go to the top of the page
 
+Quote Post
Tanya
сообщение Jul 31 2013, 16:04
Сообщение #34


Гуру
******

Группа: Модераторы
Сообщений: 8 752
Регистрация: 6-01-06
Пользователь №: 12 883



Цитата(KRS @ Jul 31 2013, 18:19) *
Но как видно минимально возможное время все таки у алгоритма Xenia! И во многих случаях для реальных измерений с АЦП он будет давать лучшее время!

"Мой" алгоритм еще быстрее. Вот только имеет ли тут смысл говорить о скорости...
Цитата(KRS @ Jul 31 2013, 18:19) *
А физический смысл для рандомного массива применять медианный фильтр не очень понятный!

Вот по этой самой причине. Тем более, что массив вполне себе не массив, а последовательность с вполне понятными особенностями.
Go to the top of the page
 
+Quote Post
AlexandrY
сообщение Aug 1 2013, 07:45
Сообщение #35


Ally
******

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



Цитата(Tanya @ Jul 31 2013, 19:04) *
"Мой" алгоритм еще быстрее. Вот только имеет ли тут смысл говорить о скорости...

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


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

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

Медиану с успехом применял для декодеров частотно модулированных сигналов.
Вот ниже график сигналов с которыми я имею дело при проектировании силовой электроники.
Вверху исходный сигнал. посередине зашумленный, внизу сравнение 3-х фильтров: желтый - скользящее среднее (31 отсчет), голубой - ФНЧ (Fpass=80Hz, Fstop=1000Hz, Astop=60dB ), розовый - медиана (31 отсчет)
Как видно медиана ближе всех повторяет исходный вид сигнала.
[attachment=78573:Median_vs_Average.png]
Go to the top of the page
 
+Quote Post
Tanya
сообщение Aug 1 2013, 08:13
Сообщение #36


Гуру
******

Группа: Модераторы
Сообщений: 8 752
Регистрация: 6-01-06
Пользователь №: 12 883



Цитата(AlexandrY @ Aug 1 2013, 11:45) *
Говорить не надо, надо только измерить. wink.gif
Если имеете исходник, то выкладывайте.

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


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

Какой исходник? Одно сравнение и один же инкремент или декремент.
Медиана у Вас разве не задерживает сигнал?
Go to the top of the page
 
+Quote Post
AlexandrY
сообщение Aug 1 2013, 09:09
Сообщение #37


Ally
******

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



Цитата(Tanya @ Aug 1 2013, 11:13) *
Какой исходник? Одно сравнение и один же инкремент или декремент.
Медиана у Вас разве не задерживает сигнал?


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

Еще раз повторил тесты уже с нормальной библиотекой 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
Go to the top of the page
 
+Quote Post
demiurg_spb
сообщение Aug 1 2013, 09:17
Сообщение #38


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

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



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


--------------------
“Будьте внимательны к своим мыслям - они начало поступков” (Лао-Цзы)
Go to the top of the page
 
+Quote Post
Tanya
сообщение Aug 1 2013, 09:44
Сообщение #39


Гуру
******

Группа: Модераторы
Сообщений: 8 752
Регистрация: 6-01-06
Пользователь №: 12 883



Цитата(AlexandrY @ Aug 1 2013, 13:09) *
Мы наверно о разном. laughing.gif

Может...
Я про то, что
1. Человек не может крутить ручку с бесконечно большой скоростью.
2. Скрип возникает только в процессе верчения.
3. Выбросы поэтому будут только на фронтах и всегда в одну сторону (к нулю?).
4. Не имеет смысла часто и много все это измерять.
5. Выбросы нужно просто выбрасывать.
Ваш медианный алгоритм может давать очень большую ошибку, если более половины выборки придется на фронт с большими выбросами.
Go to the top of the page
 
+Quote Post
neiver
сообщение Aug 1 2013, 10:20
Сообщение #40


Местный
***

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



Цитата(AlexandrY @ Aug 1 2013, 13:09) *
Алгоритм Neiver-а не прошел верификацию, там получается явно не медиана.

Кстати да, там у меня ошибочка была, при не четном колчестве элементов. Сейчас поправил и чуть оптимизировал.
Код
int Median(const int *data, size_t size)
{
    int min, max, newMax, newMin;
    size_t middle = (size+1) / 2;
    min = max = data[0];
    for(size_t i = 0; i < size; i++)
    {
        max = std::max(data[i], max);
        min = std::min(data[i], min);
    }

    int medianCandidate;
    size_t less, more;
    do
    {
        medianCandidate = (min + max) / 2;
        more = less = 0;
        newMax = min;
        newMin = max;

        for(size_t i = 0; i < size; i++)
        {
            if(data[i] > medianCandidate)
            {
                more ++;
                if (data[i] < newMin)
                    newMin = data[i];
            }else if(data[i] < medianCandidate)
            {
                less ++;
                if (data[i] > newMax)
                    newMax = data[i];
            }
        }

        if (less <= middle && more <= middle)
            break;
        else if(less > middle)
        {            
            max = newMax;
        }
        else
        {
            min = newMin;
        }
    }
    while(max - min > 1);
    size_t equal = size - more - less;
    if (less >= middle) return newMax;
    else if (less + equal >= middle)
        return medianCandidate;
    else return newMin;
}
Go to the top of the page
 
+Quote Post
AlexandrY
сообщение Aug 1 2013, 10:23
Сообщение #41


Ally
******

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



Цитата(Tanya @ Aug 1 2013, 12:44) *
Может...
Я про то, что
1. Человек не может крутить ручку с бесконечно большой скоростью.
2. Скрип возникает только в процессе верчения.
3. Выбросы поэтому будут только на фронтах и всегда в одну сторону (к нулю?).
4. Не имеет смысла часто и много все это измерять.
5. Выбросы нужно просто выбрасывать.
Ваш медианный алгоритм может давать очень большую ошибку, если более половины выборки придется на фронт с большими выбросами.


А..., все про тот же резистор.
Честно говоря уже много лет на моих платах нет переменных резисторов. Ничего по существу в этой теме сказать не могу. sad.gif
Разве что предложил бы поставить сенсорный скроллер. biggrin.gif

Любой фильтр даст ошибку если его применить вне его полосы.
Кстати, знаете что у медианного фильтра можно измерить частотную характеристику?
Go to the top of the page
 
+Quote Post
Altemir
сообщение Aug 1 2013, 10:36
Сообщение #42


Местный
***

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



AlexandrY
Могу скромно попросить на той же платформе проверить медианку на основе сортировки Шелла? rolleyes.gif Заранее благодарю.

CODE

#define compGT(a, B ) (a > B )

//------------------------------------------------------------------------------
//! Медианный фильтр с сортировкой Шелла
//------------------------------------------------------------------------------
int Med_Shell_Sort(int *a, int ub){
int n, h, i, j;
int t;

// compute largest increment
n = ub + 1;
h = 1;
if (n < 14) h = 1;
else
if (sizeof(tblIndex) == 2 && n > 29524) h = 3280;
else
{
while (h < n) h = 3*h + 1;
h /= 3;
h /= 3;
}

while (h > 0)
{
// sort-by-insertion in increments of h
for (i = h; i < ub; i++)
{
t = a[i];
for (j = i-h; (j >= 0) && compGT(a[j], t); j -= h) a[j+h] = a[j];
a[j+h] = t;
}
// compute next increment
h /= 3;
}

return ((ub+1)/2);
}
Go to the top of the page
 
+Quote Post
AlexandrY
сообщение Aug 1 2013, 10:54
Сообщение #43


Ally
******

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



Цитата(neiver @ Aug 1 2013, 13:20) *
Кстати да, там у меня ошибочка была, при не четном колчестве элементов. Сейчас поправил и чуть оптимизировал.


Да теперь правильно.
Обновил в таблице результатов в предыдущем посте.

Ого! Уже и алгоритм Neiver обошел (по некоторым параметрам) алгоритмы из Рецептов .
Значит рецепты пора выбрасывать в утиль. wink.gif

Цитата(Altemir @ Aug 1 2013, 13:36) *
Могу скромно попросить на той же платформе проверить медианку на основе сортировки Шелла? rolleyes.gif Заранее благодарю.


Да, мог бы если скажете где взять sizeof(tblIndex) и что там за магические числа типа 29524 и 3280
Go to the top of the page
 
+Quote Post
Altemir
сообщение Aug 1 2013, 11:09
Сообщение #44


Местный
***

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



Цитата(AlexandrY @ Aug 1 2013, 14:54) *
Да, мог бы если скажете где взять sizeof(tblIndex) и что там за магические числа типа 29524 и 3280


Код
typedef int tblIndex;


Упс. Прошу прощения. sizeof(tblIndex) - это я не всё заменил, чтобы выложить сюда. В моём примере следует использовать 4 - размерность элемента входного массива, байт. 29524 и 3280 - это действительно, магические числа. Брал пример очень давно из какого-то учебника, не вспомнить. Если правильно понимаю, сии числа ограничивают элементарный блок сортировки и шаг при работе с очень большими массивами, чтобы оптимизировать эффективность. Вроде, в книжке был целый ряд магических чисел для ещё бОльших длин, но я не использую у себя массивы, превышающие 1000 элементов, так что указанная часть кода неактуальна. Если найду источник этого примера, скину ссылку.

Ага. Нашёл: http://www.info-system.ru/library/algo/sortsearch.pdf
И отдельные файлы реализации, откуда выдирал:
http://www.cs.auckland.ac.nz/~jmor159/PLDS...emann/s_shl.txt
http://epaperpress.com/sortsearch/txt/shl.txt
Видно, что h может быть получен разными способами.
Go to the top of the page
 
+Quote Post
AlexandrY
сообщение Aug 1 2013, 11:22
Сообщение #45


Ally
******

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



Цитата(Altemir @ Aug 1 2013, 14:09) *
Упс.


Добавил. См. пост №37
Go to the top of the page
 
+Quote Post
Altemir
сообщение Aug 1 2013, 11:24
Сообщение #46


Местный
***

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



AlexandrY
Огромное спасибо! Как и ожидал - у меня стабильный середнячок sm.gif Буду теперь смотреть альтернативу. Благо, они есть в этой ветке
Go to the top of the page
 
+Quote Post
VAI
сообщение Aug 1 2013, 18:45
Сообщение #47


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

Группа: Модераторы
Сообщений: 1 120
Регистрация: 17-06-04
Пользователь №: 37



Я опоздал, все результаты посчитаны, да и алгоритм медианного фильтра у меня поддерживает размер буфера от 3 до 31. 65 никак.
Когда-то выкладывал на сахаре и здесь, но здесь не нашел, вот ссылки на сахару.
http://caxapa.ru/15518.html
http://caxapa.ru/170626.html
В обоих вариантах алгоритм получения медианы одинаков, только способ работы с буфером разный.
AlexandrY Вы обещали выложить проект, хочу скомпилить с буфером размером 31 и сравнить скорость моего алгоритма с другими.
Или Вы добавьте мой алгоритм и выложите результат. :-)


--------------------
Если зайца бить, его можно и спички научить зажигать
Сколько дурака не бей - умнее не будет. Зато опытнее
Go to the top of the page
 
+Quote Post
AlexandrY
сообщение Aug 5 2013, 08:18
Сообщение #48


Ally
******

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



Цитата(VAI @ Aug 1 2013, 21:45) *
... и выложите результат. :-)


Да, вот результат.
Уменьшил массив до 31 отсчета, под алгоритм VAI. Количество итераций довел до 10000.
И сделал сравнение для двух компиляторов Keil (ARM compiler) и IAR.
Чтобы еще было больше доверия к тесту и виднее различие компиляторов вначале идет тест CoreMark.
CODE
INVCPU v1.00 PCB Test Suit. INVCPU v1.00 PCB Test Suit.
System Core Clock = 120000000 Hz System Core Clock = 120000000 Hz
Now. System ticks = 750122l, Sysytem time = 0.006251 s Now. System ticks = 707843l, Sysytem time = 0.005899 s
---------------------------------------------- ----------------------------------------------
Coremark test Coremark test
---------------------------------------------- ----------------------------------------------
2K performance run parameters for coremark. 2K performance run parameters for coremark
CoreMark Size : 666 CoreMark Size : 666
Total ticks : 2090014320 Total ticks : 1619555180
Total time (secs): 17.416786 Total time (secs): 13.496293
Iterations/Sec : 287.079373 Iterations/Sec : 370.472095
Iterations : 5000 Iterations : 5000
Compiler version : ARM CC 5030069 Compiler version : IAR CC 6050006
Compiler flags : -o3 Compiler flags : -o3
Memory location : STACK Memory location : STACK
seedcrc : 0xe9f5 seedcrc : 0xe9f5
[0]crclist : 0xe714 [0]crclist : 0xe714
[0]crcmatrix : 0x1fd7 [0]crcmatrix : 0x1fd7
[0]crcstate : 0x8e3a [0]crcstate : 0x8e3a
[0]crcfinal : 0xbd59 [0]crcfinal : 0xbd59
Correct operation validated. Correct operation validated.
CoreMark 1.0 : 287.079373 / ARM CC 5030069 -o3 / STACK CoreMark 1.0 : 370.472095 / IAR CC 6050006 -o3 / STACK
---------------------------------------------- ----------------------------------------------
Algorithm 'Qsort' (standart C library ): Algorithm 'Qsort' (standart C library ):
---------------------------------------------- ----------------------------------------------
Results CRC = 81AD Results CRC = 7CBB
Number of samplles = 31 Number of samplles = 31
Number of iterations = 10000 Number of iterations = 10000
Measured max. number of cycles = 6452, time(us)=53.77 Measured max. number of cycles = 12865, time(us)=107.21
Measured min. number of cycles = 3722, time(us)=31.02 Measured min. number of cycles = 6340, time(us)=52.83
Measured aver.number of cycles = 4788, time(us)=39.90 Measured aver.number of cycles = 8508, time(us)=70.90

---------------------------------------------- ----------------------------------------------
Algorithm 'Shell' Algorithm 'Shell' (http://electronix.ru/forum/index.php?s=&showtopic=114436&view=findpost&p=1181687):
---------------------------------------------- ----------------------------------------------
Results CRC = 81AD Results CRC = 7CBB
Number of samplles = 31 Number of samplles = 31
Number of iterations = 10000 Number of iterations = 10000
Measured max. number of cycles = 2372, time(us)=19.77 Measured max. number of cycles = 2610, time(us)=21.75
Measured min. number of cycles = 1355, time(us)=11.29 Measured min. number of cycles = 1565, time(us)=13.04
Measured aver.number of cycles = 1949, time(us)=16.24 Measured aver.number of cycles = 2098, time(us)=17.48

---------------------------------------------- ----------------------------------------------
Algorithm 'Select' Algorithm 'Select' (Numerical Recipes in C, 8.5 Selecting the Mth Largest p.341):
---------------------------------------------- ----------------------------------------------
Results CRC = 81AD Results CRC = 7CBB
Number of samplles = 31 Number of samplles = 31
Number of iterations = 10000 Number of iterations = 10000
Measured max. number of cycles = 1453, time(us)=12.11 Measured max. number of cycles = 1530, time(us)=12.75
Measured min. number of cycles = 472, time(us)=3.93 Measured min. number of cycles = 495, time(us)=4.13
Measured aver.number of cycles = 825, time(us)=6.88 Measured aver.number of cycles = 858, time(us)=7.15

---------------------------------------------- ----------------------------------------------
Algorithm 'Selip' Algorithm 'Selip' (Numerical Recipes in C, 8.5 Selecting the Mth Largest p.343):
---------------------------------------------- ----------------------------------------------
Results CRC = 81AD Results CRC = 7CBB
Number of samplles = 31 Number of samplles = 31
Number of iterations = 10000 Number of iterations = 10000
Measured max. number of cycles = 3325, time(us)=27.71 Measured max. number of cycles = 3590, time(us)=29.92
Measured min. number of cycles = 1090, time(us)=9.08 Measured min. number of cycles = 1170, time(us)=9.75
Measured aver.number of cycles = 2957, time(us)=24.64 Measured aver.number of cycles = 3266, time(us)=27.22

---------------------------------------------- ----------------------------------------------
Algorithm 'Xenia' Algorithm 'Xenia' (http://electronix.ru/forum/index.php?s=&showtopic=114436&view=findpost&p=1180687):
---------------------------------------------- ----------------------------------------------
Results CRC = 81AD Results CRC = 7CBB
Number of samplles = 31 Number of samplles = 31
Number of iterations = 10000 Number of iterations = 10000
Measured max. number of cycles = 5375, time(us)=44.79 Measured max. number of cycles = 5340, time(us)=44.50
Measured min. number of cycles = 469, time(us)=3.91 Measured min. number of cycles = 474, time(us)=3.95
Measured aver.number of cycles = 2947, time(us)=24.56 Measured aver.number of cycles = 2766, time(us)=23.05

---------------------------------------------- ----------------------------------------------
Algorithm 'Neiver' Algorithm 'Neiver' (http://electronix.ru/forum/index.php?s=&showtopic=114436&view=findpost&p=1180944):
---------------------------------------------- ----------------------------------------------
Results CRC = 81AD Results CRC = 7CBB
Number of samplles = 31 Number of samplles = 31
Number of iterations = 10000 Number of iterations = 10000
Measured max. number of cycles = 3014, time(us)=25.12 Measured max. number of cycles = 2540, time(us)=21.17
Measured min. number of cycles = 754, time(us)=6.28 Measured min. number of cycles = 681, time(us)=5.68
Measured aver.number of cycles = 761, time(us)=6.34 Measured aver.number of cycles = 695, time(us)=5.79

---------------------------------------------- ----------------------------------------------
Algorithm 'Ekstrom' Algorithm 'Ekstrom' (http://embeddedgurus.com/stack-overflow/tag/median-filter/):
---------------------------------------------- ----------------------------------------------
Results CRC = 81AD Results CRC = 7CBB
Number of samplles = 31 Number of samplles = 31
Number of iterations = 10000 Number of iterations = 10000
Measured max. number of cycles = 694, time(us)=5.78 Measured max. number of cycles = 637, time(us)=5.31
Measured min. number of cycles = 632, time(us)=5.27 Measured min. number of cycles = 563, time(us)=4.69
Measured aver.number of cycles = 637, time(us)=5.31 Measured aver.number of cycles = 574, time(us)=4.78

---------------------------------------------- ----------------------------------------------
Algorithm 'Klen' Algorithm 'Klen' (http://electronix.ru/forum/index.php?s=&showtopic=114436&view=findpost&p=1180679)
---------------------------------------------- ----------------------------------------------
Results CRC = 81AD Results CRC = 7CBB
Number of samplles = 31 Number of samplles = 31
Number of iterations = 10000 Number of iterations = 10000
Measured max. number of cycles = 1967, time(us)=16.39 Measured max. number of cycles = 7120, time(us)=59.33
Measured min. number of cycles = 822, time(us)=6.85 Measured min. number of cycles = 2580, time(us)=21.50
Measured aver.number of cycles = 1169, time(us)=9.74 Measured aver.number of cycles = 4260, time(us)=35.50

---------------------------------------------- ----------------------------------------------
Algorithm 'VAI' (http://caxapa.ru/170626.html): Algorithm 'VAI' (http://caxapa.ru/170626.html):
---------------------------------------------- ----------------------------------------------
Results CRC = 81AD Results CRC = 7CBB
Number of samplles = 31 Number of samplles = 31
Number of iterations = 10000 Number of iterations = 10000
Measured max. number of cycles = 5707, time(us)=47.56 Measured max. number of cycles = 4627, time(us)=38.56
Measured min. number of cycles = 5650, time(us)=47.08 Measured min. number of cycles = 4549, time(us)=37.91
Measured aver.number of cycles = 5652, time(us)=47.10 Measured aver.number of cycles = 4562, time(us)=38.02


С шаблонами на C++ в обоих компиляторах наблюдаются проблемы.
Keil при наличии шаблонов отказывается делать кроссмодульную оптимизацию. IAR как видно с шаблонами тоже жутко затормозил быстродействие.
Вообщем шаблоны не наш путь. wink.gif

Проекты для микроконтроллера MK60FN1M0VLQ12 здесь:
[attachment=78643:Test_Suite_IAR.zip]
[attachment=78644:Test_Suite_Keil.zip]

Сообщение отредактировал IgorKossak - Aug 1 2018, 17:18
Причина редактирования: [codebox] для длинного кода, [code] - для короткого!
Go to the top of the page
 
+Quote Post
VAI
сообщение Aug 5 2013, 11:59
Сообщение #49


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

Группа: Модераторы
Сообщений: 1 120
Регистрация: 17-06-04
Пользователь №: 37



Херовый у меня алгоритм.. :-(


--------------------
Если зайца бить, его можно и спички научить зажигать
Сколько дурака не бей - умнее не будет. Зато опытнее
Go to the top of the page
 
+Quote Post
p_v
сообщение Aug 1 2018, 16:07
Сообщение #50


Участник
*

Группа: Участник
Сообщений: 68
Регистрация: 16-06-18
Из: СПб
Пользователь №: 105 099



Тема старая, но очень годная. Спасибо за тесты, очень удобно, когда есть на что ориентироваться.

Решил пойти немного другим путем - вместо median считать truncated mean (как в школе на лабораторках). То есть:
  1. Считаем среднее.
  2. Считаем дисперсию.
  3. Считаем среднее только по тем элементам, которые укладываются в заданную дисперсию относительно среднего.
Каких-то особых целей не преследовал, просто захотелось посмотреть как еще можно покоцать резкие отклонения.

CODE
#define ADC_BUFFER_SIZE 32
#define ADC_BUFFER_MASK 0x1F

uint32_t truncated_mean(uint16_t *src, int head, int count)
{
int i = 0;
int idx = 0;

// Collect mean
i = count;
idx = head;
int s_mean = 0;
while (i)
{
i--;
idx--;
idx &= ADC_BUFFER_MASK;
s_mean += src[idx];
}

// add (count >> 1) for better rounding
int mean = (s_mean + (count >> 1)) / count;

// Collect sigma
i = count;
idx = head;
int s_sigma = 0;
while (i)
{
i--;
idx--;
idx &= ADC_BUFFER_MASK;
int val = src[idx];
s_sigma += (mean - val) * (mean - val);
}

int sigma_square_4 = s_sigma * 4 / count;

// Drop big deviations and count mean for the rest
i = count;
idx = head;
int s_mean_filtered = 0;
int s_mean_filtered_cnt = 0;

while (i)
{
i--;
idx--;
idx &= ADC_BUFFER_MASK;
int val = src[idx];

if ((mean - val) * (mean - val) < sigma_square_4)
{
s_mean_filtered += val;
s_mean_filtered_cnt++;
}
}

// Protection from zero div. Should never happen
if (!s_mean_filtered_cnt) return mean;

return (s_mean_filtered + (s_mean_filtered_cnt >> 1)) / s_mean_filtered_cnt;
<div> }


На M3 31 отсчет 1460 циклов по дебагеру. Не знаю на сколько врет, он вообще странный немного.
  • Это медленнее Ekstrom, но зато тут есть встроенный усреднятор (теоретически должен дрожать меньше медианы).
  • Можно переписать вычисление дисперсии и среднего на 1 проход вместо 2, но тогда целочисленная арифметика потянет только 15 12-битных отсчетов
  • Тут еще дополнительные накладные расходы, чтобы читать из кольцевого буфера, лень было вырезать.
  • Кольцевой буфер делался чтобы не возиться с локами и т.п. - размер сделан с запасом, поэтому пока мы обрабатываем данные из одной части, АЦП пишет дальше.


Сообщение отредактировал IgorKossak - Aug 1 2018, 17:07
Причина редактирования: [codebox] для длинного кода, [code] - для короткого!
Go to the top of the page
 
+Quote Post
p_v
сообщение Aug 3 2018, 15:16
Сообщение #51


Участник
*

Группа: Участник
Сообщений: 68
Регистрация: 16-06-18
Из: СПб
Пользователь №: 105 099



https://pastebin.com/xpLbyWN6 - двухпроходный вариант, ~ на 20% быстрее предыдущего, + настраиваемый порог срабатывания. 2*sigma многовато оказалось. 1.1*sigma поправдоподобнее.
Go to the top of the page
 
+Quote Post

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

 


RSS Текстовая версия Сейчас: 22nd July 2025 - 10:23
Рейтинг@Mail.ru


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