Методы и алгоритмы разделения смеси сигналов. Ii. Применение м-градиента к анализу независимых компонент Текст научной статьи по специальности «Электротехника, электронная техника, информационные технологии»

ISSN 0868-5886

НАУЧНОЕ ПРИБОРОСТРОЕНИЕ, 2009, том 19, № 4, с. 83-95 ОБРАБОТКА И АНАЛИЗ СИГНАЛОВ =

УДК 681.51; 621.391; 519.21 © В. М. Малыхин, А. В. Меркушева

МЕТОДЫ И АЛГОРИТМЫ РАЗДЕЛЕНИЯ СМЕСИ СИГНАЛОВ. II. ПРИМЕНЕНИЕ М-ГРАДИЕНТА К АНАЛИЗУ НЕЗАВИСИМЫХ

КОМПОНЕНТ

Анализируется задача разделения смеси сигналов (с восстановлением вида ее компонент) при отсутствии информации о пропорциях и типе смешивания. Метод основан на использовании информационного критерия и адаптивного алгоритма для обучающейся нейронной сети. Рассмотрены несколько видов распределения первичных сигналов, поступающих на сенсоры информационно-измерительной сети, и соответствующие им изменения функции преобразования нейронов. Подход к разделению смеси сигналов включает применение модифицированного градиента (м-градиента) в схеме анализа независимых компонент. (Статью I цикла см. "Научное приборостроение", 2009, т. 19, № 2).

Кл. сл.: смесь сигналов, методы разделения, адаптивный алгоритм, нейронная сеть, м-градиент, анализ независимых компонент

ИСПОЛЬЗОВАНИЕ М-ГРАДИЕНТА И ОБУЧЕНИЕ НЕЙРОННОЙ СЕТИ (НС)

Простая схема линейного смешивания сигналов (рис. 1) в векторно-матричной форме описывается соотношением (1):

x(k) = Hs(k) + v(k),

(1)

где х^) = [х^),...,хт^)]т — зашумленный сигнал на сенсорах; s(k) = [¿^),..., sn ^ )]т — первичный сигнал (ПС); у^) = ),..., ут ^ )]т — вектор шума; Н — неизвестная смешивающая матрица размера т х п и полного ранга. При этом предполагается, что доступным является только вектор-сигнал х(^), регистрируемый на сенсорах информационно-измерительной сети (ИИС).

Цель — разделение сигналов смеси (РСС) — достигается построением НС с прямым распространением сигнала или с обратными связями (так называемой рекуррентной НС) и соответствующего этим сетям адаптивного обучающего алгоритма, который позволяет оценивать компоненты первичного сигнала s(k), идентифицировать смешивающую матрицу (СМ) Н или разделяющую матрицу W.1 Свойство адаптивности алгоритма состоит в возможности отслеживать изменяющуюся во времени систему сигналов.

Сначала выполнен анализ обучающих алгоритмов для простого (мгновенного) смешивания сигналов при известном и одинаковом числе п ПС и

1 Более полно описание задачи разделения смеси сигналов содержится в [1].

сенсоров (т). При этом Н и W — несингулярные матрицы размера п х п . Предполагается, что первичные сигналы ¿¡(к) взаимно независимы и имеют нулевые средние значения (т. е. Е^(£)}=0), а аддитивный шум \(к) пренебрежимо мал или становится малым после предварительной обработки. На следующих этапах эти ограничения будут постепенно сняты.

Критерий Кульбака—Лейблера как мера статистической независимости

