Пролог

В программировании микроконтроллеров порой приходится прибегать к такому программному компоненту как медианный фильтр. Вот вам яркий пример. Вы измеряете расстояние UWB трансиверами между радиопередатчиками и вдруг у вас образуются единичные редкие отбросы, которые никак не характеризует измеряемую величину. Или вы опрашиваете емкостную кнопку и от скача питания она тоже может дать редкие, но высокие значения семпла. Классическое решение для исключения таких отбросов - это медианный фильтр.

Ситуация осложняется тем, что на курсах ЦОС(DSP) в ВУЗах обычно про медианный фильтр даже не вспоминают, отводя всё время на обзор FIR(КИХ) и немного IIR(БИХ) фильтров. Хотя достоинство МФ в то, что у него всего один конфиг: порядок фильтра и медианный фильтр (МФ) предназначен для того, чтобы как раз исключать одиночные отбросы сигнала.

В этом тексте я произвел разбор решения LeetCode задачи 480. Sliding Window Median в контексте реализации на языке программирования Си.

Постановка задача

Реализовать оптимизированный по быстродействию медианный фильтр на языке программирования Си, чтобы его можно было применять в микроконтроллерных прошивках. Фильтр должен работать в потоковом режиме. На каждом отсчете обрабатывать семпл на входе и выдавать семпл на выходе.

Определения

медиана — это центральный элемент отсортированного массива нечетной длинны (или среднее арифметическое двух центральных элементов в случае массива четной длинны). Медиана определяется, как среднее значение упорядоченного списка.

порядок фильтра (размер окна, число k) - количество отсчетов, которые фильтр анализирует на каждом такте. Проше говоря сколько памяти в этом фильтре. Порядок фильтра еще называют размером окна.

семпл - целое число. Отсчеты ADC. В конечном счете семпл - это какая-либо физическая величина: расстояние, напряжение, ток, освещенность и т.п.

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

нечетное число - это число, которое не делится на два в множестве целых чисел. Пример нечетных чисел: 1; 3; 5; 7; 9; 11; и т д

бинарная куча - Это абстрактная структура данных. У кучи есть только три легальные операции: добавить элемент в кучу, получить корень кучи (максимальный или минимальный элемент), сбалансировать кучу. Вот и всё. Бинарная куча это двоичное дерево.

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

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

Реализация

Текстовое описание алгоритма

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

Основная идея в том чтобы использовать две кучи: максимальную и минимальную.

Обслуживая эти две кучи мы всегда можем получать быстрый доступ к средему элементу. Схематически это выглядит так.

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) для отложенного удаления чисел, которые были логически удалены, но не удалены физически из кучей.

Какой вместительности должны быть эти две кучи? В идеале бесконечного. Однако так не бывает. В случае сильного разрастания кучей можно вручную удалить ненужные более элементы согласно данным из hash таблицы (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ле удаления все равно делать балансировку, которая всё восстановит (закопает обратно). В результате не будет потребности в динамической памяти для куч, а хеш-таблицу и вовсе можно будет выбросить за ненадобностью. Однако это крайне продолжительная операция. По сути та же самая сортировка, что и в классической реализации на сортировке элементов окна.

Важные замечания к алгоритму

--Надо сделать виртуальные размеры куч. int smallSize и largeSize. При добавлении и извлечении элементов увеличивать и уменьшать надо именно эти отдельные переменные. Реальные же размеры кучей не должны участвовать, как управляющие данные для алгоритма.

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

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

Исходный код

На языке программирования Си это особенно трудно, так как в 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 размера окна, то фильтр станет работать в 1106 раз медленнее, чем реализация на сортировке.

5--Если фильтровать случайный поток чисел, то размер хэш-таблицы тоже бесконтрольно увеличивается.
6--Если фильтровать случайный поток 32-битных чисел, то в хэш-таблице происходят постоянные коллизии, которые тоже просаживают производительность. Ситуация сильно упрощается, если мы имеем дела с 8-битными семплами. Тогда в качестве ключа можно использовать сами значения и ликвидировать коллизии как таковые в принципе. Цена решения - 1kByte RAM. Однако для чисел 32 бит такая реализация хэш-таблици потребовала бы от 4GByte до 16GByte RAM.
7--Сумма элементов в верхней и нижней кучах может оказаться больше, чем размер окна (порядок фильтра). Элемент, который надо удалить может уйти далеко от корней кучи или вовсе оказаться крайним листом в дереве.

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

Боюсь, что больше плюсов просто нет.

Итог

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

Без процедуры flush иногда медианный фильтр на кучах работает медленнее. Иногда быстрее. Иногда и вовсе зависает. Все зависит от того как генератор случайных чисел сформирует входные данные.

Реализация МФ на двух кучах не имеет смысла так, как та же qsort на PC реализована с учетом аппаратных возможностей конкретного CISC процессора, чтобы быть максимально производительной.

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

Словарь

Ссылки

Вопросы:

1--Существуют ли аппаратные реализации медианного фильтра в ASIC исполнении? В виде отдельных микросхем с SPI-интерфейсом или вовсе с параллельной шиной.

2--Как быть если добавляемый элемент оказался между верхушками куч? Куда его добавлять? В LargeHeap или в SmallHeap?

3--А что если поставить последовательно друг за другом два медианный фильтра длинны 3? Заменят ли они собой один медианный фильтр порядка 6 ?

4--Существуют ли микроконтроллеры с аппаратным блоком сортировки ? Это может пригодится для потоковой реализации медианной фильтрации .