в интервале (0, fmax ). Сначала получим стандартные равномерно

распределённые

случайные величины

xi R(0,1),i =

,

по

1,n +1

этим

исходным

данным

определим

величины

ξ(k ) = a

+

(b −a

)x

,i =

,ξ(k )

= f

x

.

i

i

1,n

max

n+1

i

i

i

n+1

Пусть

(k )

= (ξ1(k ),ξ2(k ),...,ξn(k ) )T

– n-мерный случайный

век-

ξ

тор, причём k – номер его реализации, зависящий от t . Точка ξ(k ) принимается в качестве реализации случайного процесса ξ(t) с

n -мерной плотностью вероятности

f (x1,t1, x2 ,t2 ,..., xn ,tn ), если

f

(k )

≥ ξ

(k )

, и отбраковывается, если

(k )

< ξ

(k )

. При бра-

ξ

ξ

n+1

f

n+1

ковке точки ξ(k ) происходит переход к новому шагу и индекс k

увеличивается на единицу. Таким образом, здесь моделируются случайные точки n +1-мерного пространства, равномерно распре-

делённые в объёме под гиперповерхностью z = f ξ(k ) .

Как и все без исключения методы моделирования случайных процессов, метод отбора обладает достоинствами и недостатками. Достоинствами являются простота реализации и программирования, а также отсутствие всяких предварительных расчётов в противоположность, например, методу условных распределений. Основной недостаток метода – работа «вхолостую» в том случае, ко-

				(k )	≥ ξ	(k )	. Это приводит к избыточ-
		ξ
	гда нарушается условие f					n+1

ному обращению к датчикам случайных чисел и увеличивает время работы программы.

4.4. Моделирование случайных процессов с заданными корреляционными свойствами

С практической точки зрения получение возможных значений случайного процесса ξ(t) в рамках заданной корреляционной тео-

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

120

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

4.4.1. Метод формирующего фильтра. Формирующим фильт-

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

заданными статистическими характеристиками. На рис. 4.2 представлена схема аналогового формирующего фильтра, когда из бе-

лого шума η(t) на входе получают случайный процесс ξ(t) с заданной ковариационной матрицей Kξ(τ) и спектральной плотностью Sξ(ω). Характеристикой формирующего фильтра является передаточная функция K(p), которая может быть найдена по ковариационной функции Kξ(τ) и (или) по спектральной плотности Sξ(ω) нужного случайного процесса ξ(t).

η(t)

ξ(t)

K(p)

Рис. 4.2. Схема формирующего фильтра

Например, так как характеристики белого шума (см. подразд. 2.4) Kξ (t)= cδ(τ), Sξ (ω)= 2cπ = S0 , то, предполагая процесс

η(t) гауссовским с M [η(t)]= 0 , получим [22] K(p)= S(p).

S0

121

В современной литературе известно несколько разных методов вычисления передаточной функции. В соответствии со спосо-

бом вычисления K(p) метод моделирования случайного процесса ξ(t) имеет то или иное название (см., например, п. 4.4.2).

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

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

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

Пусть ξ(n) – последовательность значений случайного процесса ξ(t) в точках tn = n∆t . Оператор скользящего суммирования есть линейный оператор преобразования:

ξ(n)= ∑N ck x(n −k ),	(4.4)
k=0

где x(k ) N (0,1) – отcчëты белого шума, ck – коэффициенты циф-

рового фильтра. В силу линейности преобразования (4.4) последовательность ξ(n) будет также нормальным случайным процессом.

Название «метод скользящего суммирования» отражает суть формулы (4.4), так как значения ξ(n) вычисляются как взвешенная

сумма	входных отсчётов x(k ) в сдвигающемся окне		отсчётов
(n − N ),(n − N +1),...,n . Коррелированность		случайных	величин
ξ(n),	ξ(n + k ) обеспечивается за счёт того,	что в их образовании

участвуют k общих случайных величин последовательности x(n).

122

При k = N значения ξ(n) и ξ(n + k ) становятся некоррелирован-

ными. Характер корреляционной связи зависит только от значений коэффициентов ck .

Так как ковариационная функция последовательности x(n) имеет вид

	Kx (n)= M [x(k )x(k + n)]= δn					1,n = 0,					(4.5)
						=			≠ 0,
						0,n
