#include #include #include /* Ne znam zasto ali koliko sam upoznat ova konstatna treba da je deo ansi C-a medjutim kada se kod kompajlira sa opcijom -ansi kompajler ne prepoznaje konstantu */ #ifndef M_PI #define M_PI 3.14159265358979323846 /* pi */ #endif /* Sluzi kao flag u f-ji shuffle da ukaze koju implementaciju funkcije da koristi pri prevodjenju */ #define SHUFFLE /* Ovu konstanu koristim da bih odredjivao frekvenciju kojom cu pozivati funkcije sinus i kosinus i biblioteke math.h Da bih smanjio gresku u racunanju sinusa i kosinusa preko formula za zbir sinusa i kosinusa. Broj koji se koristi za vrednost frekvencije je oblika 2^frekvencija. Posto u postavci zadatka nije bilo zahteva za odredjenom tacnoscu uzeo sam da je optimalno na svakih 128=2^7 racunanja potrebno ponovo izracunati sinus i kosinus. Kako bi anulirao gresku u racunanju */ #define FREKVENCIJA 7 /* Funkcija shuffle() transformise dati broj u broj koji se dobija kada se odredjen broj bitova sa desne strane procita naopacke. Parametri funkcije su: n - broj koji treba transformisati p - broj bitova koje treba procitati naopacke radi opisane transformacije Dati broj mora biti u opsegu [0,2^p), a broj bitova koji se citaju naopacke mora biti u opsegu [1,15]. */ static int shuffle (int n, int p) { /* Ovaj algoritam je brzi ali koristi promenjljivu vise */ #ifdef SHUFFLE register int rez; /* Rezultat funcije, cesto koriscena promeljiva pa je zbog optimizacije postavljeno kompajleru da ju postavi u registar procesora prilikom racunanja */ int i; /* Brojac u petlji */ /* Ispravnije bi bilo zbog portabilnosti koda koristiti sizeof(int) kao pokazatelj gornje granice, sa druge strane p je mogao da bude deklarisan kao unsigned short int */ assert (p >= 1 && p <= 15); /* Provera da li je "n" u opsegu [0,2^p] */ assert (0 <= n && n < (1 << p)); /* Deo bitova koji se ne menja */ rez = n >> p; for (i = 0; i < p; i++) { rez <<= 1; /* Pravim mesto za upisivanje sledeceg bita */ rez |= n & 1; /* Maska na poslednji bit i dodavanje promeljivoj rez */ n >>= 1; /* Pozicioniranje na sledeci bit */ } return rez; #endif /* Ovaj algoritam je sporiji ali ne koristi pomocnu promenljivu vec racunanje se vrsi u ulaznoj promenljivoj */ #ifndef SHUFFLE int i; /* Brojac u petlji */ assert (p >= 1 && p <= 15); assert (0 <= n && n < (1 << p)); /* Ukoliko je p parno sve je u redu, a ukoliko je neparno for ne petlja ne dolazi do bita koji se nalazi u sredini intervala pa njega ne treba ni menjati */ for (i = 0; i < p / 2; i++) /* XOR daje tacno kao rezultat akko je samo jedan od bitova tacan a drugi ne => bitovi na datim pozicijama su razliciti pa im treba zameniti mesta */ if (((n >> i) & 1) ^ ((n >> (p - i - 1)) & 1)) { /* Zamena vrednosti bitovima */ n ^= 1 << i; n ^= 1 << (p - i - 1); } return n; #endif } /* Funkcija complex_powers() izracunava stepene zadatog kompleksnog broja W sa jedinicnim modulom. Parametri funkcije su: theta - argument (fazni ugao) koji odredjuje kompleksan broj n - broj stepeni koje treba izracunati re, im - polja u koja ce biti smjesteni realne, odnosno imaginarne komponente stepena kompleksnog broja Broj stepeni koje treba izracunati mora biti veci od 0. */ static void complex_powers (double theta, int n, double *re, double *im) { short int epsilon; /* Promeljiva kojom podesavam frekvenciju pozivanja funkcija sinus i kosinus */ int i; /* Brojac u petlji */ assert (n > 0); /* Ukoliko je ukupan broj elemenata manji od frekvencije ponovnog racunanja nema potrebe za testiranjima pa onda podesavam epsilon na 0 */ epsilon = ((n >> FREKVENCIJA) > 0 ? FREKVENCIJA : 0); /* Racunam stepene kompleksnog broja pomucu formule za zbir sinusa i kosinusa,prednost ovog resenja je veca brzina izracunavanja nego da se pozivaju f-je sin,cos, a mana je smanjena racunska tacnost resenja */ re[0] = 1.0; /* theta^0 */ im[0] = 0.0; /* theta na nulti stepen */ /* Ukoliko trebam da racunam samo dva elementa u FFT-u da ne napravim prelivanje buffera */ if(n==1) return; /* Racunam kompleksni broj theta */ re[1] = cos (theta); im[1] = sin (theta); for (i = 2; i < n; i++) { /* Proveravam da li je "i" deljivo sa 2^epsilon( koristim shiftovanje u desno i levo jer su bitske operacije najbrze), tj. proveravam da li je ostatak deljenja jednak 0. Mogla bi se koristiti i finija podela frekvencije ako nam je podela 2^epsilon previse gruba tada bi sa operatorom % ispitivali dali je ostatak jednak nuli.Takodje mora da se proverava dali je epsilon razlicito od nule jer u suprotnom ako bi imali broj elemenata manji od frekvencije svaki put bi program racunao sin i cos koristeci standardne funkcije, takodje prvo se uporedjuje epsilon zbog optimizacije prilikom kompjaliranja i ako je on jednak 0 ostatak izraza se ne proverava */ if (epsilon && (i == ((i >> epsilon) << epsilon))) { re[i] = cos (i * theta); im[i] = sin (i * theta); } else { /* cos((n+1)*theta)=cos(n*theta)*cos(theta)-sin(n*theta)*sin(theta) */ re[i] = re[i - 1] * re[1] - im[i - 1] * im[1]; /* sin((n+1)*theta)=sin(n*theta)*cos(theta)+cos(n*theta)*sin(theta) */ im[i] = im[i - 1] * re[1] + re[i - 1] * im[1]; } } } void fft (int direct_fft, int p, double *re_in, double *im_in, double *re_out, double *im_out) { int duzinaVektora; /* Duzina vektora furijeve transformacije, koristim ju kao pomocnu promeljivu da ne bih svaki put shiftovao p ulevo */ double *wReal, *wImag; /* Niz realnih i imaginarnih stepena promenljive W, ovo su pomocne promenjljive koje koristim da bih smanjio visestruko racunanje stepena broja W */ /* Promenljive gornjiIndex, donjiIndex, stepen da bih izbegao visestruko racunanje u for petlji a sa druge strane povecavam citljivost koda */ int gornjiIndex; /* Gornji index leptira */ int donjiIndex; /* Donji index leptira */ int stepen; /* Stepen broja W */ double real, imag; /* Pomocne promeljive koje koristim da bih smestio medjurezultat koji racunam u ciklusu FFT-a */ int i, j, k; /* Brojaci u petljama */ /* Ispravnije bi bilo zbog portabilnosti koda koristiti sizeof(int) kao pokazatelj gornje granice, sa druge strane p je mogao da bude deklarisan kao unsigned short int */ assert (p >= 1 && p <= 15); /* Racunam duzinu vektora, shiftovanje u levo je ekvivalentno mnozenju broja sa dva dok je shiftovanje u desno ekvivalentno deljenju broja sa dva */ duzinaVektora = 1 << p; /* Pripremam vektor za racunanje tako sto obrcem indeke i kopiram vrednosti u izlazni vektor, i ne koristim dodatne pomocne promeljive vec racunaje direktno koristim izlazne rezultate */ for (i = 0; i < duzinaVektora; i++) { j = shuffle (i, p); /* Stedim jedan poziv f-je shuffle */ re_out[i] = re_in[j]; im_out[i] = im_in[j]; } /* Alokacija memorije za kompleksan broj W i za stepene broja W, izraz duzinaVektora>>1 je evivalentan izrazu duzinaVektora/2, shiftovanje koristim jer je brza operacija nego aritmeticko deljenje */ wReal = (double *) malloc ((duzinaVektora >> 1) * sizeof (double)); wImag = (double *) malloc ((duzinaVektora >> 1) * sizeof (double)); /* Racunam stepene broja W */ complex_powers (2 * M_PI / duzinaVektora, duzinaVektora >> 1, wReal, wImag); /* Ukoliko se radi o inverznoj FFT, flag direct_fft podesavam na vrednost -1 kako bi mogao racunam leptire uz pomoc tog flaga da ne pravim novu promeljivu */ if (direct_fft == 0) direct_fft = -1; /* Racunanje FFT-a, koristeci semu leptira a c \ _____ / | | | W^s | |_____| / \ b d c=a+b*W^s d=a-b*W^s Novo izracunate vrednosti c i d smestam nakon izracunavanja u vrednosti a i b jer, a i c imaju istu vrednost indeksa(gornjiIndex), a promenljive b i d imaju vrednost indexa donjIndex */ for (i = 1; i <= p; i++) /* Korak u FFT-u */ for (j = 0; j < duzinaVektora; j += 1 << i) /* Podskup leptira u koraku */ for (k = 0; k < (1 << (i - 1)); k++) /* Konretan clan iz podskupa leptira */ { gornjiIndex = j + k; donjiIndex = gornjiIndex + (1 << (i - 1)); stepen = (1 << (p - i)) * k; /* Ovde racunam b*W^s */ real = wReal[stepen] * re_out[donjiIndex] - direct_fft * wImag[stepen] * im_out[donjiIndex]; imag = wReal[stepen] * im_out[donjiIndex] + direct_fft * wImag[stepen] * re_out[donjiIndex]; /* Ovde racunam d=a-b*W^s */ re_out[donjiIndex] = re_out[gornjiIndex] - real; im_out[donjiIndex] = im_out[gornjiIndex] - imag; /* Ovde racunam c=a+b*W^s */ re_out[gornjiIndex] += real; im_out[gornjiIndex] += imag; } /* Oslobadjam dinamicki alociranu memoriju */ free (wReal); free (wImag); /* Ukoliko se radi o direktnoj brzoj furijevoj transformaciji potrebno je podeliti izlazni vektor sa duzinom vektora. Nazalost ovde ne mogu da shiftujem u desno */ if (direct_fft == 1) for (i = 0; i < duzinaVektora; i++) { re_out[i] /= duzinaVektora; im_out[i] /= duzinaVektora; } }