Методы оценки энергетического спектра случайных процессов и их сравнительные характеристики Текст научной статьи по специальности «Математика»

2005

НАУЧНЫЙ ВЕСТНИК МГТУ ГА 93

сер. Радиофизика и радиотехника

№

УДК 629.735

Методы оценки энергетического спектра случайных процессов и их сравнительные характеристики

А.А. БУТЕНКО

Статья представлена доктором технических наук, профессором Шахтариным Б.И.

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

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

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

Е{8(та)} = 8х(та). (1)

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

1 к-1

8Б (®) = - У

У х(п + Ш)в

п=О

, где К=КЬ , (2)

т.е. х(п) - исходная последовательность, разделенная на К неперекрывающихся последовательностей, длинной Ь.

Рассмотрим следующий пример. Пусть у нас есть корреляционная функция вида:

Рп = + 2.0сов[2тс(0,2)п] + гп, п = 0,1,...10,

где 5п - символ Кронекера; гп - нормально распределенная величина с нулевым средним значением и стандартным отклонением 0,1. Оценим спектр по методу Бартлетта. Все построения проведем в пакете МАТЬАВ. Наложив 3 реализации, получим следующий график (рис. 1).

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

Представим известные значения АКП в виде:

2

Рис.1. Оценка спектра по методу максимума Бартлетта. 1-я кривая: оценка по методу Бартлетта для случая Гп =0; 2-я и 3-я кривые: оценка по методу Бартлетта для случая, когда Гп — нормально распределенная величина с нулевым средним значением и стандартным отклонением 0,1

г,(к) = — /х,(шУ , (|к| < N), (3)

-п

где N - количество отсчетов.

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

р

х(п) = акх(п -к) + Ь0и(п) , (4)

к=1

где и (п) и х(п) - соответственно входное воздействие и реакция фильтра на это воздействие.

Если обозначить 2-преобразования входного и выходного сигналов как и (2) и X (2), то передаточная функция этого фильтра Н (2) описывается выражением:

Н(2) = ^ =-----------Ь-------= -Ь0-. (5)

и(2) 1 ^ -1 А(2) V 7

1 + Е ак2

к=1

Пусть и (п) - стационарный в широком смысле дискретный случайный процесс с нулевым средним значением. Выход х(п) фильтра также будет процессом, стационарным в широком смысле. Обозначим 2-преобразования корреляционных функций ги (к) и гх (к) как 8и (2) = 2 {ги (к)} , Хх (2) = 2\гх(к)} . Известно, что связь между 2-преобразованиями 5и (2) и (2) имеет следующий вид: (2) = 8и (2)Н(2)Н*(1/ 2*).

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

_.2 7 2

а А

А( 2) А" (1/2* )

Без потерь общности можно положить Ь0 = 1, так как входное воздействие можно соответ-

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

_2

Я (I) =-------------. (6)

х А(і)А (1/)

Оценка спектра методом максимальной энтропии для той же последовательности: Рп = + 2.0ео8[2п(0,2)и] + гп, п = 0,1,...10.

S (т ) , дБ

У / П

Рис. 2. Оценка спектра по методу максимума энтропии. 1-я кривая: оценка по методу максимума энтропии для случая Гп =0; 2-я и 3-я кривые: оценка по методу максимума энтропии для случая, когда Гп — нормально распределенная величина с нулевым средним значением и стандартным отклонением 0,1

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

L (*>•••> xn; 0 ) = p (xi; 0 )••• p (xn; 0), (7)

где p(t;q) в случае непрерывного распределения интерпретируется как плотность вероятности случайной величины X, а в дискретном случае - как вероятность того, что случайная величина Х примет значение t. Функцию L(Xi, . . ., Xn; q) от случайных величин Xi, . . ., Xn называют функцией правдоподобия, а оценкой максимального правдоподобия параметра q называют такое значение оценки 0 = 0 (Xi, . . ., Xn) (являющееся случайной величиной), при котором функция правдоподобия достигает наибольшего возможного значения. Так как точка максимума для log L та же, что и для L, то для нахождения оценок максимального правдоподобия следует решить так называемое уравнение правдоподобия

dLogL( Xi, X2,... Xn; 0) = Q d 0 .

Метод максимального правдоподобия не всегда приводит к приемлемым результатам, однако можно утверждать, что, если для параметра q существует несмещенная эффективная оценка q* по выборке объёма N, то уравнение правдоподобия имеет единственное решение ® ш ^ . Проделав аналогичные операции и наложив 3 реализации на один график, получаем оценку энергетического спектра по методу максимального правдоподобия для автокорреляционной последовательности.