то значения ковариационной функции						Kξ(t)			в точках	tn = n∆t ,
вычисленные с использованием формулы (4.4), будут равны:
			N			N
	Kξ (0)= M [ξ(k)ξ(k)]= M ∑ck x(k)∑ck x(k) =
			k=0			k=0
= M{[c0 x(0)+ c1x(1)+... + cN x(N )] [c0 x(0)+ c1x(1)+... + cN x(N )]}=
= M c2	x(0)x(0)+ c	0	x(0)c x(1)+... + c	0	x(0)c		N	x(N )+ c x(1)c		0	x(0)+
0			1						1

+c12 x(1)x(1)+... + c1x(1)cN x(N )+... + cN x(N )c0 x(0)+

+cN x(N )c1x(1)+... + cN2 x(N )x(N )]} = c02 M [x(0)x(0)]+

+c0c1M [x(0)x(1)]+... +cN2 M [x(N )x(N )]= c02 1+c0c1 0 +

N

+... +cN2 1 = ∑ck2 .

k =0

Аналогично другие индексы:

Kξ(1)= c0c1 + c1c2 +... + cN −1cN ,

Kξ (2)= c0c2 + c1c3 +... + cN −3cN −1, …,

Kξ (k )= c0ck + c1ck+1 +... + cN +1cN −k+1, …,

Kξ(N −1)= c0cN −1 + c1cN , Kξ(N )= c0cN , Kξ(N +1)= 0,

где Kξ (n)= Kξ (n∆t).

Итак, между коэффициентами линейного фильтра (4.4) ck и значениями ковариационной функции Kξ(τ) в точках отсчётов tn = n∆t существуют следующие соотношения:

123

K

N

ξ (0)= ∑ck2 = c02 + c12 +... + cN2 ,

K

(1)= c c + c c +... + c c ,

k

=0

ξ

0 1

1 2

N −1 N

K

ξ

(2)= c

0

c

2

+ c c

3

+... + c

N −3

c

N −1

,

1

....................................................

(4.6)

K

ξ

(k)= c

0

c

k

+ c c

k+1

+... + c

N −k

+1

c

N

+1

,

1

.....................................................

Kξ (N −1)= c0cN −1 + c1cN ,

K

ξ

(N )= c

0

c

N

.

Таким образом, очевидным способом синтеза дискретного цифрового фильтра, моделирующего случайный процесс ξ(t) с

заданной ковариационной матрицей Kξ(τ) и математическим

ожиданием M [ξ(t)]= 0 , является получение коэффициентов ck из

решения нелинейной системы алгебраических уравнений (4.6). После того как коэффициенты ck вычислены, реализации

случайного процесса в моменты tn = n∆t осуществляются по

формуле (4.4).

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

ность, которая зависит от числа учитываемых коэффициентов ck в

формуле (4.4). Его целесообразно применять при известной ковариационной функции Kξ(τ) нормального случайного процесса ξ(t).

4.4.3. Метод авторегрессии-скользящего среднего. Для моде-

лирования нормальных случайных процессов можно использовать и рекуррентные алгоритмы. Этот метод основан на уравнении ав- торегрессии-скользящего среднего:

ξ(n)= ∑l ak x(n −k )+ ∑m bk ξ(n −k ).		(4.7)
k=0	k=1

В данном случае вид ковариационной функции случайного процесса определяется набором значений параметров ak и bk ,

124

а также их количеством. Параметры ak и bk определяются на

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

Коэффициенты ak и bk можно определить методом факто-

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

Метод факторизации применяют и при определении коэффициентов ck линейного формирующего фильтра, когда уравнение (4.7) имеет более простую структуру, например:

ξ(n)= ∑l ck x(n −k ), xk N (0,1).	(4.8)
k=0

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

	Sξ(ω)= ∞∫K(τ)e−iωτdτ =	S1	(ω)	,	(4.9)
		S2	(ω)
	−∞
	где S1(ω) и S2 (ω) – многочлены степеней l			и m, (l < m)	соответ-

ственно.

Случайные процессы со спектральной плотностью вида (4.9) имеют также дробно-рациональную передаточную функцию [3]

	K(iω)=	K1	(iω)		,	(4.10)
		K2	(iω)

где K1(iω) и K2 (iω) – многочлены степени l и m(l < m). Спек-

тральная плотность случайного процесса ξ(t),							получаемого по
формуле (4.7) на выходе системы, будет иметь вид
Sξ(ω)=		K(iω)		2 = K(iω)K(− iω)=	K1	(iω)K1	(− iω)	.	(4.11)
					K2	(iω)K2	(− iω)

