Kuli-Tukijev algoritam

Kuli-Tukijev algoritam (енгл. Cooley-Tukey algorithm) je najčešće korišćen algoritam za izračunavanje brze Furijeove transformacije (engl. The Fast Fourier Transformation). Algoritam je prvi put objavljen 1969. od strane Džejmsa Kulija i Džona Tukija u članku “Brza Furijeova transformacija I njene primene'” (engl. The Fast Fourier Transformationa and its Applications). Algoritam je zasnovan na principu “zavadi pa vladaj” (engl. divide–and–conquer) i funkcioniše tako što izračunava diskretnu Furijeovu transformaciju cele ulazne sekvence koristeći diskretnu Furijeovu transformaciju podsekvenci ulazne sekvence.

Istorija Cooley-Tukey algoritma

Ovaj algoritam, uključujući njegovu rekurzivnu primenu, izmislio je Karl Fridrih Gaus oko 1805. godine, koji ga je koristio za interpolaciju trajektorija asteroida Palas i Juno, ali njegov rad nije bio široko prepoznatljiv (objavljen je samo posthumno i na neo-latinskom jeziku). ^[1]^[2] Međutim, Gaus nije analizirao asimptotsku vremensku složenost algoritma. Razne ograničene forme su takođe otkrivene nekoliko puta tokom 19. i početka 20. veka. FFT je postao popularan nakon što su Džejms Kuli iz IBM-a i Džon Tjuki sa Prinstona objavili je rad 1965. godine u kome je algoritam detaljno opisan, kao i postupak kako ga praktično izvoditi na računaru. ^[3]

Tjuki je, navodno, došao do ideje o algoritmu tokom sastanka Savetodavnog Odbora za naučne studije predsednika Kenedija o načinu otkrivanja testova nuklearnog oružja u Sovjetskom Savezu koristeći seizmometre locirane izvan zemlje. Ovi senzori bi generisali seizmološke vremenske serije. Međutim, analiza ovih podataka zahtevala bi brze algoritme za izračunavanje DFT zbog broja senzora i dužine vremena. Ovaj zadatak je bio kritičan za ratifikaciju predložene zabrane nuklearnog testiranja, tako da bi se bilo kakva kršenja mogla otkriti bez potrebe za posetom sovjetskih objekata. ^[4]^[5] Još jedan učesnik na tom sastanku, Ričard Garvin iz IBM-a, prepoznao je potencijal metoda i stavio Tjukija u kontakt sa Kulijem, iako je bio siguran da Kuli nije znao prvobitnu svrhu. Umesto toga, Kuliju je rečeno da je to potrebno za određivanje periodičnosti spin orijentacija u 3-D kristalu Helijum-3. Kuli i Tjuki su kasnije objavili svoj zajednički članak, a široko usvajanje je brzo usledilo zbog istovremenog razvoja Analogno-digitalnih pretvarača sposobnih za uzorkovanje brzinom do 300 kHz.

Činjenica da je Gaus opisao isti algoritam (iako bez analize njegovog asimptotskog vremenskog troška) nije prepoznata do nekoliko godina nakon izdanja rada Kulija i Tjukija 1965-e. Njihov rad je citiran kao inspiracija pri radu I.J. Guda na onome što se danas naziva FFT zasnovan na prostim činiocima (eng. Prime-factor FFT algoritam, PFA); iako je algoritam Guda prvobitno smatran ekvivalentom Kuli-Tjuki algoritmu, brzo se shvatilo da je PFA sasvim drugačiji algoritam (radi samo za veličine koje imaju uzajamno proste faktore i oslanja se na Kinesku teoremu o ostacima, za razliku od Kuli-Tjukijeve koja radi sa bilo kojim veličinama). ^[2] Their paper cited as inspiration only the work by I. J. Good on what is now called the prime-factor FFT algorithm (PFA);^[3] although Good's algorithm was initially thought to be equivalent to the Cooley–Tukey algorithm, it was quickly realized that PFA is a quite different algorithm (only working for sizes that have relatively prime factors and relying on the Chinese Remainder Theorem, unlike the support for any composite size in Cooley–Tukey).^[6]

O algoritmima za brzu Furijeovu transformaciju

