Реферат: Матричные операции в вейвлетном базисе

Курсовая работа студентки 3 курса ГромовойМарии Михайловны

Белорусский государственный университет

Факультет прикладной математики иинформатики

Кафедра математической физики

Минск 2003

Введение

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

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

Базовая идея восходит квременам 200-летней давности и принадлежит Фурье: аппроксимировать сложнуюфункцию взвешенной суммой простых функций, каждая из которых, в свою очередь,получается из одной функции-прототипа. Эта функция-прототип выполняет рольстроительного блока, а искомая аппроксимация получается комбинированиемодинаковых по структуре блоков. При этом, если «хорошая»аппроксимация получается при использовании небольшого числа блоков, то темсамым достигается значительное уплотнение информации. В качестве таких блоковФурье использовал синусоиды с различными периодами.

Что прежде всего отличаетвейвлет-анализ от анализа Фурье? Основным недостатком Фурье-преобразованияявляется его «глобальная» чувствительность к «локальным»скачкам и пикам функции. При этом модификация коэффициентов Фурье (например,обрезание высоких гармоник с целью фильтрации шума) вносит одинаковые измененияв поведение сигнала на всей области определения. Это особенность оказывается полезнойдля стационарных сигналов, свойства которых в целом мало меняются со временем.

При исследовании женестационарных сигналов требуется использование некоторых локализованных вовремени компактных волн, коэффициенты разложения по которым сохраняют информациюо дрейфе параметров аппроксимируемой функции. Первые попытки построения такихсистем функций сводились к сегментированию сигнала на фрагменты(«окна») с применением разложения Фурье для этих фрагментов.Соответствующее преобразование — оконное преобразование Фурье — было предложенов 1946-47 годах Jean Ville и, независимо, Dennis Gabor. В 1950-70-х годахразными авторами было опубликовано много модификаций времени-частотныхпредставлений сигналов.

В конце 70-хинженер-геофизик Морли (Jean Morlet) столкнулся с проблемой анализа сигналов,которые характеризовались высокочастотной компонентой в течение короткогопромежутка времени и низкочастотными колебаниями при рассмотрении большихвременных масштабов. Оконные преобразования позволяли проанализировать либовысокие частоты в коротком окне времени, либо низкочастотную компоненту, но необа колебания одновременно. В результате был предложен подход, в котором дляразличных диапазонов частот использовались временные окна различнойдлительности. Оконные функции получались в результате растяжения-сжатия исмещения по времени гаусиана. Morlet назвал эти базисные функции вейвлетами(wavelets) — компактными волнами. В дальнейшем благодаря работам Мейера (YvesMeyer), Добеши (Ingrid Daubechies), Койфмана (Ronald Coifman), Маллы (StephaneMallat) и других теория вейвлетов приобрела свое современное состояние.

Среди российских ученых,работавших в области теории вейвлетов, необходимо отметить С.Б. Стечкина, И.Я.Новикова, В.И. Бердышева.

1. Многомасштабный анализ и вейвлеты

Определение 1.Многомасштабный анализ (multiresolutional analysis) – разложение гильбертовапространства L2(Rd), d³1, в последовательностьзамкнутых подпространств

  />, (1.1)

обладающих следующимисвойствами:

1. />, и /> полно в L2(Rd),  

2. Для любого fÎ L2(Rd),для любого jÎ Z, f(x)ÎVj тогда и только тогда, когда

f(2x) ÎVj-1,

3. Для любого fÎ L2(Rd),для любого kÎ Zd, f(x)ÎV0тогда и толькотогда, когда f(x-k)ÎV0,

4. Существует масштабирующая(scaling) функция jÎV0, что {j(x-k)}kÎZdобразует 

базис Ритца в V0.

Для ортонормальныхбазисов можно переписать свойство 4 в виде:

4’. Существуетмасштабирующая функция jÎV0, что {j(x-k)}kÎZdобразует ортонормальный базис в V0.

Определим подпространствоWj как ортогональное дополнение к Vj в Vj-1,

   />,  (1.2)

и представим пространствоL2(Rd) в виде прямой суммы

   />  (1.3)

Выбирая масштаб n, можемзаменить последовательность (1.1) следующей последовательностью:

  />  (1.4)

и получить

   />  (1.5)

