Основы цифровой обработки сигналов
06. Быстрое преобразование Фурье (БПФ)
БПФ — это алгоритм, который переводит сигнал из временной области в частотную (или наоборот). Он ускоряет вычисления дискретного преобразования Фурье (ДПФ), сокращая количество операций с O(N²) до O(NlogN).
Данный алгоритм применяется повсеместно — от цифровой обработки сигналов, звука и изображений до радиолокации. В аудиотехнике и спектральном анализе БПФ позволяет наглядно увидеть, какие именно частоты (гармоники) присутствуют в сигнале
Расчёт коэффициентов (поворачивающих множителей)
Для восьми и четырёх отсчётов дискретного преобразования Фурье (ДПФ):
wvec8 = exp.(-1*im*2pi.*(0:7)./8)
wvec4 = exp.(-1*im*2pi*(0:3)./4)
Отобразим коэффициенты на комплексной плоскости:
θ = range(0, 2*π, length=100);
x = cos.(θ);
y = sin.(θ);
plot(x, y, lw=1)
scatter!(wvec8,size=(400,400))
scatter!(wvec4,leg=false)
Далее идёт "хитрый" код для формирования матрицы, в столбцах которой будут коэффициенты для вычисления последовательных точек ДПФ:
wvec_long = repeat(wvec8,8);
wmatrix = zeros(ComplexF64,8,8);
wmatrix[:,1] = wvec8[1] .* ones(8);
for n = 1:7
wmatrix[:,n+1] = wvec_long[1:n:8*n];
end
Наконец, создадим тестовый сигнал во временной области, который мы будем преобразовывать:
x = cos.(LinRange(0,2pi,8))
plot(x,m=:c,leg=false)
Вычисление ДПФ для восьми отсчётов сигнала
"Честно" перемножим сигнал поэлементно восемь раз на значения в столбцах матрицы коэффициентов:
multmatrix_8 = wmatrix .* x;
И посчитаем сумму по каждому столбцу. Это даст нам результат вычисления ДПФ по известной формуле:
DFT8 = vec(sum(multmatrix_8, dims=1));
plot(abs.(DFT8),l=:stem,m=:c,leg=false)
Спектр симметричен, что ожидаемо. Но присутствует постоянная составляющая - ненулевой результат на первом отсчёте ДПФ. Это связано с тем, что исследуемая синусоида не совсем точная, и отрицательные значения не доходят до минус единицы. Сумма отсчётов во времени не равна нулю.
Алгоритм "бабочка" вычисления БПФ для сигнала из 8-ми отсчётов
Попробуем по действиям вычислить схему на иллюстрации и сравнить результат с ДПФ:
Вычисление значений вектора на выходе первой секции:
a = zeros(ComplexF64,8);
a[1] = x[1] + x[5];
a[2] = x[1] - x[5];
a[3] = x[3] + x[7];
a[4] = x[3] - x[7];
a[5] = x[2] + x[6];
a[6] = x[2] - x[6];
a[7] = x[4] + x[8];
a[8] = x[4] - x[8];
a
И на выходе второй секции:
b = zeros(ComplexF64,8);
b[1] = a[1] + a[3];
b[2] = a[2] - im*a[4];
b[3] = a[1] - a[3];
b[4] = a[2] + im*a[4];
b[5] = a[5] + a[7];
b[6] = a[6] - im*a[8];
b[7] = a[5] - a[7];
b[8] = a[6] + im*a[8];
b
Результирующий вектор БПФ и его сравнение с ДПФ:
X8 = zeros(ComplexF64,8);
X8[1] = b[1] + b[5];
X8[2] = b[2] + b[6] * wvec8[2];
X8[3] = b[3] - im*b[7];
X8[4] = b[4] + b[8] * wvec8[4];
X8[5] = b[1] - b[5];
X8[6] = b[2] - b[6] * wvec8[2];
X8[7] = b[3] + im*b[7];
X8[8] = b[4] - b[8] * wvec8[4];
plot(abs.(X8),l=:stem,m=:c)
plot!(abs.(DFT8),l=:stem,m=:c,leg=false)
Проверим максимальное значение отклонения между двумя алгоритмами:
err = DFT8 - X8;
maxerr = maximum(abs.(err))
Сравнение со специализированной функцией БПФ из состава EngeeDSP
Мы можем вызывать алгоритм БПФ в одну команду:
myFFT = EngeeDSP.Functions.fft(x)
plot(abs.(myFFT),l=:stem,m=:c,leg=false)
plot!(abs.(X8),l=:stem,m=:c)
Заключительное сравнение результатов:
errFFT = myFFT - X8;
maxerrFFT = maximum(abs.(errFFT))