125

Множитель K1((iω)) в формуле (4.11) и будет передаточной функ-

K2 iω

цией формирующего фильтра (4.9). Разложение (4.11) возможно, хотя и не является однозначным.

По теории [19], всякая положительная дробно-рациональная функция относительно аргумента может быть представлена свои-

ми нулями и полюсами:

S1	(ω)			∏l (iω− iω0k )			2	,	(4.12)
			= Ck	k =1
S2 (ω)				m
				∏(iω− iωk )
				k =1

где ω0k – нули передаточной функции (нули многочлена K1(iω)), ωk – полюсы передаточной функции (нули многочлена K2 (iω)), Ck – константа. В соответствии с этим разложением сама передаточная функция K(p) может быть представлена формулой

K(iω)=			∏l1 (iω− iω1k )
	C		k =1		,	(4.13)
			m1
			∏(iω− iω2k )

k =1

где ω1k и ω2k – нули и полюсы передаточной функции, лежащие в верхней полуплоскости (имеющие положительную мнимую

часть), а множитель С выбирается из условия K(iω)2 = Sξ(ω).

Для получения формулы скользящего суммирования по теории необходимо получить функцию импульсной переходной характеристики формирующего фильтра [1, 3, 13]:

k ri −1	t	j
h(t)= ∑ ∑Cij			e pit ,	(4.14)
		j!
i=1 j=0

где pi – полюсы передаточной функции, т.е. корни знаменателя формулы (4.13) кратности ri каждый (r1 + r2 +... + rk = m), а

	1			d ri − j−1	k
	Cij =				[K(p)(p − pi ) i ]	p=pi .	(4.15)
		(ri − j −1)!		dpri − j−1

126

Если на вход фильтра с импульсной переходной характеристикой h(t) воздействует белый шум с ковариационной функцией

Kx (τ)= δ(τ), то на выходе фильтра случайный процесс ξ(t) выражается интегралом Дюамеля

ξ(t)= ∞∫h(τ)xδ(t − τ)dτ .	(4.16)
0

Белый шум xδ(t) с дельтаобразной ковариационной функцией

имеет бесконечную дисперсию (см. также формулу (2.24)). Чтобы этого избежать, при вычислении значений случайного процесса уравнение (4.16) заменяется на

ξ(t)= ∫t	h(τ)x0 (t − τ)dτ ,	(4.17)
0

где x0 (t) – белый шум с ограниченной частотой ωГ , причём в полосе (−ωГ ,ωГ ) должна находиться основная часть мощности

процесса ξ(t). x0 (t) имеет дисперсию D0 = ωπГ и не коррелиро-

ванные в точках tn = n∆t = n	π	значения.
	ωГ

Если заменить интеграл (4.17) конечной суммой с шагом ∆t , получим алгоритм реализации дискретных значений случайного процесса ξ(t) в виде

ξ(n)= ξ(n∆t)= ∆t ∑∞ h(k )x0 (n −k )= ∑∞ ck x(n −k ),		(4.18)
k=0	k=0

где сk = ∆th(k ), x(m) N (0,1).

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

127

Пример 17 [3]. Найдём формулу для моделирования случайного процесса с ковариационной функцией Kξ (τ)= e−ωk τ .

Найдём прежде всего спектральную плотность мощности заданного случайного процесса:

1

∞

1

∞

−ω

τ

−iωτ

S(ω)=

2π

Kξ(τ)e−iωτdτ

=

e

k

dτ

=

−∞∫

2π −∞∫

1

0

−ω

τ

−iωτ

∞

−ω

τ

−iωτ

=

∫e

k

dτ+ ∫e

k

dτ

=

2π −∞

0

=

1

1

e

−ω

τ

−iωτ

0

+

1

e

−ω

τ

−iωτ

∞

=

k

k

−iω

−ωk −iω

2π

ωk

−∞

0

=

1

1

+

1

=

ωk

.

2π

−iω

π(ω2k + ω2 )

ωk

ωk +iω

Здесь разделение интервала интегрирования связано со свой-

ствами функции τ = − τ,τ < 0,

τ,τ ≥ 0.

Корни спектральной функции S(ω) равны ±iωk . Передаточная функция формирующего фильтра, вычисленная по формуле

(4.12), будет равна: K(p)=			∏l (p − p1k )
	C		k=1			, где p	= iω	– корни
			∏m	(p − pk )		1k	1k

