CC BY
27
УДК 004
Вестник СПбГУ. Сер. 10, 2008, вып. 1
Н. А. Лебединская, Д. М. Лебединский
МНОГОПОТОЧНЫЙ АЛГОРИТМ СПЛАЙН-ВЕЙВЛЕТНОГО СЖАТИЯ ЦИФРОВОГО ПРЕДСТАВЛЕНИЯ СИГНАЛА *>
В последнее время персональные компьютеры часто оснащают процессорами, поддерживающими почти одновременное выполнение двух потоков (одновременное, если потоки используют разные узлы процессора), а также двухъядерными. Возможны и комбинации, сочетающие в себе и то, и другое. Для эффективного использования возможностей таких процессоров нужны многопоточные алгоритмы, число потоков в которых соответствует возможностям процессора.
Данная статья посвящена разработке многопоточных алгоритмов для адаптивного сплайн-вейвлетного сжатия цифрового представления сигнала.
1. Формулы декомпозиции и реконструкции. Рассматриваемые в работе алгоритмы основаны на формулах декомпозиции и реконструкции сплайн-вейвлетного преобразования на неравномерной сетке, которое было предложено Ю. К. Демьяновичем в работе [1]. В отличие от традиционных вейвлетов локальный базис, представленный в ней, позволяет использовать неравномерную сетку. Это означает, что те точки на вещественной прямой, в которых записываются значения интересующей нас функции (сигнала), в дальнейшем называемые узлами (в традиции литературы по численным методам), могут не образовывать арифметическую прогрессию.
Весьма существенной платой за возможность использования неравномерной сетки является потеря ортогональности функций базиса. Однако в работах Ю. К. Демьяновича для предложенного им локального базиса строится двойственный базис, т. е. базис в пространстве функционалов, применение элемента которого к функции / дает соответствующую координату / в локальном базисе. Это позволяет получить формулы декомпозиции и реконструкции, но затрудняет оценку погрешности при удалении узлов в стандартной норме (той самой, которая давалась скалярным произведением, относительно которого были ортогональны стандартные вейвлеты).
В настоящей работе используются два следствия из возможности пользоваться неравномерной сеткой: во-первых, при сжатии узлы можно удалять по одному, причем при удалении конкретного узла необходимо подправлять значения только нескольких (их число зависит от порядка участвующих в локальном базисе сплайнов и функционалов в двойственном базисе и никак не зависит ни от сетки, т. е. набора узлов, ни от сигнала) соседей справа и слева.
В данной работе применяются В-сплайны второго порядка. Формулы декомпозиции и реконструкции для таких сплайнов получены Ю. К. Демьяновичем и А. В. Зиминым в работе [2], однако для удобства они приводятся ниже. Для указанных сплайнов при удалении узла необходимо подправлять коэффициенты, соответствующие двум предыдущим узлам, а также вычислять функцию, добавляемую к сигналу в качестве погрешности при удалении узла. Эта функция представляет собой одну из базисных функций (до удаления узла), коэффициент при которой дается одной из формул декомпозиции. Носитель данной функции лежит от левого соседа удаляемого узла до второго соседа удаляемого узла справа.
Работа выполнена при финансовой поддержке Российского фонда фундаментальных исследований (гранты № 04-01-00692 и 04-01-00026).
© Н. А. Лебединская, Д. М. Лебединский, 2008
Перейдем теперь к краткому изложению формул декомпозиции и реконструкции, вводя по мере необходимости соответствующие понятия.
.В-сплайном второго порядка называется кусочно-квадратичная функция одной переменной Ь с 4 параметрами х\ < х2 < хз < Х4, определяемая формулой
'(„51,). ье[хих2]
«-Х1)2 , (Х1-Х4)(^-Х2)2 . г- г 1
о С/4) — ^ (х1 —хг)(х1 —хз) (Х1-Х2)(Х2-ХЗ)(Х4-Х2) 1 с 1^2,
"Х1,Х2,ХЗ,ХДЧ - (¡-х4)2 . сГ 1 •
(Х4-Х3)(Х4-Х2) ' Г (= [Хз, ДГ4]
О, £ ^ [11,14]
Указанная функция является финитной (ограниченный носитель) и гладкой.
Пусть, далее, имеется отрезок [а, 6] (на котором рассматривается интересующий нас сигнал), и введена сетка < Хо < = а < х2 < ■ ■ ■ < хп = Ь < жп+] < жп+2- В качестве базисных функций рассматриваютсяограничения на отрезок [а, 6] В-сплайнов
(0 = ^ = 1,0,..., 71 1.
В качестве биортогоналыюй системы функционалов могут выступать два набора функционалов, определяемых формулами ^(/) = /(24) + (эц+1 — ж»)/'(£*)/2, ЩО./) = ¿1,3+1 и Fi(f) ~ /(хг) - (х{ - х^1)}'{хг)/2, = 5^+2- Поскольку не-
возможно брать значение сигнала и его производных вне отрезка [а, 6], приходится сочетать функционалы из обеих систем; это можно делать многими разными способами, здесь используется такой:
п-1
¿=о
где / - приближение сигнала в виде разложения по базису из 5-сплайнов, которое подвергается сжатию. Коэффициенты в этом разложении (в порядке следования) в дальнейшем обозначаются о_1,ао,... ,ап-1- Также под двойственным базисом к системе функций в дальнейшем будем понимать систему функционалов, использованную в указанном разложении /.
Пусть теперь нужно удалить узел х1у / ^ 4. Обозначим новые базисные функции (соответствующие сетке без узла х{) через Г^(£); имеем
Uj(t) = Q.j{t) при 3 = -1,0,... ,г - 4, г + 1,... ,п - 1;
^г — 1 = ,х< + 1 ,х,+2,х1+з 1
отсутствует.
Новое представление того же сигнала / запишем в виде
п-1
/ = 6-10-1(4) + XI +
3=
где последнее слагаемое - это погрешность, возникающая при удалении узла ж;. Действительно, все функционалы двойственного к flj(t) базиса, соответствующие узлам новой сетки, равны 0 на функции П;_1(4).
Итак, формулы декомпозиции выглядят следующим образом: bj = а] при з = — 1,0,..., I — 3, г + 1,..., п — 1,
£¿+1 Х{ Х{—2 Яа+1 Ог-2 = ---Йг-3 Н--~-Щ-2,
Ьг-1 = СЧ,
ь _ (Х1+! - Хг)(Хг+2 - Хг) д (^¿+2 - Ж{)(з;г-2 ~ ж' + 1)д +
(Хг+2 - - Ж.) г~3 (Хг+2 ~ Х{-г)(х{-2 ~ Х{) 1~~2
Х{ — 1 Хг
+ аг_ 1--аг.
Хг-1 - Х1+2
Формулы реконструкции выглядят так:
Ь2 при з = — 1,0,..., г — 3, г + 1,... ,п — 1,
Хг+1
-Ог-_з Н--Ог_2,
Хг+1 - Хг-2 Яг-2 ~ £¿+1
-Ьг-2 + -»1 — 1 + О»,
Хг+2 - £¿-1 ~
Ьг-1-
Строго говоря, в работе [2] формулы декомпозиции и реконструкции выписаны для случая периодического сигнала, здесь же рассматривается непериодический случай. Однако, ввиду того, что первые три узла не удаляются, в нашем случае формулы будут те же самые, с точностью до обозначений (мы придерживаемся других обозначений).
2. Алгоритм сжатия. Второе обстоятельство, вытекающее из возможности использования неравномерной сетки и возможности удалять узлы по одному, состоит в возможности использования адаптивного алгоритма. Вообще говоря, в самой требовательной постановке задача сжатия сигнала может выглядеть так: при заданном допустимом уровне погрешности (например, чтобы заранее выбранная норма разности исходного и полученного сигналов не превосходила заранее заданного положительного числа) удалить максимальное возможное число узлов, или двойственная задача - при заданном числе удаляемых узлов (заданной степени сжатия) достичь наименьшей погрешности. Однако в такой постановке задача сжатия сигнала кажется авторам безнадежной.
Несмотря на вышеуказанные трудности, требования задачи-максимум можно ослабить таким образом, чтобы в итоге получить несколько вполне решаемых за приемлемое время задач. Поскольку базис является локальным, произведение функций, добавляемых к сигналу в качестве погрешности при удалении достаточно далеких друг от друга узлов, тождественно равно нулю, т. е. они ортогональны при любом разумном скалярном произведении. Достаточно далекими являются узлы, удаление которых приводит к пересчету коэффициентов, отвечающих непересекающимся множествам узлов. Для рассматриваемых в данной работе £?-сплайнов второго порядка между такими узлами должны лежать по крайней мере два узла.
В формулах декомпозиции при удалении х, используется значение коэффициента (¡¿—з, что на первый взгляд может приводить к взаимному влиянию при удалении Хг и
«г-2 =
аг-1 = а, =
з, однако, во-первых, в этих формулах не участвует а во-вторых, при удале-
нии узла его коэффициент просто переписывается в предыдущую позицию, т. е. нам повезло и вне зависимости от того, надо ли удалять Xj_3, коэффициент в третьем узле слева от Xi будет одним и тем же.
Таким образом, если ограничиться вопросом об удалении узлов с номерами, образующими арифметическую прогрессию с первым элементом 4 и разностью 3, погрешность сжатия (вычисляемая по ¿2-норме) будет вычисляться как корень из суммы квадратов погрешностей, вносимых удалением отдельных узлов. Это позволяет, задавшись заранее максимально допустимым уровнем погрешности, выбрать в качестве удаляемых узлов такие, удаление которых приводит к минимальным погрешностям, в количестве, определяемом тем обстоятельством, чтобы сумма их квадратов не превосходила квадрата общей допустимой погрешности. Такая процедура может быть затем использована повторно, возможно с другим допустимым уровнем погрешности. К сожалению, прогнозировать общую погрешность в результате нескольких применений такой процедуры трудно, однако она в любом случае не превосходит суммы допустимых уровней погрешности, установленных на каждом из прогонов.
У рассмотренного метода есть и недостатки. В качестве первого можно отметить плохую его приспособленность к удалению длинных последовательностей подряд идущих узлов. Второй недостаток состоит в том, что для применения этого метода следует иметь цифровое представление сигнала (разложение по исходному базису) все целиком. Однако это - плата за глобальную оптимизацию погрешности. Если это невозможно, можно организовать буферизацию и разделить представление сигнала на части, расположенные на разных отрезках, и затем применять метод на каждом отрезке отдельно.
Теперь встает вопрос о действиях, необходимых для однократного применения вышеописанной процедуры к результатам разложения сигнала по локальному базису. Первым действием является пересчет коэффициентов разложения, соответствующий удалению рассматриваемых узлов, с вычислением погрешностей, вносимых удалением каждого из них. Ввиду того обстоятельства, что «зоны влияния», т. е. те наборы узлов, коэффициенты которых надо пересчитывать в связи с удалением узлов 4, 7, ..., не пересекаются, такой процесс может быть распараллелен почти полностью. Новые коэффициенты попадают во второй массив.
Далее вычисляется набор погрешностей, возводим их в квадрат и должны выбрать набор минимальных из них, сумма которых не превосходит заданной величины lim. Первое, что приходит в голову, - отсортировать множество погрешностей и затем выбрать подходящее количество минимальных. Однако, воспользовавшись алгоритмом, лежащим в основе пирамидальной сортировки (см., например, [3]), можно несколько сэкономить.
Представляется разумным объединить в пары погрешности и номера соответствующих узлов и устроить очередь по приоритетам из пар, соответствующих узлам, назначенным для удаления; в качестве приоритета в этой очереди будет использоваться величина погрешности. Потребуются также две вспомогательные переменные, min из которых хранит минимум из просмотренных, но не попавших в очередь, погрешностей на данный момент (если таких нет, min = max, т. е. максимальной погрешности из попавших в очередь) и sum - сумму имеющихся в очереди погрешностей.
Перебирая последовательно погрешности, поддерживается следующее условие: очередь содержит пары, погрешности которых минимальны из всего множества просмотренных к данному моменту погрешностей, и любая просмотренная погрешность из числа не попавших в очередь больше разности предельного значения суммы lim и
суммы погрешностей в очереди sum. Изначально очередь пуста. Также, естественно, поддерживаются утверждения о правильности значений дополнительных переменных. Рассматривая очередную пару, разбираем следующие случаи (здесь х - погрешность из рассматриваемой пары):
1) х ^ min. Не делаем ничего;
2) min > х ^ max, х > lim — sum. Заносим значение х в min.
3) min > х, х ^ lim — sum. Заносим текущую пару в очередь, добавляем х к sum.
4) max > х, х > lim—sum. Удаляем пару с максимальной погрешностью из очереди, заносим в очередь текущую. Соответствующим образом подправляем min и sum.
Экономия по сравнению с полной сортировкой погрешностей здесь состоит в том, что максимальное количество пар в очереди по приоритетам ограничено числом пар, составляющих ответ, тогда как при полной сортировке оно достигает всего количества сортируемых элементов.
В итоге этого процесса в очереди оказываются пары, номера узлов в которых и дают удаляемые узлы.
Только что описанный процесс распараллелить сложнее, чем предыдущий (вычисление новых коэффициентов, соответствующих удалению узлов). Однако при большом числе узлов и небольшом числе потоков можно просто разделить набор погрешностей на две (для двух потоков), примерно равные по количеству части и обработать их независимо (с тем же значением lim, что и для всего исходного набора), и затем объединить результаты и еще раз применить к ним только что выписанный алгоритм.
Заключительный этап состоит в том, чтобы получить третий, окончательный, массив коэффициентов. В него попадают исходные коэффициенты из первого массива для узлов, которые остаются, и новые коэффициенты из второго массива для «зон влияния» удаленных узлов. Очевидно, данный процесс также может быть распараллелен почти полностью.
3. Алгоритм восстановления. Вообще говоря, если все, что нам нужно, это получение значений сигнала в определенных точках, то никакого явного восстановления, т. е. вычисления разложения по тому локальному базису, по которому сигнал был разложен до сжатия, не требуется. Тем более, что процедура восстановления сигнала его почти не меняет (с точностью до ошибок округления). Но если по каким-то причинам восстановление все же требуется, кажется разумным проводить его в порядке, обратном тому, в котором происходило сжатие, за такое же число «обратных прогонов», на каждом из которых восстанавливаются те узлы, которые были удалены на соответствующем прогоне сжатия. При этом тот факт, что при сжатии на каждом конкретном прогоне между удаляемыми узлами лежали по крайней мере два оставляемых узла, обеспечивает почти полное распараллеливание восстановления на конкретном прогоне так же, как он обеспечивал почти полное распараллеливание первого этапа сжатия.
Summary
Lebedinskaya N. A., Lebedinski D. М. A multithreaded algorithm of the spline-wavelet compressing of a digital representation of a signal.
A multithreaded algorithms of digital signal processing is established for the compression and decompression, using the decomposition and reconstruction formulas of the spline-wavelet transformation in the case of an irregular net.
Литература
1. Демьянович Ю. К. Всплесковые разложения в пространствах сплайнов на неравномерной сетке // Докл. РАН. 2002. Т. 382, № 3. С. 313-316.
2. Демьянович Ю. К., Зимин А. В. Всплесковое (вейвлетное) разложение пространств периодических 5-сплайнов второй степени на неравномерной сетке // Вестн. С.-Петерб. унта. Сер. 1: Математика, механика, астрономия. 2006. № 3. С. 72-83.
3. Седжвик Р. Фундаментальные алгоритмы на С++ / Пер. с англ.; Под ред. Ю. Н. Ар-теменко. СПб.: ООО «ДиаСофтЮП», 2002. 688 с.
Статья рекомендована к печати проф. Л. А. Петросяном.
Статья принята к печати 11 ноября 2007 г.
