Algoritmul FFT

Se foloseste des termenul FFT = Fast Fourier Transform (Transformata Fourier Rapi­da). De fapt, este vorba despre un algoritm care foloseste proprietatile functiilor peri­o­di­ce pentru reducerea numarului de operatii aritmetice si nu este o transformata cu proprietati matematice diferite de cele ale transformatei DFT (Transformata Fourier Discreta).

Acest capitol reprezinta forma matriceala a transformatei Fourier discrete, o demon­stra­tie a algoritmului FFT si modul in care sunt organizate datele si calculele intr-un pro­gram care implementeaza algoritmul FFT.

Forma matriceala a transformatei Fourier discrete

Primii pasi in scrierea unui program care implementeaza transformata Fourier discreta sunt organizarea datelor sub forma matriceala si gasirea unei metode eficiente pentru calculul exponentialelor complexe.

Transformata DFT directa

Fie , o secventa prelevata din semnalul

, (3.1)

prin esantionarea cu perioada T. Transformata Fourier discreta transforma un semnal esan­tionat ,   intr-un spectru esantionat ,

. (3.2)

Daca se noteaza:

, (3.3)

unde se numeste nucleul transformatei Fourier discrete de dimensiune N. Cu notatia (3.3) esantioanele spectrului complex dat de transformata Fourier discreta sunt:

. (3.4)

Daca se scrie functia (3.4) pentru este indexul pentru linii iar n este in­de­xul pentru coloane, se obtine forma matriceala a transformatei Fourier discrete

, (3.5)

Pentru ca matricea transformatei Fourier discrete este simetrica fata de diago­nala principala. Functia vectoriala (3.5) se poate pune si sub forma:

, (3.6)

unde este vectorul esantioanelor semnalului analizat, este vectorul esanti­oanelor spectrului esantionat dat de transformata Fourier discreta iar matricea patrata , se numeste matricea transformatei DFT de dimensiune N.

Calculul exponentialelor

De exemplu pentru transformata DFT sub forma matriceala este :

. (3.7)

In cazul primei linii a matricei transformatei DFT k = 0 si . Prin inmultirea primei linii a matricei cu vectorul se obtine , care este valoarea me­die a sem­na­lului analizat.

Figura 3.1 Calculul exponentialelor

In figura 3.1 se prezinta cercul trigonometric si exponentialele , , care divi­zea­za cercul trigonometric in 8 puncte echidistante. Exponentialele se calculeaza la in­ce­putul programului si se memoreaza intr-un vector de numere complexe

. (3.8)

In cazul celei de-a doua linii a matricei k = 1 si cercul din figura 3.1 este parcurs o sin­gu­ra data. Prin inmultirea celei de-a doua linii a matricei transformatei cu vectorul se obtine , corespunzator frecventei fundamentale a transformatei Fourier discre­te de­fi­nita in sectiunea 2.2.3.

In cazul celei de-a treia linii k = 2, cercul din figura 3.1 este parcurs de doua ori, deci se calculeaza coeficientul Fourier c₂ corespunzator frecventei , armonica a doua a frecventei fundamentale a transformatei Fourier discrete. Exponentialele se iau tot din vectorul (3.8) pentru ca

. (3.9)

Calculul celorlalti coeficienti Fourier nu mai are importanta pentru ca nu are sem­ni­fi­catie fizica , , (v. sectiunea 3.2.4).

In cazul general, exponentialele se iau dintr‑un vector construit asemanator cu vectorul (3.8) folosind regula de indexare

. (3.10)

Transformata DFT inversa

Daca se cunoaste spectrul esantionat al unui semnal, un esantion al semnalului original se calculeaza cu formula:

, (3.11)

care este forma discreta a formulei (2.3). Similar ca in sectiunea 3.1.1 se noteaza:

(3.12)

si se obtine

, (3.13)

din care rezulta forma matriceala a transformatei Fourier inverse :

, (3.14)

unde si au aceeasi semnificatie ca (3.6), iar este matricea transformatei Fourier discreta inversa.

Fie o functie armonica cu perioada principala (v. definitiile din sectiunea 2.1.1).

. (3.15)

Armonicele acestei functii au urmatoarea proprietate

. (3.16)