Для получения хорошей оценки у = W х ПС s используется "целевая функция" р^,1^ . Ее математическое ожидание, называемое функцией риска R( = Е{р(у,, представляет меру взаимной независимости выходного сигнала у(^). Функция риска R(W) достигает минимума при условии, что компоненты сигнала у становятся независимыми, а матрица W (с точностью до масштабирования и перестановки) становится равной обратной СМ Н-1. Для достижения такого минимума используется критерий Кульбака—Лейблера [2], который служит мерой независимости. Если Ру (у^) — плотность вероятности (ПВ) стохастической ("случайной") переменой у = W х = = W Н s и q(y) — ПВ переменной у, у которой все компоненты статистически независимы и, слеп

довательно, q(y) = ^ q¡ (у1) (т. е. представима

1=1

произведением ПВ компонент qi (у{)), то q(y) — служит в критерии плотностью для сопоставления.

Неизвестная часть структуры

Рис. 1. Блок-схема, иллюстрирующая анализ независимых компонент и разделение сигналов смеси

Критерий Кульбака—Лейблера (ККЛ) дает меру отличия ПВ р (у^) (т. е. ПВ, где у получено с

действительной величиной W) от ПВ ч(у), предполагающей независимость компонент.

Таким образом, функция риска R(W) =

= Е{р(у^)} представима в виде ККЛ Крч (W):

Я( W) = Е{р(у^)} = Ки (W) =

= К

Ру (у; W)||

¿у

(2)

11«(У )1=ГРу (у^^

■> ч (у)

Следовательно, Я(W) = Крч (W) показывает, насколько ПВ ру (y;W) далека от сопоставляемой ПВ ч(у). Если ПВ ч(у) равна истинному распределению р,, (s) первичных сигналов, Ру (У^) = р% (у) при W = Н-1, поскольку в этом

случае у = W Н s = s.2

Для задачи анализа независимых компонент принимается, что ч(у) является произведением ПВ независимых переменных уд(у) = р (у) =

= ПР'(у>), где р,(у,) — ПВ у (/ = 1, ..., п).3 По-

¿=1

этому мера независимости принимает форму выражения (3):

Крч = E{р(y,W)} = | ру (у)^

р(у) ГГч , (у,)

-¿у,

(3)

где р(у^) = log(ру (у) / ч(у)).

ККЛ имеет также трактовку взаимной информации:

п ™

Крч = -н (у) - X} ру (y)log ч, (у, )ау,

(4)

где дифференциальная энтропия выходного сигнала у = W х определяется соотношением

Н(у) = -1 ру (у) ^ ру (у^у .

(5)

2 Использование Крд (W) в качестве функции риска R(W) имеет дополнительное преимущество. Даже если д(у) не равно р, , минимизация R(W) все равно определяет W = Н-1. Кроме того, большую часть обучающих алгоритмов, не использующих ККЛ явно, можно трактовать на основе функции риска в виде Крч (W) (отличие состоит только в выборе ПВ для сопоставления д(у)). Среди них — алгоритмы на основе максимизации энтропии [3], анализа независимых компонент [4] и анализа максимального правдоподобия [5].

3 Строго говоря, р. (у,) — это маргинальная ПВ для у,:

р. (у,) = 1-™ ру (у) •¿у, где интегрирование осуществля-

ется по у, [ у1,..., у,-1, ум,..., уп ]т (т. е. из у удаляется пе-

ременная у,).

Когда qi (yi) = р,( yi), то (учитывая, что ¿у = ¿у iдLyi) второе слагаемое в (4) может быть выражено через маргинальную энтропию Н(у):

+да +да +да

| ру ^^ р,(у, Му = | log р,(у,) | ру ДОФ '¿уг =

-да -да —да

+да

= | р,(у,^р,(у,щ=щ^ср(у,))}=-н(у,). (6)

—да

Следовательно, критерий КЛ может быть выражен разностью между энтропией Н(у) и маргинальной энтропией Н(у):

Кр9 = Е{р(у^)} = -Н(у) + £Н, (у,) .

(7)

Для у = W х дифференциальная энтропия выражается соотношением4

Н(у) = Н(х) +

На основании использования правила соответствия ПВ линейному преобразованию W.

1=1

где H(x) = -1 px (x) log px (x)dx не зависит от мат-

—да

рицы W. С использованием этих соотношений функция риска для общего вида ПВ q приобретает форму

R W) = Kpq = E{p(y,W)} =

= — H(x) — log|det(W)| — ]Г E{log q, (y,)}; (8)

1=1

причем H(x) можно опустить, т. к. это слагаемое не включает зависимости от W.

Основные правила обучения, основанные на использовании м-градиента

Для получения матрицы W (обратной относительно СМ H-1) требуется минимизация функции риска (8). Поэтому обучающий алгоритм реального времени (РВ) основывается на методе градиентного спуска:

AW(k) = W(k +1) — W(k) = — ri(k)

dp[y(k), W]

aw

где ц(к) — скорость обучения, др/дW — (п х п) -

матрица градиента с элементами дрдwij .

Вычисление матрицы градиента (с учетом (8)) приводит к соотношению

AW(k) = ^(к) |^-т (к) - f [у(к )]хт (к)],

где f (У) = [/1( Jl), /2(У2), . . ., /п (Уп)]т — вектор-столбец с 7-й компонентой, которая равна

f (Уг) = —

= — g,(У г )

g г (Уг У

dl0g g г (Уг ) d g г (Уг ) / 4Уг

dy

g г (Уг )

где gi(Уi), г = 1,...,n — приближенные (обычно параметрические) модели ПВ первичных сигналов {s,}. Обычно —dp/aw представляет направление максимального уменьшения функции p, но если пространством оптимизируемых параметров является множество матриц {W}, то, как показано Амари, Кардосо, Ченом, Сичоки (Amari, Cardoso, Chen, Cichocki) [6, 7], такое направление определяет модифицированный градиент (м-градиент), который представляется соотношением5

dp (y, W)

aw

WT W = [I — f (y)y T]W.

Следовательно, алгоритм обучения представляется в следующей форме:

AW(k) = — r(k)

dp[y(k), W] wt

aw

wt w =

= r(k)[ I — f [y(k)] y T(k )]W(k).

(9)

Другой подход к вычислению м-градиента для р(у^), выраженного соотношением (8), состоит в использовании дифференциала р по W:

dp = p(y,W + dW) — p(y, W) = Хг

3 dw 3

Строго говоря, множество матриц {W}является мультипликативной группой (с единицей I). В работах [6] и

[7] для такого пространства матриц {W} введена "естественная" метрика, которая порождает м-градиент. Иногда м-градиент называют также естественным.

Это выражение представляется в виде dр = - tr(dWW-1) + f т(у^у, где ^ — след матрицы; Д(у) — вектор-столбец с компонентами

/ ■7 (У г ) = -Чг( У г )/ Ч 7 (У г ) [7] Для оценки ПС

у = Wx справедливо соотношение dy = dW • х = = dW • W-1y. Следовательно, можно обозначить dX = dW • W-1, и компоненты йг.. являются линейными комбинациями базиса { dw . } .

В базисе dX м-градиент представлен теперь выражением dр = - ^ ^Х) + Дт (у^Х • у, которое приводит к алгоритму обучения

АХ(к) = -п(к)%Р = Л(к)[ I - Д (у (к ))у т(к)]. (10) dX

С учетом АХ(к) = AW(к) • W-1 (к) (10) преобразуется в алгоритм обучения для матрицы W:

W(k +1) = W(k) + П(к)[I - Д[у(к)]ут (к)]W(k). (11)

Применение м-градиента в соотношении (8) приводит к алгоритму обучения с усреднением (символ ( )):

дК

AW(k) = -П(кWтW = П(к)[I - Д(у)ут

а представление его в РВ (опуская операцию усреднения) показывает эквивалентность этого алгоритма с полученным ранее выражением (9).

ОБОБЩЕНИЕ АЛГОРИТМА, ОСНОВАННОГО НА М-ГРАДИЕНТЕ

Чтобы уменьшить элемент неопределенности масштаба (свойственной задаче РСС), в алгоритме

Правило обучения, представленное выражением (9), введено в [8], а м-градиент использован несколько ранее в [9].

обучения, основанном на м-градиенте, используются ограничения на величину восстановленных сигналов смеси: ЕЩуг)}=1. Однако если первичные сигналы (ПС) нестационарные и их средние значения изменяются достаточно быстро, эти ограничения требуют столь же быстрого изменения значения "разделяющей" матрицы W. Избежать этой трудности позволяет использование особой формы ограничений: направление изменения матрицы W должно быть ортогонально классу эквивалентности матриц, образованному за счет масштабирования. Такой алгоритм, определенный соотношением (12), адаптируется к быстрым изменениям величины ПС:

AW(k) = ф) [A(k) - f [y(k)] yT (k )]W(k), (12)

где Л = diag{^,X1,..,Xn} — положительно определенная матрица и Xi = f (yi)yi.

В другой модификации алгоритм обучения имеет вид

AW(k) = ^(k)^(k) - y(k) gT( y(k))] W(k), (13)

где использована видоизмененная Атиком и Ред-лихом (Atick, Redlich) [10] форма м-алгоритма. В этом случае нелинейность gt (yi) является обратной относительно функции f (yi): gi = f_1.

Еще более общая форма алгоритма обучения определяется соотношением (14) [11]:

AW(k) = V(k) • [Л(k) - а • (f [y (k)] y T(k)) +

+ßiy(k)gT[y(k)]yT(k)) ] W(k), (14)

где Л — диагональная матрица, которая предназначена для нейтрализации диагональных элементов всей матрицы в правой части соотношения (12), а а и ß — адаптивно определяемые параметры [7].

М-градиент при ограничениях на ортогональность

Можно считать, что над вектором-сигналом x(k) = Q H s(k) = A s(k) уже выполнена предварительная обработка отбеливания и сигнал нормализован [1], так что

Rxx = E{x(k)xT(k)} = Im ; Rss = E{s(k)sT(k)} = In .

(15)

При равном числе ПС и сенсоров (п = т) из (15) следует, что ААт = 1п и СМ А = QHeRпхп является ортогональной. В этом случае для восстановления ПС у = Wx можно рассматривать только ортогональные матрицы W (такие, что W-1 = Wт).

В рассмотренном выше соотношении ¿X = = dW • W-1 = dW • Wт базисные функции dW являются косо-симметричными (т. е. ¿X = -dXT). Это видно из цепочки равенств

0 = ¿I = d(WWT) = dW • Wт + W • dWT = ¿X + ¿Хт.

Отсюда следует, что м-градиент бр/ бХ = = (бр/бW )• Wт также косо-симметричный, и поэтому

(бр/бХ)• Wт = %)ут -yfт(у).

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

W(k +1) =

dp

= W(k)-ц(к) • ^^р(у,^ = W(k)-ц(к) • б^. (16)

Матрица W является (приблизительно)6 ортогональной на каждом шаге итерации. В этом случае х = Wт у, и алгоритм (14) приобретает более простую форму:

W(k +1) = W(k) - ц(к) х

х\_/(у (к)) • у т (к) - у(к) • / т (у(к ))]• W(k). (17)

Практически благодаря косо-симметричности члена [/(у(к)) • ут(к) - у(к) • /т(у(к))] в (17) выполнение обработки смеси сигналов, связанной с отбеливанием (декорреляцией), может производиться одновременно с РСС7, а сам алгоритм принимает вид

ДЩк) =ц(к) х

