Оптимальные коэффициенты при разложении сигнала по B-сплайнам Текст научной статьи по специальности «Математика»

Сер. 10. 2012. Вып. 2

ВЕСТНИК САНКТ-ПЕТЕРБУРГСКОГО УНИВЕРСИТЕТА

УДК 519.651

Н. А. Лебединская, Д. М. Лебединский

ОПТИМАЛЬНЫЕ КОЭФФИЦИЕНТЫ

ПРИ РАЗЛОЖЕНИИ СИГНАЛА ПО В-СПЛАЙНАМ*)

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

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

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

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

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

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

Лебединская Наталия Александровна — кандидат физико-математических наук, доцент кафедры параллельных алгоритмов математико-механического факультета Санкт-Петербургского государственного университета. Количество опубликованных работ: 10. Научное направление: локальная аппроксимация. E-mail: [email protected].

Лебединский Дмитрий Михайлович — кандидат физико-математических наук, доцент кафедры параллельных алгоритмов математико-механического факультета Санкт-Петербургского государственного университета. Количество опубликованных работ: 6. Научное направление: сплайн-вейвлеты. E-mail: [email protected].

*) Работа частично выполнена при финансовой поддержке Российского фонда фундаментальных исследований (гранты № 10-01-00-245 и 10-01-00-297).

(¡5 Н. А. Лебединская, Д. М. Лебединский, 2012

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

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

Второй подход [1] - использовать двойственную систему функционалов, оценивая отсутствующие значения производных рассматриваемой функции как значения производных интерполяционных многочленов. Такой подход дает лучшие результаты, чем предыдущий.

Третий подход [2] - применить ^-преобразование и сверточные соотношения для дискретных В-сплайнов.

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

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

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

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

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

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

..., а-п = 0,..., а-т-1 = 0, а-т > 0,...,ат > 0, ат+1 = 0,...,а„ = 0,...

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

Элементы базиса, по которому раскладывается вектор значений исходной функции, будем обозначать вг, так что вг(]) = Наша цель состоит в том, чтобы подобрать

такие коэффициенты сг, чтобы у отличалось как можно меньше от ^ сз вз.

зеъ

Подход расширенной свертки предполагает искать коэффициенты е в виде е =

кг

Iв^уг+з. Поскольку базисные функции отличаются друг от друга сдвигом, числа

з = -кг

кг и коэффициенты будем искать не зависящими от г.

Для дальнейшего фиксируем к > т и будем искать коэффициенты вз, —к ^ ] ^ к. Подставляя имеющуюся формулу для е в разложение для у, получаем

к

У1 - вг(1) X взУг+з ■

гех з = -к

Вспоминая определение вг и меняя порядок суммирования, находим

к т

У1 - Ш И а-гвзУг+з ■

3 = — к г=-т

Обозначая г + ] = в и заменяя переменные суммирования соответствующим образом, имеем

1+т+к / ш'т(8+т,к)

У1 ~ ^ ( У3 X аз-*вз

8=1-т-к \ з=шах(8-т,-к)

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

т+ к

шт тах

в 8 = -т-к

ш{п(з+т,к)

¿80 — аз-8 вз

з=шах(8-т,-к)

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

к

^ аз-звз = 8 — ( — 1)4 —(к +1)< в < к + 1. (1)

з=-к

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

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

Теорема. Система (1) однозначно разрешима.

Доказательство. Прибавим к каждому уравнению (кроме первого) системы (1) предыдущее, начиная с последнего. Это приведет к тому, что из всех уравнений, кроме первого, исчезнет неизвестная а; в первое уравнение основной системы она входит

с ненулевым коэффициентом, так что утверждение теоремы о системе (1), очевидно, эквивалентно аналогичному утверждению о системе

Теперь воспользуемся тем, что полученная система уравнений симметрична - если поменять местами в] и в—], будем иметь ту же самую систему.

Рассмотрим вспомогательную систему полученную из (2) подстановкой в—] = в] для ] > 0 и удалением уравнений, соответствующих —к ^ в ^ 0.

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

