Адаптивные вейвлеты. (Алгоритм спектрального анализа сигналов)

An adaptive transform technique dim is suitable for spectral analysis of observational data with irregular sampling gaps is considered. Wc present some artificial examples of data with gaps analysed by both .standart and adaptive techniques. We show that common aliasing problems due to finite length of observations and data gap can be suppressed in both large and small scales in favor of the adaptive wavelets algorithm.

Традиционный спектральный анализ, основанный на преобразовании Фурье, эффективен при обработке периодических сигналов. На практике нередко возникает необходимость знать спектральные свойства сигналов, не являющихся периодическими и, к тому же, известными лишь на ограниченном отрезке времени (или ряде отрезков, разделенных пробелами). Зачастую интерес представляют сигналы, спектральный состав которых меняется со временем. В этом случае необходим некий ‘'локальный” анализ спектра.С последней задачей хорошо справляется вейвлет-анаииз. сформировавшийся в самостоятельный ра<дел прикладной математики в последнее десятилетие (см., например, [3.7,9]).Используя разложение по осцилирующим функциям, локализованным как в физическом, так и в фурье-sпространствах, вейвлет-преобразование отображает исходный одномерный сигнал на полуплоскость время-частота, характеризуя спектральный состав сигнала в каждый момент времени.

В данной работе рассматриваются проблемы, связанные с обработкой сигналов, известных с пробелами, когда как фурье-. гак и вейвлет-спектры искажаются разрывами в данных. 'Идея предлагаемого алгоритма возникла при исследовании астрофизических данных [4.5.8]. которые по условиям наблюдения содержат неизбежные пропуски во временных рядах.Первое описание метода дано в работе [4] в приложении к обработке данных по хромосферной активности звезд. Цель данной статьи - привести подробное описание алгоритма с иллюстрациями на искусственных примерах. Предлагаемый подход может быть полезен при обработке самых различных видов данных.

Вейвлеты. Семейство вейвлет-функций генерируется из одной "материнской" функции w(t), называемой также анализирующим вейвлетом, при помощи растяжения (сжатия) и сдвига,

^..Н^^а V -------- •                          (1)

V а 7

где а - масштабный множитель, характеризующий растяжение, а b - сдвиг вейвлета.

В качестве вейвлета y(t) выбирается осцилирующая функция, удовлетворяющая следующим условиям. Во-первых, среднее значение функции должно равнялся нулю

J y/(t)dt O.                                              (2)

Во-вторых, функция должна быть локализована и в физическом, и в фурье-пространствах (время и частота), т.е. области определения функции в физическом пространстве At и в пространстве частот Аю должны удовлетворять условию AtAo = const > 2тг.

Вейвлет преобразование w(a.t) сигнала f(t) определяется как wfa^CX jvf--W)dt ,(3)

v а V где ц/(1) действительный или комплексный анализирующий вейвлет, а знак * означает комплексное сопряжение, с,. = jk4!|v7(6v)l'd

и («) есть фурье-образ вейвлета чЧО у>(<д) = J ^/(t)e ""'dt.(5)

Если С v <ое, то для вейвлет-преобразования существует формула обращения [6]

fit)" С;1'2 J JV'^f^'Lfa.t')1^1(6)

„       ' 8а"

и справедливо соотношение, являющееся аналогом теоремы Парсеваля.

j f,(t)fk(t)dt = С^ J w,(a,t)w*(a,t)^,(7)

о обеспечивающее, в частности, равенство энергии в физическом и вейвлет-пространствах.

Вейвлет-разложение можно выразить через фурье-образы вейвлета >;'(,<») и сигнала f(

w(a.t) = С4а‘" (2т)'1 ^ уЧам^со^йго.(8)

ЙИ = С;,: J J 4т^юЩ! К" —” •

О -со

Это дает возможность связать фурье-спектр сигнала E(a :

М(а) = JlHa.t))2 dt.(10)

С помощью (8) получаем соотношение

М(а)- aj£(w)|7(a69)j"dft>.(11)

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

Подробное описание свойств непрерывного вейвлет-преобразования алгоритмов и разнообразные примеры можно найти в книгах [7.9]. краткий очерк на русском языке - в препринте [2] и статье [1].

Приведем в качестве примеров два наиболее часто применяемых вейвлета, действительный вейвлет, называемый "мексиканской шляпой".

-//(D-fl-t^e"'2.                                  (12)

и комплексный вейвлет Морле (Morlet)

^(t)-е-2,2е^"*.                                         (13)

В наших расчетах используется вейвлет Морле с ш0 = 2,т .

Адаптивные вейвлегы. Пусть имеется стохастический или квазипериодический сигнал f(t), который регистрируется только на некоторых интервалах времени. Иначе говоря, вместо функции f(t) нам известна функция

f(t) = f(t) * G(t) .                                          (14)

