Материал: ТАУ Лекция ч2

© К.Ю. Поляков, 2009

SX (ω)

Белый шум «содержит» все частоты, но они по-разному преобразуются. Постоянный сигнал (имеющий частоту ω = 0 ) передается на выход системы без изменений. На низких частотах искажения достаточно малые, а высокие частоты подавляются (фильтруются) системой. Это типичный фильтр низких частот (он пропускает низкочастотные сигналы и блокирует высокочастотные). Отметим очень важный факт: поскольку высокие частоты подавляются, отклонения спектра входного сигнала от равномерного спектра белого шума в этой области не будут существенно влиять на спектр выхода. На этой идее основано компьютерное моделирование случайных процессов (см. далее).

3.4. Моделирование случайных сигналов

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

3.4.1. Случайные числа

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

математической формуле на основе предыдущего xk (или нескольких предыдущих). Например, во многих системах программирования используется линейный конгруэнтный метод:

xk+1 = (axk +c) mod m .

Здесь a , c и m – некоторые целые числа. По этой формуле вычисляются псевдослучайные числа, равномерно распределенные8 на интервале от 0 до m −1. Для краткости будем далее называть псевдослучайные числа просто случайными.

Параметры a , c и m нужно выбрать так, чтобы полученная последовательность содержала как можно меньше закономерностей и как можно дольше не повторялась. Это отдельная математическая проблема, которая до сих пор не имеет однозначного решения. Например, в функции rand среды MATLAB используется генератор с параметрами9

a = 75 =16807, c = 0, m = 231 −1 = 2147483647 .

8Теоретически, конечно. Реальное распределение всегда немного неравномерное.

9См. http://www.mathworks.com/company/newsletters/news_notes/pdf/Cleve.pdf.

26

© К.Ю. Поляков, 2009

Последовательность начинается с некоторого начального числа (англ. seed – семя, зерно). Чтобы получить в точности такую же последовательность (например, чтобы повторить полученные ранее результаты моделирования), нужно начать последовательность с того же числа.

Как показано в разд. 1.6, при сложении равномерно распределенных чисел распределение их суммы стремится к нормальному. Поэтому простейший способ получения нормального распределения – сложить некоторое количество случайных чисел, полученных с помощью стандартного датчика. Часто используют суму 12 чисел, равномерно распределенных на интервале [0;1], потому что это сразу дает случайную величину с единичной дисперсией.

В MATLAB есть встроенная функция randn, которая генерирует последовательность случайных чисел с нормальным распределением (с нулевым средним и единичной дисперсией) на основе более сложного алгоритма10.

3.4.2. Формирующий фильтр

На практике обычно известна спектральная плотность SX (ω) и требуется смоделировать процесс, имеющий такую спектральную плотность.

Вспомним, что спектральная плотность неотрицательна для любой частоты. Тогда функция SX (s) , полученная при подстановке ω2 = −s2 в SX (ω) , неотрицательна на мнимой оси, то есть при s = jω для всех ω . Можно доказать, что в этом случае ее можно представить в виде произведения SX (s) = F(s)F(−s) , то есть в форме (7). При этом всегда можно выбрать передаточную функцию F(s) так, чтобы она была устойчивой (не имела полюсов в правой полуплоскости) и минимально-фазовой (не имела нулей в правой полуплоскости). Такой переход от

SX (s) к F(s) называется факторизацией (англ. разложение на множители).

Как следует из (7), если на вход звена с передаточной функцией F(s) подать единичный белый шум, процесс на выходе будет иметь заданную спектральную плотность SX (ω) . Функ-

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

Проще всего моделировать процессы с дробно-рациональной спектральной плотностью. Например, одна из моделей морского волнения (подробности см. в главе 4) описывается спектральной плотностью

S(ω) = 4Drαω2 ,

ω4 + 2(α2 − β2 )ω2 +(α2 + β2 )2

где Dr – дисперсия волновой ординаты, α – коэффициент затухания и β – частота максимума спектра. Заменяя ω2 на −s2 , получаем

S(s) = −4Drα s2 .

s4 −2(α2 −β2 )s2 +(α2 + β2 )2

Очевидно, что формирующий фильтр будет иметь передаточную функцию вида

F(s) = γ s ,