Algoritmi za izračunavanje brze Furijeove transformacije se koriste za izračunavanje Furijeove transformacije za datu ulaznu sekvencu. Diskretna furijeova transformacija transformiše diskretni impulsni signal iz vremenskog domena u domen frekvencije. Brza Furijeova transformacija se koristi u mnogo različitih sfera, ukjučujući digitalnu obradu signala, telekomunikacije, i analizu zvučnog signala. Vremenska složenost algoritma brze Furijeove transformacije na standardnim mikroprocesorskim sistemima sa fon Nojmanovom arhitekturom je $O(NlogN)$ . Ubrzana varijanta algoritma radi u vremenskoj složenosti $O(logN)$ . Iako je vremenska složenost ubrzanog algoritma značajno manja nego složenost standardnog algoritma, treba imati u vidu činjenicu da aplikacije za ubrzanje koriste Maxelex DFE tehnologiju, koja ima neke razlike u poređenju sa standardnim fon Nojmanovim modelom programiranja.^{[појаснити]} Takođe, treba uzeti u obzir da je, prilikom korišćenja Maxeler DFE tehnologije za ubrzanje aplikacije, potrebno određeno vreme za prenos podataka na i iz DFE-a. Sa druge strane, Maxeler DFE sistemi imaju prednost da mogu da izračunavaju rezultate pri svakom ciklusu sata, nakon perioda pocetne latentnosti, stvaranjem hardverskih puteva. Ova osobina omogućuje adekvatnim aplikacijama da postignu znacajna ubrzanja koriscenjem Maxeler DFE sistema.

Pregled različitih algoritama za brzu Furijeovu transformaciju

Postoji mnogo drugih algoritama koji računaju brzu Furijeovu transformaciju (FFT), a koji se razlikuju od Kuli-Tjuki algoritma, kao što su: Brunov FFT algoritam, Riderov FFT algoritam, Vinogradov FFT algoritam, FFT algoritam zasnovan na prostim činiocima, I Blustajnov FFT algoritam. Ovi algoritmi koriste različite matematičke metode da izračunaju brzu Furijeovu transformaciju. Korišćene matematičke metode variraju od teorije brojeva, preko numeričke matematike, do teorije grafova. Svi ovi algoritmi imaju jednu zajedničku osobinu: da im je vremenska složenost $O(NlogN)$ kada su pokrenuti na procesoru. Iako ne postoje dokazi da je nemoguće konstruisati brzi algoritam sa manjom vremenskom složenošću od $O(NlogN)$ , takav algoritam ipak do sada još nije konstruisan.

Upoređivanje performansi velikog broja algoritama koji računaju brzu Furijeovu transformaciju je napravljeno od strane Matea Frigoa i Stivena G. Džonsona sa Masačusetskog tehnološkog instituta u SAD-u. U njihovim istraživanjima korišćene su implementacije različitih autora, napisane u periodu dužem od 35 godina. Eksperimenti su izvrsavani na različitim arhitekturama racunara. Treba naglasiti da su Frigo i Džonson u eksperimentima koristili samo jednoprocesorsko jezgro u multiprocesorskim sistemima. Slika 1 je grafička reprezentacija podataka sa Frigovog i Džonsonovog eksperimenta, i prikazuje poređenje između performasi različitih algoritama za računanje brze Furijeove transformacije za kompleksnu ulaznu sekvencu u zavisnosti od dužine ulazne sekvence. Eksperiment je izveden na Intel Xeon 3.60 GHz Pentium 4 (Prescott) mikroprocesoru. Vrednosti na Y-osi pokazuju brzinu algoritama predstavljenje u Mflops-ima. Vrednost Mflops-a se izracunava iz srednje vrednosti vremena izvršavanja algoritma koristeći sledeću formulu:

mflops = 5Nlog_{2}(N)/(time for one FFT in microseconds)

Slika 1. Grafički prikaz podataka iz jednog od eksperimenata koje su uradili Matteo Frigo i Steven G. Johnson.

Radix 2 Cooley-Tukey algoritam

Radix 2 Cooley-Tukey algoritam se sastoji od sledećih koraka:

Ulazna sekvenca se deli na dve podsekvence. Elementi ulazne sekvence sa parnim indeksima čine jednu sekvencu, a elementi sa neparnim indeksima čine drugu sekvencu.
Izračunaj diskretnu Furijeovu transformaciju za obe sekvence
Zasnovano na diskretnoj Furijeovoj transformaciji podsekvenci, izračunaj diskretnu Furijeovu transformaciju cele ulazne sekvence

