П3. Вейвлет-преобразование

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

Вейвлет-преобразование, как и преобразование Фурье, состоит в вычислении корреляций между анализируемым временным рядом и базисной функцией преобразования. Так, преобразование Фурье направлено на выявление гармонических составляющих временного ряда. Для этой цели применяется бесконечно-осциллирующая гармоническая функция, которая, по-существу, "накладывается" на анализируемую реализацию процесса. Затем проводится сравнение поведения гармонической функции и изучаемой реализации путем вычисления корреляции. Если в результате сравнения выяснено, что они линейно зависимы, т.е. коррелированы между собой, то это означает, что в составе процесса имеются гармонические составляющие выбранной частоты. Затем частота гармонической функции изменяется и процедура сравнения повторяется. Результатом является спектральная функция, которая отражает исходный процесс из временной области в частотную.

Вейвлет-преобразование временного ряда состоит в разложении ряда по базису, сконструированному из обладающей определенными свойствами функции, называемой вейвлетом^П1, посредством ее масштабных изменений и переносов. Каждая вейвлет-функция базиса характеризуется определенным масштабом (частотой) и локализацией во времени. В отличие от преобразования Фурье вейвлет-преобразование дает двумерную развертку одномерного процесса, при этом частота и время рассматриваются как независимые переменные. В результате появляется возможность анализировать свойства процесса одновременно во временной и частотной областях.

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

Операция огрубления исходной информации осуществляется путем сглаживания исходного ряда с помощью оконной функции -- вейвлета $\psi(t)$. Термином вейвлет обозначают локализованные во временной и частотной областях солитоноподобные функции, обладающие следующими свойствами:

-- нулевым средним $\langle \psi(t)\rangle =\int_{-\infty}^{\infty}\psi (t)dt=0$;

-- ограниченностью: функция $\psi(t)$ быстро убывает при $t\rightarrow \pm \infty$; $\int_{-\infty}^{\infty}\psi(t) dt< \infty$;

-- автомодельностью: при масштабных преобразованиях вейвлета количество осцилляций функции $\psi(t)$ не меняется.

Данные свойства определяют большой класс действительных и комплексных функций, которые являются вейвлетами.

Итак, вейвлет локализован сразу в двух областях -- временной и частотной, и для осуществления вейвлет-преобразования произвольной реализации процесса $x(t)$ необходимо предусмотреть возможность сдвигов вейвлет-функций вдоль временной оси и масштабных преобразований в частотной области путем сжатия или растяжения исходного вейвлета. Такую возможность реализует базисная функция следующего вида:

$$\psi_{ab}(t)=\vert a\vert^{-1/2}\psi \left(\frac{t-b}{a}\right),$$ (П25)

в которой параметры $a$ и $b$ являются действительными числами и определяют масштаб (величину, обратно пропорциональную частоте) и временной сдвиг, соответственно. На основе этой базисной функции вейвлет-преобразование непрерывной функции $x(t)$, определенной на всей временной оси $-\infty <t<\infty$, записывается в виде:

$$W_\psi (a,b)=\vert a\vert^{-1/2}\int_{-\infty}^{\infty}x(t) ... ...t-b}{a}\right) dt=\int_{-\infty}^{\infty}x(t)\psi_{ab}^*(t)dt;$$ (П26)

здесь $W_\psi (a,b)$ -- коэффициенты вейвлет-преобразования; символом "*" отмечена комплексно-сопряженная функция. Параметр $b$ меняется в интервале $-\infty <b<\infty$, пробегая всю временную ось, т.е. всю временную область, на которой определена функция $x(t)$. Параметр $a$ меняется от $-\infty$ до $\infty$, но обычно рассматривается только положительная область. Индекс $\psi$ у коэффициентов $W$ подчеркивает, что значение $W_\psi (a,b)$ зависит от выбора функции $\psi$. При использовании действительного вейвлета получается двумерный массив коэффициентов $W(a,b)$, а при применении комплексного -- двумерные массивы значений модуля $\vert W(a,b)\vert$ и фазы $\Phi(a,b)$:

$$W(a,b) = \vert W(a,b)\vert e^{j\Phi(a,b)} ,$$

