Математические операции над двумерными сигналами в базисах сферических гармоник

Куликова Л. И. (Kulikova@impb.ru), Махортых С. А. Институт математических проблем биологии РАН

Сферические функции являются важным классом специальных функций, тесно связанных с классическими ортогональными полиномами. Они возникают при решении широкого круга проблем, например, при решении уравнения Лапласа в сферических координатах. Так как непрерывные решения уравнения Лапласа называют гармоническими функциями, то сферические функции называют также сферическими гармониками. Нахождение решения уравнения Лапласа Аи = 0 в сферических координатах r,6 ,ф , определение явного выражения для сферических функций, свойства последних подробно описано в работах [1, 2]. Сферические гармоники и их свойства достаточно широко представлены и в других публикациях зарубежных и отечественных авторов [3-6].

Сферические функции имеют общий вид: где функция ©1т определена следующим образом:

Ylm (6 ,Ф) = = етф© m (cose), (1)

© m (x)

2/ +1 (l - m)!pm

2 (l + m)!

здесь if (x) = (1 - x2)m/2 .

l dxm

Функции p/ (x)- функции или полиномы Лежандра, называются еще зональными гармониками. Функции p/4 (x) называются присоединенными функциями Лежандра. Так как е1ф = COSф + i sinф , выражение (1) перепишется в виде

----- p (cos6) cos /ф + ipm (cos6 )sin ,

4n (l + m)!

где функции рп (cos тф,sin тф ) = pm (cos6 )cos тф, pr/ (cos6 )sin тф называются тессеральными гармониками, по-видимому, по названию одного из видов игры в кости,

известного в древнем Риме [3]. При m = l тессеральные гармоники называются секторальными гармониками.

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

2 l

исходной функции. Поскольку функция (x - 1) имеет l - кратные нули в точках ± 1, ее первая производная будет иметь (l - 1) - кратные нули в этих точках и один нуль, лежащий между ними, вторая производная будет иметь два нуля между точками -1 и +1, а функция Pl будет иметь l действительных нулей, лежащих между точками -1 и +1. Степени sin 0 равны нулю в рассматриваемом интервале только при x = ±1, поэтому функция p (l > 0) имеет нули при 0 = 0 u П , и еще m - l действительных нулей.

Таким образом, зональные гармоники сохраняют знак в каждом из l +1 поясов, заключенных между двумя параллелями; области, содержащие полюса, также считаются поясами. При увеличении индекса m на 1 число параллелей уменьшается на 1, а число меридианов увеличивается на 2, где соответствующая гармоника обращается в нуль. Из вида функции

