Алгоритм оптимальной обработки широкополосньгх сигналов в пространстве состояний Текст научной статьи по специальности «Математика»

Алгоритм оптимальной обработки широкополосных сигналов в пространстве состояний

В.М. Советов, д.т.н., профессор Российского государственного университета туризма и сервиса

(РГУТиС), г. Москва, e-mail: [email protected]

В.А. Коекин, аспирант РГУТиСа, г. Москва, e-mail: [email protected]

Проведен синтез алгоритма оптимального приема сигнала; показано, что алгоритм Калмана для рассматриваемой модели может быть представлен в виде рециркулятора с преобразованием.

Synthesis of optimal algorithm of receiving of signals is performed. It is shown Kalman algorithm for given model can be represented as recirculator with transformation.

Ключевые слова: пространство состояний, отношение правдоподобия, фильтр Калмана.

Key words: state space model, likelihood ratio, Kalman filter.

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

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

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

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

Как известно, для представления сигнала можно использовать модель динамической системы в пространстве состояния. В частности, дискретный сигнал можно представить в виде уравнения перехода состояний (1) и уравнения выхода (2):

х (к + 1) = Ах ( к ), (1)

у (к) = Сх(к), (2)

где А - матрица перехода состояний (МПС),

при умножении на которую вектора состояния х(к)

на к-м шаге получается вектор состояния на (к+1)-м шаге; С - матрица выхода, определяющая связи переменных состояний с выходом системы.

Например, модель динамической системы в пространстве состояний, образующая гармоническое колебание будет иметь вид

X,

1 (к +!)

0

1

Х2 ( к + l)J L-1 2cos«r _||^ х2 (к )

X

1(к)

У(к ) = [1 °]

X1(к)' х2 (к )_

При этом МПС A =

0

1

-1 2 cos oT

находится из

сигнала

Pr(Y(к) | xод(к)) =

Pr[xо і (к), Y (к)]

‘и (0 ) =

Z -преобразования синусоидального (здесь T - интервал дискретизации).

Начальная фаза и амплитуда колебания будут зависеть от вектора начального состояния. Для образования дискретного гармонического колебания с начальной фазой и в качестве начального вектора состояний необходимо задать вектор 8ш(и) sin (юТ + и)

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

Рассмотрим задачу оптимального приема такого сигнала в общем виде. Пусть начальный вектор x(0) представляет случайный вектор сообщения с детерминированными, но неизвестными значениями элементов. При этом для простоты рассмотрим случай проверки двух гипотез:

Н : х0(£ +1) = Ax0(k); y0(k) = Cx0(k) + v(k), к = 1,2,..., K ,

H: xi(k +1) = АхДк); yt(k) = СхДк) + v(k), к = 1, 2,. , K .

Аддитивная помеха v(k) представляет белый шум с нормальным распределением и известными статистическими характеристиками

E{ v (к)} = E { v 0( к)} = E { v 1( к)} = 0; cov{ v(к), v(j)} = Vv (к)Sk (к - j) . Решение о выборе той или иной гипотезы принимается по K наблюдениям:

Y(к) = y(1), y(2),. ,y(к).

Если векторный параметр х(к) не случаен и его значение неизвестно, то можно воспользоваться обобщенным отношением правдоподобия [1] maxPr(Y(к) | x(к), Н1)

Л (Y (к)) = -^(к)

(3)

maxPr( Y (к) | x (к), H 0)

Хо( к)

Рг( х од( к))

Представим Y(k) как объединение нового наблюдения у(к) и У (к— 1) предыдущих наблюдений:

Рг[ х 01( к), у (к), У (к -1)] (4)

рГ( у (к )|х оХ( к)) = „ —-• (4)

Рг( х од( к))

Числитель выражения (4) можно представить как Рг[х0,1 (к),у(к),У (к - 1)] = Рг[у(к) | хоД(к)] х

X Рг[ход(к) | У (к - 1)]Рг[У (к - 1)].

Отсюда Р г( У (к )|х од( к)) =

= Рг[ у (к) | ход(к )]Рг[ х од(к) | У (к - 1)]Рг[У(к -1)]

Рг(хоД(к)) •

При равной вероятности гипотез обобщенное отношение правдоподобия (3) будет иметь вид т ахРг[ у (к )| хх(к)]

Л (У (к )) = -^(к) х

maxPr[y (к )| xо(к)]

Хо( к)

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

Вычислим правдоподобие для какой-либо гипотезы Но или Н при получении к наблюдений. По теореме умножения

.. Pr[x 1(к) | Y (к - 1)] , (5)

X * ' '

Pr[ x °(к) | Y (к - 1)]

Найдем плотность вероятности Рг[у(к)| х0д(к)]. В зависимости от гипотезы наблюдение у(к) связано с состоянием системы х(к) линейным уравнением

