LINUX.ORG.RU

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

 , median,


2

2

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

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

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

★★★★★

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

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

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

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

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

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.

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

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

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

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

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

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

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

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

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

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

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

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...

anonymous
()

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

Harald ★★★★★
()

5-9 отсчетов

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

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

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

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

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

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

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

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

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

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

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

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

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

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

quickquest ★★★★★
()

5-9 отсчетов

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

xaizek ★★★★★
()

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

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

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

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

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

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

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

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

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

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

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

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

под 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 будет быстрее работать.

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

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

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

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

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

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

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

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

Vit ★★★★★
() автор топика

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

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

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

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

missxu
()

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

quickquest ★★★★★
()

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

Vit ★★★★★
() автор топика
Вы не можете добавлять комментарии в эту тему. Тема перемещена в архив.