Pl (x)={ h й t( x-1) l (x+1)l f

применяя формулу Лейбница, легко установить, что все члены, кроме одного, обратятся в нуль при x = ±1 , в силу чего получим

pl (1) = 1, pl (-1) = (- Г/.

В действительности функция (x) < 1.

Имеют место и соотношения ортогональности. Каждая зональная гармоника pi ортогональна полиному, зависящему от аргумента x = cos0 более низкой степени.

Ортогональность определяется формулой

J Pk (x) Pn (x)dx = 5 , -1 2n +1

где 8!(п символ Кронекера.

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

J Pnm (х) PH (x)dx = m

-i (n - m)! 2n +1

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

т т s s

pl cos тф, pl sin тф, pl cos sф, pl sin sф ортогональна, интеграл от произведения их по ф обращается в нуль, за исключением т = s , когда берется интеграл от квадрата гармоники, беря в обоих множителях либо косинус, либо синус. Линейная независимость стандартных гармоник У1т (6 ,ф ) (1) легко доказывается [3].

Что касается вычисления присоединенных полиномов , можно воспользоваться явной формулой [3, 5]:

Pk(cos6) = (n + k)! (1 - cos2 6)-k/2 j(2)!c0s2-n+k6 =

п 2п(п - к)! =0!(п - )!(2 - п + к)!

= НГ(1-cos26)k/2 j(2)!cos2-n-k 6 , (2)

2 =0!(п - )!(2 - п - к)!

или

ss ал (П - S)l 2а л*12 Д (n + s)! n!(cos6 - 1)П s u n!(cos6 + 1)

psn(cos6) -V(1 - cos26)s/2 j---------.

2n(n!)2 U=0(s + u)!(n - u)! (n - s - u)! u!

Однако при больших значениях n приходится оперировать очень большими числами, что приводит к потере точности, переполнениям регистра и прочим неудобствам. Поэтому имеет смысл воспользоваться рекуррентными соотношениями для присоединенных полиномов Лежандра [7,8]:

(l - т)Рг = x(2l -1)Plт1 - (l + т -1)P-2, Рт = (-1)т(2т -- x2)т/2,

Ртт+1 = х(2т +1) Pz

Процедура вычисления присоединенных полиномов Лежандра имеется [9].

Система сферических гармоник является полной в пространстве интегрируемых с квадратом функций f (6 ,ф). Ранее Джеффрис Г. и Свирлс Б. [3] доказали, в частности, для гравитационного потенциала на сфере, когда притягивающее вещество находится

со n

f (0 ,ф) = II (anspsn C0S + bnspsn Sin X (3)

n=0s=0

где ans , bns - коэффициенты разложения, которые требуется определить.

Необходимо заметить, что это разложение обладает свойством, аналогичным свойству разложения Фурье или любого другого разложения по ортогональным функциям [11]. Именно, коэффициенты разложения подбираются так, чтобы I - сумма конечного

числа членов ряда гармоник удовлетворяет выражению: JJ[/(0 ,ф) - I]-(1\\йф - 0. Таким образом, как и в случае разложения по классическим ортогональным полиномам, решается задача минимизации функционала Jj[/(0 ,ф) - I]-й\\Лф. Имеется также аналог теоремы Парсеваля

JJ[/(0,ф)]2= I Iain(pn) C0s- S(P + I Ibn-\\{psn) sin- ФФ ,

n=0s=0 n=0s=0

означающий, что среднее значение квадрата f на сфере равно сумме средних значений

квадратов составляющих гармоник.

Принимая во внимание все изложенное, получим:

n n

B(0 ,ф) II (ankpk C0S kф + bnkpk Sin kф)

n=0k=0

С учетом свойства ортогональности коэффициенты разложения ank , bnk вычисляются по формулам

1 2n

ank JJ(pnk )-C0S- kф ФФ

-1 0

1 2п п 2п

JJb(0 ,ф)pkn C0s kфd\\dф =JJb(0 ,ф)pkn C0skф sin0d0dф

-1 0 0 0

1 2n

bnk JJ(pn )-sin- k(P ФФ

-1 0

1 2п п 2n (4)

JJb(0,ф)pkn sinkфdцdф =JJb(0 ,ф)pkn sinkф sin0d0dф

-1 0 0 0

либо целиком вне сферы, либо целиком внутри нее. Обобщение на случай, когда функция f (0 ,ф ) имеет непрерывные вторые производные доказали Курант и Гильберт [10].

Итак, если некоторая функция f (0 ,ф) интегрируема с квадратом на сфере, она может быть разложена в ряд по сферическим гармоникам

1 2п 1 2п п 2п

an0\\{pn) d\xdp = JJB(0,ф)d)idp =JJB(0,ф)sin0d0dp.

-10 -10 0 0

Коэффициенты Jn0 равны нулю в силу того, что в формулу вычисления входит sin /ар .

Отсюда видно, что коэффициенты разложения выражаются через определенные интегралы.

Значения нормировочных множителей вычисляются по формулам

J j (Pn )2 sin0d0dp =1Г+-Л ,

JJ [Pj (cos кф ,sin кф )sin 0 d\\dp=---)-.

Что касается интегралов, стоящих справа в выражениях для коэффициентов разложения (4), предлагается при вычислении их воспользоваться квадратурной формулой Гаусса [1, 8], которая дает возможность интеграл заменить суммой

j f(у)dy ~ J-a iwif ( yiX

a 2 i=1

J - a J + a

где у = -2- xi +--2-, Xi - i-тый нуль многочлена Лежандра pn (x), wi - весовые

коэффициенты. Процедуры вычисления нулей полиномов и весовых коэффициентов представлены в [9].

Погрешность аппроксимации функции определяется по формуле

i iai Jj(psn) cos2 sp d\xdp + i i J2m Jj(psn) sin2 sp d\xdp

n=0s=0 n=0s=0

5 = 1 -

п 2п

JJ B 2(0 ,p )sin 0d0dp

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

В силу того, что обобщенный спектрально-аналитический метод (ОСАМ) предполагает проведение полной обработки информационных массивов в пространстве коэффициентов [13, 14], возникает задача вывода аналитического соотношения в

пространстве коэффициентов при типовых операциях с сигналами, например, для

ф дВ(0 ,ф) дВ(0 ,ф)

получения производных от функции двух переменных-, -.

Запишем

дВ(0 ,ф) n-1 n д k д

II {Clnti- (P>* kф ) + bnki- (Pnk sin kф )).

д0 n=0k=0 д0 д0

Найдем производную по 0 от присоединенного полинома Лежандра (2):

dpm -(-1)n+m sinm0i (-1)k(2k)!C0s2k-n-m0 / 2n iTk!(n - k)!(2k - n - m)! 0 ~

(-1) m-1Q ov(-1) (2k)!C0s 0 ----msin 0 C0s0I------+

2n k k!(n - k)!(2k - n - m)!

+ (-1)n+m (- m+10)I (-1)k(2k)!C0s2k-n-m-10(2k - n - m) 2n ( * k!(n - k )!(2k - n - m)!

-n-m

C0s0 (-1)n+m m (-1)k(2k)!C0s2k-n-m 0

= m-----sin m 0I------

sin0 2n k k!(n - k)!(2k - n - m)!

n+m+1 k 2k-n-m-1 (-1) m+1 0 I (-1) (2k)!C0s ()

2n(-1) It k!(n - k)!(2k - n - m - 1)!

m m +1

= m Ctg 0- С + Pn

С другой стороны,

дВ(0, ф) n-1 n

д0 n=0k=0

II (cnkpkcos k(P + dnkpk sin kф ).

следовательно, можно записать

>m+1

/ + a Pnm+1(C0s0),

cnk = (m Ctg 0 + -n----) ank

nk v б (C0s 0 )> nk

d ( t 0 + Pm+1(C0s0)) b

nk V 5 p (C0s0 nk

здесь F = a(ar + r2 - (r0, r)),

VF = (a2r - + a- (a, r) + 2a + 2r)r - (a + 2r + a- (a, r)){

1 + a 4a, r) + 2a + 2rjr -la + 2r + ч ул, i jp0, a = r - r0, a = a, r = r, \n\ = 1.

Источники магнитной активности головного мозга моделируются точечными токовыми диполями. Каждый такой диполь характеризуется двумя векторами: r0 - радиус-вектор диполя (положение диполя) и Q - вектор, задающий направление и силу диполя (момент диполя). Для вычисления магнитной индукции на поверхности головы используется модель токового диполя в проводящей сфере.

На рис. 1 и 2 приводятся результаты аппроксимации пространственного распределения биомагнитного поля сферическими гармониками.

Результаты тестирования предлагаемых алгоритмов описания функций на сфере

Ниже в качестве примера рассматриваются модельные данные, полученные в расчетах из магнитной энцефалографии (МЭГ) [15]. Данное направление исследований -интенсивно развивающаяся область экспериментального изучения высшей нервной деятельности человека, функциональных областей мозга и диагностики различных патологий. Это обусловлено возможностью неинвазивного получения данных о процессах, происходящих как в коре головного мозга, так и в глубоких его отделах. Чрезвычайная слабость как спонтанных, так и вызванных магнитных полей, возбуждаемых токовыми источниками в мозге, накладывает повышенные требования к используемой измерительной аппаратуре и методам обработки полученных данных. Несмотря на возникающие технические трудности и высокую стоимость необходимого экспериментального оборудования, использование данного подхода как дополнения или альтернативы электрической энцефалографии оправдано следующим обстоятельством. Магнитное поле, в сравнении с электрическим, испытывает значительно меньшие искажения на внутричерепных неоднородностях и покрывающих тканях, что существенно повышает точность локализации источников и снижает требования к знанию структуры внутричерепной среды.

Основные теоретические аспекты данной проблемы изложены в работе [16]. Величина измеряемого магнитного поля, если датчик находится в точке r и имеет направление нормали n, определяется следующим образом:

10 п

Аппроксимация функции BB(fi,teta) сферическими гармониками, (fi= pi)

, Km ФГ №1

4 -2 -0 -

1,57

3,14

BBfi,teta)+noise(l0%) -vost BB(fi,teta)

Аппроксимация функции BB(fi,teta) сферическими гармониками, (fi= 0)

0 T-

1 57

.3,14

BBf,teta)+noise(10%) -vostBBfi,teta)

Рис. 1. Аппроксимация функции магнитного поля сферическими гармониками. Уровень шума, наложенного на функцию равен 10%, глубина аппроксимации равна 4. Изображение кривой BB(0,ф) при фиксированной переменной ф : а) ф=п ;

б) ф = 0 .

Ошибка аппроксимации составляет 0,045 и 0,038, соответственно.

Аппроксимация функции магнитного поля 2i сферическими гармониками. (teta= pi/25)

I 1,57

6,28

BBfi, teta)+noise(10%) vost BBfi,teta)

0--гт -2 У

-4-У

Аппроксимация функции BB(fi,teta) сферическими гармониками. (teta= pi/6)

А<* v.

1,57

3,14

4,71

6, 28

\i kt

BBfi,teta)+noise(10%) -vostBB(fi,teta)

Рис. 2. Аппроксимация функции BB(0 ,ф ) сферическими гармониками.

Уровень шума, наложенного на функцию равен: а) 20%, б) 10%. Изображение кривой BB(0 ,ф) при фиксированной переменной 0 : а) 0= П /25;