х[1 - y(k )yT(k) - f (y(k)) • yT(k) + y(k) • f T(y(k)) ]x xW(k).

(18)

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

Эта матрица точно ортогональная для алгоритма с непрерывным временем. Для дискретного времени реализуется почти точная ортогональность в том случае, если параметр скорости обучения п будет иметь достаточно

малое значение.

В связи с этим алгоритм (17) условно назван "легким" [9].

AW(k) = r(k) х

x[l — (y(k )yT(k)) — (f (y(k)) • yT(k))

+

+

(y(k) • f T(y(k ))) ]• W(k).

(19)

K (W ) = K p (y, W )

"p^'V^ " e'llq(y)

Г p(y,W)

= fРу (y,We )log Р(У; ) e) dy,

q(y)

(20)

= — r(k) •

dp(y, We)

a w.

— w.

ap(y, we)

a w.

w.

(21)

ВЫДЕЛЕНИЕ ГРУППЫ ПЕРВИЧНЫХ СИГНАЛОВ ИЗ ИХ СМЕСЕЙ, РЕГИСТРИРУЕМЫХ СЕНСОРАМИ ИИС

Алгоритмы, использующие метод м-градиента вместо "простого" последовательного РСС, позволяют осуществлять выделение произвольной группы ПС. Положим, что число сенсоров т больше (часто неизвестного) количества п первичных сигналов, сигналы сенсоров отбелены и требуется выделить группу определенного количества е первичных сигналов. Тогда сигналы сенсоров удовлетворяют условию Rxx = Е{ххт} = \п (где

х = х1 = Qx ); общая смешивающая (п х п)-матрица А = QH является ортогональной (т. е. А Ат = Ат А = In), и для РСС при е = п разделяющая (п х п) -матрица имеет вид W = А-1 = Ат .

