Преобразование Гильберта-Хуанга

ЦИФРОВАЯ ОБРАБОТКА СИГНАЛОВ

Тема 24. ПРЕОБРАЗОВАНИЕ ГИЛЬБЕРТА-ХУАНГА

Судьба новой истины такова: в начале своего существования она всегда кажется ересью.

Томас Генри Гексли.

Эту стадию Хуанг уже прошел. Вытирать об него ноги математики прекратили и скопом ринулись обосновывать новый метод. А практикам понравилось: работает хорошо, интеграл всего один, и тот позаимствован у Гильберта. Наши, правда, опять лопухнулись, мир вопит от восторга, а мы интересуемся – о чем галдеж?

Игорь Бреднев. Уральский  геофизик.

Содержание.

Введение.

Рекомендуемые материалы

1. Эмпирический метод декомпозиции (EMD) сигналов. Функции внутренних мод сигналов.  Процесс отсева функций IMF.  Критерии останова процесса декомпозиции.  Ортогональность базиса декомпозиции. Примеры практического применения EMD.

2. Спектральный анализ гильберта (HSA).

3. EMD шумовых сигналов. EMD  «белого шума».  Частотная избирательность EMD. 

Введение

            Под преобразованием Гильберта-Хуанга (Hilbert-Huang transform – HHT)  понимается эмпирический метод декомпозиции (EMD) нелинейных и нестационарных процессов и Гильбертов спектральный анализ (HSA). HHT представляет собой частотно-временной анализ данных (сигналов) и не требует априорного функционального базиса. Функции базиса получаются адаптивно непосредственно из данных процедурами отсеивания EMD. Мгновенные частоты вычисляются от производных фазовых функций Гильбертовым преобразованием функций базиса. Заключительный результат представляется в частотно-временном пространстве /1/.

            EMD-HSA был предложен Норденом Хуангом в 1995 в США (NASA) для изучения поверхностных волн тайфунов, с обобщением на анализ произвольных временных рядов коллективом соавторов  в 1998 г. /2/. В последующие годы, по мере расширения применения EMD-HSA для других отраслей науки и техники, вместо термина EMD-HSA был принят более короткий термин преобразования HHT.

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

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

24.1. ЭМПИРИЧЕСКИЙ МЕТОД ДЕКОМПОЗИЦИИ (EMD) СИГНАЛОВ /1,2,3/

Функции внутренних мод сигналов.  Декомпозиция сигналов основана на предположении, что любые данные состоят из различных режимов ­колебаний. В любой момент времени данные могут иметь много различных сосуществующих режимов колебаний, ­нанесенных одно на другое. Каждый режим, линейный или нелинейный, представляет простое колебание, которое имеет экстремумы и нулевые пересечения. Кроме того, колебание будет в определенной степени «симметрично» относительно локального среднего значения. Результат – конечные  сложные данные.

Каждый из этих колебательных режимов может быть представлен функцией внутренней моды (intrinsic mode function - IMF) со ­следующим определением:

1. Число экстремумов и число нулевых пересечений функции должны быть равными или отличаться не более чем на 1.

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

IMF представляет собой колебательный режим, как часть простой гармонической функции, но вместо постоянной амплитуды и частоты, как в простой гармонике, у IMF могут быть переменная амплитуда и частота, как функции независимой переменной (времени, координаты, и пр.). Любую функцию и любой произвольный сигнал можно разделить на семейство функций IMF, придерживаясь изложенной ниже методики. Для наглядности методику реализации EMD рассмотрим на примере разложения цифрового массива модельного сигнала y(k), представленного на рис.24.1.1. Сигнал смоделирован суммой трех нестационарных по амплитуде гармоник различной частоты на интервале отсчетов по k от 0 до 200, и продлен на начальном и конечном участках на интервалы tp=4 для задания начальных и конечных условий преобразования и устранения ошибок преобразования на концевых интервалах обрабатываемого массива данных.  

