Участник: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] рекурсивных вызовов.

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

Псевдокод:

fft(X, N)
begin
  (m, n) = factor(N)

  if (n == 1) // Предполагается, что n <= m
    direct_fft(X, N) // Через умножение на матрицу

  // 1
  transpose(X, m, n) //Транспонирует X как матрицу m x n, хранящуюся по строкам
   
  // 2
  do i = 0, n - 1
    fft(X(i * m : (i + 1) * m - 1), m)

  // 3
  do i = 0, n - 1
    do j = 0, m - 1
      X(i * m + j) *= eps(i * j, N)

  // 4
  transpose(X, n, m)

  // 5
  do j = 0, m - 1
    fft(X(j * n : (j + 1) * n - 1), n)

  // 6
  transpose(X, m, n)
end

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

Кроме того, шаги 3 и 4 можно переставлять произвольным образом (с очевидными изменениями в шаге 3).

Особо следует отметить, что если порядок элементов в выходном массиве не важен, шаг 6 можно пропустить.

Так, например, если преобразование Фурье используется для вычисления свёртки векторов, то конкретный порядок элементов в их образах Фурье совершенно не важен, лишь бы он был одинаков для обоих векторов. В этом случае в прямом преобразовании можно опустить шаг 6, а в обратном — 1.

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

Непосредственно из описания алгоритма получаем, что если сложность вычисления преобразования Фурье обозначить [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].

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

Рисунок 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].

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

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

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

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

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

Для описания этих случаев нам понадобится матрица

[math]\displaystyle{ Q = \begin{bmatrix} 1 & & & & & \\ & & & & & 1 \\ & & & & 1 & \\ & & & \cdots & & \\ & & 1 & & & \\ & 1 & & & & \end{bmatrix} }[/math]

1. Если [math]\displaystyle{ X }[/math] — вещественный, то [math]\displaystyle{ Y = Q\bar{Y} }[/math].
2. Если [math]\displaystyle{ X }[/math] — чисто мнимый, то [math]\displaystyle{ Y = -Q\bar{Y} }[/math].
3. Если [math]\displaystyle{ X = Q\bar{X} }[/math], то [math]\displaystyle{ Y }[/math] — вещественный.
4. Если [math]\displaystyle{ X = -Q\bar{X} }[/math], то [math]\displaystyle{ Y }[/math] — чисто мнимый.

В каждом случае, [math]\displaystyle{ X }[/math] и [math]\displaystyle{ Y }[/math] требуют для своего хранения в два раза меньше памяти, чем в общем случае.

Таким образом, при наличии дополнительной информации о входном векторе, возможно сократить объём памяти используемый для хранения входных и выходных данных.