Для выделения е первичных сигналов используется подходящая форма ККЛ:

При этом если в начале вычислений выполнено ограничение We (0) • Weт (0) = !е, то оно сохраняется на каждом шаге итерации, т. е. [11]

^^е (к) • ^^е^ (к) = I, .

Соотношение (21) соответствует правилу обучения в виде (22) с начальным значением We (0), удовлетворяющим условию We(0) • Weт (0) = !е :9

где ч(у) = ^ е= Ч (У') представляет ПВ группы соответствующих независимых сигналов на выходе структуры РСС. Следовательно, с учетом (8) минимизируемая "функция стоимости" имеет вид

р(у ) =-Х Ч (У) при ограничении WeWeт = Ie (или эквивалентно: wг.wт = 8.), где We eRехп — разделяющая матрица с е < п , w i — 7-я вектор-строка матрицы We .8 Для того чтобы эти ограничения выполнялись в процессе обучения, форма м-градиента в алгоритме выделения группы ПС должна использоваться в следующем виде:

А^^е (к) = ^^е (к + 1) - ^^е (к) =

We (k +1) = We (k) = — r(k) x

x { f [y (k)] • xT (k) — y (k) • f T [y (k)] • We (k) }. (22)

МОДЕЛЬ ОБОБЩЕННОГО ГАУССОВА РАСПРЕДЕЛЕНИЯ ВЕРОЯТНОСТЕЙ ДЛЯ ПС И РАЗДЕЛЕНИЕ СМЕСИ СИГНАЛОВ ПО СХЕМЕ АНАЛИЗА НЕЗАВИСИМЫХ КОМПОНЕНТ (АНК)

Оптимальная нелинейная функция Д (у), входящая в обучающий алгоритм (9) и определенная ранее выражением Д (у) = [/1( уД /,(У2),..., / (Уп )]т, где / (У г) = - £'(У г)/£ г (У г), требует знания распределения ПВ первичных сигналов, которое обычно бывает неизвестным. Каждая модель ПВ соответствует некоторому виду функции Д(у), которая в терминах нейронной сети, реализующей алгоритм обучения, называется активационной функцией. Для некоторых типов моделей ПВ выражение для активационной функции приведено в таблице. Например, для косинус-гиперболического распределения (имеющего "надгауссову" форму ПВ)10 активационная функция имеет вид / ( У ) =

= tgh(yiу1), где ^ обычно имеет значение 1/с2 .

Вид ограничения связан с тем, что смешивающая матрица А = QH является квадратной ортогональной, а разделяющая ("размешивающая") матрица We после выделения е первичных сигналов должна удовлетворять соотношению We А =[!е, 0п-е ].

При e = n (т. е. при "полном" РСС) разделяющая матрица W = W является ортогональной, обучающее

правило упрощается и сводится к алгоритму (17), предложенному Кардосо и Лахэльдом (Cardoso, Laheld ) [9]:

W(k +1) = W(k) = — r(k) x

x { f [y(k)] • yT (k) — y(k) • f T [y(k)] } • W(k).

10 Термин "надгауссова" ПВ используется для распределения с эксцессом, большим чем у гауссовой ПВ, поскольку вершина такого распределения выше гауссова пика. В таком же смысле далее говорится о "подгауссо-вой" ПВ, если величина ее эксцесса меньше, чем у гауссовой ПВ, и, следовательно, вершина этого распределения ниже, чем у гауссовой ПВ.

T

Типовые распределения ПВ д(у) и соответствующие нормализованные активационные функции

/ (у) = - dlog ч( у)! ¿у

Название ПВ

ПВ-функция ч(у)

Активационная функция /у)

Гаусса

(V ч/2жст)- ехр (-| у|2 ¡2а2)

у/ст2

Лапласа

(12ст)-ехР (-| у\/ст)

sign( у)/ст

Коши

(1 ст

1 + (у/ст)

2 у

ст2 + у2

Треугольное

ст(1 -1 у\/ст)

sign( у)

ст(1 -1 у\/ст)

Косинус-гиперболическое

1/ (ж cosh( у / ст2))

tgh (у/ст2)

Экспоненциально-степенное (или обобщенное Гаусса)

2ст Г(1 / г)

• ехр

Г 1 У_ г Л

г V ст /

г > 0

у

ст'

• sign( у)

1

2

1

г -1

г

0.8

0.7

0.6

0.5

0.4

0.3

0.2

0.1

1 \

г = 1

\

/ г = 2

Г/ \\ г = 6

№

¡у

у

Экспоненциально-степенная ПВ (при различных соотношениях величины ее параметров о и г ) представляет целое семейство распределений, отличающихся формой и включающих распределения Лапласа (при г = 1), Гаусса (г = 2) и равномерное (при г >5 ^ 6) [12], [13]. Таким образом, семейство экспоненциально-степенных ПВ (ЭС_ПВ) включает над- и подгауссовы формы распределения (рис. 2).

Обучающий алгоритм (9), включающий нелинейные функции / (у,,) = - g'(у,)/g, (у,), при использовании ЭС_ПВ для ПС адаптируется к форме распределения и определяется его параметрами. Для ЭС_ПВ в виде

ч,(у,) = -

-ехр

А г Л

V г ' /

(23)

Рис. 2. Экспоненциально-степенная ПВ при различных значениях параметра г

2ст, Г(1/ г)

Íда

^ у-1 ехз(-у)ф—

гамма-функция, стг = Е{| у,|г} — обобщенная мера дисперсии ЭС_ПВ, параметр г, принимает значение в диапазоне (0 ^ да).

Оптимальные нормированные нелинейные ак-тивационные функции выражаются соотношением

/ (у,) = - d log ((у,) )/d у, = I у Г ^п( у,),

0

2

1

0

1

2

и т. к

■. к. sign(y) = У/|уI то f (У) = У/|У

|2—r 11

Для оценки значений f (y,) = y/|y |2 r = = yt (k)/cf(2 r) (k) может использоваться также

скользящее среднее, т. е. в качестве ст берется оценка

(2—Г) _

= тп

?(2—r)

(k +1) = (1 — ) + 1 У (k)

I (2 r )

= 1 (k) -

r 2(k)

ция f (У) = у/(|У|

cc,2(k)

f; =

где

\—f (У) Уj для Y,J <1,

[—У if J (У;) для других случаев, у, = cfJ Ef y )}E{f;( y)} и

/(У) = У /(1У|2-Г +е)

при значении г г между нулем и единицей.

В качестве другого подхода после некоторых упрощений можно использовать правило обучения (24) с диагональной матрицей Л(к) =

= diag{F (у (к)) • gт (у(к))} . Активационные функции этой матрицы адаптивно меняются со временем, устойчивы к спайк-шуму и имеют вид:

Параметр r, может быть фиксированным (если доступно знание статистики ПС) или получен в процессе обучения нейронной сети, реализующей РСС. Правило подстройки параметра r, основанное на градиенте, имеет вид:

Ar(k) = —i,(k) Ф/ a r, = i(k) •a log (q( у г Va r, =

0.1r (k) + | y, (k )|r(k) (1 — log (I y (k )|r(k)))

f, (У) =

g, (У г ) =

для k4( У{) > S0,

для других случаев;

для k4( У, ) >— S0, для других случаев,

(25)

В случае предельно обостренной ПВ (т. е. для ПС в виде так называемых спайк-сигналов) и Г = 0 в качестве оптимальной получается функ-

У (к)

предложенная ра-

нее Матсуокой (Matsuoka) [14]. Аналогично для ПС с плотностью распределения Лапласа (т. е. для

г =1) /(у,) = У/(ы)=у,(кV(к), а для больших значений г, (г = 5 ^ 10) и (почти) равномерной ПВ — / (У) = У (к )(| У (к)|8) = У (к) • с?(к).

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

AW(k) = Е(у^(к), (24)

где элементы матрицы Е(у) определены соотношениями

где к4(yi) = Е{уг4}/Е2{уг2} -3 — это нормализованное (относительно гауссова) значение эксцесса распределения плотности вероятности ПС, 80 — малое пороговое значение.

Алгоритм обучения (24) и (25) наблюдает и оценивает статистики каждого выходного сигнала, зависящие от знака или величины нормализованного эксцесса (который является мерой отличия ПВ сигнала от гауссовой ПВ). Затем алгоритм автоматически отбирает (переключает) подходящие нелинейные активационные функции так, чтобы успешно разделять все негауссовы ПС.

Подобные методы применимы к другим ПВ. Так, для обобщенного распределения Коши, определяемого тремя параметрами г > 0, vi > 0 и с плотность вероятности имеет вид:

2 1 '

Я,(У) = -

B(r, V,)

+1/r,

(26)

¡1 + (1/у.)

У A(r),

где

A(r) = [с,2 Г(1 / r )/Г(3/r )]1/2;

B(r,v,) = r v—Уr Г(г, +1 / r) / 2A(ri)Г(у,) Г(1 / rt).

При этом активационная функция выражается соотношением

f, (У) =

(rv +1)

(v, | A(r) |r, +1У |r)

lylr—1 • sign(y). (27)

В особом случае импульсивных (или "спайк-") сигналов параметр г может принимать значение между нулем и единицей. В этом случае целесообразно принимать несколько модифицированную активационную функцию в виде /(у) = у/[|У |2~г + е], 0 < г < 1, где е — малое положительное число (~10-3).

Статистические моменты и эксцесс ЭС_ПВ

Более полное представление об ЭС_ПВ дают статистические моменты этого распределения, особенно второй и четвертый, которые определяют оценку эксцесса. Общее выражение п-го мо-

^да

упр(у; г^у , и в силу

-да

симметрии ПВ при нечетном п тп = 0 (в том числе среднее т1 = 0). Вычисление четных моментов, которые полностью характеризуют ЭС_ПВ, основано на формуле [12], [15]:

Г уп-1е^у° ¿у =1 ц-Уп Г[ - |.

Второй момент ЭС_ПВ определяется выражением

^да да

у2 Р( y, г )4у = 21 у

-да J 0

= 21 у2-г-

0 2ст Г(1 / г)

Iу

е ¿у =

ст Г(1 / г)

р да -[у)

]0 у2е ¿у

которое подстановкой у = г • ст приводится к вы-

ражению т2 = -

гст

-Г г Vгг ¿г = ст2^^)

I .10

Г(1/ г) Г(1/ г)

Подобным путем определяется момент четвер-.4 Г(5 / г)

того порядка т4 = ст

Г(1/ г) „ Г((2к +1)/ г)

и момент 2к-й степе-

2к Г(1/ г)

Безразмерная (в отличие от моментов) величина эксцесса является показателем уплощенности или обостренности ПВ. "Нормализованный" эксцесс к(у) показывает его превышение над гауссовым и определяется выражением: к (у) = = т4/(т2)2 -3 .

Метод РСС на основе анализа независимых компонент (АНК) при критерии эксцесса ЭС_ПВ ПС

Модели ЭС_ПВ (с учетом данных таблицы и комментария к формуле (22)) соответствует нелинейная функция алгоритма (18), имеющая вид

/(у,) = -ч, (у, Vdу, =1 у,|г-1 ^ёпСу). при г = 4

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

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

несколько видов нелинейных функций (активаци-онных функций системы РСС, реализуемой на основе нейронной сети). Кроме того, при низких значениях к(г) (для подгауссовой ПВ) его величина очень слабо зависит от г, она практически постоянна для г в пределах от 4 до 12. Поэтому если оценка для к(г) отрицательна, то можно использовать г, = 4.

При положительном значении к(г) зависимость от параметра г достаточно значительна, так что можно предполагать несколько возможных значений этого параметра. Однако опыт применения метода показывает, что для успешности РСС бывает достаточно использовать 2-3 значения параметра г для всего практического диапазона положительных значений эксцесса (т. е. у сигналов с надгауссовой ПВ).

Из числа применяемых вариантов нелинейных функций для алгоритма РСС по методу АНК [15] целесообразно проанализировать устойчивость алгоритма (18) для различных случаев: 1) г, = 4 при к, < 0; 2) г, = 1; 3) г, = 0.8 для к, > 0 [15].

1) г, = 4 — подгауссовы распределения. Выбор г, = 4 предполагает у ПВ ПС наличие отрицательного эксцесса (к, < 0) и соответствует кубической нелинейной (активационной) функции, т. е.

/ (у,) =1 у,|2 у,.

2) г, = 1 — распределение Лапласа. При выборе г, = 1 ЭС_ПВ переходит в распределение Лапласа, т. е. qi (yi) = (1/2ст,.)ехр -1 yi / ст{ |, которое соответствует "ступенчатой" нелинейной функции / (уг) = sign(уг) = уг /1 уг | (активационная функция с жестким ограничением).

3) г, < 1 — распределение с отчетливо выраженным пиком. При резко выраженном пике ПВ (большом значении эксцесса, к >> 1) выбирается г, < 1. Это соответствует невозрастающей нелинейной функции с сингулярностью в окрестности нуля. Для практических приложений / (у{) принимается постоянной (и ограниченной) величиной в пределах интервала [-8, е], который определяется малым числом е.

Дисперсия у, для ЭС_ПВ определяется соотношением Е {у,2} = ст,2 • Г(3 / г) / Г(1/ г), и при ограничении нормализацией Е{ у,2} = 1 дисперсия имеет величину ст1 = [ Г(1/ г) / Г(3 / г )]1/2. Условием устойчивости алгоритмов обучения (11) и (18), предложенным Амари, Ченом и Сичоки [7] и Кардосо, Лахэльдом [9], является проверка условия + X, > ^ где = ЕШ у г )} - Е{/ ( Уг ) у г }.

Для каждой отдельной компоненты у, восстановленного ПС достаточным условием устойчивости является х^ > 0. Вне области у{ е [-£,£] вычис-

а

а

г

ношениям:

ЕЩ Я)} =

