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

Немного из другой оперы, но может будет интересно - онлайн подсчет оценки медианы, без сортировки значений.

Может скользящее среднее: s=(7*s+x)/8 подойдет?

Нет, уровень ардуины и младшей дошкольной группы мне не нужен.

Там для больших объемов. Не то.

Теоретически есть select в std. Он безусловно весьма быстр, но у меня есть сомнения в его оптимальности конкретно для моей архитектуры и размеров.

Хотелось бы услышать тех, кто имел дело с практикой.

а можно ссылку на описание метода, а то что то не найду

https://en.wikipedia.org/wiki/Quickselect
https://stackoverflow.com/a/810905/1031804

Если ничего не перепутал, это оно. И там внутри скорее всего еще партиционирование есть (но не заглядывал, врать не буду).

На SO еще выше расписано через HEAP на 27 элементов и бенчмарки, но глядя на стиль программирования автора, у меня есть некоторые сомнения в результах.

Для эмбедов много ссылок есть сюда https://embeddedgurus.com/stack-overflow/2010/10/median-filtering/, но как-то тоже не довелось на гитхабе увидеть библиотеки с реализацией от Phil Ekstrom.

Спасибо, сейчас почитаю

Будешь смотреть - обрати внимание, что чаще фильтр делают для картинок (2D), и самые выгодные оптимизации касаются именно пересчета соседей всякими скользящими методами.

Но у меня АЦП (1D), выборки не пересекаются (чтобы фазу не колбасило), и надо просто каждый раз по новой выборке заново посчитать 1 число.

я уже реализовывал этот фильтр для одномерного случая, но для меня скорость работы фильтра была не важна и использовал я обычный std::sort. получалось примерно 0.11 -0.13 с на 50000 точек

Мне надо 5-9 точек, но бюджет порядка микросекунды :)

не 50000 точек - это весь обрабатываемый объем, среднее окно было 256 точек, т.е было 2,2*10^-6 c на точку, при процессоре intel i5. Если выше описанный метод шустрее стандартного sort - должен быть прирост

Но у меня АЦП (1D)

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

3.3.1 Averaging samples ... In STM32L0 and STM32L4 microcontrollers, averaging can be performed by using the hardware oversampling feature: the ADC performs built-in hardware averaging according to configurable parameters (number of samples to average and final right bit shift of result). The advantage of averaging is to improve ADC precision without any hardware changes. The drawback is that the conversion speed is lower as well as the frequency response (it is equivalent to decreasing effective sampling frequency).

https://www.st.com/content/ccc/resource/technical/document/application_note/g...

апаратно вычисляется среднее арифметическое

Насколько я понял, нужна медиана, а не среднее.

сортируешь, потом выбираешь средний элемент из массива

5-9 отсчетов

Вангую, что в этом случае быстрее и проще тупо сортировкой и взятием среднего элемента или среднего значения 2х центральных обойтись. Все эти решения через медианы медиан это для сильно бôльших массивов данных. А с хипом так вообще решения для потока.

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

если схемотехника - говно, никакая медиана и прочие цифровые фильтры не помогут

Мне надо 5-9 точек, но бюджет порядка микросекунды :)

Для 5-9 точек могут быть полезны аппаратные сортировочные сети с особыми свойствами ©. У меня была разработка для ранговой фильтрации 9 точек ТВ сигнала в реальном времени 13 МГц на основе аналогичной сети.
Но это был «большой железный ящик» :)

Глянь: оптимальная 9-элементная сортировочная сеть ©. Там в тексте ссылки на Кнута и его соратников.

Ещё: поиск k-ой порядковой статистики за линейное время © для фильтрации гибче и полезнее медианы.

Ежели по времени допустимо, то аппаратно сеть сортировки/медианы можно сократить до 2-х параллельных слоёв работающих «вход < — > выход».

если схемотехника - говно, никакая медиана и прочие цифровые фильтры не помогут

Схемотехника может быть идеальна, но сам исходный сигнал сложный, например, селекция цели на фоне импульсных помех.

А ранговая фильтрация — самое эффективное средство супротив импульсных помех и борьбы с тяжёлыми хвостами спектра мощности шума.

5-9 отсчетов

Сгенерировать if-else? Кажется, что на таком маленьком количестве никакие хитрые алгоритмы выигрыша не дадут, а только замедлят процесс.

Въехал в тему. Штука хорошая. Но тут вспомнил как мы на лабах по физике тупо считали среднее потом sigma, выкидывали все данные abs(x - avg) < 3 * sigma и снова пересчитывали среднее. Все делается за O(n) но нужен FPU или на крайняк считать в fixed point ах.
Принципе на микроконтроллерах главное стабильная скорость обработки какая бы она ни была. А то на удачных данных ты посчитаешь за 2мкс а на неудачных за 10мкс. Итого неровная квантизация выборки данных. Так что лучше от алгоритмов где есть много условных операций отказаться. Или выполнять их по прерыванию таймера но тогда будет не максимальная скорость обработки.

http://we.easyelectronics.ru/STM32/filtraciya-zvuka-na-stm32-s-ispolzovaniem-...

А есть чем сгенерить? Я только за.

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

Угу, почитав compare networks я тоже стал подозревать, что должны быть генераторы сишечки, но тоже не нашел.

А есть чем сгенерить?

Проще руками, например, компараторные модули сети заменяешь на программные if-else, а каждая линия — поток данных.

Примерчик с кодом: сеть обменной сортировки со слиянием Бэтчера ©.

