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

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

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

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

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

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

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

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

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

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

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

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

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

В этом случае

[math] \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} [/math]

Это даёт следующий алгоритм вычисления преобразования:

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

[math] h(mn) = 1 + h(m) + h(n) [/math]

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

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

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

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

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

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

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