ление Е{/¡(yi)} иЕ{^(yi)yi} проводится по соот- где d = 10т2т4 -18т3 -12т32; т2, т3,т4 — вто-12 рой, третий и четвертый статистические моменты

у. Преимущество в использовании системы ПВ Пирсона в том, что она позволяет эффективно реализовать РСС при несимметричных распределениях.

^+да

(Т - 2)|у,Г

-да

2СТ Г(1/ Т)

еХР- (| У, Г /п,Т) • ^г =

( Vл) \ Г[(г-1)/г]; Г(1/ г)

Е{£ (У,) У,} =

У,1 У, Г1^У,) •

-да

П Г(1/Г)

ехР- (| У,\Т /пт) • 4у =

ГП

Г(1/Т )

Г[(г +1)/гг ].

Модель ПВ Пирсона

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

Широкий класс ПВ (симметричных и несимметричных) может моделироваться системой распределений Пирсона, которая описывается дифференциальным уравнением

Ч( У) =

(У - а)д(у) Ь0 + Ь1 у + Ь2 У2

где а, Ь0, Ь1, Ь2 — параметры, зависящие от ПВ оцененных ПС. Оптимальная активационная функция для системы ПВ Пирсона имеет вид

/ (у) = - Я (у v я(у) = (а - У)/(ьо + Ъ У + Ь2 У2). К разновидностям форм модели Пирсона относятся широко используемые ПВ, в том числе бета, гамма, нормальное (Гаусса), Стьюдента и некоторые другие специальные виды распределений. Параметры а, Ь0, Ь1, Ь2 оценивают непосредственно методом моментов:

а = Ь = -

т3(т4 + 3т2)

Ьо =-

т2(4т2т4 - 3т^)

d

Ь2 = -

(2т2т4 - 3т32 - 6т3) d :

АЛГОРИТМЫ НА ОСНОВЕ М-ГРАДИЕНТА ДЛЯ РСС НЕСТАЦИОНАРНЫХ ПЕРВИЧНЫХ СИГНАЛОВ

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

Основным предположением при РСС является статистическая независимость ПС. При этом если ПС взаимно независимы, но каждый из них имеет негауссову ПВ независимо и одинаково распределенных временных отсчетов, то для РСС требуется использование статистик высокого порядка [16]. Тогда процедура разделения смеси связана с методом анализа независимых компонент (АНК), в котором многомерные данные разлагаются в линейную сумму неортогональных базисных векторов, где коэффициенты статистически независимы. Напротив, когда первичные сигналы нестационарные, то для РСС достаточно статистик второго порядка.

Основные исходные положения модели

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

1) Смешивающая матрица Н имеет полный ранг по столбцам.

