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

Быстрое преобразование Фурье (БПФ, FFT — Fast Fourier Transform) — это алгоритм быстрого вычисления дискретного преобразования Фурье (ДПФ), который широко применяется в цифровой обработке сигналов (ЦОС) для анализа и обработки в частотной области. В ЦОС часто требуется переводить сигналы из временной области в частотную (и обратно), чтобы:

• Анализировать спектр сигнала
• Фильтровать сигналы
• Сжимать данные
• Реализовывать алгоритмы модуляции (OFDM в Wi-Fi, 5G)
• Обрабатывать радиолокационные и медицинские сигналы

Рассмотрим, как применять функцию БПФ из подключаемой библиотеки FFTW.jl для задачи анализа спектра сигнала:

using FFTW

Простой сигнал и его спектр

Познакомимся с особенностями применения функции fft для анализа спектра на примере действительного синусоидального сигнала длительностью 0.2 сек, с частотой дискретизации 2 кГц, амплитудой 0.5 и частотой повторений 100 Гц:

fs = 2000;
dt = 1/fs;
t = 0:dt:0.2-dt;
sinewave = 0.5 * sin.(2*pi*t*100);
plot(t,sinewave)

Стоит отметить длину нашего дискретного сигнала sinewave - 400 отсчётов. Это важно для дальнейшего анализа.

nsamp = length(t)

400

Применим к нему функцию fft и визуализируем результат:

A = fft(sinewave);
plot(A)

Выход функции БПФ - это комплексный вектор A с таким же количеством отсчётов, как и у исходного сигнала во временной области.

Для оценки амплитудного спектра сигнала необходимо взять отсчёты выхода БПФ по модулю:

plot(abs.(A))

Мы наблюдаем два пика, соответствующих положительной и отрицательной области частот. Но такое представление не даёт нам достоверной информации ни о частоте, ни об амплитуде этих пиков.

Для более явного отображения спектра воспользуемся функцией fftshift для сдвига выходного вектора БПФ, а также зададим вектор частот, который отложим по оси Х. Мы наблюдаем спектр сигнала в так называемой первой зоне Найквиста, которая простирается от -fs/2 до +fs/2. Разрешение по частоте (RBW, resolution bandwidth), или шаг частотной сетки, мы можем узнать, поделив весь рассматриваемый частотный диапазон (либо от -fs/2 до +fs/2, либо от 0 до fs) на количество точек БПФ:

B = fftshift(A);
println("Разрешение по частоте (RBW) = ", fs/nsamp, " Гц")
freq_vec = -fs/2:fs/nsamp:fs/2-df;
plot(freq_vec, abs.(B))

Разрешение по частоте (RBW) = 5.0 Гц

Теперь мы видим два пика на частотах -100 Гц и +100 Гц. Но их амплитуды не соответствуют амплитуде рассматриваемой синусоиды.

Для точной оценки энергетики сигнала на его спектре мы нормируем выходной вектор БПФ по количеству отсчётов, то есть делим на nsamp. Но так как половина энергии нашего действительного сигнала "ушла" в область отрицательных частот, нам надо умножить результирующий спектр на 2. Ну и для отображения воспользуемся линейчатым графиком (stem) с круглыми маркерами в области положительных частот от 0 до 400 Гц:

S = 2*B/nsamp;
println("Разрешение по частоте (RBW) = ", fs/nsamp, " Гц")
plot(freq_vec, abs.(S), 
        line=:stem,
        marker=:circle,
        xlim = (0, 400),
        title = "Амплитудный спектр синусоиды")

Разрешение по частоте (RBW) = 5.0 Гц

Мы наблюдаем пик на частоте 100 Гц амплитудой 0.5, что точно соответствует параметрам исходного синусоидального сигнала.

Изменение разрешения по частоте

Теперь попробуем изменить количество отсчётов входного сигнала ("окно" БПФ). Подадим на вход функции fft первые 80 отсчётов синусоиды. В этом случае на выходе алгоритма БПФ так же будет 80 комплексных отсчётов, а значит, изменится шаг частотной сетки. Сигнал стал короче в 5 раз, следовательно, разрешение по частоте стало 25 Гц:

winsize = 80;
rbw80 = fs/winsize;
freq80 = 0:rbw80:fs - rbw80;
S80 = (2*fft(sinewave[1:winsize]))./winsize;
println("Разрешение по частоте (RBW) = ", rbw80, " Гц")
plot(freq80, abs.(S80), 
        line=:stem,
        marker=:circle,
        xlim = (0, 400))

Разрешение по частоте (RBW) = 25.0 Гц

