Програмна реалізація основних елементів ШПФ

Алгоритм попередньої перестановки

Розглянемо конкретну реалізацію ШПФ. Нехай є N=2^T елементів послідовності x_{N} і треба одержати послідовність X_{N}. Розділимо x_{N} на дві послідовності: парні і непарні елементи і т.д з кожною послідовністю. Ітераційний процес закінчиться, коли залишаться послідовності довжиною по 2 елементи. Приклад процесу для N=16 показаний нижче:

Разом виконується (log₂N)-1 ітерацій.

Розглянемо двійкове представлення номерів елементів і займаних ними місць. Елемент із номером 0 (0000) після всіх перестановок займає позицію 0 (0000), елемент 8 (1000) - позицію 1 (0001), елемент 4 (0100) - позицію 2 (0010), елемент 12 (1100) - позицію 3 (0011). І так далі. Позиції до і після перестановок дзеркально симетричні. Двійкове представлення кінцевої позиції виходить із двійкового представлення початкової позиції перестановкою бітів у зворотному порядку. І навпаки.

Це є закономірним для всіх N. На першому кроці парні елементи з номером n перемістилися в позицію n/2, а непарні з позиції в позицію N/2+(n-1)/2, де n=0,1,..., N-1. Таким чином, нова позиція обчислюється зі старої позиції як:

ror(n,N) = [n/2] + N{n/2},

де [x] - ціла частина числа, {x} - дробова.

В асемблері ця операція - циклічний зсув вправо (ror), якщо N - ступінь двійки. Спочатку береться двійкове представлення числа n, потім всі біти, крім молодшого зміщуються на 1 позицію вправо, а молодший біт переміщається на місце, що звільнилося, найстаршого (лівого) біта (див. рис.2.2).

Рис. 2.2

Подальші розбиття виконуються аналогічно. На кожному наступному кроці кількість послідовностей подвоюється, а число елементів у кожній з них зменшується вдвічі. Операції ror піддаються вже не всі біти, а тільки кілька молодших. Старші j-1 бітів залишаються недоторканими (зафіксованими), де j - номер кроку (див. рис.2.3):

Рис. 2.3

Простежимо за довільним бітом номера позиції. Нехай біт перебував в j -му двійковому розряді. Біт буде послідовно зсуватися вправо на кожному кроці доти, поки не виявиться в крайній правій позиції (після j -го кроку). На j+ 1 -му кроці буде зафіксовано j старших бітів, а той біт, за яким стежимо, переміститься в розряд з номером T-j-1, що і необхідно для дзеркальної перестановки біт. 

Алгоритм перестановки бітів наведений нижче:

for(I = 1; I < N-1; I++)

{

J = reverse(I,T); //reverse переставляє біти в I у зворотному порядку

if (I >= J)      // пропустити вже переставлені

conitnue;

S = x[I]; x[I] = x[J]; x[J] = S; // перестановка елементів xI і xJ

}

Операція зворотної перестановки біт може бути реалізована через інші бітові операції. Нижче наведений алгоритм функції перестановки T молодших бітів у числі I:

unsigned int reverse(unsigned int I, int T)

{

int Shift = T - 1;

unsigned int LowMask = 1;

unsigned int HighMask = 1 << Shift;

unsigned int R;

for(R = 0; Shift >= 0; LowMask <<= 1, HighMask >>= 1, Shift -= 2)

R |= ((I & LowMask) << Shift) | ((I & HighMask) >> Shift);

return R;

}

У змінних LowMask і HighMask зберігаються маски, що виділяють два біти, які переставляються. Перша маска у двійковому поданні виглядає як 0000...001 і в циклі змінюється, зсуваючи одиницю щораз на 1 розряд вліво:

0000...001

0000...010

0000...100

...

Маска HighMask за кожну ітерацію зсуває одиничний біт на 1 розряд вправо.

1000...000

0100...000

0010...000

...

Ці зсуви виконуються інструкціями LowMask <<= 1 і HighMask >>= 1.

Змінна Shift показує відстань (у розрядах) між бітами, що переставляються. Спочатку вона дорівнює T-1 і кожну ітерацію зменшується на 2. Цикл припиняється, коли відстань стає менше або дорівнює нулю.