2) ПС {5,. (к)}являются статистически независимыми и имеют нулевое среднее.

3) При условии, что г,] = 1,...,п и г Ф ], отношение П2 (к) / П5 (к) при изменении времени (аргумента к) не является постоянной величиной.13

Поскольку гамма-функция Г(х) имеет много сингулярных точек (обращаясь в да), особенно при х < 0, это важно учитывать при выборе тг < 1.

Оценка дисперсии сигнала П2 (к) осуществляется методом скользящего среднего:

П (к) = (1 2 (к -1) + щ У](к).

Первые два допущения обычны для большинства подходов к РСС, третье — существенно для анализируемых ниже алгоритмов, т. к. именно оно позволяет осуществлять РСС с использованием только статистик второго порядка (СВТ).

Критерий на основе СВТ для формирования алгоритмов РСС

Критерий для РСС по методу АНК (при приведенных выше допущениях) основан на взаимной информации. Формирование этой характеристики требует знания ПВ восстановленных сигналов, но поскольку эти ПВ в начале процедуры РСС неизвестны, то алгоритмы, связанные с методом АНК, основываются на гипотетических распределениях. При этом явно или неявно в схему получения алгоритмов РСС приходится включать и статистики более высокого порядка.

Достаточно эффективным является критерий Матсуоки [14], который позволяет не использовать кросс-корреляционную матрицу E{ yi (k) yj (k)} (вычисления которой достаточно трудоемки):