б) 0=п/6.

Погрешность аппроксимации равна: а) 0, 0942 и б) 0,043. На рис. 3 представлены результаты вычисления первой производной функции магнитного поля, при чем на функцию наложен шум (10%). Производная функции получена тремя коэффициентами, ошибка составляет: а) 0,176, б) 0,164.

12 п 10 1

0 --

-2 0

-4 -

-10 -12

Вычисление производной функции BB(fi,teta), fi= pi)

i i Д AfA .iinM.ii t.M:i-y iff ifck

1,57

-BB(fi,teta)c шумом(10%)

-vost BB(fi,teta)

- вычиcл.производн.BB(fi,teta)

10 -8 -6 -4 2

0 3,14 -2 -4 -6 -8 -10 -12 j

Вычисление производной функции BB(fi,teta), (fi= 0)

1,57

3,14

teta

BB(fi,teta)c шумом(10%) vost BB(fi,teta)

вычиcл.производн.BB(fi.teta)

Рис. 3. Вычисление производной функции по выведенным формулам.

а) при ф= п ,

б) ф = 0 .

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

Развитые подходы будут реализованы в задачах спектральной диагностики и классификации биомедицинских систем [17].

Работа выполнена при поддержке РФФИ, проекты 04-02-17368, 04-01-00756.

