Уровень алгоритма

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

Материал из Алговики
Перейти к навигации Перейти к поиску
 
(не показано 114 промежуточных версий 5 участников)
Строка 1: Строка 1:
 +
{{Assignment|IgorS|Kronberg}}
 +
 
{{algorithm
 
{{algorithm
| name              = Алгоритм Кули-Тьюки одномерного преобразования Фурье для действительных чисел
+
| name              = Алгоритм одномерного преобразования Фурье для действительных чисел
 
| serial_complexity = <math>O (N \log N)</math>
 
| serial_complexity = <math>O (N \log N)</math>
 
| pf_height        = <math>O (\log N)</math>
 
| pf_height        = <math>O (\log N)</math>
Строка 7: Строка 9:
 
| output_data      = <math>\lfloor N/2 \rfloor+1</math> комплексных чисел
 
| output_data      = <math>\lfloor N/2 \rfloor+1</math> комплексных чисел
 
}}
 
}}
 +
Авторы описания [[Участник:Ivanov.kir.m | К.М.Иванов]], [[Участник:Pumpkin | А.C.Сапин]] Вклад каждого участника равномерный, поскольку каждый пункт вычитывался и исправлялся обоими участниками. При сильном приближении можно сказать,что больший вклад [[Участник:Ivanov.kir.m | К.М.Иванов]] внес в заполнение теоретической части, а [[Участник:Pumpkin | А.C.Сапин]] в заполнение практической части.
 +
 
'''Быстрое преобразование Фурье (БПФ, FFT)''' — алгоритм быстрого  вычисления дискретного преобразования Фурье (ДПФ). То есть, алгоритм вычисления за количество действий, меньшее чем <math>O(N^{2})</math>, требуемых для прямого (по формуле) вычисления ДПФ. Иногда под БПФ понимается один из быстрых алгоритмов, называемый алгоритмом прореживания по частоте/времени или алгоритмом по основанию 2, имеющий сложность <math>O(N\log(N))</math>. Cуществует несколько различных алгоритмов для вычисления ДПФ считающимся быстрым преобразование Фурье:
 
'''Быстрое преобразование Фурье (БПФ, FFT)''' — алгоритм быстрого  вычисления дискретного преобразования Фурье (ДПФ). То есть, алгоритм вычисления за количество действий, меньшее чем <math>O(N^{2})</math>, требуемых для прямого (по формуле) вычисления ДПФ. Иногда под БПФ понимается один из быстрых алгоритмов, называемый алгоритмом прореживания по частоте/времени или алгоритмом по основанию 2, имеющий сложность <math>O(N\log(N))</math>. Cуществует несколько различных алгоритмов для вычисления ДПФ считающимся быстрым преобразование Фурье:
 
