Участник:1: различия между версиями

Материал из Алговики
Перейти к навигации Перейти к поиску
 
(не показаны 3 промежуточные версии этого же участника)
Строка 1: Строка 1:
= ЧАСТЬ. Свойства и структура алгоритмов =
 
== Свойства и структура алгоритма ==
 
  
=== Общее описание алгоритма ===
 
 
Главная идея БПФ состоит в сокращении многочисленных повторов сложения и умножения, которые выполняются в ДПФ (повторяющиеся множители: <math>e^{\frac{-2kpi\pi}{n}}</math> Например при условии, что  <math>(pk)mod(n)=const</math>, при всех <math>k=\overline{0,n-1}</math> и <math>p=\overline{0,n-1}</math>
 
 
Напомним формулу ДПФ:
 
 
<math> y_{p}=\sum_{k=0}^{n-1}x_{k}e^{\frac{-2kpi\pi}{n}}, p=\overline{0,n-1}</math>
 
 
За счет использования данного факта, происходит значительное сокращение количества вычислительных операций. В данной статье используется алгоритм Cooley – Tukey, который предполагает собой разбиение исходного множества входных сигналов на два равных подмножества. При этом множество входных сигналов должно быть степенью двойки. Далее операция деления множества продолжается для образовавшихся подмножеств, и так далее. В результате над каждым из полученных подмножеств выполняется БПФ. Данные вычисления можно выполнять независимо друг от друга. Этот факт используется для реализации параллельной версии алгоритма. В общем алгоритм Cooley – Tukey строится из двух основных шагов:
 
 
* Преобразование входного массива данных (бит-реверсирование).
 
* Вычисление БПФ.
 
Приведем более детальное описание БПФ (математическое описание):
 
*Пусть <math>e^{\frac{-2kpi\pi}{n}}</math> = <math>\omega _{n}^{kp}</math>.
 
*Разобьем исходный сигнал <math>x(n)</math> на два подмножества, с четными <math>x_{1}(2n)</math> и нечетными <math>x_{2}(2n+1)</math> номерами.
 
Тогда получаем: <math>y_{p}=\sum_{k=0}^{n-1}x_{k}\omega _{n}^{kp}= \sum_{k=0}^{\frac{n}{2}-1}x_{2k}\omega _{n}^{2kp}+\sum_{k=0}^{\frac{n}{2}-1}x_{2k+1}\omega _{n}^{(2k+1)p}</math>.
 
 
Выпишем два свойства омеги:
 
*<math>\omega_{\frac{n}{2}}=e^{2(\frac{-2i\pi}{n})}=\omega_{n}^{2}</math>
 
*<math>\omega_{n}^{\frac{n}{2}}=e^{\frac{-2ni\pi}{2n}}=e^{-i\pi}=-1</math>
 
Воспользуемся первым: <math>y_{p}= \sum_{k=0}^{\frac{n}{2}-1}x_{2k}\omega _{n}^{2kp}+\omega_{n}^{p}\sum_{k=0}^{\frac{n}{2}-1}x_{2k+1}\omega _{n}^{2kp}=\sum_{k=0}^{\frac{n}{2}-1}x_{2k}\omega _{\frac{n}{2}}^{kp}+\omega_{n}^{p}\sum_{k=0}^{\frac{n}{2}-1}x_{2k+1}\omega _{\frac{n}{2}}^{kp}</math>. Или, если обозначить <math>y_{1}(p)=\sum_{k=0}^{\frac{n}{2}-1}x_{2k}\omega _{\frac{n}{2}}^{kp}</math>, а <math>y_{2}(p)=\sum_{k=0}^{\frac{n}{2}-1}x_{2k+1}\omega _{\frac{n}{2}}^{kp}</math>, то получим формулу для вычисления координат изображения БПФ при <math>p=\overline{0...n/2}</math>: <math>y_(p)=y_{1}(p)+\omega_{n}^{p}y_{2}(p)</math>.
 
 
Заметим, что <math>y(p); y_{1}(p); y_{2}(p)</math> обладают свойством периодичности, с периодом <math>\frac{n}{2}</math>. Покажем это для <math>y_{1}(p)</math> : <math>y_{1}(p+\frac{n}{2})=\sum_{k=0}^{\frac{n}{2}-1}x_{2k}\omega_{\frac{n}{2}}^{k(p+\frac{n}{2})}=\sum_{k=0}^{\frac{n}{2}-1}x_{2k}\omega_{\frac{n}{2}}^{kp}\omega_{\frac{n}{2}}^{\frac{kn}{2}}=\sum_{k=0}^{\frac{n}{2}-1}x_{2k}\omega_{\frac{n}{2}}^{kp}</math>
 
 
Так как: <math>\omega_{\frac{n}{2}}^{\frac{kn}{2}}=e^{\frac{-4kni\pi}{2n}}=e^{-2ki\pi}=1</math>
 
 
Далее необходимо определить <math>y(p)</math> при <math>p\geq \frac{n}{2}</math>. Теперь воспользуемся свойством периодичности обнаруженными нами ранее: <math>y(p+\frac{n}{2})=y_{1}(p+\frac{n}{2})+\omega_{n}^{(p+\frac{n}{2})}y_{2}(p+\frac{n}{2})</math>. Воспользуемся вторым свойством омеги: <math>y(p)=y_{1}(p)-\omega_{n}^{p}y_{2}(p)</math>.
 
 
В конечном итоге мы пришли к конечной формуле для первого этапа вычисления БПФ. Дальнейшие этапы делят подмножества полученные на предыдущим этапе вычисления БПФ.
 
 
=== Математическое описание алгоритма ===
 
 
'''Бит-реверсирование'''
 
 
Бит-реверсирование – это преобразование двоичного числа путем изменения порядка следования бит в нем на противоположный. Бит-реверсирование применимо только к индексам элементов массива и предназначено для изменения порядка следования этих элементов, при этом значения самих элементов не изменяются. Так, к примеру, четырехмерный вектор <math>(x_{0},x_{1},x_{2},x_{3})</math> изменится на <math>(x_{0},x_{2},x_{1},x_{3})</math>. Данная процедура необходима для разделения исходного сигнала на элементы подмножеств с четными и нечетными номерами, которые необходимы для вычисления <math>y_{1}(p)</math> и <math>y_{2}(p)</math>.
 
 
'''Вычисление БПФ'''
 
 
После выполнения этапа бит-реверсирования выполняются вычисления проиллюстрированные ''на рисунке 2'' в разделе [[#Информационный граф|Информационный граф]]. Рисунок 2 приведен для случая четырехмерного входного сигнала. Он предполагает совершение двух этапов. В общем же количество этапов равно <math>log_{2}n</math>. Каждый этап содержит в себе вычисление одной и той же операции "бабочки". Где на вход подаются два числа <math>a</math> и <math>b</math>, на выходе же получаются <math>a+bU</math> и <math>a-bU</math>. Где <math>U</math>- это ''коэффициент поворота''. Общая формула для коэффициента поворота: <math>e^{\frac{-ji\pi}{k}}</math>. Где <math>k</math>- размер шага "бабочки", <math>j</math>- номер итерации вычислений над группой элементов, размером <math>2k</math> <math>(0\leq j < k )</math>.
 
 
=== Вычислительное ядро алгоритма ===
 
 
Из двух вышеописанных этапов БПФ, бит-реверсирование является самым быстровыполнимым. Поэтому вычислительное ядро алгоритма будет лежать именно во втором этапе. А именно, им будет являться вычисление "бабочек". Количество этапов выполнения БПФ, как описывалось выше, равно <math>log_{2}n</math>. Учитывая количество этапов, размерность входного сигнала и идею деления на подмножества, всего "бабочек" получается <math>\frac{nlog_{2}n}{2}</math>. При этом операция умножения выполняется <math>\frac{n}{2}</math> раз. Это связано с тем, что результаты вычисления "бабочки" <math>a+bU</math> и <math>a-bU</math>.
 
 
=== Макроструктура алгоритма ===
 
 
Из всего вышесказанного заключаем, что макроструктуру алгоритма легче всего описать при помощи рекурсии:
 
*<math>\frac{n}{2}</math> преобразований Фурье порядка 2
 
*умножение <math>\frac{n}{2}</math> пар комплексных чисел
 
*<math>n</math> комплексных сложений
 
*снова преобразование Фурье
 
 
=== Схема реализации последовательного алгоритма ===
 
 
Программный код (описан на языке С++) строится из ряда функций:
 
*бит-реверсирование
 
*инициализация размерности входного сигнала
 
*инициализация входного сигнала (при помощи генератора случайных чисел)
 
*процесс очистки памяти
 
*выполнение операции "бабочка"
 
*выполнение БПФ
 
*функция вывода результата (для тестирования)
 
----
 
'''главный программный код БПФ:'''
 
<source lang="С++">
 
int main()
 
{
 
complex<double> *inputSignal = NULL; 
 
complex<double> *outputSignal = NULL;
 
int size = 0;
 
cout << "FFT. Serial version." << endl;     
 
ProcessInitialization(inputSignal, outputSignal, size);
 
bit(inputSignal, outputSignal, size);
 
SerialFFTCalculation(outputSignal, size);
 
PrintSignal(outputSignal, size); //в случае тестирования
 
ProcessTermination(inputSignal, outputSignal);
 
cin.get();
 
}
 
</source>
 
----
 
==== функция бит-реверсирования ====
 
<source lang="С++">
 
void bit(complex<double> inputSignal[], complex<double> *outputSignal, int size)
 
{ // процесс бит-реверсирования
 
int j = 0, i = 0;
 
while (i < size)
 
{
 
if (j > i)
 
{
 
outputSignal[i] = inputSignal[j];     
 
outputSignal[j] = inputSignal[i];
 
}
 
else     
 
if (j == i)         
 
outputSignal[i] = inputSignal[i];
 
 
int m = size >> 1;   
 
while ((m >= 1) && (j >= m))
 
{
 
j -= m;     
 
m = m >> 1;
 
}   
 
j += m;   
 
i++;
 
}
 
}
 
</source>
 
==== инициализация размерности входного сигнала ====
 
<source lang="С++">
 
void ProcessInitialization(complex<double>* &inputSignal, complex<double>* &outputSignal, int &size)
 
{
 
// инициализация размера входного сигнала 
 
do
 
{
 
cout << "Enter the input signal length: ";
 
cin >> size;
 
if (size < 4)
 
cout << "Input signal length should be >= 4" << endl;
 
else
 
{
 
int tmpSize = size;
 
while (tmpSize != 1)
 
{
 
if (tmpSize % 2 != 0)
 
{
 
cout << "Input signal length should be powers of two" << endl;
 
size = -1;
 
break;
 
}
 
tmpSize /= 2;
 
}
 
}
 
}
 
while (size < 4); 
 
cout << "Input signal length = " << size << endl;
 
inputSignal = new complex<double>[size]; 
 
outputSignal = new complex<double>[size];
 
DummyDataInitialization(inputSignal, size);
 
}
 
</source>
 
==== инициализация входного сигнала ====
 
<source lang="С++">
 
void DummyDataInitialization(complex<double>* mas, int size)
 
{
 
for (int i = 0; i < size; i++)
 
{
 
mas[i] = int(rand() % 100);
 
}
 
}
 
</source>
 
==== процесс очистки памяти ====
 
<source lang="С++">
 
void ProcessTermination(complex<double>* &inputSignal, complex<double>* &outputSignal)
 
{
 
//очистка памяти
 
delete[] inputSignal; 
 
inputSignal = NULL; 
 
delete[] outputSignal; 
 
outputSignal = NULL;
 
}
 
</source>
 
==== выполнение операции "бабочка" ====
 
<source lang="С++">
 
void Butterfly(complex<double> *signal, complex<double> &u,int offset, int butterflySize)
 
{
 
complex<double> tem = signal[offset + butterflySize] * u;
 
signal[offset + butterflySize] = signal[offset] - tem; 
 
signal[offset] += tem;
 
}
 
</source>
 
==== выполнение БПФ ====
 
<source lang="С++">
 
void SerialFFTCalculation(complex<double> *signal, int size)
 
{
 
int m = 0; 
 
for (int tmp_size = size; tmp_size > 1; tmp_size /= 2, m++);
 
for (int p = 0; p < m; p++)
 
{
 
int butterflyOffset = 1 << (p + 1);
 
int butterflySize = butterflyOffset >> 1;
 
double coeff = M_PI / butterflySize;
 
 
for (int i = 0; i < size / butterflyOffset; i++)
 
for (int j = 0; j < butterflySize; j++)
 
Butterfly(signal,complex<double>(cos(-j*coeff),sin(-j*coeff)),j+i*butterflyOffset,butterflySize);
 
}
 
}
 
</source>
 
==== функция вывода результата (для тестирования) ====
 
<source lang="С++">
 
void PrintSignal(complex<double> *signal, int size)
 
{
 
cout << "Result signal" << endl; 
 
for (int i = 0; i < size; i++)   
 
cout << signal[i] << endl;
 
}
 
</source>
 
=== Последовательная сложность алгоритма ===
 
Рассмотренный алгоритм  Cooley – Tukey выполняется за <math>nlog_{2}n</math> операций комплексного сложения и <math>\frac{n}{2}log_{2}n</math> операций комплексного умножения. Таким образом последовательная сложность алгоритма является '''линейно-логарифмической'''.
 
=== Информационный граф ===
 
[[file:Граф.jpg|center|thumb|800px|Рисунок 1. Алгоритм Cooley – Tukey для n=4, при предварительно сделанном бит-реверсировании]]
 
=== Ресурс параллелизма алгоритма ===
 
[[file:Парал.jpg|center|thumb|800px|Рисунок 2. Распараллеливание алгоритма Cooley – Tukey для n=8]]
 
Рисунок 2 получен из информационного графа, главной идеей которого является процедура деления множеств полученных на предыдущем шаге на два равных по мощности подмножества и применение к ним БПФ. То есть процесс вычисления распределяется между <math>\frac{n}{2}</math> процессорами. И так как при последовательной реализации алгоритма было <math>\frac{n}{2}\log_{2}n</math> операций комплексного умножения и <math>n\log_{2}n</math> комплексного сложения/вычитания. Получается, что при параллельной реализации алгоритма мы будем иметь <math>\log_{2}n</math> операций комплексного умножения и <math>2\log_{2}n</math> комплексного сложения/вычитания. При этом, если используется функция HT(Hyper threading), то параллельная сложность быстрого преобразования Фурье (базовый алгоритм Кули-Тьюки) для векторов с длиной, равной степени двойки составляет <math>\log_{2}n</math>.
 
=== Входные и выходные данные алгоритма ===
 
*входные (вектор): <math>\overline{a}=(a_{1},a_{2},...,a_{n})</math>
 
*выходные данные (вектор): <math>\overline{b}=(b_{1},b_{2},...,b_{n})</math>
 
Данные могут быть как действительными так и комплексными.
 
=== Свойства алгоритма ===
 
*Соотношение последовательной и параллельной сложности является ''линейным''. Это хорошо видно из пункта [[#Ресурс параллелизма алгоритма|ресурс параллелизма алгоритма]].
 
*Вычислительная мощность алгоритма является ''логарифмической''.
 
*Устойчивость.
 
**Рассмотрим вопрос об устойчивости поближе. Так как БПФ отличается от ДПФ простым лишь сокращением одинаковых операций, то вопрос устойчивости можно рассматривать на примере ДПФ. Так, к примеру, возьмем формулу ДПФ и добавим возмущение: <math> y_{p}=\sum_{k=0}^{n-1}(x_{k}+\varepsilon )e^{\frac{-2*\pi*k*p*i}{n}}, p=\overline{0,n-1}</math>. После чего преобразуем ее: <math>y_{p}=\sum_{k=0}^{n-1}x_{k}e^{\frac{-2*\pi*k*p*i}{n}}+\varepsilon\sum_{k=0}^{n-1} e^{\frac{-2*\pi*k*p*i}{n}}</math> и оценим вторую сумму, которая является геометрической прогрессией со знаменателем <math>q=e^{\frac{-2*\pi*p*i}{n}}</math>. Таким образом оцениваемая сумма будет равняться: <math>\varepsilon\frac{1-e^{\frac{-2*\pi*p*i*(n-1))}{n}}}{1-e^{\frac{-2*\pi*p*i}{n}}}</math>. При достаточно больших размерностях входного сигнала данное выражение сводится к: <math>\varepsilon (1-e^{-2*\pi*p*i})</math>. В тригонометрической форме данного выражения: <math>\varepsilon (cos(\frac{2*\pi*p*(n-1)}{n})-i*sin(\frac{2*\pi*p*(n-1)}{n}))</math> видно, что принимая во внимание целочисленность индексов <math>p</math> получается 0. Таким образом, мы можем говорить об устойчивости алгоритма в целом, и абсолютной устойчивости алгоритма в асимптотике.
 
*исходя из информационного графа, можно сделать вывод о сбалансированности алгоритма.
 
*алгоритм полностью детерминирован.
 
*Степень исхода вершины информационного графа не превышает двух.
 
 
= ЧАСТЬ. Программная реализация алгоритма =
 
== Масштабируемость алгоритма и его реализации ==
 
[[file:масш.jpg|center|thumb|1000px|Рисунок 3. Параллельная реализация алгоритма. Изменение эффективности(%)в зависимости от числа процессоров и размера области.]]
 
Как описывалось выше, каждый процессор получает достаточно низкий объем работы, что приводит к неэффективному распараллеливанию алгоритма (высокие затраты на передачу данных между процессорами по сравнению с вычислениями на них). Поэтому необходимо находить оптимальное количество процессоров для каждой отдельно взятой задачи. В данном примере, при размерности входного сигнала менее 262144, оптимальное количество процессоров составляет 1. При размерности 262144-4194304 оптимум достигался при 8 процессорах. И при размерности более 4194304- 16 процессоров. Вычисления проводились на суперкомпьютере «Ломоносов».
 
== Существующие реализации алгоритма ==
 
Самой распространенной и известной реализацией дискретного БПФ является библиотека [http://www.fftw.org/ FFTW], разработанная Маттео Фриго и Стивеном Г. Джонсоном в Массачусетском технологическом институте. Данная библиотека существует для языков C и FORTRAN и позволяет производить вычисления в одном или нескольких измерениях. Входной сигнал может быть как действительным так и комплексным, а его размерность не обязательно должна быть степенью двойки. На данный момент последняя версия FFTW 3.3.5. Другие особенности библиотеки:
 
*Скорость. (Поддержка SSE / SSE2 / Altivec, начиная с версии 3.0. Версия 3.3.1 поддерживает AVX и ARM Neon.)
 
*Параллельные вычисления (OpenMP, MPI).
 
*[http://www.fftw.org/fftw3_doc/ Документация] в форматах HTML, PDF.
 
*Данное программное обеспечение является свободным и распространяется по лицензии GNU (General Public License).
 
Здесь приводится пример реализации прямого БПФ для комплексных данных входного сигнала на языке C++.
 
<source lang="С++">
 
#include <fftw3-mpi.h>
 
# include <stdlib.h>
 
# include <stdio.h>
 
#include <sys/stat.h>
 
#include <fcntl.h>
 
# include <time.h>
 
#include <math.h>
 
 
int main(int argc, char **argv)
 
{
 
const ptrdiff_t N0 = 4194304;
 
fftw_plan planForw;
 
fftw_complex *data, *dataOut;
 
ptrdiff_t alloc_local, local_ni, local_i_start, i, j, local_no, local_o_start;
 
int index, size;
 
double startwtime, endwtime;
 
MPI_Init(&argc, &argv);
 
fftw_mpi_init();
 
MPI_Comm_rank(MPI_COMM_WORLD, &index);
 
MPI_Comm_size(MPI_COMM_WORLD, &size);
 
 
/* get local data size and allocate */
 
alloc_local = fftw_mpi_local_size_1d(N0, MPI_COMM_WORLD, FFTW_FORWARD, FFTW_ESTIMATE,
 
&local_ni, &local_i_start, &local_no, &local_o_start);
 
data = fftw_alloc_complex(alloc_local);
 
dataOut = fftw_alloc_complex(alloc_local);
 
/* create plan  */
 
planForw = fftw_mpi_plan_dft_1d(N0, data, dataOut, MPI_COMM_WORLD,
 
FFTW_FORWARD, FFTW_ESTIMATE);
 
/* initialize data to some function my_function(x,y) */
 
for (i = 0; i < local_ni; ++i)
 
{
 
data[i][0] = rand() / (double)RAND_MAX;
 
data[i][1] = rand() / (double)RAND_MAX;
 
}
 
if (index == 0){
 
startwtime = MPI_Wtime();
 
 
}
 
 
fftw_execute(planForw);
 
if (index == 0){
 
endwtime = MPI_Wtime();
 
printf("wall clock time = %f\n",
 
endwtime - startwtime);
 
 
 
}
 
 
fftw_destroy_plan(planForw);
 
MPI_Finalize();
 
}
 
</source>
 
[http://www.fftw.org/download.html Скачав] библиотеку и установив ее на рабочий компьютер, можно поэкспериментировать, используя вышеописанный исходный код. При этом, на стадии компиляции необходимо указывать директорию, где находятся заголовочные файлы и директорию, где находятся библиотечные файлы. Пример: mpicxx name.cpp -I/path/to/include/files -L/path/to/lib/files  -lfftw3_mpi -lfftw3 -lm. Для версии 2.1.4 описание на русском языке приводится [https://parallel.ru/cluster/fftw.html здесь].
 
 
Остальные реализации:
 
*[http://www.gnu.org/software/gsl/ GSL - GNU Scientific Library]
 
*[https://developer.nvidia.com/cufft cuFFT]
 
*[http://www-03.ibm.com/systems/power/software/essl/ Engineering and Scientific Subroutine Library (ESSL) and Parallel ESSL]
 
= Литература =
 
[1] Бахвалов Н. С.,  Жидков Н. П., Кобельков. Г. М. — 6-е изд. — М. : БИНОМ. Лаборатория знаний
 
 
[2] Г. Нуссбаумер. Быстрое преобразование Фурье и алгоритмы вычисления сверток
 
 
[3] А.О. Лапаев Параллельное вычисление дискретного преобразования Фурье  полинома в простом поле
 

Текущая версия на 01:05, 15 октября 2016