Литература

1. Бейтмен Г., Эрдейи Э. Высшие трансцендентные функции. Справочная математическая библиотека. Том 2. - М: Наука, 1974.

2. Никифоров А.Ф., Уваров В. Б. Специальные функции математической физики. Учебное пособие, Гл. ред. Физ.-мат. лит.,- М.: Наука, 1984, 344с.

3. Джеффрис Г., Свирлс Б. Методы математической физики. - М.: Мир, вып. 3, т. 3,

345с.

4. Стейн И., Вейс Г. Введение в гармонический анализ на Евклидовых пространствах.

- М.: Мир, 1974, 336с.

5. Виленкин Н.Я. Специальные функции и теория представления групп. - М.: Наука, 1991, 576с.

6. Никольский С.М. Приближение функций многих переменных и теоремы вложения.

- М.: Наука, Гл. ред. Физ.-мат. лит., 1969, 480с.

7. Magnus W., Oberhettinger F. F0rmulas and The0rems 0f MathematiCal PhysiCs, 1949, pp.54.

8. Abram0witz M., Stegun I.A. Handb00k 0f MathematiCal FunCti0ns, Applied MathematiCs Series, V. 55, Chapter 8. (Специальные функции. Под ред. Абрамовица и Стиган. - М.: Мир, 1982).

9. Press W.H., Teuk0lsky S.A., Vetterling W.T., Flannery B.P. NumeriCal ReCipes in C. The Art 0f SCientifiC C0mputing. - Cambridge University Press, 1992.

10. Курант Р., Гильберт Д. Методы математической физики, т. 1, М.-Л., Гостехиздат,

1951.

11. Л.И. Куликова, С. А.Махортых. Выбор оптимальных условий аналитического описания контурных объектов. Электронный журнал Исследовано в России , 144, стр. 1501-1511, 2005г. http: zhurnal.ape.relarn.ru/artiCles/2005/144.pdf

12. Дедус Ф.Ф. Аналитическое представление экспериментальных данных и их обработка. Кибернетика и вычислительная техника. Вып.74, Наукова думка ,

Киев, 1987.

13. Dedus F.F., Makh0rtykh S.A., and Ustinin M.N. AppliCati0n 0f the Generalized SpeCtral-AnalytiC Meth0d in Inf0rmati0n Pr0blems Pattern ReC0gniti0n and Image Analysis,

V0l. 12, N0. 4, 2002, pp.429-437.

14. Дедус Ф.Ф., Махортых С.А., Устинин М.Н., Дедус А.Ф. Обобщенный спектрально-аналитический метод обработки информационных массивов. Задачи анализа изображений и распознавания образов. М.: Машиностроение, 1999, 357с.

15. Sarvas J. BasiC mathematiC and eleCtromagnetiC C0nCepts 0f the bi0magnetiC inverse problem. Phys. Bi0l., V0l. 32, N 1, 1987, pp. 11-22.

16. Устинин М.Н., Махортых С.А., Молчанов А.М. и др. Задачи анализа данных магнитной энцефалографии. В кн. Компьютеры и суперкомпьютеры в биологии. М.: Институт компьютерных технологий, 2002, С. 327-349.