Slika 2. Grafička ilustracija algoritma radix 2. Ova ilustracija pokazuje kako funkcioniše radix 2 algoritam za ulaznu sekvencu dužine 8. Ulazna sekvenca podijeljena je na dve podsekvence. Prva podsekvenca sastoji se od elemenata sa parnim indeksima; druga podsekvenca sastoji se od elemenata sa neparnim indeksima. Ilustracija prikazuje kako se diskretna Fourierova transformacija celokupne sekvence izračunava iz diskretnih Fourierovih transformacija njenih podsekvenci.

Matematička osnova algoritma

Diskretna Furijeova transformacija cele ulazne sekvence može biti izračunata pomoću diskretne Furijeove transformacije podsekvenci parnih i neparnih elemenata na sledeći način:

{\begin{matrix}X_{k}&=&\sum \limits _{n=0}^{N-1}x_{n}e^{-{\frac {2\pi i}{N}}nk}\\&=&\sum \limits _{m=0}^{N/2-1}x_{2m}e^{-{\frac {2\pi i}{N}}(2m)k}+\sum \limits _{m=0}^{N/2-1}x_{2m+1}e^{-{\frac {2\pi i}{N}}(2m+1)k}\\&=&\sum \limits _{m=0}^{N/2-1}x_{2m}e^{-{\frac {2\pi i}{N/2}}2mk}+e^{-{\frac {2\pi i}{N}}k}\sum \limits _{m=0}^{N/2-1}x_{2m+1}e^{-{\frac {2\pi i}{N/2}}mk}\end{matrix}}

gde indeks k ide od 0 do n-1. Prva suma u poslednjem izrazu predstavlja diskretnu Furijeovu transformaciju parne podsekvence, dok druga suma predstavlja diskretnu Furijeovu transformaciju neparne podsekvence. Ako označimo prvu sumu sa $E_{k}$ , drugu sumu sa $O_{k}$ , i $e^{-2\pi /N}$ kao $W_{n}^{k}$ , dobijamo sledeću formulu:

{\begin{matrix}X_{k}&=&E_{k}+W_{N}^{k}O_{k}\end{matrix}}

Zbog osobina periodičnosti diskretne Furijeove transformacije, važe sledeće dve jednačine:

{\begin{matrix}E_{k+{\frac {N}{2}}}&=&E_{k}\\O_{k+{\frac {N}{2}}}&=&O_{k}\end{matrix}}

Zbog osobina Ojlerove formule, sledeće dve jednačine su tačne:

{\begin{matrix}e^{-{\frac {2\pi i}{N}}(k+N/2)}&=&e^{-{\frac {2\pi i}{N}}k}\\W_{N}^{k}&=&-W_{N}^{k}\end{matrix}}

Kombinacijom prethodnih jednačina sa formulama za izračunavanje diskretne Furijeove transformacije, dobijamo konačnu formulu koja nam omogućuje implementaciju algoritma koji se koristi za izračunavanje diskretne Furijeove transformacije ulazne sekvence pomoću podsekvenci parnih i neparnih članova:

{\begin{matrix}X_{k}&=&E_{k}+W_{N}^{k}O_{k}&,k=0...N/2\\X_{k+{\frac {N}{2}}}&=&E_{k}-W_{N}^{k}O_{k}&,k=0...N/2\end{matrix}}

Implementacije Cooley-Tukey algoritma

Originalna implementacija algoritma Cooley-Tukey je iterativna radix 2 implementacija. Implementacija ima optimalnu vremensku složenost $O(NlogN)$ , gde je $N$ broj podataka u ulaznoj sekvenci. Prostorna složenost prvobitne implementacije je Math > O(1) </math>. Zasnovano na formulama iz prethodnog odeljka, može se konstruisati sledeći rekurzivan algoritam za izračunavanje diskretne Furijeove transformacije:

Pseudo kod Radix 2 Cooley-Tukey Algoritma

U pseudokodu, ova procedura bi mogla biti napisana na sledeći način: ^[7]

