«Ленивые» вейвлеты эрмитовых кубических сплайнов и алгоритм с расщеплением Текст научной статьи по специальности «Математика»

ВЕСТНИК ТОМСКОГО ГОСУДАРСТВЕННОГО УНИВЕРСИТЕТА

2011

Управление, вычислительная техника и информатика

№ 1(14)

УДК 519.6

Б.М. Шумилов

«ЛЕНИВЫЕ» ВЕЙВЛЕТЫ ЭРМИТОВЫХ КУБИЧЕСКИХ СПЛАЙНОВ И АЛГОРИТМ С РАСЩЕПЛЕНИЕМ

В статье изучен неизвестный ранее тип «ленивых» вейвлетов для эрмитовых кубических сплайнов. Получен алгоритм вейвлет-разложения в виде двух независимых трехдиагональных систем линейных уравнений со строгим диагональным преобладанием.

Ключевые слова: эрмитовы кубические сплайны, вейвлеты, разложение

Вейвлетом называется малая, то есть короткая или быстро затухающая волна, множество сжатий и смещений которой порождает некоторое пространство ограниченных функций на всей числовой оси [1 - 3]. За счет сжатия вейвлеты выявляют с разной степенью подробности различие в характеристиках измеренного сигнала, а путем сдвига способны проанализировать свойства сигнала в разных точках на всем изучаемом интервале. Свойство локальности вейвлетов обеспечивает им известное преимущество при анализе нестационарных сигналов, например, по сравнению с преобразованием Фурье. К недостаткам ортонормальных вейвлетов [1] относится то, что они не имеют аналитического представления и графически похожи на фрактальные кривые. Недостатком сплайн-вейвлетов [2] является то, что для них не существует явных конечных формул разложения. Поэтому при вычислениях используют приближенные соотношения для главных коэффициентов разложения [2] либо решают системы линейных алгебраических уравнений [3]. В данной работе для случая эрмитовых кубических сплайнов [4] рассмотрен неизвестный ранее тип «ленивых» вейвлетов, для которых предложен эффективный алгоритм вейвлет-анализа на основе конечных неявных соотношений разложения.

1. Построение систем «ленивых» эрмитовых сплайн-вейвлетов

Основой для построения вейвлет-преобразования является набор вложенных пространств ... ¥Ь-1с ¥Ь с ¥ь+1 .... В данном случае пространство ¥Ь является пространством эрмитовых кубических сплайнов на отрезке [а, Ь] с равномерной сеткой узлов АЬ : иг = а + (Ь - а) / / 2Ь, / = 0, 1,., 2Ь, Ь > 0, и базисными функциями N^,-,0 (V) = ф3(у - /), МЬ,;1 (у) = у(у - /) V/, V = 2Ь (и - а) / (Ь - а) + 1, с центрами в целых числах, порожденными сжатиями и сдвигами двух функций вида

(1 - г)2(1 + 2/), 0 < г < 1,

г (1 - г)2, г2(1 - г), 0,

0 < г < 1,

1 < г < 2, /г [0,2].

Фз (г) = ] г2(3 - 2г), 0,

.2

1 < г < 2, у(г) = <

г г [0,2];

На любой сетке АЬ, Ь > 0, интерполяционный кубический эрмитов сплайн может быть представлен как

SL (u) = £ CL’° NL0 (u) + X (u), a < u < b, (1)

i=0 i=0

где коэффициенты CLk, k = 0, 1, являются узловыми значениями и соответственно производными аппроксимируемой функции. Если записать базисные сплайн-

функции в виде единой матрицы-строки, qL = |^ N 1, Nf0, N^,..., 1J, и упоря-

дочить коэффициенты сплайна в виде вектора, СL =^С0'1, С^0, С1£'1,..., ^L1 J ,

то уравнение (1) переписывается как SL(u) = ф (u) CL.

