В этой статье наши инженеры хотели бы поделиться с Хабром достаточно интересным инструментом, который можно эффективно применять для фильтрации зашумленных сигналов, пользуясь априорным знанием об унимодальности сигнала.
Задача оффлайновой фильтрации сигналов в случае, когда ожидаемая форма сигнала известна с точностью до нескольких неизвестных параметров, сводится к задаче аппроксимации. Например, если известно, что сигнал линейно растет на рассматриваемом промежутке, задача сведётся к линейной регрессии, а если можно предположить, что шум — нормален, то правильным методом будет МНК. Но однажды мы столкнулись с задачей оценки формы профиля рентгеновского микрозонда (пучка), про которую априори было достоверно известно только одно: профиль унимодален, а именно имеет ровно один максимум. Оказывается, и в этом случае можно наилучшим (в смысле, например, L2 метрики) образом приблизить экспериментальный сигнал функцией, принадлежащей известному множеству (множеству унимодальных функций). Причём — с приемлемой ассимптотикой вычислительной сложности.
===> 
===> 
Итак, речь пойдёт о задаче аппроксимации некоторого зашумлённого сигнала 
Вначале, для простоты изложения, имеет смысл рассмотреть ту же задачу для монотонной функции. Описание приведено для монотонно нестрого возрастающей функции, хотя, разумеется, для монотонно убывающей можно проделать всё то же самое.
Если предположить, что шум нормально распределен, то задача сведётся к минимизации квадратичного функционала 
Существует альтернативное решение, не обладающее указанными недостатками, но при его использовании шкалу сигнала 
Рис. 1: a) Модельный монотонный сигнал до налолжения шума. b) Модельный монотонный сигнал после наложения шума.
Оказывается, в дискретном приближении фильтрация нормально зашумленного монотонного сигнала сводится к динамическому программированию с двумя параметрами, причём, пользуясь этим решением, удаётся построить вычислительную схему той же сложности и для унимодального случая.
О том, что такое динамическое программирование, можно почитать здесь, а как и к каким задачам его применяют — здесь.
Фильтрация монотонного сигнала методом динамического программирования
В рассматриваемой задаче таблица динамического программирования — двумерная, её столбцам i соответствуют возрастающие слева направо значения дискретного аргумента 
Внешний цикл ведётся по координате i слева направо. Для каждого значения 
причём за пределами таблицы функционал считается равным нулю.
При этом, чтобы не пересчитывать каждый раз минимум на меньших значениях j, его полезно сохранять в процессе заполнения вместе со значением функционала, а поскольку в конце требуется восстановить путь, сохранять удобнее значение m, на котором этот минимум был достигнут:
Тогда вычисление самого функционала можно выполнять за O(1) для каждой ячейки:
После завершения прямого прохода остаётся получить аппроксимацию в явном виде. Для этого ищется минимум среди значений функционала для крайнего правого столбца (столбца с наибольшей координатой i), а затем выполняется обратный проход по i, с восстановлением индекса значения отфильтрованного сигнала j* в каждом столбце:
Окончательные значения 
На рис. 2 показан пример результата описанного алгоритма для монотонного сигнала.
Рис. 2: Фильтрация монотонно возрастающего сигнала динамическим программированием. Синим цветом показан входной зашумлённый сигнал, зелёным — результат работы алгоритма, красным — идеальный незашумлённый сигнал.
На этом рассмотрение вспомогательного монотонного случая закончено.
Унимодальный сигнал
Используя полученный результат, можно решить задачу об унимодальном сигнале, заполнив 2 таблицы ДП: одна заполняется как для монотонной функции, вторая — таким же образом, но в обратную сторону. В результате вычисляется 2 функционала: 
На рисунках 3 и 4 показан пример исходного унимодального сигнала и результата его обработки предложенным алгоритмом.
Рис. 3: Зашумлённый унимодальный сигнал.
Рис. 4: Фильтрация унимодального сигнала динамическим программированием. Синим цветом показан входной зашумлённый сигнал, зелёным — результат работы алгоритма, красным — идеальный незашумлённый сигнал.
Ограничение на производную
В случае, если известно, что сигнал не может меняться быстрее некоторой фиксированной скорости 
Правда, чтобы не проиграть в трудоёмкости, нужно иметь возможность получать минимальное значение функционала на отрезке [j — dj, j] за константное время.
Это можно сделать с помощью алгоритма Ван Херка — Гила — Вермана, описанного, например, здесь.
Результат применения алгоритма с ограничением на производную 
Рис. 5: Фильтрация унимодального сигнала динамическим программированием с ограничением на скорость изменения сигнала. Синим цветом показан входной зашумлённый сигнал, зелёным — результат работы алгоритма, красным — идеальный незашумлённый сигнал.
Разоблачение чёрной магии
Пытливый читатель может заподозрить, что таких же или лучших результатов можно достигнуть более простыми способами. Конечно, при определённых параметрах сигнала и шума описанная техника будет излишней. Кроме того, короче та дорога, которую знаешь. Поэтому на каждый хитрый случай у специалиста наверняка найдётся любимый метод, который именно в его руках даст приемлемый результат.
Чтобы рассеять ощущение возможного мошенничества, на рис. 6 и 7 приведены примеры обработки сигнала гауссовским сглаживанием и медианной фильтрацией с зашкаливанием в нуле. Параметры фильтров подбирались с целью получения оптимального результата, но без фанатизма.
Рис. 6: Фильтрация унимодального сигнала гауссовским сглаживанием. Синим цветом показан входной зашумлённый сигнал, зелёным — результат фильтрации, красным — идеальный незашумлённый сигнал.
Рис. 7: Медианная фильтрация унимодального сигнала. Синим цветом показан входной зашумлённый сигнал, зелёным — результат фильтрации, красным — идеальный незашумлённый сигнал.
Как видно, на достаточно произвольных данных метод динамического программирования дает результат получше, чем простейшие общеизвестные фильтры. При этом сравнивать их на модельных примерах подробно — абсолютно бессмысленно. Изложенный выше метод ценен ровно тем, что является примером “костюма, сшитого по мерке”: формализованная постановка задачи была положена в саму основу метода. В рассматриваемом случае это было условие унимодальности. Аналогично, линейные методы “по мерке” сигналам, частотно разделимым с шумом. Еще один инструмент в наборе — что может быть лучше для инженера?
Автор: Smart Engines Ltd.