Операція I & LowMask виділяє перший біт, потім він зсувається на місце другого (<<Shift). Операція I & HighMask виділяє другий біт, потім він зсувається на місце першого (>>Shift). Після чого обидва біти записуються в змінну R операцією "|".

Замість того щоб переставляти біти позицій місцями, можна застосувати і інший метод. Для цього треба вести відлік 0,1,2,...,N/ 2-1 уже зі зворотним проходженням бітів. Алгоритм збільшення на одиницю можна реалізувати програмно, приклад для T=4 наведений нижче.

I = 0;

J = 0;

for(J1 = 0; J1 < 2; J4++, J ^=1)

for(J2 = 0; J2 < 2; J3++, J ^=2)

for(J4 = 0; J4 < 2; J4++, J ^=4)

for(J8 = 0; J8 < 2; J8++, J ^=8)

{

if (I < J)

{

S = x[I]; x[I] = x[J]; x[J] = S; // перестановка елементів xI і xJ

}

I++;

}

У алгоритмі враховано, що при збільшенні числа від 0 до нескінченності (зі збільшенням на одиницю) кожний біт міняється з 0 на 1 і назад з певною періодичністю: молодший біт - щораз, наступний - кожний другий раз, наступний - кожний четвертий і так далі. Ця періодичність реалізована у вигляді T вкладених циклів, у кожному з яких один з бітів позиції J перемикається туди і назад за допомогою операції XOR (В C/C++ записується як ^=). Позиція I використовує звичайний інкремент I++.

Даний алгоритм залежить від кількості вкладених циклів залежно від T. Оскільки T звичайно обмежено деякою розумною межею (16..20), то можна написати декілька варіантів алгоритму. Нижче наведений варіант цього алгоритму, що емулює вкладені цикли через стеки Index і Mask.

int Index[MAX_T];

int Mask[MAX_T];

int R;

for(I = 0; I < T; I++)

{

Index[I] = 0;

Mask[I] = 1 << (T - I - 1);

}

 

J = 0;

for(I = 0; I < N; I++)

{

if (I < J)

{

S = x[I]; x[I] = x[J]; x[J] = S; // перестановка елементів xI і xJ

}

for(R = 0; R < T; R++)

{

J ^= Mask[R];

if (Index[R] ^= 1) // еквівалентно Index[R]^=1; if (Index[R]!= 0)

break;

}

}

Величина MAX_T визначає максимальне значення для T і в найгіршому випадку дорівнює розрядності цілочисельних змінних. Даний алгоритм повільнішим від попереднього, але економніший за обсягом коду.

І ще один алгоритм, який використовує класичний підхід до багаторозрядних бітових операцій: треба розділити 32-біта на 4 байти, виконати перестановку в кожному з них, після чого переставити самі байти.

Перестановку бітів в одному байті можна робити по таблиці, для якої потрібно мати масив reverse256 з 256 елементів (записані числа від 0 до 255 з переставленими в кожному порядок бітів).

Даний масив можна застосувати для реалізації функції reverse:

unsigned int reverse(unsigned int I, int T}

{

unsigned int R;

unsigned char *Ic = (unsigned char*) &I;

unsigned char *Rc = (unsigned char*) &R;

Rc[0] = reverse256[Ic[3]];

Rc[1] = reverse256[Ic[2]];

Rc[2] = reverse256[Ic[1]];

Rc[3] = reverse256[Ic[0]];

R >>= (32 - T);

Return R;

}

Звертання до масиву reverse256 переставляють у зворотному порядку біти в кожному байті. Покажчики Ic і Ir дозволяють звернутися до окремих байтів 32-бітних чисел I і R і переставити у зворотному порядку байти. Зсув числа R вправо наприкінці алгоритму усуває розходження між перестановкою 32 біт і перестановкою T біт. Геометрична ілюстрація алгоритму, де стрілками показані перестановки бітів, байтів і зсув наведена на рис.2.4.

Рис. 2.4

У всіх випадках є N перестановок з порівнянням I і J, що запобігає повторній перестановці деяких елементів.