Изменение частотной сетки на оценку амплитуды и частоты сигнала не повлияло. Дискрет в 25 Гц позволяет нам "точно попасть" в спектральный компонент на 100 Гц. Но что произойдёт, если шаг спектра не будет кратен частоте оцениваемой синусоиды?

winsize = 110;
rbw110 = fs/winsize;
freq110 = 0:rbw110:fs - rbw110;
S110 = (2*fft(sinewave[1:winsize]))./winsize;
println("Разрешение по частоте (RBW) = ", rbw110, " Гц")
plot(freq110, abs.(S110), 
        line=:stem,
        marker=:circle,
        xlim = (0, 400))

Разрешение по частоте (RBW) = 18.181818181818183 Гц

Мы наблюдаем два пика на частотах 91 Гц и 109 Гц. Это явление утечки спектра (spectral leakage): энергетика сигнала на 100 Гц была захвачена основным и боковыми лепестками (см. теорию ЦОС) оконной функции на соседних дискретах частотной сетки.

Дополнение нулями при БПФ

Zero padding (добавление нулей) при БПФ — это техника, которая заключается в дополнении исходного сигнала нулями до нужной длины перед выполнением БПФ. Добавление нулей увеличивает кол-во точек БПФ, уменьшая шаг частотной сетки и делая спектр визуально более гладким. Но это не добавляет новой информации о сигнале — просто интерполирует существующий спектр.

Дополним нашу синусоиду нулями до длины вектора в 2000 отсчётов:

winsize = 2000;
rbw2000 = fs/winsize;
freq2000 = 0:rbw2000:fs - rbw2000;
numofzeros = winsize - length(sinewave);
println("Добавлено нулей = ", numofzeros)
S2000 = (2*fft([sinewave; zeros(numofzeros)]))./winsize;
println("Разрешение по частоте (RBW) = ", rbw2000, " Гц")
plot(freq2000, abs.(S2000), 
        line=:stem,
        marker=:circle,
        xlim = (0, 400))

Добавлено нулей = 1600
Разрешение по частоте (RBW) = 1.0 Гц

Можно заметить, что пик спектра находится на уровне 0.1, а не 0.5. Это связано с тем, что сигнал стал длиннее в пять раз за счёт нулей, и его энергетика, соответственно, также была ослаблена в пять раз.

Классический алгоритм работает эффективнее, когда число точек БПФ является степенью 2, и добавление нулей полезно для выравнивания длины и визуализации, но оно не увеличивает истинное разрешение по частоте. Для реального улучшения спектрального анализа нужно увеличивать длительность записи сигнала.

Интерактивный анализ спектра

В заключении, применим техники оценки спектра для сигнала, представляющего собой сумму трёх синусоид с частотами 80 Гц, 120 Гц и 245 Гц, и амплитудами 10, 5 и 2.4 соответственно. Отобразим форму результирующего вектора во временной области:

three_sines = [10 5 2.4] .* sin.(2*pi*t*[80 120 245]);
sig = sum(three_sines, dims = 2);
sig = vec(sig);
plot(t,sig)

А теперь воспользуемся функционалом скриптов Engee для создания маски кодовой ячейки. Мы можем изменять значение разрешения по частоте (RBW) при помощи инструмента slider.

В коде ячейки происходит расчёт количества выходных точек БПФ, сравнение его с количеством точек дискретного сигнала во времени, и в случае превышения одного над другим выполняется или же дополнение нулями, или же усечение входного "окна" алгоритма. Затем вычисляется и отрисовывается амплитудный спектр.

RBW = 10 # @param {type:"slider",min:1,max:20,step:1}
nfft = Int64(round(fs/RBW));
println("Точек БПФ:  ", nfft,
"    Отсчётов сигнала:  ", nsamp);

if nsamp >= nfft
        win = sig[1:nfft];
        Y = 2*fft(win)./nfft;
        println("Усечение входного окна")
else
        win = [sig; zeros(nfft - nsamp)];
        Y = 2*fft(win)./nsamp;
        println("Дополнение нулями: интерполяция спектра")
end

ff = (0:nfft-1)*fs/nfft;
plot(ff, abs.(Y), 
        line=:stem, 
        xlim = (0, 400), 
        marker =:circle,
        lab = "Сумма синусоид",
        title = "Амплитудный спектр сигнала")

Точек БПФ:  200    Отсчётов сигнала:  400
Усечение входного окна

Попробуйте поэкспериментировать с различными значениями RBW, от 1 до 20 Гц, и посмотреть на точность оценки частоты и амплитуды отдельных спектральных компонент сигнала.

Вывод

Мы познакомились с особенностями применения функции fft библиотеки FFTW.jl для оценки спектра дискретных сигналов, а также узнали, как применять макси кодовых ячеек и выполнять задачи цифровой обработки сигналов в интерактивном режиме.