У °,1(к) = Cx °,1(к) + v(к). (6)

Ясно, что плотность вероятности Рг[у(к)| Х0,х(к)] будет нормальной, поэтому для записи этой плотности вероятности необходимо найти лишь математическое ожидание и дисперсию. Среднее значение процесса

Е{У(к) | x°,1(к)} = E{[Cx°,1(к) +

+v(к)]| x0,1(к)} = Cx°,1(к).

Дисперсия случайного процесса по определению равна

уаг[у(к) | х°,1(к)] = Е{[Сх°д(к) + v(k)] х

x[Cx °д( к) + v (к )]т} = Cx °д(к) x т д( к )Ст + Vv (к).

Отсюда плотность вероятности Рг[у(к)| xo,!^)] будет иметь вид

Pr[у(к)| x°,1(к)] =

1

=---------------------------------------------V2 х

(2л)”/2 {det[Cx°д(к)xтд(к)Ст + Vv (к)]}V

х ехр {-\[ у (к) - Cx °,1 (к )]т x[Cx °д(к )x °д(к )Ст + + Vv (к )]-1[у (к) - C^)]}.

Найдем плотность вероятности Рг[х0д(к)| У(к)]. Она тоже будет нормальной. Среднее значение

вычисляется как

Е[хо1(к) | У (к - 1)] = Ахо,1 (к - 1) = хо,1(к | к - 1).

То есть среднее значение вектора состояния х0; х(к) при получении (к-1) наблюдений представляет собой предсказание оценки ход(к | к - 1) •

Дисперсия х(к) при полученной последовательности У(к — 1) совпадает с дисперсией ошибки оценки Vо,(к | к - 1):

уаг[х(к) | У (к -1)] = Уход (к | к -1) = АУход (к -1) А

Таким образом, плотность вероятности Рг[ход(к)|У(к-1)] будет иметь вид

Рг[ход(£)1 У(к - !)] = ■

1

(2п)п> 2[ае1 V* (к|к -1)]

1/2

х ехр ^ - 2[х(к) - хо,1(к 1 к - 1)]т ^ (к 1 к -1) х к -1)

х[х(к) - хо,1(к | к -1)]}.

Используя полученные выражения, можно записать обобщенное отношение правдоподобия (5) в виде Л(У(к)) =

аег[Сх0 (к)*о (к)ст + vv (к)] аег vXо (к | к -1)

ае1[Схх (к)*т (к)Ст + % (к)] ёе1 VX1 (к | к -1) ехр^- ^2[у(к) - Сх1(к)]т [схо,1(к)хіт (к)с

112

тах

х1(к)

+

тах

хо(к)

+

+

«хр-^- >) - &о(к )Г [сх«(к)х0 (к )с’ V,(к)]-‘[у(к) - Сх,(к)]}

V, (к) ]-‘[ у (к) - Сх о( к )]}Х 1

+

ехр і-у[ х (к) - хі( к | к - 1)]Т

х-----і-------------------------х

ехР і - ^[х(к) - хо (к | к - 1)]т

Vх11 (к | к - 1)[х(к) - х1(к | к - 1)]}

V5-С)1 (к | к - 1)[х(к) - хо(к | к - 1)]}

Используя определения шума измерения у(к)=[у(к)-Схо1(к)] и ошибки оценки состоя-

хо1 (к | к -1) = [хо1 (к) - хо,1 (к | к -1)], упростим запись и представим обобщенное отношение правдоподобия как

тах

Л (У (к)) =

Ао х1( к)

тах

хо( к)

ехр ^- ^ ¥1(к )Т V“(к) “1(к)

V о (к)т V,-; (к), о( к)

ехр|- ;2 х1(к | к - 1)т ^1(к | к - 1)х1(к | к -1)'|

ехр|-1 хо(к | к - 1)т ^,1(к | к - 1)хо(к | к - О?

(7)

т ГДЄ ~~ =

Ао {аеі[Схо(к)хо(к)Ст +

А ~ {аег[Сх1(к)хт (к)Ст +

+% (к )]аег Vxо(k |к -1)}

- общий множитель.

+% (к)]аеі Vxl(k|k -1)}2

Взяв логарифм с каждой стороны (7), получаем: 1п Л (У (к)) = 1п Ао - 1п А1 +

+тах

х1( к)

1

-1 “1(к) Т VV11 (к) “1(к) -

- - х1(к|к - 1)т Vx11 (к|к - 1)х#|к -1)

- тах

хо(к)

- 22 хо(к | к -1) т vx"01 (к | к - 1)хо(к | к -1)

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

- -2 “1(к) т Vvl1 (к) “1(к) -

- - хД к|к - 1)т к | к - 1) хД к|к - 1)

= о,

- 2 х о (к | к -1)т Vxо1 (к|к - 1)х о (к | к -1) = о.

