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

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

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

3.2.3. Использование дискретного преобразования Фурье

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

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

F(ω) = ∞∫ f (t) e− jωt dt .

−∞

Оно имеет смысл для любой детерминированной (неслучайной) функции f (t) , которая абсолютно интегрируема, то есть интеграл от ее модуля на всей оси сходится:

∞∫ f (t) dt < ∞ .

−∞

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

Однако если рассмотреть усеченный процесс xT (t) , равный реализации x(t) случайного процесса на интервале [−T;T ] и нулю вне этого интервала, для него можно найти преобразование Фурье:

FX (ω) = ∫xT (t) e− jωt dt = ∫x(t) e− jωt dt .

−∞ −T

Квадрат модуля этой функции, деленный на ширину интервала 2T , характеризует среднюю мощность сигнала на частоте ω . В пределе, при T → ∞, мы должны получить спектральную плотность мощности. Так как x(t) – это только одна реализация случайного процесса, в окончательной формуле нужно использовать усреднение по ансамблю (математическое ожидание):

SX (ω) = lim E{FX (ω) 2 }.

T →∞ 2T

21

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

При реальных измерениях мы знаем только одну реализацию случайного процесса x(t) на интервале [0;T ] , поэтому усреднение по ансамблю чаще всего невозможно. Тогда для оценки спектральной плотности можно использовать предыдущую формулу без усреднения::

Предположим, что известны его значения xk = x(k ∆t) при t = k ∆t для k = 0,1,..., N −1, так что интервал [0;T ] разделен на N подынтервалов шириной ∆t (поэтому T = N ∆t ). Тогда интегрирование можно приближено заменить суммой:

Для оценки спектра в теории обработки сигналов обычно используют сетку частот (в герцах) fm = Nm∆t , m = 0, 1, K, N −1

с шагом ∆f = N 1∆t . В теории управления принято строить спектры как функции угловой час-

тоты (в радианах в секунду), которая получается из «обычной» частоты умножением на 2π :

ωm = N2π∆t m, m = 0, 1, K, N −1.

Для частоты ωm получаем

k=0

где через X m обозначена сумма, называемая дискретным преобразованием Фурье (ДПФ):

Заметим, что эта величина – комплексная, содержащая как вещественную, так и мнимую части. Легко подсчитать, что при расчете ДПФ по этим формулам для N частот количество операций сложения и умножения будет пропорционально N 2 (обозначается O(N 2 ) ). Это значит, что если N увеличивается, скажем, в 10 раз, то количество операций – примерно в 100 раз. Для больших N , особенно при анализе сигналов в реальном времени, такие расчеты выполняются

недопустимо долго.

Для быстрого вычисления ДПФ были разработаны специальные алгоритмы, которые называются быстрым преобразованием Фурье (БПФ). Они позволили сократить количество операций с O(N 2 ) до O(N log N) . В функции fft среды MATLAB используется модификация алгоритма БПФ, предложенного Дж. Кули и Дж. Тьюки. Этот алгоритм наиболее эффективен, если число отсчетов N представляет собой степень двойки ( N = 2p при целом p > 0 ). Заметим, что если это не так, всегда можно дополнить ряд нулями до ближайшей степени двойки.

22

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

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

Согласно теореме Котельникова-Шеннона, по дискретным измерениям с периодом ∆t

можно восстановить частотные свойства сигнала только до частоты fmax = 21∆t (или до соот-

ветствующей угловой частоты ωmax = ∆πt , которая называется частотой Найквиста6). Поэтому только оценка спектра на частотах ω0 , K, ωN / 2 дает нам практически полезную информацию7.

Для сглаживания спектральной плотности так же, как и в методе Блэкмана-Тьюки, используются окна. Только теперь на весовую функцию умножается не оценка корреляционной функции, а сама реализация на интервале [0;T ] :

Далее дискретное преобразования Фурье вычисляется для отсчетов взвешенной функции, то есть, вместо (6) получаем

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

6Иначе говоря, синусоиду нужно измерять более 2 раз за период.

7Можно доказать, что X N −k = X k* , где звездочка обозначает комплексно-сопряженную величину.

23

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

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

T T

∫q2 wh2 (t) dt = ∫1 dt = T .

0 0

Несложно подсчитать, что для окна Хэмминга из этого условия следует

3.3. Прохождение случайных сигналов через линейные системы

Существует два подхода к исследованию систем управления при случайных возмущениях:

1)вероятностный – на основе плотностей распределения вероятностей;

2)статистический – с помощью усредненных характеристик: математического ожидания, дисперсии, корреляционной функции, спектральной плотности.

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

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

стью SX (ω) .

Прежде всего, отметим, что при стационарном случайном входе выход y(t) линейной стационарной системы (у которой все характеристики не зависят от времени) – тоже стационарный случайный процесс. Для процесса y(t) требуется найти

•математическое ожидание y ;

•дисперсию Dy ;

•корреляционную функцию RY (τ) ;

•спектральную плотность SY (ω) .

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

y = kF x, kF = lim F(s) .

s→0

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

24

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

Остальные характеристики удобнее определять с помощью спектральной плотности выхода. Вспомним, что спектральная плотность – это плотность распределения мощности сигнала по частотам. Сначала рассмотрим, как изменяется мощность гармонического сигнала x(t) = Asinω0t , когда он проходит через линейную систему. Из классической теории автомати-

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

y(t) = Bsin(ω0t +ϕ) ,

причем его амплитуда B определяется по частотной характеристике F( jω0 ) :

B = A F( jω0 ) .

Мощность гармонического сигнала (средний квадрат) пропорциональна квадрату амплитуды:

B2 = A2 F( jω0 ) 2 = A2 F( jω0 ) F(− jω0 ) .

В последнем равенстве использовано свойство комплексного числа F( jω0 ) – квадрат его модуля равен произведению этого числа на комплексно-сопряженное, то есть на F(− jω0 ) .

Таким образом, мы с помощью простых рассуждений вышли на очень важный результат: при гармоническом входе с частотой ω0 мощность сигнала на выходе линейной системы равна мощности входного сигнала, умноженной на F( jω0 )F(− jω0 ) . Учитывая, что это свойство

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

SY (ω) = SX (ω) F( jω) F(− jω) .

Соответствующая корреляционная функция RY (τ) может быть найдена как обратное преобразование Фурье от SY (ω) . Для вычисление среднего квадрата процесса y(t) нужно проинтегрировать SY (ω) :

y2 = 1 ∞∫SY (ω) dω .

π 0

Если процесс центрированный, средний квадрат совпадает с дисперсией, то есть Dy = y2 .

В общем случае для вычисления дисперсии нужно использовать равенство Dy = y2 −(y)2 . Выделим один важный случай, когда входной сигнал – это единичный белый шум с посто-

янной спектральной плотностью SX (ω) =1 (белый шум единичной интенсивности). Тогда получаем

график которой показан на рисунке:

25