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

Авторы: Чачба А.Н., Костоев Р.С.

Быстрое преобразование Фурье (БПФ)
Последовательный алгоритм
Последовательная сложность	$O(N \log N)$
Объём входных данных	$N$
Объём выходных данных	$N$
Параллельный алгоритм
Высота ярусно-параллельной формы	$O(\log N)$
Ширина ярусно-параллельной формы	$O(N)$

1 ЧАСТЬ. Свойства и структура алгоритмов

1.1 Общее описание алгоритма

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

$\hat{f}(\omega)=\frac{1}{\sqrt{2\pi}}\int\limits_{-\infty}^{\infty}f(x)e^{-ix\omega}\,dx$

Дискретное преобразование Фурье, в свою очередь, если аналог непрерывного преобразования Фурье, но для дискретного сигнала содержащего $N$ отсчетов. Широко применяет в цифровой обработке сигналов, теории вероятностей, криптографии и акустике. Преобразование Фурье обратимо, причем обратное преобразование имеет практически ту же форму, что и прямое. Преобразование Фурье имеет сложность $O(N^2)$, но существует быстрый вариант преобразование Фурье со сложностью $O(N\log{N})$.

1.2 Математическое описание алгоритма

Пусть исходный сигнал имеет значения $x_n,\quad n = 0,\dots,N-1$, тогда дискретное прямое преобразование Фурье (ДПФ) имеет вид:

$X_k = \sum_{n=0}^{N-1}x_ne^{-\frac{2\pi i}{N}kn},\quad k = 0, \dots, N-1$

Обозначим $\varepsilon_{N} = e^{-\frac{2\pi i}{N}}$, тогда ДПФ можно перезаписать в матричной форме:

$\bar X = A\bar x$

где матрица $A = \{e^{-\frac{2\pi i}{N}(i - 1)(j - 1)}\}_{i,j=1}^{N}$

1.3 Вычислительное ядро алгоритма

Пусть, для простоты $N = km$, тогда рекурсивная реализация преобразования Фурье, за счет ${k}$ рекурсий на первом этапе и $m$ на последнем этапе, имеет суммарную сложность $O(Nk + Nm + km) = O(N(k + m))$.

В случае, например $N = 2^n$ сложность БПФ составляет $O(N\log{N})$.

В общем же случае, когда $N = \prod_{i=1}^np_i$ сложность БПФ составляет $O(N(\sum_{i=1}^np_i))$.

1.4 Макроструктура алгоритма

Макроструктура БПФ для случая $N = km$ описывается рекурсивно:

1. $k$ независимых преобразований векторов меньшей размерности $m$
2. Умножение элементов на поворотные коэффициенты ($N$ умножений)
3. $m$ обратных преобразований векторов размерностей $k$

1.5 Схема реализации последовательного алгоритма

Рекурсивный метод для случая $N = 2^k$, без оптимизации на C++:

#include <vector>
#include <complex>

using namespace std;

typedef complex<double> cd;
typedef vector<cd> vcd;

vcd fft(const vcd &as) { 
  int n = as.size();
  if (n == 1) return vcd(1, as[0]);

  vcd w(n); // Calculate roots
  for (int i = 0; i < n; i++) {
    double alpha = 2 * M_PI * i / n;
    w[i] = cd(cos(alpha), sin(alpha));
  }

  vcd A(n / 2), B(n / 2);
  for (int i = 0; i < n / 2; i++) {
    A[i] = as[i * 2];
    B[i] = as[i * 2 + 1];
  }
  vcd Av = fft(A);
  vcd Bv = fft(B);
  vcd res(n);
  for (int i = 0; i < n; i++)
    res[i] =   Av[i % (n / 2)] +
        w[i] * Bv[i % (n / 2)];
  return res;
}

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

1.6 Последовательная сложность алгоритма

Алгоритм состоит из трех этапов, следовательно если $N = \prod_{i = 1}^np_i$, то общая сложность составляет порядка $O(N\sum_{i=1}^np_i)$ операций.

Пусть существует некоторое фиксированное число $P$, такое что:

$p_i \leq P \quad \forall i = 1...n,$

тогда, учитывая, что $n \leq \log_{P}{N}$, сложность алгоритма можно записать в виде $O(N\log{N})$.

1.7 Информационный граф

Рисунок 1. Информационный граф преобразования Кули-Тьюки для $N = km = 25, k = 5, m = 5$. Желтые вершины - умножение на поворотные множители.

1.8 Ресурс параллелизма алгоритма

В силу независимости всех рекурсивных вызовов и того факта, что путь максимальной длины в графе (так называемый критический путь) имеет порядок $O(\log{N})$, можно утверждать, что параллельная сложность алгоритма составит $O(\log{N})$.

Замечание: формально основание логарифма на каждом рекурсивном шаге зависит от множителя в факторизации $N$, потому использована О-символика.

При условии ограниченности количества доступных вычислительных узлов рекомендуется выбирать множители в разложении числа $N$ близкими к числу доступных узлов, таким образом доступные ресурсы будут использованы максимально эффективно.

Итого алгоритм относится к логарифмическому классу алгоритмов по параллельной сложности.

1.9 Входные и выходные данные алгоритма

Входные данные: Чаще всего для обработки сигналов в качестве входных данных для БПФ подается вектор размерности $N$ вещественных элементов. Но БПФ работает и для случая элементов над комплексым полем. Таким образом, например, можно экономить на количестве применений БПФ в некоторой конкретной задаче путем приведения двух векторов вещественных чисел к одному вектору комплексных чисел с вещественной частью равной первому вектору и комплексной частью равной второму вектору. Выходные данные: Вектор размерности $N$ комплексных чисел.

1.10 Свойства алгоритма

Матрица Вандермонда преобразования $A = \{e^{-\frac{2\pi i}{N}(i - 1)(j - 1)}\}_{i,j=1}^{N}$ такая, что:

$A^{-1} = \frac{1}{N}A^*$

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

2 ЧАСТЬ. Программная реализация алгоритма

2.7 Существующие реализации алгоритма

• Самая популярная библиотека для работы с преобразованиями Фурье это FFTW
• Реализация от Intel в рамках Math Kernel Library
• Реализация ДПФ на графических картах от NVidia - cuFFT

3 Литература

• Бахвалов Н. С., Жидков Н. П., Кобельков. Г. М. — 6-е изд. — М. : БИНОМ. Лаборатория знаний, 2008. — 636 с.
• Описание алгоритма на Вкипедии: Быстрое преобразование Фурье
• Материалы лекций по "Алгоритмам и структурам данных" Школы Анализа Данных Яндекса