Также любое решение 5 без труда превращается в решение (2), таким образом достаточно доказать, что однородная система, соответствующая (2), не имеет нетривиальных решений. Такое утверждение, в свою очередь, следует из невырожденности матрицы, полученной из матрицы системы (2) удалением последней строки.

Матрица, о которой идет речь, является теплицевой ленточной. Поэтому ее невырожденность вытекает из того, что некоторый определитель [4], в который входят степени корней многочлена Эйлера и —1, отличен от нуля. Однако хорошо известно [5], что корни многочленов Эйлера простые, вещественные, отрицательные. Потому данный определитель раскладывается в произведение определителя Вандермонда от корней многочлена Эйлера и —1, что отлично от нуля, и многочлена Шура от тех же корней и —1, что также отлично от нуля, поскольку многочлен Шура однороден и все его ненулевые коэффициенты положительны, а корни многочлена Эйлера вещественны и отрицательны, так что в полученной сумме все слагаемые имеют один и тот же знак.

3. Решение задачи на минимакс. После доказательства того, что система (1) однозначно разрешима, логично было бы доказать, что ее решение дает решение задачи на минимакс. Но пока это можно доказать лишь для небольших значений т и произвольных значений к, так как на данном пути встречаются определенные трудности.

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

Второй шаг - завершение доказательства гипотезы. Используя критерий для ми-нимакса из работы [6], можно утверждать справедливость этой гипотезы, если найти набор неотрицательных коэффициентов, которые не все равны нулю, и такие, что линейная комбинация строк матрицы системы (1), умноженных на ( — 1)номер строки, с этими коэффициентами равна нулевой строке. Такое утверждение будем называть второй частью доказательства основной гипотезы.

Рассмотрим сначала первый шаг. Начнем с леммы, для формулировки и доказательства которой нам понадобятся некоторые факты и обозначения из [4], а также несколько собственных обозначений, которые здесь будут приведены.

В дальнейшем понадобятся обозначения и некоторые факты из [7].

(2)

Э = — к

Пусть —Л1 < —Л2 < ... < —\2т < 0 - корни многочлена Эйлера степени 2т. Определим разбиения Rkmpq = 1рц(2к + 1)т+1, а также многочлены Рт = ,т,т—1—1,—

, т, т—1—1, к+2. Рассмотрим также числа Лт = РтЛ,..., Лт, 1,Лт+ь ... ,Л2т).

Лемма 1. Утверждение первой части основной гипотезы следует из следующих соотношений:

1) Лкт1 > 0;

2) Лк,т,г+1 > Лкт1;

3) Лк ,т,1+1 Лкт1 > Лкт1 Лк,т,1 —1.

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

q-l

-г - -^аг-г(а„(% + -) + ап(1\д+ - + 1)), г = 0, .. ., то - 1

1=0

Для дальнейшего введем следующее обозначение:

В

(—Л1)0

( — Лт)0

( —1)0 0 ( —Лт+1)

(—Л2т )0

(—Л1)т-1 (—Л1)2к+1+т

( — Лт )

( — 1)

( \ \т— 1 ( —Лт+1)

т— 1 т — 1

( —Л2т)

т— 1

( —Лт) (—1) ( — Лт+1 )2к+1+т

2к+1+т 2к+1+т

( —Л2т)2к+1+т

(-Л1)2к+1+2т

(—Лт) ( — 1)

(—Лт+1)2к+1+2т

2к+1+2т 2к+1+2т

(-Л2т)2к+1+2т

Затем также понадобятся определители Вц, полученные из В заменой степени I на т + к + г — 1, где I = 0,...,т — 1 и г = 1, 2.

Вычисляя интересующие нас значения функций, после некоторых преобразований определителей (уже в приведенном выше обозначении) находим

/о =

Е(—1)

I Вц + В2

1=0

В

где /о - те значения, про которые нам нужно доказать, что они ^ а. Осуществим это так: покажем, что в предположениях леммы числа /о возрастают, поскольку /т—1 = а. Данный факт не изменится, если домножить числа /о на одно и то же положительное число. Пользуясь таким утверждением, достаточно доказать, что числа

9о

Е(—1)1Л

кт1

1=0

также возрастают, а это уже очевидным образом вытекает из условия леммы.

Лемма 2. Соотношения из условия леммы 1 следуют из того, что корни многочлена Эйлера удовлетворяют таким неравенствам:

2т

1) £ Л < 1;

