Участник:Ivanov.kir.m/Быстрое дискретное преобразование Фурье

Авторы описания К.М.Иванов, А.C.Сапин Вклад каждого участника равномерный, поскольку каждый пункт вычитывался и исправлялся обоими участниками. При сильном приближении можно сказать,что больший вклад К.М.Иванов внес в заполнение теоретической части, а А.C.Сапин в заполнение практической части.

Быстрое преобразование Фурье (БПФ, FFT) — алгоритм быстрого вычисления дискретного преобразования Фурье (ДПФ). То есть, алгоритм вычисления за количество действий, меньшее чем [math]O(N^)[/math] , требуемых для прямого (по формуле) вычисления ДПФ. Иногда под БПФ понимается один из быстрых алгоритмов, называемый алгоритмом прореживания по частоте/времени или алгоритмом по основанию 2, имеющий сложность [math]O(N\log(N))[/math] . Cуществует несколько различных алгоритмов для вычисления ДПФ считающимся быстрым преобразование Фурье:

Алгоритм Кули-Тьюки [1]
Алгоритм Гуда-Томаса [2]
Алгоритм Бруна [3]
Алгоритм Блюштейна [4]

Содержание

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

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

Быстрое дискретное преобразование Фурье (БПФ) - популярный в современном мире математический метод преобразования вектора/матрицы комплексных/действительных чисел — сигнала — в вектор/матрицу комплексных чисел — спектр. С математической точки зрения преобразование сигнала [math]s[/math] (как функции по времени [math]t[/math] ) в спектр [math]\hat[/math] описывается формулой

В отличие от своего классического аналога, БПФ имеет сложность [math]O(N\log)[/math] , а не [math]O(N^2)[/math] , это достигается благодаря разделению исходного вектора на несколько частей так, чтобы, применив к этим частям БПФ и объединив результаты путем умножения на поворотные множители и дополнительного БПФ, можно было получить преобразование Фурье исходного вектора — принцип «разделяй и властвуй».

Область применения данного метода простирается от обработки аудио информации до собственно спектрального анализа. Поскольку в общем случае алгоритм БПФ применим к произвольным векторам чисел, выделяют отдельные подзадачи, для которых можно ввести дополнительную оптимизацию. Одной из таких подзадач является быстрое преобразование Фурье вектора действительных чисел, которому и посвящена данная статья. В качестве тестируемой реализации был выбран алгоритм Кули-Тьюки в библиотеке FFTW, написанной на языке C.

Замечание: Поскольку общие описание алгоритма БПФ весьма сложно, в некоторых пунктах данная статья будет ссылаться на простой случай — алгоритм Кули-Тьюки быстрого преобразование Фурье для векторов с размерностью равной степени двойки. Отличительной особенностью данного алгоритма является то, что он обходится без использования специфических приемов, использующихся именно для степеней четверки, восьмерки и т.п. Однако благодаря тому, что на вход данному алгоритму подается вектор чисто вещественных чисел, выходной вектор удовлетворяет эрмитовой избыточности (Hermitian redundancy) , т.е. [math]out[i][/math] является сопряженным с [math]out[n-i][/math] . Это обстоятельство позволяет достичь роста скорости и снижения затрат памяти примерно в 2 раза по сравнению с комплексным аналогом алгоритма.

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

Входные данные: вектор действительных чисел [math]x = (x_1,x_2. x_N)[/math] .

Выходные данные: вектор комплексных чисел [math]X = (X_1,X_2. X_N)[/math] .

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

В простейшем случае(для степеней двойки) алгоритм Кули-Тьюки рассчитывает ДПФ для четных [math](x_=x_0, x_2, \ldots, x_)[/math] и нечетных элементов отдельно [math](x_=x_1, x_3, \ldots, x_)[/math] . Затем результаты рассчетов объединяются для получения ДПФ всей последовательности:

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

1.2.1 Рекурсивное описание

Входной вектор [math]a = (a_0,a_2. a_)[/math] преобразуется в матрицу [math]A[/math] размера [math]n_1 \times n_2 [/math] , где [math]N=n_1 \cdot n_2[/math] и [math]n_1 \lt n_2[/math]