X_0,...,N−1 ← fft_radix2(x):      
    ;specijalan slučaj kada je dužina ulazne sekvence jednaka jedan
    if N = 1 then
        X₀ ← x₀                                      
    else
     
     ;diskretna Furijeova transformacija parne podsekvence
        X_{0,...,N/2−1} ← fft_radix2(even_elements(x))             
     
     ;diskretna Furijeova transformacija neparne podsekvence
       X_N_/2,...,N−1 ← fft_radix2(odd_elements(x))           
        
      ;izračunava diskretnu Furijeovu transformaciju sekvence
      ;pomoću diskretne Furijeove transformacije neparnih i parnih podsekvenci
         
          for k = 0 to N/2−1                           
            t ← X_k
            X_k ← t + exp(−2πi k/N) X_k_+N/2
            X_k_+N/2 ← t − exp(−2πi k/N) X_k_+N/2
        endfor
    endif

Implementacija Cooley-Tukey algoritma u C++

Implementacija algoritma u programskom jeziku C++ ima optimalnu vremensku složenost $O(NlogN)$ , gde je N broj ulaznih tačaka ulazne sekvence. Prostorna složenost je $O(1)$

/* fft.cpp
 * 
 * Implementacija 
 * Cooley-Tukey rekurzivnog  FFT algoritma.
 * Kompilacija se vrši pomoću GNU/GCC kompilatora na sledeći način:
 *
 * g++ -o fft fft.cpp -lm
 *
 * Pokretanje izvršnog fajla se vrši sledećom naredbom:
 * ./fft
 *
 */

#include <complex>
#include <stdio.h>

using namespace std;

// odvaja parne/neparne elemente u donju/gornju polovinu niza redom.
// Usled "leptir" kombinacija, ispostavlja se da je ovo najjednostavniji način
// da se obavi posao bez pretvaranja pogrešnih elemenata

void separate (complex<double>* a, int n) {
    complex<double>* b = new complex<double>[n/2];  // get temp heap storage
    for(int i=0; i<n/2; i++)    // kopira sve neparne elemente na heap memoriju
        b[i] = a[i*2+1];
    for(int i=0; i<n/2; i++)    // kopira sve parne elemente na nižu polovinu niza a[
        a[i] = a[i*2];
    for(int i=0; i<n/2; i++)    // kopira sve neparne elemente (sa heap-a) na višu polovinu niza a[]
        a[i+n/2] = b[i];
    delete[] b;                 // oslobađa heap
}

// N mora biti deljivo sa dva, u suprotnom stvari mogu poći po zlu.
// U algoritmu nema provere ovog uslova.
//
// Na N ulaznih podataka u nizu X[] je primenjena brza Furijeova transformacija, 
//i rezultati ostaju u X[] 
//
// Samo prvih N/2 rezultata u X[] su rešenje brze Furijeove transformacije ulazne sekvence.

void fft2 (complex<double>* X, int N) {
    if(N < 2) {
        // Izlaz iz rekurzije.
        // Ovde se ništa ne dešava, jer je već X[0] = X[0]
    } else {
        separate(X,N);      // svi parni u nižu polovinu, svi neparni u višu polovinu niza
        fft2(X,     N/2);   // rekurzivno se obrađuju parni elementi
        fft2(X+N/2, N/2);   // rekurzivno se obrađuju neparni elementi
        // kombinuju se rezultati obe polovine niza rekurzivno
        for(int k=0; k<N/2; k++) {
            complex<double> e = X[k    ];   // parni
            complex<double> o = X[k+N/2];   // neparni
                         
            complex<double> w = exp( complex<double>(0,-2.*M_PI*k/N) );
            X[k    ] = e + w * o;
            X[k+N/2] = e - w * o;
        }
    }
}