Процесс отсева функций IMF.  Алгоритм эмпирической декомпозиции сигнала складывается из следующих операций его преобразования.

Операция 1. Идентифицируем по координатам и амплитудам все локальные экстремумы (максимумы и минимумы) сигнала (рис. 24.1.2). Группируем раздельно массивы векторов координат (номеров отсчетов) х_max(k) и соответствующих амплитудных значений у_max(k) максимумов, и аналогичные массивы векторов x_min(k) и y_min(k) минимумов всех выделенных экстремумов.

            Операция 2.  Кубическим (или каким либо другим) сплайном вычисляем верхнюю и нижнюю огибающие сигнала по выделенным максимумам и минимумам, как это показано на рис. 24.1.3 (красный и синий цвет соответственно). Определяем функцию средних значений m₁(k) между огибающими (черный цвет) и находим первое приближение к первой функции IMF:

h₁(k) = y(k) – m₁(k).           (24.1.1) 

Операция 3.  Повторяем операции 1 и 2, принимая вместо y(k) функцию h₁(k), и находим второе приближение к первой функции IMF – функцию h₂(k).

h₂(k) = h₁(k) – m₂(k).                                                (24.1.2)

            Аналогично находим третье и последующие приближения к первой функции IMF. По мере увеличения количества итераций функция m_i(k), равно как и функция h_i(k), стремится к неизменяемой форме. С учетом этого, естественным критерием останова итераций является задание определенного предела по нормализованной квадратичной разности между двумя последовательными операциями приближения, определяемой как

d = S_k |h_i-1(k) - h_i-1(k)|² / S_k h_i-1(k)².                               (24.1.3)

Пример изменения значений d в процессе итераций приведен на рис. 24.1.4. При пороге d = 0.0001 количество итераций, как правило, не превышает 6-8.

Останов итераций по нормализованной квадратичной разности (24.1.3) исторически был первым по применению.  Однако для сложных и объемных (по количеству отсчетов) сигналов в процессе итераций может изменяться количество выделяемых экстремумов, при этом, как правило, наблюдаются скачки значения d в большую сторону, начиная с которых снова начинается процесс уменьшения d, и общее количество итераций может увеличиваться до 20-30 без существенного повышения качества отсева. Опыт показывает, что для оптимальных отсеиваний число итераций порядка 6-8 является вполне достаточным, и все большее применение находит останов по максимальной итерации, при которой число экстремумов и нулевых пересечений функции h(k) не изменяется по сравнению с предыдущими итерациями.

Последнее значение h_i(k) итераций принимается за наиболее высокочастотную функцию с₁(k) = h_i(k) семейства IMF, которая непосредственно входит в состав исходного сигнала y(k). Это позволяет вычесть с₁(k) из состава сигнала и оставить в нем более низкочастотные составляющие (показано на рис. 5):

r₁(k) = y(k) – c₁(k).           (24.1.4)   

Функция r₁(k) обрабатывается как новые данные по аналогичной методике с нахождением второй функции IMF – c₂(k), после чего процесс продолжается:

r₂(k) = r₁(k) – c₂(k),   и т.д.                                              (24.1.5)  

            Таким образом, достигается декомпозиция сигнала в n – эмпирическом приближении:

y(t) =   c_n(t)+r_n(t).                                                  (24.1.6)

Критерии останова процесса декомпозиции сигнала могут быть следующими:

1. Остаток r_n(k) во всем интервале задания сигнала становятся несущественными по своим значениям по сравнению с сигналом.
2. Остаток r_n(k) становится монотонной функцией, из которой больше не может быть извлечено функций IMF.
3. Так как в конечном итоге суммирование всех функций IMF (реконструкция сигнала) должно давать исходный сигнал, то можно останавливать разложение заданием относительной погрешности среднеквадратической реконструкции (без учета остатка r_n(k)) .
4. По мере увеличения количества функций IMF относительная среднеквадратическая погрешность реконструкции достаточно сложных и протяженных сигналов уменьшается, но, как правило, имеет определенный минимум. По-видимому, это определяется попытками алгоритма разложить остаток на функции, частично компенсирующие друг друга. Соответственно, останов программы может выполняться, если следующая выделенная функция IMF увеличивает погрешность реконструкции.