Перевагу необхідно надати останньому алгоритму - він має усього N переходів.

Оптимізований останній алгоритм наведено нижче.

static unsigned char reverse256[]=

{

0x00, 0x80, 0x40, 0xС0, 0x20, 0xA0, 0x60, 0xE0,

 0x10, 0x90, 0x50, 0xD0, 0x30, 0xB0, 0x70, 0xF0,

0x08, 0x88, 0x48, 0xC8, 0x28, 0xA8, 0x68, 0xE8,

0x18, 0x98, 0x58, 0xD8, 0x38, 0xB8, 0x78, 0xF8,

0x04, 0x84, 0x44, 0xC4, 0x24, 0xA4, 0x64, 0xE4,

0x14, 0x94, 0x54, 0xD4, 0x34, 0xB4, 0x74, 0xF4,

0x0C, 0x8C, 0x4C, 0xCC, 0x2C, 0xAC, 0x6C, 0xEC,

0x1C, 0x9C, 0x5C, 0xDC, 0x3C, 0xBC, 0x7C, 0xFC,

0x02, 0x82, 0x42, 0xC2, 0x22, 0xA2, 0x62, 0xE2,

0x12, 0x92, 0x52, 0xD2, 0x32, 0xB2, 0x72, 0xF2,

0x0A, 0x8A, 0x4A, 0xCA, 0x2A, 0xAA, 0x6A, 0xEA,

0x1A, 0x9A, 0x5A, 0xDA, 0x3A, 0xBA, 0x7A, 0xFA,

0x06, 0x86, 0x46, 0xC6, 0x26, 0xA6, 0x66, 0xE6,

0x16, 0x96, 0x56, 0xD6, 0x36, 0xB6, 0x76, 0xF6,

0x0E, 0x8E, 0x4E, 0xCE, 0x2E, 0xAE, 0x6E, 0xEE,

0x1E, 0x9E, 0x5E, 0xDE, 0x3E, 0xBE, 0x7E, 0xFE,

0x01, 0x81, 0x41, 0xC1, 0x21, 0xA1, 0x61, 0xE1,

0x11, 0x91, 0x51, 0xD1, 0x31, 0xB1, 0x71, 0xF1,

0x09, 0x89, 0x49, 0xC9, 0x29, 0xA9, 0x69, 0xE9,

0x19, 0x99, 0x59, 0xD9, 0x39, 0xB9, 0x79, 0xF9,

0x05, 0x85, 0x45, 0xC5, 0x25, 0xA5, 0x65, 0xE5,

0x15, 0x95, 0x55, 0xD5, 0x35, 0xB5, 0x75, 0xF5,

0x0D, 0x8D, 0x4D, 0xCD, 0x2D, 0xAD, 0x6D, 0Xed,

0x1D, 0x9D, 0x5D, 0xDD, 0x3D, 0xBD, 0x7D, 0xFD,

0x03, 0x83, 0x43, 0xC3, 0x23, 0xA3, 0x63, 0xE3,

0x13, 0x93, 0x53, 0xD3, 0x33, 0xB3, 0x73, 0xF3,

0x0B, 0x8B, 0x4B, 0xCB, 0x2B, 0xAB, 0x6B, 0xEB,

0x1B, 0x9B, 0x5B, 0xDB, 0x3B, 0xBB, 0x7B, 0xFB,

0x07, 0x87, 0x47, 0xC7, 0x27, 0xA7, 0x67, 0xE7,

0x17, 0x97, 0x57, 0xD7, 0x37, 0xB7, 0x77, 0xF7,

0x0F, 0x8F, 0x4F, 0xCF, 0x2F, 0xAF, 0x6F, 0xEF,

0x1F, 0x9F, 0x5F, 0xDF, 0x3F, 0xBF, 0x7F, 0xFF,

}

 

unsigned int I, J;

unsigned char *Ic = (unsigned char*) &I;

unsigned char *Jc = (unsigned char*) &J;

 

for(I = 1; I < N - 1; I++)

