Синусоидальное моделирование и опечатки в Калтехе
Этот пост про относительно новый метод обработки сигналов, описанный в статье Adaptive data analysis via sparse time-frequency representation, а также про крохотную, но сбившую лично меня с толку, ошибку. Сею статью опубликовали в 2011 году профессора прикладной математики Калифорнийского Технологического института Томас И. Хоу и Зуокьянг Ши, и, вероятно, к моменту, как вы это читаете, они уже её поправили.На эту статью я наткнулся в поиске различных методов частотно-временного анализа нелинейных и нестационарных сигналов — в моем случае ультразвуковых сигналов от передвигающихся форменных элементов крови в сосудах человека. Суть такого анализа состоит в отслеживании изменений характеристик сигнала, иначе говоря, мы хотим знать зависимость составляющих сигнал частот от времени. За исключением широко распространенных методов — спектрального и вейвлет-анализа, были найдены такие методы как EMD (разложение на эмпирические моды) и синусоидальное моделирование, о котором далее пойдет здесь речь.Метод эмпирических мод довольно прост в применении, однако не особо развит с точки зрения обоснованности полученных результатов. Томас Хоу и Зуокьянг Ши пошли дальше в развитии математического аппарата и предложили свой метод синусоидального моделирования сигнала. Его идея заключается в разреженной декомпозиции сигнала на гармоники с гладкими амплитудами. Какой результат мы ожидаем получить — на картинке выше. В данном случае раскладывался сигнал, полученный функцией f (t) = 6t + cos (8πt) + 0.5 cos (40πt). Разложение сигнала, естественно, не уникально, поэтому был введен критерий минимума составляющих гармоник, и задача сформировалась следующим образом:
Решением задачи P0 будут гармоники c амплитудами a и фазовыми функциями theta, чьи производные будут соответствовать частотам сигнала. Авторы предлагают решать эту задачу рекурсивно: извлекаем одну гармоническую составляющую, получаем остаток (медиану) — работаем с ним, извлекаем из него следующую гармонику, получаем остаток и так далее. Так как медиана должна быть гораздо более гладкой, чем сигнал, нахождение гармоники сводится к решению следующей задачи: 
где a1*cos (theta1(t)) — полученная гармоника, a0 — медиана. Гладкость амплитуды определяется общей вариацией: 
Чем больше гладкость функции g, тем ближе к нулю её (n+1) производная, тем, соответственно, меньше общая вариация n-го порядка. Можно сказать еще так: чем меньше порядок полинома, которым мы можем интерполировать, тем глаже функция. Логично предположить, что лучше всего будет искать минимум вариации нулевого порядка, однако, в этом случае мы рискуем получить кусочно-постоянную функцию (staircase effect) и потерять некоторую нужную нам информацию более высокого порядка (например, кривизну). Авторы статьи предлагают уменьшать вариацию третьего порядка, устремляя четвертые производные к нулю, тем самым делая нечто похожее на аппроксимацию кубическими сплайнами. В результате, мы ищем min TV^3(a0). Но нам также важно, чтобы и a1 была достаточно гладкой, поэтому в лучшем случае мы должны найти min TV^3(a0) + min λ*TV^3(a1), где λ — множитель Лагранжа, играющий в данном случае роль весового коэффициента. Для упрощения авторы оставили λ = 1 и пришли к следующей задаче: 
Алгоритм
Для решения этой задачи предложен итеративный алгоритм. В условии добавляется еще одна гармоника — b1*sin (theta1(t)), её гладкость также становится нашей целью. Инициализируем фазовую функцию theta^0 таким образом, чтобы она была возрастающей функцией, удовлетворяющей условию cos (theta^0) = f (t), и запускаем основную итерацию: 
На первом шаге минимизируем общие вариации амплитуд. На следующем — обновляем theta по вышеозначенной формуле, затем проверяем, насколько она изменилась. К концу алгоритма предполагается, что амплитуда b1 будет уже пренебрежимо мала, соответственно, мал будет artcg (b1/a1), и theta сойдется к какому-то значению. Это и будет наша фазовая функция. Получив theta^n, мы выделяем первую гармонику a1*cos (theta1(t)), где theta1 = theta^n. Далее запускаем этот алгоритм на остатке, то есть работаем со значениями f (t) = a0 + b1 * sin (theta1(t)). Выделяем следующую гармонику a2*cos (theta2(t)), продолжаем с остатком. И так до тех пор, пока норма остатка не будет пренебрежимо мала. Я брал первую норму. Для тех, кто не помнит её определение: 
Для применения задачу нахождения минимума общей вариации авторы свели к задаче L1-минимизации, заменив дифференциальное уравнение разностным. У нас имеется сигнал, а точнее упорядоченный набор из N точек (узлов), который мы хотим разложить, то есть представить как сумму векторов a0, a1*cos (theta1(t)) и b1*sin (theta1(t)). И в переформулированной задаче была замечена злосчастная опечатка: 
Предполагая, что a0^n, a1^n, b1^n — сеточные функции, мы аппроксимируем их четвертые производные умножением их на матрицу D^4. Затем, минимизируя сумму модулей полученных значений, мы уменьшаем соответствующие вариации. Все логично и правдоподобно, но, реализовав сей алгоритм, я ужаснулся, увидев результаты — получалось нечто совершенно хаотическое. Долго пытаясь обнаружить ошибку в коде, я, наконец, осознал, что ошибка была в статье.Дело в том, что при умножении Phi на x мы получаем вектор, где элементами являются суммы производных a0^n, a1^n, b1^n в соответствующих точках. То есть мы минимизируем общую вариацию от (a0^n + a1^n + b1^n), а не сумму трех вариаций, как было указано сверху.Ошибка есть ошибка, нашел, исправил и забыл на несколько месяцев. И вот на майских праздниках меня все-таки дернуло написать письмо профессорам с предложением заменить в статье min ||Phi*x|| на min ||D^4*a0^n|| + ||D^4*a1^n|| + ||D^4*b1^n||. Удивительно, уже через час пришел ответ от профессора Хоу с фразой «Thanks for your email. We will study your note and get back to you soon» и простой подписью «Tom».Еще через 6 часов пришло письмо от профессора Ши:
Dear Aleksandr,
Thanks a lot for your email. We did make a typo here.\Phi should be diag (D4, D4, D4) which was also what we used in thealgorithm. Thanks.
Best, Zuoqiang Shi
Естественно, моя поправка и то, что они хотели написать, идентичны. Единственное, должен заметить, что, если не хранить матрицу Phi в предположении, что она разряженная, то есть хранить кучу нулевых элементов, память может скоро закончиться.Как решается задача L1-минимизации
Для своего алгоритма профессора предлагают методы split Bregman iteration и L1-regularized least square method. Однако, эти методы достаточно свежи и, вероятно, зашиты далеко не во все пакеты. Кому лениво (как мне), реализовывать какой-нибудь из алгоритмов нахождения минимума первой нормы, тому будет полезно узнать, как эта задача сводится к более распространенной задаче линейного программирования. Исходная задача выглядит следующим образом: 
при условии
где А — матрица NxM, Phi — матрица KxM, x — вектор размерности M.Она эквивалентна следующей задаче линейного программирования:
при условиях
где phi (i, j) — элементы матрицы Phi.Формально переход можно сделать следующим образом. К вектору x добавляется K переменных, тем самым он становится размерности K+M:
К матрице А добавляется нулевая матрица размерностью NxK: 
Вектор с (коэффициенты в задаче ЛП) состоит из M нулей и K единиц: 
Матрица F составляется из матрицы Phi и единичной матрицы I. Размерность F — (2*K)x (M+K): 
Таким образом, формализованная задача выглядит так:
при условиях
Результаты
Вот каковы результаты работы алгоритма для ранее указанного примера f (t) = 6t + cos (8πt) + 0.5 cos (40πt). На графиках полученные этим алгоритмом значения (синие) сравниваются с аналитическими значениями (красные) и значениями, полученными разложением на эмпирические моды (черные).
На следующем примере видно, как этот метод помогает отслеживать частотные изменения сигнала.
В статье также указаны и другие примеры, а также результаты работы алгоритма на реальных данных.Итог
Алгоритм весьма интересен и, на мой взгляд, требует дальнейшего анализа, благо простора для раздумий здесь полно. Его явными минусами являются большая вычислительная сложность, а также чувствительность к шуму и, как следствие, нестабильная работа при зашумленных данных. Авторы предлагают разрешать проблему чувствительности при помощи низкочастотной фильтрации и преобразования Фурье. Но об этом в другой раз.Напоследок хочется лишний раз похвастаться и сказать, что, видимо, и такое случается, что профессора в Калтехе допускают пускай и ничтожные, но ошибки, и русский студент, не имеющий пока и бакалавра, может им на это указать.
