Анотація: У цій лекції розглянута теорія та реалізація алгоритму дискретного перетворення Фур'є в контексті паралельного програмування.

Нагадаємо коротко основні поняття і визначення, що відносяться до дискретного предобразованію Фур'є (ДПФ). Більш докладно про це див., Наприклад, у розділі 1 книги [ [4] ].

нехай є вектором з речовими або комплексними компонентами. Дискретним перетворенням Фур'є вектора називається вектор довжини з комплексними компонентами, які визначаються рівностями:

де і .

Обчислення перетворення Фур'є у вигляді (1), як воно записано вище, вимагає порядку умножений і складань.

У книзі [ [5] ], Розділ 32.3, наведено алгоритм обчислення ДПФ, що вимагає порядку операцій, і тому названий алгоритмом швидкого перетворення Фур'є (ШПФ). Однак, цей алгоритм важко піддається распараллеливанию, а тому, використовуючи його, важко отримати подальше зниження складності алгоритму БПФ.

Однак, існує інший алгоритм - так званий алгоритм Кулі-Тьюки (книга [ [4] ], Глава 4), який можна застосувати коли число є складовим. Зокрема, якщо для деяких і , То алгоритм Кулі-Тьюки вимагає порядку умножений, але який ефективно паралелі. Особливості паралельної реалізації алгоритму Кулі-Тьюки полягають у тому, що

1. вхідний і вихідний вектори розглядаються як двовимірні таблиці, стовпці і рядки яких можуть оброблятися незалежно в різних потоках,
2. в рамках одного потоку, обробка окремих стовпців і рядків є обчислення ДПФ за допомогою алгоритму БПФ.

Суть алгоритму Кулі-Тьюки (додаткові подробиці див. В книзі [ [4] ], Глава 4) складається в перетворенні виразу ~ (\ ref {travial}) відповідно до таких равенствами для індексів і :

Тоді, вираз (1) може бути переписано у вигляді

в котрому

У цьому виразі, при кожному значенні внутрішня сума являє собою -точковий перетворення Фур'є, а зовнішня сума є -точковий перетворення Фур'є. Відповідно, блок-схема послідовного алгоритму Кулі-Тьюки може бути представлена ​​наступним чином:

Зауважимо, що реального відображення вихідного вектора в двовимірний масив не відбувається - ця операція замінюється відповідним перерахунком індексів.

Нижче представлена ​​базова частина БПФ-алгоритму Кулі-Тьюки, яка включає в себе 4 етапи:

1. застосування БПФ до кожному колонку,
2. поелементне множення на ,
3. порядкове застосування БПФ,
4. відновлення результуючого вектора з двовимірної таблиці.

// Cooley-Tukey Algorithm public Complex [] FFT (Complex [] X, int n1, int n2) {int n = X.Length; Complex [,] Y = new Complex [n1, n2]; Complex [,] P = new Complex [n1, n2]; Complex [] T = new Complex [n]; for (int k = 0; k <n1; k ++) ColumnFFT (X, Y, n1, n2, k); for (int k = 0; k <n1; k ++) TwiddleFactorMult (Y, n1, n2, k); for (int i = 0; i <n2; i ++) RowFFT (Y, P, i, n1); // Sort elements in right order to create output vector for (int j = 0; j <n2; j ++) for (int i = 0; i <n1; i ++) T [i * n2 + j] = P [i, j]; return T; }

Тут, ColumnFFT і RowFFT - функції, які здійснюють БПФ для одного стовпчика і, відповідно, рядки матриці, функція TwiddleFactorMult - функція, що здійснює домноженіе стовпця матриці на додатковий множник, а функція TL_idx реалізує транспонування линеаризованной таблиці.

Паралельна реалізація ШПФ-алгоритму Кулі-Тьюки складається з двох наступних один за одним циклів Parallel .For, де в першому циклі виконуються кроки (1) і (2) послідовного алгоритму, а в другому циклі - кроки (3) і (4).

// Cooley-Tukey Algorithm. Parallel Implementation public Complex [] FFT (Complex [] X, int n1, int n2) {int n = X.Length; Complex [,] Y = new Complex [n1, n2]; Complex [,] P = new Complex [n1, n2]; Complex [] T = new Complex [n]; Parallel.For (0, n1, k => {ColumnFFT (X, Y, n1, n2, k); TwiddleFactorMult (Y, n1, n2, k);}); Parallel.For (0, n2, q => {RowFFT (Y, P, q, n1); // Sort elements in right order to create output vector for (int i = 0; i <n1; i ++) T [i * n2 + q] = P [i, q];}); return T; }

Повністю код паралельної версії БПФ-алгоритму Кулі-Тьюки наведено нижче:

// PFX Parallel Fast Fourier Transform (pFFT) // This implementation is based on the Cooley-Tukey algorithm using System; using System.Text; using System.Threading; public class Complex {public double Re = 0,0; public double Im = 0,0; public Complex () {} public Complex (double re, double im) {Re = re; Im = im; } Public override string ToString () {return Re + "" + Im; }} Class Program {public static void Main (string [] args) {if (args.Length! = 3) {Console.WriteLine (Usage ()); return; } // n - input vector length (must be power of two) // n1 - number of Cooley-Tukey's matrix columns // n2 - number of Cooley-Tukey's matrix rows int n = 0, n1 = 0, n2 = 0; n = (int) Math.Pow (2, Int32.Parse (args [0])); n1 = Int32.Parse (args [1]); // input vector generation Complex [] X = new Complex [n]; Random r = new Random (); for (int i = 0; i <n; i ++) X [i] = new Complex (r.NextDouble () * 10, 0); if ((n1> n) || (n1 <= 0)) {Console.WriteLine (Usage () + "Param n1 is invalid: n1 =" + n1.ToString () + ". Vector length =" + X. Length.ToString ()); return; } N2 = n / n1; Console.WriteLine ( "* RUN *"); DateTime dt1 = DateTime.Now; Complex [] pY = (new Program ()). FFT (X, n1, n2); DateTime dt2 = DateTime.Now; Console.WriteLine ( "Parallel FFT:"); Console.WriteLine ( "n =" + n + "n1 =" + n1 + "n2 =" + n2 + "Elapsed time is" + (dt2 - dt1) .TotalSeconds); } // Cooley-Tukey Algorithm. Parallel Implementation public Complex [] FFT (Complex [] X, int n1, int n2) {int n = X.Length; Complex [,] Y = new Complex [n1, n2]; Complex [,] P = new Complex [n1, n2]; Complex [] T = new Complex [n]; Parallel.For (0, n1, k => {ColumnFFT (X, Y, n1, n2, k); TwiddleFactorMult (Y, n1, n2, k);}); Parallel.For (0, n2, q => {RowFFT (Y, P, q, n1); // Sort elements in right order to create output vector for (int i = 0; i <n1; i ++) T [i * n2 + q] = P [i, q];}); return T; } Private void TwiddleFactorMult (Complex [,] Y, int n1, int n2, int k) {// Column Twiddle Factor Multiplication double wn_Re = 0, arg = 0, wn_Im = 0, tmp = 0; for (int q = 0; q <n2; q ++) {arg = 2 * Math.PI * k * q / (n2 * n1); wn_Re = Math.Cos (arg); wn_Im = Math.Sin (arg); tmp = Y [k, q] .Re * wn_Re - Y [k, q] .Im * wn_Im; Y [k, q] .Im = Y [k, q] .Re * wn_Im + Y [k, q] .Im * wn_Re; Y [k, q] .Re = tmp; }} Public void ColumnFFT (Complex [] a, Complex [,] A, int n1, int n2, int k) {int q, m, m2, s; double wn_Re, wn_Im, w_Re, w_Im; double arg, t_Re, t_Im; double u_Re, u_Im, tmp; int logN = 0; m = n2; while (m> 1) {m = m / 2; logN ++; } Int temp; for (q = 0; q <n2; q ++) {temp = bit_reverse (q, logN); A [k, temp] = a [k + q * n1]; } For (s = 1; s <= logN; s ++) {m = 1 <<s; arg = 2.0 * Math.PI / m; wn_Re = Math.Cos (arg); wn_Im = Math.Sin (arg); w_Re = 1.0; w_Im = 0,0; m2 = m>> 1; for (int i = 0; i <m2; i ++) {for (int j = i; j <n2; j + = m) {t_Re = w_Re * A [k, j + m2] .Re - w_Im * A [ k, j + m2] .Im; t_Im = w_Re * A [k, j + m2] .Im + w_Im * A [k, j + m2] .Re; u_Re = A [k, j] .Re; u_Im = A [k, j] .Im; A [k, j] .Re = u_Re + t_Re; A [k, j] .Im = u_Im + t_Im; A [k, j + m2] .Re = u_Re - t_Re; A [k, j + m2] .Im = u_Im - t_Im; } Tmp = w_Re * wn_Re - w_Im * wn_Im; w_Im = w_Re * wn_Im + w_Im * wn_Re; w_Re = tmp; }}} Public void RowFFT (Complex [,] a, Complex [,] A, int q, int n1) {int j, k, m, m2, s; double wn_Re, wn_Im, w_Re, w_Im; double arg, t_Re, t_Im; double u_Re, u_Im, tmp; int logN = 0; m = n1; while (m> 1) {m = m / 2; logN ++; } For (k = 0; k <n1; k ++) A [bit_reverse (k, logN), q] = a [k, q]; for (s = 1; s <= logN; s ++) {m = 1 <<s; arg = 2.0 * Math.PI / m; wn_Re = Math.Cos (arg); wn_Im = Math.Sin (arg); w_Re = 1.0; w_Im = 0,0; m2 = m>> 1; for (j = 0; j <m2; j ++) {for (k = j; k <= n1; k + = m) {t_Re = w_Re * A [k + m2, q] .Re - w_Im * A [k + m2, q] .Im; t_Im = w_Re * A [k + m2, q] .Im + w_Im * A [k + m2, q] .Re; u_Re = A [k, q] .Re; u_Im = A [k, q] .Im; A [k, q] .Re = u_Re + t_Re; A [k, q] .Im = u_Im + t_Im; A [k + m2, q] .Re = u_Re - t_Re; A [k + m2, q] .Im = u_Im - t_Im; } Tmp = w_Re * wn_Re - w_Im * wn_Im; w_Im = w_Re * wn_Im + w_Im * wn_Re; w_Re = tmp; }}} Public int bit_reverse (int k, int size) {int right_unit = 1; int left_unit = 1 <<(size - 1); int result = 0, bit; for (int i = 0; i <size; i ++) {bit = k & right_unit; if (bit! = 0) result = result | left_unit; right_unit = right_unit <<1; left_unit = left_unit>> 1; } Return (result); } Static string Usage () {string pname = System.Reflection.Assembly.GetEntryAssembly (). GetName (). Name; StringBuilder s = new StringBuilder (); for (int i = 0; i <Console.WindowWidth; i ++) s.Append ( "-"); return (s.ToString () + Environment.NewLine + "Usage:" + pname + "p n1 np" + Environment.NewLine + "where" + Environment.NewLine + "p - n = 2 ^ p - input vector length" + Environment.NewLine + "n1 - width of block" + Environment.NewLine + "np - number of processes" + Environment.NewLine + s.ToString () + Environment.NewLine); }}