Въехал в тему. Штука хорошая. Но тут вспомнил как мы на лабах по физике тупо считали среднее потом sigma, выкидывали все данные abs(x - avg) < 3 * sigma и снова пересчитывали среднее. Все делается за O(n) но нужен FPU или на крайняк считать в fixed point ах.

На truncaded mean смахивает. Только не очень понятно, как выбирать sigma чтобы обеспечить нужную сходимость, и какой будет профит по сравнению с медианой.

Может есть какая-то «методичка», где написаны рекомендации по количеству отсчетов и т.п.?

под sigma тут имеется ввиду https://en.wikipedia.org/wiki/Standard_deviation#Discrete_random_variable
И по графику https://upload.wikimedia.org/wikipedia/commons/8/8c/Standard_deviation_diagra... 3 * sigma это многовато. Лучше 2 * sigma
Корни тут кстати очевидно можно не считать а сравнивать (x - avg) ** 2 < 4 * sigma ** 2
Сходимость тут хорошая тк это вообще стандартный научный метод отсева данных экспериментов тогда как медиана все равно будет скакать как по кочкам.
По профиту не знаю. Смотря что на stm32 будет быстрее работать.

Блин, сказал бы сразу что это среднеквадратичное отклонение.

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

Короче, в кортексе M3 самая жопа в доступе к памяти - 10 циклов. Поэтому твой вариант лучше, т.к. там чистых 2 прохода чтения и все.

Причем сейчас я пытался сначала ring buffer копировать в линейный, а надо сразу из кольцевого буфера считать.

Математика там целочисленная, т.к. данные 12 бит и даже и если возводить в степень, то на буфер из 16 элементов хватит без переполнений.

Завтра будем тестить.

Набрехал, надо 3 прохода.

малость пооффтоплю:
недавно столкнулся с такой проблемой - нужно считать перцентили (или хотябы медиану) по массиву данных, разбитому на произвольные куски. склеить все вместе невыйдет - данных дофига. для большинства кусков медиана уже посчитана заранее. думал это довольно банальная задача - но готового решения или хотябы внятной методологии не нашел.
мой лучший вариант расчета медианы дает погрешность в 1.6% по верхнему квартилю (0.75) и 8.5% максимальную. алгоритм предельно топорный - медианы всех частей принимаем за новый массив, считаем по нему медиану и среднее и берем среднее от этих двух значений. тест проводился по 10кк наборам данных длинной от 10 до 10к элементов случайных целых чисел от 0 до 1000.
может есть у кого идеи как улучшить?

Eigen не хватает?

https://ru.wikipedia.org/wiki/Eigen

Перебор. Пока осиливаем truncated mean.

в eigen там тупо шаблоны, нечего компилять даже(тупо инклюд добавить), и куча туториалов в инете

если не пользовался попробуй

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

Один резистор и один конденсатор на входе АЦП. Сложность алгоритма - O(0). Время на расчёты - 0мкс.

склеить все вместе

Глянул по диагонали. Слишком наворочено и неудобно как внешнюю зависимость в PlatformIO интегрировать.

Все равно из-за оптимизаций приходится ручками шаманить, иначе медленно совсем.

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

надо скорость и контроль — пиши на асме. или кишка тонка?

Медианы это для отсева шума экстремумов, ты задачу не понимаешь

может есть у кого идеи как улучшить?

Покопался немного в своём архиве:

Хуанг Т.С. Быстрые алгоритмы в цифровой обработке изображений. Преобразования и медианные фильтры ©.
Там был алгоритм скользящей модификации медианы по большой выборке в прямоугольном окне изображений, иногда можно применить и для одномерного случая, типа медиана из медиан.

Ещё глянь: Fast median search: an ANSI C implementation ©.

Примеры кода для самостоятельной модификации под задачу:
mediator.c ©,
quick_median.c ©,
nth_element ©.

Откуда 3 прохода? Сигма(с.к.о.) считается за один (в две переменные).

Я может чего не понимаю, но cortex m0-m3 не содержат ни аппаратной плавучки, ни DSP-инструкций, как m4 в stm32f4 соответственно. Поэтому вопрос об оптимизации под архитектуру звучит странно. Просто берете и реализуете подходящий алгоритм на C.

тебе как раз подойдет эта штука. Ссылка на статью там же есть. Сам активно использую для поиска медианы по матрице с 10^8-10^10 значений. Просто скармливая элементы по рядам.

Дык надо сначала среднее посчитать?

1. Считаем среднее
2. Считаем СКО
3. Считаем среднее на тех числах, которые вписались в 2*сигма

У меня три прохода получилось.

Вот ссыль: http://ndevilla.free.fr/median/median/index.html

Выбирай себе любой метод, смотря на задачи.

Я, как и анон, тоже что-то упустил. Откуда три прохода? Речь же шла только о медиане, чем подсчет только авеража не угодил? Почему и как ты пришел к именно таким кол-ву и порядку проходов?

За один проход считаем s1-сумма и s2-сумма квадратов, потом sigma=sqrt(n*s2-s1*s1)/n;

Вдогонку: Надо в формулу для сигма подставить формулу для среднего и раскрыть скобки с квадратами.

Да, нашел: https://en.wikipedia.org/wiki/Algorithms_for_calculating_variance. Я так и не понял, почему базовую формулу корректно приводить к такой, но так действительно делают.

Правда тогда для 12-битного АЦП целочисленная 32-битная арифметика потянет около 15 отсчетов без переполнения. С другой стороны, мне больше и не надо.

Надо будет попробовать если руки дойдут, годная тема. Спасибо.