// test program i main funkcija
int main () {
    const int nSamples = 64;
    double nSeconds = 1.0;                      // ukupno vreme izvršavanja programa
    double sampleRate = nSamples / nSeconds;    // n Hz = n / second 
    double freqResolution = sampleRate / nSamples; // frekvencija izvršavanja koraka u FFT resultatu
    complex<double> x[nSamples];                // lokacija za smeštanje ulaznih podataka
    complex<double> X[nSamples];                // lokacija za smeštanje obradjenih FFT podataka
    const int nFreqs = 5;
    double freq[nFreqs] = { 2, 5, 11, 17, 29 }; // poznate frekvencije za testiranje
    
    // generiše podatke za test primer
    for(int i=0; i<nSamples; i++) {
        x[i] = complex<double>(0.,0.);
        // postavi vrednosti sinusoida u x[]
        for(int j=0; j<nFreqs; j++)
            x[i] += sin( 2*M_PI*freq[j]*i/nSamples );
        X[i] = x[i];        // kopiraj u X[] podatke obradjene brzom Furijeovom transformacijom i rezultate
    }
    // izracunaj brzu Furijeovu transformaciju za ove podatke
    fft2(X,nSamples);
    
    printf("  n\tx[]\tX[]\tf\n");       // stampanje rezultata
    // loop to print values
    for(int i=0; i<nSamples; i++) {
        printf("% 3d\t%+.3f\t%+.3f\t%g\n",
            i, x[i].real(), abs(X[i]), i*freqResolution );
    }
}

// eof

Analiza ubrzane verzije Cooley-Tukey FFT algoritma

Slično kao i originalne CPU implementacije Cooley-Tukey algoritma, ubrzana verzija sastoji se od preuređivanja reda elemenata ulazne sekvence na okrenuti bitovski redosled.

Ubrzanje se računa sledećom formulom:

$\eta ={\frac {T_{cpu}}{T_{max}}}$

$T_{cpu}$ predstavlja vreme koje je potrebno da se izračunavanje izvrši na CPU-u, a $T_{max}$ predstavlja vreme koje je potrebno da se izračunavanje izvrši na Maxeler DFE-u.

Detalji implementacije

Kod kernela za računanje FFT prima parametar N koji predstavlja broj čvorova ulazne sekvence.^[8] Ovaj kernel ima dve ulazne i dve izlazne sekvence podataka. Te dve izlazne sekvence se koriste za transfer realnog i imaginarnog dela ulazne sekvence tačaka. Ove sekvence su tipa ArrayType koji je definisan na sledeći način:


DFEType floatType = dfeFloat(8, 24);
DFEArrayType<DFEVar> arrayType =

   new DFEArrayType<DFEVar>(floatType, n);

Koeficijenti $W_{n}^{k}$ se izračunavaju tokom kompilacije, i ne koriste se kao hardverske konstante.

Korišćenje DFE za izračunavanje FFT

MaxCompiler generiše zaglavlje FFT.h koje sadrži sledeći metod:


FFT(

  int number_of_input_sequences,

  float *imag_in,

  float *real_

  in,

  float *imag_out,

  float '''*real_out'''
  
 );

Pozivanjem ovog metoda iz C/C++ koda, šalju se podaci sa CPU-a na DFE, vrši se izračunavanje FFT-a na DFE-u i rezultat se šalje nazad CPU-u. Za MaxCompiler, da bi generisao zaglavlje FFT.h, SliC opcija mora biti uključena.

MaxJ kod Kernela koji izračunava FFT


import com.maxeler.maxcompiler.v2.kernelcompiler.Kernel;
import com.maxeler.maxcompiler.v2.kernelcompiler.KernelParameters;
import com.maxeler.maxcompiler.v2.kernelcompiler.types.base.DFEType;
import com.maxeler.maxcompiler.v2.kernelcompiler.types.base.DFEVar;
import com.maxeler.maxcompiler.v2.kernelcompiler.types.composite.DFEArray;
import com.maxeler.maxcompiler.v2.kernelcompiler.types.composite.DFEArrayType;

class FFTKernel extends Kernel {

   private static final DFEType floatType = dfeFloat(8, 24);