где G(t) есть "функция дырок", которая равна единице, если сигнал регистрируется, и нулю во всех остальных точках (т.е. внутри пробелов и вне сигнала). В результате, при вычислении вейвлет-коэффициентов вместо искомых значений w(a.t) (3) получаются коэффициенты w"(a, I) — С'^а'1'3 I )dt',                                (15)

V а / где у/1 — । = у/ — *G(t).

v а ) V а /

Попадая на дыры, функция д-1 перестает удовлетворять требованиям, предъявляемым к вейвлетам, в частности, условию (2), т.е.

L/C^dr^O.                        (16)

Д а )

Идея адаптивных вейвлетов состоит в перенесении проблемы пробелов с неизвестной функции f(t) на хорошо известную функцию \р. и замене "испорченного" вейвлета ^‘ на "адаптированный" вейвлет у , который должен, по меньшей мере, удовлетворять условию < у >=0 и стремиться к исходному вейвлету ц/ при исчезновении пробелов.

На языке фурье-представлений можно дать следующую интерпретацию предлагаемой идеи. Особенностью вейвлетов является то. что их фурье-образ локализован в ограниченной полосе частот. Нарушение условия (2) приводит к появлению в спектре |^'(a,)|' низких частот вплоть до ы-0. а разрывы на краях дыр производят высокочастотный шум. Необходимо указать алгоритм, подавляющий как низкочастотный, так и высокочастотный шум. обусловленный дырами и краями (край будем в дальнейшем рассматривать как полубесконечную дыру).

Для реализации этой идеи представим исходный вейвлет ц/ в виде

V(t)=h(t)^(t) ,                                       (17)

где Ф(Г) есть положительно определенная масштабная функция ("оболочка", в качестве которой обычно используется гауссиан Ф(1) exp(-t2/2)), a h(t) -"заполняющая" функция (парабола для мексиканской шляпы (12), комплексная гармоническая функция для вейвлета Морле (13)).

Искомый адаптивный вейвлет ^ будем искать в виде ^(l)-[h(t)-C^(t)G(t) .                              (5 8)

Константа (?. зависящая от масштаба а и положения центра вейвлета t. выбирается из условия < у >=0 и равна

(г (       7 т’д\

С(а,1)= I Ф — Gft’^t' I h -- Ф!            .                 (19)

(^ v a 7                v a 7 V a 7

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

Рис. 1.

Примеры. Для иллюстрации свойств алгоритма рассмотрим простые сигналы, анализируемые двумя видами вейвлетов (вещественным и комплексным) по стандартной и адаптивной методике (Рис.2-3). В нижней части каждой из восьми панелей, представленных на этих рисунках, показан обрабатываемый сигнал, а в верхней - фазовая плоскость время-частота, на которой интенсивность черного цвета пропорциональна модулю вейвлет-коэффипиентов. Горизонтальная ось на фазовой плоскости - это ось времени в том же масштабе, что и приведенный под ней сигнал. Вертикальная ось - это ось временных масштабов а (или связанных с ними частот а>=2я/а, где a, -a_nda", а а_п выбирается из условия, что период., соответствующий вейвлету наименьшего масштаба, должен охватывать две точки сигнала), причем большим масштабам (малым частотам) соответствует низ плоскости. По вертикальной оси выбран логарифмический масштаб. его является традиционным для представления частотных характеристик.

Первый пример (Рис.2) - это простой гармонический сигнал, наблюдаемый на коротком интервале времени с одним дополнительным пропуском. На рис.2а представлены результаты стандартного вейвлет-преобразования с помощью мексиканской шляпы (121. а на рис. 1с тог же сигнал обработан с помощью вейвлета Морле (13). В первом случае фазовая плоскость имеет типичную для вещественных вейвлетов структуру с вертикальными белыми полосами. которые возникают при смене знака коэффициентов разложения. При использовании комплексного вейвлета структура ноля модулей становится монотонной, а информация о положении максимумов содержится в фазе коэффициентов. которая тдесь не приводится, так как не играет принципиального значения в контексте тех задач, о которых идет речь. Сравнивая рисунки 2а и 2с, следует отметить и более высокое спектральное разрешение вейвлета Морле, которое выражается в том. что темная горизонтальная полоса соответствующая основной частоте сигнала, шачительно сужается на втором рисунке.

Правые панели показывают результат обработки с помощью адаптивных вейвлетов (рис.2Ь - мексиканская шляпа, рис.2d - вейвлет Морле) . В обоих случаях наблюдается значительное ослабление сигнала в нижней части фазовой плоскости.

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

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

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

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

вейвлет-область только

крупномасштабные функции, то на узких дырах (и на широких вблизи краев)

будет воспроизводиться и мелкомасштабная структура сигнала.

Рис. 2

л>ЛЯГИЯ Д.Н ФрИЕ НЛ'

Рис. 4.

Работа выполнена при финансовой поддержке Российского фонда фундаментальных исследований (гранты 94-01-00951-а, 96-02-16252-а). Компьютерная реализация алгоритма на рабочих станциях SUN осуществлена благодаря поддержке Королевского Технологического Института в Стокгольме и гранту Шведской Академии наук.