Demonstratia este imediata: pentru , , iar pentru , integrala calculata pe un numar intreg de perioade dintr-o functie armonica este nula. Daca exponentialele complexe (3.16) sunt esantionate in N puncte si daca se folosesc notatiile (3.3) si (3.12) rezulta

. (3.17)

Demonstratia este similara. Pentru , deci suma da valoarea 1. Cercul din figura 3.1 este divizat uniform in N puncte simetrice fata de axa reala.

Pentru , sunt parcurse toate cele N puncte iar cercul este parcurs de ori folosind regula de indexare (3.10). In consecinta, partea imaginara a sumei (3.17) este nula pentru ca se aduna partile imaginare ale unor puncte simetrice fata de axa reala.

Daca N este par atunci punctele de pe cerc sunt simetrice si fata de axa imaginara. In consecinta si partea reala a sumei (3.17) este nula. Proprietatea (3.17) poate fi demonstrata si pentru N impar, dar demonstratia este laborioasa si rezultatul nu are importanta practica. De exemplu, pentru suma (3.17) da

. (3.18)

Utilizand (3.17), din (3.6) si (3.14) rezulta ca

. (3.19)

Algoritmul FFT

Transformata Fourier rapida (sau algoritmul FFT) este o metoda de calcul pentru o tran­sfor­mata DFT cu , esantioane ale semnalului , din care se obtin esantioane independente ale spectrului esantionat. Numarul natural se numeste dimensiunea transformatei FFT, iar numarul este numarul de iteratii necesare in algoritmul FFT.

Principiul algoritmului FFT

In sens matematic, algoritmul este o metoda sistematica de rezolvare a unei categorii de probleme. In sens informatic, algoritmul este un set de instructiuni prestabilit prin care se rezolva o anumita problema intr‑un numar finit de pasi.

Algoritmul FFT are h iteratii, . In continuare se prezinta o varianta a algo­rit­mului FFT in care se parcurg urmatorii pasi:

iteratia 0 se aplica vectorului de numere complexe de dimen­si­u­ne , indexat dupa regula, , din care se obtin doua tran­sfor­ma­te Fourier rapide de dimensiune . Aceste rezultate intermediare sunt depuse in vectorul ;

iteratia 1 a algoritmului se aplica separat celor doua jumatati ale vectorului , indexate dupa regulile si respectiv din care se obtin patru transformate Fourier discrete de dimensiune ;

algoritmul se incheie cu iteratia , care se aplica celor 2^h transformate Fourier Rapide de dimensiune din vectorul din care se obtin esan­ti­oane de frecventa de forma .

Observatie: in iteratiile , vectorii de numere complexe nu au sem­ni­­fi­catie fizica si din acest motiv, vor fi numiti "vectorii variabilelor intermediare" ;

Demonstratia algoritmului este in stil constructivist. In principiu, daca se demonstreaza ca fiecare iteratie se dubleaza numarul de transformate Fourier rapide si se injumatateste dimen­si­u­nea lor si daca se stabileste semnificatia fizica a numerelor complexe din vectorul , atunci algoritmul este demonstrat.

Demonstratia algoritmului FFT

Se porneste de la formula transformatei Fourier discrete (3.4) pusa sub forma

, (3.20)

unde indicii zero din si indica numarul iteratiei curente, iar este indexul sumei din iteratia curenta. La fiecare iteratie numerele complexe din vectorul , , isi schimba pozitiile datorita metodei de calcul prezentata in continuare.

Se divizeaza secventa in doua jumatati

. (3.21)

In a doua suma se face substitutia

, (3.22)

unde se tine seama ca formula (2.20) a fost dedusa in cazul primei ipoteze din sectiunea 2.2.2, deci in cazul indexului si a limitei de indexare calculele se fac intr-o clasa de resturi modulo

. (3.23)

Dupa ce se desface paranteza de la exponentul celei de-a doua sume din (3.22) si se regrupeaza termenii asemenea se obtine

, (3.24)

unde .

Se considera separat cazurile k par si k impar. Se noteaza si se obtine

, (3.25)

iar pentru k impar se noteaza si se obtine

,

. (3.26)

Se observa ca paranteza rotunda din (3.25) si paranteza patrata din (3.26) sunt niste nu­me­re complexe usor de evaluat. Se fac notatiile:

, . (3.27)

