Об одном забавном подходе к фильтрации унимодальных сигналов

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

Задача оффлайновой фильтрации сигналов в случае, когда ожидаемая форма сигнала известна с точностью до нескольких неизвестных параметров, сводится к задаче аппроксимации. Например, если известно, что сигнал линейно растет на рассматриваемом промежутке, задача сведётся к линейной регрессии, а если можно предположить, что шум — нормален, то правильным методом будет МНК. Но однажды мы столкнулись с задачей оценки формы профиля рентгеновского микрозонда (пучка), про которую априори было достоверно известно только одно: профиль унимодален, а именно имеет ровно один максимум. Оказывается, и в этом случае можно наилучшим (в смысле, например, L2 метрики) образом приблизить экспериментальный сигнал функцией, принадлежащей известному множеству (множеству унимодальных функций). Причём — с приемлемой ассимптотикой вычислительной сложности.

===> ===>

Итак, речь пойдёт о задаче аппроксимации некоторого зашумлённого сигнала $Y$ как функции дискретного аргумента $X$ , причем известно, что исходный сигнал — унимодальный. В рассматриваемом случае дополнительно известно, что на краях сигнал ассимптотически приближается к нулю. Следует отметить, что рассматриваемый ниже метод сам по себе не требует такого допущения.

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

Если предположить, что шум нормально распределен, то задача сведётся к минимизации квадратичного функционала $\hbox{\sum_i(\widetilde{Y}(X_i) - Y(X_i))^2}$ в присутствии неравенств вида $\widetilde{Y}(X_{i+1}) \ge \widetilde{Y}(X_i)$ . Это — задача квадратичного программирования, и для нее на просторах интернета можно найти некоторое количество готовых решателей. Но, во-первых, задача квадратичного программирования в общем случае обсчитывается довольно медленно, а во-вторых, задача аппроксимации унимодального сигнала не выписывается в подобном виде, если положение экстремума не известно заранее.

Существует альтернативное решение, не обладающее указанными недостатками, но при его использовании шкалу сигнала $\widetilde{Y}$ придётся считать дискретной. Тут стоит отметить, что задача не интересна в случае слабого шума, поскольку обычная линейная фильтрация, скорее всего, снимет все практические проблемы. Что же такое интересный, сильный шум? На рис. 1 можно увидеть модельный сигнал до и после такого зашумления. Шум здесь задан нормальным распределением $N(0, 0.8^2)$ . В подобных условиях дискретизация аппроксимирующего сигнала не выглядит существенной потерей. Кроме того, как будет видно далее, при увеличении числа отсчетов сложность вычислений будет возрастать линейно.

Рис. 1: a) Модельный монотонный сигнал до налолжения шума. b) Модельный монотонный сигнал после наложения шума.

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

О том, что такое динамическое программирование, можно почитать здесь, а как и к каким задачам его применяют — здесь.

Фильтрация монотонного сигнала методом динамического программирования

В рассматриваемой задаче таблица динамического программирования — двумерная, её столбцам i соответствуют возрастающие слева направо значения дискретного аргумента $X$ экспериментальных точек, а строкам j — возрастающие снизу вверх возможные значения аппроксимирующей функции $\widetilde{Y}$ . Минимизируемый функционал — сумма квадратов отклонений $\widetilde{Y}$ от $Y$ вдоль путей с допустимыми переходами. Очевидно, что допустимыми будут переходы на 1 вправо с одновременным неотрицательным смещением вверх.

Внешний цикл ведётся по координате i слева направо. Для каждого значения $X$ выполняется внутренний цикл по координате j снизу вверх. Минимальное значение функционала определяется по формуле

$W_{i, j} = \min\limits_{k \le j}(W_{i-1, k}) + (Y_j - \widetilde{Y}_j)^2,$

причём за пределами таблицы функционал считается равным нулю.

При этом, чтобы не пересчитывать каждый раз минимум на меньших значениях j, его полезно сохранять в процессе заполнения вместе со значением функционала, а поскольку в конце требуется восстановить путь, сохранять удобнее значение m, на котором этот минимум был достигнут:

$m_{i, j} = \begin{cases} m_{i, j-1}, & W_{i, m_{i, j - 1}} < W_{i, j} \\ j, & \text{otherwise} \end{cases}$

Тогда вычисление самого функционала можно выполнять за O(1) для каждой ячейки:

$W_{i, j} = W_{i-1, m_{i-1, j}} + (Y_j - \widetilde{Y}_j)^2.$

После завершения прямого прохода остаётся получить аппроксимацию в явном виде. Для этого ищется минимум среди значений функционала для крайнего правого столбца (столбца с наибольшей координатой i), а затем выполняется обратный проход по i, с восстановлением индекса значения отфильтрованного сигнала j* в каждом столбце:

$j^*_{i} = \begin{cases} \operatornamewithlimits{argmin}\limits_{j}(W_{n-1, j}), & i = n - 1 \\ m_{i,j^*_{i+1}}, & i < n - 1 \end{cases}$