   protected FFTKernel(KernelParameters parameters, int n) {
     super(parameters);

     DFEArrayType < DFEVar > arrayType =

        new DFEArrayType < DFEVar > (floatType, n);
     
     DFEArray < DFEVar > in_real = io.input(“in_real”, arrayType);
     DFEArray < DFEVar > in_imag = io.input(“in_imag”, arrayType);



     int log_n = Integer.numberOfTrailingZeros(n);


     DFEVar[][] real = new DFEVar[log_n + 1][n];
     DFEVar[][] imag = new DFEVar[log_n + 1][n];


     for (int i = 0; i < n; i++) {
        int j = Integer.reverse(i) >>> (32 - log_n);

        real[0][i] = in_real[j];
        imag[0][i] = in_imag[j];
     }


     for (int i = 1, log_i = 0; i < n; i *= 2, log_i++) {
        for (int k = 0; k < i; k++) {
           float w_real = (float) Math.cos(-2 * Math.PI * k / (i * 2));
           float w_imag = (float) Math.sin(-2 * Math.PI * k / (i * 2));


        for (int j = 0; j < n; j += i * 2) {
           DFEVar temp_real =
              real[log_i][j + k + i] * w_real - imag[log_i][j + k + i] * w_imag;
           DFEVar temp_imag =
              real[log_i][j + k + i] * w_imag + imag[log_i][j + k + i] * w_real;

           real[log_i+ 1][j + k + i] = real[log_i][j + k] - temp_real;
           imag[log_i+ 1][j + k + i] = imag[log_i][j + k] - temp_imag;


           real[log_i+ 1][j + k] = real[log_i][j + k] + temp_real;
           imag[log_i+ 1][j + k] = imag[log_i][j + k] + temp_imag;
        }
     }
   }

   DFEArray < DFEVar > out_real = arrayType.newInstance(this);

   for (int i = 0; i < n; i++)
     out_real[i] <= = real[log_n][i];

   io.output(“out_real”, out_real, arrayType);

   DFEArray < DFEVar > out_imag = arrayType.newInstance(this);


   for (int i = 0; i < n; i++)
     out_imag[i] <== imag[log_n][i];

   io.output(“out_imag”, out_imag, arrayType);
   }
}

MaxJ kod Cooley-Tukey Radix-2 Algoritma


package fft;

import com.maxeler.maxcompiler.v2.build.EngineParameters;
import com.maxeler.maxcompiler.v2.managers.custom.CustomManager;
import com.maxeler.maxcompiler.v2.managers.custom.DFELink;
import com.maxeler.maxcompiler.v2.managers.custom.blocks.KernelBlock;
import com.maxeler.maxcompiler.v2.managers.custom.stdlib.Max3RingConnection;
import com.maxeler.maxcompiler.v2.managers.custom.stdlib.MaxRingBidirectionalStream;
import com.maxeler.maxcompiler.v2.managers.engine_interfaces.CPUTypes;
import com.maxeler.maxcompiler.v2.managers.engine_interfaces.EngineInterface;
import com.maxeler.maxcompiler.v2.managers.engine_interfaces.InterfaceParam;


public class FFTManager extends CustomManager {

public FFTManager(EngineParameters engineParameters, int n) {
super(engineParameters);
KernelBlock bitreversePhase = addKernel(
new BitReverseKernel(makeKernelParameters(“BitreversePhase”), n));


KernelBlock fftPhase = addKernel(
new FFTKernel(makeKernelParameters(“FFTPhase”), n, 0, Integer.numberOfTrailingZeros(n)/2, false, true));


DFELink in_real = addStreamFromCPU(“in_real”);
DFELink in_imag = addStreamFromCPU(“in_imag”);


bitreversePhase.getInput(“in_real”) <== in_real;
bitreversePhase.getInput(“in_imag”) <== in_imag;


fftPhase.getInput(“in_real”) <== bitreversePhase.getOutput(“out_real”);
fftPhase.getInput(“in_imag”) <== bitreversePhase.getOutput(“out_imag”);


MaxRingBidirectionalStream real_ring = addMaxRingBidirectionalStream(“realRingStream”,
Max3RingConnection.MAXRING_A);
MaxRingBidirectionalStream imag_ring = addMaxRingBidirectionalStream(“imagRingStream”,
Max3RingConnection.MAXRING_A);


fftPhase.getInput(“in_real_dummy”) <== real_ring.getLinkFromRemoteDFE();
fftPhase.getInput(“in_imag_dummy”) <== imag_ring.getLinkFromRemoteDFE();


real_ring.getLinkToRemoteDFE() <== fftPhase.getOutput(“out_real”);
imag_ring.getLinkToRemoteDFE() <== fftPhase.getOutput(“out_imag”);
}



static EngineInterface interfaceDefault (int n) {

EngineInterface ei = new EngineInterface();
InterfaceParam size = ei.addParam(“dataSize”, CPUTypes.UINT64);


ei.setTicks(“FFTPhase”, size);
ei.setTicks(“BitreversePhase”, size);


ei.setStream(“in_real”, CPUTypes.FLOAT, size * n * CPUTypes.FLOAT.sizeInBytes());
ei.setStream(“in_imag”, CPUTypes.FLOAT, size * n * CPUTypes.FLOAT.sizeInBytes());
return ei ;
}



public static void main(String[] args) {
EngineParameters params = new EngineParameters(args);
int n = Integer.parseInt(System.getenv(“n”));
System.out.println(n);


FFTManager manager = new FFTManager(params, n);

  //manager.setIO(IOType.ALL_CPU);
manager.addMaxFileConstant(“n”, n);
manager.createSLiCinterface(interfaceDefault(n));

  //manager.getBuildConfig().setBuildEffort(Effort.LOW);
manager.getBuildConfig().setMPPRCostTableSearchRange(1, 3);
manager.getBuildConfig().setMPPRParallelism(3);


manager.build();
}
}