Если имеем конечное числомасштабов, то, не нарушая общности, можно положить j=0 и рассматривать

  />, V0Î L2(Rd) (1.6)

вместо (1.4). В числовойреализации подпространство V0конечномерно.

Функция j — такназываемая масштабирующая (скейлинг-) функция. С ее помощью можно определитьфункцию y — вейвлет — такую, что набор {y(x-k)}kÎZ образуетортонормальный базис в W0. Тогда

   />,m=0..M-1.  (1.7)

Из свойства 4’непосредственно следует, что, во-первых, функция j может быть представлена ввиде линейной комбинации базисных функций пространства V-1. Так какфункции {jj,k(x)=2-j/2j(2-jx-k)}kÎZобразуют ортонормальный базис в Vj, то имеем

   />. (1.8)

Вообще говоря, сумма ввыражении (1.8) не обязана быть конечной. Можно переписать (1.8) в виде

   />, (1.9)

где

   />, (1.10)

а 2p-периодическаяфункция m0определяется следующим образом:

   />. (1.11)

Во-вторых,ортогональность {j(x-k)}kÎZ подразумевает, что

  /> (1.12)

и значит

  />  (1.13)

и   />.  (1.14)

Используя (1.9), получаем

  />  (1.15)

и, рассматривая сумму в(1.15) по четным и нечетным индексам, имеем

 />. (1.16)

Используя 2p-периодичностьфункции m0и (1.14), после замены x/2 на x, получаем необходимоеусловие

   />  (1.17)

для коэффициентов hkв (1.11). Заметив, что

   />  (1.18)

и определив функцию yследующим образом:

  />, (1.19)

где

  />, k=0,…,L-1,  (1.20)

или преобразование Фурьедля y

  />,  (1.21)

где

   />,  (1.22)

можно показать, что прикаждом фиксированном масштабе jÎZ вейвлеты

{yj,k(x)=2-j/2y(2-jx-k)}kÎZобразуют ортонормальный базис пространства Wj.

Равенство (1.17)определяет пару квадратурных зеркальных фильтров (quadrature mirror filters,QMF) H и G, где /> и />. Коэффициенты QMF H и G вычисляютсяс помощью решения системы алгебраических уравнений. Число L коэффициентовфильтра в (1.11) и (1.22) связано с числом исчезающих моментов М, и всегдачетно.

Выбранный фильтр Нполностью определяет функции j и y и, таким образом, многомасштабный анализ.Кроме того, в правильно построенных алгоритмах значения функций j и y почтиникогда не вычисляются. Благодаря рекурсивному определению вейвлетного базиса,все операции проводятся с квадратурными зеркальными фильтрами H и G, даже еслив них используются величины, связаные с j и y.

2. Быстрое вейвлет-преобразование

После того, как вычисленыкоэффициенты hk и gk, т.е. выбран определенный вейвлет,можно проводить вейвлет-преобразование сигнала f(x), поскольку заданортонормальный базис (yj,k, jj,k). Любая функцияf(x)ÎL2(R) полностью характеризуется ее вейвлет-коэффициентамиразложения по этому базису и потому может быть представлена формулой

  />. (2.1)

Зададим все пределысуммирования в формуле (2.1). Функцию f(x) можно рассматривать на любом n-муровне разрешения jn. Тогда разделение между ее усредненнымизначениями на этом уровне и флуктуациями вокруг них выглядят как

  />. (2.4)

На бесконечном интервалепервая сумма может быть опущена, и в результате получается «чистое»вейвлет-разложение.

Коэффициенты sj,k иdj,k содердат информацию о составе сигнала на разных масштабах ивычисляются по формулам:

   />,  (2.2)

   />.  (2.3)

Однако при этомкомпьютерные расчеты занимают довольно длительное время, т.к. при вычисленииприходится проводить O(N2) операций, где N – число имеющихсязначений функции. Опишем более быстрый алгоритм.

В реальных ситуациях социфрованным сигналом мы всегда имеем дело с конечным набором цифр (точек).Поэтому всегда существует наилучший уровень разрешения, когда каждый интервалсодержит по одному числу. Соответственно и суммирование по k будет идти вконечных пределах. Удобно изменить шкалу разрешения (или шкалу f), приписавзначение j=0 этому наилучшему уровню разрешения. В этом случае легко вычислитьвейвлет-коэффициенты для более усредненных уровней j³1. Многомасштабныйанализ приводит естественным путем к иерархической и быстрой схеме вычислениявейвлет-коэффициентов заданной функции.