cu notatiile (3.27) se obtin doua transformate Fourier rapide de dimensiune :

, (3.28)

, (3.29)

pentru ca din (3.3) rezulta ca

. (3.30)

Raman de rezolvat cateva probleme legate de organizarea datelor si calculelor. In (3.27) se observa variabilele intermediare si pot fi calculate simultan din si si ca rezultatele raman "in situ".

In (3.28) se observa ca variabilele intermediare , care intra in for­mu­la de calcul a esantioanelor pare de frecventa sunt ordonate in sens crescator in prima jumatate a vectorului. Daca se fa­ce substitutia in (3.29), ti­nand cont de (3.23), se obtine formula de calcul a esantioanelor impare de frecventa

, , (3.31)

unde variabilele intermediare , sunt ordonate in sens crescator in a doua jumatate a vectorului.

Formulele (3.28), (3.30) si (3.31) demonstreaza ca o iteratie a algoritmului dubleaza nu­ma­rul de transformate Fourier rapide si le injumatateste dimensiunea si ca sunt necesare iteratii pentru aflarea valorii esantioanelor de frecventa.

Ordonarea esantioanelor de frecventa

Din (3.28) si (3.31) rezulta dupa prima iteratie esantioanele de frecventa sunt ames­te­cate. Acest proces de amestecare continua si in ite­ratiile urmatoare deci demonstratia din sectiunea anterioara tre­bu­ie completata cu un algoritm de ordonare finala a esantioa­ne­lor de frecventa.

In practica nu se memoreaza vectorul esantioanelor de frecventa pentru ca, dupa ultima iteratie se obtine . Este suficient sa se gaseasca o functie de forma:

, , (3.32)

care sa indice pozitiile esan­ti­oa­nelor de frecventa.

In sectiunea precedenta s‑a aratat ca dupa fiecare iteratie a unei transformate Fourier rapide esantioanele de frecventa sunt sortate: intai cele pare in prima jumatate a vectoru­lui, apoi cele impare sunt asezate in ordine crescatoare in a doua jumatate a vectorului.

In tabelul figura 3.2 se da un exemplu pentru . In prima coloana se gaseste n indexul pentru vectorul , in a doua coloana se gaseste scris cu cifre binare iar in a doua coloana se gaseste scris tot cu cifre binare.

Figura 3.2. Metoda de calcul a indecsilor ,

Ordonarea descrisa anterior consta in luarea primei cifre binare din fiecare numar si asezarea ei la sfarsitul numarului. In sectiunea 3.2.1 s‑a aratat ca dupa prima iteratie se obtin doua transformate Fourier rapide de dimensiune . Coloana a doua a tabelului contine numerele . Primele trei cifre binare din numerele sunt indexul variabilei intermediare care intra in transformata Fou­ri­er rapida de dimensiune .

In coloana a treia a tabelului se gasesc numerele . Numerele se obtin prin luarea primei cifre binare din grupul de trei cifre binare si asezarea ei in ultima pozitie a gru­pu­lui. Dupa doua iteratii se obtin patru transformate Fourier rapide de dimensiune . Primele doua cifre binare din numerele sunt indexul variabilei intermediare care intra in transformata Fou­ri­er rapida de dimensiune .

Dupa trei iteratii are cifrele binare ale numarului asezate in ordine inversa. Daca numerele si se scriu cu cifre binare atunci pozitia finala a esantionului de frecventa se afla inversand ordinea cifrelor binare si functia (3.32) devine:

, . (3.33)

Programarea algoritmului FFT

Pentru scrierea unui program care sa implementeze algoritmul FFT, demonstratia pre­zen­tata anterior trebuie completata prin precizarea urmatoarelor probleme:

care este memoria necesara pentru rularea unui program de dimensiune ;

cum se calculeaza simultan variabilele intermediare cu formulele (3.27) astfel incat calculele sa se faca "in situ";

cum se calculeaza indecsii necesari pentru parcurgerea vectorului variabilelor in­ter­mediare si vectorul exponentialelor in iteratiile stiind ca la fiecare ite­ra­tie dimensiunea transformatei se injumatateste si se dubleaza numarul de tran­sfor­mate Fourier rapide;

cum se face ordonarea finala a esantioanelor de frecventa tinand seama de ope­ra­tiile de sortare a vectorului variabilelor intermediare (3.29) si (3.31).