* Алгоритм Кули-Тьюки [https://en.wikipedia.org/wiki/Cooley%E2%80%93Tukey_FFT_algorithm]
 
* Алгоритм Кули-Тьюки [https://en.wikipedia.org/wiki/Cooley%E2%80%93Tukey_FFT_algorithm]
Строка 16: Строка 20:
  
 
== Общее описание алгоритма ==
 
== Общее описание алгоритма ==
Одним из вариантов '''быстрого преобразования Фурье''' для вектора '''действительных чисел''' с размерностью равной степени двойки является '''алгоритм Кули-Тьюки'''. Отличительной особенностью данного алгоритма является то, что он обходится без использования специфических приемов, использующихся именно для степеней четверки, восьмерки и т.п. Однако благодаря тому, что на вход данному алгоритму подается вектор чисто вещественных чисел, выходной вектор удовлетворяет '''эрмитовой избыточности''' (''Hermitian redundancy'') , т.е.  <math>out[i]</math> является сопряженным с <math>out[n-i]</math>. Это обстоятельство позволяет достичь роста скорости и снижения затрат памяти примерно в 2 раза по сравнению с ''комплексным'' аналогом алгоритма.
+
'''Быстрое дискретное преобразование Фурье''' (''БПФ'') - популярный в современном мире математический метод преобразования вектора/матрицы комплексных/действительных чисел — ''сигнала'' — в вектор/матрицу комплексных чисел — ''спектр''. С математической точки зрения преобразование сигнала <math>s</math> (как функции по времени <math>t</math>) в спектр <math>\hat{s}</math> описывается формулой
 +
:<math>
 +
\hat{s}(\omega)=\frac{1}{\sqrt{2\pi}}\int\limits_{-\infty}^{\infty}s(t)e^{-it\omega}\,dt
 +
</math>
 +
 
 +
В отличие от своего классического аналога, БПФ имеет сложность <math>O(N\log{N})</math>, а не <math>O(N^2)</math>, это достигается благодаря разделению исходного вектора на несколько частей так, чтобы, применив к этим частям БПФ и объединив результаты путем умножения на поворотные множители и дополнительного БПФ, можно было получить преобразование Фурье исходного вектора — принцип «разделяй и властвуй».
 +
 
 +
Область применения данного метода простирается от обработки аудио информации до собственно спектрального анализа.
 +
Поскольку в общем случае алгоритм БПФ применим к произвольным векторам чисел, выделяют отдельные подзадачи, для которых можно ввести дополнительную оптимизацию. Одной из таких подзадач является быстрое преобразование Фурье вектора действительных чисел, которому и посвящена данная статья. В качестве тестируемой реализации был выбран алгоритм Кули-Тьюки в библиотеке FFTW, написанной на языке C.
 +
 
 +
 
 +
'''Замечание''': Поскольку общие описание алгоритма БПФ весьма сложно, в некоторых пунктах данная статья будет ссылаться на простой случай — алгоритм Кули-Тьюки быстрого преобразование Фурье для векторов ''с размерностью равной степени двойки''. Отличительной особенностью данного алгоритма является то, что он обходится без использования специфических приемов, использующихся именно для степеней четверки, восьмерки и т.п. Однако благодаря тому, что на вход данному алгоритму подается вектор чисто вещественных чисел, выходной вектор удовлетворяет ''эрмитовой избыточности'' (''Hermitian redundancy'') , т.е.  <math>out[i]</math> является сопряженным с <math>out[n-i]</math>. Это обстоятельство позволяет достичь роста скорости и снижения затрат памяти примерно в 2 раза по сравнению с ''комплексным'' аналогом алгоритма.
  
 
== Математическое описание алгоритма ==
 
== Математическое описание алгоритма ==
  
Входные данные: вектор действительных чисел <math>a = (a_1,a_2,...,a_N)</math>.
+
Входные данные: вектор действительных чисел <math>x = (x_1,x_2,...,x_N)</math>.
  
Выходные данные: вектор комплексных чисел <math>b = (b_1,b_2,...,b_{\lfloor N/2 \rfloor+1})</math>.
+
Выходные данные: вектор комплексных чисел <math>X = (X_1,X_2,...,X_N)</math>.
  
'''''Замечание:''''' В простейшем случае '''алгоритм Кули-Тьюки''' применяется к векторам '''размерности степени двойки''', поэтому на практике вектора иной размерности часто '''дополнять''' до ближайшей степени двойки. Такой подход делает '''алгоритм Кули-Тьюки''' не самым эффективным алгоритмом БПФ, поскольку дополнение до степени двойки может сильно усложнить задачу.
+
Дискретное преобразование Фурье задается следующей формулой:
  
== Рекурсивное описание ==
+
<math>X_k = \sum_{n=0}^{N-1} x_n e^{-\frac{2\pi i}{N} nk}, k = \overline{0,N-1}</math>
 +
 
 +
В простейшем случае(для степеней двойки) алгоритм Кули-Тьюки рассчитывает ДПФ для четных <math>(x_{2m}=x_0, x_2, \ldots, x_{N-2})</math> и нечетных элементов отдельно <math>(x_{2m+1}=x_1, x_3, \ldots, x_{N-1})</math>. Затем результаты рассчетов объединяются для получения ДПФ всей последовательности:
 +
 
 +
<math>
 +
  \begin{matrix} X_k & =
 +
& \sum \limits_{m=0}^{N/2-1} x_{2m}e^{-\frac{2\pi i}{N} (2m)k}  +  \sum \limits_{m=0}^{N/2-1} x_{2m+1} e^{-\frac{2\pi i}{N} (2m+1)k}
 +
  \end{matrix}
 +
</math>
 +
 
 +
Таким образом, получается рекурсивный алгоритм, работающий по принципу [https://ru.wikipedia.org/wiki/%D0%A0%D0%B0%D0%B7%D0%B4%D0%B5%D0%BB%D1%8F%D0%B9_%D0%B8_%D0%B2%D0%BB%D0%B0%D1%81%D1%82%D0%B2%D1%83%D0%B9_(%D0%B8%D0%BD%D1%84%D0%BE%D1%80%D0%BC%D0%B0%D1%82%D0%B8%D0%BA%D0%B0) разделяй и влавствуй].
 +
 
 +
=== Рекурсивное описание ===
 
Алгоритм:
 
Алгоритм:
 
# Входной вектор <math>a = (a_0,a_2,...,a_{N-1})</math> преобразуется в матрицу <math>A</math> размера <math>n_1 \times n_2 </math>, где <math>N=n_1 \cdot n_2</math> и <math>n_1 < n_2</math>
 
# Входной вектор <math>a = (a_0,a_2,...,a_{N-1})</math> преобразуется в матрицу <math>A</math> размера <math>n_1 \times n_2 </math>, где <math>N=n_1 \cdot n_2</math> и <math>n_1 < n_2</math>
Строка 36: Строка 63:
 
  \end{pmatrix} </math>
 
  \end{pmatrix} </math>
 
# К каждой строке полученной матрицы применяется быстрое дискретное преобразование Фурье (БПФ) порядка <math>n_1</math>
 
# К каждой строке полученной матрицы применяется быстрое дискретное преобразование Фурье (БПФ) порядка <math>n_1</math>
# Каждый элемент полученный после применения БПФ умножается на '''''поворотные множители''''' <math>exp (2 \pi i(m-1)(j-1)/N)</math>, где <math>m</math> - номер строки, а <math>j</math> - номер столбца
+
# Каждый элемент полученный после применения БПФ умножается на ''поворотные множители'' <math>exp (2 \pi i(m-1)(j-1)/N)</math>, где <math>m</math> - номер строки, а <math>j</math> - номер столбца
 
# Полученная после шагов 1-2 матрица <math>A</math> транспонируется
 
# Полученная после шагов 1-2 матрица <math>A</math> транспонируется
 
# К каждой строке матрицы <math>A^T</math> применяется БПФ порядка <math>n_2</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> оно дополняется(как правило единицей) до составного числа. В процессе слияния добавленные значения отбрасываются.
  
 
== Вычислительное ядро алгоритма ==
 
== Вычислительное ядро алгоритма ==
В случае размерности входа равной степени двойки, вычислительным ядром алгоритма является, так называемая, "''бабочка''". В простейшем случае "''бабочка''" представляет из себя двухточечное преобразование. Рассмотрим этот случай:
+
В случае размерности входа равной степени двойки, вычислительным ядром алгоритма является так называемая, "''бабочка''". В простейшем случае "''бабочка''" представляет из себя двухточечное преобразование. Рассмотрим этот случай:
  
 
На вход алгоритму подается двухэлементный вектор ‒ <math> v = (v[0], v[1]) </math>. Тогда для вычисления будут происходить по следующим формулам:
 
На вход алгоритму подается двухэлементный вектор ‒ <math> v = (v[0], v[1]) </math>. Тогда для вычисления будут происходить по следующим формулам:
Строка 73: Строка 106:
 
== Макроструктура алгоритма ==
 
== Макроструктура алгоритма ==
 
Для исходного вектора <math>a = (a_1,a_2,...,a_N)</math> размерности <math>N = n_1 \cdot n_2</math>, <math>n_1 < n_2</math> БПФ представляется как:
 
Для исходного вектора <math>a = (a_1,a_2,...,a_N)</math> размерности <math>N = n_1 \cdot n_2</math>, <math>n_1 < n_2</math> БПФ представляется как:
 
+
#<math>n_2</math> БПФ порядка <math>n_1</math>  
<math>n_2</math> БПФ порядка <math>n_1</math> <math>\Rightarrow</math> <math>n_1 \cdot n_2</math> умножение комплексных чисел <math>\Rightarrow</math> <math>n_1</math> БПФ порядка <math>n_2</math>
+
#<math>n_1 \cdot n_2</math> умножение комплексных чисел  
 +
#<math>n_1</math> БПФ порядка <math>n_2</math>
  
 
== Схема реализации последовательного алгоритма ==
 
== Схема реализации последовательного алгоритма ==
Не рекурсивная схема организации состоит в том, что на каждом шаге <math>i</math> для выполнения "''бабочки''" все элементы разбиваются на <math>n_{i_1}</math> векторов по <math>n_{i_2}</math> элементов, причем <math>n_{i_1} \cdot n_{i_2} = n_i</math>, где <math>n_i</math> длина входного вектора на текущем шаге <math>i</math>. В случае если такое разбиение невозможно, по причине того что <math>n_i</math> простое число, исходный вектор на текущем шаге дополняется элементом. В зависимости от номера шага, разница координат для каждой пары элементов увеличивается соразмерно разбиению  <math>(n_{i_1},n_{i_2})</math>.
+
Схема организации состоит в том, что на каждом шаге <math>i</math> для выполнения "''бабочки''" все элементы разбиваются на <math>n_{i_1}</math> векторов по <math>n_{i_2}</math> элементов, причем <math>n_{i_1} \cdot n_{i_2} = n_i</math>, где <math>n_i</math> длина входного вектора на текущем шаге <math>i</math>. В случае если такое разбиение невозможно, по причине того что <math>n_i</math> простое число, исходный вектор на текущем шаге дополняется элементом. В зависимости от номера шага, разница координат для каждой пары элементов увеличивается соразмерно разбиению  <math>(n_{i_1},n_{i_2})</math>.
 
При этом результат суммы записывается в элемент с меньшим номером, а результат вычитания с последующим умножением - в элемент с большим.
 
При этом результат суммы записывается в элемент с меньшим номером, а результат вычитания с последующим умножением - в элемент с большим.
  
 
== Последовательная сложность алгоритма ==
 
== Последовательная сложность алгоритма ==
Быстрое дискретное преобразование Фурье выполнимо за <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_{N-1})</math> в вектор чисел <math>b = (b_0, \dots, b_{N-1})</math>, такой, что <math>b_i=\sum_{j=0}^{N-1}a_j\varepsilon^{ij}</math>, где <math>\varepsilon^n=1</math> и <math>\varepsilon^k\neq 1</math> при <math>0<k<N</math>.
+
Быстрое дискретное преобразование Фурье выполнимо за <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_{N-1})</math> в вектор комплексных чисел <math>b = (b_0, \dots, b_{N-1})</math>, такой, что <math>b_i=\sum_{j=0}^{N-1}a_j\varepsilon^{ij}</math>, где <math>\varepsilon^n=1</math> и <math>\varepsilon^k\neq 1</math> при <math>0<k<N</math>.
  
 
Основной шаг алгоритма состоит в сведении задачи для <math>N</math> чисел к задаче для <math>n_1=N/n_2</math> числам, где <math>n_2</math> — делитель <math>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_{n_2+i}, \dots, a_{n_2(n_1-1)+i}</math> для <math>i=0,1,\dots,n_2-1</math>. Покажем теперь, что за <math>O(Np)</math> действий можно решить исходную задачу. Заметим, что <math>b_i=\sum_{j=0}^{n_2-1} \varepsilon^{ij} (\sum_{k=0}^{n_1-1}a_{kn_2+j}\varepsilon^{kin_2})</math>.
 
Выражения в скобках нам уже известны — это <math>(i \pmod p)</math>-тое число после преобразования Фурье <math>j</math>-го вектора. Таким образом, для вычисления каждого <math>b_i</math> нужно <math>O(n_2)</math> действий, а для вычисления всех <math>b_i</math> всего <math>O(Nn_2)</math> действий, что и требовалось.
 
  
В '''общем''' случае. Пусть <math>4N>2^k\ge2N</math>. Заметим, что тогда <math>b_i=\varepsilon^{-i^2/2}\sum_{j=0}^{N-1}\varepsilon^{(i+j)^2/2}\varepsilon^{-j^2/2}a_j</math>. Обозначим <math>\bar{a}_i=\varepsilon^{-i^2/2}a_i</math><math>\bar{b}_i=\varepsilon^{i^2/2}b_i</math>,   <math>c_i=\varepsilon^{(2N-2-i)^2/2}</math>. Тогда <math>\bar{b}_i=\sum_{j=0}^{2N-2-i}\bar{a}_jc_{2N-2-i-j}</math>, если положить <math>\bar{a}_i=0</math> при <math>i\ge N</math>.
+
Пусть мы уже умеем решать задачу для <math>N/n_2</math> чисел. Применим преобразование Фурье к векторам <math>a_i,a_{n_2+i}, \dots, a_{n_2(n_1-1)+i}</math> для <math>i=0,1,\dots,n_2-1</math>. Покажем теперь, что за <math>O(Np)</math> действий можно решить исходную задачу. Заметим, что <math>b_i=\sum_{j=0}^{n_2-1} \varepsilon^{ij} \left(\sum_{k=0}^{n_1-1}a_{kn_2+j}\varepsilon^{kin_2}\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{a_0}, \cdots, \bar{a}_{2^k-1})</math> и <math>(c_1,\cdots,c_{2^k-1})</math>, перемножаем поэлементно результаты и выполняем обратное преобразование Фурье.
+
В общем случае:
  
Вычисления всех <math>\bar{a}_i</math> и <math>c_i</math> требуют <math>O(N)</math> действий, три преобразования Фурье требуют <math>O(N\log(N))</math> действий, перемножение результатов преобразований Фурье требует <math>O(N)</math> действий, вычисление всех <math>b_i</math> зная значения свертки требует <math>O(N)</math> действий. '''Итого для дискретного преобразования Фурье требуется''' <math>O(N\log(N))</math> '''действий для любого''' <math>N</math>.
+
Пусть <math>4N>2^k\ge2N</math>. Заметим, что тогда <math>b_i=\varepsilon^{-i^2/2}\sum_{j=0}^{N-1}\varepsilon^{(i+j)^2/2}\varepsilon^{-j^2/2}a_j</math>. Обозначим <math>\bar{a}_i=\varepsilon^{-i^2/2}a_i</math>,  <math>\bar{b}_i=\varepsilon^{i^2/2}b_i</math>,  <math>c_i=\varepsilon^{(2N-2-i)^2/2}</math>. Тогда <math>\bar{b}_i=\sum_{j=0}^{2N-2-i}\bar{a}_jc_{2N-2-i-j}</math>, если положить <math>\bar{a}_i=0</math> при <math>i\ge N</math>.
 +
 
 +
Таким образом задача сведена к вычислению ''свёртки'', а это можно сделать с помощью трёх преобразований Фурье для <math>2^k</math> элементов. Выполняем прямое преобразование Фурье для <math>(\bar{a_0}, \cdots, \bar{a}_{2^k-1})</math> и <math>(c_1,\cdots,c_{2^k-1})</math>, перемножаем поэлементно результаты и выполняем обратное преобразование Фурье.
 +
 
 +
Вычисления всех <math>\bar{a}_i</math> и <math>c_i</math> требуют <math>O(N)</math> действий, три преобразования Фурье требуют <math>O(N\log(N))</math> действий, перемножение результатов преобразований Фурье требует <math>O(N)</math> действий, вычисление всех <math>b_i</math> зная значения свертки требует <math>O(N)</math> действий. Итого для дискретного преобразования Фурье требуется <math>O(N\log(N))</math> действий для любого <math>N</math>.
 +
 
 +
== Информационный граф ==
 +
[[Файл:InfoGraph Ivanov Sapin.png |center|thumb|800px| Информационный граф алгоритма Кули-Тьюки для <math>N = 21</math>. Исходные данные обозначены зеленым, результат — синим.  Преобразования Фурье для <math>n_1 = 7</math> и <math>n_2 = 3 \text{  } (N = n_1 \cdot n_2)</math> представлены как «чёрные ящики». Умножение на поворотные коэффициенты (множители) представлено коричневыми ромбами.]]
 +
Как видно из рисунка, размер этого графа линейно-логарифмический, а формулы дуг имеют экспоненциальные компоненты.
  
 
== Ресурс параллелизма алгоритма ==
 
== Ресурс параллелизма алгоритма ==
Если считать только главные члены выражений, то простой алгоритм Кули-Тьюки имеет критический путь, состоящий из <math>log(N)</math> операций комплексного сложения/вычитания и <math>log(N)</math> операций комплексного умножения (основание <math>log</math> целиком зависит от выбранного на каждом шаге алгоритма разбиения). Таким образом, простой алгоритм БПФ в общем случае может быть отнесён к ''логарифмическому'' классу по параллельной сложности.
+
Если считать только главные члены выражений, то простой алгоритм Кули-Тьюки имеет критический путь, состоящий из <math>log(N)</math> операций комплексного сложения/вычитания и <math>log(N)</math> операций комплексного умножения (основание <math>log</math> целиком зависит от выбранного на каждом шаге алгоритма разбиения). Здесь стоит отметить, что все рекурсивные вызовы БПФ на каждом шаге выполняются независимо, это приводит к тому, что их можно распределить по доступным вычислительным узлам. Аналогичная ситуация обстоит и с умножением на поворотные множетели, что позволяет независимо присоединять его к любому шагу. Таким образом алгоритм БПФ имеет высокий ресурс параллелизма, а алгоритм БПФ в общем случае может быть отнесён к ''логарифмическому'' классу по параллельной сложности.
 +
 
 +
Кроме того, параллелизм БПФ можно увеличить, если брать <math>n_1</math> и <math>n_2</math> по возможности близкими к кратным доступному количеству вычислительных узлов (а не брать просто один из множителей близким к нулю, как в традиционной реализации), это позволит равномерно распределить работу между вычислительными узлами, что сильно повысит эффективность алгоритма благодаря лишь одной внутренней пересылке данных.
  
 
== Входные и выходные данные алгоритма ==
 
== Входные и выходные данные алгоритма ==
Входные данные: вектор действительных чисел <math>a = (a_1,a_2,...,a_N)</math>.
+
'''Входные данные:''' вектор действительных/комплексных чисел <math>a = (a_1,a_2,...,a_N)</math>.
 +
 
 +
'''Выходные данные:''' вектор комплексных чисел <math>b = (b_1,b_2,...,b_{\lfloor N/2 \rfloor+1})</math>.
 +
 
 +
Есть несколько ситуаций, при которых количество потребляемой памяти для хранения входного вектора <math>a</math> и выходного вектора <math>b</math> тможет быть снижено вдвое:
 +
# Если <math>a</math> — вектор действительных чисел, тогда <math>b = E\bar{b}</math>;
 +
# Если <math>a</math> — вектор чисто-мнимых чисел, тогда <math>b = -E\bar{b}</math>;
 +
# Если <math>a = E\bar{a}</math>, тогда <math>b</math> — вектор действительных чисел;
 +
# Если <math>a = -E\bar{b}</math>, то <math>b</math> — вектор чисто-мнимых чисел.
  
Выходные данные: вектор комплексных чисел <math>b = (b_1,b_2,...,b_{\lfloor N/2 \rfloor+1})</math>.
+
Где матрица <math>E</math>:
 +
:<math>
 +
  E = \begin{bmatrix}
 +
    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{bmatrix}
 +
</math>
  
 
== Свойства алгоритма ==
 
== Свойства алгоритма ==
Таким образом в общем случае у алгоритма БПФ последовательная и параллельная сложности в случае неограниченных ресурсов - ''линейна''.
+
Таким образом в общем у алгоритма БПФ в случае неограниченных ресурсов:
 
+
*последовательная сложность - ''линейно-логарифмическая''
При этом вычислительная мощность алгоритма БПФ, как отношение числа операций к суммарному объему входных и выходных данных – ''логарифмическая''.
+
*параллельная сложность  - ''логарифмическая''.  
  
 
При этом алгоритм полностью детерминирован.
 
При этом алгоритм полностью детерминирован.
  
 
= ЧАСТЬ. Программная реализация алгоритма =
 
= ЧАСТЬ. Программная реализация алгоритма =
 +
В качестве тестируемой реализации была выбрана библиотека FFTW, а именно функция, выполняющая БПФ одномерного вектора действительных чисел.
 +
==  Особенности последовательной реализации алгоритма ==
 +
Простейший вариант алгоритма для степеней двойки на языке C:
 +
<source lang="C++">
 +
#include <stdio.h>
 +
#include <math.h>
 +
#include <complex.h>
  
 +
typedef double complex cplx;
 +
 +
void _fft(cplx *buf, cplx *out, int n, int step)
 +
{
 +
    if (step < n) {
 +
_fft(out, buf, n, step * 2); //рекурсивный рассчет
 +
_fft(out + step, buf + step, n, step * 2);
 +
 +
for (int i = 0; i < n; i += 2 * step) {
 +
    cplx t = cexp(-I * PI * i / n) * out[i + step];
 +
    buf[i / 2]    = out[i] + t;
 +
    buf[(i + n)/2] = out[i] - t;
 +
        }
 +
    }
 +
}
 +
 +
void fft(cplx *buf, int n)
 +
{
 +
    cplx out[n];
 +
    for (int i = 0; i < n; i++) out[i] = buf[i];
 +
 +
    _fft(buf, out, n, 1);
 +
}
 +
</source>
 +
 +
== Локальность данных и вычислений ==
 +
== Возможные способы и особенности параллельной реализации алгоритма ==
 
== Масштабируемость алгоритма и его реализации ==
 
== Масштабируемость алгоритма и его реализации ==
 +
Тестирование проводилось на вычислительном комплексе IBM Blue Gene/P
 +
 +
=== Описание платформы 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-картам через сеть коллективных операций
 +
 +
=== Оценка масштабируемости алгоритма ===
 +
Масштабируемость алгоритма проверялась с помощью тестовой программы расположенной по [https://github.com/alesapin/parallel_fftw ссылке]. Для эксперимента использовались:
 +
* Компилятор mpicc
 +
* Библиотека OpenMPI версии 1.68
 +
* Библиотека FFTW 3.3.5
 +
Алгоритм тестировался на входных векторах размера <math>2^{15}</math> до <math>2^{26}</math>. При этом использовались от 8 до 128 процессоров с шагом степени двойки.
 +
 +
[[File:Fftw time good cr.png  |center|thumb|800px| Замеры времени работы параллельной реализации FFTW для входного вектора действительных чисел. По оси Data Size указано <math>N</math> -- число входных элементов. По оси Procs Num указано <math>p</math> -- число задействованных процессоров. По оси Elapsed Time -- время затраченное на выполненине алгоритма в секундах.]]
 +
 +
[[File:Fftw efficency good cr.png |center|thumb|800px| Эффективность параллельной реализации FFTW для входного вектора действительных чисел. По оси Data Size указано <math>N</math> -- число входных элементов. По оси Procs Num указано <math>p</math> -- число задействованных процессоров. По оси Effeciency -- эффективность алгоритма в процентах.]]
 +
 +
Из графиков можно сделать следующие выводы:
 +
* Значительное увеличение числа процессоров не дает значимого уменьшения времени исполнения, что отражено и на графике эффективности. Это означает, что время, в основном, расходуется на пересылки данных, а не на рассчет.
 +
* При большом количестве элементов и маленьком числе процессоров наблюдается сверхлинейное ускорение, что связано с маленьким число пересылок и активным использованием кеша.
 +
 +
== Динамические характеристики и эффективность реализации алгоритма ==
 +
== Выводы для классов архитектур ==
 
== Существующие реализации алгоритма ==
 
== Существующие реализации алгоритма ==
 +
Наиболее известной библиотекой для выполнения различных вариантов быстрого преобразования Фурье является [http://www.fftw.org/ FFTW], разработанная Матео Фриго и Стивеном Джонсоном в Массачусетском технологическом институте. Для преобразования небольших объемов данных используется алгоритм [https://en.wikipedia.org/wiki/Cooley%E2%80%93Tukey_FFT_algorithm Кули - Тьюки], в случае входа больших размеров используется алгоритм [https://en.wikipedia.org/wiki/Rader%27s_FFT_algorithm Радера] или [https://en.wikipedia.org/wiki/Chirp_Z-transform#Bluestein.27s_algorithm Блюштейна]. Библиотека считается самой быстрой реализацией БПФ и на входах любых размеров и размерностей сохраняет ассимптотическую сложность <math>O(N\cdot log(N))</math>. Также, существует версия FFTW с поддержкой [http://www.fftw.org/doc/Distributed_002dmemory-FFTW-with-MPI.html#Distributed_002dmemory-FFTW-with-MPI MPI], которая позволяет с объемами данных не умещающимися в память.
 +
 +
FFTW изначально написана на языке C, однако имеет API для нескольких языков, в том числе Fortran и Ada. Открытая версия библиотеки распространяется под лицензией [https://en.wikipedia.org/wiki/GNU_General_Public_License GNU General Public License]. Также существует реализация для коммерческого использования распространяемая под лицензией [https://en.wikipedia.org/wiki/MIT MIT].
 +
 +
Существует несколько реализаций алгоритма, для рассчета на GPU:
 +
* [https://developer.nvidia.com/cufft cuFFT] для рассчета на графических ускорителях от NVIDIA.
 +
* [http://gpuopen.com/compute-product/clfft/ clFFT] для рассчета на GPU и CPU.
  
 
= Литература =
 
= Литература =
 
[1] Википедия [Электронный ресурс]. Тема: Быстрое преобразование Фурье – Электрон. дан. – URL [https://ru.wikipedia.org/wiki/%D0%91%D1%8B%D1%81%D1%82%D1%80%D0%BE%D0%B5_%D0%BF%D1%80%D0%B5%D0%BE%D0%B1%D1%80%D0%B0%D0%B7%D0%BE%D0%B2%D0%B0%D0%BD%D0%B8%D0%B5_%D0%A4%D1%83%D1%80%D1%8C%D0%B5 Быстрое преобразование Фурье] (дата обращения 17.09.2016)
 
[1] Википедия [Электронный ресурс]. Тема: Быстрое преобразование Фурье – Электрон. дан. – URL [https://ru.wikipedia.org/wiki/%D0%91%D1%8B%D1%81%D1%82%D1%80%D0%BE%D0%B5_%D0%BF%D1%80%D0%B5%D0%BE%D0%B1%D1%80%D0%B0%D0%B7%D0%BE%D0%B2%D0%B0%D0%BD%D0%B8%D0%B5_%D0%A4%D1%83%D1%80%D1%8C%D0%B5 Быстрое преобразование Фурье] (дата обращения 17.09.2016)
  
 +
[2] Бахвалов Н. С.,  Жидков Н. П., Кобельков. Г. М. — 6-е изд. — М. : БИНОМ. Лаборатория знаний, 2008. — 636 с.
 
[[en:Fast Discrete Fourier Transform]]
 
[[en:Fast Discrete Fourier Transform]]
 +
 +
[3] FFTW [Электронный ресурс]. Тема: One-Dimensional DFTs of Real Data – Электрон. дан. – URL [http://www.fftw.org/doc/One_002dDimensional-DFTs-of-Real-Data.html#One_002dDimensional-DFTs-of-Real-Data One-Dimensional DFTs of Real Data] (последняя дата обращения 14.19.2016)
 +
 +
[4] FFTW [Электронный ресурс]. – Электрон. дан. – URL [http://www.fftw.org/ FFTW Lib] (последняя дата обращения 14.19.2016)

Текущая версия на 17:09, 29 ноября 2016

Symbol confirmed.svgЭта работа успешно выполнена
Преподавателю: в основное пространство, в подстраницу

Данное задание было проверено и зачтено.
Проверено Kronberg и IgorS.



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

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

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

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

Содержание

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

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

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

[math] \hat{s}(\omega)=\frac{1}{\sqrt{2\pi}}\int\limits_{-\infty}^{\infty}s(t)e^{-it\omega}\,dt [/math]

В отличие от своего классического аналога, БПФ имеет сложность [math]O(N\log{N})[/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_k = \sum_{n=0}^{N-1} x_n e^{-\frac{2\pi i}{N} nk}, k = \overline{0,N-1}[/math]

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

[math] \begin{matrix} X_k & = & \sum \limits_{m=0}^{N/2-1} x_{2m}e^{-\frac{2\pi i}{N} (2m)k} + \sum \limits_{m=0}^{N/2-1} x_{2m+1} e^{-\frac{2\pi i}{N} (2m+1)k} \end{matrix} [/math]

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

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

Алгоритм:

  1. Входной вектор [math]a = (a_0,a_2,...,a_{N-1})[/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] A = \begin{pmatrix} a_0 & a_1 & \cdots & a_{n_1-1} \\ a_{n_1} & a_{n_1+1} & \cdots & a_{2n_1-1} \\ \vdots & \vdots & \ddots & \vdots \\ a_{(n_2-1)\cdot n_1} & a_{(n_2-1)\cdot n_1+1} & \cdots & a_{n_2\cdot n_1-1} \end{pmatrix} [/math]

  1. К каждой строке полученной матрицы применяется быстрое дискретное преобразование Фурье (БПФ) порядка [math]n_1[/math]
  2. Каждый элемент полученный после применения БПФ умножается на поворотные множители [math]exp (2 \pi i(m-1)(j-1)/N)[/math], где [math]m[/math] - номер строки, а [math]j[/math] - номер столбца
  3. Полученная после шагов 1-2 матрица [math]A[/math] транспонируется
  4. К каждой строке матрицы [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]

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

Рис. 1. "Бабочка" для двумерного входа

Для 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]

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

Рис. 2. "Бабочка" для четырёхмерного входа

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

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] БПФ представляется как:

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

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

Схема организации состоит в том, что на каждом шаге [math]i[/math] для выполнения "бабочки" все элементы разбиваются на [math]n_{i_1}[/math] векторов по [math]n_{i_2}[/math] элементов, причем [math]n_{i_1} \cdot n_{i_2} = n_i[/math], где [math]n_i[/math] длина входного вектора на текущем шаге [math]i[/math]. В случае если такое разбиение невозможно, по причине того что [math]n_i[/math] простое число, исходный вектор на текущем шаге дополняется элементом. В зависимости от номера шага, разница координат для каждой пары элементов увеличивается соразмерно разбиению [math](n_{i_1},n_{i_2})[/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_{N-1})[/math] в вектор комплексных чисел [math]b = (b_0, \dots, b_{N-1})[/math], такой, что [math]b_i=\sum_{j=0}^{N-1}a_j\varepsilon^{ij}[/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_{n_2+i}, \dots, a_{n_2(n_1-1)+i}[/math] для [math]i=0,1,\dots,n_2-1[/math]. Покажем теперь, что за [math]O(Np)[/math] действий можно решить исходную задачу. Заметим, что [math]b_i=\sum_{j=0}^{n_2-1} \varepsilon^{ij} \left(\sum_{k=0}^{n_1-1}a_{kn_2+j}\varepsilon^{kin_2}\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]4N\gt 2^k\ge2N[/math]. Заметим, что тогда [math]b_i=\varepsilon^{-i^2/2}\sum_{j=0}^{N-1}\varepsilon^{(i+j)^2/2}\varepsilon^{-j^2/2}a_j[/math]. Обозначим [math]\bar{a}_i=\varepsilon^{-i^2/2}a_i[/math], [math]\bar{b}_i=\varepsilon^{i^2/2}b_i[/math], [math]c_i=\varepsilon^{(2N-2-i)^2/2}[/math]. Тогда [math]\bar{b}_i=\sum_{j=0}^{2N-2-i}\bar{a}_jc_{2N-2-i-j}[/math], если положить [math]\bar{a}_i=0[/math] при [math]i\ge N[/math].

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

Вычисления всех [math]\bar{a}_i[/math] и [math]c_i[/math] требуют [math]O(N)[/math] действий, три преобразования Фурье требуют [math]O(N\log(N))[/math] действий, перемножение результатов преобразований Фурье требует [math]O(N)[/math] действий, вычисление всех [math]b_i[/math] зная значения свертки требует [math]O(N)[/math] действий. Итого для дискретного преобразования Фурье требуется [math]O(N\log(N))[/math] действий для любого [math]N[/math].

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

Информационный граф алгоритма Кули-Тьюки для [math]N = 21[/math]. Исходные данные обозначены зеленым, результат — синим. Преобразования Фурье для [math]n_1 = 7[/math] и [math]n_2 = 3 \text{ } (N = n_1 \cdot n_2)[/math] представлены как «чёрные ящики». Умножение на поворотные коэффициенты (множители) представлено коричневыми ромбами.

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

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_{\lfloor N/2 \rfloor+1})[/math].

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

  1. Если [math]a[/math] — вектор действительных чисел, тогда [math]b = E\bar{b}[/math];
  2. Если [math]a[/math] — вектор чисто-мнимых чисел, тогда [math]b = -E\bar{b}[/math];
  3. Если [math]a = E\bar{a}[/math], тогда [math]b[/math] — вектор действительных чисел;
  4. Если [math]a = -E\bar{b}[/math], то [math]b[/math] — вектор чисто-мнимых чисел.

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

[math] E = \begin{bmatrix} 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{bmatrix} [/math]

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

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

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

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

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

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

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

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

#include <stdio.h>
#include <math.h>
#include <complex.h>

typedef double complex cplx;
 
void _fft(cplx *buf, cplx *out, int n, int step)
{
    if (step < n) {
	_fft(out, buf, n, step * 2); //рекурсивный рассчет
	_fft(out + step, buf + step, n, step * 2);
 
	for (int i = 0; i < n; i += 2 * step) {
	    cplx t = cexp(-I * PI * i / n) * out[i + step];
	    buf[i / 2]     = out[i] + t;
	    buf[(i + n)/2] = out[i] - t;
        }
    }
}
 
void fft(cplx *buf, int n)
{
    cplx out[n];
    for (int i = 0; i < n; i++) out[i] = buf[i];
 
    _fft(buf, out, n, 1);
}

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^{15}[/math] до [math]2^{26}[/math]. При этом использовались от 8 до 128 процессоров с шагом степени двойки.

Замеры времени работы параллельной реализации FFTW для входного вектора действительных чисел. По оси Data Size указано [math]N[/math] -- число входных элементов. По оси Procs Num указано [math]p[/math] -- число задействованных процессоров. По оси Elapsed Time -- время затраченное на выполненине алгоритма в секундах.
Эффективность параллельной реализации FFTW для входного вектора действительных чисел. По оси Data Size указано [math]N[/math] -- число входных элементов. По оси Procs Num указано [math]p[/math] -- число задействованных процессоров. По оси Effeciency -- эффективность алгоритма в процентах.

Из графиков можно сделать следующие выводы:

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

2.5 Динамические характеристики и эффективность реализации алгоритма

2.6 Выводы для классов архитектур

2.7 Существующие реализации алгоритма

Наиболее известной библиотекой для выполнения различных вариантов быстрого преобразования Фурье является FFTW, разработанная Матео Фриго и Стивеном Джонсоном в Массачусетском технологическом институте. Для преобразования небольших объемов данных используется алгоритм Кули - Тьюки, в случае входа больших размеров используется алгоритм Радера или Блюштейна. Библиотека считается самой быстрой реализацией БПФ и на входах любых размеров и размерностей сохраняет ассимптотическую сложность [math]O(N\cdot log(N))[/math]. Также, существует версия FFTW с поддержкой MPI, которая позволяет с объемами данных не умещающимися в память.

FFTW изначально написана на языке C, однако имеет API для нескольких языков, в том числе Fortran и Ada. Открытая версия библиотеки распространяется под лицензией GNU General Public License. Также существует реализация для коммерческого использования распространяемая под лицензией MIT.

Существует несколько реализаций алгоритма, для рассчета на GPU:

  • cuFFT для рассчета на графических ускорителях от NVIDIA.
  • clFFT для рассчета на GPU и CPU.

3 Литература

[1] Википедия [Электронный ресурс]. Тема: Быстрое преобразование Фурье – Электрон. дан. – URL Быстрое преобразование Фурье (дата обращения 17.09.2016)

[2] Бахвалов Н. С., Жидков Н. П., Кобельков. Г. М. — 6-е изд. — М. : БИНОМ. Лаборатория знаний, 2008. — 636 с.

[3] FFTW [Электронный ресурс]. Тема: One-Dimensional DFTs of Real Data – Электрон. дан. – URL One-Dimensional DFTs of Real Data (последняя дата обращения 14.19.2016)

[4] FFTW [Электронный ресурс]. – Электрон. дан. – URL FFTW Lib (последняя дата обращения 14.19.2016)