В общем случаеитерационные формулы быстрого вейвлет-преобразования имеют вид:

   />,  (2.4)

  />  (2.5)

с

  />.  (2.6)

Эти уравненияобеспечивают быстрые (или пирамидальные) алгоритмы вычислениявейвлет-коэффициентов, поскольку требуют только O(N) операций для своегозавершения. Начав с s0,k, мы вычислим все другие вейвлет-коэффициенты,если параметры вейвлета hm и gm известны. Явный видвейвлета при этом не используется. Простая форма полученных итерационныхуравнений служит единственным оправданием введения множителя /> в функциональноеуравнение (1.8). В принципе, коэффициенты hm и gm можнобыло бы перенормировать. Однако, уравнения (2.4), (2.5) используются напрактике значительно чаще других, и поэтому эту нормировку не изменяют. Любыедополнительные сомножители в них могут привести лишь к усложнению численныхрасчетов.

Остающиеся проблемысвязаны с начальными данными. Если известен явный вид функции f(x), токоэффициенты s0,k можно вычислить, используя формулу (2.6). Носитуация отличается от этой, если доступны только дискретные значения f(x).Чтобы достичь высокой точности, хорошо бы задать очень малые интервалы (плотнуюрешетку), но это зачастую недоступно из-за конечности интервалов сбораинформации. В таком случае простейшее принимаемое решение состоит внепосредственном использовании величин f(k) из доступного набора данных в видекоэффициентов s0,k и применении быстрого вейвлет-преобразования сиспользованием формул (2.4), (2.5). Это безопасная операция, т.к. пирамидальныйалгоритм обеспечивает полную реконструкцию сигнала, а коэффициенты s0,kпо сути представляют собой локальные средние значения сигнала, взвешенные соскейлинг-функцией.

В общем случае можновыбрать

  />.  (2.7)

Рассмотренная ситуацияотвечает условию s0,k=f(k), что соответствует cm=d0m.

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

3. Двумерные вейвлеты

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

Тривиальный путьпостроения двумерного ортонормального базиса исходя из одномерногоортонормального вейвлет-базиса yj,k(x)=2j/2y(2jx-k)состоит в том, чтобы путем тензорного произведения образовать соответствующиефункции из двух одномерных базисов:

  />.  (3.1)

В этом базисе двепеременных x1 и x2 сжимаются по-разному.

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

  />,j,k,lÎZ,  (3.2)

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

/>, />, />.

Тогда двумерные вейвлетызапишутся в виде

/>, />, />.

На двумерной плоскостипроисходит анализ по горизонталям, вертикалям и диагоналям с одинаковымразрешением в соответствии с тремя выписанными выше вейвлетами.

4. Матричные операции

4.1 Матричное умножение

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

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

  /> при />, (4.1.1)

  /> при />, (4.1.2)

Топология распределенияэтих матричных элементов внутри матрицы может оказаться весьма запутанной.

Рассметрим действиеоператора Т на функцию f, которое превращает ее в функцию g.

  />  (4.1.3)

Как g, так и f могут бытьпредставлены в виде вейвлет-рядов с вейвлет-коэффициентами (f sj,k;fdj,k) и (g sj,k;g dj,k).На наиболее детальном уровне разрешения jn отличны от нуля толькоs-коэффициенты, и преобразование имеет вид

   />.  (4.1.4)

На следующем уровнеполучаем

 />, (4.1.5)

 />, (4.1.6)

где

/>

и замена нижних индексовS®D соответствует подстановке j®y под знаком интеграла.

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

На конечном этапе мыимеем дело с вейвлет-представлением, описываемым формулой (2.1), в которой ввекторах остается только один s-коэффициент, представляющий взвешенное среднеефункции по всему интервалу ее задания, а SS-переход от f к g описываетсяверхним левым квадратиком на этом рисунке. В то же время на пути к этой формулеот скейлинг-представления нам приходилось иметь дело со средними величинами напромежуточных уровнях, разлагая их затем на каждом этапе на части, s и d,последующих уровней разрешения. Эти промежуточные s-коэффициенты были опущены,потому что мы заменяли их на s- и d-коэффициенты поледующих уровней. Именнопоэтому окончательная матрица при стандартном подходе приобретает такой сложныйвид.

