Моделирование: 2. Получение весовых коэффициентов путем разложения функции спектральной плотности в ряд Фурье

2. Получение весовых коэффициентов путем разложения функции спектральной плотности в ряд Фурье

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

. (2.12)

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

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

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

, (2.13)

где — спектральная плотность белого шума (рис. 2.3). Соотношение (2.13) выражает известный из теории случайных процессов факт: энергетический спектр шума на выходе линейной системы равен произведению энергетического спектра входного шума на квадрат модуля передаточной функции (комплексной частотной характеристики) системы.

Рис. 2.3.

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

. (2.14)

Импульсная переходная характеристика, соответствующая передаточной функции (2.14), равна

. (2.15)

Последнее равенство в формуле (2.15) написано в силу четности функции . Заметим, что система физически не осуществима, так как ее импульсная переходная характеристика, определяемая формулой (2.15), отлична от нуля не только при положительных, но и при отрицательных значениях , причем . Однако в данном случае это обстоятельство не является ограничением.

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

Запишем случайный процесс на входе системы в виде ряда Котельникова:

, (2.16)

где

;

— последовательность независимых нормально распределенных случайных величин с параметрами ; - дисперсия шума .

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

Выразим реакцию системы на воздействие в виде интеграла Дюамеля:

. (2.17)

Отсюда

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

, (2.8)

где

Подставив в (2.17) ряды (2.16) и (2.18), получим

Поскольку , то на основании соотношения ортогональности функций :

где

окончательно будем иметь

(2.19)

Здесь

(2.20)

где

;

- безразмерная частота.

Отметим, что формула (2.19) совпадает с известной формулой прямоугольников для приближенного вычисления интеграла, если шаг дискретизации подынтегральной функции выбрать равным . Приведенный выше вывод показывает, что формула прямоугольников с шагом в применении к интегралу (2.17) дает точный результат, если функции и имеют спектры, ограниченные частотой .

Коэффициенты , как это следует из формулы (2.20), совпадают с коэффициентами Фурье в разложении функции на интервале . Следовательно, при . Как правило, коэффициенты достаточно быстро убывают. Так, например, если функция непрерывна, то для имеется оценка [76]:

где — некоторое положительное число.

Поэтому обычно достаточно ограничиться в (2.19) небольшим количеством членов и тогда можно записать

. (2.21)

Поскольку исходная последовательность стационарна, то статистические свойства последовательности не изменятся, если алгоритм (2.21) записать в виде, аналогичном (2.1), т. е.

, (2.22)

где .

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

, (2.23)

где — дисперсия моделируемого случайного процесса; — некоторое малое число (погрешность).

Неравенство (2.23) основано на том, что сумма квадратов весовых коэффициентов должна быть равна дисперсии моделируемого случайного процесса [см. (2.9)].

В рассматриваемом методе моделирования подготовительная работа состоит в вычислении интеграла (2.20) по известной функции спектральной плотности моделируемого процесса. Подготовительная работа может быть проделана без применения ЦВМ, если интеграл в (2.20) берется в явном виде. Если же интеграл (2.20) не является табличным и не выражается в элементарных или ранее табулированных специальных функциях, то при нахождении коэффициентов может оказаться целесообразным аппроксимировать некоторой функцией, разложение которой в ряд Фурье заранее известно.

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

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

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

Представляет интерес оценка погрешности данного метода моделирования.

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

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

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

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

(2.24)

где — свертка коэффициентов , т. е.

, (2.25)

Операция свертки коэффициентов образуется при возведении в квадрат ряда Фурье (2.24).

Величина относительного среднеквадратического отклонения функции от функции по определению равна

, (2.26)

где

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

Выражение (2.26) легко преобразуется к виду

где

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

— погрешность, обусловленная заменой бесконечного ряда конечным рядом .

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

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

В силу четности функции коэффициенты будут вещественными и равными

где — коэффициенты Фурье в разложении функции на интервале , а — уже используемые ранее коэффициенты [формула (2.24)].

Согласно равенству Парсеваля среднее значение квадрата функции по некоторому промежутку равно сумме квадратов модулей коэффициентов разложения функции в ряд Фурье на этом промежутке, тогда можно записать

Учитывая, что при , окончательно получим

. (2.27)

Нетрудно найти абсолютную величину отклонения функции корреляции , получаемой при данном методе моделирования, от заданной функции корреляции . В самом деле, значения корреляционной функции совпадают со значениями коэффициентов , что следует из сравнения свертки формула (2.25) с (2.10). Поэтому

. (2.28)

Вычисление свертки сводится к перемножению матриц, аналогичных матрицам (2.11).

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

Пример 1. Пусть требуется имитировать на ЦВМ стационарный нормальный случайный процесс, корреляционная функция и энергетический спектр которого имеют вид

. (2.29)

Используя формулу (2.20) при произвольной пока , получим

Тогда

. (2.30)

Интеграл (2.30) является табличным [25]. После несложных преобразований получим

(2.31)

При достаточно большой и выражение для упрощается

. (2.32)

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

Величина погрешности метода будет зависеть при этом от выбора значений и . Для примера выберем , частоту будем отсчитывать на уровне 0,01 от максимума , что равносильно выбору . Значения , , и , рассчитанные при выбранных , и по формулам (2.31), (2.29), (2.25) и (2.28) соответственно, сведены в табл. 2.1.

Таблица 2.1

0	1	2	3	4	5	6	7	8	9	10
0,839	0,358	0,099	0,058	0,027	0,022	-	-	-	-	-
0,988	0,688	0,342	0,192	0,112	0,068	0,024	0,008	0,003	0,001	0,000
1,000	0,682	0,349	0,192	0,118	0,079	0,056	0,042	0,032	0,026	0,021
0,012	0,006	0,007	0,000	0,006	0,001	0,032	0,034	0,029	0,025	0,021
0,990	0,689	0,345	0,194	0,117	0,080	0,056	0,042	0,032	0,026	0,021

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

Для расчета относительной среднеквадратической погрешности при найдем:

Значения при для помещены в табл. 2.1. Используя табличные значения и , получим

Величину второй суммы в формуле (2.27) оценим следующим образом

Следовательно, погрешность за счет ограничения ряда согласно (2.27) оценивается величиной

Общая погрешность .

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

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

Заканчивая рассмотрение этого метода моделирования, заметим, что в основу его положена идея формирующего фильтра, включающая три момента: 1) пропускание белого шума через линейный непрерывный фильтр; 2) подбор такой передаточной функции фильтра, которая обеспечивает энергетический спектр шума на выходе, равный энергетическому спектру моделируемого процесса; 3) дискретизация процессов с целью воспроизведения фильтрации на ЦВМ. Особенностью при этом было то, что формирующий непрерывный фильтр отыскивался для простоты в классе физически неосуществимых линейных фильтров с четной импульсной переходной характеристикой.

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