S (т ) , дБ

У / П

Рис. 3. Оценка спектра по методу максимального правдоподобия. 1-я кривая: оценка по методу максимального правдоподобия для случая rn =0, 2-я и 3-я кривые: оценка по методу максимального правдоподобия для случая,

когда Гп - нормально распределенная величина с нулевым средним значением и стандартным отклонением 0,1

Рассмотрим еще один пример. Пусть автокорреляционная функция задана в виде: p(t) = [e^(cosP1t + a1/P1sinP1t)] + [[P2e~“21(cosP2t + a2/P2sinP2t)], где t=0,1,2... 10. Для P2 = 2, P1 = 2п(0,15) и P2 = 2п(0,3) получим оценки энергетического спектра при:

a) a1 = 0,02, a 2 = 0,04;

b) a1 = 0,04, a 2 = 0,04;

c) a1 = 0,12, a2 = 0,0б .

Предварительно получим теоретическое значение энергетического спектра, используя теорему Винера-Хинчина. Графики изображены на рис. 4-б соответственно.

S (т), дБ

Рис.4. Различные оценки энергетического спектра для а1 = 0,02; а2 = 0,04. 1-я кривая - теоретическое значение; 2-я кривая - оценка по максимуму энтропии; 3-я кривая - оценка спектра по методу максимума правдоподобия; 4-я кривая - оценка спектра по методу Бартлетта

Ш / П

Аналогичные отображения кривых введены на рис. 5 и рис. 6.

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

Возьмем два процесса АВ.(2) и АВ.(4), полученные фильтрацией белого гауссовского шума с единичной дисперсией фильтрами, имеющими следующие вектора коэффициентов:

а2 = [1 - 0,75 0,5]г, й0 = 1, (9)

а4 = [1 -2,7607 3,8106 -2,6535 0,9238]г, й0 = 1. (10)

£ (то), дБ

і / \ і

Л і /уГ V і

уф \\7Т\ : і.

к\ \ і і д/ д. \ : \ і

\ і / \ \і

■'Vі \ 2 \ : 1 /\ і

^ 1

і / Ні У і V ..і..

: : 1 л

: І : і — \1

Рис.5. Различные оценки энергетического спектра для а1 = 0,04; а2 = 0,04. 1-я кривая - теоретическое значение спектра; 2-я кривая - оценка спектра по методу максимума энтропии; 3-я кривая - оценка спектра по методу максимума правдоподобия; 4-я кривая - оценка по методу Бартлетта

Ш / П

8(хп),дБ

10 5 0

-10 -15 -20 -25 -ЗО

3

/

0 0.05 0.1 0.15 0.2 0.25 0.3 0.35 0.4 0.45 0 5

Рис. 6. Различные оценки энергетического спектра для а1 = 0,12; а2 = 0,06 . 1-я кривая - теоретическое значение спектра; 2-я кривая - оценка по максимуму энтропии; 3-я кривая - оценка спектра по методу максимума правдоподобия; 4-я кривая - оценка спектра по методу Бартлетта

Ш / П

Будем генерировать N = 1024 отсчета каждого процесса. Используя периодограммный метод и метод модифицированной периодограммы с окном Хэмминга, получим оценки энергетического спектра этих процессов и сравним полученные результаты.

Все построения проведем в пакете Matlab. На рис. 7, 8 показаны реализации процессов AR(2) и AR(4). В целях наглядности мы выведем графики только первых 200 отсчетов.

Для AR(2) процесса программа в Matlab примет вид:

>> N=1024;

>> X0=randn(1,N);

>> a=[1 -0.75 0.5];

>> X=f^lter(1,[a], X0);

>> plot(X(1:200));

ГП

ІЮ

Рис. 7. Реализация процесса ЛЯ(2)

Аналогично для AR(4) процесса:

>> N=1024;

>> X0=randn(1,N);

>> a=[1 -2.7607 3.8106 -2.6535 0.9238]; >> X=filter(1,[a], X0); >> plot(X(1:200));

Рис.8. Реализация процесса AR(4)

Построим теоретические энергетические спектры данных процессов, чтобы в дальнейшем было с чем сравнивать получаемые результаты. Поскольку мы имеем дело с авторегрессионы-ми процессами второго и четвертого порядков, спектральная плотность мощности будет пропорциональна квадрату модуля коэффициента функции передачи фильтра. Воспользуемся формулой (11), положив в ней сначала N = 2, а потом N = 4. Кроме того, дисперсию белого шума а2п и частоту дискретизации fd считаем равными единице. Функции спектрального анализа Signal Processing пакета Matlab для вещественного сигнала рассчитывают односторонний спектр плотности мощности, поэтому умножим результат вычислений по формуле (11) на два.