MaxJ kod za okretanje bitova

package fft;
import com.maxeler.maxcompiler.v2.kernelcompiler.Kernel;
import com.maxeler.maxcompiler.v2.kernelcompiler.KernelParameters;
import com.maxeler.maxcompiler.v2.kernelcompiler.types.base.DFEType;
import com.maxeler.maxcompiler.v2.kernelcompiler.types.base.DFEVar;
import com.maxeler.maxcompiler.v2.kernelcompiler.types.composite.DFEArray;
import com.maxeler.maxcompiler.v2.kernelcompiler.types.composite.DFEArrayType;


class BitReverseKernel extends Kernel {


private static final DFEType floatType = dfeFloat(8,24);


protected BitReverseKernel(KernelParameters parameters, int n) {
super(parameters);


DFEArrayType<DFEVar> arrayType = new DFEArrayType<DFEVar>(floatType, n);


DFEArray<DFEVar> in_real = io.input(“in_real”, arrayType);
DFEArray<DFEVar> in_imag = io.input(“in_imag”, arrayType);


int log_n = Integer.numberOfTrailingZeros(n);


DFEVar[] real = new DFEVar[n];
DFEVar[] imag = new DFEVar[n];


for(int i = 0; i < n; i++) {
int j = Integer.reverse(i) >>> (32 - log_n);

  real[i] = in_real[j];
  imag[i] = in_imag[j];
}


DFEArray<DFEVar> out_real = arrayType.newInstance(this);


for (int i = 0; i < n; i++)
out_real[i] <== real[i];

io.output(“out_real”, out_real, arrayType);

DFEArray<DFEVar> out_imag = arrayType.newInstance(this);


for (int i = 0; i < n; i++)
out_imag[i] <== imag[i];
io.output(“out_imag”, out_imag, arrayType);

}
}

Zaključak

Izvođeni eksperimenti^{[који? ]} su pokazali da ubrzana implementacija ima očekivane performanse i da to daje ekvivalentne rezultate originalnoj implementaciji. Performanse ubrzane aplikacije su prevazišle performanse druge javno dostupne implementacije brze Furijeove transformacije. Zbog ograničenog broja hardverskih resursa na Maxeler kartici koja se koristi za naše eksperimente, nismo uspeli da generišemo bitstream za izračunavanje brzog Fourierovog transformiranja ulaznih sekvenci sa više od 32 tačke. Na boljim Maxeler mašinama, kao što je MAKS4 (Maia ili Coria), moguće je generisati bitstream koji izračunava brzu Fourierovu transformaciju duže ulazne sekvence.

Zbog svojstava algoritma radix 2, moguće je napraviti hibridno rješenje koje će koristiti i CPU i Maxeler DFE za izračunavanje brze Fourierove transformacije. CPU deo aplikacije bi podelio problem na manje probleme koji bi, zavisno od njihove veličine, bili rešeni na DFE-u ili ponovo podeljeni.

Očekuje se, a izvedeni eksperimenti potvrđuju, da postoji povećanje postignutog ubrzanja s povećanjem broja elemenata ulazne sekvence. Ubrzanje se povećava jer se inicijalna latencija implementacije DFE-a povećava sporije od latencije implementacije CPU-a. Takođe, prilikom izvršavanja višestrukih proračuna, nakon perioda početne latencije, implementacija DFE-a daje rezultate na svakom taktnom ciklusu dok implementacija CPU ima istu latenciju za izračunavanje svakog rezultata.