Другими словами, остановка декомпозиции сигнала должна происходить при максимальном «выпрямлении» остатка, т.е. превращения его в тренд сигнала по интервалу задания с числом экстремумов не более 3. Даже для данных с нулевым средним значением конечный остаток может отличаться от нуля. Чтобы применять метод EMD, центрирования данных не требуется, метод нуждается только в локализациях экстремумов. Нулевая линия для каждого компонента декомпозиции будет сформирована процессом отсеивания.

На рис. 24.1.6 приведен пример полной декомпозиции сигнала с остановом по критерию 2. На верхнем графике рисунка приведен входной сигнал преобразования (красным) и сигнал обратной реконструкции (пунктиром) суммированием функций разложения c_i (c1-c5).

Компоненты EMD обычно физически значимы, поскольку ­характеристические параметры функций IMF определяются материальными данными.

Ортогональность базиса декомпозиции. Таким образом, входной сигнал y(k) в соответствии с выражением (24.1.6) раскладывается по базису, который, не определен аналитически, но удовлетворяет всем традиционным требованиям базиса. На основании проверки на модельных и опытных данных он является:

            - законченным и сходящимся (сумма всех функций IMF и остатка равна исходному сигналу и не зависит от критериев останова итераций),

            - ортогональным (все IMF и остаток ортогональны друг другу),

            - единственным.

            И, что самое главное – он является адаптивным, так как получен непосредственно из анализируемых данных эмпирическим методом.

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

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

Примеры практического применения EMD.  В качестве примера в работе /3/ приводится EMD-анализ данных девиации периода вращения Земли. В результате исследований всем выделенным функциям IMF сопоставлены определенные физические процессы, которыми и вызвано их формирование (влияние штормов и тайфунов, месячных вариаций мощности приливов, явления Эль-Ниньо и прочие факторы). Ниже, без комментариев, приводятся выборки из результатов данного анализа, демонстрирующие свойства базиса.

Рис. 24.1.8. Данные.

Рис. 24.1.9. Семейство IMF.

Рис. 24.1.10. Данные и с12 IMF.

Рис. 24.1.11. Данные и сумма с10+с11+с12.

Рис. 24.1.12. Детализация данных и сумма с8+с9+с10+с11+с12.

Рис. 24.1.13. Детализация данных и сумма с7+с8+с9+с10+с11+с12

            Еще один пример из этой же работы /3/.

Глобальная температурная аномалия. Ежегодные данные с 1856 до 2003

Средние значения IMF за 10 просеиваний: СС(1000, I)

Данные и тренд C6

24.2. СПЕКТРАЛЬНЫЙ АНАЛИЗ ГИЛЬБЕРТА (HSA) /1,2/

            IMF, определенные вышеприведенным способом, допускают вычисление физически значимых мгновенных частот, что дает возможность создать частотно-временное представление сигнала на основе преобразования Гильберта.

Для действительного сигнала s(t),преобразование Гильберта определяется главным значением  (PV) интеграла

                                  (24.2.1)

Это приводит к аналитическому сигналу

    z(t) = s(t) + j Hs(t) = a(t) exp(jq(t)),                                   (24.2.2)

                       (24.2.3)

где a(t) и q(t)- соответственно, амплитуда и фаза сигнала во время t. Мгновенная частота может быть вычислена по формуле

w(t) = dq(t)/dt.                                                  (24.2.4)

Заметим, что вышеприведенные и амплитуда и мгновенная частота - функции времени. Мы надеемся создать частотно-временное представление сигнала  с использованием преобразования Гильберта. После выполнения преобразования Гильберта на каждой компоненте IMF первоначальные данные могут быть выражены как действительные в следующей форме:

                         (24.2.5)

