|
|
| (не показаны 4 промежуточные версии этого же участника) |
| Строка 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] А.О. Лапаев Параллельное вычисление дискретного преобразования Фурье полинома в простом поле
| |