LINUX.ORG.RU

Посоветуте оптимальный алгоритм поиска медианы

 , ,


2

2

есть stm32 c ацп, к которому хочется привернуть фильтр шума. Обычно советуют медианный.

Ну ок, допустим мне надо искать медиану для 5-9 отсчетов (uint16_t). Как это сделать наиболее быстрым образом? В идеале - просто ссылка на годную библиотеку, ну и какие-то бенчмарки минимальные, желательно для ARM Cortex M0-M3.

PS Вики и stackoverflow в общих чертах смотрел.

★★★★★

Ответ на: комментарий от Vit

Там просто алгебра. Где затык? (м.б. это: sum(x*sum(x[j],j,1,n),i,1,n)=(sum(x,i,1,n))^2; )

12бит в квадрате =24, 32-24=8бит осталось на сумму, 8бит=256 слагаемых, так?

А почему интервалы должны не пересекаться, и как это мешает уползанию фазы? (тогды бы можно было считать «скользящую сигму», это еще быстрее)

anonymous ()
Ответ на: комментарий от silw

спасибо, на livestats я тестил.
вопервых у него погрешность получается таже что и у меня +- (в точной цифре не уверен, я ее тестил на меньшем наборе данных).
вовторых это всеже скользящая медиана, что в моем случае крайне не оптимально - я ведь уже имею медиану по большинству чанков, зачем же все поновой пересчитывать.
или ручь о том, что ей можно подсовывать медианы чанков? такой алгоритм не пробовал, проверю в выходные. моих знаний в матане нехватает, чтобы без теста оценить насколько такой подход будет рабочим

genryRar ★★ ()
Ответ на: комментарий от genryRar

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

silw ★★★★★ ()

У меня был очень простой алгоритм, я все значения складывал в массив, потом на контрольном значении сортировал его вставкой и брал медиану. Бенчмарков нет, но мне хватало, сенсор опрашивал каждые три секунды, отсылка раз в минуту, массив был 60/3=20 элементов. Хитрость была, если сенсор возвращал нереальное запредельное значение, то я писал в текущую ячейку медиану с прошлой минуты.

anonymous ()
Ответ на: комментарий от Vit

У меня так написалось:

int n = 8;	//ширина скользящего окна
int kk = 2 * 2;	//квадрат множителя порога отсечки (2*sigma)
int ns1 = 0;	//n*сумма (биты в хвосте не теряем)
int ns2 = 0;	//n*сумма_квадратов (n=[1-16], т.к. n*n квадратов от 12бит)
int xb = 0;	//предыдущий х
int xm = 0;	//среднее обрубленное

int filter(int x)
{
    int d, s1;
    ns1 += x - ns1 / n;
    ns2 += x * x - ns2 / n;
    d = (x - xm) * n;
    s1 = ns1 / n;
    if (d * d < kk * (ns2 - s1 * s1)) //отклонение^2 < (k*сигма_с_шумом)^2
        xb = x;
    return xm += (xb - xm) / n;
}
М.б. надо считать сигму без шума? Не застрянет ли на старте?

Данных с сенсора нет?

anonymous ()
Ответ на: комментарий от anonymous

Куча ошибок. 1) Действительно 15 слагаемых, а не 256 ((4бит_сумма + 12бит_отсчет)^2=32бит, по тексту в if-е) 2) Перепутал сумму и среднее. 3) Если считать сугму_без_шума, то будет застревать, и не только на старте. С сигмой_с_шумом тоже может застревать на сигнале типа ШИМ (при окне < ширины меандра). Заплатка медленно подтягивает xb к зашумленному_среднему.

int n = 8;	//ширина скользящего окна
int kk = 2 * 2;	//квадрат множителя порога отсечки (2*sigma)
int s1 = 0;	//скользящая сумма
int s2 = 0;	//скользящая сумма_квадратов
int xb = 0;	//предыдущий х
int xm = 0;	//среднее обрубленное

int filter(int x)
{
    int d;
    s1 += x - s1 / n;
    s2 += x * x - s2 / n;
    d = (x - xm) * n;
    if (d * d < kk * (n * s2 - s1 * s1))	//(n*отклонение)^2 < (k*n*сигма_с_шумом)^2
        xb = x;
    else
        xb += (s1 / n - xb) / (3 * n);	//заплатка от застревания
    return xm += (xb - xm) / n;
}
Итог похож на RC-фильтр.

anonymous ()
Ответ на: комментарий от tailgunner

Казахский.

Так вы спец. по ТАУ. (?) А выдаете критику тупизма, а не решения.

Если считать, что «s+=...s...;» - два обращения, таких строк 6, значит на операции «+-*/» ничего не осталось.

М.б. несколько измерений за раз считать.

anonymous ()
Ответ на: комментарий от anonymous

Я как-то пропустил момент переключения. Ну, если не надо проходить по массиву - можно попробовать распределить переменные по регистрам и за счет этого сэкономить на доступах к памяти.

tailgunner ★★★★★ ()
Ответ на: комментарий от tailgunner

Ты все проспал. Проехали медиану уже.

Для эмбедов по медиане замеры есть тут https://electronix.ru/forum/index.php?s=&showtopic=114436&view=findpo... Раза в полтора быстрее чем truncated mean. Зато тут с усреднятором и алгоритм прямой как рельса.