При дифференцировании используем правило [2]: V0(атЬ) = V0(ат)Ь + V0(Ьт)а, а, Ь є ^*0 х1 . В результате получим СтVv-о1l(к)[у(к) - СхоД(к)] -

-VхД (к|к - 1)хоД(к|к - 1).

х

X

X

х

X

Приравняв полученное выражение к нулю и решив его относительно хо1(к), придем к выражению

х тах о,1( к ) = [ V,-! ( к\к - 1) + С т V-1,, ( к )С ]-‘ х

х[Vx01l (к | к - 1)х о,1 (k\k -1) + Ст V,-! (к )у(к)]. (8)

Выражение можно упростить с использованием леммы обращения [2]:

(А + ББт)-1 = А-1 - А- 1В(1 + БтА- 1В)-1 БтА-1.

Определим термины леммы как А = V-1 (к |к-1), Б = СХ"1 (к), а Бт = С

хо,1 4 1 ' “ о,1 4 '

и в результате получим:

[vs-01l(к|к -1) + Ст^(к)С]-1 =

= vSoj (k \k -1) - vSoj (k I k - 1)Cту;! (k) x

x[I+CVxo i (k I k - 1)CTV-011 (k)]-1 CVxo i(k I k -1). (9) Представим выражение в квадратных скобках как

[I + cvSo (k\k- 1)Cту;011(k)]-1 =

т-і-1

= У^(к)[УуоД(к) + СУхоД(к|к -1)С]

Выражение в квадратных скобках справа соответствует дисперсии процесса обновления

ео,1 (к) = [у (к) - Схо,1(к | к - 1)]. Поскольку

Е{у (к )|У(к -1)} =

= СЕ{хоД (к) | У (к -1)} + Е{у(к) | У (к -1)} = = СЕ{хо1(к) | У(к - 1)} = Схо,1(к | к - 1),

то дисперсия процесса обновления уаг[у (к )| У (к -1)] =

= Е{[ у (к) - Сх о,1 (к | к -1)][ у (к) - Сх о,1 (к | к - 1)]т }.

Подставим в выражение для дисперсии вектор наблюдения (6) и получим ущ-[у (к)| У (к -1)] =

= Е{[Схо1(к) + V(к) - Схо,1 (к | к -1)] х х[Схо1(к) + V(к) - Схо,1(к | к - 1)]т} = = Е{ [Схо,1(к | к -1) + v(к)][Cxо,l(к | к -1) + v(k)]т} = = Vv (к) + СУх о1( к|к - 1)Ст = уе о1( к). Таким образом,

[I + СУ^ (к | к - 1)СтУ-о11 (к)]-1 = Уч, (к)Уео1 (к),

и (9) принимает вид

[ ухо11( к|к -1) + с т у;о11( к )С ]-1 = = Ух оД( к | к -1)- Ух оД( к|к -1)С т X хУе-11(к)СУхоД(к|к -1) = = [I - Ухо,1 (к | к - 1)СтУе-оХ1 (к)С]Ухо,1 (к | к -1) = = [ I - К од( к )С ]Ух о1( к|к -1),

где Код(к) = Ухо1(к|к-1)СтУе-о11(к) - весо-

вая матрица или коэффициент усиления.

Обозначим

Уход (к | к) = [Ухо11 (к | к -1) + Сту;1 (к)С]-1 .(1о)

и запишем (8) как

хтаход(к) =

=Уход(к | к)[Ухо11(к | к-1)хо ,1(к|к-1)+СгУ-о11(к)у(к)]. Из (10) найдем У--1 (к | к -1) и подставим

х о ,1 4 1 7

в полученное выражение:

хтах о,1 (к) = Ух (к | к) X

[Vx0^1 (k \ k) - CTV;011 (k)C]xo,1(k I k -1)

+

+CVv°1 (к )у(к)

Перемножая и группируя члены этого выражения, представим его как

:хmax °,1 (к) = x °,1 (к | к - 1) +

+Vs°,1 (к|к)CтVv-°11 (к) х

х у (к) - Cx°,1 (к | к - 1) , или окончательно:

xmax °,1 (к) = :х°,1 (к | к - 1) + K°д(к) х

у(к) - Cx°,1 (к | к - 1)

Таким образом, обобщенное отношение правдоподобия (7) можно записать как ln Л (Y (к)) = ln A° - ln A1 +

+x1( к | к -1) + к 1( к) Г у (к) - cx1( к | к -1) --x°^ | к-1) + K°(к) у(к)-С:х°(к | к-1)

= ln A ° - ln A1 + x m ax ,1 ( к ) - x max,2 (к ) Решающее правило будет иметь вид

H1

>

x max,1 (к)- x max,2 (к) < ln^ ,

H °

x

где п определяется стоимостью решений и значениями Ао1; если Ао= А1 и отношение стоимостей решений равно 1, то решающее правило сводится к вычислению