J(W, k) = (1/ 2)]Г log (с?2 (k)) - log det R£, (28)

i=1

где с2y(k) = (y2(k)) — оценка переменной во времени дисперсии в реальном времени (РВ); det(-) — обозначает определитель матрицы;

Ryy1 = (^(k )УT (k)) — автокорреляционная матрица

сигналов y(k), которая оценивается тоже методом скользящего среднего в РВ:

Ryy=(i-^o)iRyy-1) +% y(k )y T(k).

Критерий (28) представляет неотрицательную функцию, достигающую минимума, только если выполняется условие E{ yi (k) yj (k)} = 0

(i, j = 1,...,n и i Ф j), т. е. если корреляционная матрица сигналов является диагональной.14

Алгоритм обучения с использованием м-градиента

Алгоритмы обучения ориентированы на два типа структуры нейронной сети, обеспечивающей РСС (рис. 3): 1) нейронная сеть (НС) с прямым распространением сигналов, 2) полносвязная рекуррентная нейронная сеть.

1) Выходной сигнал у (к) нейронной сети с прямым распространением (рис. 3, а) определяется выражением у(к) = W х(к). Поэтому полный дифференциал МС^ критерия оптимизации параметров НС (матрицы W ее синаптических весов) определяется соотношением:

М (W) = J ^ + dW) - J (W) =

= \d {! >в< у2(к))}- (29)

- ^d{logde^ у(к )у т(к))}.

Для второго слагаемого в (29) справедлива оценка:

d{logdet( у(к)ут(к))} = d{logdet(WRxxWT)} =

= 2d{log det(W-1)} + d{log det и} =

= 2tr{dWW-1} + d{log det 11„}. При этом можно опустить 11 хх, т. к. у него нет зависимости от матрицы весов W. Поэтому, принимая матричный дифференциал dX (при т = п) в виде dX = dWW-1, получаем оценку

d{log det( у (к )у т(к))} = 2tr{dX}.

Подобным образом оценивается и первое слагаемое в (29):

d {! >g< y2(k)> }=!

2( y, (k )dy, (k))

Это следует из так называемого неравенства Адамара

det(R) <ПП=1 Г , где Я = [г ] — неотрицательно опре-

ч

деленная (п х п)-матрица , и равенство достигается

только при условии, что для / Ф . г. = 0. Поэтому

ч

определитель принимает максимальное (а (-det) — минимальное) значение как раз при условии, что корреляционная матрица сигналов 1уу — диагональная, т. е.

Е{у(к)у.(к)} = 0 (/,. = 1,...,п и /Ф.).

(У2(к))