Суть вейвлет-преобразования состоит в том, что оно позволяет иерархически разложить заданную функцию на серию все более грубых приближенных представлений и локальных уточняющих подробностей. Пусть «более грубый» уровень представления в VL-1 получается из «более подробного» уровня представления в Vl посредством удаления каждого второго отсчета. Тогда базисными функциями для VL-1 будут функции NL-1ik, в два раза большие по ширине с центрами в четных целых числах. Пространство вейвлетов WL-1 определяется как дополнение VL-1 до VL, так что любая функция в VL может быть записана в виде суммы какой-то функции в VL-1 и какой-то функции в WL-1. В [3] на примере сплайнов первой степени было предложено в качестве базисных функций в WL-1 использовать базисные функции в VL с центрами в нечетных целых числах («ленивые» вейвлеты). Тогда формулы для вейвлет-коэффициентов приобретают локальный вид и для данного случая независимо получены в [5]. А в [6] эти формулы были применены к более подходящей схеме вычисления производных для анализа негладких данных.

Мы предлагаем использовать в качестве вейвлетов для WL-1 функции NLik в VL с центрами в четных целых числах. Поскольку WL-1 должно являться дополнением VL-1 в VL, размерности этих пространств должны удовлетворять соотношению Dim (VL) = Dim (VL-1) + Dim (WL-1). Иногда для этого достаточно ограничиться рассмотрением периодического случая, когда первая и последняя контрольные точки совпадают. Это соответствует аппроксимации замкнутых кривых и поверхностей. Для незамкнутых кривых предлагается вычитать из исходных координат уравнение прямой, соединяющей первую и последнюю точки. Тогда две базисные функции по краям отрезка [a, b] удаляются из базиса, и размерности полученных пространств V°L, W°L-1 равны 2(2L+1) - 2 = 2L+1 и 2L соответственно. Следовательно, Dim (VL) = Dim (V°L-1) + Dim (W°L-1).

Обозначим базисные сплайн-функции и коэффициенты эрмитового кубического сплайна с отсутствующими элементами по концам отрезка аппроксимации как 90l и C0L'. Аналогично, обозначим базисные вейвлет-функции как ML-1i0 = 93(v - 2i), ML-1i;1 = y(v - 2i), i = 0, 1,..., 2l-1, и запишем их в виде

матрицы-строки, ^ = [^М^, М["0, М^,..., mLl 1 J . Соответствующие вейвлет-

коэффициенты на уровне разрешения L будем собирать в вектор, D0L = D^1, DjL,°, D^1,..., ^L’l J . Тогда с использованием обозначения для блочных матриц вейвлет-преобразование может быть записано как [3]

с0і =[ PL I Q0l ]

L-1

С0

Dl-1

0

(2)

Ниже представлен пример матрицы [Р0Ь | QoL], соответствующий Ь = 3:

[ Ро3|бо3 ] =

х

2

3

4

1

0

1

2

_ 3

4

1

2

3

4

1

0

1

2

_ 3

4

1

Здесь и далее пустые позиции представляют собой нулевые элементы. Блоки матрицы РЬ составлены из коэффициентов двухмасштабных соотношений для эрмитовых сплайнов 3-й степени [7]:

Ф3(/) =1 Ф3(2/) + ф3(2/ _ 1) +1 ф3(2/ _ 2) + 3 ((2/) _ у (2/ _ 2)),

У(0 = 2 ¥(2/ _ 1) _ 1 у(2/) _ 8 ¥(2/ _ 2) _ 1 (ф3 (2/) _ Ф3 (2/ _ 2)),

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

Обратный процесс разбиения коэффициентов СЬ на более грубую версию СЬ-1 и уточняющие коэффициенты БЬ-1 состоит в решении системы линейных уравнений (2). При этом для облегчения вычислений матрицу [Р0Ь | Q0L] можно сделать ленточной, изменив порядок неизвестных так, чтобы столбцы матриц Р0Ь и Q0L перемежались. Тем не менее, хотя разрешимость полученной системы и гарантирована линейной независимостью базисных функций, вопрос ее хорошей обусловленности остается открытым.

2. Алгоритм с применением расщепления

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

Теорема 1. Пусть величины ЕЬ = І^’1, І^0, ,..., ] , Ь > 1, получены из

соотношений вида

|Ь’0 = сЬ’0, і = 2,4,...,2Ь -2; ^Ь’1 = С^1, і = 0,2,...,2Ь;

-11 1

1 -10

-9 1

1 -10 1

1 -11

1 -9

£

Ь,0

:Ь,0

:Ь,0

Ь,1

:Ь,1

3

ї,0

X--1 Ь,1

Ь,1

(3)

^к = С1 к, к = 0,1.

Тогда коэффициенты вейвлет-разложения представляют собой результат

матричного умножения

где

С

ї-1

= Я

ЬггЬ

Я1 =

ЯЬ =

4

-4

-2

2

"0

-4

-4

2

2

-48

32 0

-4 -4 0 0 -4 4 0

24 16 0 0 -24 16 0

0 -4 -4 0 ' 0 -4 4

0 24 16 0 ' 0 -24 16

0 48 32 0

24 -16 0

4 4 1 0 4 -4 0

-12 -8 0 1 12 -8 0

0 4 4 1 ' 0 4 -4

0 -12 -8 0 ' 1 12 -8

0 -24 -16 1

(4)

Доказательство. Согласно построению, на каждых трех внутренних шагах сетки А перекрываются по две пары широких базисных функций и вейвлетов и пять пар узких базисных функций. Поэтому после приведения к сетке с единичным шагом на отрезке [-1, 2] можно записать конечное неявное соотношение разложения

2 0 0

X а Фз (2х-1 )= X Ь] Фз (2х-1 -2/)+ X Ь; V(2x-1 -2/) +

г=~2 /=-1 ]=— (5)

