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

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

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

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

Из сниппетов.

Ну и общий случай.

спасибочки, ознакомлюсь

За микросекунду можно успеть ручками отсортировать и выбрать центральный элемент

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

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

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

Да, не точно посчитал длину фильтра. Больше влезет.

Переписал, процентов на 25 быстрее стало https://github.com/speedcontrols/ac_sc_grinder/commit/b727033e0cd0ff2e659e596....

Что ты все выдумываешь что-то? Уже давно в сниппетах все необходимое есть, я примеры кода выше кидал.

P.S. Как видно по первой ссылке, я opt_med9 вообще на тупом STM8 запускаю! И ничего, довольно шустро обновляет показания АЦП на экранчике...

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

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;
}

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

Куча ошибок. 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;
}

За микросекунду можно успеть ручками отсортировать и выбрать центральный элемент

При тактовой частоте 120МГц и при цене обращения к памяти 10 тактов?

стандартный sort на маленьких выборках может быть медленнее самописной процедуры

Т.е. можно всего 12 чтений/записей? Чето жестко слушком.

Офф. А у вас Приозерск казахский или карельский?

Т.е. можно всего 12 чтений/записей? Чето жестко слушком.

О чем и речь. Я подозреваю, что в микросекунду влезть просто невозможно.

Казахский.

Казахский.

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

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

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

Так вы спец. по ТАУ. (?)

Я системный программист.

...

ЯННП.

ЯННП.

int filter(int) (текст выше) - зашкаливает за 1мкс. при грубом подсчете. Что делать - не понимаю.

int filter(int) (текст выше)

Это точно медианный фильтр?

зашкаливает за 1мкс. при грубом подсчете.

Ну я и говорю - мне кажется, что задача нерешаема.

Это точно НЕ медианный фильтр. Вит, переключился на truncated_mean, это, как понимаю, среднее из значений внутри интервала 3*sigma.

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

На М3 цена кажись 5 тактов. 10 - это если копирование.

В общем, 2-проходная версия 4 отсчета перебирает за 160 тактов.

Мне кажется, что с 5 тактами на обращение можно сделать и тупо вставку в отсортированный список.

Может и можно. А смысл в чем?

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

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

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

Проехали медиану уже.

Я так и не понял, чем тебе плоха медиана, если она вписывается во время, но тебе виднее.

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

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

А давайте забаним этого топикстартера, чисто из любви к искусству системного программирования! :)

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

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

Каюсь не совсем медиана получается, а вернее даже совсем не она.

Ну если медиана еще актуальна - вот для 5 элементов (дело было вечером, делать было нечего)

https://wandbox.org/permlink/M5HWjcElceiXShAn

Раскуроченная сортировка слиянием

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

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

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

Какое, нахрен, среднее? Откуда у тебя возьмется гауссов шум? Дробовой + пуассонов != гауссов!

Только медиана и спасет!

Ну да 9! больше, чем 2^12 поэтому за 12 физически нельзя. 19 получается как раз минимум.

А аналогичный онлайн из С в асм для АРМа есть?

https://godbolt.org/

Спасибо.

Почему то, все версии gcc (кроме арм64) генерят

bl      __aeabi_idiv

До 8 увеличили пока что. 4 маловато оказалось.

До 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;
}

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

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

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

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

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

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

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

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

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

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

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

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

#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);
}

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

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

как и обещал, провел тесты.

arrMedianLivestats1: q=(0.2%  0.5%  0.9%)	max=5.0%
    arrMedianParts1: q=(0.4%  0.9%  1.6%)	max=8.0%