Участник:Timokhinivan/Быстрое преобразование Фурье

Данная страница в настоящее время активно редактируется участником Timokhinivan (обсуждение). Пожалуйста, не вносите в неё никаких изменений до тех пор, пока не исчезнет это объявление. В противном случае могут возникнуть конфликты редактирования.

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

1 Свойства и структура алгоритмов

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

Преобразование Фурье переводит сигнал в его спектр и обратно, и широко применяется в самых различных областях вычислительной математики — собственно спектральный анализ, сжатие данных, вычисление свёрток и т.д. В связи с этим, особый интерес представляют быстры алгоритмы для вычисления этого преобразование. На сегодняшний день, наилучшие из них имеют сложность $O(n \log n)$.

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

Преобразование Фурье задаётся следующей формулой:

$Y_k = \sum_{l = 0}^{N-1} X_l \epsilon^{kl}_N,\qquad k = \overline{0,N-1}$

где $\epsilon_N = e^{\frac{2i\pi}{N}}$ — первообразный корень из 1 степени $N$.

Таким образом, преобразование Фурье является линейным и задаётся матрицей $F = \left\{ \epsilon^{kl}_N \right\}_{k\,l = 0}^{N-1}$.

1.2.2 Сведение к двумерному преобразованию

Если $N$ — составное, т.е. $N = mn$, то возможно существенно сократить вычислительные расходы за счёт сведение к двумерному преобразованию Фурье. А именно, положим $X(l_1, l_2) = X_{l_1 n + l_2}$, $Y(k_1, k_2) = Y_{k_1 m + k_2}$, где $l_1, k_2 = \overline{0, m-1}$, $l_2, k_1 = \overline{0, n-1}$.

В этом случае можно показать, что

$\begin{align} Y(k_1, k_2) &=& \sum_{l_2=0}^{n-1} (\epsilon^{k_2 l_2}_N \hat{X}(k_2, l_2)) \epsilon^{k_1 l_2}_n \\ \hat{X}(k_2, l_2) &=& \sum_{l_1=0}^{m-1} X(l_1, l_2) \epsilon^{k_2 l_1}_m \end{align}$

Таким образом, получаем алгоритм вычисления преобразования:

1. Записываем исходный вектор в матрицу $m\times n$ по строкам.
2. Применяем к каждому столбцу преобразование Фурье.
3. Умножаем элемент в позиции $(i,j)$ на $\epsilon^{ij}_N$.
4. Применяем к каждой строке преобразование Фурье.
5. Результат записан в получившейся матрице по столбцам.

Данный алгоритм уже даёт существенный выигрыш по сравнению с обычным умножением матрицы на вектор: $O(m^2 n + n^2 m)$ против $O(m^2 n^2)$. Однако наилучших результатов можно добиться, если применять этот алгоритм рекурсивно на этапах 2 и 4.

Так, если $N = \prod_{i=1}^{K} p_i$, то целесообразно на каждом уровне рекурсии «отщеплять» одно $p_i$. В этом случае задача сводится к вычислению $\prod_{i\neq j} p_i$ преобразований Фурье порядка $p_j$ для всех $j$. При небольших $p_i$ (например, 2), это можно проделывать «в лоб».

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

На каждом уровне рекурсии наиболее дорогостоящими этапами являются рекурсивные вызовы преобразования Фурье: они требуют в сумме $O(mn (\log m + \log n))$ операций против $O(mn)$ для умножения на поправочные коэффициенты (шаг 3).

То же верно и для алгоритма в целом в случае $N = \prod_{i=1}^{K} p_i$; а именно, суммарно наибольшие вычислительные затраты связаны с вычислением преобразований Фурье порядка $p_i$.

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

Фактически, макроструктура алгоритма уже описана в разделе математического описания. На каждом уровне рекурсии, алгоритм состоит из

1. $n$ рекурсивных вызовов алгоритма.
2. Умножение всех элементов рабочего вектора на поправочные коэффициенты.
3. Ещё $m$ рекурсивных вызовов.

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

Непосредственно из описания алгоритма получаем, что если сложность вычисления преобразования Фурье обозначить $f_N$, то $f_{mn} = n f_{m} + m f_{n} + mn$.

Если $N = \prod_{i=1}^K p_i$, и на каждом шаге «отщеплять» от него одно $p_i$, а преобразование Фурье для $p_i$ вычислять «в лоб», то общее количество операций составит $O\left(N\sum_{i=1}^K p_i\right)$.

В частности, если $p_i \leq P$, то $K \leq \log_P N$ и для сложности получаем $O(NP\log_P N)$. Полагая $P$ фиксированным, получаем заявленную сложность $O(N \log N)$.

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

Рисунок 1. Информационный граф алгоритма для $N = 15$. Исходные данные обозначены фиолетовым, результат — красным. Преобразования Фурье для $m = 5$ и $n = 3$ представлены как «чёрные ящики». Умножение на дополнительные коэффициенты представлено оранжевыми узлами.

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

$h(mn) = 1 + h(m) + h(n)$

Для «элементарных» БПФ высота ЯПФ будет та же, что и для умножения матрицы на вектор.

Таким образом, с учётом разложения $N$, $h(N) = O\left(\sum_{i=1}^K p_i\right) = O(P \log_P N) = O(\log N)$.

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

Поскольку все преобразования Фурье на шагах 2 и 4 алгоритма совершенно независимы, кажется естественным распределить их по доступным вычислительным узлам. Шаг 3 при этом и вовсе выполняется независимо на каждом элементе рабочего вектора, и может быть беспрепятственно присоединён к любому из них.

При этом, в отличие от традиционной реализации типа Кули-Тьюки, в которой на каждом этапе один из множителей берётся малым, при параллельной реализации целесообразно взять и $m$ и $n$ по возможности близкими к кратным доступному количеству вычислительных узлов, поскольку в этом случае возможно равномерно распределить работу между ними и реализовать весь алгоритм всего с одной внутренней пересылкой.

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

В общем случае на входе и на выходе имеются комплексные векторы порядка $N$.