0 0 ( )

+ Хс/ Фз (х - /)+ Xс) V (х - /)

/=-1 /=-1

где ф3 (2 х- /) - эрмитовы базисные сплайны на густой сетке, ф3 (х - /), V (х -/) -эрмитовы базисные сплайны на прореженной сетке, ф3 (2х -1- 2/), V (2х -1- 2/) -базисные вейвлеты на прореженной сетке.

Тогда для определения неизвестных коэффициентов согласно табл. 1 имеем соответственно при

1 11 Г 1 Л п 3 ! ( 1

х = —: а 2 = с 1 — + с, I — 1, 0 = с,— + с, I —

2 2 1 2 11 8 ) 1 2 11 4

х = 0: а-1 = Ь_ 1 + с -1, 0 = Ь-1 • 2 + с-1,

1 1 1 ! 1

х = —: а0 = с_ . — + с0 — + с_. — + с01 —

2

0="_■ (_ 2)+с°'3+"_■ (_1)+(_ \

х = 1: а = Ь0 + с0, 0 = Ь • 2 + с0,

3 1 11 ( 3Л 1(1

х = —: а2 = с0 + с0 , 0 = с01—1 + с01 —

2 2 0 2 0 8 0 1 2) 01 4

Решение полученной системы уравнений имеет вид

а-1 = а1 = 0, Ь0 = 4, Ь-1 = 4, с-1 = с0 = -4, Ь1-1 = 12, Ь10 = -12,

с1-1 = -24, с'о = 24, а-2 = а2 = 1, а0 = -10.

Таким образом, из соотношения (5) выделяется неявное трехчленное разложение вида

Ф3 (2х + 2)-10 ф3 (2х) + ф3 (2х - 2) = 4 ф3 (2х +1) + 4 ф3 (2х -1) +12 V(2х +1)--12V(2х-1)-4ф3 (х +1)-4ф3 (х)-24V(х +1) + 24V(х).

Т аблица 1

Значения базисных функций и их производных в точках отрезка [0, 2]

х Ф3(х) ф'3(х) ¥(» УМ Ф3(2х) Ф'3(2х) у(2х) ^(2х)

0 0 0 0 0 0 0 0 0

1/2 1/2 3/2 - 1/8 - 1/4 1 0 0 2

1 1 0 0 1 0 0 0 0

3/2 1/2 - 3/2 1/8 - 1/4 0 0 0 0

2 0 0 0 0 0 0 0 0

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

у (2х + 2)-10у (2х) + у (2х - 2) = -4ф3 (2х +1) + 4 ф3 (2х -1)- 8у (2х +1)-

-8у(2х-1) + 4ф3 (х +1)-4ф3 (х) +16у(х +1) +16у(х).

В нечетных узлах выполняются тождества для базисных сплайн-вейвлетов ф3 (2х ± 1), у (2х± 1).

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

2 0 0

Xаф3(2х-1 ) = й0ф3(2х-1)+ Xь)^(2х-1 -2})+с0ф3(х)+ X с°(х-}); (6)