Здесь, остаток r_n был не учтен, поскольку это - или монотонная функция или постоянная.

Те же самые данные в представлении Фурье:

                                           (24.2.6)

Рис. 24.2.1.

На рис. 24.2.1 приведено сопоставление частотно-временного представления модельного сигнала в трех представлениях.

24.3. EMD ШУМОВЫХ СИГНАЛОВ

Шумы, сопровождающие полезную информацию в сигнале, в принципе, не относятся к типу колебательных в прямом смысле этого понятия. Но они полностью удовлетворяют приведенным выше определениям функций IMF. При распределении во всем частотном диапазоне входного сигнала и выполнении EMD, они распределяются по всем функциям IMF. Так как информация в главном частотном диапазоне дискретных сигналов обычно является низкочастотной, то шумы «отсеиваются», в основном, в высокочастотные функции IMF. Но в эти функции могут "просачиваться" и высокочастотные спектральные составляющие информационной части сигнала в зависимости от их положения в главном частотном диапазоне. Соответственно, на первый план выдвигается задача формирования определенных критериев отбора только шумовых функций IMF (для исключения их при последующей реконструкции сигнала) и влияния на этот отбор как статистических и спектральных характеристик самих шумов, так и спектрального состава полезной информации в сигнале.

EMD  «белого шума».  Белый шум является стационарным случайным процессом q(t), автокорреляционная функция которого описывается дельта - функцией Дирака и, соответственно, спектральная плотность мощности шумов не зависит от частоты и имеет постоянное значение Wq(f) = s², равное дисперсии значений q(t). По существу, это идеализированный случайный процесс с бесконечной энергией. Но в случае постоянства спектральной плотности мощности случайного процесса в конечном диапазоне частот это существенно упрощает анализ сигналов. Многие помехи в радиотехнике, в технике связи и в других отраслях техники рассматривают как белый шум, если эффективная ширина спектра сигналов много меньше эффективной ширины спектра шумов и спектральная плотность мощности шумов слабо изменяется в интервале спектра сигнала.

Понятие "белый шум" определяет только спектральную характеристику случайного процесса, а, следовательно, под это понятие подпадают любые случайные процессы, имеющие равномерный энергетический спектр и различные законы распределения. На рис. 24.3.1. Приведена гистограмма единичной реализации модельного «белого шума» y(k) в системе Mathcad (5000 отсчетов) с равномерным распределением отсчетов от -0.55 до +0.55 и дисперсией 0.1. Спектральную плотность мощности модельной реализации шума можно посмотреть на рис. 24.3.5.

Рис. 24.3.2.

На рис. 24.3.2 приведен результат EMD модели y(k) шума на первых 750 отсчетах. Шумовой сигнал данной реализации был разложен на 12 функций IMF, первые девять из которых приведены на рисунке (4 последних функции в увеличенном масштабе). Останов процесса декомпозиции выполнен по минимуму погрешности реконструкции без учета остатка, процесс снижения погрешности при увеличении количества функций IMF приведен на рис. 24.3.3. Количество функций IMF в различных реализациях случайного сигнала изменяется от 8 до 14. Останов итераций при вычислении каждой функции IMF был установлен по относительному расхождению между последовательными итерациями с порогом 0.01%, при этом количество итераций для первых функций IMF достигает 30 и, как правило, постепенно снижается. Погрешность реконструкции с учетом остатка практически равна нулю. Вычислением скалярного произведения любых двух функций IMF можно убедиться в их взаимной ортогональности.