k=1

числителя, pk = iω2k – корни знаменателя. Числитель корней не имеет (не зависит от ω), корни знаменателя ±iωk . Согласно теории выберем корень с положительной мнимой частью, это корень

+iωk . Тогда K(p)=

1

(в выражение для K(p) входит

C

p −ωk

коэффициент при мнимой части).

K(iω)

2 = S(ω) даёт

Найдём теперь константу C . Условие

ωk

K

(p)

2

K(iω)

2

= C

1

2

C

π(ω2k + ω2 )=

=

iω−iωk

=

ω2 + ωk2

.

128

Здесь использовано свойство модуля комплексного выражения, именно, если z = a +ib , то z = a2 +b2 . Отсюда С = ωπk .

Вычислим

теперь

импульсную переходную

характерис-

тику

формирующего

фильтра

по

формулам

(4.14) и (4.15).

1 0

t

j

Так

как k =1,r1 =1, p1 = −ωk ,

то

h(t)= ∑∑Cij

e pit = C10e p1t ,

j!

i=1 j=0

1

d1−0−1

1

C10 =

[K(p)(p − p1 )

]

p=p1

= [K(p)(p − p1 )]

p=p1 =

(1−0 −1)!

dp1−0−1

= C p +1ωk (p + ωk )p=ωk = C . Тогда h(t)= Ce−ωkt = ωπk e−ωkt .

Наконец, получим формулы, моделирующие случайный процесс с данной ковариационной функцией (см. формулы (4.18)):

		ωk				и ξ(n)= ∑∞ ck x(n −k), xk N (0,1).
ck = ∆t			e−ωk∆tk =	a	e−ak ,a = ωk ∆t
		π		π		k=0

Пример 18. Осуществим моделирование случайного процесса ξ(t) с ковариационной функцией Kξ(τ)= De−α τ [2δ(τ−α(sign(τ))2 )].

Найдём вначале спектральную плотность. Так как Kξ(τ) – действительная чётная функция, применим формулу (2.14).

При этом	используем свойства	функций					δ(τ) = ∞, τ = 0,
									0, τ ≠ 0
1,τ > 0,					2	∞K
и sign(τ)=	0,τ = 0, В этом случае	S		(ω)=					(τ)cosωτdτ =
			ξ		π 0∫			ξ
−1,τ < 0.

=2 ∞∫De−α τ [2δ(τ)− α(sign(τ))2 ]cosωτdτ =

π0

=2D ∞∫[2e−α τ δ(τ)− αe−α τ (sign(τ))2 ]cosωτdτ =

π0

=4D ∞∫e−ατδ(τ)cosωτdτ − 2αD ∞∫e−ατ(sign(τ))2 cosωτdτ =

π0 π 0

129

=

2D

e

−α 0

cosω 0 −

2αD ∞

−ατ

cosωτdτ,

так как по свойству дель-

π

π

∫e

0

та-функции

0+∫ε f (τ)δ(τ)dτ =

1

f (0). Последний интеграл в выраже-

0

2

нии

для

Sξ(ω)

легко

берётся

по

частям,

именно

∫e

−ατ

cosωτdτ =

α2

e

−ατ

1

ω

−

cosωτ+

sin ωτ .

Тогда

α

2

+ ω

2

α

α

2

Sξ(ω)=

2D

− 2αD

α

=

2Dω2

.

Графики

функций

π

α2 + ω2

π(α2 + ω2 )

π

Kξ(τ) и Sξ(ω) приведены на рис. 4.3.

K(τ)

S(ω)

τ

ω

Рис. 4.3. Графики ковариационной функции и функции спектральной

плотности мощности примера 18

Корень числителя спектральной функции равен нулю, корни знаменателя ± iα. Передаточная функция формирующего фильтра

равна (формула (4.13)): K(p)= C p −piα , так как для знаменателя выбирается один корень с положительной мнимой частью.

Определим

значение

константы

С. По усл

овию

Sξ(ω)=

K(iω)

2 ,

2D

ω2

iω

2

ω2

т.е.

= C

=

.

π

α2 + ω2

iω− iα

α2 + ω2

130

Отсюда C = 2πD . Импульсная переходная характеристика формирующего фильтра аналогична таковой же в примере 17. Именно,

k =1,r =1, p = −α,

h(t)= C e p1t ,

C

= [K(p)(p − p

)]