/=0 ] =-1 ]=-1

2 Л 0 0

X ау(2х-1) = Ь0ф3 (2х-1)+ X ь)у(2х-1 -2]) + с0ф3 (х)+ X с°(х-}). (7)

г=-1 ] =-1 ] =-1

Тогда для отыскания коэффициентов соотношения (6), согласно табл. 1 имеем соответственно при

х = 0: 0 = Ь- • 2 + с-1,

1 1 111 ( 1 | А 3 1 ( 1 | 1 ( 1

х = —: а0 = с0— + с,- + с01 — 1, 0 = с0— + с-11 — 1 + с01 —

2 0 0 2 1 8 0 1 8 ) 0 2 11 4) 01 4

х = 1: а! = Ъ() + с0, 0 = Ь0 • 2 + с0,

х = -: а9 = сп- + с0°, 0 = с01--1 + с01-°|.

л2 — '-'0 ' '-'0 > — с'0 I I ~ с'0 .

2 2 0 2 0 8 0 ^ 2) 0 ^ 4

Решение полученной системы уравнений имеет вид а1 = 0, с1_1 = -48, Ь1^ = Л = 24, Ь10 = -12, с0 = -4, Ь0 = 4, а2 = 1, а0 = -11. Таким образом, из соотношения (6) выделяется неявное двухчленное разложение вида

-11ф3 (2х) + ф3 (2х - 2) = 4 ф3 (2х -1) + 24у (2х +1)-12у (2х -1)-

- 4 ф3 (х)- 48у (х +1) + 24у (х).

Аналогично, для отыскания коэффициентов соотношения (7) имеем х = 0: 2<з-1 = Ь- • 2 + с-,

х - 2: 0 - + С-і8 + 4 [ 8 ] , Ч - [ 4 ) + < [ 4

х - 8: 0 - Ь0 + с0, 2а8 - • 2 + <0,

3 Л - 8 -8 8 - ( 3 А „8 8

х -“ : 0 - с0~ + C0-, 2а2 - С0 1~ 1+ С0 I

2 2 8 - ^ 2) \ 4,

Решение полученной системы уравнений имеет вид

а-1 = а1 = 0, Ь0 = 4, с0 = -4, с;-1 = 32, Ь°1 = -16, Ь0 = -8, с0 = 16, а2 = 1, а0 = -9.

Таким образом, из соотношения (7) выделяется неявное двухчленное разложение вида

-9у(2х)+у(2х - 2) = 4ф3 (2 х-1)-16у(2 х+1)-8у(2 х-1)-4ф3 (х)+32у(х+1)+16у(х).

На крайних справа шагах сетки А0 разложения отражены полученным выше: ф3 (2х+2)-11ф3 (2х)=4ф3 (2 х+1)+12у(2 х+1)-24у(2 х-1)-4ф3 (х+1)-24у(х+1)+48у(х), у(2х+2)-9у(2х)=-4ф3 (2 х+1)-8у(2 х+1)-16у(2 х-1)+4ф3 (х+1)+16у(х+1)+32у(х).

И на нулевом уровне разрешения разложения имеют вид

0 0

аоф3 (2х)= X Ь1 ^(2х -1 - 2І)+ X СІ-^(х - І)

І=-1 І=-1

X аі ^(2 х -1 )= Xь) ^(2 х -1 - 2 І)+ Xг) ^(х - І)

і=-1 І=-1 І=-1

откуда вытекают соотношения

ф3 (2х) = -2у(2х +1) + 2у (2х -1) + 4у (х +1)- 4у (х), у (2х) = 2 у (2 х +1) + 2 у (2 х -1)- 4 у (х +1)- 4 у (х).

Введем матрицы О0, блоки которых составлены из коэффициентов левых частей полученных соотношений:

О1 =

, О0 =

-11 0 0 0 1 0 0

0 -9 0 0 0 1 0

0 0 1 0 0 0 0

