Пролог
В программировании микроконтроллеров порой приходится прибегать к такому программному компоненту как медианный фильтр. Вот вам яркий пример. Вы измеряете расстояние UWB трансиверами между радиопередатчиками и вдруг у вас образуются единичные редкие отбросы, которые никак не характеризует измеряемую величину. Или вы опрашиваете емкостную кнопку и от скачка питания она тоже может дать редкие, но высокие значения отдельного семпла. Классическое решение для исключения таких отбросов - это медианный фильтр (МФ). Метафорично выражаясь медианный фильтр - это острая бритва для отсечения выпирающих семплов.
Ситуация осложняется тем, что на курсах ЦОС(DSP) в ВУЗах обычно про медианный фильтр даже не вспоминают, отводя всё время на обзор FIR (КИХ) и немного IIR (БИХ) фильтров. Хотя достоинство МФ в том, что у него всего один конфиг: порядок фильтра. Плюс медианный фильтр (МФ) как раз и предназначен для того, чтобы исключать эти самые одиночные отбросы сигнала (промахи).
Это настолько важная в программировании задача, что её можно даже встретить на LeetCode. В связи с этим я произвел детальный разбор решения LeetCode задачи 480. Sliding Window Median в контексте реализации МФ на языке программирования Си.
Постановка задача
Реализовать оптимизированный по быстродействию медианный фильтр на языке программирования Си, чтобы его можно было применять в микроконтроллерных прошивках. Фильтр должен работать в потоковом режиме. На каждом отсчете обрабатывать семпл на входе и выдавать семпл на выходе.
Определения
медиана — это центральный элемент отсортированного массива нечетной длинны (или среднее арифметическое двух центральных элементов в случае массива четной длинны). Медиана определяется, как среднее значение упорядоченного списка.
порядок фильтра (размер окна, число k) - количество отсчетов, которые фильтр анализирует на каждом такте. Проше говоря сколько памяти в этом фильтре. Порядок фильтра еще называют размером окна.
семпл - целое число. Отсчеты ADC. В конечном счете семпл - это какая-либо физическая величина: расстояние, напряжение, ток, освещенность и т.п.
абстрактная структура данных - простыми словами это способ организации хранения, помещения и извлечение данных в памяти компьютера.
нечетное число - это число, которое не делится на два в множестве целых чисел. Пример нечетных чисел: 1; 3; 5; 7; 9; 11; и т д
Двоичное дерево (Бинарное дерево) — иерархическая структура данных, в которой каждый узел имеет не более двух потомков (детей). Как правило, первый называется родительским узлом, а дети называются левым и правым наследниками.
Бинарная куча (она же приоритетная очередь) - это абстрактный тип данных, двоичное дерево. У кучи есть только три легальные операции: добавить элемент в кучу, посмотреть на корень кучи (максимальный или минимальный элемент), вытащить максимальный (или минимальный) элемент, сбалансировать кучу. При работе с бинарной кучей удалять можно только то, что лежит в корне кучи. Вот, пожалуй, и всё.
Балансировка медианного фильтра - это перекладывание из одной кучи в другую через вершину.
Хэш-таблица - это абстрактная структура данных с мгновенной вставкой элемента в хранилище и мгновенным извлечением или проверкой элемента из хранилища.
Сдвиговый регистр (он же звено задержки)- FIFO постоянного размера.
Реализация
Текстовое описание алгоритма
Интуитивно понятно, что на каждом шаге надо добавить новый элемент (семпл), удалить старый элемент, выполнить вычисления и выдать результат.
Основная идея в том чтобы использовать две кучи: максимальную и минимальную.
Название кучи |
для чего |
max_heap |
хранить половину с меньшими элементами |
min_heap |
хранить половину с большими элементами |
Обслуживая эти две кучи мы всегда можем получать быстрый доступ к средему элементу. Схематически это выглядит так.
Max_heap позволяет нам получить доступ к наибольшему элементу из меньшей половины элементов. Min_heap позволяет нам получить доступ к наименьшему элементу в половине , где все элементы больше чем в нижней куче. Для простоты min-heap назовем large, a max-heap small
Учитывая, что в куче large все элементы больше или равны всем элементам, чем в куче small, то можно даже перерисовать бинарные кучи прямоугольными треугольниками. Вот так.
Когда размер окна нечетный (3; 5; 7) медиана - это верхушка max_heap. Когда k-четное , медиана - это среднее между верхушками(коренями) двух куч.
На каждом шаге надо делать балансировку куч. Надо следить чтобы количество элементов в верхней куче и в нижней куче было одинаково или отличалось на 1 (в случае нечетного порядка фильтра).
Инициализация
До начала работы с фильтром надо проинициализировать max_heap, чтобы хранить меньшую половину, min_heap - чтобы хранить большую половину. Надо определить hash таблицу (delayed) для отложенного удаления чисел, которые были логически удалены, но не удалены физически из куч.
Какой вместительности должны быть эти две кучи? В идеале бесконечного. Однако так не бывает. В случае сильного разрастания куч можно вручную удалить ненужные более элементы согласно данным из хэш-таблицы (delayed). Для этого будет отдельная функция flush.
Ради определенности перед первым пуском надо загнать в фильтр k нулей. Нулей должно быть столько же сколько составляет размер фильтра. Это нужно , чтобы поведение фильтра было предсказуемым при первом пуске.
Добавление элемента х
Вот поступило число. В какую кучу его положить? Когда новое число добавляется, оно сравнивается с наименьшим числом из кучи large. Если x меньше или равен root.large, то этот элемент x должен принадлежать куче small и добавляется в small. Если root.large < x, то х в large. Как только элемент добавлен вызывается алгоритм rebalance, чтобы убедиться, что обе кучи содержат примерно одинаковое количество элементов, допуская максимальную разницу в размере только в один элемент.
Перебалансировка медианного фильтра
Метод перебалансировки гарантирует, что количество элементов в small и large кучах равно либо
small.size = large.size, при k четном
small.size + 1 = lagre.size , при k нечетном
Это достигается путем переноса верхнего элемента из одной кучи в другую при нарушении условия разности размеров.
Балансировка куч выполняется после каждой операции добавления или удаления семпла, чтобы соблюдалось свойство половинного размера двух куч относительно размера окна (порядка фильтра).
Вычисление медианы
Метод определения медианы вычисляет медиану на основе размера куч. Если K нечетное, медиана равна верхнему(корневому) элементу smal.root. Если K четное, что медиана равна среднему арифметическому от корней куч
median = small.root, при k нечетном ( 1; 3; 5; 7 ...)
median = 0.5*(small.root+large.root), при k четном (2; 4; 6; 8 ... )
Удаление чисел
На каждом шаге нам надо добавить новый элемент и удалить старый. Выявлять самый старый элемент придется при помощи программного сдвигового регистра размера k+1. После добавления семпла мы также спрашиваем сдвиговый регистр какой элемент вывалился из регистра. Его и предстоит удалить. Размер сдвигового регистра должен быть на 1 больше, чем размер окна.
Чтобы предотвратить разрастания кучи применяется подход отложенного удаления с использованием hash таблиц. При использовании этого метода элементы, которые надо удалить, условно, помечаются словом "откладываются", а фактическое удаление из кучи происходит лениво, то есть мы удаляем элемент только тогда, когда он оказывается наверху той или иной кучи.
При перемещении окна надо удалять числа, выходящие за его пределы. Чтобы избежать дорогостоящих по времени операций удаления, используется подход ленивого удаления, при котором число помечается для удаления в hash таблице с отсрочкой, но не удаляется немедленно из кучи. Фактическое удаление происходит в методе prune, который запускается, когда удаляемое число оказывается наверху одной из куч.
Можно конечно же на каждой итерации удалять элемент принудительно вручную из куч методом откапывая числа в нужном направлении. Просто последовательно перекладывая элементы из кучи в кучи пока этот удаляемый элемент внезапно не окажется корнем в откапываемой куче. Поcле удаления все равно делать балансировку, которая всё восстановит (закопает обратно). В результате не будет потребности в динамической памяти для куч, а хеш-таблицу и вовсе можно будет выбросить за ненадобностью. Однако это крайне продолжительная операция. По сути та же самая сортировка, что и в классической реализации на сортировке элементов окна.
В конце концов можно удалять элементы прямо из середины бинарных куч. Поиск по значению - это O(log(N)) поменяем значения в массива и затем сделает балансировку поддерева за те же O(log(N)). Таким образом и удалим элемент из кучи. Тогда тоже не нужна хеш-таблица, плюс размер куч всегда будет ограничен сверху заранее известным значением.
Важные замечания к алгоритму
1--Надо сделать виртуальные размеры куч: small_size и large_size. При добавлении и извлечении элементов увеличивать и уменьшать надо именно эти отдельные переменные. Реальные же размеры куч не должны участвовать, как управляющие данные для алгоритма. Они живут своей жизнью.
2--Мы всегда должны внимательно наблюдать за верхушкой куч. Как только там окажется число из реестра на ликвидацию, так сразу это число надо вытащить и выбросить, обновив при этом сам список на удаление (т.е. удалить элемент из хэш-таблицы).
3--По хорошему кучи должны быть бездонными. Со временем любой размер в конце концов переполнится.
Исходный код
На языке программирования Си это реализовать особенно трудно, так как в Cи нет готовых реализаций абстрактных типов данных: бинарные кучи, хэш-таблицы и сдвиговые регистры. Всё это предварительно надо написать или найти и крепко протестировать. И уже на основе этого как из кирпичей строить медианный фильтр.
Скрытый текст
#include "median_filter_fast.h"
#include <stdlib.h>
#include "circular_buffer_dword.h"
#include "code_generator.h"
#include "float_utils.h"
#include "hash_table_s8.h"
#include "log.h"
#include "max_heap.h"
#include "min_heap.h"
#include "utils_math.h"
COMPONENT_GET_NODE(MedianFilterFast, median_filter_fast)
COMPONENT_GET_CONFIG(MedianFilterFast, median_filter_fast)
static bool median_filter_fast_init_custom(void) {
bool res = true;
return res;
}
static bool median_filter_form_min_heap_config(const MedianFilterFastConfig_t* const Config,
BinHeapConfig_t* const ConfigMin) {
bool res = false;
if(ConfigMin) {
if(Config) {
ConfigMin->array = Config->tempLarge;
ConfigMin->capacity = Config->bin_heap_size;
ConfigMin->name = "MinHeapSmall";
ConfigMin->num = 10;
ConfigMin->valid = true;
res = true;
}
}
return res;
}
static bool median_filter_form_hash_table_s8_config(const MedianFilterFastConfig_t* const Config,
HashTableS8Config_t* const ConfigHash) {
bool res = false;
if(ConfigHash) {
if(Config) {
// ConfigHash->size = Config->hash_table_s8_size;
// ConfigHash->Memory = Config->HashTableS8Memory;
ConfigHash->name = "HashTableS8";
ConfigHash->num = 12;
ConfigHash->valid = true;
res = true;
}
}
return res;
}
static bool median_filter_form_max_heap_config(const MedianFilterFastConfig_t* const Config,
BinHeapConfig_t* const ConfigMax) {
bool res = false;
if(ConfigMax) {
if(Config) {
ConfigMax->array = Config->tempSmall;
ConfigMax->capacity = Config->bin_heap_size;
ConfigMax->name = "MaxHeapSmall";
ConfigMax->num = 11;
ConfigMax->valid = true;
res = true;
}
}
return res;
}
static bool prune_large(MedianFilterFastHandle_t* const Node) {
bool res;
bool loop = true;
bool ok = false;
while(loop) {
uint32_t cnt = min_heap_size(&Node->Large);
if(cnt) {
int32_t root_large = 0;
res = min_heap_peek(&Node->Large, &root_large);
if(res) {
loop = hash_table_s8_check_ll(&Node->ToDelete, (int8_t)root_large);
if(loop) {
res = min_heap_delete_root(&Node->Large);
res = hash_table_s8_pull_ll(&Node->ToDelete, (int8_t)root_large);
ok = true;
}
}
} else {
loop = false;
}
}
return ok;
}
static bool prune_small(MedianFilterFastHandle_t* const Node) {
bool res;
bool ok = false;
bool loop = true;
while(loop) {
uint32_t cnt = max_heap_size(&Node->Small);
if(cnt) {
int32_t root_small = 0;
res = max_heap_peek(&Node->Small, &root_small);
if(res) {
loop = hash_table_s8_check_ll(&Node->ToDelete, (int8_t)root_small);
if(loop) {
res = max_heap_delete_root(&Node->Small);
res = hash_table_s8_pull_ll(&Node->ToDelete, (int8_t)root_small);
ok = true;
}
}
} else {
loop = false;
}
}
return ok;
}
bool prune_both(MedianFilterFastHandle_t* const Node) {
bool res = false;
bool res1 = prune_small(Node);
bool res2 = prune_large(Node);
res = res1 || res2;
return res;
}
static MedianFilterBalanceDir_t median_filter_fast_get_balance_dir_even(MedianFilterFastHandle_t* const Node) {
MedianFilterBalanceDir_t balance_dir = MED_FILT_BALANCE_DIR_UNDEF;
/* k is even devided by 2 (4 6 8)*/
if(Node->small_size < Node->large_size) {
balance_dir = MED_FILT_BALANCE_DIR_DOWN;
} else {
if(Node->large_size < Node->small_size) {
balance_dir = MED_FILT_BALANCE_DIR_UP;
} else {
// large_size == small_size
balance_dir = MED_FILT_BALANCE_DIR_NONE;
}
}
return balance_dir;
}
static MedianFilterBalanceDir_t median_filter_fast_get_balance_dir_odd(MedianFilterFastHandle_t* const Node) {
MedianFilterBalanceDir_t balance_dir = MED_FILT_BALANCE_DIR_UNDEF;
/* k is odd (1; 3; 5; 7; 9) */
if(Node->small_size == (Node->large_size + 1)) {
balance_dir = MED_FILT_BALANCE_DIR_NONE;
} else {
if(Node->small_size < Node->large_size) {
balance_dir = MED_FILT_BALANCE_DIR_DOWN;
} else {
// large_size<(small_size+2)
balance_dir = MED_FILT_BALANCE_DIR_UP;
}
}
return balance_dir;
}
static MedianFilterBalanceDir_t median_filter_fast_get_balance_dir(MedianFilterFastHandle_t* const Node) {
MedianFilterBalanceDir_t balance_dir = MED_FILT_BALANCE_DIR_UNDEF;
switch(Node->k_parity) {
case MATH_PARITY_EVEN: {
balance_dir = median_filter_fast_get_balance_dir_even(Node);
} break;
case MATH_PARITY_ODD: {
balance_dir = median_filter_fast_get_balance_dir_odd(Node);
} break;
default: {
balance_dir = MED_FILT_BALANCE_DIR_UNDEF;
} break;
}
return balance_dir;
}
static bool median_filter_fast_shift_down(MedianFilterFastHandle_t* const Node) {
bool res = false;
int32_t top_val = 0;
res = min_heap_pull(&Node->Large, &top_val);
if(res) {
Node->large_size--;
res = max_heap_push(&Node->Small, top_val);
Node->small_size++;
prune_large(Node);
}
return res;
}
static bool median_filter_fast_shift_up(MedianFilterFastHandle_t* const Node) {
bool res = false;
int32_t top_val = 0;
res = max_heap_pull(&Node->Small, &top_val);
if(res) {
res = min_heap_push(&Node->Large, top_val);
Node->small_size--;
Node->large_size++;
prune_small(Node);
}
return res;
}
uint32_t median_filter_fast_get_tolal_size(const MedianFilterFastHandle_t* const Node) {
uint32_t total_size = 0;
uint32_t large_size = min_heap_size(&Node->Large);
uint32_t small_size = max_heap_size(&Node->Small);
total_size = small_size + large_size;
return total_size;
}
static bool median_filter_fast_try_delete_old_top(MedianFilterFastHandle_t* const Node,
const int32_t sample_to_delete,
MedianFilterBalanceDir_t* dir) {
bool res = false;
bool delete_done = false;
MedianFilterBalanceDir_t balance_dir = MED_FILT_BALANCE_DIR_UNDEF;
int32_t large_root = 0;
res = min_heap_peek(&Node->Large, &large_root);
if(res) {
if(large_root == sample_to_delete) {
delete_done = min_heap_delete_root(&Node->Large);
Node->large_size--;
} else {
if(large_root < sample_to_delete) {
balance_dir = MED_FILT_BALANCE_DIR_DOWN;
}
}
}
if(!delete_done) {
int32_t small_root = 0;
res = max_heap_peek(&Node->Small, &small_root);
if(res) {
if(sample_to_delete == small_root) {
delete_done = max_heap_delete_root(&Node->Small);
Node->small_size--;
} else {
if(sample_to_delete < small_root) {
balance_dir = MED_FILT_BALANCE_DIR_UP;
}
}
}
}
*dir = balance_dir;
return delete_done;
}
static bool median_filter_fast_shift(MedianFilterFastHandle_t* const Node, const MedianFilterBalanceDir_t balance_dir) {
bool loop = false;
switch(balance_dir) {
case MED_FILT_BALANCE_DIR_DOWN: {
loop = median_filter_fast_shift_down(Node);
} break;
case MED_FILT_BALANCE_DIR_UP: {
loop = median_filter_fast_shift_up(Node);
} break;
case MED_FILT_BALANCE_DIR_NONE: {
loop = false;
} break;
default:
loop = false;
break;
}
return loop;
}
bool median_filter_fast_delete_old_force_ll(MedianFilterFastHandle_t* const Node, const int32_t old_elememt) {
bool res = false;
bool loop = true;
while(loop) {
MedianFilterBalanceDir_t seek_dir = MED_FILT_BALANCE_DIR_UNDEF;
res = median_filter_fast_try_delete_old_top(Node, old_elememt, &seek_dir);
if(res) {
loop = false;
} else {
res = median_filter_fast_shift(Node, seek_dir);
}
}
return res;
}
static bool median_filter_fast_remove_num_ll(MedianFilterFastHandle_t* const Node, const int32_t old_elememt) {
bool res = false;
bool delete_done = false;
res = hash_table_s8_push_ll(&Node->ToDelete, (int8_t)old_elememt);
int32_t small_root = 0;
res = max_heap_peek(&Node->Small, &small_root);
if(res) {
if(old_elememt <= small_root) {
Node->small_size--;
delete_done = true;
if(old_elememt == small_root) {
prune_small(Node);
}
}
}
if(false == delete_done) {
int32_t large_root = 0;
res = min_heap_peek(&Node->Large, &large_root);
if(res) {
if(large_root <= old_elememt) {
Node->large_size--;
delete_done = true;
if(old_elememt == large_root) {
prune_large(Node);
}
}
}
}
return delete_done;
}
static bool median_filter_fast_rebalance_ll(MedianFilterFastHandle_t* const Node) {
bool res = false;
bool loop = true;
while(loop) {
MedianFilterBalanceDir_t balance_dir = MED_FILT_BALANCE_DIR_UNDEF;
balance_dir = median_filter_fast_get_balance_dir(Node);
loop = median_filter_fast_shift(Node, balance_dir);
}
res = true;
return res;
}
static float median_filter_fast_calc_median(const MedianFilterFastHandle_t* const Node) {
float median = 0.0;
bool res = false;
int32_t small_top = 0;
int32_t large_top = 0;
switch(Node->k_parity) {
case MATH_PARITY_ODD: {
/* not devided by 2*/
res = max_heap_peek(&Node->Small, &small_top);
if(res) {
median = (float)small_top;
}
} break;
case MATH_PARITY_EVEN: {
/*devided by 2*/
res = min_heap_peek(&Node->Large, &large_top);
if(res) {
res = max_heap_peek(&Node->Small, &small_top);
if(res) {
median = ((float)(small_top + large_top)) / 2.0;
}
}
} break;
default:
break;
}
return median;
}
static bool median_filter_fast_detete_num_ll(MedianFilterFastHandle_t* Node , int32_t del){
bool res=false;
res = median_filter_fast_remove_num_ll(Node, del);
return res;
}
static bool median_filter_flush_dir(MedianFilterFastHandle_t* Node, MedianFilterBalanceDir_t dir) {
bool res = false;
bool loop = true;
while (loop) {
loop = median_filter_fast_shift(Node, dir);
if (loop) {
res = prune_both(Node);
}
}
return res;
}
static bool median_filter_flush(MedianFilterFastHandle_t* Node){
bool res = false;
uint32_t tolal_size = median_filter_fast_get_tolal_size(Node);
if ((Node->size*2) <= tolal_size) {
bool res1 = median_filter_flush_dir(Node, MED_FILT_BALANCE_DIR_DOWN);
bool res2 = median_filter_flush_dir(Node, MED_FILT_BALANCE_DIR_UP);
res = res1 || res2;
Node->flush_cnt++;
}
return res;
}
static bool median_filter_fast_push_small(MedianFilterFastHandle_t* Node, const int32_t x) {
bool res = false;
res = max_heap_push(&Node->Small, x);
Node->small_size++;
return res;
}
static bool median_filter_fast_push_large(MedianFilterFastHandle_t* Node, const int32_t x) {
bool res = false;
res = min_heap_push(&Node->Large, x);
Node->large_size++;
return res;
}
static bool median_filter_fast_insert_sample(MedianFilterFastHandle_t* Node, const int32_t x) {
bool res = false;
bool is_insert = false;
int32_t small_root = 0;
res = max_heap_peek(&Node->Small, &small_root);
if (res) {
if (x<=small_root ) {
is_insert = median_filter_fast_push_small(Node, x);
}
}
if (!is_insert) {
int32_t large_root = 0;
res = min_heap_peek(&Node->Large, &large_root);
if (res) {
if (large_root <= x) {
is_insert = median_filter_fast_push_large(Node, x);
}
}
}
if (!is_insert) {
is_insert = median_filter_fast_push_small(Node, x);
}
return res;
}
bool median_filter_fast_proc_in_out(uint8_t num, const int32_t x, float* const out_val) {
bool res = false;
MedianFilterFastHandle_t* Node = MedianFilterFastGetNode(num);
if(Node) {
res = median_filter_fast_insert_sample(Node,x);
bool pull_res = circular_buffer_push_pull(&Node->SlidingWindow, x, &Node->old_elememt);
if(pull_res) {
res = median_filter_fast_detete_num_ll(Node, Node->old_elememt);
median_filter_fast_rebalance_ll(Node);
}
Node->proc_cnt++;
*out_val = median_filter_fast_calc_median(Node);
median_filter_flush(Node);
}
return res;
}
static bool median_filter_fast_init_window(MedianFilterFastHandle_t* const Node) {
bool res = false;
uint32_t i = 0;
for(i = 0; i < Node->size; i++) {
float out_val = 0.0;
res = median_filter_fast_proc_in_out(Node->num, 0, &out_val);
}
return res;
}
static bool median_filter_fast_depensencies(const MedianFilterFastConfig_t* const Config,
MedianFilterFastHandle_t* const Node) {
bool res = false;
BinHeapConfig_t ConfigLarge = {0};
res = median_filter_form_min_heap_config(Config, &ConfigLarge);
res = min_heap_init_one_ll(&ConfigLarge, &Node->Large) && res;
BinHeapConfig_t ConfigSmall = {0};
res = median_filter_form_max_heap_config(Config, &ConfigSmall);
res = max_heap_init_one_ll(&ConfigSmall, &Node->Small) && res;
HashTableS8Config_t ConfigToDel = {0};
res = median_filter_form_hash_table_s8_config(Config, &ConfigToDel);
res = hash_table_s8_init_one_ll(&ConfigToDel, &Node->ToDelete) && res;
res = circular_buffer_dword_init(&Node->SlidingWindow, Config->x, Config->size + 1) && res;
if(res) {
Node->init = true;
}
return res;
}
bool median_filter_fast_init_one(uint8_t num) {
bool res = false;
const MedianFilterFastConfig_t* Config = MedianFilterFastGetConfig(num);
if(Config) {
res = MedianFilterFastIsValidConfig(Config);
if(res) {
MedianFilterFastHandle_t* Node = MedianFilterFastGetNode(num);
if(Node) {
res = median_filter_fast_init_common(Config, Node);
Node->small_size = 0;
Node->large_size = 0;
Node->k_parity = math_calc_parity(Node->size);
res = median_filter_fast_depensencies(Config, Node);
res = median_filter_fast_init_window(Node);
} else {
res = false;
}
} else {
res = false;
}
} else {
res = false;
}
return res;
}
COMPONENT_INIT_PATTERT(MEDIAN_FILTER_FAST, MEDIAN_FILTER_FAST, median_filter_fast)
Отладка
Тестировал этот медианный фильтр я другим медианным фильтром, который реализован просто на основе сортировки элементов в окне. Вот буквально пара итераций из жизни медианного фильтра с окном равным 6 семплов
69.227,2322 N,[MedianFilterFast]
Proc,in:-640,N:4,K:6,Ssz:3,Lsz:3,Small:[;-571;-789;-645;-875;-839;-967;-848;-919;-989],Large:[;-231;0;261;536;589;663;449;550;882;986;840;937;790;994;839;648],ToDel:[648,663,-989,-967,-919,790,-875,839,840,-848,261,882,-789,937,986,994,449,550,589,],Procnt:73,CbX:3/7,Init:On,
69.275,2323 D,[MedianFilterFast] x:-640,small.root:-571
69.289,2324 D,[MedianFilterFast] Push,Small:x:-640
69.298,2325 N,[MedianFilterFast] AfterPush,N:4,K:6,Ssz:4,Lsz:3,Small:[;-571;-640;-645;-875;-789;-967;-848;-919;-989;-839],Large:[;-231;0;261;536;589;663;449;550;882;986;840;937;790;994;839;648],ToDel:[648,663,-989,-967,-919,790,-875,839,840,-848,261,882,-789,937,986,994,449,550,589,],Procnt:73,CbX:3/7,Init:On,
69.329,2326 D,[MedianFilterFast] ElementToDel:0
69.334,2327 D,[MedianFilterFast] pushToDel,Ok
69.342,2328 D,[MedianFilterFast] small_root:-571
69.347,2329 D,[MedianFilterFast] large_root:-231
69.359,2330 D,[MedianFilterFast] Rebalane,N:4,K:6,Ssz:4,Lsz:2,Small:[;-571;-640;-645;-875;-789;-967;-848;-919;-989;-839],Large:[;-231;0;261;536;589;663;449;550;882;986;840;937;790;994;839;648],ToDel:[0,648,663,-989,-967,-919,790,-875,839,840,-848,261,882,-789,937,986,994,449,550,589,],Procnt:73,CbX:4/7,Init:On,
69.395,2331 D,[MedianFilterFast] BalanceDir:Up
69.401,2332 D,[MedianFilterFast] MoveUp:-571
69.408,2333 D,[MedianFilterFast] prune_small
69.421,2334 D,[MedianFilterFast] root_small:-640
69.435,2335 N,[MedianFilterFast] Rebalane,N:4,K:6,Ssz:3,Lsz:3,Small:[;-640;-789;-645;-875;-839;-967;-848;-919;-989],Large:[;-571;-231;261;0;589;663;449;536;882;986;840;937;790;994;839;648;550],ToDel:[0,648,663,-989,-967,-919,790,-875,839,840,-848,261,882,-789,937,986,994,449,550,589,],Procnt:73,CbX:4/7,Init:On,
69.462,2336 D,[MedianFilterFast] prune
69.468,2337 D,[MedianFilterFast] prune_small
69.473,2338 D,[MedianFilterFast] root_small:-640
69.493,2339 D,[MedianFilterFast] prune_large
69.499,2340 D,[MedianFilterFast] root_large:-571
69.514,2341 D,[MedianFilterFast] BalanceDir:None
69.520,2342 N,[MedianFilterFast] Rebalane,N:4,K:6,Ssz:3,Lsz:3,Small:[;-640;-789;-645;-875;-839;-967;-848;-919;-989],Large:[;-571;-231;261;0;589;663;449;536;882;986;840;937;790;994;839;648;550],ToDel:[0,648,663,-989,-967,-919,790,-875,839,840,-848,261,882,-789,937,986,994,449,550,589,],Procnt:73,CbX:4/7,Init:On,
69.552,2343 D,[MedianFilterFast] prune
69.556,2344 D,[MedianFilterFast] prune_small
69.566,2345 D,[MedianFilterFast] root_small:-640
69.581,2346 D,[MedianFilterFast] prune_large
69.586,2347 D,[MedianFilterFast] root_large:-571
69.604,2348 D,[MedianFilterFast] ReBalance,Ok
69.611,2349 N,[MedianFilterFast] N:4,K:6,Ssz:3,Lsz:3,Small:[;-640;-789;-645;-875;-839;-967;-848;-919;-989],Large:[;-571;-231;261;0;589;663;449;536;882;986;840;937;790;994;839;648;550],ToDel:[0,648,663,-989,-967,-919,790,-875,839,840,-848,261,882,-789,937,986,994,449,550,589,],Procnt:74,CbX:4/7,Init:On,
69.641,2350 D,[MedianFilterFast] Even
69.653,2351 D,[MedianFilterFast] PeakMin,Ok
69.658,2352 D,[MedianFilterFast] PeakMax,Ok
69.666,2353 D,[MedianFilterFast] SmallRoot:-640,LargeRoot:-571,median:-605.5
69.675,2354 D,[MedianFilterFast] Proc,in:-640->out:-605.5
69.689,2355 D,[MedianFilterFast] i:67,Exp:-605.5,Real:-605.5:
69.696,2356 N,[MedianFilterFast]
Proc,in:930,N:4,K:6,Ssz:3,Lsz:3,Small:[;-640;-789;-645;-875;-839;-967;-848;-919;-989],Large:[;-571;-231;261;0;589;663;449;536;882;986;840;937;790;994;839;648;550],ToDel:[0,648,663,-989,-967,-919,790,-875,839,840,-848,261,882,-789,937,986,994,449,550,589,],Procnt:74,CbX:4/7,Init:On,
69.743,2357 D,[MedianFilterFast] x:930,small.root:-640
69.755,2358 D,[MedianFilterFast] Push,Large:x:930
69.765,2359 N,[MedianFilterFast] AfterPush,N:4,K:6,Ssz:3,Lsz:4,Small:[;-640;-789;-645;-875;-839;-967;-848;-919;-989],Large:[;-571;-231;261;0;589;663;449;536;882;986;840;937;790;994;839;648;550;930],ToDel:[0,648,663,-989,-967,-919,790,-875,839,840,-848,261,882,-789,937,986,994,449,550,589,],Procnt:74,CbX:4/7,Init:On,
69.795,2360 D,[MedianFilterFast] ElementToDel:-645
69.800,2361 D,[MedianFilterFast] pushToDel,Ok
69.810,2362 D,[MedianFilterFast] small_root:-640
69.816,2363 D,[MedianFilterFast] Rebalane,N:4,K:6,Ssz:2,Lsz:4,Small:[;-640;-789;-645;-875;-839;-967;-848;-919;-989],Large:[;-571;-231;261;0;589;663;449;536;882;986;840;937;790;994;839;648;550;930],ToDel:[0,648,663,-989,-967,-919,790,-875,839,840,-848,261,882,-789,937,986,994,449,-645,550,589,],Procnt:74,CbX:5/7,Init:On,
69.850,2364 D,[MedianFilterFast] BalanceDir:Down
69.856,2365 D,[MedianFilterFast] min_heap_pull,Ok
69.871,2366 D,[MedianFilterFast] MoveDown:-571
69.876,2367 D,[MedianFilterFast] max_heap_push,Ok
69.884,2368 D,[MedianFilterFast] PushSmall:-571 Ok
69.890,2369 D,[MedianFilterFast] prune_large
69.904,2370 D,[MedianFilterFast] root_large:-231
69.919,2371 N,[MedianFilterFast] Rebalane,N:4,K:6,Ssz:3,Lsz:3,Small:[;-571;-640;-645;-875;-789;-967;-848;-919;-989;-839],Large:[;-231;0;261;536;589;663;449;550;882;986;840;937;790;994;839;648;930],ToDel:[0,648,663,-989,-967,-919,790,-875,839,840,-848,261,882,-789,937,986,994,449,-645,550,589,],Procnt:74,CbX:5/7,Init:On,
69.946,2372 D,[MedianFilterFast] prune
69.952,2373 D,[MedianFilterFast] prune_small
69.957,2374 D,[MedianFilterFast] root_small:-571
69.978,2375 D,[MedianFilterFast] prune_large
69.985,2376 D,[MedianFilterFast] root_large:-231
69.999,2377 D,[MedianFilterFast] BalanceDir:None
70.006,2378 N,[MedianFilterFast] Rebalane,N:4,K:6,Ssz:3,Lsz:3,Small:[;-571;-640;-645;-875;-789;-967;-848;-919;-989;-839],Large:[;-231;0;261;536;589;663;449;550;882;986;840;937;790;994;839;648;930],ToDel:[0,648,663,-989,-967,-919,790,-875,839,840,-848,261,882,-789,937,986,994,449,-645,550,589,],Procnt:74,CbX:5/7,Init:On,
70.037,2379 D,[MedianFilterFast] prune
70.044,2380 D,[MedianFilterFast] prune_small
70.059,2381 D,[MedianFilterFast] root_small:-571
70.079,2382 D,[MedianFilterFast] prune_large
70.084,2383 D,[MedianFilterFast] root_large:-231
70.100,2384 D,[MedianFilterFast] ReBalance,Ok
70.108,2385 N,[MedianFilterFast] N:4,K:6,Ssz:3,Lsz:3,Small:[;-571;-640;-645;-875;-789;-967;-848;-919;-989;-839],Large:[;-231;0;261;536;589;663;449;550;882;986;840;937;790;994;839;648;930],ToDel:[0,648,663,-989,-967,-919,790,-875,839,840,-848,261,882,-789,937,986,994,449,-645,550,589,],Procnt:75,CbX:5/7,Init:On,
70.143,2386 D,[MedianFilterFast] Even
70.154,2387 D,[MedianFilterFast] PeakMin,Ok
70.159,2388 D,[MedianFilterFast] PeakMax,Ok
70.170,2389 D,[MedianFilterFast] SmallRoot:-571,LargeRoot:-231,median:-401
70.177,2390 D,[MedianFilterFast] Proc,in:930->out:-401
70.187,2391 D,[MedianFilterFast] i:68,Exp:-401,Real:-401:В целом, реализация на кучах работает, но есть целый ворох проблем
Недостатки реализации медианного фильтра на двух кучах
1--Не каждая итерация фильтра приводит к удалению старого элемента.
2--На Си алгоритм получился весьма сложным. Его трудно отлаживать, настраивать и ремонтировать.
3--Требуются внешние вспомогательные программные компоненты: бинарные кучи и hash таблицы и сдвиговый регистр.
4--Если фильтровать случайный поток чисел, то размер максимальной кучи и размер минимальной кучи бесконтрольно увеличивается. Даже внутри балансировки prune не улучшает ситуацию. Со временем любая статическая память переполняется и фильтр просто заклинивает. К примеру, если же делать flush каждый раз, когда память куч превысит, условно, 2 размера окна (2k), то фильтр станет работать аж в 1106 раз медленнее, чем реализация на сортировке qsort.
5--Если фильтровать случайный поток чисел, то размер хэш-таблицы тоже бесконтрольно увеличивается.
6--Если фильтровать случайный поток 32-битных чисел, то в хэш-таблице происходят постоянные коллизии, которые тоже просаживают производительность. Ситуация сильно упрощается, если мы имеем дела с 8-битными семплами. Тогда в качестве ключа можно использовать сами значения семплов и ликвидировать коллизии как таковые в принципе. Цена решения - 1kByte RAM. Однако для чисел 32 бит такая реализация хэш-таблици потребовала бы уже от 4GByte до 16GByte RAM памяти.
7--Сумма элементов в верхней и нижней кучах может оказаться больше, чем размер окна (порядок фильтра). Элемент, который надо удалить может уйти далеко от корней кучи или вовсе оказаться крайним листом в дереве.
Достоинства реализации медианного фильтра на двух кучах
++ Можно использовать медианный фильтр на кучах только, как один большой интеграционный тест для бинарных куч, хэш-таблицы и сдвигового регистра.
Боюсь, что больше плюсов просто нет.
Что можно улучшить?
Если требуется найти не точную медиану, а лишь её приблизительное значение, то можно всё сделать на одном AVL дереве. Добавил, удалил и прочитал получившийся корень AVL дерева.
Итог
Удалось реализовать медианный фильтр на основе двух бинарных куч. Тем не менее я бы не рекомендовал использовать такую реализацию фильтра в микроконтроллерных прошивках. Уж очень она сильно пожирает RAM память. На случайном наборе входных данных кучи могут сильно разрастаться.
Без процедуры flush иногда медианный фильтр на кучах работает медленнее. Иногда быстрее. Иногда и вовсе зависает. Все зависит от того, как генератор случайных чисел сформирует входные семплы.
Реализация МФ на двух кучах не имеет смысла так, как та же qsort на PC реализована с учетом аппаратных возможностей конкретного CISC процессора, чтобы быть максимально производительной.
Если вам удавалось реализовать на чистом Си медианный фильтр, который работает быстрее, чем сортировка qsort-ом окна на каждой итерации, то напишите про это в комментариях
Словарь
Акроним |
Расшифровка |
AVL-дерево |
Дерево имени Адельсона-Вельского Георгия Максимовича и Ландиса Евгения Михайловича. (Georgy Adelson-Velsky and Evgenii Landis) |
RAM |
Random Access Memory |
FIFO |
First in, first out |
ЦОС |
Цифровая обработка сигналов |
ВУЗ |
Высшее учебное заведение |
SPI |
Serial Peripheral Interface |
FIR |
Finite impulse response |
IIR |
Infinite impulse response |
CISC |
Complex instruction set computer |
DSP |
Digital signal processor |
ASIC |
Application-specific integrated circuit |
Ссылки
Название |
URL |
вычисление медианы массива (или произвольной моды) |
http://electronix.ru/forum/index.php?showtopic=114436&hl=median&st=0 |
Визуализация деревьев |
https://www.cs.usfca.edu/~galles/visualization/Algorithms.html |
Функция qsort |
|
Структуры данных: двоичная куча (binary heap)@awolf |
|
Leetcode 480 Sliding Window Median Heaps PriorityQueue Problem Solving DSA Greedy |
|
480. Sliding Window Median - Detailed Explanation |
https://www.designgurus.io/answers/detail/480-sliding-window-median-detailed-explanation |
разбор задачи 480. Sliding Window Median |
|
480. Sliding Window Median |
|
Median filter |
|
480. Sliding Window Median |
|
Медианный фильтр |
https://chipenable.ru/index.php/embedded-programming/203-mediannyy-filtr.html |
@Mihahanya |
|
@AndreyIvanoff |
|
@vadimr |
Вопросы:
1--Существуют ли аппаратные реализации медианного фильтра в ASIC исполнении? В виде отдельных микросхем с SPI-интерфейсом или вовсе с параллельной шиной.
2--Как быть если добавляемый элемент оказался между верхушками куч? Куда его добавлять? В LargeHeap или в SmallHeap?
3--А что если поставить последовательно друг за другом два медианный фильтра длинны 3? Заменят ли они собой один медианный фильтр порядка 6 ?
4--Существуют ли микроконтроллеры с аппаратным блоком сортировки ? Это может пригодится для потоковой реализации медианной фильтрации .
Комментарии (92)
apcs660
07.09.2025 00:54торможу с утра, рано встал но эта задачка напомнила мне старые схемы, типа "интегратор - линия задержки - вычитатель":
cumsum - интегратор
линия задержки - на выход идет cumsum[i-n] где n ширина фильтра
вычитатель вычитает задержанный сигнал cumsum[i-n] из текущего сигнала с интегратора cumsum[i] : window_sum = cumsum[i] - cumsum[i-n]
делитель на n - moving_average[i] = window_sum / n.
остается проверить наличие ближайщих целых значений в буфере по вычисленному среднему из 4го шага (хеш мап значения элемента на массив/лист индекса элемента в буфере использовать - reverse index типа) или просто округлить (в реальной задаче, в этой элемент так понял нужно брать из существующих а не просто округление)? Пойду кофе попью
randomsimplenumber
07.09.2025 00:54Как то сложно. Вставка в отсортированный массив не решает ту же задачу?
wataru
07.09.2025 00:54Решает, но медленнее. Тут делается O(log(W)) операций при сдвиге окна, а для вставки в отсортированный массив надо O(W), где W - размер окна. Если же вместо массива делать какое-нибудь балансированное бинарное дерево, то там константа чуть хуже, потому что эти деревья чуть-чуть медленнее куч.
apcs660
07.09.2025 00:54С утра не сразу проснулся когда статью открыл. Почитал, хороший алгоритм:
Добавление нового значения: O(log N)
Удаление старого значения: O(log N)
Получение медианы: O(1)
Амортизированная стоимость на операцию: O(log N)
Две кучи: 2-4 операции O(log N) + балансировка
Простая сортировка: O(1) для малых N
На 15 элементах почти в 3 раза быстрее чем алгоритм со вставкой в сортированный массив но и требует 3 раза больше памяти (обмен память на скорость).
aamonster
07.09.2025 00:54В сбалансированное дерево тогда уж, чтоб вставка/удаление дорогими не были.
aabzel Автор
07.09.2025 00:54Сбалансированные деревья скорее всего дадут не точную медиану, а значение медианы с погрешностью.
aamonster
07.09.2025 00:54С чего бы? Сбалансированное дерево позволяет получить элемент с заданным номером за O(log N).
С погрешностью – это если на гистограммах делать или приближать гистограмму моментами.
aabzel Автор
07.09.2025 00:54Нужно поставить решающий эксперимент.
aamonster
07.09.2025 00:54Зачем? Контейнер для хранения чисел внутри окна – чёрный ящик. Мы в него на каждом шаге кладём одно число, убираем одно число и извлекаем медиану.
Референсная реализация – отсортировать массив и взять из средней позиции. Легко доказать, что результат будет тот же.
Эксперимент был бы интересен для приближения гистограммы моментами и тому подобных приближённых методов (позволяющих считать за O(1) независимо от размера окна)
Или вы имели в виду эксперимент по скорости?
wataru
07.09.2025 00:54Нет, в сбалансированных деревьях можно найти элемент по его номеру. Правда, для этого надо хранить в каждой вершине размер поддерева и стандартные реализации в большинстве языков так не умеют.
wataru
07.09.2025 00:54Проблема со сбалансированным деревом в том, что там для поиска медианы надо писать его руками и хранить в каждой вершине размер поддерева. Стандартные реализации в большинстве языков так не умеют. Так что реализация с двумя кучами или двумя деревьями получается проще и даже быстрее.
aamonster
07.09.2025 00:54Ну так-то да. Реализация с деревом – просто база, от которой начинаешь думать, уже зная, что уложишься в O(log N).
f45d07
07.09.2025 00:54Вот вам яркий пример. Вы измеряете расстояние UWB трансиверами между радиопередатчиками и вдруг у вас образуются единичные редкие отбросы, которые никак не характеризует измеряемую величину. Или вы опрашиваете емкостную кнопку и от скача питания она тоже может дать редкие, но высокие значения семпла.
С примером графика сырых и отфильтрованных данных получившимся фильтром было бы интереснее)
apcs660
07.09.2025 00:54да, это бритва для среза волосков. А усреднитель это совочек который выброс размажет и несущую исказит
mozg37
07.09.2025 00:54Допустим, окно из 16 элементов. Делаем сдвиговый регистр(линию задержки) на 16 элементов. С каждым новым элементом считаем добавляем к некому числу текущий элемент и вычитаем 16й. Это число есть сумма последних 16 элементов. Добавляйте свою логику для контроля единичных выбросов. Так у меня сделано на fpga. Окно при использовании RAM блоков может быть весьма большим, на одном из проектов 2*1024 элемента. Без всяких куч и прочего.
aamonster
07.09.2025 00:54Скользящее среднее – совсем другой фильтр. Лучше убирает гауссовский шум, хуже удаляет выбросы. Разница очень заметна при обработке изображений: медиана убирает одиночные точки/пятнышки на изображении, оставляя границы резкими, бегущее среднее – размывает границы, а точки полностью не удаляет.
Это не в том смысле, что медианный фильтр лучше – он просто решает другие задачи.
rukhi7
07.09.2025 00:54Неужели так сложно догадаться как с помощью скользящего среднего решить задачу удаления выбросов? Без всяких куч или даже просто массивов?
Это, наверно, из тоже же серии задач про Метод инверсии или методика проектирования «от противного», абсурдной перестановки. или проще?
aabzel Автор
07.09.2025 00:54Да, сложно.
rukhi7
07.09.2025 00:54в ВУЗах обычно про медианный фильтр даже не вспоминают,
Вы же знаете почему про него не вспоминают в ВУЗах? Почему на нем нельзя построить никакой математики цифровой обработки сигналов? Почему с точки зрения ЦОС это изобретение, вообще говоря, не является фильтром?
Хотя в качестве упражнения по программированию вполне себе интересно, пока вы не пытаетесь строить ЦОС на этой основе.
aamonster
07.09.2025 00:54Как человек, 10 лет занимавшийся разработкой софта для анализа изображений, могу ответить: да, сложно.
Конечно, скользящее среднее – это база, если в вашей задаче можно им воспользоваться – ура, мы сэкономили кучу ресурсов. Но оно (как и дешёвые решения на его основе) подходит не всегда. Могут понадобиться оба фильтра – и медиана, и скользящее среднее (возможно, повторённое 2-3 раза, чтобы приблизить к гауссовскому фильтру).
Ещё, кстати, в копилку дешёвых методов – open и close фильтры. Вначале выполняем бегущий минимум (эрозия), потом бегущий максимум (дилатация). Ну или наоборот. Убирает светлые или тёмные пятна соответственно. Для изображений есть усовершенствование – open/close with reconstruction. Но тоже особо не используется отдельно.
Кстати, бегущий максимум/минимум, не зависящий от размера окна – сам по себе интересный алгоритм. Изобрели, если я правильно помню, только в 80-х. Попробуйте до него додуматься сами, не подглядывая – это правда интересно (делается без сложных структур данных).
Alexandroppolus
07.09.2025 00:54бегущий максимум/минимум, не зависящий от размера окна
https://leetcode.com/problems/sliding-window-maximum - этот?
aamonster
07.09.2025 00:54Да, это оно. Конечно, для leetcode задачка так себе, т.к. она больше на знание алгоритма (не думаю, что многие смогут придумать его сами). Интересно, уложится ли в тайминги литкода решение с логарифмической сложностью.
Alexandroppolus
07.09.2025 00:54Там по факту два разных решения, одно из которых легко придумать по мотивам реализации очереди на двух стеках и хранении максимумов для стека (у меня так и вышло)
aamonster
07.09.2025 00:54Не понял. Короче, у вас решение за O(1) работает или за O(log N)? Для O(1) есть хитрость, просто хранить две кучи и обновлять их на каждом шаге – получите O(log N).
Alexandroppolus
07.09.2025 00:54У меня работает за амортизированное O(1) на элемент, то есть за O(N) на всё про всё. Вот, опубликовал. Раз в k итераций собираю стек максимумов по окну справа налево, далее на следующих k итерациях благодаря этому под рукой имеется максимум среди тех элементов, которые покидают окно
aamonster
07.09.2025 00:54Угу, это тот самый алгоритм, который изобрели удивительно поздно, несмотря на его простоту.
Alexandroppolus
07.09.2025 00:54Кстати, если не ошибаюсь, тут можно добиться не амортизированного, а "истинного" O(1) на элемент, если вместо базовой идеи очереди на двух стеках, использовать великую и ужасную очередь на 6 стеках (и соответственно к ним стеки с максимумами). Насколько это будет оптимально, не берусь судить, но факт любопытный
Melirius
07.09.2025 00:54Я б скользящее окно размером порядка фильтра сохранял и удалял последний старый элемент на каждой итерации из соответствующей кучи, добавляя новый.
aabzel Автор
07.09.2025 00:54При работе с бинарной кучей удалять можно только то, что лежит в корне кучи.
Если число закопано, то ничего быстро не сделаешь.Alexandroppolus
07.09.2025 00:54Долго искать элемент внутри кучи. Но если знаем его позицию, то можно быстро удалить (переносим на его место последний элемент, потом просеивание вниз или вверх)
ripatti
07.09.2025 00:54Это в базовой куче. Если завести массив отслеживания позиций элементов в куче (в код кучи вносится модификация - при каждом свапе элементов этот массив корректируем) - то можно вполне быстро.
VT100
07.09.2025 00:54Спасибо. Но... проверьте правописание и пропущенные фрагменты предложений.
aabzel Автор
07.09.2025 00:54Где ещё остались опечатки?
Nemoumbra
07.09.2025 00:54При использовании этого метода элементы которые надо удалить помечаются словом откладываются
Я в этот момент решил, что статья - перевод без тега "перевод".
Mike-M
07.09.2025 00:54Перебалансировка куч выполняется после каждой операции добавления или удаления, чтобы свойство половинного размера двух куч.
Как насчет использования бесплатных или дешевых сервисов проверки текстов, таких как Яндекс Спеллер, Text.ru, Орфограммка? А то обидно получается: элементарные ошибки затрудняют восприятие достойного контента.
VT100
07.09.2025 00:54По оборванным предложениям - ответили @Nemoumbra и @Mike-M.
"но не удалены физически из куч
ей."
"В случае сильного разрастания кучейможно вручную удалить ненужные более элементы"
"Реальные же размеры кучейне должны участвовать, как управляющие данные для алгоритма."
Sulerad
07.09.2025 00:54Строго говоря, алгоритм довольно несложно доработать так, чтобы устаревшие элементы удалялись из куч в тот же момент, когда они вышли из окна. Но для этого потребуется заменить кучу на ваше любимое дерево поиска — тогда удаление станет возможно за O(log N) на каждый новый элемент, где N — размер окна. А реализация какого-нибудь КЧД это уже далековато от "несложно" =)
Однако же альтернативный вариант можно подглядеть здесь: https://stackoverflow.com/a/10696252 — там предлагается аналогичное решение, но с заменой кучи/дерева на список с пропусками, который обладает примерно теми же преимуществами и без необходимости реализовывать балансировку. Но в худшем случае он потребует O(N), что становится практически аналогично quicksort.
В целом, задача как ни крути сводится к поддержанию сортированного списка и быстрее чем O(log N) её решить в общем случае вряд ли выйдет. Учитывая что N (размер окна) довольно мало, то боюсь что любой выигрыш от сложных структур данных запросто будет нивелирован прямолинейностью quicksort.
aabzel Автор
07.09.2025 00:54Если просить не точно медиану, а ее приблизительное значение, то можно всё сделать на одном AVL дереве.
aamonster
07.09.2025 00:54Из него не обязательно извлекать корень. Например, держите в узле количество его подэлементов. Тогда можете сразу определить, искомый элемент в корне, левом или правом поддереве, и рекурсивно повторить.
LaRN
07.09.2025 00:54Если числа целые и ограничены например 00-FF, то ту наверное обычный массив который в ячейке с индексом равным этому числу храниться количество раз, когда это число встретилось в потоке может быть проще и быстрее для вычисления медианы и по памяти не очень затратно.
aabzel Автор
07.09.2025 00:54У меня так хеш-таблица работает.
Но как из хэш-таблицы посчитать медиану?
aamonster
07.09.2025 00:54Не хэш-таблицы, гистограммы. Просто считать частичную сумму, пока не доберётесь до половины длины окна. Ну, естественно, можно оптимизировать, храня и обновляя некоторые из частичных сумм.
Имеет смысл только для очень больших размеров окна, и то не факт.
Jijiki
07.09.2025 00:54интересно. я не понял, а зачем сортировку делать. Gap_buffer мне это вспомнилось
есть еще емуляция буферного окна, это когда просто есть мгновенный доступ к нужным позициям-курсор в окне, тоесть установить курсор(и проделать фильтр)
тогда получается есть 1 heap-array(cм. ниже) для всех семплов например песня, если надо его можно поделить на части более маленькие, тогда окно будет выборка тоесть get[i] или нужного количества семплов из общего количества семплов.
если это перводить на рельсы сленга gap-buffer, тогда мы получаем на выходе следующее:
есть абстракция array-heap кусочек - строка
есть абстракция buffer-heaps которая знает сколько кусочков и где они, в ней есть курсор и определение строки(в примере это размер array-heap)
тогда МФ будет просто работа с кусочками[n] из buffer(по типо как в окно программы помещается 40 строк эти строки мы рисуем и управляем курсором и позицией)
тоесть, вычисляем общее количество семплов, потом делим, добавляем семплы в array-heap-это 1 кусок, а буфер хранит эти строки так называемые и даёт доступ к адресам этих строк-семпловых
typedef struct { char *data; size_t length; size_t capacity; } String; typedef struct { String **line; //pointer to lines size_t nlines; //number of lines size_t capacity; // size_t currLine; size_t totalSizeChars; int stateFlag;//0 scratch,1 openFile } Buffer;
aamonster
07.09.2025 00:54Удаление или замена элемента в куче – O(log N), можно не хранить удалённые элементы.
aabzel Автор
07.09.2025 00:54Хорошо. А как реализация этого на си выглядит?
aamonster
07.09.2025 00:54Уффф... Спросите у чатгпт, я на память не воспроизведу. Вкратце – заменяем элемент (для удаления – на последний из кучи) и восстанавливаем порядок в куче операциями siftUp/siftDown, там не надо всю кучу перелопачивать.
aabzel Автор
07.09.2025 00:54static bool min_heap_make_down(BinHeapHandle_t* Node, int32_t index) { bool res = true; while(index < Node->size) { int32_t left_child, right_child, smallest; left_child = 2 * index + 1; right_child = 2 * index + 2; smallest = index; if(left_child < Node->size && Node->array[left_child] < Node->array[smallest]) { smallest = left_child; } if(right_child < Node->size && Node->array[right_child] < Node->array[smallest]) { smallest = right_child; } if(smallest == index) { break; } swap_i32(&Node->array[index], &Node->array[smallest]); index = smallest; } return res; } static int32_t min_heap_check_node(BinHeapHandle_t* const Node, const int32_t value, uint32_t my_index) { int32_t index = -1; if (my_index < Node->size) { if (value == Node->array[my_index]) { index = my_index; } else { if (Node->array[my_index] < value) { uint32_t left = 2 * my_index + 1; index = min_heap_check_node(Node, value, left); if (-1 == index) { uint32_t right = 2 * my_index + 2; index = min_heap_check_node(Node, value, right); } } else { index = -1; } } } return index; } bool min_heap_delete(BinHeapHandle_t* const Node, const int32_t value){ bool res = false; int32_t del_index = min_heap_check_node(Node, value, 0); if(0 <= del_index) { swap_i32(&Node->array[del_index], &Node->array[Node->size-1]); Node->size--; res = min_heap_make_down(Node, del_index); } return res; }
wataru
07.09.2025 00:54В такой же хеш таблице надо для каждого элемента в куче хранить его позицию (при просеивании элементов не забывать обновлять позиции у обоих). Потом при удалении вы уже знаете, где элемент лежит - меняете его с последним в куче и оттуда удаляете. А то, что на его место поставили надо просеять или вниз или вверх.
krotos139
07.09.2025 00:54На мой взгляд сделать можно сильно проще.
1) Сделайте 1 сложносвязанный список размером с размер окна и связывающий каждый элемент в 2 списка. Далее эти элементы свяжите в 1 список отсортированный по значению величини в ноде и во 2 список по порядку добавления.
Далее при добавлении нового элемента в кучу сделайте следующее: возьмите последний элемент из списка по добавлению и удалите его(поменяв указатели которые формируют список по порядку и список по значению). Далее в этом элементе поменяйте значение на новое и вставьте этот элемент в начало списка по порядку и в список отсортированный по значению в нужное место.
Медианный фильтр готов.. Не нужно ни внешние библиотеки, ни ребалансировки ни хэш таблицы и потребление памяти фиксированное, какие бы значения туда ни попадали.
Для получения среднего значения возьмите размер окна поделите на 2 и возьмите этот элемент из списка отсортированного по величине. Если вам нужны значения для серидины верхней и нижней кучи из вашего алгоритма - просто возьмите элементы 1/4 и 3/4 размера окна.
rukhi7
07.09.2025 00:54настоящее решение никому не интересно на хабре :). Как это без
библиотеки; ребалансировки; хэш таблицы; просеивании элементов ; поОкать - О(амен) О(м); ... ?
тут сразу становится не о чем поговорить :) !
wataru
07.09.2025 00:54А потом пользователи жалуются, что все тормозит. А просто программисту было лень "Окать" и он впендюривал наивные очень медленные решения повсюду. О-большое - это база программирования. Это знать надо.
rukhi7
07.09.2025 00:54когда нам 25+ лет назад надо было HD видео играть на пентиумах нас про всякие О-разноразмерных никто не спрашивал, всех интересовали фреймы в секунду. А теперь все учатся правильно писать О-большое, про фреймы в секунду никто не вспоминает! По моему это просто такой вид моды, мода пройдет, я думаю. Хотя теоретическая математика тоже важна, конечно, но она достаточно ограничена.
aabzel Автор
07.09.2025 00:54Я не понял. Есть ли возможность нарисовать эту структуру данных?
krotos139
07.09.2025 00:54
Вот нарисовал такой сложносвязанный список из 5 элементов. Он отображает какая будет структура связей при получении данных 35, 8, 17, 500, 21.
Если обходить используя указатели Head1, Tail1 - то это будет список по порядку добавления. Если его обходить используя указатели Head2, Tail2 - то это будет список отсортированный по значениям (8, 17, 21, 35, 500).
wataru
07.09.2025 00:54вставьте этот элемент ... в список отсортированный по значению в нужное место
Вот тут проблема. А как найти это место? Бинарный поиск в списке не работает. Просто пройтись по списку - это линейная сложность O(размер окна). Когда как решение в статье рабоатет за O(log(размер окна)), что на порядки быстрее.
Все эти хитрые структуры данных исключительно чтобы быстрее работало.
Если вам неважна скорость, то зачем извращаться со списками - просто берете последние W элементов из цикличиского буфера и сортируете каждый раз при запросе медианы. Два массива и sort.
rukhi7
07.09.2025 00:54Бинарный поиск в списке не работает.
это прям не решаемая проблема сделать чтобы бинарный поиск работал на упорядоченном списке?
wataru
07.09.2025 00:54Вот тут явно видно нежелание задуматься о сложности алгоритма.
Бинарный поиск основан на том, что можно быстро взять элемент по его индексу, что в связных списках не работает.
Бинарный поиск смотрит на элемент в середине и решает в какую сторону дальше идти: влево или вправо. Но как посмотреть на средний элемент в списке? Это не массив, где можно по индексу в одну арифметическую операцию получить адрес элемента. Тут придется пройтись по всему списку до середины. И потом опять также в нужной половине. Потом опять в четвертине и т. д. Каждый раз надо проходится до середины.
В итоге такой бинарный поиск выполнит n/2 + n/4 + n/8 +... = n-1 операций, т.е. оказывается ничем не лучше простого наивного линейного поиска. Только сложнее код и всякие лишние операции: линейный поиск может остановится не дойдя до середины, если наткнется на нужное место, бинарный же - всегда выполнит эти самые n-1 операций.
rukhi7
07.09.2025 00:54быстро взять элемент по его индексу, что в связных списках не работает.
он там выразился не совсем корректно, а вы не стали вникать, потому что просто не верите в более простое решение (без всяких умных слов), мне кажется.
Сделайте 1 сложносвязанный список размером с размер окна
поскольку размер списка жестко задан - это массив, а не список.
Верить полезно - английские ученые доказали :) !
wataru
07.09.2025 00:54> поскольку размер списка жестко задан - это массив, а не список.
Ну подумайте чуть глубже и посмотрите на весь алгоритм. Это все еще список и элементы хоть и хранятся в масиве, они там не упорядочены по возрастанию. Третий по порядку элемент может быть хоть в начале, хоть в конце массива и не проходясь по ссылкам вы его никак не найдете.Если же вы сделаете хранение элементов в массиве по возрастанию, то зачем там вообще сложно-связный список, если все указатели будут на следующий элемент в массиве? А главная проблема, вставка в такой массив все равно будет линейной, ибо надо будет почти все элементы сдвигать, чтобы новый вставить в середину.
Так что никак вы бинарный поиск со списками так просто не скрестите, вы или убиваете скорость бинарного поиска, или скорость изменения списка. В итоге получается более сложное и медленное решение чем с простым массивом.
Вообще, если в эту сторону думать, то в итоге вы придете все к тем же балансированным деревьям поиска или скип-листам, которых и пытались избежать.
rukhi7
07.09.2025 00:54Ну подумайте чуть глубже и посмотрите на весь алгоритм.
дело в том что я смотрю не на алгоритм, а на задачу и думаю как ее решить самым простым способом.
Сначала надо доказать что выбранный алгоритм является самым простым.
wataru
07.09.2025 00:54Тут всегда есть выбор между эффективностью и простотой: самый простой способ обычно самый медленный, а самый быстрый - сложнее. Иногда есть что-то среднее балансирующее между простотой и скоростью.
Так вот, предложенный @krotos139 алгоритм и медленный и сложный. Он все такой же линейный как наивный способ вставки в сортированный массив, но при этом усложнен двусвязными списками.
Сначала надо доказать что выбранный алгоритм является самым простым.
Самый простой еще медленее: Вы просто берете последние W элементов и сортируете их каждый раз. Куда уж проще? Хотя, возможно, кому-то более простым методом покажется метод за O(W^2) - берете каждый элемент из окна и считаете, сколько там меньших или равных ему, и берете минимальный с нужным количеством. Тут нет сложной концепции сортировки, хотя если использовать библиотечную функцию, этот метод будет содержать больше кода и переменных.
Но вообще, формально доказать, что какой-то метод самый простой очень сложно, ибо простота - вещь сложно формализуемая. И вообще, оценки снизу доказывать сильно сложнее оценок сверху.
Jijiki
07.09.2025 00:54сдвигать? всмысле memcpy есть же, выделенное окно или та структура где будут фильтры можно представлять в виде блока данных просто вставили участок байт же нет?
просто вы сейчас описали тривиальный момент вставка определенного количества байт в центр строки
потом если смотреть по вики что такое МФ это ну картинка самое простое, а картинка это строки хотим мы того или нет
тоесть семпл это елемент строки, а строка это последовательность семплов и такая строка может считаться в байтах
тоесть самое ближайшее это uint32 codepoints строки
19 UTF-32 code units: [ 0x20 0x20 0x2f 0x2f 0x70 0x72 0x69 0x6e 0x74 0x5f 0x73 0x74 0x72 0x28 0x26 0x73 0x29 0x3b 0xa ] 25 UTF-32 code units: [ 0x20 0x20 0x2f 0x2f 0x70 0x72 0x69 0x6e 0x74 0x66 0x28 0x22 0x25 0x7a 0x75 0x5c 0x6e 0x22 0x2c 0x6d 0x61 0x78 0x29 0x3b 0xa ] 19 UTF-32 code units: [ 0x20 0x20 0x53 0x74 0x72 0x69 0x6e 0x67 0x5f 0x66 0x72 0x65 0x65 0x28 0x26 0x73 0x29 0x3b 0xa ] 25 UTF-32 code units: [ 0x20 0x20 0x2f 0x2f 0x70 0x72 0x69 0x6e 0x74 0x66 0x28 0x22 0x25 0x7a 0x75 0x5c 0x6e 0x22 0x2c 0x6d 0x61 0x78 0x29 0x3b 0xa ] 12 UTF-32 code units: [ 0x20 0x20 0x72 0x65 0x74 0x75 0x72 0x6e 0x20 0x30 0x3b 0xa ] 2 UTF-32 code units: [ 0x7d 0xa ]тоесть не окно идёт по семплам а курсор и в месте курсора подставляется медиана тоесть в месте курсора фильтрация
тогда получается из потока семплы попадают в буффер, и там курсор
in [buffer[cursor]] out
если произошел out скинуть курсор в начало и получить новый блок данных из входа
wataru
07.09.2025 00:54сдвигать? всмысле memcpy есть же
И за какую сложность, по-вашему, работает memcpy?
Cтроки тут вообще нипричем. Дальше просто несвязный бред написан, как у вас обычно и просиходит.
Jijiki
07.09.2025 00:54просто выделить 3*1024, и через функцию в буфере адресно можно конструировать выводной поток дописывать в конец поблочно по 3*1024, помимо memcpy есть прямой доступ к елементам 3х3(например ширина курсора внутри буфера)
wataru
07.09.2025 00:54Ничего непонятно, но очень интересно (нет).
Jijiki
07.09.2025 00:54да, потомучто у дерева есть недостатки в буферной обработке(лишние итерации, которые можно опустить прямым доступом)
а деревом пользоваться, а обращаться к лепесткам как по массиву, зачем если есть итак массив
пришло 3 строки в буфер курсор 3х3 считаем пишем чем хуже
если нет строк последовательность можно поделить на 3 и использовать offset
или не делить и идти построчным курсором до конца последовательности(вобщем деление последовательности даст улучшение, это все как с текстом, мы можем иметь 1 большую строку на гигабайт, а можем делить на строки )
wataru
07.09.2025 00:54Остановитесь, глубоко вдохните и выдохните. Осознайте, что другие люди не имеют доступа к вашим мыслям, у нас не распределенный разум-рой и нет телепатии. Перечитайте статью и ветку комментариев, потом составьте на бумаге всю цепочку ассоциаций от этого общего контекста к вашим мыслям и укажите главные из них в комментарии. Тогда и только тогда будут шансы что вас поймут.
Каждый раз перед тем как запостить комментарий, убедитесь, что он опрерирует теми же коцепциями, что и тот, на который вы отвечаете. А если нет, то вы должны составить логическую цепочку от уже упомянутых концепций к вашим.Вот memcpy вы даже к цитате привязали и на эту часть я ответить что-то разумное смог. Остальное как с неба свалилось: какие-то строки, курсоры, буфера.
Jijiki
07.09.2025 00:54я выдыхая осознал когда вы заминусили существование использования в таком роде структур memcpy впринципе всё понял
а что не так, человек пишет про последовательность байт, эту последовательность при получении надо поделить на 3 равных, знать эти адреса, указать ширину курсора(3х3) и вуаля есть все нужные адреса, если надо записать в центр 1/3 есть memcpy блочный
wataru
07.09.2025 00:54а что не так, человек пишет про последовательность байт, эту последовательность при получении надо поделить на 3 равных
Какой человек? Откуда взялось число 3?
Jijiki
07.09.2025 00:54ну семпл это же мерная величина измеряется в байтах
потомучто в 1 проход дороже чем по трём полосам на своих адресах
в 1 полосе 1 адрес и величина блока больше на 1/3, если поделить 1 большой на 3 равных по семплам выравнивая это будет логичнее, тогда надо позаботиться чтобы на входе было равное поступление семплов за такт или отрезок времени
wataru
07.09.2025 00:54Хорошо, семпл измеряется в байтах. Вообще все в компьютере в байтах представлено. Но при чем тут это, чем это помогает в поиске медианы?
Какой человек писал:
человек пишет про последовательность байт, эту последовательность при получении надо поделить на 3 равных
Откуда взялось 3, почему равные части?
Вы опять выдали какую-то свою мысль никак ее не привязав к упомянутому раньше и я опять не понял, что вы хотели сказать.
Jijiki
07.09.2025 00:54ну МФ в последовательности семплов так? последовательность семплов это и строка семплов, а значит если настроить равное поступление семплов например на 1/3 больше или просто чтобы равное поступало поделив например на 3, и распределив эти куски на свои адреса алгоритм попросторнее, пошустрее будет работать, это как со строками, вы же не работаете с 1 строкой на 3килобайта(аски например) я надеюсь, обычно нет и первая оптимизация в работе со строками, иметь на каждую строку свой адрес
а в сумме они могут достигнуть 3килобайта например, а курсор в этом контексте ходит по 1 символу и строке
wataru
07.09.2025 00:54ну МФ в последовательности семплов так?
Пока да.
а значит если настроить равное поступление семплов например на 1/3 больше или просто чтобы равное поступало поделив например на 3,
А тут я вас потерял. Что за "поступление" семплов? Поэтому неясно, как оно может быть равным и на 1/3 больше, ведь я не знаю, что за поступление.
Jijiki
07.09.2025 00:54всё зависит от того как вы читаете поток предположу что он состоит из семплов
в любом случае поступает адрес согласны и длинна?
вот я выше пишу что надо прочитать из 1 адреса в трёх адресное равномерное пространство, а писать как будто была бы последовательность этих трёх адресов
123456789 fn(+1) 123 456 789 2345678910поидее же 1 итерация 3 вычисления и запись уже по знакомому адресу
Jijiki
07.09.2025 00:54из вики пример 1
x = [2 50 6 3] -- init *x1 [2 2 50] y[0] *x2 [2 50 6] y[1] *x3 [50 6 3] y[2] *x4 [6 3 3 ] y[3] calculate for int=0;i<1 fn(*x1) y[i] fn(*x2) y[i+1] fn(*x3) y[i+2] fn(*x4) y[i+3]
wataru
07.09.2025 00:54Все еще ничего не понятно. Вы так и не остановились и не выдали никаких связей по шагам. Вы продолжаете развивать свою гениальную мысль, но оно вообще далека от понимания.
в любом случае поступает адрес согласны и длинна?
Куда поступает? На вход программе решающей задачу? Поступает одномерный массив и число W - размер окна. Можно считать что массив это адрес да. Если вы это имелли ввиду адресом и длиной, то да. Если нет, то разъясните.
Jijiki
07.09.2025 00:54ну и получается ждём завершение всех потоков и читаем по 1024
в результате будет выход с медианами
прям как в курле, есть пример на просторах интернета
а вы хотите по очереди выполнять?
wataru
07.09.2025 00:54Я заканчиваю эту ветку. Не буду дальше разбираться - надоело. Вы каждый раз вместо того чтобы разъяснить вашу прошлую мысль и ответить на мои вопросы выдаете что-то новое и не понятное.
Jijiki
07.09.2025 00:54значит надо разобраться возможно дерево это проигрыш, есть источники как работать с последовательностями на С, почему 1 большой поток это тяжелее чем что-то поделённое и почему это легко запустить через пул потоков например
можно в разные адреса загнать, можно на пул потоков(по 4 потока) и есть TimSort там косвенно об этом же
извините я тоже перестану
krotos139
07.09.2025 00:54Вот тут проблема. А как найти это место?
Обычно фильтры применяются к значениям которые плавно меняются. По этому мы можем взять последний добавленный элемент и сравнить новое значение с ним. Далее используя связи на следующий и предыдущий элемент последовательно сравнивая каждый элемент, перемещая указатель в нужном направлении, найти нужное место. В лучшем случае это будет О(1), в худшем случае это будет О(размера окна).
wataru
07.09.2025 00:54Если вы применяете фильтры, то значения уже зашумленные. Они могут чуть-чуть колебаться в окресности истинного гладкого значения, но при этом идти в полную перемешку, а значит новый элемент может попасть в любое место в отсортированном списке. Вот это самое "В лучшем случае это будет О(1)" - это на практике будет весьма редкая вещь.
nickolaym
07.09.2025 00:54Если последовательность монотонно растёт, то мы постоянно добавляем в левую кучу новые элементы (на вершину), а старые помечаем как выбывшие. Но они никогда не окажутся наверху кучи, и поэтому никогда не будут удалены. В результате у нас разрастается и куча, и мусорная корзина. И мы получаем и утечку, и логарифм от маляра Шлемюэля.
Или нет?
nickolaym
07.09.2025 00:54От логарифма на каждом шаге мы никак не избавимся, потому что heappush. Мы хотели бы избавиться от n или даже nlogn при удалении произвольного элемента из кучи.
Но давайте-ка сделаем совсем по-другому.
Пусть у нас есть упорядоченный контейнер с логарифмическим поиском-удалением-вставкой. То есть, любое ДДП (хоть АВЛ, хоть КЧД, неважно). И пусть у этого контейнера есть двунаправленные итераторы! Да-да, std::map.
И вот пусть мы знаем не только значение медианного элемента (или пары смежных), но и итератор на медианный элемент.
Теперь, добавляем новое значение в мап. Пусть мап сам занимается своей балансировкой, нам пофигу. Но мы ведь и понимаем, добавили в левое или правое крыло относительно значения медианы. Поэтому сдвигаем итератор медианы в нужную сторону, и это амортизированная O(1).
А потом удаляем старое значение из мапа. И мы также понимаем, из левого или правого крыла. И опять сдвигаем итератор.
Alexandroppolus
Если массив уже отсортирован по возрастанию, то куча "max_heap" будет постоянно расти, потому что "удаленные" элементы внутри неё никогда не станут её корнями.
Навскидку, можно попробовать хранить в кучах не сами значения элементов, а их индексы в массиве. В хэш-таблице для каждого индекса хранить его позицию в куче (обновлять эти данные при просеивании), потому что зная эту позицию, можно быстро удалить элемент из кучи, даже если он не корневой. Причем вместо хэш-таблицы тогда достаточно массива длиной k, ведь "хэш" для индекса можно считать как i % k. У нас в любой момент времени всего k элементов в работе, и когда вываливается i-й элемент с хэшем h1 = i % k, то добавляется элемент h2 = (i + k) % k = h1, то есть заезжает в ту же ячейку