Organizarea datelor

In memoria calculatorului se rezerva:

, , un vector de numere complexe folosit pen­tru memorarea variabilelor inter­me­­di­are;

, , un vector pentru memorarea expo­nen­tialelor complexe;

, , un vector de numere intregi fara semn in care se memo­rea­za functia (3.33) tabelata. In sectiunea 3.3.4 se prezinta algoritmul prin care se gene­rea­za acest vector.

Figura 3.3. Alocarea memoriei

In figura 3.3 se prezinta declaratiile prin care se aloca memoria necesara pentru vectorii , si . Daca se face alocare statica a memoriei, de regula se declara spatii mai mari decat cele necesare in rularea respectiva. Daca se face alocare dinamica atunci, in timpul rularii programului se pot aloca exact spatiile necesare.

Vectorii si sunt formati din obiecte din clasa . Clasa are doua cam­puri si si metode care implementeaza operatiile de adunare, scadere, inmul­ti­re, calculul modulului si a fazei. In sectiunea 3.3.2 se va prezenta modul in care se folosesc metodele clasei .

La pornirea algoritmului FFT, , in partea reala a numerelor complexe din vectorul varia­bilelor intermediare se incarca in ordine esanti­oa­ne­le semnalului analizat (v. for­mu­la (3.1)) iar partea imaginara este anulata

. (3.34)

In iteratiile, numerele din vectorul nu au semnificatie fizica. La in­che­­ierea algoritmului, , vectorul contine valorile esantioa­ne­lor de frec­ven­ta care sunt amestecate datorita ordonarilor care se fac in fiecare iteratie.

La pornirea algoritmului vectorul exponentialelor complexe se initializeaza cu formula

. (3.35)

(v. notatia (3.3) din sectiunea 3.1.1). In figura 3.4 se prezinta procedura de initializare a vectorului exponentialelor complexe.

Figura 3.4. Initializarea vectorului exponentialelor complexe

Din formula (3.26) rezulta ca in prima iteratie sunt necesare doar , iar in fiecare dintre iteratiile urmatoare numarul exponentialelor complexe folosite in calcule se injumatateste dar exponentialele se iau tot din vectorul . In sectiunea 3.3.3 se va deduce o regula pentru calculul indexului cu care se adre­sea­za vectorul exponentialelor complexe.

Organizarea calculelor

Un graf de fluenta se compune din puncte si arce orientate. Un punct din care pornesc unul sau mai multe arce este un operand initial. Un punct in care converge cel putin un arc si din care porneste cel putin un arc este un operand intermediar. Un punct din care nu porneste nici un arc este un rezultat. Un nume scris langa un punct desemneaza o da­ta initiala, o variabila intermediara sau un rezultat.

Daca intr‑un punct converg mai multe arce atunci in punctul respectiv se face o insu­ma­re. Daca intr‑un punct converge un singur arc atunci se face doar operatia de atribuire. Numarul scris in dreptul unui arc este ponderea cu care participa operandul din punctul initial in operatia din punctul final.

Figura 3.5. Graful de fluenta al unui fluture

Graful de fluenta din figura 3.5 care implementeaza calculele din formulele (3.27) poar­ta numele sugestiv de fluture. Un fluture calculeaza "in situ" variabilele intermediare si din variabilele intermediare si . In figura 3.6 se prezinta procedura care implementeaza calculele dintr‑un fluture.

Figura 3.6. Implementarea unui fluture

Clasa este organizata ca o biblioteca aritmetica cu un singur operand. Cam­pu­ri­le si au rolul de acumulator. Metodele , si au ca parametru tot un numar complex transmis prin valoare si implementeaza operatiile aritmetice de adunare, scadere si respectiv inmultire a numerelor complexe. Rezultatul operatiei ra­ma­ne in acumulator.

Indecsi

In figurile 3.5 si 3.6 este indexul pentru prima linie a fluturelui, indexul pentru a doua linie a fluturelui iar este indexul pentru exponentiala complexa. Calculul indecsilor , si va fi detaliat in sectiunea urmatoare.

Pentru unificarea notatiilor, in afara indexului pentru iteratii, , se introduc doi indecsi noi ale caror limite maxime de indexare se exprima in functie de :