0 0 0 1 0 0 0

1 0 0 0 -10 0 0 '

0 1 0 0 0 -10 0 '

0 0 0 1 ' 0 1 0

0 0 0 0 ' 0 0 1

0 1 0 0 ' 0 0 0

0 0 1 0 ' 1 0 0

0 -11 0

0 0 -9 0

1

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

ф0° Оь = [фоі-1 | уоі-1] Я°, 0 > 1.

Отсюда, используя свойство дополнения пространства вейвлетов, находим

К-|

с,

0-1

в0

= ф0 С = [

ф°1

] я° (о0)

Вычислим вспомогательные переменные 2 из системы линейных уравнений Оь 2 = С0Ь. При этом матрица Оь представляет собой хорошо определенную ленточную матрицу со строгим диагональным преобладанием. Более того, после

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

3. Пример

Рассмотрим в качестве тестовой функцию Хартена [6]:

Ї (х) =

-28іп(3лх), х < І,

^іп(4лх)|, 3 < х < р --^єіПЗлх), х > у.

Эта кусочно-гладкая функция имеет разрывы первого рода в точках х=1/3 и 2/3 и угол (разрыв первой производной) в точке х=1/2. Ее производная - также кусочно-гладкая функция.

На рис. 1 представлены результаты реконструкции вейвлет-анализа интерполяционного эрмитового кубического сплайна 53(0 с числом разбиений п = 2Ь на интервале 0 < х < 1 при условии точного вычисления производных. Здесь и далее длина шага Ах = 1/п, где верхний уровень разрешения Ь = 5, то есть п = 32, сплошной линией обозначается исходная функция. На рис. 2 представлены аналогичные результаты реконструкции вейвлет-анализа эрмитового кубического сплайна 53(/) при условии численного вычисления производных. При этом выброс на месте разрыва производной приобретает симметричный вид.

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

Рис. 1. Графики вейвлет-реконструкции эрмитового кубического сплайна при условии точного вычисления производных

Рис. 2. Графики вейвлет-реконструкции эрмитового кубического сплайна при условии численного вычисления производных

Заключение

Представленные в статье схемы построения «ленивых» эрмитовых кубических сплайн-вейвлетов и получения для них неявных соотношений разложения могут быть распространены и на сплайны более высокой степени [8] и предоставляют широкие возможности для создания эффективных параллельных алгоритмов построения и использования эрмитовых сплайн-вейвлетов.

ЛИТЕРАТУРА

1. Добеши И. Десять лекций по вейвлетам: пер. с англ. Ижевск: НИЦ «Регулярная и хаотическая динамика», 2001. 332 с.

2. Чуи Ч. Введение в вейвлеты: пер. с англ. М.: Мир, 2001. 412 с.

3. Столниц Э., ДеРоуз Т., Салезин Д., Вейвлеты в компьютерной графике: пер. с англ. Ижевск: НИЦ «Регулярная и хаотическая динамика», 2002. 272 с.

4. Завьялов Ю.С., Квасов Б.И., Мирошниченко В.Л. Методы сплайн-функций. М.: Наука, 1980. 352 с.

5. Warming R., Beam R. Discrete multiresolution analysis using Hermite interpolation: Bior-thogonal multiwavelets // SIAM J. Sci. Comp. 2000. V. 22:1. P. 269-317.

6. Arandiga F., Baeza A., Donat R. Discrete multiresolution based on hermite interpolation: computing derivatives // Communications in Nonlinear Science and Numerical Simulation. 2004. V. 9. P. 263-273.

7. Heil С., Strang G., Strela V. Approximation by translate of refinable functions // Numer. Math. 1996. V. 73. P. 75-94.

8. Турсунов Д.А., Шумилов Б.М., Эшаров Э.А., Турсунов Э.А. Новый тип эрмитовых мультивейвлетов пятой степени // Пятая Сибирская конференция по параллельным и высокопроизводительным вычислениям: Тез. докл. (1 - 3 декабря 2009 года). - Томск: Изд-во Том. ун-та, 2009, С. 57-58.

Шумилов Борис Михайлович Томский государственный университет E-mail: [email protected]

Поступила в редакцию 26 октября 2010 г.