Представлена конструкция мультивейвлетов на основе эрмитовых сплайнов произвольной нечетной степени на конечном отрезке. Мультивейвлеты ортогональны многочленам той же степени и имеют наименьший носитель. Предложен алгоритм мультивейвлет преобразования с применением метода блочной матричной прогонки. Исследована устойчивость на больших сетках. Дан численный пример аппроксимации и сжатия данных для случая мультивейвлетов седьмой степени.
Multiwavelets of odd degree, orthogonal to polynomials informatics and programming.pdf Вейвлетом называется короткая или быстро затухающая волновая функция (всплеск), множество сжатий и смещений которой порождает пространство измеримых функций на всей числовой оси [1-3]. Основой для построения вейвлетов является наличие набора аппроксимирующих пространств ... Vl-1 с Vl с Vl+1 ..., таких, что каждая базисная функция в V- может быть выражена в виде линейной комбинации базисных функций в VL. В частности, таким свойством обладают сплайны - гладкие функции, склеенные из кусков многочленов степени m, на вложенной последовательности сеток. Если порядок склейки равен m - 1, то классические полуортогональные вейвлеты (элементы пространства Vl, ортогональные пространству Vl-1) имеют довольно большой носитель из 2m + 1 шагов сетки. В случае аппроксимации на конечном интервале ситуация усугубляется, так как приходится учитывать краевые эффекты. Это препятствует их широкому использованию для решения задач математического моделирования. В отличие от этого эрмитовы сплайны нечетной степени m = 2r + 1 (соответствующие склейке порядка r) приводят к вейвлетам с носителем из трех шагов сетки, что, несомненно, предпочтительнее. Поскольку в базисе таких функций несколько, они называются мультивейвлетами [4-8]. В [9, 10] были предложены методы построения вейвлетов, ортогональных многочленам, с уменьшенными носителями. Идея уменьшения носителя вейвлета за счет замены свойства ортогональности пространству сплайнов на прореженной сетке ортогональностью многочленам представляется привлекательной. Действительно, с точки зрения скорости приближения гладких функций [3] данные типы вейвлетов эквивалентны, а ортогональность многочленам обеспечивает локально максимальную «похожесть» на наилучшее среднеквадратическое приближение. В данной работе мы излагаем результаты исследования мультивейвлетов произвольной нечетной степени с носителем [0, 2], ортогональных многочленам той же степени, на конечном отрезке. 1. Построение системы базисных мультивейвлетов произвольной нечетной степени на конечном отрезке Пусть на отрезке [a, b] задана вложенная последовательность равномерных сеток AL: xi = a + i h, i = 0, 1, ..., 2L, h = (b - a)/2L, L > 0. Мультивейвлеты возникают тогда, когда с каждым узлом сетки связано более одной базисной функции [4]. В частности, если базисные функции NLiXx) = - i), k = 0, 1, ..., r Vi, где v = (x - a)/h, с центрами в целых числах, порождены сжатиями и сдвигами r + 1 функций вида [11. С. 82]: (-1)к юк (-t) при -1 < t < 0, юк (t) при 0 < t < 1, 0, t е[-1,1], где сок (V) = (l - tY+ Z -- //l+''', /с = 0,1.....г, то при условии отсечения выступающих за концы отр=о &!|3!г! резка половинок функций ф0(0, . • •, ф?-(0 полученное пространство Vl является пространством эрмитовых сплайнов степени 2r + 1 гладкости С. На любой сетке AL, L > 0, интерполяционный эрмитов сплайн (2r + 1)-й степени может быть представлен как г 2l SL (X) = Z Z CL,k nL", (X), a < x < b, (1) к=0 i=0 где коэффициенты CL,k, k = 0, ..., r, являются значениями и соответствующими производными аппроксимируемой функции f (х), помноженными на h в соответствующей степени: {fk)(rh)h, к = 0, ..., r, i = 0, ..., 2l}, в узлах сетки. Поскольку сетка AL-1, L > 1, получена из AL посредством удаления каждого второго узла, то соответствующее пространство VL-1 с базисными функциями NL-1i,k, с носителями, в два раза большими по ширине, и центрами в четных целых числах, вложено в Vl. Остаток Wl-1 от разности пространств Vl и Vl-1, размерности (r + 1)(2L + 1 - 2L-1 - 1) = (r + 1)2L-1, называется пространством вейвлетов. В классической теории вейвлетов базисные функции Wl-1 ортогональны всем базисным сплайнам на прореженной сетке Vl-1 по отношению к некоторому скалярному произведению. Это означает, что пространство Vl представляет прямую сумму Vl-1 и Wl-1 : Vl = Vl-1 © Wl-1 . В отличие от этого, будем искать вейвлеты MLi,k(x), к = 0, 1, ..., r Vi, как линейные комбинации базисных эрмитовых сплайнов на сетке AL+1, удовлетворяющие условиям ортогональности многочленам (2r + 1)-го порядка [12] ь {x)xmdx = 0, k = 0,,...,rVi(m = 0,,...,2r+V) (2) а и имеющие минимальные из возможных носителей. Теорема. Система функций MLi,k(x), к = 0, 1, ..., r, i = 1, 2, ..., 2L, удовлетворяющих условиям (2) с носителями не более чем из двух шагов сетки AL, существует и образует базис в пространстве Wl. Доказательство. а) Пусть L > 1, и носители вейвлетов [x2i-1, X2i+p+1], p > 1, полностью располагаются внутри отрезка аппроксимации [a, Ъ]. Приp > 1 для всех i = 1, 2, ..., 2L - 1 - [p/2] имеем разложение ML" (X) = Z Z a^i (t-Д t = (x - x2 i)/h, -1 < t < p +1. (3) l=0 j=0 Разложение (3) подставим в (2) и вычислим все необходимые интегралы, учитывая при этом, что подынтегральные выражения обращаются в нуль вне носителей вейвлетов. Условия ортогональности (2) определяют однородную систему 2r + 2 уравнений относительно коэффициентов вейвлета MLk(x). В силу линейной независимости на интервалах [x2i-1, X2i+p+1] базисных функций N^ k(x), l = 0, ., r, j = 0, 1, ., p, и поверочных функций xm, m = 0, 1, ., 2r + 1 [1], ранг полученной системы равняется min[2r + 2, (r + 1)(p + 1)]. Если предположить, что носитель вейвлета равен трем шагам сетки AL+1, т.е. p = 1, то однородная система становится квадратной и, будучи невырожденной, может иметь только тривиальное решение. Поэтому будем далее считать, что носитель вейвлета равен четырем шагам сетки AL+1, те вейвлет построен из 3r + 3 базисных сплайнов. В этом случае ранг матрицы равен 2r + 2. Поэтому однородная система 2r + 2 уравнений с 3r + 3 неизвестными имеет r + 1 линейно независимых частных решений. Длина носителей, полученных для каждого номера i вейвлетов, равна длине носителей базисных сплайнов на сетке AL. Их линейная независимость вытекает из того факта, что они представляют собой множество сдвигов ненулевых функций с компактными носителями. При этом количество вейвлетов, полностью лежащих внутри отрезка аппроксимации, равно (r + 1)(2L - 2) = (r + 1)2L - 2(r + 1), что на 2(r + 1) меньше разности между размерностями пространств Vl+i и Vl. (r + 1)(2L+1 + 1) - (r + 1)(2l + 1) = (r + 1)2L. Фк (t) = Вблизи концов отрезка интервалы интегрирования не выходят за пределы отрезка. Поэтому с учетом того, что по краям отрезка [a, b] от выступающих за концы отрезка функций фо(0, • ••, Фг(0 остается по половинке, граничные вейвлеты отличаются от вейвлетов, предложенных выше, при условии ортогональности многочленам (2r + 1)-й степени. Аналогичные рассуждения приводят к двум системам линейно независимых функций wo(t), •.., wr(t) отдельно на левом и на правом концах, дополняющих до базиса построенную выше вейвлет-систему. б) Пусть L = 0. Вычисления снова дают r + 1 линейно независимых частных решений, которые дают базис в пространстве V - V0 с размерностью (r + 1) 3 - (r + 1) 2 = (r + 1). Теорема доказана. В соответствии с правилами линейной алгебры сформулируем следующий алгоритм построения эрмитовых мультивейвлетов нечетной степени. Следствие. Введем в рассмотрение три матрицы размерности (2г + 2) х (г + 1) с элементами j+1 | f^{2t-j)tmdt, 7=0,1,2,/ = 0,1,..., г, m = 0,l,...,2r + l j-1 и составим из них блочную квадратную матрицу (2r + 2)-го порядка вида [Л0 | Л1]. Тогда при условии, что коэффициенты разложения (3) а! = {1, l = к; 0, l Ф к}, каждый из r + 1 столбцов матрицы дшпег у^дшпег = 1/' дает значения коэффициентов a1-, j = 0,2, / = 0,1,.. .,г, соответствую щего k-го базисного вейвлета, полностью лежащего внутри отрезка аппроксимации. На левом конце отрезка элементы матрицы R0 вычисляются по укороченному интервалу: 1 (2t)tmdt, I = 0,1,..., г, т = 0,1,..., 2r +1, а. 0 а коэффициенты разложения (3) а'-, j - 1,2, / = 0,1,...,г, соответствующего к-то базисного вейвлета даются значениями столбцов матрицы Га1^/Af ] - -[R11R2 ]-1 R0 при условии, что коэффициенты - {1, l = к; 0, l Ф к}. Аналогично на правом конце отрезка аппроксимации элементы матрицы R2 вычисляются по укороченному интервалу: 2 Rl,i = IФ/ (2^ - 2)tmdt, I = 0,1,..., г, m = 0,1,..., 2r +1, 1 а коэффициенты разложения (3) а', j - 0,1, / = 0,1,.. .,г, соответствующего к-то базисного вейвлета даются значениями столбцов матрицы Aq®111/A"®111 - -[R0 | R1]-1 R2 при условии, что коэффициенты а2 = {1, l = к; 0, l Ф к} . При L = 0 все интегралы вычисляются по интервалу [0, 2], а граничных вейвлетов не возникает: [ A^ /A^ ] = -[R0 | R2 ]-1 R1 при условии, что а1 = {1, l = к; 0, l Ф к}. Ф0(t) Ф0 (2t- -k) Фх(0 = ЪНк Ф1С2t~ -к) к=0 .Фг(0. _ФД21- -к)_ 1 j, матрица Uразмерности (r +1) х(r +1) задана элементами где H2 = U~x Ли, Л = diag (2 -г-1 г*-2г- 0 = SH2S U k, j Аналогично запишем базисные вейвлет-функции (2r + 1)-й степени на уровне масштабирования L fL '1,0,fL 11Д,...,- ,L 1 г'" 'L . Тогда для уровня L - 1 можно выразить в виде матрицы-строки: у = функции 1 в виде линейных комбинаций функций ф^ т.е. 1 = фLQL, где блоки матрицы QL составлены из коэффициентов разложения (3). Соответствующие коэффициенты сплайна будем собит Гь,0 Гь,1 С0 , С0 , L,0 L,r рать в вектор С = С ..Ц) , Cj а соответствующие вейвлет-коэффициенты - (здесь символ Т обозначает операцию транспонирова т DLfi вы D1 ,D1 , L,r L,r в вектор D = .А D ния). С использованием обозначений для блочных матриц процесс получения CL из CL 1 и DL 1 может быть записан как [14] CL = [ PL Ниже представлены примеры матриц [P | Q О a center A0 Hi i L-1 C Ql ] (4) dl-1 соответствующих L = 1, 2, 3: H1 H T о о о о о I о H0T о A1left о H1 о A2eft Aright H2T Hi о Aright о H1 о I о о о о о о a inner A0 о о I о о inner inner A0 о о I о о inner Agght о о Aright о о I H1 HI о [ P11Q1] = IQ2] = ,[ P Hi A^ [ P3IQ3] = H о о о о I H i H i о о о Aleft о H о о о Aleft о Hi H 0l о о о о о H1 о о о о о H i H0T о о о о о H1 о о о о о H2T H0T о о о о о H1 о Здесь O обозначает матрицу (r + 1)-го порядка с нулевыми коэффициентами, тогда как I - единичная матрица (r + 1)-го порядка. Обратный процесс разбиения коэффициентов CL на более грубую версию CL-1 и уточняющие коэффициенты DL 1 состоит в решении системы линейных уравнений (4). Разрешимость данной системы гарантирована линейной независимостью базисных функций. Процедуру разбиения CL на часть CL1, соответствующую низшему масштабу, и уточняющие коэффициенты DL-1 можно применить рекурсивно и к самой этой части CL-1. Следовательно, исходные значения можно представить в виде иерархии грубых версий c масштабами C0, C1, ..., CL-1 и уточняющих деталей D0, D1, ..., DL-. Результирующее вейвлет-разложение эрмитового сплайна (2r + 1)-й степени SL(x) может быть записано в виде: L-1 2J SL (x) = Z Z Cf,k N0,k (x) + Z Z DJ,k Mjk (x) k=01 i=0 j=0 i=0 При этом по величине вейвлет-коэффициентов D1, j = 0, 1,., L - 1, можно судить о значимости соответствующих уточняющих деталей. Незначимые убираются с целью сжатия информации. 3. Проблема устойчивости мультивейвлет-преобразования При больших L матрицу [PL | Q] следует сделать блочной трехдиагональной, изменив порядок неизвестных так, чтобы блоки матриц/^ и (У перемежались (см. [15]): Нх I 0 ••• I , a < x < b. Hi 4left н1 о 4eft ■•. ■•. I L > 1. [PlQl ] 4right Hi H 0 0 I 1 у Тогда можно применить для решения системы (4) алгоритм блочной матричной прогонки [16]. Введем для L > 1 обозначения L = blocktridiag{H2T, 0, 0; A2left, 0, 0; H2', 0, 0; A2inner, 0, 0; ... ; I, 0, 0}; U = blocktridiag{0, 0, I; 0, 0, HT; 0, 0, Acinner; ... ; 0, 0, Acnght; 0, 0, Ш}; T = blockdiag{T, i = 0, 1,. , 2L}. В результате получим [PL | Ql] = (L + T)T-1(U + T), (5) где для блоков Ti справедливы выражения [17] •/, =Я,:72 =4вй-Нт2Т~1-,Т3 =Я,-A^T-'Hl; Tt =/-Я27:1Чппег;?;+1 =Н, -АГеТХ 0 = 4,6,...,2l -2); = -я^^а;^; T2l+l =НХ-TjHl Аналогично для L = 1 имеют место ( = Щ; ( = I - нЦ*61'; T = H - A^T-Hl. Если ввести обозначения T^CU + T)^"1, Df"1,..., T = z, Cf = Cf'1,..Ct Dl = /)/"',..., /.)/"''J, то процесс решения системы (4) разбивается на два этапа: (6) •L,r (L + T) z = CL; (I + T-1U) T r-1i ^L-1 nL-1 nL-1 = Z. Покомпонентно каждый из них можно записать в следующем виде: для решения (L + T) z = CL последовательно выполняются действия z = ( (7) T = z можно воспользоваться рекуррентными а для вычисления (I + T"1U)[c^"1, ..., , C^l! формулами: 1 = TT1 (CL) ;Z2 = T2-1((cL) -HTZ1);Z3 = T"1^) -Af^); =^"1((cf_1)T=ЩС1)Т - АГЧ) (/ = 4,6,...,2L - 2); T2L ((C2Ll-1 ) - H2 Z2L-1); Z2L +1 = T2L+1 ((ClL ) - Z2L ), Z2L = 1 2l ) = Z2L +1; (D2L~1) = Z2L T2lHV {pi- ) ; (cf-1 )T = zM -т;:114ппег (Dj-1 )T (i = 2L~X -1,2L~X -2,...,2); (nr)T =Z2 -г2-чт (c^-1)1 ;(coi_1)T = -7Г1 (Al_1)t • Здесь Zi суть подвекторы порядка r + 1. Нетрудно убедиться, что при L > 1, кроме первой и последней строк, условия даже слабого диагонального преобладания [16] не выполняются. Объективное представление об устойчивости или неустойчивости алгоритма вычисления мультивейвлет-преобразования дает наблюдение за поведением чисел обусловленности в эвклидовой норме прогоночных матриц Ti (таблица). Значения стандартной функции MathCad'a cond2(7i), i = 1, 2, ..., 33, r = 1, ..., 5 1 2 3 4 5 6 i cond2(Ti) r = 1 r = 2 r = 3 r = 4 r = 5 1 2 4 8 16 32 2 64,759 2,695e+004 2,799e+007 2,142e+011 8,333e+014 3 2,513 157,2 5,832e+004 8,799e+007 3,717e+010 4 2,953 60,32 7 805 2,468e+006 8,122e+008 5 3,493 478,7 1,817e+005 7,345e+007 7,162e+010 6 2,835 52,05 7 124 1,999e+006 7,353e+008 7 3,56 492,9 1,864e+005 7,264e+007 8,126e+010 8 2,831 51,79 7 102 1,991e+006 7,328e+008 9 3,563 493,4 1,865e+005 7,263e+007 8,158e+010 10 2,831 51,79 7 102 1,991e+006 7,327e+008 11 3,563 493,4 1,865e+005 7,263e+007 8,159e+010 32 533,592 3,61e+005 1,379e+009 5,382e+011 3,713e+017 33 84,003 2,933e+004 1,965e+007 7,101e+009 1,308e+015 Таким образом, можно полагать, что для любого уровня разрешения L они являются невырожденными, что на практике означает корректность представленного выше монотонного варианта алгоритма матричной прогонки (6)-(8). Для r = 1 (случай приближения кубическими мультивейвлетами) это согласуется с тем фактом, что на бесконечной числовой оси L2-устойчивость базиса гарантирована существованием оценок констант Ритца 0,414076 < A < B < 0,414265 [10], таких что 2 r ZZZ DJ,k MJ,k (x) k=0 i DJ1 AZ ZZ k=0 J i < Dj 1 внутри отрезка аппроксимации получаем (все вычисления проводились в системе MathCad, дробные числа там, где они затрудняют восприятие результата, приводятся округленно в виде конечной десятичной дроби) 1,1168 0,082 0,0159 0,0004" -7,894 -0,2968 - 0,0915 - 0,001 -67,9076 -6,9996 -1,1658 -0,039 ' 1792,4185 129,5384 26,3043 0,6684 На левом конце отрезка аппроксимации имеем 15,9089 1,2448 0,0667 0,0019 -12,3588 -0,7126 -0,0257 -0,0004 -1984,0775 -158,6934 -8,679 - 0,2531 ' 30335,1818 2401,4609 130,1543 3,7654 На правом конце отрезка аппроксимации значения совпадают с точностью до знака и перестановки. При этом H1 = diag (8,4,2,1)/8, f 384 840 0 -5040> f384 - 840 0 5040^ 1 -132 -228 360 2520 1 132 -228 -360 2520 > H 2 =- 768 18 24 -84 -360 2 768 18 -24 -84 360 v -1 -1 6 18 V v 1 -1 -6 18 V Для xe[0, 1], полагая на верхнем уровне разрешения L = 5, находим длину шага сетки h = 2-5=1/32. В соотношениях уточнения (4) при выполнении вейвлет-преобразования в качестве исходных используются значения функции и трех производных, всего 132 числа. Рассматривая в качестве тестовой функции многочлен седьмой степени f(x) = (2 - x)3(x2 - 1)2, находим на последнем этапе рекуррентного алгоритма вейвлет-разложения восемь значений сплайна S0(x) и трех его производных в концах отрезка C = [8, -12, -20, 138, 0, 0, 8, -48]T. При этом все вейвлет-коэффициенты пренебрежимо малы и могут быть обнулены, обеспечивая в данном случае коэффициент сжатия K = 132/8 = 16,5. Заключение В работе представлена общая схема построения и вычисления мультивейвлетов на основе эрмитовых сплайнов и свойства ортогональности многочленам. Исследована устойчивость и даны результаты численного эксперимента по восстановлению функции, заданной многочленом. Дополнительные возможности для оптимизации методов обработки численной информации состоят в построении сглаживающей сплайн-аппроксимации, альтернативной методу наименьших квадратов, и в сжатии дискретных данных на участках гладкости функций.
Добеши И. Десять лекций по вейвлетам : пер. с англ. Ижевск : Регулярная и хаотическая динамика, 2001. 332 с.
Чуи Ч. Введение в вейвлеты : пер. с англ. М. : Мир, 2001. 412 с.
Новиков И.Я., Протасов В.Ю., Скопина М.А. Теория всплесков. М. : Физматлит, 2006. 616 с.
Strela V. Multiwavelets: regularity, orthogonality and symmetry via two-scale similarity transform // Stud. Appl. Math. 1997. V. 98, No. 4. P. 335-354.
Strela V., Heller P.N., Strang G., Topivala P., Heil C. The application of multiwavelet filterbanks to image processing // IEEE Trans. Signal Processing. 1999. V. 8, No. 4. P. 548-563.
Warming R., Beam R. Discrete multiresolution analysis using Hermite interpolation: Biorthogonal multiwavelets // SIAM J. Sci. Comp. 2000. V. 22, No. 1. P. 269-317.
Dahmen W., Han B., Jia R.-Q., Kunoth A. Biorthogonal multiwavelets on the interval: cubic Hermite splines // Constr. Approx. 2000. V. 16. P. 221-259.
Han B. Approximation properties and construction of Hermite interpolants and biorthogonal multiwavelets // J. Approxim. Theory. 2001. V. 110. P. 18-53.
Koro K., Abe K. Non-orthogonal spline wavelets for boundary element analysis // Engineering Analysis with Boundary Elements. 2001. V. 25. P. 149-164.
Han B., Kwon S.-G., Park S.S. Riesz multiwavelet bases // Appl. Comput. Harmon. Anal. 2006. V. 20. P. 161-183.
Завьялов Ю.С., Квасов Б.И., Мирошниченко В. Л. Методы сплайн-функций. М. : Наука, 1980. 352 с.
Шумилов Б.М. Мультивейвлеты эрмитовых сплайнов третьей степени, ортогональные кубическим многочленам // Математическое моделирование. 2013. № 4. С. 17-28.
Strang G., Strela V. Short wavelets and matrix dilation equations // IEEE Trans. Signal Processing. 1995. V. 43, No. 1. P. 108-115.
Столниц Э., ДеРоуз Т., Салезин Д. Вейвлеты в компьютерной графике : пер. с англ. Ижевск : Регулярная и хаотическая динамика, 2002. 272 с.
Шумилов Б.М. Алгоритм матричной прогонки вычисления мультивейвлетов нечетной степени, ортогональных многочленам // Автометрия. 2015. Т. 51, № 2. С. 83-92.
Самарский А.А., Николаев Е.С. Методы решения сеточных уравнений. М. : Наука, 1978. 591 с.
Васильева Е.А. Достаточные условия корректности метода матричной прогонки // Известия КГТУ. 2011. № 20. С. 103-108.
Schneider A. Biorthogonal cubic Hermite spline multiwavelets on the interval with complementary boundary conditions // Results in Mathematics. 2009. V. 53, is. 3. P. 407-416.
Cerna D., Finek V. Construction of optimally conditioned cubic spline wavelets on the interval // Adv. Comput. Math. 2011. V. 34. P. 519-552.