pentru fluturii dintr‑o transformata Fourier rapida;

pentru grupurile de fluturi. Un grup de fluturi corespunde cu o transfor­mata Fourier rapida.

Figura 3.7. Generarea indecsilor in algoritmul FFT

In figura 3.7 se prezinta rutina care genereaza indecsii , si , indexul folosit pentru adresarea vectorului exponentialelor complexe. Se observa ca sunt trei bucle controlate de indecsii , indexul pentru iteratii, , indexul pentru grupurile de fluturi si , indexul pentru fluturi.

In prima iteratie exista un singur grup de fluturi. Regula de indexare rezulta din formulele (3.27) , . In iteratia urmatoare si lucrurile se complica pentru ca sunt doua transformate Fourier Rapide de dimensiune deci sunt doua grupuri de cate fluturi.

In cazul general, fiecare transformata Fourier rapida de dimensiune are cate un grup de flu­turi, iar a doua linie a unui fluture se afla la distanta de fata de prima linie. Indecsii , si se calcu­lea­za cu urmatoarele formule:

   , , (3.36)

In prima iteratie vectorul este parcurs in ordine complet. La injumatatirea dimensi­u­nii tran­sfor­matei Fourier rapide se injumatateste si dimensi­u­nea vectorului exponen­ti­a­le­lor complexe, dar exponentialele divizeaza uniform tot o jumatate de cerc. Din (3.30) rezulta ca , deci vectorul este parcurs din doua in do­ua valori, ceea ce de­monstreaza formula .

La scrierea unui program care implementeaza algoritmul FFT cea mai complicata pro­ble­ma este calculul corect al indecsilor , si . In aceasta fata de dezvoltare a programului procedurile si nu fac decat sa scrie intr‑un fisier text indec­sii , , , , si sub forma unui tabel. Dupa verificarea corectitudinii indecsilor se inlatura apelul procedurii iar procedura va avea continutul din figura 3.6.

Algoritmul de ordonare

Fie un numar intreg fara semn oarecare, cu cifre binare . Algoritmul care implementeaza func­tia (3.3.3) este complicat pentru ca nu exista o succesiune sim­pla de in­struc­tiuni scrise in limbaj de asamblare sau in limbaj de nivel inalt care sa in­ver­­se­ze ordinea cifrelor binare.

Din acest motiv in practica fo­lo­seste un algoritm iterativ care genereaza direct nume­rele binare cu cifrele asezate in or­di­ne inversa. In memorie se declara vectorul , si se initializeaza , . Apoi, tabelul func­tiei (3.3.3) se genereaza cu formulele iterative:

   , , (3.37)

Se observa ca la fiecare iteratie algoritmul dubleaza lungimea vectorului si ca algoritmul poate continua pana la orice dimensiune. Prima dintre formulele 3.37 adauga o cifra binara 0 la sfarsitul numarului iar a doua formula adauga o cifra binara 1 la sfarsitul nu­ma­rului.

In figura 3.5 se prezinta algoritmul de generare a numerelor binare asezate in ordine in­ver­sa pentru , deci sunt necesare patru iteratii indexate .

Figura 3.8. Algoritmul de generare a numerelor binare asezate in ordine inversa

In prima coloana a tabelului se gaseste indexul , in coloanele urmatoare se gasesc numerele , , obtinute dupa fiecare iteratie iar in ultima coloana se ga­seste scris cu cifre zecimale. Functia discreta (3.3.3) de forma da po­zi­tiile adevarate a esantioanelor de frecventa. Prima si ultima coloana din tabel sunt func­tia (3.3.3) tabelata.

Figura 3.9. Procedura care genereaza vectorul

In figura 3.9 se prezinta procedura care genereaza vectorul de di­men­siune apelul procedurii are ca efect initializarea cu 0 a vectorului .

Figura 3.10. Procedura de ordonare a esantioanelor de frecventa

In figura 3.2 se prezinta procedura care ordoneaza esantioanele de frecventa dupa in­che­ie­rea rularii algoritmului FFT. Din tabelele din figurile 3.2 si 3.8 rezulta ca sunt cateva esan­­tioane de frecventa care iti pastreaza pozitia. In celelalte cazuri esantioanele de frec­­­venta cu pozitiile si trebuie schimbate intre ele.