{

Jc[0] = reverse256[Ic[3]];

Jc[1] = reverse256[Ic[2]];

Jc[2] = reverse256[Ic[1]];

Jc[3] = reverse256[Ic[0]];

J >>= (32 - T);

if (I < J)

{

 

   S = x[I];

   x[I] = x[J];

   x[J] = S;

}

}

Основний цикл алгоритму

На рис.2.5 наведений другий етап обчислення ДПФ. Лініями зверху вниз показано використання елементів для обчислення значень нових елементів. Зручно, що два елементи на певних позиціях в масиві дають два елементи на тих самих місцях. Тільки належати вони будуть вже іншим, більш довшим масивам, що розміщені на місці попередніх (коротших). Це дозволяє обійтись одним масивом даних для вихідних даних, результату і зберігання проміжних результатів для всіх T ітерацій.

Рис. 2.5

Нижче наведені дії, які необхідно виконати після первинної перестановки елементів.

#define T 4

#define Nmax (2 << T)

complex Q, Temp, j(0,1);

static complex x[Nmax];

static double pi2 = 2 * 3.1415926535897932384626433832795;

int N, Nd2, k, mpNd2, m;

for(N = 2, Nd2 = 1; N <= Nmax; Nd2 = N, N += N)

{

for(k = 0; k < Nd2; k++)

{

    W = exp(-j*pi2*k/N);

    for(m = k; m < Nmax; m += N)

    {

        mpNd2 = m + Nd2;

        Temp = W * x[mpNd2];

        x[mpNd2] = x[m] - Temp;

        x[m] = x[m] + Temp;

    }

}

}

Змінна Nmax містить повну довжину масиву даних і кінцеву кількість елементів ДПФ. В таблиці наведена відповідність між виразами в формулі (2.9) і змінними в програмі.

В алгоритмі	В формулі
N	N
Nd2	N/2
k	K
m	k з поправкою на номер послідовності
mpNd2	k+N/2 з поправкою на номер послідовності
pi2	2π
j	j (уявна одиниця)
x[k] (на початку циклу)	X{N/2}[even]k
x[mpNd2] (на початку циклу)	X{N/2}[odd]k
x[k] (в кінці циклу)	X{N}k
x[mpNd2] (в кінці циклу)	X{N}N/2+k
W

Зовнішній цикл - це основні ітерації. На кожній из них 2Nmax/N ДПФ (довжиною по N/2 елементів кожне) перетворюється в Nmax/N ДПФ (довжиною по N елементів кожне).

Наступний цикл по k є циклом синхронного обчислення елементів з індексами k і k + N/2 у всіх результуючих ДПФ.

Внутрішній цикл перебирає Nmax/N штук ДПФ одне за другим: спочатку обчислюються елементы з номерами 0 и N/2, у всіх ДПФ в кількості Nmax/N штук, потім з номерами 1 і N/2 + 1 і т.д. На рис.2.5 показана послідовність обчислень для N = 8, в якій забезпечується однократне обчислення .

x[mpNd2] обчислюється раніше від x[k], щоб x[k ] не було модифіковано перед використанням.

Алгоритм оберненого ДПФ відрізняється тим, що замість -j треба використовувати +j і після всіх обчислень поділити кожне x[k] на Nmax.

Оптимізація повертаючих множників

Для обчислення повертаючих множників можна використати формулу (10), що дозволить

.                                                           (2.10)

уникнути операцій піднесення до ступеня і обійтися одними множеннями, якщо заздалегідь обчислити W_N для N = 2, 4, 8,…,Nmax. У міру збільшення k повертаючий множник буде змінюватися, але разом з тим буде рости похибка обчислень. І після N/2 підряд множень у повертаючому множнику, може нагромадитися велике відхилення від точного значення, оскільки при множенні величин їхні відносні похибки сумуються.

Цього можна було б уникнути, будь число W_N цілим, але воно не ціле при N > 2, тому що обчислюється через значення синуса і косинуса:

Уникнути багатьох звертань до повільних функцій синуса і косинуса можна за допомогою алгоритму обчислення ступеня через багаторазові піднесення до квадрату і множення. Наприклад:

У цьому випадку необхідно всього 5 множень ( не потрібно обчислювати двічі) замість 13. У найгіршому випадку для піднесення до ступеня від 1 до N/ 2-1 потрібно log₂N множень замість N/2, що дає цілком прийнятну похибку для більшості практичних задач.