s2 +δ1s +δ0

так что

10 См. http://www.mathworks.com/company/newsletters/news_notes/clevescorner/spring01_cleve.html.

27

© К.Ю. Поляков, 2009

F(s) F(−s) = −γ 2 s2 .

s4 −(δ12 −2δ0 )s2 +δ02

Приравнивая коэффициенты числителя и знаменателя S(s) и F(s)F(−s) , сразу находим:

В более сложных случаях факторизация выполняется с помощью численных методов. Нужно разложить на простейшие сомножители числитель и знаменатель S(s) и включить в F(s) только те множители, корни которых находятся в левой полуплоскости (в числителе допускаются корни на мнимой оси). Посмотрим, как работает этот метод на том же примере. Числитель спектральной плотности S(s) имеет двойной корень в точке s = 0 , его можно представить в симметричном виде

−4Drα s2 = 2 Drα s 2 Drα (−s) .

Знаменатель имеет корни в точках α ± jβ и −α ± jβ . Учитывая, что α > 0 , определяем, что корни −α ± jβ расположены в левой полуплоскости. Поэтому знаменатель F(s) представляет собой произведение (s +α + jβ)(s +α − jβ) = s2 + 2αs +α2 + β2 . Таким образом,

Такой ж результат был получен ранее методом неопределенных коэффициентов.

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

3.4.3. Случайный ступенчатый сигнал вместо белого шума

Как мы видели, на компьютере несложно получить последовательность случайных (точнее – псевдослучайных) чисел с равномерным или нормальным распределением. По этим числам можно построить ступенчатый сигнал, фиксируя каждое значение в течение некоторого времени τk :

Теоретически эти числа некоррелированы11; при этом можно показать, что корреляционная функция RX (τ) ступенчатого сигнала – треугольная (см. рисунок справа). При τ = 0 она равна дисперсии D последовательности случайных чисел, а при τ >τk обращается в нуль (по-

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

28

© К.Ю. Поляков, 2009

тому что моменты времени t и t +τ находятся на разных интервалах и, следовательно, соответствующие значения некоррелированы). Число τk называют интервалом корреляции – так назы-

вается интервал, после которого можно считать корреляционную функцию (примерно) равной нулю.

Взяв преобразование Фурье от корреляционной функции

Вычисляя предел этой функции при ω → 0 , находим lim SX (ω) = Dτk , так что при выборе

ω→0

D =1/τk это значение равно 1 (как у белого шума). Заметим, что SX (ω) = 0 , когда cosωτk =1 , то есть ωτk = 2π m при любом целом m . Форма спектральной плотности показана на рисунке ниже (здесь и далее принимается D =1/τk ):

стоянной на всех частотах. Тем не менее, при уменьшении интервала корреляции τk «колокол» расширяется, и для низких частот можно считать, что SX (ω) ≈1. В пределе при τk → 0 спектр стремится к равномерному спектру единичного белого шума. Далее будет показано, что при грамотном выборе τk такой сигнал можно использовать в качестве источника вместо белого

шума.

Для примера предположим, что нужно получить сигнал со спектральной плотностью

SY (ω) = ω21+1 , то есть формирующий фильтр имеет передаточную функцию F(s) = s 1+1 . В ка-

честве входного сигнала для этого звена будем использовать описанный выше ступенчатый сигнал при τk = 0,5 с. На рисунке приведены графики спектральной плотности ступенчатого

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

29

© К.Ю. Поляков, 2009

ω

По графику видно, что в существенной полосе частот (где частотная характеристика звена ненулевая) спектр входного сигнала существенно неравномерный, поэтому желаемый и фактический спектры на выходе системы немного различаются в области высоких частот. Приняв τk = 0,01 c , имеем совершенно другую картину:

S(ω)

1

τk = 0,01 c

ω

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

Очевидно, что при выборе τk нужно учитывать частотные свойства формирующего фильтра, точнее, полосу частот, где его частотная характеристика достаточно отличается от нуля. Для этого используют понятие полосы пропускания системы ωb – так называется частота,

для которой амплитудная частотная характеристика уменьшается на 3 дБ (децибела) в сравнении с максимальным значением (составляет примерно 0,708 от максимума). Разработчики MATLAB рекомендуют при моделировании использовать значение

30