К каждой строке полученной матрицы применяется быстрое дискретное преобразование Фурье (БПФ) порядка [math]n_1[/math]
Каждый элемент полученный после применения БПФ умножается на поворотные множители [math]exp (2 \pi i(m-1)(j-1)/N)[/math] , где [math]m[/math] - номер строки, а [math]j[/math] - номер столбца
Полученная после шагов 1-2 матрица [math]A[/math] транспонируется
К каждой строке матрицы [math]A^T[/math] применяется БПФ порядка [math]n_2[/math]

Так как алгоритм реализует принцип «разделяй и властвуй», глубина рекурсии составляет [math]O (\log N)[/math] от числа входных элементов.

Замечание: Как правило все поворотные множители вычисляются заранее и хранятся в специальном массиве.

Замечание: В общем случае [math]N[/math] разбивается на [math]n_1 \cdot n_2[/math] так, чтобы [math]n_1[/math] и [math]n_2[/math] сами не являлись простыми числами и разлагались на множители, не являющиеся простыми числами.

Замечание: В случае простого [math]N[/math] оно дополняется(как правило единицей) до составного числа. В процессе слияния добавленные значения отбрасываются.

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

В случае размерности входа равной степени двойки, вычислительным ядром алгоритма является так называемая, "бабочка". В простейшем случае "бабочка" представляет из себя двухточечное преобразование. Рассмотрим этот случай:

На вход алгоритму подается двухэлементный вектор ‒ [math] v = (v[0], v[1]) [/math] . Тогда для вычисления будут происходить по следующим формулам:

[math]V[0] = W_2^0 v[0] + W_2^0 v[1] = v[0] + W_2^0 v[1] [/math]

[math]V[1] = W_2^0 v[0] + W_2^1 v[1] = v[0] + W_2^1 v[1] [/math]

Данный процесс удобно изобразить с помощью следующей схемы:

Для 4-х элеметного вектора [math]v=(v[0],v[1],v[2],v[3])[/math] , алгоритм строится похожим образом. Сначала создаются простейшие "бабочки", а потом их результаты соединяются с противоположеной "бабочкой":

[math]V[0]=v[0]+W_2^0 v[2]+W_4^0(v[1]+W_2^0 v[3])[/math]

[math]V[1]=v[0]-W_2^0 v[2]+W_4^1(v[1]-W_2^0 v[3])[/math]

[math]V[2]=v[0]+W_2^0 v[2]-W_4^0(v[1]+W_2^0 v[3])[/math]

[math]V[3]=v[0]-W_2^0 v[2]-W_4^1(v[1]-W_2^0 v[3])[/math]

Схема в таком случае будет выглядеть следующим образом:

Для случая, когда вход не является степенью двойки, "бабочки" будут "несимментричными", но в остальном вычисления будут проходить схожим образом.

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

Для исходного вектора [math]a = (a_1,a_2. a_N)[/math] размерности [math]N = n_1 \cdot n_2[/math] , [math]n_1 \lt n_2[/math] БПФ представляется как:

[math]n_2[/math] БПФ порядка [math]n_1[/math]
[math]n_1 \cdot n_2[/math] умножение комплексных чисел
[math]n_1[/math] БПФ порядка [math]n_2[/math]

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

Схема организации состоит в том, что на каждом шаге [math]i[/math] для выполнения "бабочки" все элементы разбиваются на [math]n_[/math] векторов по [math]n_[/math] элементов, причем [math]n_ \cdot n_ = n_i[/math] , где [math]n_i[/math] длина входного вектора на текущем шаге [math]i[/math] . В случае если такое разбиение невозможно, по причине того что [math]n_i[/math] простое число, исходный вектор на текущем шаге дополняется элементом. В зависимости от номера шага, разница координат для каждой пары элементов увеличивается соразмерно разбиению [math](n_,n_)[/math] . При этом результат суммы записывается в элемент с меньшим номером, а результат вычитания с последующим умножением - в элемент с большим.

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