На рис. 24.3.4-А приведены гистограммы первых четырех IMF в сопоставлении с гистограммой входного сигнала y(k). Как следует из этих графиков, EMD существенно изменяет плотности распределения выходных функций. Распределение первой IMF становится двумодальным с прогибом вниз на малых (близких к нулевым) амплитудах. Это объясняется тем, что для рядом расположенных однополярных импульсов при EMD выделяются экстремумы импульсов большей амплитуды, которые и отсеиваются в первую функцию IMF. При вычитании этой функции из входного сигнала распределение оставшейся части шумов становится близким к гауссовому с нулевым средним значением и резким сокращением рядом расположенных однополярных импульсов. На отборе всех последующих IMF этот фактор уже не сказывается, и они имеют распределение, близкое к гауссовому, а рядом расположенные однополярные импульсы воспринимаются, как более низкочастотные составляющие шума. Статистика последовательной реконструкции шумового сигнала частично показана на рис. 24.3.4-В.

Рис. 24.3.4.

Проверка процесса EMD на шумовых сигналах с другими законами распределения (Гаусса, Пуассона и пр.) показала, что качественный характер процесса остается практически неизменным.

На рис. 24.3.5 приведены спектры плотности мощности сигнала y(k) и первых семи функций IMF в главном частотном диапазоне сигнала 0-p (5000 отсчетов шума с Dk = 1, отсчеты по спектру 0-2500 с шагом Dw = p/2500).

Рис. 24.3.5.

            Частотная избирательность EMD.  По рис. 24.3.5 видно, что процесс EMD обладает вполне определенной частотной избирательностью на каждом уровне EMD. Но говорить о каких-либо частотных передаточных функциях EMD, по-видимому, будет некорректным, так как любая частотная составляющая w_i исходного сигнала в процессе EMD может быть расщеплена по амплитуде и фазе на составляющие разных уровней функций IMF. Это можно видеть на рис. 24.3.6, где приведены графики модулей «эквивалентных» частотных передаточных характеристик разложения для первых пяти функций IMF, полученные осреднением отношения спектров функций к спектру исходного сигнала. Для получения достаточно гладких «эквивалентных» передаточных функций входной белый шум реализовался массивом из 30000 отсчетов, а сглаживание отношения спектров выполнялось в скользящем временном окне 2000 отсчетов.

Рис. 24.3.6.

            На рис. 24.3.7 приведены графики последовательного суммирования коэффициентов «эквивалентных» передаточных функций, которые показывают процесс последовательного перекрытия всего частотного диапазона входного сигнала.

Рис. 24.3.7.

Обратным преобразованием Фурье по спектрам мощности могут быть вычислены нормированные автокорреляционные функции семейства IMF, первые 5 из которых приведены на рис. 24.3.8. Как следует из графиков, статистическая независимость отсчетов в какой-то мере сохраняется только для первой IMF. Но даже в ней появляется отрицательная (знакопеременная) корреляция между последовательными отсчетами. Во всех остальных функциях четко прослеживается появление затухающей косинусоидальной зависимости между отсчетами с последовательным увеличением интервала корреляции по мере увеличения номера IMF.

Рис. 24.3.8.

Информация в лекции "71 Революция 1848-1849 гг. в чешских землях" поможет Вам.

Литература.

1. The Hilbert-Huang transform and its applications / editors, Norden E. Huang, Samuel S.P. Shen. - World Scientific Publishing Co. Pte. Ltd. 5 Toh Tuck.   Link, Singapore 596224

2. Huang, N. E., Z. Shen, S. R. Long, M. C. Wu, H. H. Shih, Q. Zheng, N.-C. Yen, С. C. Tung, and H. H. Liu, 1998: The empirical mode decomposition and the Hilbert spectrum for nonlinear and non-stationary time series analysis. Proc. R. Soc. London, Ser. A, 454, 903-995.

3. An Introduction to Hilbert-Huang Transform: A Plea for Adaptive Data Analysis. Norden E. Huang. Research Center for Adaptive Data Analysis.  National Central University.

О замеченных опечатках, ошибках и предложениях по дополнению:  davpro@yandex.ru.

Copyright © 2010 Davydov А.V.