Окончательные значения $\widetilde{Y}$ определяются как

$\widetilde{Y}(X_i) = \widetilde{Y}_{j^*_i}$

На рис. 2 показан пример результата описанного алгоритма для монотонного сигнала.

Рис. 2: Фильтрация монотонно возрастающего сигнала динамическим программированием. Синим цветом показан входной зашумлённый сигнал, зелёным — результат работы алгоритма, красным — идеальный незашумлённый сигнал.

На этом рассмотрение вспомогательного монотонного случая закончено.

Унимодальный сигнал

Используя полученный результат, можно решить задачу об унимодальном сигнале, заполнив 2 таблицы ДП: одна заполняется как для монотонной функции, вторая — таким же образом, но в обратную сторону. В результате вычисляется 2 функционала: $W^{lr}$ и $W^{rl}$ , а также значения $m^{lr}$ и $m^{rl}$ для восстановления путей по таблицам. После этого среди всех значений i и j ищется такоая пара $(i_0, j_0)$ , для которой сумма функционалов минимальна — это будут индексы максимума аппроксимированной функции. Значения функции с каждой стороны от него восстанавливаются по соответствующей таблице.

$(i_0, j_0) = \operatornamewithlimits{argmin}\limits_{i, j}(W^{lr}_{i, j} + W^{rl}_{i, j})$

$j^*_{i} = \begin{cases} j_0, & i = i_0 \\ m^{lr}_{i, j^*_{i+1}}, & i < i_0 \\ m^{rl}_{i, j^*_{i-1}}, & i > i_0 \end{cases}$

$\widetilde{Y}(X_i) = \widetilde{Y}_{j^*_i}$

На рисунках 3 и 4 показан пример исходного унимодального сигнала и результата его обработки предложенным алгоритмом.

Рис. 3: Зашумлённый унимодальный сигнал.

Рис. 4: Фильтрация унимодального сигнала динамическим программированием. Синим цветом показан входной зашумлённый сигнал, зелёным — результат работы алгоритма, красным — идеальный незашумлённый сигнал.

Ограничение на производную

В случае, если известно, что сигнал не может меняться быстрее некоторой фиксированной скорости $v$ (что вполне вероятно), точность регрессии можно улучшить, введя это ограничение в алгоритм в явном виде. Для этого достаточно вычислить максимально возможное изменение индекса значения сигнала $dj_{i}$ для каждого i и запретить переходы в таблице ДП с изменением j более, чем на эту величину:

$dj_{i} = \lceil v \cdot (X_{i} - X_{i-1}) / (\widetilde{Y}_{1} - \widetilde{Y}_{0})\rceil$

$W_{i, j} = \min\limits_{j-dj \le k \le j}(W_{i-1, k}) + (\widetilde{Y}_{j} - Y_{j})^2$

Правда, чтобы не проиграть в трудоёмкости, нужно иметь возможность получать минимальное значение функционала на отрезке [j — dj, j] за константное время.
Это можно сделать с помощью алгоритма Ван Херка — Гила — Вермана, описанного, например, здесь.

Результат применения алгоритма с ограничением на производную $v = 1.2$ показан на рис. 5.

Рис. 5: Фильтрация унимодального сигнала динамическим программированием с ограничением на скорость изменения сигнала. Синим цветом показан входной зашумлённый сигнал, зелёным — результат работы алгоритма, красным — идеальный незашумлённый сигнал.

Разоблачение чёрной магии

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

Чтобы рассеять ощущение возможного мошенничества, на рис. 6 и 7 приведены примеры обработки сигнала гауссовским сглаживанием и медианной фильтрацией с зашкаливанием в нуле. Параметры фильтров подбирались с целью получения оптимального результата, но без фанатизма.

Рис. 6: Фильтрация унимодального сигнала гауссовским сглаживанием. Синим цветом показан входной зашумлённый сигнал, зелёным — результат фильтрации, красным — идеальный незашумлённый сигнал.

Рис. 7: Медианная фильтрация унимодального сигнала. Синим цветом показан входной зашумлённый сигнал, зелёным — результат фильтрации, красным — идеальный незашумлённый сигнал.

Как видно, на достаточно произвольных данных метод динамического программирования дает результат получше, чем простейшие общеизвестные фильтры. При этом сравнивать их на модельных примерах подробно — абсолютно бессмысленно. Изложенный выше метод ценен ровно тем, что является примером “костюма, сшитого по мерке”: формализованная постановка задачи была положена в саму основу метода. В рассматриваемом случае это было условие унимодальности. Аналогично, линейные методы “по мерке” сигналам, частотно разделимым с шумом. Еще один инструмент в наборе — что может быть лучше для инженера?

ссылка на оригинал статьи https://habrahabr.ru/post/279427/

Об одном забавном подходе к фильтрации унимодальных сигналов

Фильтрация монотонного сигнала методом динамического программирования

Унимодальный сигнал

Ограничение на производную

Разоблачение чёрной магии

Комментарии

Добавить комментарий Отменить ответ