$$\vert W(a,b)\vert = \sqrt{ Re [W(a,b)]^2+Im [W(a,b)]^2},$$ (П27)

$$\Phi(a,b) = {\rm arctg} \frac{Im [W(a,b)]}{Re [W(a,b)]}.$$

Процедура вычисления коэффициентов вейвлет-преобразования по формуле (2.26) для каждой пары параметров $a$ и $b$ выглядит следующим образом:

-- растянуть вейвлет $\psi$ в $a$ раз по горизонтали и в $1/a$ раз по вертикали;

-- сдвинуть вейвлет в точку $t^0=t-b$; получаем вейвлет $\psi_{ab}^*$;

-- усреднить значение функции $x(t)$ c помощью временного окна $\psi_{ab}^*$.

Далее процедура повторяется для другой пары параметров $a$ и $b$.

Обратное вейвлет-преобразование записывается с помощью той же базисной функции, что и прямое:

$$x(t)= C_\psi^{-1}\int_{-\infty}^{\infty} \int_{-\infty}^{\infty}W_\psi(a,b) \psi_{ab}(t)\frac{da~db}{a^2};$$ (П28)

здесь  $C_\psi$  -- нормализующий коэффициент,  который определяется фурье-образом $\Psi(f)$ вейвлета:

$$C_\psi =\int_{-\infty}^{\infty} \vert\Psi(f)\vert^2\vert f\vert^{-1} df < \infty.$$ (П29)

Коэффициенты вейвлет-преобразования временного ряда $x(i)$ длины $N$ вычисляются по формуле (шаг дискретизации равен единице):

$$W_\psi (a,b)= \vert a \vert ^{-1/2} \sum_{i=0}^{N-1}x(i)\psi ^*\left(\frac{i-b}{a}\right),$$ (П30)

в которой параметр сдвига $b$ принимает целые значения из области $[0; N-1]$, а параметр масштаба $a$ -- любое значение из области $[2;N]$.

Коэффициенты вейвлет-преобразования содержат информацию об анализируемом процессе и используемом вейвлете. Поэтому выбор анализирующего вейвлета определяется тем, какую информацию необходимо извлечь из процесса. Каждый вейвлет имеет характерные особенности во временной и частотной областях, поэтому иногда с помощью разных вейвлетов можно полнее выявить и подчеркнуть те или иные свойства анализируемого процесса. Если провести аналогию с микроскопом, а вейвлет часто называют математическим микроскопом, то параметр сдвига $b$ фиксирует точку наведения микроскопа, параметр масштаба $a$ -- его увеличение и, наконец, базисный вейвлет $\psi$ определяет оптические свойства микроскопа.

Рассмотрим несколько примеров  вейвлетообразующих  функций. Простейшим является HAAR-вейвлет, названный так по имени предложившего его Хаара, и определяемый соотношением:

$$\psi(t)= \left\{ \begin{array}{cc} 1,&0\le t < 1/2, \\ -1,& 1/2 \le t < 1, \\ 0,& t<0 ,~ t\ge 1. \end{array} \right.$$ (П31)

Этот вейвлет имеет резкие границы во временной области и, как следствие, бесконечно осциллирующие, но убывающие хвосты в частотной области. Кроме того, HAAR-вейвлет является несимметричным. Часто применяется очень похожий, но симметричный FHAT-вейвлет (French HAT -- французская шляпа):

$$\psi(t)= \left\{ \begin{array}{cc} 1,& \vert t\vert \le 1/3, \\ -1/2,& 1/3 < \vert t\vert \le 1, \\ 0,& \vert t\vert >1. \end{array} \right.$$ (П32)

Фурье-образ $\Psi(f)$ этого вейвлета имеет вид:

$$\Psi (f)= 3 \Theta(f) \left( \frac{\sin (f)}{f }-\frac{\sin (f/3)}{3f} \right) ;$$ (П33)

здесь $\Theta(f)$ -- функция Хевисайда. Если правые чпсти формул (2.32) и (2.33) поменять местами, то получим вейвлет Литтвуда-Пели (LP-вейвлет), который в противоположность FHAT-вейвлету имеет резкие границы в частотной области и медленно спадает во временной области. FHAT- и LP-вейвлеты можно считать предельными случаями, между которыми находятся практически все вейвлет-функции.

Наиболее часто среди вещественных вейвлетов используются вейвлеты, сконструированные на основе производных функций Гаусса:

$$\psi_m(t) = (-1)^m \partial^m_t \left[ \exp \left( \frac{-t^2}{2}\right) \right] ,$$

$$\Psi_m (f) = m (j f)^m \exp \left( \frac{-f^2}{2}\right);$$ (П34)

здесь $\partial^m_t[...] = \partial^m[...]/ \partial t^m, m\ge 1$. Более высокие производные позволяют извлечь информацию об особенностях высокого порядка, содержащихся во временном ряде. Вейвлеты на основе функций Гаусса используются для анализа масштабных свойств временных рядов, в том числе мультифрактальных: типа особенности (степенная, импульсная, ступенька и т. п.), ее интенсивности, распределения особенностей по масштабам и т. д. На рис. 19,  а и рис. 19,  б показаны вейвлеты и их фурье-образы, полученные для $m=1$ и $m=2$. По внешнему виду первый из них называют обычно WAVE-вейвлетом (WAVE-волна), второй -- мексиканской шляпой, или MHAT-вейвлетом (Mexican HAT).

Среди комплексных вейвлетов наибольшее распространение получил вейвлет Морле:

$$\psi(t) = \exp(j2 \pi k_0 t) \exp \left[\frac{-t^2}{2}\right],$$

$$\Psi(f) = \Theta(f) \exp\left( \frac{ -(f-k_0)^2}{2}\right) ,$$ (П35)

который представляет собой плоскую волну, промодулированную гауссианом единичной ширины (здесь $k_0$ -- параметр). На рис. 19,  в вейвлет Морле и его преобразование Фурье показаны для $k_0=6$. С увеличением $k_0$ растет частотная избирательность базиса, но ухудшается временная. Данный вейвлет используется, в основном, для частотно-временного анализа процессов, у которых спектральный состав меняется во времени.

: Рис. 19. Временное и спектральное представления WAVE-вейвлета ( а), MHAT-вейвлета ( б) и вейвлета Морле ( в).

Результаты вейвлет-преобразований во многом зависят от используемого вейвлета, однако вейвлет-преобразования обладают некоторыми универсальными свойствами, которые не зависят от вида вейвлета. Рассмотрим некоторые из этих свойств. Будем использовать обозначение $W(a,b)=W[x(t)]$.

1. Линейность:

$$W[\alpha x_1(t)+\beta x_2(t)]=\alpha W_1(a,b)+\beta W_2(a,b).$$ (П36)

2. Инвариантность относительно сдвига:

$$W[x(t-b_0)]=W(a,b-b_0).$$ (П37)

Из этого свойства следует коммутативность дифференцирования, в частности, $\partial_t W[x(t)]=W[\partial_t x(t)]$. Поскольку коэффициенты вейвлет-преобразования являются гладкими функциями параметров даже для негладкой функции $x(t)$, то это свойство можно использовать для численного дифференцирования негладких функций с разрывами производной.

3. Инвариантность относительно растяжения (сжатия):

$$W\left[x\left(\frac{t}{a_0} \right)\right]=\frac{1}{a_0}W\left( \frac{a}{a_0},\frac{b}{a_0} \right).$$ (П38)

Это свойство позволяет выявлять особенности ряда и их характер путем анализа поведения коэффициентов вейвлет-преобразования при $a\rightarrow 0$ вместо анализа функции $x(t)$.

4. Частотно-временная локализация и наличие частотно-временного окна. Как известно, чем лучше функция сконцентрирована во времени, тем больше она "размазана" в частотной области. При переходе от одного масштаба к другому произведение временного и частотного разрешения (площадь частотно-временного окна) остается постоянным. Данное ограничение учитывается вейвлетным базисом, который извлекает высокочастотную информацию из относительно малых временных интервалов и наоборот, низкочастотную спектральную информацию из относительно широких интервалов. В результате вейвлет-преобразование позволяет проследить изменение частотного состава временного ряда с течением времени.

5. Дифференцирование:

$$W[\partial^m_tx(t)]= (-1)^m \int_{-\infty}^{\infty} x(t) \partial^m_t[\psi^*_{ab}(t)] dt.$$ (П39)

Это свойство означает, что вместо дифференцирования исходного временного ряда, например, для удаления полиномиального тренда, можно продифференцировать нужное число раз вейвлет.

6. Энергетическое свойство:

$$\int_{-\infty}^{\infty}x^2(t) dt = C_\psi^{-1} \int_{-\infty}^{\infty}\int_{-\infty}^{\infty}W(a,b)W^*(a,b)\frac {da~db}{a^2}=$$

$$~ = \int_{-\infty}^{\infty}\int_{-\infty}^{\infty} \vert W(a,b)\vert ^2\frac {da~db}{a^2},$$ (П40)

из которого следует, что энергия процесса на каком-либо масштабе может быть рассчитана через коэффициенты вейвлет-преобразования; $\vert W(a,b)\vert ^2$ можно интерпретировать как плотность энергии процесса в частотно-временном пространстве. Если провести усреднение вейвлет-коэффициентов по времени, то получим распределение энергии по масштабам или так называемый глобальный вейвлет-спектр:

$$E_W(a)=\int_{-\infty}^{\infty}\vert W(a,b)\vert ^2 db.$$ (П41)

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

=.45     =.45

=.45

=.45     =.5

=.9

=.9

Рис. 20. Временной ряд $x(i)$ с меняющейся во времени частотой ( а), его спектр мощности, рассчитанный с помощью преобразование Фурье ( б), глобальные вейвлет-спектры ( в), рассчитанные с помощью MHAT-вейвлета (сплошная линия) и вейвлета Морле (пунктирная линия), скелетоны коэффициентов вейвлет-преобразования, рассчитанные с помощью вейвлета Морле ( г) и MHAT-вейвлета ( д), частотно-временные вейвлет-спектры, полученные с помощью вейвлета Морле ( е) и MHAT-вейвлета ( ж). Более светлый цвет на рисунках ( е) и ( ж) соответствует б $\acute {\rm о}$льшим значениям коэффициентов вейвлет-преобразования.

Способы графического представления результатов вейвлет-преобразования могут быть самыми различными. Вейвлет-спектр $W_\psi (a,b)$ представляет собой поверхность в трехмерном пространстве. Вместо изображения поверхности часто используется, во-первых, проекция на плоскость с изображением изоуровней, которые позволяют проследить изменения интенсивностей коэффициентов вейвлет-преобразования на разных временных масштабах, и, во-вторых, картины локальных экстремумов поверхностей, так называемый скелетон, четко выделяющий структуру анализируемого процесса.

Рассмотрим возможности вейвлет-преобразования в выявлении характерных особенностей процессов. В качестве реализации СП рассмотрим нестационарный временной ряд $x(i)=\sin (\omega _i i)$, частота $\omega _i$ которого сначала линейно растет, а потом скачком выходит на постоянное значение. На рис. 20 представлен ряд $x(i)$, его спектр мощности $S_{xx}(f)$ ^П1, рассчитанный через преобразование Фурье, глобальные вейвлет-спектры $E_W(f)$ и скелетоны $W_m(f,b)$, полученные с помощью MHAT-вейвлета и вейвлета Морле, а также частотно-временные спектры $\vert W(f,b)\vert$. Видно (рис. 20,  а), что спектр мощности, вычисленный через преобразование Фурье, дает неверную информацию о частотном составе реализации: энергия высокочастотных компонент значительно превышает истинную. Глобальные вейвлет-спектры (рис. 20,  б) показывают правильное соотношение между энергиями низко- и высокочастотных составляющих ряда. Скелетон (рис. 20,  г) и коэффициенты вейвлет-преобразования (рис. 20,  е), рассчитанные с помощью вейвлета Морле, дают локализацию гармонических составляющих ряда во времени. Изображения коэффициентов $\vert W(f,b)\vert$ и их скелетона, полученные с помощью MHAT-вейвлета (рис. 20,  д, ж), более непривычные, но отражают ту же информацию. MHAT-вейвлет позволяет также оценить характер "скачка" частоты, реагируя на изменение формы реализации, а не ее частотного состава.