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

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

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

.

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

.

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

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

.

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

.

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

, где .

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

.

Для оценки спектра в теории обработки сигналов обычно используют сетку частот (в герцах)

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

.

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

                                           ,

           (6)

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

.

Заметим, что эта величина – комплексная, содержащая как вещественную, так и мнимую части.

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

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

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

Согласно теореме Котельникова-Шеннона, по дискретным измерениям с периодом  можно восстановить частотные свойства сигнала только до частоты  (или до соответствующей угловой частоты , которая называется частотой Найквиста [1]). Поэтому только оценка спектра на частотах  дает нам практически полезную информацию[2].

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

1) с помощью БПФ (функция fft в Matlab) найти массив ;

2) взяв первую половину этого массива, рассчитать соответствующие значения  для частот, не превышающих частоту Найквиста ;

3) для каждой частоты  найти оценку спектральной плотности мощности по формуле: .

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

Для этого случая окно Хэмминга на интервале  принимает вид

.

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

                              , где .

               

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

.

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

.