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

Файл:Warning sign font awesome.svg

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

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

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

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

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

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

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

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

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

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

Многомерное преобразование Фурье

Естественным обобщением описанного одномерного преобразования Фурье является многомерный вариант:

[math]\displaystyle{ X(l_1,\cdots,l_d) = \sum_{k_1 = 0}^{N_1} \cdots \sum_{k_d = 0}^{N_d} Y(k_1,\cdots,k_d) \prod_{j=1}^{d} \epsilon_{N_j}^{l_j k_j} }[/math]

Многомерное преобразование Фурье эквивалентно последовательному выполнению одномерных преобразований по каждому измерению; при использовании для этого рекурсивного алгоритма, общая сложность вычислений также составляет [math]\displaystyle{ O(N\log N) }[/math], [math]\displaystyle{ N = \prod_{j=1}^d N_j }[/math].

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

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

В этом случае

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

[math]\displaystyle{ h(mn) = 1 + h(m) + h(n) }[/math]

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

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

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

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

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

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

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