г=т+1 2т

2) £

г=т+2

Доказательство. Воспользуемся представлением многочлена Шура как суммы одночленов, соответствующих полустандартным таблицам Юнга.

Начнем с первого утверждения о том, что Лит > 0. Это число получено как разность двух многочленов Шура, причем диаграмма Юнга, соответствующая второму, получена из диаграммы, отвечающей первому, пририсовыванием одной клетки. Поэтому можно сгруппировать одночлены, входящие в указанную разность, следующим образом: отдельные одночлены, у которых в соответствующей таблице Юнга в клетке сразу над местом, куда пририсовывалась дополнительная клетка, стоит максимально возможное число - это те одночлены из первого многочлена Шура, у которых нет соответствующих во втором многочлене; далее, идут разности - одночлены из первого многочлена Шура, у которых в соответствующей таблице Юнга в клетке сразу над местом, куда пририсовывалась дополнительная клетка, стоит число, меньшее максимально возможного, и из него вычитаются те одночлены из второго многочлена Шура, которые соответствуют таблицам Юнга, полученным пририсовыванием одной клетки и заполнением ее некоторым числом, из таблицы Юнга уменьшаемого одночлена. Будем утверждать, что все такие группы одночленов (одиночные одночлены - очевидно, и разности тоже) неотрицательны.

Для разностей это следует из следующих утверждений: число клеток в столбце, куда надо пририсовывать новую, равно то+1. Учитывая строгое возрастание чисел в столбце полустандартной таблицы Юнга, в новой клетке может появиться число не меньше то + 2 - следовательно, сумма вычитаемых одночленов группы получается из уменьшаемого при помощи умножения на сумму модулей корней многочлена Эйлера, которые меньше 1, а она (см. условие леммы) предполагается меньше 1.

Аналогичные рассуждения позволяют доказать и второе, и третье утверждения леммы 2. Для доказательства второго утверждения мы рассматриваем разность Ли,т,1+1 — Лит и объединяем группы одночленов из обоих сумм Л в большие группы способом, аналогичным применявшемуся ранее. Как и в предыдущем случае, утверждаем, что полученные укрупненные группы разностей неотрицательны.

Третье утверждение - самое сложное. Здесь мы опять пририсовываем клетку и по этому принципу группируем разности разностей одночленов. Если в клетке над тем местом, куда пририсовывалась вторая из дополнительных клеток, стояло число п, то соответствующая разность в Ли,т,1+1 — Лит имела вид х(1 — Е2=П Л»), а в Лит1 — Ли,т,1-1 - вид X Е2=п Л*(1 — Е Л-), где X - разность одночленов, полученная при пририсовывании первой клетки. Вычитая одно из другого, имеем

2т 2т 2т

1 — 2£Лг + £ Е лл

г=п г=п -=г+1

Положительность этого выражения следует из второго неравенства в условии леммы, если принять во внимание тот факт, что п ^ то + 2, поскольку длина столбца (первого в соответствующей диаграмме Юнга), к которому мы собираемся пририсовывать клетку, не менее то+2.

Условия леммы 2 для небольших то могут быть проверены при помощи (довольно громоздких, но точных) вычислений следующим образом.