=

1

1

10

10

1

p=p1

=

p

(p − α)

=

α =

2D

α ,

h(t)=

2D

αe−αt .

C

C

p − α

π

π

p=α

Таким образом, формулы, моделирующие случайный про-

цесс

ξ(t) с заданной

ковариационной

функцией,

имеют

вид

ξ(n)= ∆t ∑∞ h(k )x0 (n −k )= ∑∞ ck x0 (n −k ),

где ck =

∆th(k)=

k=0

k=0

e−ak , a = α∆t, x0 N (0,1).

=

2D

αe−α∆tk =

2D∆t

α2 e−α∆tk =

2Dαa

∆t

π

π

π

4.4.4. Метод канонических разложений. Каноническое пред-

ставление случайного процесса основано на выражении модели случайного процесса ξ(t) в виде детерминированной функции

случайных величин, т.е. в виде

ξ(t)= mξ(t)+ ∑∞ uk xk (t),	(4.19)
k=1

где uk – коэффициенты разложения – случайные величины, математические ожидания которых равны нулю, xk (t) – детерминированные функции, образующие систему функций. Функции xk (t) называются координатными функциями канонического разложе-

ния. Коэффициенты разложения uk		изменяются от реализации к
реализации; необходимо, чтобы uk		были некоррелированными,
т.е. M [uiu j ]= 0	при i ≠ j .
Согласно теории [18] некоррелированность случайных коэффи-
циентов uk	обеспечивается при	выборе в качестве системы

xk (t),k =1,2,...,∞ функций, являющихся решениями интегрального уравнения типа Фредгольма относительно функции x(t):

131

T∫MK(t1,t2 )xk (t2 )dt2 = λk uk (t1 ),	(4.20)
0

где K(t1,t2 ) – заданная ковариационная функция– ядро интегрального уравнения, TM – интервал моделирования случайного процесса, λk – собственные числа данного интегрального уравнения.

Основным достоинством представления (4.20) является возможность моделирования случайного процесса ξ(t) для любого

момента времени t , так как ξ(t) моделируется как функция непре-

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

1. Решение интегрального уравнения (4.20) аналитически найдено лишь для весьма ограниченного набора ковариационных

функций K(t1,t2 ).

2. В соответствии с (4.19) необходимо использовать бесконечное число функций xk (t). Это не может быть реализовано при

практическом моделировании. Если же взять число членов в правой части формулы (4.19) конечным, возникает методическая ошибка.

В работе [18] В.С. Пугачёв предложил метод, в котором случайный процесс ξ(t) моделируется таким образом, что истинная

ковариационная функция K(t1,t2 ) совпадает с ковариационной функцией моделируемого случайного процесса K Д (t1,t2 ) лишь

для заданных дискретных моментов времени. Таким образом, K(t1,t2 )= K Д (t1,t2 ) для конечного множества моментов, а само

моделирование происходит по формуле

ξ(t)= mξ(t)+ ∑N uk xk (t),	(4.21)
k=1

где N – число точек, в которых K(t1,t2 )= K Д (t1,t2 ). Для вычисления по формуле (4.20) необходимо определить ортонормированную систему функций xk (t),k =1, N и некоррелированные случай-

ные величины uk ,k =1, N . При этом накладываются условия лишь на первый и второй моменты случайных величин uk [18]:

132

M [u

]= 0,k

=

, D

= K(t ,t

), x (t)=

1

K(t,t

),

k

1, N

21

21

1

11

1

D1

Dk = K(t1,k ,t2,k )− k∑−1Di xi (tk ),k = 2,3,..., N,

(4.22)

i=1

x

k

(t)=

1

K(t,t

2,k

)− k−1D

x

(t

k

) ,k = 2,3,..., N.

∑ i

i

Dk

i=1

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

уравнениями (4.22).

4.4.5. Метод канонических разложений для стационарных случайных процессов. У стационарных случайных процессовкова- риационная функция, зависящая от разности аргументов

Kξ(t1,t2 )= Kξ(t2 −t1 )= Kξ(τ), связана со спектральной плотностью

мощности парой преобразования Фурье (2.18)-(2.19), поэтому при моделировании стационарных случайных процессов можно ис-

пользовать функцию Sξ(ω).

При условии Kξ(t1,t2 )= Kξ(τ) собственными функциями

xk (t) интегрального уравнения (4.20) (координатными функциями

уравнения (4.19)) являются гармонические функции, т.е. ряд (4.19) превращается в ряд Фурье:

ξ(t)= mξ(t)+ ∑∞ (vk coskωt +uk sin kωt),0 ≤ t ≤ TM ,	(4.23)
k=0

где T	– интервал	моделирования,	ω =	2π	, T	П
	M			TП
M [vk ]= M [uk ]= 0 ,		M [viu j ]= M [uiu j ]= M [viv j ]= 0

i ≠ j, M [viui ]= 0,TП – интервал, равный периоду формулы выбранный так, чтобы ковариационная функция процесса значимо не отличалась от заданной Kξ(τ).

Определим ковариационную функцию процесса (4.23).

Kξ (t1,t2 )= M [ξ(t1 )ξ(t2 )]=

∞	+ uk sin kωt1 ) (vk coskωt2		=
= M ∑(vk coskωt1		+ uk sin kωt2 )
k=0

≥ TM ,

при

(4.23),

(4.23)

133

= ∑∞ (Dv	k	coskωt1 coskωt2 + Du	sin kωt1 sin kωt2 ). Для того чтобы
k=0			k

случайный процесс ξ(t) был стационарным в широком смысле,

необходимо, чтобы Dv	k	= Du	= Dk . Тогда
		k
Kξ(t1,t2 )= ∑∞ Dk (coskωt1 coskωt2				+sin kωt1 sin kωt2 )=
k=0				(4.24)
∞			∞
= ∑Dk coskω(t2			−t1 )= ∑Dk coskωτ.
k=0			k=0

Выражение (4.24) есть разложение ковариационной функции в ряд Фурье с периодом TП = 2ωπ , коэффициентами которого явля-

ются дисперсии Dk . В этом случае, так как Kξ (τ) – чётная функция, то на интервале [−TП ,TП ] коэффициенты ряда Фурье равны:

	D0 =	1	T∫ПKξ (τ)dτ,	Dk =	2	T∫ПKξ (τ)coskωτdτ,	(4.25)
		TП			TП
			0			0

где TП >> τk ,τk – интервал корреляции случайного процесса ξ(t). Если TП >> τk , то с небольшой методической погрешностью уравнения (4.25) можно переписать в виде

D0	=	1		∞∫Kξ (τ)dτ =	π	Sξ (0),
		TП 0			2TП			(4.26)
Dk =	2		∞∫Kξ (τ)coskωτdτ =				π	S(kω).
							TП
	TП 0

Таким образом, дисперсия гармоник формул (4.23) и (4.24) определяется по заданной спектральной плотности мощности

S(ω) с точностью до постоянных множителей.

Обычно функция S(ω) резко убывает и становится малой по достижении некоторой верхней частоты ωВ , поэтому число членов усечённого ряда в (4.23) может быть приближённо оценено

как M = ωωВ = ω2πВ TП . С другой стороны, поскольку общая дис-

134

персия моделируемого случайного процесса равна Dξ = Kξ(0), дисперсия же при усечении формулы (4.23) определяется соотно-

	M
	шением DM = ∑Dk	,	то критерием выбора M			может быть выра-
	k=0
	жение				∞
			DM =	1
	1−				∑Dk < ε ,	(4.27)
			Dξ	Kξ(0)k=M +1

где ε – заданная погрешность при моделировании.

При таких условиях и ограничениях моделирование случайного процесса ξ(t) происходит по формуле

ξ(t)= mξ(t)+ ∑M (vk coskωt +uk sin kωt),			(4.28)
	k=0
где vk и uk –	некоррелированные	случайные	величины с
M [vk ]= M [uk ]= 0 ,	D[vk ]= D[uk ]= Dk .	Распределение случайных

величин vk и uk может быть произвольным, если процесс ξ(t) нормальный, то vk и uk – также нормально распределённые случайные

величины с определёнными первыми и вторыми моментами.

Так как каждое слагаемое в правой части (4.28) представляет собой простую гармонику с частотой kω, то при условии нор-

мальности случайного процесса ξ(t) эта гармоника имеет фазу, равномерно распределённую на интервале [− π,π], и амплитуду, распределённую по закону Рэлея с параметром DR = 2Dk , т.е.

	f (A)=	A	−	A2
				4D
			e	k . В этом случае моделирование по		формуле
		2Dk
	(4.28) может быть заменено моделированием по формуле
				ξ(t)= mξ(t)+ ∑M Ak cos(kωt +ϕk ),		(4.29)
					k=0

где вместо двух случайных величин vk и uk с одинаковым (нормальным) распределением используются две случайные величины с разными законами распределения: амплитуда Ak распределена

135

по закону Рэлея с параметром 2Dk , фаза ϕk распределена равномерно на R[− π,π].

Формулу (4.29) можно ещё упростить, исключив случайные

амплитуды гармоник Ak .	Так как мощность амплитуды каждой
гармоники в (4.28) равна:	M [vk2 ]+ M [uk2 ]= D[vk ]+ D[uk ]= 2Dk , то
формулу (4.29) можно переписать в виде
ξ(t)= mξ(t)+ ∑M Ck cos(kωt +ϕk ),		(4.30)
	k=0

где Ck – неслучайный параметр и Ck = 2Dk .

Формулу (4.30) используют при большом значении M , так как в силу центральной предельной теоремы процесс ξ(t), получаемый по формуле (4.30), сходится к нормальному.

4.4.6. Моделирование марковских случайных процессов (дис-

кретных марковских цепей). Моделирование марковских случай- ных процессов осуществляется на основе метода условных распределений (см. подразд. 4.2). Сам процесс моделирования зависит от порядка марковского процесса. Рассмотрим вначале процесс первого порядка. По определению марков-

ской	цепи,	P{ξ(tn )= in	ξ(t1 )= i1,ξ(t2 )= i2 ,...,ξ(tn−1 )= in−1}=
= P{ξ(tn )= in		ξ(tn−1 )= in−1}. Кроме того, вероятность перехода из
состояния i в состояние j			определяется переходной вероятно-

стью pij , которая должна быть вычислена или задана. Для дис-

кретных однородных марковских цепей с дискретным временем справедлива формула (2.45).

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

f0 (t0 ) начального значения ξ(t0 ) в начальный момент времени t0 . Пусть задана начальная плотность f0 (t0 )= f0 (ξ(t0 ),t0 ).

По этой плотности соответствующим датчиком случайных чисел реализуется значение первого отсчёта ξ1(t0 ) первой реализа-

136

ции случайного процесса ξ(t). Для следующего, второго, отсчёта определяется условная одномерная плотность вероятности

f (ξ ξ )= f (ξ(t ),t						ξ(t		),t		)=	f01 (ξ(t0 ),ξ(t1 ))	=	f01 (ξ(t0 ),t0 ,ξ(t1 ),t1 )
											f0 (ξ(t0 ),t0 )
1	1	0	1	1	1		0		0					f0 (ξ(t0 ),t0 )
и по этой плотности – значение второго отсчёта														ξ1(t1 ) первой

реализации. Третий отсчёт будет иметь условную одномерную плотность вероятности f1 (ξ2 ξ0 ,ξ1 )= f2 (ξ(t2 ),t2 ξ(t1 ),t1 )=

=

f012 (ξ(t0 ),ξ(t1 ),ξ(t2 ))

=

f012 (ξ(t0 ),t0 ,ξ(t1 ),t1

,ξ(t2 ),t2 )

, по которой

f

01

(ξ(t

0

),ξ(t ))

f

01

(ξ(t

0

),ξ(t

))

1

1

реализуется значение третьего отсчёта ξ1 (t2 ).

Величины

f01,

f012 ,..., f01...k в общем случае рассчитываются

по формулам (1.11). В результате получается последовательность чисел ξ1(t0 ),ξ1(t1 ),ξ1(t2 ),...,ξ1(tk ), изображающих первую реализацию марковского случайного процесса ξ(t). Для получения второй

и последующих реализаций повторяются те же операции. Простая дискретная цепь Маркова представляет собой после-

довательность случайных величин ξ(ωi ),i =1,2,... с возможными значениями i0 ,i1,...,ik . Условные вероятности p(ξk ξk−1 ) для каждого значения процесса ξ(tk ) в момент tk , т.е. вероятности пере-

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

Само моделирование дискретных марковских однородных це-

пей происходит по формуле (2.45): p j (k )= ∑n pi (k −1)pij ,

i=1

k =1,2,..., j =1,2,...,n , где n – число состояний цепи. Для запуска процесса моделирования необходимо задать начальное распределение вероятностей pi (0),i =1,2,...,n и матрицу переходных веро-

ятностей pij .

Тогда вероятности всех состояний p1(1), p2 (1),..., pn (1) на

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

137