Da bi se postiglo ubrzanje korišćenjem ubrzane implementacije potrebno je izvršiti više uzastopnih proračuna brze Furijeove transformacije. Broj uzastopnih proračuna potrebnih za izvođenje, u cilju postizanja ubrzanje u poređenju sa originalnom implementacijom CPU-a, zavisi od broja tačaka u ulaznoj sekvenci.

Reference

^ Gauss, Carl Friedrich, "Theoria interpolationis methodo nova tractata^{[мртва веза]}", Werke, Band 3, 265–327 (Königliche Gesellschaft der Wissenschaften, Göttingen, 1866)
^ ^а ^б Heideman, M. T., D. H. Johnson, and C. S. Burrus, "Gauss and the history of the fast Fourier transform," IEEE ASSP Magazine, 1, (4), 14–21 (1984)
^ ^а ^б Cooley, James W.; Tukey, John W. (1965). „An algorithm for the machine calculation of complex Fourier series”. Math. Comput. 19: 297—301. doi:10.2307/2003354.
^ Cooley, James W.; Lewis, Peter A. W.; Welch, Peter D. (1967). „Historical notes on the fast Fourier transform”. IEEE Trans. on Audio and Electroacoustics. 15 (2): 76—79. doi:10.1109/tau.1967.1161903.
^ Rockmore, Daniel N. , Comput. Sci. Eng. 2 (1), 60 (2000). The FFT — an algorithm the whole family can use Special issue on "top ten algorithms of the century "„Archived copy” (PDF). Архивирано из оригинала (PDF) 2009-04-07. г. Приступљено 2009-03-31.
^ James W. Cooley, Peter A. W. Lewis, and Peter W. Welch, "Historical notes on the fast Fourier transform," Proc. IEEE, vol. 55 (no. 10), p. 1675–1677 (1967).
^ S. G. Johnson and M. Frigo, "Implementing FFTs in practice," in Fast Fourier Transforms (C. S. Burrus, ed.), ch. 11, Rice University, Houston TX: Connexions, September 2008.
^ Milutinovic, V., Salom, J., Trifunovic, N., Giorgi, R., "Guide to DataFlow SuperComputing," Springer, 2015.

Spoljašnje veze

Radix-2 Decimation in Time FFT Algorithm
Radix-2 Decimation in Frequency FFT Algorithm
a simple, pedagogical radix-2 Cooley–Tukey FFT algorithm in C++.
KISSFFT: a simple mixed-radix Cooley–Tukey implementation in C (open source)

[1] Gauss, Carl Friedrich, "Theoria interpolationis methodo nova tractata^{[мртва веза]}", Werke, Band 3, 265–327 (Königliche Gesellschaft der Wissenschaften, Göttingen, 1866)

[Heideman84-2] а ^б Heideman, M. T., D. H. Johnson, and C. S. Burrus, "Gauss and the history of the fast Fourier transform," IEEE ASSP Magazine, 1, (4), 14–21 (1984)

[CooleyTukey65-3] а ^б Cooley, James W.; Tukey, John W. (1965). „An algorithm for the machine calculation of complex Fourier series”. Math. Comput. 19: 297—301. doi:10.2307/2003354.

[4] Cooley, James W.; Lewis, Peter A. W.; Welch, Peter D. (1967). „Historical notes on the fast Fourier transform”. IEEE Trans. on Audio and Electroacoustics. 15 (2): 76—79. doi:10.1109/tau.1967.1161903.

[5] Rockmore, Daniel N. , Comput. Sci. Eng. 2 (1), 60 (2000). The FFT — an algorithm the whole family can use Special issue on "top ten algorithms of the century "„Archived copy” (PDF). Архивирано из оригинала (PDF) 2009-04-07. г. Приступљено 2009-03-31.

[6] James W. Cooley, Peter A. W. Lewis, and Peter W. Welch, "Historical notes on the fast Fourier transform," Proc. IEEE, vol. 55 (no. 10), p. 1675–1677 (1967).

[7] S. G. Johnson and M. Frigo, "Implementing FFTs in practice," in Fast Fourier Transforms (C. S. Burrus, ed.), ch. 11, Rice University, Houston TX: Connexions, September 2008.

[8] Milutinovic, V., Salom, J., Trifunovic, N., Giorgi, R., "Guide to DataFlow SuperComputing," Springer, 2015.

[1]

[2]

[3]

[4]

[5]

[6]

[7]

[8]