Можна вдвічі зменшити кількість множень на кожному кроці, якщо використовувати результати минулих обчислень

,

для зберігання яких потрібно додатково log₂(Nmax) комплексних комірок пам'яті:

Якщо чергове не було обчислено попередньо, то береться двійкове представлення k і аналізується. Кожному одиничному біту відповідає рівно один множник. У загальному випадку одиниці в біті з номером b відповідає множник , що зберігається в b-ій комірці масиву.

Для зменшення кількості множень при обчисленні до одного на два цикли необхідно відвести N/2 комплексних комірок для зберігання всіх попередніх . Непарні елементи обчислюються за формулою (2.10). Парні обчислюються за формулою , тобто, нічого не обчислюється, а береться одне зі значень, обчислених на попередньому кроці, коли N було вдвічі менше. Щоб не потрібно було копіювати величини на нові позиції досить їх відразу розташовувати в тій позиції, що вони займуть при N = Nmax і вводити просте виправлення Skew.

Пояснення до лістингу, що наведений нижче.

1. Реалізовані найпростіші операції над комплексними числами (класи Complex і ShortComplex), оптимізовані під дану задачу.

2. Масив W2n містить заздалегідь обчислені коефіцієнти W₂, W₄, W₈,...,W_Nmax.

3. Для обчислень використовуються найточніші представлення чисел із плаваючою крапкою в C++: long double розміром в 10 байт. Для зберігання результатів у масивах використовується тип double довжиною 8 байт.

/*

Fast Fourier Transformation

====================================================

*/

 

#ifndef FFT_H_

#define FFT_H_

 

struct Complex;

struct ShortComplex;

 

/*

Fast Fourier Transformation: direct (complement= false)

and complement (complement = true). 'x' is data source.

'x' contains 2^T items.

 

*/

 

extern void fft(ShortComplex *x, int T, bool complement);

struct ShortComplex

 

{

double re, im;

inline void operator=(const Complex &y);

};

struct Complex

{

long double re, im;

inline void operator= (const Complex &y);

inline void operator= (const ShortComplex &y);

};

 

inline void ShortComplex::operator=(const Complex &y) { re = (double)y.re; im = (double)y.im; }

inline void Complex::operator= (const Complex &y)  { re = y.re; im = y.im};

inline void Complex::operator= (const ShortComplex &y) { re = y.re; im = y.im};

#endif

 

 

Fft.cpp

/*

Fast Fourier Transformation

=====================================================

*/

#include "fft.h"

// This array contains values from 0 to 255 with reverse bit order