1

^ |ч - jюT - j2юГ -iNюТ I

Jд 1 - ае - @2^ -... - @е

Строим график спектра плотности мощности для AR(2) процесса (рис. 9): >> Т=1;

>> уаг=1;

>> £=Н^расе(0,1/(2*Т),100);

>> P=2*var*T./abs(1+0.75*exp(-i*2*pi*f*T)-0.5*exp(-i*2*2*pi*f*T)).Л2;

>> psdplot(P,f,'Hz').

5 (ш ) , дБ

(11)

Рис.9. Теоретический спектр процесса AR(2)

Для AR(4) соответственно (рис. 10):

>> T=1;

>> var=1;

>> f=linspace(0,1/(2*T),100);

A

ш / П

>>2*уаг*Т./аЬ8(1+2.7607*ехр(-і*2*рі*РТ)-3.8106*ехр(-і*2*2*рі*РТ)+2.6535* ехр(-і*2*3*рі*РТ)-0.9238*ехр(-і*2*4*рі*РТ)).л2; >> psdplot(P,f,'Hz').

£ (то), дБ

50

40 ---

30 ---

20 ---

1 О ---

О ---

-1 О ---

-20

-30 С

0.05 0.1 О. 15 0.2 0.25 0.3 0.35 О. А 0.45 О 5

Рис.10. Теоретический спектр процесса ЛЯ(4)

Теперь получим оценки энергетического спектра периодограммным методом. Для этого воспользуемся следующими процедурами:

Для AR(2):

>> N=1024;

>> X0=randn(1,N);

>> a=[1 -0.75 0.5];

>> X=filter(1,[a], Х0);

>> periodogram(X,[],[],1);

Для АЕ.(4): — / п

>> N=1024;

>> X0=randn(1,N);

>> a=[1 -2.7607 3.8106 -2.6535 0.9238];

>> X=fllt.er(1,[a], Х0);

>> periodogram(X,[],[],1);

Рис. 11 и 12 демонстрируют полученные результаты.

Ш / П

Рис.11. Наложение теоретического спектра процесса ЛЯ(2) и его периодограммной оценки

Для получения модифицированной периодограммы воспользуемся функцией МаЙаЬ w=hamming(n,,sflag,), реализующей окно Хемминга. Мы остановимся на симметричном случае, тогда параметр можно опустить, а процедуру periodogram(X,[],[],1) заменим на

periodogram(X,w,[],1). Результаты отражены на рис. 13 и рис.14.

8{тп),дБ

Ш / П

Рис.14. Наложение теоретического спектра процесса ЛЯ(4) и его модифицированной периодограммной оценки

ВЫВОДЫ

Рассмотренные примеры позволяют сделать следующие выводы:

1. Методы максимума энтропии и метод максимального правдоподобия позволяют существенно улучшить разрешающую способность получаемых оценок энергетического спектра (в отличие от метода Бартлетта).

2. Метод максимума энтропии позволяет получить более качественную оценку (с меньшей дисперсией) энергетического спектра в отличие от метода максимального правдоподобия.

3. Периодограмма является некорректной оценкой спектра для AR(4) процесса. Наклон графика обнаружен не был.

4. В случае AR(2) имеет место меньший динамический диапазон спектра (около 15 дБ).

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

6. В обоих случаях сильные случайные флуктуации введением окна уменьшить не удалось.

ЛИТЕРАТУРА

1.Шахтарин Б.И. Случайные процессы в радиотехнике. М.: Радио и связь, 2000.

2.Шахтарин Б.И. Случайные процессы в радиотехнике. 2-е изд. Ч.1. Линейные системы. М.: Радио и связь, 2002.

3.Сергиенко А.Б. Цифровая обработка сигналов. СПб.: Питер, 2003.

Л.А. Butenko

Estimation methods for power spectrum of random processes and their comparison

The article is devoted to the research of perspective methods of estimation of a power spectrum of random processes. It is a method of maximum likelihood and maximum entropy and periodogram estimations. On the concrete examples it is shown, that the given methods can improve the quality of received estimations. Comparison of these methods is carried out.

Сведения об авторе

Бутенко Алексей Анатольевич, 1978 г.р., окончил МГТУ им. Баумана (2001), аспирант кафедры автономных информационных и управляющих систем МГТУ им. Н.Э. Баумана, автор 7 научных работ, область научных интересов - нейрокомпьютерная техника, прогнозирование и фильтрация случайных процессов.