Быстрое дискретное преобразование Фурье выполнимо за [math]O(N(n_1+\cdots+n_m))[/math] действий при [math]N=n_1n_2\cdots n_m[/math] (в простом случае, при [math]N=2^m[/math] необходимо [math]O(N\log_2(N))[/math] действий). Дискретное преобразование Фурье преобразует вектор [math]a = (a_0, \dots, a_)[/math] в вектор комплексных чисел [math]b = (b_0, \dots, b_)[/math] , такой, что [math]b_i=\sum_^a_j\varepsilon^[/math] , где [math]\varepsilon^n=1[/math] и [math]\varepsilon^k\neq 1[/math] при [math]0\lt k\lt N[/math] .

Основной шаг алгоритма состоит в сведении задачи для [math]N[/math] чисел к задаче для [math]n_1=N/n_2[/math] числам, где [math]n_2[/math] — делитель [math]N[/math] .

Пусть мы уже умеем решать задачу для [math]N/n_2[/math] чисел. Применим преобразование Фурье к векторам [math]a_i,a_, \dots, a_[/math] для [math]i=0,1,\dots,n_2-1[/math] . Покажем теперь, что за [math]O(Np)[/math] действий можно решить исходную задачу. Заметим, что [math]b_i=\sum_^ \varepsilon^ \left(\sum_^a_\varepsilon^\right)[/math] . Выражения в скобках нам уже известны — это [math](i\mod p)[/math] -тое число после преобразования Фурье [math]j[/math] -го вектора. Таким образом, для вычисления каждого [math]b_i[/math] нужно [math]O(n_2)[/math] действий, а для вычисления всех [math]b_i[/math] всего [math]O(Nn_2)[/math] действий, что и требовалось.

Таким образом задача сведена к вычислению свёртки, а это можно сделать с помощью трёх преобразований Фурье для [math]2^k[/math] элементов. Выполняем прямое преобразование Фурье для [math](\bar, \cdots, \bar_)[/math] и [math](c_1,\cdots,c_)[/math] , перемножаем поэлементно результаты и выполняем обратное преобразование Фурье.

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

Как видно из рисунка, размер этого графа линейно-логарифмический, а формулы дуг имеют экспоненциальные компоненты.

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

Если считать только главные члены выражений, то простой алгоритм Кули-Тьюки имеет критический путь, состоящий из [math]log(N)[/math] операций комплексного сложения/вычитания и [math]log(N)[/math] операций комплексного умножения (основание [math]log[/math] целиком зависит от выбранного на каждом шаге алгоритма разбиения). Здесь стоит отметить, что все рекурсивные вызовы БПФ на каждом шаге выполняются независимо, это приводит к тому, что их можно распределить по доступным вычислительным узлам. Аналогичная ситуация обстоит и с умножением на поворотные множетели, что позволяет независимо присоединять его к любому шагу. Таким образом алгоритм БПФ имеет высокий ресурс параллелизма, а алгоритм БПФ в общем случае может быть отнесён к логарифмическому классу по параллельной сложности.

Кроме того, параллелизм БПФ можно увеличить, если брать [math]n_1[/math] и [math]n_2[/math] по возможности близкими к кратным доступному количеству вычислительных узлов (а не брать просто один из множителей близким к нулю, как в традиционной реализации), это позволит равномерно распределить работу между вычислительными узлами, что сильно повысит эффективность алгоритма благодаря лишь одной внутренней пересылке данных.

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

Входные данные: вектор действительных/комплексных чисел [math]a = (a_1,a_2. a_N)[/math] .

Выходные данные: вектор комплексных чисел [math]b = (b_1,b_2. b_)[/math] .

Есть несколько ситуаций, при которых количество потребляемой памяти для хранения входного вектора [math]a[/math] и выходного вектора [math]b[/math] тможет быть снижено вдвое:

Если [math]a[/math] — вектор действительных чисел, тогда [math]b = E\bar[/math] ;

Если [math]a[/math] — вектор чисто-мнимых чисел, тогда [math]b = -E\bar[/math] ;

Если [math]a = E\bar[/math] , тогда [math]b[/math] — вектор действительных чисел;

Если [math]a = -E\bar[/math] , то [math]b[/math] — вектор чисто-мнимых чисел.

Где матрица [math]E[/math] :

[math] E = \begin 1 & 0 & 0 & \cdots & 0 & 0 \\ 0 & 0 & 0 & \cdots & 0 & 1 \\ 0 & 0 & 0 & \cdots & 1 & 0 \\ 0 & 0 & 0 & \ddots & 0 & 0 \\ 0 & 0 & 1 & \cdots & 0 & 0 \\ 0 & 1 & 0 & \cdots & 0 & 0 \end [/math]

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

