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

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

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

.

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

.

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

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

.

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

.

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

, где .

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

.

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

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

.

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

,

(4)

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

.

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

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

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

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

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

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

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

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

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

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

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

.

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

, где .

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

.

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

.