static unsigned char reverse256[]={

0x00, 0x80, 0x40, 0xС0, 0x20, 0xА0, 0x60, 0xE0,

0x10, 0x90, 0x50, 0xD0, 0x30, 0xB0, 0x70, 0xF0,

0x08, 0x88, 0x48, 0xС8, 0x28, 0xA8, 0x68, 0xE8,

0x18, 0x98, 0x58, 0xD8, 0x38, 0xB8, 0x78, 0xF8,

0x04, 0x84, 0x44, 0xC4, 0x24, 0xA4, 0x64, 0xE4,

0x14, 0x94, 0x54, 0xD4, 0x34, 0xB4, 0x74, 0xF4,

0x0C, 0x8C, 0x4C, 0xCC, 0x2C, 0xAC, 0x6C, 0xEC,

0x1C, 0x9C, 0x5C, 0xDC, 0x3C, 0xBC, 0x7C, 0xFC,

0x02, 0x82, 0x42, 0xC2, 0x22, 0xA2, 0x62, 0xE2,

0x12, 0x92, 0x52, 0xD2, 0x32, 0xB2, 0x72, 0xF2,

0x0A, 0x8A, 0x4A, 0xCA, 0x2A, 0xAA, 0x6A, 0xEA,

0x1A, 0x9A, 0x5A, 0xDA, 0x3A, 0xBA, 0x7A, 0xFA,

0x06, 0x86, 0x46, 0xC6, 0x26, 0xA6, 0x66, 0xE6,

0x16, 0x96, 0x56, 0xD6, 0x36, 0xB6, 0x76, 0xG6,

0x0E, 0x8E, 0x4E, 0xCE, 0x2E, 0xAE, 0x6E, 0xEE,

0x1E, 0x9E, 0x5E, 0xDE, 0x3E, 0xBE, 0x7E, 0xFE,

0x01, 0x81, 0x41, 0xC1, 0x21, 0xA1, 0x61, 0xE1,

0x11, 0x91, 0x51, 0xD1, 0x31, 0xB1, 0x71, 0xF1,

0x09, 0x89, 0x49, 0xC9, 0x29, 0xA9, 0x69, 0xE9,

0x19, 0x99, 0x59, 0xD9, 0x39, 0xB9, 0x79, 0xF9,

0x05, 0x85, 0x45, 0xC5, 0x25, 0xA5, 0x65, 0xE5,

0x15, 0x95, 0x55, 0xD5, 0x35, 0xB5, 0x75, 0xF5,

0x0D, 0x8D, 0x4D, 0xCD, 0x2D, 0xAD, 0x6D, 0xED,

0x1D, 0x9D, 0x5D, 0xDD, 0x3D, 0xBD, 0x7D, 0xFD,

0x03, 0x83, 0x43, 0xC3, 0x23, 0xA3, 0x63, 0xE3,

0x13, 0x93, 0x53, 0xD3, 0x33, 0xB3, 0x73, 0xF3,

0x0B, 0x8B, 0x4B, 0xCB, 0x2B, 0xAB, 0x6B, 0xEB,

0x1B, 0x9B, 0x5B, 0xDB, 0x3B, 0xBB, 0x7B, 0xFB,

0x07, 0x87, 0x47, 0xC7, 0x27, 0xA7, 0x67, 0xE7,

0x17, 0x97, 0x57, 0xD7, 0x37, 0xB7, 0x77, 0xF7,

0x0F, 0x8F, 0x4F, 0xCF, 0x2F, 0xAF, 0x6F, 0xEF,

0x1F, 0x9F, 0x5F, 0xDF, 0x3F, 0xBF, 0x7F, 0xFF,

};

 

void fft(ShortComplex *x, int T, bool complement

{

unsigned int I, J, Nmax, N, Nd2, k, m, mpNd2, Skew;

unsigned char *Ic = (unsigned char*) &I;

unsigned char *Jc = (unsigned char*) &J;

ShortComplex S;

ShortComplex *Wstore, *Warray;

Complex WN, W, Temp, *pWN;

Nmax = 1 << T;

//first interchanging

for(I = 1; I < Nmax - 1; I++)

{

   Jc[0] = reverse256[Ic[3]];

   Jc[1] = reverse256[Ic[2]];

   Jc[2] = reverse256[Ic[1]];

   Jc[3] = reverse256[Ic[0]];

   J >>= (32 - T);

   if (I < J)

   {

       S = x[I];

       x[I] = x[J];

       x[J] = S;

   }

}

//rotation multiplier array allocation

Wstore = new ShortComplex[Nmax / 2];

Wstore[0].re = 1.0;

Wstore[0].im = 0.0;

 

//main loop

for(N = 2, Nd2 = 1, pWN = W2n, Skew = Nmax >> 1; N <= Nmax; Nd2 = N, N+= N, pWN++, Skew >>= 1)

{

   //WN = W(1, N) = exp(-2*pi*j/N)

  WN= *pWN;

   if (complement)

       WN.im = -WN.im;

   for(Warray = Wstore, k = 0; k < Nd2; k++, Warray += Skew)

   {

       if (k & 1)

       {

           W *= WN;

           *Warray = W;

       }

       else

           W = *Warray;

       for(m = k; m < Nmax; m += N)

       {

           mpNd2 = m + Nd2;

           Temp = W;

           Temp *= x[mpNd2];

           x[mpNd2] = x[m];

           x[mpNd2] -= Temp;

           x[m] += Temp;

       }

   }

}

 

delete [] Wstore;

if (complement)

{

   for(I = 0; I < Nmax; I++)

       x[I] /= Nmax;

}

 

}