Vit ★★★★★ ()
Ответ на: комментарий от tailgunner

Я так и не понял, чем тебе плоха медиана,

Да ничем особо не плоха. Просто не было полностью готовых библиотек под platformio, а вникать в код Ekstrom чтобы его нормально переписать - мне было влом. 2-проходный truncated mean почти такой же быстрый, но тупо проще в несколько раз.

Vit ★★★★★ ()

Я, конечно, чайник, но нельзя как нибудь так Что самое простое? Найти медиану для трёх элементов. Как найти медиану для 9 элементов? Надо три раза найти медиану для 3х элементов, а потом найти медиану из трёх медиан. Медиана из 3х находится за 3 сравнения, следовательно нам нужно 4*3 = 12 сравнений, чтобы найти медиану из 9 элементов.

nguseff ()
Ответ на: комментарий от nguseff

При этом используя результаты предыдущих шагов получается ещё меньше сравнений. На каждое следующее значение нужно сделать всего 6 сравнений, так как две медианы уже посчитаны.

nguseff ()
Ответ на: комментарий от nguseff

Неа, 19 сравнений нужно.

Потому что медиана девяти элементов никак не равна медиане трех групп по три элемента!

Пример. Ряд 1,9,8,3,5,8,1,4,4. Сортируем: 1,1,2,3,4,5,8,8,9, медиана - 4. В группах же медианы: 8, 5 и 4, т.е. медиана групп равна пяти, а не четырем!

anonymous ()
Ответ на: комментарий от Vit

До 8 увеличили пока что

Можно сделать самонастраивающееся. Как-то так:

int m = 0, a = 1, b = 1, c = 30;
int fa(int x)
{
    int y, z;
    y = a + b;
    z = y + c;
    m = (m * c + x * y) / z;
    a = c * y / z;
    return m;
}
Чтений из памяти 4, записей 2, сложений,умножений 8. b и c - const, зависят от параметров сигнала и шума, a - подстраивается.

anonymous ()
Ответ на: комментарий от Vit

Просто не было полностью готовых библиотек под platformio, а
вникать в код Ekstrom чтобы его нормально переписать - мне было
влом.

Помойму проще:

http://ndevilla.free.fr/median/median/src/optmed.c

Это если не надо менять кол-во элементов на ходу.

А вообще, можно сделать накопление сигмы «RC-фильтром», и накопление среднего - другим «RC фильтром» с управляемой «постоянной времени»: типа, если отклонение > 3 сигма - длинное время, если меньше - короткое.

Так накапливают «фон» в детекторах движения.

shkolnick-kun ★★★★ ()
Последнее исправление: shkolnick-kun (всего исправлений: 1)
Ответ на: комментарий от shkolnick-kun

Ты про кольцевой буфер забыл. С фиксированным размером еще попадаешь на копирование.

От кроилова на скольжении профит только если интервалы пересекаются. То есть, экономии на бюджете проца не будет.

Vit ★★★★★ ()
Ответ на: комментарий от Vit

Ты про кольцевой буфер забыл. С фиксированным размером еще попадаешь на копирование.

В таком случае можно попробовать:

Torben's method

This method was pointed it out to me by Torben Mogensen.

It is certainly not the fastest way of finding a median, but it has the very interesting property that it does not modify the input array when looking for the median. It becomes extremely powerful when the number of elements to consider starts to be large, and copying the input array may cause enormous overheads. For read-only input sets of several hundred megabytes in size, it is the solution of choice, also because it accesses elements sequentially and not randomly. Beware that it needs to read the array several times though: a first pass is only looking for min and max values, further passes go through the array and come out with the median in little more time that the pixel_qsort() method. The number of iterations is probably O(log(n)), although I have no demonstration of that fact.

http://ndevilla.free.fr/median/median/src/torben.c

shkolnick-kun ★★★★ ()
Последнее исправление: shkolnick-kun (всего исправлений: 2)
Ответ на: комментарий от Vit

Еще можно попробовать вот так:

#include <sdtint.h>
#include <stdbool.h>

/*Используем с99 и знаковые типы, так проще*/

const int32_t fast = 100; /*Постоянная времени быстрого фильтра*/
const int32_t med  = 10;  /*Постоянная времени промежуточного фильтра*/
const int32_t slow = 1;   /*Постоянная времени медленного фильтра*/

const int32_t denom = 1000;

int32_t mean = 0;  /*Оценка среднего*/
int32_t noise = 0; /*Оценка шума*/

int16_t tr_mean(int16_t input)
{
    int32_t delta = (int32_t)input - (mean / denom);
    int32_t nz = delta * delta; /*Оцениваем квадрат отклонения*/
    
    alpha = (nz > noise * 9)?med:fast; /*Изменения < 3s будем копить быстрее*/
    mean  += delta * alpha;

    noise += (nz - (noise / denom)) * slow; /*Оценку шума делаем медленнее всего*/
    return (int16_t)(mean / denom);
}

Кольцевой буфер тут не нужен. Но есть неилюзорный шанс переполнения...

shkolnick-kun ★★★★ ()
Последнее исправление: shkolnick-kun (всего исправлений: 4)
Ответ на: комментарий от shkolnick-kun

Забей. Уже все написали, закоммитили, и уже переписываем заново из-за переделки логики прерываний :). Ковырять глубже уже просто не интересно - дальше двинулись.

Vit ★★★★★ ()