Таким образом в общем у алгоритма БПФ в случае неограниченных ресурсов:

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

параллельная сложность - логарифмическая.

При этом алгоритм полностью детерминирован.

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

В качестве тестируемой реализации была выбрана библиотека FFTW, а именно функция, выполняющая БПФ одномерного вектора действительных чисел.

2.1 Особенности последовательной реализации алгоритма

Простейший вариант алгоритма для степеней двойки на языке C:

2.2 Локальность данных и вычислений

2.3 Возможные способы и особенности параллельной реализации алгоритма

2.4 Масштабируемость алгоритма и его реализации

Тестирование проводилось на вычислительном комплексе IBM Blue Gene/P
2.4.1 Описание платформы IBM Blue Gene/P
IBM Blue Gene/P — массивно-параллельная вычислительная система, которая состоит из двух стоек, включающих 8192 процессорных ядер (2 x 1024 четырехъядерных вычислительных узлов), с пиковой производительностью 27,9 терафлопс (27,8528 триллионов операций с плавающей точкой в секунду).

две стойки с вычислительными узлами и узлами ввода-вывода

1024 четырехъядерных вычислительных узла в каждой из стоек

16 узлов ввода-вывода в стойке (в текущей конфигурации активны 8, т.е. одна I/O-карта на 128 вычислительных узлов)

выделенные коммуникационные сети для межпроцессорных обменов и глобальных операций

программирование с использованием MPI, OpenMP/pthreads, POSIX I/O

высокая энергоэффективность: ∼ 372 MFlops/W (см. список Green500)

система воздушного охлаждения

Стойка (rack, cabinet) состоит из двух midplane'ов. В midplane входит 16 node-карт (compute node card), на каждой из которых установлено 32 вычислительных узла (compute card). Midplane, 8 x 8 x 8 = 512 вычислительных узлов, — минимальный раздел, на котором становится доступна топология трехмерного тора; для разделов меньших размеров используется топология трехмерной решетки. Node-карта может содержать до двух узлов ввода-вывода (I/O card). Вычислительный узел включает в себя четырехъядерный процессор, 2 ГБ общей памяти и сетевые интерфейсы.

модель: PowerPC 450

рабочая частота: 850 MHz

адресация: 32-битная

кэш инструкций 1-го уровня (L1 instruction): 32 KB

кэш данных 1-го уровня (L1 data): 32 KB

кэш предвыборки (L2 prefetch): 14 потоков предварительной выборки (stream prefetching): 14 x 256 байтов

два блока 64-битной арифметики с плавающей точкой (Floating Point Unit, FPU), каждый из которых может выдавать за один такт результат совмещенной операции умножения-сложения (Fused Multiply-Add, FMA)

пиковая производительность: 2 FPU x 2 FMA x 850 MHz = 3,4 GFlop/sec per core

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

четыре микропроцессорных ядра PowerPC 450 (4-way SMP)

пиковая производительность: 4 cores x 3,4 GFlop/sec per core = 13,6 GFlop/sec

пропускная способность памяти: 13,6 GB/sec

2 ГБ общей памяти

2 x 4 МБ кэш-памяти 2-го уровня (в документации по BG/P носит название L3)

легковесное ядро (compute node kernel, CNK), представляющее собой Linux-подобную операционную систему, поддерживающую значительное подмножество Linux-совместимых системных вызовов

асинхронные операции межпроцессорных обменов (выполняются параллельно с вычислениями)

операции ввода-вывода перенаправляются I/O-картам через сеть коллективных операций

2.4.2 Оценка масштабируемости алгоритма
Масштабируемость алгоритма проверялась с помощью тестовой программы расположенной по ссылке. Для эксперимента использовались:

Компилятор mpicc

Библиотека OpenMPI версии 1.68

Библиотека FFTW 3.3.5

Алгоритм тестировался на входных векторах размера [math]2^[/math] до [math]2^[/math] . При этом использовались от 8 до 128 процессоров с шагом степени двойки.
📎📎📎📎📎📎📎📎📎📎