= 2Е{ут (к )Л-1 (к )ёу(к)} = 2Е{у т (к) Л-1 (к ^Ху(к)}, где Л(к) — диагональная матрица, у которой элемент / на диагонали Xi = ау = (у2(к. Использование матричного дифференциала dX позволяет представить критерий обучения НС в виде соотношения М(W)/d = Е{Л-1 (к)у (к)ут(к)} -1. С учетом того, что ЛW(к) = ДХ(к)W(k), применение градиентного метода непосредственно дает алгоритм обучения в РВ: ДW(k) = -^(к) х

х[М X] W(k) = ^(к) [I - Л-1(к )у (к )ут(к )]W(k), который можно также представить в виде

W(k +1) = W(k) + ц{к)[ Л(к) - 11уу) ] W(k), (30)

а

б

х1(к)

Х2(к)

хга(к)

\Упт

Рис. 3. Две структуры нейронных сетей.

а — с прямым распространением сигнала; б — полносвязная рекуррентная нейронная сеть

где Л — диагональная матрица с элементом г на диагонали « Е{у,2}, который можно оценивать в РВ: (к) = (1 -ц0)Л{(к -1) + ^0у2(к) . Так же оценивается корреляционная матрица выходных сигналов НС: К<? = (1 - ^ )Яуу-1) + ^ у (к )ут (к).

Алгоритм (30) отражает одну из форм метода АНК и приспособлен к реализации РСС при нестационарных ПС [7], он устойчив и практически не зависит от вида распределения ПС [17].

2) Структурой полносвязной рекуррентной нейронной сети (рис. 3, б), предназначенной для реализации РСС, выполняются преобразования, в соответствии с которыми выходной вектор-сигнал у(к) определяется соотношением у (к) = х(к) -

-\\у(к) = (I + \\)-1 х(к) (при естественном предположении, что матрица (I + \) не сингулярная). Путем анализа, аналогичного проведенному для НС с прямым распространением сигнала, получается алгоритм обучения для полносвязной рекуррентной нейронной сети. Матрица \ весовых коэффициентов НС, которая одновременно дает пре-

образование, обратное к смешивающей матрице, удовлетворяет соотношению

д\(к) = -7 (к) [I + \ (к)] [ Л(к) - ]. (31)

ЗАКЛЮЧЕНИЕ

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

- Адаптивные алгоритмы обучения НС по-

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

— Даны критерий взаимной информации и критерий Кульбака—Лейблера, определяющие меру различия распределений.

— Рассмотрены схемы разделения сигналов при различных неизвестных распределениях первичных сигналов (ПС) и количестве ПС.

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

— Активационные функции для нейронных структур выведены вполне корректно для распределения ПС в виде ЭС_ПВ.

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

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

Построению алгоритмов с использованием м-градиента посвящены работы Амари и Кардосо. Основные подходы к РСС нестационарных сигналов исследованы Матсуокой, Охиа и Кавамотой.

СПИСОК ЛИТЕРАТУРЫ

1. Меркушева А.В., Малыхина Г. Ф. Методы и алгоритмы разделения смеси сигналов. I. Применение декорреляции и статистик второго порядка // Научное приборостроение. 2009. Т. 19, № 2. С. 90-103.

2. Кульбак С. Теория информации и статистика. М.: Изд. Наука, 1967. 408 с.

3. Bell A.J., Sejnowski T.J. An Information Maximization Approach to Signal Mixture Separation and Deconvolution // Neural Computation. 1995. V. 7, N 6. P. 1129-1139.

4. Comon P. Independent Component Analysis, a New Conception? // Signal Processing. 1994. N 3. April. P. 287-314.

5. Moulines E., Duhamel P., Cardoso J.-F., Mayrar-gue S. Subspace Method for the Identification of Multichannel FIR Filters // IEEE Transaction on Signal Processing. 1995. V. 43, N 2. P. 516-525.

6. Amari S., Cardoso J.-F. Source Separation -Semi-Parametric Statistical Approach // IEEE Transaction on Signal Processing. 1997. V. 45, N 11. P. 2692-2700.

7. Amari S., Chen T.-P., Cichocki A. Stability Analysis of Adaptive Blind Source Separation // Neural Networks. 1997. V.10, N 8. P. 1345-1351.

8. Cichocki A., Unbehauen R. Robustneural Networks on-Line Learning for Blind Identification and Separation of Sources // IEEE Transaction on Circuits and Systems I: Fundamental Theory and Applications. 1996. V. 43, N 11. P. 894-906.

9. Cardoso J.-F., LaheldB.H. Equivariant Adaptive Source Separation // IEEE Transaction on Signal Processing. 1996. V. 14, N 12. P. 3017-3030.

10. Atick J.J., Redlich A.N. Convergent Algorithm for Sensory Receptive Field Development // Neural Computation. 1993. V. 5, N 1. P. 45-60.

11. Amari S. Natural Gradient Works Efficiently in Learning // Neural Computation. 1998. V. 10. P.271-276.

12. Малыхина Г. Ф., Меркушева А.В. Детектирование речевого сигнала и фильтрация с адаптивным порогом // Микропроцессорные средства измерений. Сб. трудов ФТК, № 2. Изд-во СПб ГПТУ, 2001. С. 26-35.

13. Меркушева А.В. Фильтрация нестационарного сигнала (речи) в вейвлет-области с адаптацией к виду и динамике шума // Научное приборостроение. 2003. Т. 13, № 2. С. 73-87.

14. Matsuoka K., Ohia M., Kawamoto M. A Neural Net for Separation of non-Stationary Signals // Neural Networks. 1995. V. 8, N 2. P. 411-419.

15. Choi S., Cichocki A., Amari S. Flexible Independent Component Analysis // Journal of VLSI Signal Processing. 2000. T. 26, N 1 / 2. P 25-38.

16. Малыхина Г.Ф., Меркушева А.В. Сеть с симметричной функцией преобразования нейронов на основе статистик высокого порядка (СВП) для снижения шума нестационарного сигнала // Научное приборостроение. 2008. Т. 18, № 4. С. 129-136.

17. Choi S., Cichocki A., Amari S. Equivariant Un-stationary Source Separation // Neural Networks. 2002.V. 15, N 1.

Санкт-Петербург

Материал поступил в редакцию 23.06.2009.

SIGNAL MIXTURE DECOMPOSITION METHODS AND ALGORITHMS.

II. APPLICATION OF M-GRADIENT TO ANALYSIS OF INDEPENDENT

COMPONENTS

V. M. Malykhin, A. V. Merkusheva

Saint-Petersburg

The problem of signal mixture decomposition (with reconstruction of its component forms) is analyzed without any information concerning the proportions and type of mixing. The method is based on using information criteria and adaptive algorithm for learning neural network. Several forms of distributions were considered for original signals reaching the sensors of information-measurement system. Variations of neural transfer functions corresponding to several types of distributions are given. The approach to signal mixture separation includes application of m-gradient in the analysis scheme for independent component analysis.

Keywords: signals mixture, separation methods, adaptive algorithms, neural network, m-gradient, independent component analysis