/>

Рис.1. Матричноепредставление при стандартном подходе к вейвлет-анализу.

 Части матрицы сненулевыми вейвлет-коэффициентами заштрихованы.

/>

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

Рис.2. Нестандартноематричное умножение при вейвлет-анализе.

Различные уровниоказались полностью развязанными, потому что в матрице теперь полностьюотсутствуют блоки, которые ранее перепутывали их. Блок с SS-элементамиизвлечен, а на его место вставлена нулевая матрица. Полная матрицасоответстваенно искусственным образом увеличилась. Вместе с ней увеличились ивекторы, характеризующие функции f и g. Теперь здесь удерживаются всепромежуточные s-коэффициенты вейвлет-разложения функции f. Каждый блок Sj+1получается из Sj и Dj. В матрице преобразования равнынулю все SS-элементы за исключением их величин на низшем уровне S0S0.Все остальные SD-, DS-, DD-матрицы почти диагональны вследствие конечностиобласти задания вейвлетов и скейлинг функций. Приведенная на рис. 2 формафункции g преобразуется в ее обычное вейвлет-представление из рис. 1 путемразделения каждого Sj в Sj-1 и Dj-1стандартным методом. Затем эти Sj-1 и Dj-1 добавляются всоответствующие компоненты вектора. Эта процедура итерируется, начиная теперьуже с Sj-1, вполоть до S0, когда мы приходим к обычномувейвлет-представлению функции g. Таким способом мы избавляемся от всехs-коэффициентов за исключением s0. Вычисления можно теперь проделатьочень быстро.

4.2 Обращение матрицы

Утверждение 1.Последовательность матриц Xk такова, что

   Xk+1=2Xk-XkАXk,   (4.2.1)

   X0=aА*,  (4.2.2)

где А* — сопряженная матрица и a выбирается таким образом, чтобы наибольшее собственноезначение матрицы aА*А меньше двух. Тогда последовательность сходитсяк обобщенной обратной матрице А-1.

Если это утверждениескомбинировать с алгоритмом быстрого матричного умножения, то получаетсяалгоритм для построения обратной матрицы в стандартной форме с трудоемкостью /> и внестандартной форме с трудоемкостью />, где R – число обусловленностиматрицы. С помощью числа R можно оценить соотношение между наибольшим инаименьшим сингулярными числами выше порога точности.

4.3 Вычислениеэкспоненты, синуса и косинуса от матрицы.

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

Алгоритм вычисленияэкспоненты матрицы основывается на тождестве

  />.  (4.3.1)

Во-первых, exp(2-LA)может быть посчитана, например, с помощью ряда Тейлора. Число L выбираетсятаким образом, чтобы наибольшее сингулярное число матрицы 2-LA быломеньше единицы. На втором шаге алгоритма для достижения результата матрица 2-LAвозводится в квадрат L раз.

Аналогично, синус икосинус от матрицы могут быть посчитаны с исподьзованием формул двойного угла.

   />  (4.3.2)

   />,  (4.3.3)

при l=0,…,L-1

   />  (4.3.4)

   />,  (4.3.5)

где I – тождество. Сновавыбираем L таким образом, чтобы наибольшее сингулярное число матрицы 2-LAбыло меньше единицы, вычисляем синус и косинус матрицы 2-LA, спомощью рядов Тейлора, а затем используем формулы (4.3.4) и (4.3.5).

Обычно такие алгоритмытребуют по меньшей мере O(N3) операций, так как должне бытьвыполнено достаточно много операций по умножению густых матриц. Быстрыйалгоритм для умножения матриц в стандартной форме уменьшает сложность до неболее чем /> операций,а быстрый алгоритм для умножения матриц в нестандартной форме – до O(N)операций.

Списоклитературы

Beylkin G.Wavelets and Fast Numerical Algorithms.

Beylkin G.Wavelets, Multiresolution Analysis and Fast Numerical Algorithms.

Дремин И.М., Иванов О.В.,Нечитайло В.А. Вейвлеты и их использование // Успехи физических наук – 2001,№5. – С.465-500