Быстрое преобразование Фурье (БПФ)

Для процесса с дискретным временем { x_k }, x_k = x (k Δ t), k =1,…, n сохраняются все введенные выше определения и результаты с той разницей, что интегралы заменяются суммами, а частоту (аргумент f спектральной плотности) имеет смысл рассматривать только на интервале [0 F_N ], где F_N = 1/(2Δ t) - частота Найквиста. В частности, прямое преобразование Фурье последовательности { x_k } – непараметрическую оценку спектральной плотности F (f) - можно вычислить по формуле

обратное – по формуле

Вычисления по каждой из этих формул требуют n ² операций. В 1966 г. Мак-Кован предложил геометрическую интерпретацию преобразования Фурье как некоторого поворота вектора x в n – мерном пространстве и показал, что его можно разложить на последовательность поворотов вокруг каждой из осей, что требует всего 4 n операций. Этот подход получил название быстрого преобразования Фурье, БПФ. Сейчас такие быстрые версии разработаны для всех основных интегральных преобразований. Благодаря технике БПФ при корреляционном анализе оказалось выгодным находить корреляционную функцию как обратное преобразование Фурье спектральной плотности, а не наоборот, как это делалось раньше. Пример вычисления БПФ в ИМС MatLab приведен в Документе 7.1 и проиллюстрирован на рис.7.1.

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

При анализе сигнала s (t) на фоне аддитивного шума ξ(t) большое значение имеет безразмерный параметр snr (signal to noise ratio) – отношение сигнал/шум. В случае дискретного времени его определяют как отношение дисперсий:

Эту величину часто выражают в децибелах (дБ):

дБ.

Документ 7.1. Построение спектральной плотности с помощью БПФ

Спектральная плотность гармонического сигнала с частотой f ₀ на фоне аддитивного шума. Файл sp_my_0.m clear; clc; del=0.005; Nyq=1/2/del; n=200; T=n*del; t=0:del:T; f0=30; %Несущая частота (Гц) fi=pi/3/f0; s=0.5; %СКО аддитивного шума y0=sin(2*pi*f0*t+fi); y=y0+s*randn(size(y0)); snr=sqrt(2)/0.5; %отношение сигнал/шум SNR=20*log(snr); %отношение сигнал/шум в Дб S0=abs(fft(y)).^2; %спектральная плотность мощности f=linspace(0,Nyq,n/2+1); subplot(2,1,1); plot(t,y,'LineWidth',2); grid; title('Синусоида с шумом','FontName','Courier New Cyr'); xlabel('Время (с)','FontName','Courier New Cyr'); ylabel('Сигнал (мВ)','FontName','Courier New Cyr'); axis([0 T min(y)*1.1 max(y)*1.1]); subplot(2,1,2); plot(f,S0(1:n/2+1)*22'LineWidth',2); grid; title('Спектральная плотность мощности','FontName','Courier New Cyr'); xlabel('Частота (Гц)','FontName','Courier New Cyr'); ylabel('СПМ (мВ^2)','FontName','Courier New Cyr'); axis([0 Nyq 0 max(S0)*1.1]);

Для случайного ряда длиной n БПФ дает оценку СПМ S (f), значения которой представляют собой практически не коррелированные случайные величины с дисперсией порядка S ²(f), которая, очевидно, не убывает с ростом n. Основной выход состоит в том, чтобы вычислить оценки для набора интервалов, а затем усреднить. В ИМС MatLab широко используется функция spectrum, в которой реализован алгоритм Уэлча. В этой процедуре данный отрезок временного ряда разбивается на M (по умолчанию M =8) частично пересекающихся отрезков. В каждом отрезке удаляется линейный тренд, для каждого из остатков вычисляется БПФ с использованием одного из сглаживающих окон (по умолчанию используется окно Хеннинга) и результаты усредняются. Удаление линейного тренда приводит к небольшому смещению, поэтому результаты такого анализа могут несколько отличаться от результатов работы других спектроанализаторов.