Н1 > .

,1 (к ) ^ хтах.

Н о

2 (к)

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

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

Представим оптимальный фильтр Калмана для модели (1), (2) в виде

х (к) = Ах (к -1) + К (к)[ у (к) - САх (к -1)], (11)

Ух (к + 11 к) = АУх (к) А т ,

(12)

К(к) = Ух (к | к - 1)Ст[СУх(к | к - 1)Ст + Vv ]-1 ,(13) Ух (к) = [ 1-К (к) С ] Ух (к|к -1). (14)

Рассмотрим также весовую матрицу К на (к+1)-м шаге:

К(к +1) = Ух (к +1| к)Ст X

х[СУх (к + 1|к )Ст + Vv Г1. (15)

Подставим в (12) выражение для матрицы ковариации ошибки оценки (14) и получим предсказание матрицы ковариации ошибки оценки на шаг вперед:

Ух (к +1| к) = А[1 - К (к)С]Ух (к | к -1) Ат. (16)

Теперь получим предсказание весовой матрицы на (к+1)-м шаге при известной весовой матрице на к-м шаге, для чего подставим (16) в (15) и получим

К(к +1) = А[1 - К(к)С]Ух (к | к -1) АтСт х

х[СА[1 - К (к )С]Ух (к | к - 1)А тСт + У Г1. (17)

Из (13) имеем

Ух (к | к-1)=[1-К(к )С]-1 K(k)Vv (Ст)-1.(18)

Подставив (18) в (17), получим выражение весовой матрицы на (к+1)-м шаге через весовую матрицу на к-м шаге и априорные параметры модели:

К (к +1) = АК (к) у (Ст )-1 х

хАтСт[САК(к)у (Ст)“1АтСт + У ]-1. (19) Выразим из (19) весовую матрицу К(к) через К(к+1):

К (к ) = А-1 [1-К (к + 1)С ]-1 х

х К (к + 1) Vv ( С т )“1 ( А т )-1 С т Vv-1. (2о)

Рассмотрим работу алгоритма (11) - (14) оценки вектора состояния модели (1), (2) в пошаговом режиме.

При первом наблюдении, когда к = о, по формуле (11) получим оценку

х (1) = ах (о) + К (1)[ у (1) - САх (о)].

Не нарушая общности, можно принять х(о) = о , в результате чего получим оценку х(1) в виде

х (1) = К (1) у (1). (21)

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

х(2) на шаге к = 1:

х(2) = Ах(1) + К(2)[у(2) - САх(1)]. (22)

Подставив в (22) оценку (21), получаем

х(2) = [1 - К(2)С] АК(1)у(1) + К(2)у (2). (23)

Используя (20), выразим весовую матрицу К(1) через К (2):

К(1)=А-1 [I -К(2)С]-1 К(2)УЛ,(Сг)"1(Ат)“1СтУЛ,“1.

Подставив найденное выражение в (23), получим

х (2) = К(2) [ Vv (Ст)-Х(Ат)- 1CгVv-ху(1) + у (2)].

На шаге к = 2, проделав аналогичные вычисления, получим

х (3) = К (3) ( Vv (Ст)-1 (Ат)-1Ст Vv-1 )х х( Vv (Ст)-1 (Ат)-1Ст Vv-1) у (1) +

+К(3)( Vv (Ст)-х(Ат)-хСХ-1) у (2) + К(3)у(3) =

= К (3) [( Vv (Ст)-1 (Ат)-1Ст Vv-1 )2 у (1) +

+ ( Vv (С т)-1 (А т)-1 С т Vv“1) у (2) + у (3). Очевидно, что

(у (Ст)-'(Ат)-1 CтVv -1)" = Vv (Ст)- '(Ат)-" Ст^-1

при любом п и оценку на произвольном шаге можно вычислить по формуле

х (к +1) = К (к + 1)У (Ст)-1 х

хХ ( А т )гк С^ V ^ + 1). (24)

г=о

Если помеха некоррелированная, то матрица У диагональная. И если матрица выхода имеет вид С = [1 о ... о] или диагональная (сигнал передается без искажений), то выражение (24) можно записать как

к

х (к +1) = К (к +1)£ (Ат У~к у (/ +1). (25)

i=0

Интересен случай, когда матрица А ортогональна и, следовательно, АТ=А-1. Тогда выражение (25) можно записать как

к

х (к +1) = К (к + 1)^ А к~ i у (/ +1). (26)

i=о

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

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

ЛИТЕРАТУРА

1. Сейдж Э., Мелс Дж. Теория оценивания и ее приме-

нение в связи и управлении: Пер. с англ. - М.: Связь, 1976.

2. Candy, J. KModel-based signal processing. A John Wiley & Sons, inc., publication, 2°°6.

Поступила 12.09.2009 г.