По многочлену Эйлера можно построить другой многочлен f, корнями которого будут всевозможные суммы корней исходного многочлена Эйлера в количестве m штук. Это утверждение есть тривиальное следствие хорошо известной теоремы о том, что всякий симметрический многочлен является многочленом от элементарных симметрических многочленов. В системе компьютерной алгебры MAXIMA в пакете SYM имеется функция direct для такого преобразования многочленов. Далее, по полученному многочлену f можно построить еще один многочлен g, корнями которого будут величины, обратные корням f. Теперь первое условие леммы 2, очевидно, будет верно тогда и только тогда, когда минимальный корень g (а это будет величина, обратная сумме тех корней многочлена Эйлера, которые больше —1) сам будет < —1. Поскольку такой корень в принципе может быть только один, данное условие равносильно тому, что знак g( —1) отличается от знака g (—ж).

Второе условие проверяется чуть сложнее, чем первое. Как и для первого, построим многочлен g, корнями которого будут величины, обратные суммам корней многочлена Эйлера, но уже в количестве не m, а m — 1, и наше условие будет состоять в том, что у этого многочлена есть корни < —2. Поскольку таких корней может быть несколько, для проверки данного условия недостаточно знать знак только g(—2). Однако, поскольку все корни многочлена Эйлера вещественны, все корни g тоже вещественны, так что для проверки того, что у g нет корней, меньших -2, достаточно знать знак всех производных g в точке -2.

Эти вычисления были успешно проведены для m = 1, 2, 3, 4 с использованием системы компьютерной алгебры MAXIMA версии 5.23.2 для Linux (Fedora 14). Здесь нужно только указать, что для m = 4 размера стека управления по умолчанию в 2 Мб недостаточно, и нужно установить его в 20 Мб, задав соответствующую опцию. Для m = 4 один из многочленов имеет степень 70 и вычисляется за несколько минут.

Заметим, что приближенные вычисления можно провести гораздо быстрее и для значительно больших значений m, хотя такие значения m в реальных задачах маловероятны.

Наконец, приведем вторую часть доказательства основной гипотезы. Она следует из того, что искомые коэффициенты (по первым 2k + 2 строкам) выражаются как определители, полученные из соответствующей матрицы вычеркиванием соответствующей строки, умноженные на знак, чередующийся при переходе от одной строки к следующей. Но у наших строк знаки тоже чередуются, поэтому итоговый знак будет такой же, как у определителей подматриц исходной (с положительными элементами) матрицы (или везде другой - тогда изменим его у всех коэффициентов одновременно). А у исходной матрицы, являющейся подматрицей матрицы значений B-сплайнов в узлах сетки, все миноры неотрицательны [8].

Литература

1. Демьянович Ю. К. Гладкость пространств сплайнов и всплесковые разложения // Докл. РАН. 2005. Т. 401, вып. 4. С. 1-4.

2. Unser M., Aldroubi Akram, Eden Murray. Fast B-spline transforms for continuous image representation and interpolation // IEEE Transactions on Pattern Analysis and Machine Intelligence. 1991. Vol. 13, N 3. P. 277-285.

3. Renhong Wang, Yan Xu, Zhiqiang Xu. A spline interpretation of Eulerian numbers / / arXiv:0808.2349v2. 2008.

4. Trench W. F. Explicit inversion formulas for Toeplitz band matrices // SIAM J. on Algebraic and Discrete Methods. 1985. Vol. 6, N 4. P. 546-554.

5. Frobenius G. Über Bernoullische Zahlen und Eulersche Polynomen // Sitzungsberichte der Preussischen Akademic der Wissenschaften. 1910. P. 809-847.

6. Иванов В. К. Задача о минимаксе системы линейных функций // Матем. сб. 1951. Т. 28(70), вып. 3. С. 685-706.

7. Стенли Р. Перечислительная комбинаторика. Деревья, производящие функции и симметрические функции / пер. с англ. Д. М. Лебединского, Н. В. Цилевич, М. А. Всемирнова; под ред. А. М. Вер-шика. М.: Мир, 2005. 767 с.

8. De Boor C., deVore R. A Geometric proof of total positivity for spline interpolation // Mathematics of Computation. 1985. Vol. 45, N 172. P. 497-504.

Статья рекомендована к печати член-кор. РАН, проф. Г. А. Леоновым. Статья принята к печати 28 февраля 2012 г.