La parcurgerea com­pleta a vectorului perechea de indecsi si este intalnita de doua ori. Esan­tioanele de frecventa si se sunt schimbate intre ele da­ca . In acest fel se evita schimbarea esantioanelor de frecventa care isi pas­trea­za pozitia si se evita ca esan­tioanele de frecventa si sa fie schim­ba­te intre ele de doua ori.

Concluzii

Din acest capitol se desprind doar doua concluzii in comparatie cu transformata Fourier discreta algoritmul FFT este eficient dar este complicat.

Eficienta calculelor

Daca calculele se fac cu transformata Fourier discreta pusa sub forma matriceala (v. (3.5), sectiunea 3.1.1) atunci pentru calculul unui esantion de frecventa sunt necesare inmultiri si adunari. In consecinta, pentru calculul celor esantioane de frec­ven­ta sunt necesare inmultiri si adunari.

Fiecare transformata Fourier rapida de dimensiune are cate un grup de flu­turi. In iteratia respectiva sunt transformate Fourier rapide, deci in fiecare iteratie se calculeaza fluturi.

Pentru calculul unui fluture sunt necesare doua adunari si o inmultire, deci intr‑o iteratie sunt necesare adunari si inmultiri. Algoritmul are iteratii deci in total sunt necesare adunari si inmultiri.

Figura 3.11. Eficienta algoritmului FFT

In tabelul din figura 3.11 sunt comparate numarul de operatii aritmetice necesare pentru tran­sfor­mata Fourier discreta de dimensiune si pentru algoritmul FFT. Se ob­ser­va ca algoritmul reduce mult numarul operatiilor aritmetice necesare pentru calculul spec­­trului discret.

De ce FFT este un algoritm complicat?

Argumentatie:

in subcapitolul 3.4.1 se precizeaza notiunile si notatiile cu care opereaza tran­sfor­mata Fou­ri­er Discreta. Subcapitolul Sectiunea 3.4.1 este simpla in com­pa­ratie cu cele urmatoare;

in subcapitolul 3.4.2 se prezinta demonstratia unei variante a algoritmului FFT cu­­nos­cu­ta sub nu­mele decimation‑in‑frequency FFT Algorithm. In de­mon­stra­tie in­ter­vin multe schimbari de variabila si multe notatii;

in subcapitolul 3.4.3 se definesc structurile de date ale programului care imple­men­teaza algoritmul FFT si modul de organizare a calculelor (v. sectiunea 3.3.2). Se folosesc trei indecsi independenti , si si inca trei indecsi , si calculati din primii trei cu formulele (3.37).

Din aceste motive, in timp, s‑au incercat mai multe variante de demonstrare si de im­ple­­men­tare ale algoritmului FFT. Toate eforturile s‑au soldat cu acelasi rezultat: com­ple­xi­ta­tea al­go­ritmului se pastreaza in oricare dintre variantele lui.

De exemplu, in algoritmul cunoscut sub nu­mele decimation‑in‑time FFT Algorithm, se demonstreaza ca prin agre­ga­rea a doua transformate Fourier rapide de dimensiune N, se obtine o tran­sfor­mata Fourier rapida de dimensiune . Demonstratia decurge in sens invers fata de ceea prezentata in subcapitolul 3.4.2. Pentru ca toate functiile care apar in implementarea unei transformate DFT discrete sunt liniare, atunci este evident ca aceste functii sunt inversabile si ca algoritmul poate redactat in ordine inversa. Deci rezultatele sunt reversibile.

Bibliografie

Popescu D. Prelucrarea digitala a semnalelor. Editura ICPE Bucuresti 2000.

Pop E., I. Nafornita, V. Tiponut, A. Mihaescu, L. Toma. Metode in prelucrarea sem­na­le­lor. Editura Facla, Timisoara, 1986.

Proakis J.G., D. G. Manolakis. Digiral Signal Processing, Principles, Aghorithms and Applications. Thire Edittion. Prince / Hall International, Inc., 1996.

Oppenhiem A.V., Schaffer R. V. Digiral Signal Processing. Prince / Hall International, Inc.

Gade S., H. Herlufsen: Use of Weighting Functions in DFT / FFT Analysis (Part I).

Brüel & Kjær Teckhnical Rewiew, No. 3 1997;