EÖTVÖS LORÁND TUDOMÁNYEGYETEM TERMÉSZETTUDOMÁNYI KAR
Lengyel Richárd FFT algoritmus
Szakdolgozat Matematika BSc Alkalmazott Matematikus Szakirány Témavezet®:
Schipp Ferenc Numerikus Analízis Tanszék
Budapest, 2015
Köszönetnyilvánítás Köszönettel tartozom Schipp Ferenc professzornak, aki elvállalta a konzulensi feladatkört és az utolsó pillanatig is segítette munkámat.
1
Tartalomjegyzék
1. A Fourier-transzformáció 1.1.
4
Fourier-sor . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
1.2.
Fourier-transzformáció
. . . . . . . . . . . . . . . . . . . . . . . . . . . .
2. FFT 2.1. 2.2. 2.3.
4 6
9
Jelölések, deníciók . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . Az általános algoritmus Speciális esetek
9
. . . . . . . . . . . . . . . . . . . . . . . . . . . .
13
. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
16
2.3.1.
A Cooley Tukey algoritmus
. . . . . . . . . . . . . . . . . . . . . .
16
2.3.2.
A Walsh-rendszerek . . . . . . . . . . . . . . . . . . . . . . . . . . .
20
3. Alkalmazások 3.1.
Párhuzamosítás tesztelése
3.2.
Sz¶rés
22 . . . . . . . . . . . . . . . . . . . . . . . . . . .
22
. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
24
3.3.
Polinomok szorzása . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
26
3.4.
NIST . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
28
3.4.1.
28
Az algoritmus . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
4. Forráskód
31
2
Bevezetés FFT (Fast Fourier Transform), azaz a gyors Fourier-transzformáció az egyik ismert legjelent®sebb és legfontosabb algoritmus. Ezt az is alátámasztja, hogy 2000-ben az
Guest Editors Introduction to the top 10 algorithms
IEEE:
publikációban megválasztották a 20.
század egyik legnagyobb hatással rendelkez® algoritmusának, olyan algoritmusok között mint a szimplex módszer vagy a quicksort algoritmus. Legalapvet®bb felhasználása a digitális jelfeldolgozásban van (frekvencia spketrum analízis), hiszen egy jelet az id® tartományból a frekvencia tartományba transzformál, azaz komplex sinusoidok kombinációjává alakít. Jelfeldogozásban sz¶rhetünk, konvolálhatunk vele. Ugyanakkor a legkülönfélébb helyeken is el®fordul, egészen az orvostechnikai alkalmazásoktól, a közgazdasági területeken át a kriptográáig. Polinomokat szorozhatunk segítéségével, dierenciálegyenleteket oldhatunk meg vagy akár prímtesztet is végezhetünk. Ha pedig tekintjük az algoritmus egy általánosítását, még szélesebb körben válik alkalmazhatóvá mind az alkalmazott mind az elméleti matematikában. A dolgozatban ezt az általános algoritmust fogjuk különböz® rendszereken vizsgálni. Konkrétan a trigonometrikus rendszeren és a Walsh rendszereken. Megnézzük még az algoritmus párhuzamosítását és ki is próbáljuk egy R7 m szériás grakus kártyán, az ATI-SDK implementációját használva. A Walsh rendszereken python implementációt fogunk tekinteni. A dolgozat harmadik fejezetében különböz® alkalmazásokra nézünk példákat, azaz egy képfeldolgozást, véletlenszámgenerátorok tesztelését és polinomok szorzását.
3
1. fejezet A Fourier-transzformáció
Mint a legtöbb matematikai fogalmat a Fourier-transzformációt is sokféleképp lehet bevezetni, el®ször bemutatom a Fourier-sorokat majd ennek általánosításaként a Fouriertranszformációt. Végül pedig azt még általánosabban különböz® terekre, rendszerekre is deniáljuk. Így ki fog derülni, hogy nemcsak a Fourier-sor, de a Fourier-együttható is, egy speciális esete a Fourier-transzformációnak.
1.1. Fourier-sor A Fourier-sor hasonlóan a Fourier-transzformációhoz az alkalmazott matematika egyik legfontosabb eszköze. Segítségével egy periódikus függvényt egyszer¶bb függvényekre (ortogonális sinus, cosinusok) rendszerére bonthatjuk és állíthatjuk azt vissza, így a függvény
H
vizsgálatát jelent®sen egyszser¶bbé téve. A deníciókban a
h., .i
pedig annak egy skalárszorzatát jelöli.
1.1.1. Deníció
(Ortonormált rendszer)
ortonormáltnak hívünk ha Ahol
egy tetsz®leges Hiblert-teret,
δij
.
Egy
(en ) ⊂ H,
(n ∈ N)
vektorsorozatot
hei , ej i = δij
a Kronecker-delta szimbólum, azaz:
δij = 1
ha
(i = j)
és
0
egyébként.
1.1.2. Deníció (Teljes rendszer). Egy (en ) ⊂ H, (n ∈ N) vektorsorozat teljes rendszer, ha:
∀x ∈ H : hx, en i = 0 (∀n ∈ N),
1.1.3. Deníció
akkor
x = 0.
.
(Teljes ortonormált rendszer (TONR))
Egy
(en ) ⊂ H,
torsorozat teljes ortonormált rendszer, ha teljes és ortonormált.
4
(n ∈ N)
vek-
1.1.4. Deníció (Fourier-sor, Fourier-együttható). Az x =
∞ X hx, ei iei
szummát Fourier-
i=1 sornak, az (x
∈ H , (en ) n ∈ N
hx, ei i skalárszorzatokat pedig Fourier-együtthatóknak
)A
nevezzük. Jogosan merülhet fel a kérdés hogy mikor is létezik ilyen el®állítása egy függvények. A következ® tétel erre ad választ.
1.1.1. Tétel (Fourier-sorok f®tétele). H Hilbert-tér, (en ) egy teljes ortonormált rendszer, ekkor: ∞ X
x=
hx, ei iei
(∀x ∈ H)
i=1
1.1.1. Megjegyzés (speciális esetek). Az alkalmazásokban az egyik leggyakrabban használt rendszer a komplex sinus, cosinus rendszer, azaz: Legyen
H = L2 (0, 2π),
a teljes ortonormált rendszer pedig:
en (x) =
√1 einx . Ekkor a 2π
Fourier-sorok f®tétele szerint:
Z ∞ 1 X inx 2π f (x) = e f (t)e−int dt 2π n=1 0 Ugyanis:
f (x) = Ahol
cn
∞ X
1 cn √ einx 2π n=1
(∀f ∈ H)
a Fourier-együttható, azaz
Z
2π
cn = hf, en i =
Z
0
1.1.2. Megjegyzés
2π
f en = 0
(nem periódikus a függvény)
.
1 f (t) √ e−int dt 2π El®z® esetben
2π
szerint periódikus
függvényeket néztünk, persze ez tetsz®leges periódussal is megtehet®. Ha pedig nem periódikus az intervallum, azaz
L2 (R),
akkor
∀y ∈ R : ey (x) =
√1 eiyx nem különböztethet®k 2π
meg, ellenben az el®z® esettel. Ekkor lesz szükség a Fourier-transzformációra, azaz az szerint futó szummát egy
y
szerinti integrál vált fel. Azaz:
1.1.5. Deníció (Fourier-transzformáció). Z fb(x) :=
1 f (y) √ eixy dy, 2π R
5
(x ∈ R, f ∈ L1 (R)
n
1.2. Fourier-transzformáció A következ® szakaszban még tovább terjesztjük ki a deníciót, ezzel megmutatva a kapcsolatot a Fourier-sor, Fourier-együttható és Fourier-transzformáció között mind diszkrét mind folytonos esetben. Így a transzformáció nem csak
L1 (R)-en,
de sok más téren is
használható lesz. Ehhez az alábbi fogalmak szükségesek.
1.2.1. Deníció
(normált csoport)
.
Legyen
G
egy kommutatív csoport. Normált cso-
portnak hívjuk, ha értelmezve van rajta egy norma. Azaz:
G → [0, ∞), x 7→ kxk i) kxk = 0 ⇔ x = 0 ii) k − xk = kxk (∀x ∈ G) iii) kx + yk ≤ kxk + kyk (∀x, y ∈ G)
1.2.2. Deníció (topologikus vagy folytonos csoport). Ha
∀xn , yn ∈ G (n ∈ N)
lim ρ(xn , x) = 0
sorozatra
n→∞
és
lim ρ(yn , y) = 0
n→∞
esetén teljesül:
lim ρ(xn + yn , x + y) = 0, lim ρ(−xn , −x) = 0
n→∞
1.2.3. Deníció
n→∞
(lokálisan kompakt)
.
Egy
(G, ρ)
metrikus teret lokálisan kompaktnak
nevezünk, ha:
{x ∈ G : kxk ≤ r, r > 0} Azaz a zárt gömbjei kompakt halmazok.
1.2.4. Deníció (Diszkrét csoport).
Egy
G
csoportot diszkrétnek nevezünk, ha:
inf{kxk : x 6= 0} > 0
1.2.5. Deníció
(Borel-halmaz)
dukált metrika. Ekkor a
(G, ρ)
. (G, +) egy normált Abel-csoport, ρ a norma által in-
metrikus tér nyílt részhalmazai által generált
elemei a Borel-halmazok.
6
B σ -algebra
1.2.6. Deníció (Transzláció invariáns halmaz és mérték).
A
B
Borel halmazok osztálya
transzláció invariáns, azaz
∀H ∈ B -re m
teljesül :
x + H := {x + y : y ∈ H} ∈ B
mérték transzláció invariáns, ha
∀x ∈ G, ∀H ∈ B
Borel-halmazra
m(x + H) = m(H)
1.2.1. Tétel (Haar Alfréd). Minden lokálisan kompakt topológikus csoporton van nem triviális transzláció invariáns Borel-mérték. 1.2.7. Deníció.
Az el®z® (1.2.1) tételben szerepl® mértéket Haar-mértéknek nevezzük.
A Haar-mérték konstans faktor erejéig egyértelm¶en meghatározott.
1.2.8. Deníció (M, Mm és T). M := [0, 1) intervallum a mod 1 összeadással. (Ez izomorf a R/Z faktorcsoporttal) k : k = 0, 1, ..., m − 1 az M egy m-edrend¶ ciklikus részcsoportja. Mm := m T := {z ∈ C : kzk = 1}
a szorzással az egydimenziós tórusz nev¶ csoportot alkotja
1.2.1. Megjegyzés. A komplex trigonometrikus függvény egy izomorzmust ad meg M és T között. Vagyis : M → T, (t) := ei2πt , ahol a bal oldalon
1.2.9. Deníció
M-beli
(t ∈ R) egy bijekció, amelyre: (x+y) = (x)(y)
összeadás, a jobb oldalon pedig
.
(Karakter)
A
γ :G→T
leképezés a
T-beli
(G, +)
szorzás szerepel.
csoport karaktere, ha
γ
folytonos és:
γ(x + y) = γ(x)γ(y) (∀x, y ∈ G) Azaz
γ
egy homomorzmus
G
és
T
között. Látható, hogy ez az el®z® megjegyzésben
szerepl® komplex trigonometrikus függvény tulajdonságai alapján lett deniálva. Az eddig deniált fogalmak segítségével mostmár bevezethetjük a Fourier-transzformáció fogalmát általánosabban is.
1.2.10. Deníció (Fourier-transzformált). f ∈ L1m (G) és γ karaktere a G normált, lokálisan kompakt csoportnak. Az
m
pedig a
G
csoport Haar-mértéke.
Z (Ff )(γ) := fb(γ) :=
f (t)γ(t)dm(t) G
7
Láthatjuk hogy ez valóban általánosítása az eddig bevezetett Fourier-sornak, Fourieregyütthatónak és Fourier-transzformációnak is mivel:
(Ff )
speciális esetekben:
Z T F T : fb(x) =
f (t) exp(−2πixt)dt R
Z f (t) exp(−2πint)dt
T F E : fb(n) = M
T F S : fb(x) =
X
∞ X
f (n) exp(−2πinx) =
! xn e
−2πinξ
,
xn := f (n)
n=1
n∈Z
N −1 1 X xn e−2πink/N , m n=0
1 X DF T : fb(n) = f (t) exp(−2πint) = m t∈M m
TFT: trigonometrikus Fourier-transzformáció,
k∈Z
f ∈ L1 (R), x ∈ R
f ∈ L1 (R), n ∈ Z
TFE: trigonometrikus Fourier-együttható, TFS: trigonometrikus Fourier-sor,
!
f ∈ `1 (Z), x ∈ M f : Mm → C, n ∈ Zm
DFT: diszkrét Fourier-transzformáció,
(A zárójelben a mérnöki matematikában elterjedt jelölést használva) Ahol az esetben
exp(−2πixt) m
függvény az adott tér karakterének konjugáltja. A
mérték a Lebesgue-mértéket, a
G ∈ {Z, Mm }
G ∈ {R, M}
esetben a diszkrét mértéket
jelenti.
1.2.2. Megjegyzés
.
(Inverzió)
A dolgozatban csak a DFT-nek fontos az inverze, de a
fent említett 4 transzformációt egyszerre kezelhetjük a következ® tétellel:
f, fb ∈ L1 (R) (G = R) vagy f ∈ L1 (M), fb ∈ l1 (Z) ezen feltételek mellett ha: f ∈ L1m (G) 1 b b a G csoport karaktereinek halmazát jelenti, akkor: és Ff ∈ L (G) , ahol G m
F ∗ (Ff ) = f , ahol
(F ∗ f )(x) := (Fbf )(−x),
b (x ∈ G
és
F, Fb
a
G
illetve
b G
Fourier-transzformáltja)
Példa a legfontosabb esetre:
DF T : f (x) :=
X
fb(n) exp(2πinx),
n∈Zm
8
(f ∈ L1 (Mm ))
2. fejezet FFT
Ebben a fejezetben ortogonális rendszerek egy fontos osztályával fogunk foglalkozni. Ez az osztály magába foglalja a harmonikus analízis legfontosabb rendszereit, mind folytonos és diszkrét esetben, s®t egy- és többváltozós esetben is. Ezen a rendszeren mutatunk egy általános gyors algoritmust, amely az összes ide tartozó rendszeren a Fourier-analízist és Fourier-szintézist is gyorsan el®állítja. Ahhoz hogy ezt az algoritmust megismerjük, bevezetünk néhány fogalmat, jelölést.
2.1. Jelölések, deníciók 2.1.1. Jelölés (Jn ). Jn := {A ⊂ I : |A| = 2−n }, ahol I := [0, 1) Tehát
Jn
a
[0, 1)
intervallum kett®ad részekre osztott szakaszai.
2.1.2. Jelölés (L). Ln := {f : I → C, f (Tehát
Ln
egy
2n
az
dimenziós lineáris tér, és
L :=
∞ [
Ln ,
Jn
-beli intervallumokon konstans}
L0
pedig a konstans függvények osztálya )
és
L1 := L1 (I)
n=0 Ekkor nyilvánvaló hogy:
L0 ⊂ L1 ⊂
...
Ln ⊂ L ⊂ L1
2.1.1. Deníció (En -feltételes várható érték). 1 (En f )(x) := |In (x)| Ahol
f ∈ L1 , x ∈ I
és
In ∈ Jn . 9
Z f (t)dt In (x)
2.1.1. Megjegyzés. En f (x)
In intervallumra vett 1 átlaga. Mivel a fenti denícióban In eleme az J -nek, ezért |In (x)| = . Az operátor f -et 2n Ln -be képezi. az
f
függvény
x
pontot tartalmozó
F®bb tulajdonságai:
Z
1
f (t)dt (f ∈ L)
0) E0 f = 0
i) En : L1 → Ln , lineáris operátor R R ii) (En f )(t)dt = I f (t)dt, (f ∈ L1 , I ∈ Jn ) I (λ ∈ Ln , f ∈ L1 )
iii) En (λf ) = λEn f,
(azaz az operátor
1
(f ∈ L , m ≤ n,
iv) En (Em f ) = Em (En f ) = Em f,
Ln -homogén)
m, n ∈ N)
Ezen tulajdonságok a denícióból könnyen igazolhatóak. A
0)
minden fogalom, amelyben integrál szerepel kiterjeszthet® az
tulajdonságot kihasználva
En
operátorral, hiszen
E0
esetben a szokásos integrál funkcionált kapjuk.
2.1.2. Deníció (En -ortonormált).
ϕi ∈ L2 (I)
Két
függvény
En -ortonormált,ha:
En (ϕk ϕl ) = δk,l
2.1.2. Megjegyzés.
Ha két függvény
lemben is ortonormált. Mivel a
Z
0)
és
En
iv)
értelemben ortonormált akkor, a szokásos érte-
tulajdonságokat kihasználva:
1
ϕk (t)ϕl (t)dt = E0 (En (ϕk ϕl )) = E0 (δk,l ) 0
2.1.3. Deníció (En -Fourier-együttható). ϕk ∈ L2 (I) ekkor: En (f ϕk )
2.1.3. Megjegyzés.
Látható, hogy ezek a deníciók
E0
esetben valóban a 1.1.1 és 1.1.4
deníciókat adják vissza. Algoritmikusan diadikus martingál dierenciákból
⇒ En -ortonormált
rendszerek
⇒
kö-
zönséges ortonormált rendszerek képezhet®k.
2.1.4. Deníció (UDMD). ψ = (ϕn )n∈N , ϕ ∈ L sorozat egy diadikus martingál dierencia (DMD), ha: ϕn ∈ Ln+1 Ha még a
|ϕn (x)|= 1 (x ∈ I, n ∈ N)
és
En (ϕn ) = 0 (n ∈ N)
is teljesül, akkor UDMD az angol
Matringale Dierence szavakból. 10
Unitary Dyadic
2.1.5. Deníció (Rademacher függvény). 1, 0 ≤ x < 12 r(x) := −1, 1 ≤ x < 1 2 A periódusa pedig 1, azaz : Legyen
r(x) = r(x + 1) (x ∈ R)
rn (x) := r(2n x) (x ∈ I, n ∈ N)
vagyis:
2.1.4. Megjegyzés (Rademacher-rendszer). Az (rn )n∈N rendszert Rademechar-rendszernek hívjuk. Ez a legegyszer¶bb UDMD-rendszer, ami a denícióból azonnal adódik. Az alábbi ábrán egy Rademacher-függvény látható.
Az
2.1.6. Deníció ∞ X
mn 2n ,
n=0 és
xn
mn
.
(Bináris felírás)
(m ∈ N, mn ∈ B)
bitek pedig az
2.1.5. Megjegyzés. k vel. Ha x = , 2m
rn (x)
n.
Az
függvény grafikonja
Az
∞ X xn x= , n+1 2 n=0
szummák az
x
és
m
(x ∈ I, xn ∈ B := {0, 1})
rn (x)
számok bináris/diadikus el®állításai,
függvények értékei kifejezhet®ek az
m−1
x
szám bináris jegyei-
}, m ∈ N), akkor válasszuk a végtelensok 0-ra végz®d®
el®állítást. Ekkor a Rademacher-függvények a következ® alakban állíthatóak el®:
rn (x) = (−1)xn x
m=
bináris jegyek.
(k ∈ {0, 1..., 2
Azaz a függvény csak az
és
szám bináris jegyét®l függ.
11
2.1.6. Megjegyzés.
Minden UDMD-rendszer el®áll mint:
1
I, n ∈ N, ϕ∗n ∈ Ln ) és |ϕ∗n (x)|= 1, mivel ϕn minden n−1 2 1 hosszúságú diadikus intervallumra vett átlaga 0. 2n
ϕn (x) = (−1)xn ϕ∗n (x),
(x ∈
intervallumon konstans és minden
2.1.7. Deníció (Szorzatrendszer). Vegyük a (ϕn )n∈N függvényrendszer összes lehetséges véges szorzatait. Ezt könnyen megtehetjük a bináris felírást használva.
Y
Φm :=
ϕn =
∞ Y
n ϕm k ,
(m ∈ N)
n=0
k:mn =1
Látni fogjuk, hogy igensok rendszer el®állítható mint UDMD rendszer szorzatrendszere. Ez többek között a következ® tétel miatt is hasznos.
2.1.1. Tétel. Minden UDMD-rendszer szorzatrendszere, TONR. (Ebb®l rögtön következik, hogy az UDMD rendszerek nem lehetnek teljesek, hiszen a szorzatrendszerben benne van az eredeti UDMD rendszer is.)
2.1.7. Megjegyzés. A komplex trigonometrikus rendszer m (x) := e2πimx el®áll az (2k )k∈N rendszer szorzatrendszereként. Ugyanis
m (x) := e
m-et
P k 2πi( ∞ k=0 mk 2 )x
=
bináris jegyeivel felírva:
∞ Y
2πimk 2k x
e
k=0
=
∞ Y
k m 2k (x)
k=0
A Walsh rendszer egy átrendezése pedig a Rademacher rendszer szorzatrendszereként áll el®, amelyet Walsh-Paley rendszernek szokás nevezni.
W := wn
n∈N
∞ Y rnmn , wn :=
(m ∈ N)
n=0
Az eredeti Walsh rendszer pedig az :
φ0 (x) = (−1)x0 φn (x) = (−1)xn +xn−1 , rendszer szorzatrendszere.
12
(n = 1, 2, 3...)
2.1.8. Deníció (Bitfordító transzformációk). ∞ X xN −1 xN −2 x0 xn πN (x) := + 2 + ... + N + 2 2 2 2n+1 n=N
N −1
π bN (m) := mN −1 + mN −2 2 + ... + m0 2
+ ... +
∞ X
(x ∈ I)
mn 2n
(m ∈ N)
n=N Ahol
mn
és
xN
az
x
és
2.1.8. Megjegyzés.
m
szám bináris jegyei.
Legyen
vn := (−1)xn exp(2πi( Ekkor a
(vn )n∈N
x0 xn−1 + ... + n+1 )) 2 2 2
rendszer UDMD rendszer, és az általa generált szorzatrendszer
(Vm )m∈NN
bitfordító transzformációkkal kifejezhet® a diszkrét trigonometrikus rendszerrel:
m (πN (x)) = VπbN (m) (x)
2.2. Az általános algoritmus Végig valamely UDMD rendszer által generált
Ψ
gozunk. A 2.1.1. tétel miatt ekkor a
LN -re
nézve (ahol
TONR így
NN := {0, 1, 2..., 2N − 1}),
Ψ szorzatrendszer részrendszerével dol a ΨN := ψm részrendszer teljes m∈N N
tehát a Fourier-sorok f®tétele alkalmaz-
ható és így a függvényeket jellemezhetjük a Fourier-együtthatóival. A most leírt algoritmus mind a függvény Fourier-együtthatóiból vett el®állításást (Fourier-szintézis), mind a Fourier-együtthatók kiszámítását (Fourier-analízis) a denícióval vett számolásnál sokkal gyorsabban adja meg.
2.2.1. Jelölés.
(f ∈ LN ). Azaz: 1 X k k b f (m) := hf, ψm i = N f ψm , (m ∈ NN ) N 2 k∈N 2 2N Jelöljük
fb(m)-el
az
m.
Fourier-együtthatót
N
Ekkor deníció szerint:
n f N 2 Az
LN
=
X k∈NN
beli függvényekhez összesen
összesen
2N 2N
darab szorzás és
fb k ψk
2N
n , 2N
(n ∈ NN )
darab Fourier-együtthatót kell kiszámolni. Ehhez
2N (2N −1) darab összeadás kellene. A Fourier-szintézishez 13
is ugyanennyi m¶velet lenne szükséges. Denícióval való számolás helyett használjuk ki hogy
Ψ
egy UDMD-rendszer szorzatrendszere. Tudjuk hogy
goritmus során
E0 (f ψm )-on
kívül kiszámítjuk az alábbi
fb(m) = E0 (f ψm ).
Az al-
EN (f ψm )-Fourier-együtthatókat
is:
(m ∈ [0, 2N ), 1|m)
E0 (f ψm ), E1 (f ψm ),
(m ∈ [0, 2N −1 ), 2|m)
E2 (f ψm ),
(m ∈ [0, 2N −2 ), 4|m)
E3 (f ψm ),
(m ∈ [0, 2N −3 ), 8|m) ...
En (f ψm ), Azaz kompaktabb jelöléssel:
(m ∈ [0, 2N −n ), 2n |m)
En (f ψm ),
(m = 2n l∗ , 0 ≤ l∗ < 2N −n )
Az algoritmus lényege az alábbi tulajdonságon múlik:
2.2.1. Megjegyzés. (2.1.1)
A
Ψ
UDMD (2.1.4) rendszer 1. tulajdonságát és az
iii) (Ln -homogenitás)
és
iv)
En
operátor
tulajdonságát kihasználva:
En (f ψm ) = En ϕn l0 En+1 (f ψl2n+1 )
(m = 2n l∗ , l∗ = 2l + l0 , 0 ≤ l < 2N −n−1 , l0 ∈ B)
2.2.2. Megjegyzés. En operátor deníciója miatt pedig minden f ∈ Ln+1 -re: (En+1 (f ))(x) + (En+1 (f ))(x0 ) , (En (f ))(x) = 2 Ahol
x
és
x0
a következ®k:
x :=
k , 2n
x0 := x +
k ∈ Nn , n = 0, 1..N − 1
1 2n+1
Így az általános gyors algoritmus rekurziója a következ®:
En (f ψ2n (2l+l0 ) ) (x) = ϕn l0 (x) En+1 (f ψ2n+1 l ) (x) + ϕn l0 (x0 ) En+1 (f ψ2n+1 l ) (x0 ) = 2 (l ∈ NN −n−1 , n = 0, 1..N − 1) 14
Ezt kompaktabb módon is felírhatjuk, ha bevezetjük a következ® jelöléseket:
2.2.3. Megjegyzés.
ϕn
, ahol
ϕ∗n
k 1 + n+1 n 2 2
=
−ϕ∗n
k 2n
a (2.1.6) megjegyzés szerinti függvény.
2.2.2. Jelölés.
Ekkor
k 1 k ∈ NN , n = 0, 1, ..., N − 1 αn,k := ϕn n 2 2 k (n) fk,l := En f ψ2n l k ∈ Nn , l ∈ NN −n , n = 0, 1, 2, ..., N 2n (N ) (0) fk,0 = f 2kN (k ∈ NN ) és f0,l = E0 f ψ l (l ∈ NN )
(Azaz az eredeti függvényt és a Fourier-együtthatót kapjuk vissza.) A 2.2.2 megjegyzésben szerepl® tulajdonságot és a rekurziós formulát felhasználva, azt a következ®vé egyszer¶síthetjük:
(n) fk,2l+l0
= αn,k
(n+1) f2k,l
+ (−1)
l0
(n+1) f2k+1,l
k ∈ Nn , l ∈ NN −n−1 , l0 ∈ B, n = NN −1
Ezt a rekurziót használva, a kezdeti függvényértékekb®l kiindulva, azaz
(N )
fk,0
-ekb®l elké-
(n) szítjük a fk,l szorzást és
(k ∈ NN −n , l ∈ Nn ) értékeket. Ehhez minden lépésben 2N −n 2n = 2N darab (0) (k ∈ NN ) Fourier-együtthatók kiösszeadást végzünk. Így tehát fb(j) = f 0,j
számításához összesen: ritmus
N
N2
INPUT ja az f ∈ LN
fb(j) (j ∈ NN )
darab szorzás és ugyanennyi összeadás kellett. Tehát az algofüggvényértékek. Az
Fourier-együtthatók. Ezeket az
(n) segítségével kapjuk meg, fk,l
OUTPUT
O(1).
(n = N − 1, N − 2, ..., 1, 0)
Tároljuk
A
k A(k) := f N 2
rendszer szerinti kiszámításának
sorrendben.
inplace
algoritmusként is. Ekkor
tömbben az inputot, mégpedig a következ®
képpen:
ΨN
En -Fourier-együtthatók
Az algoritmus elkészíthet® úgynevezett helybenjáró az algoritmus tárigénye
pedig a
(N )
= fk,0 15
(k ∈ NN )
Az
En+1 -Fourier-együtthatókat
pedig az alábbi módon tároljuk:
(n+1)
A(2N −n−1 k + j) := fk,j Ekkor a rekurzió szerint az a következ®
(n)
n fk,2l , fk,2l+1
(k ∈ Nn+1 , j ∈ NN −n−1 )
A(2k2N −n−1 + l), A((2k + 1)2N −n−1 + l) értékekkel kiszámoljuk
értékeket és
2k2N −n−1 + l = k2N −n + l
és
(2k + 1)2N −n−1 + l = k2N −n + 2N −n−1 + l
alapján az éppen felhasznált két adat helyen tároljuk. Ekkor az felel®en:
N . lépésben el®állnak a Fourier-együtthatók a bitfordítótranszformációnak megA(b πN (l)) = fb(l) (l ∈ NN )
2.3. Speciális esetek 2.3.1. A Cooley Tukey algoritmus A diszkrét trigonometrikus rendszeren fogunk dolgozni,most legyen
f m.
f ∈ CN . f (m)
koordinátáját jelöli. Ekkor deníció szerinti analízis és szintézis:
fb(m) := F (m) =
N −1 X
f (n)e−i2πnm/N =
n=0
=
N −1 X
f (n)(m)n1/N = hf, N (m)i,
(m = 0, 1, ..N − 1)
n=0
Ahol
N (m)
az alábbit jelenti:
−0 1/N
WN−0
−1 1/N WN−1 = 2.3.1. Deníció. WN := N := ... , így tehát: ... −(N −1) −(N −1) W 1/N N −0m 1/N 01/N (−m) WN−0m −1m WN−1m 11/N (−m) 1/N m WN = WN (m) = N (m) = = = ... ... ... (N −1) −(N −1)m −(N −1)m 1/N (−m) WN 1/N Ahol pedig
WNn -nek
a twiddle factorok:
16
az
2.3.2. Deníció (Twiddle factor). WNn := exp(2πin/N ) Látható, hogy így N-1 darab összeadást és N darab szorzást végzünk. Mivel
f
és a
karakter is komplex, így egy ilyen szorzás 4 szorzást és 2 összeadást igényel.
(a + ib)(c + id) = (ac − bd) + i(bc + ad) Így összesen
4N 2
szorzás, és
4N 2 − 2N
összeadás történt. Mivel már léteznek olyan
CPU, GPU-k amelyek 1 m¶veletként végzik el a szorzást és összeadást (a kett®t egyszerre
SIMD,FMA), így a m¶velet igény tekinthet® 4N 2 -nek.
Ezzel a jelöléssel:
F (m) :=
N −1 X
f (n)WN−nm
n=0 Az általános algoritmusból látható:
N −1 X
F (m) :=
f (n)WN−nm
+
WN−m
n:n=2k
N −1 X
K
(k = 0, 1, ..N − 1)
n:n=2k+1
És tudjuk hogy ekkor a 2 szumma, éppen egy periódikus
−(n−1)m f (n)WN
szerint, és
WKKk = 1,
N/2 pontú FFT. Viszont egy K
pontú DFT
azaz kihasználjuk hogy:
2.3.1. Jelölés. Jelölje
F01..2N −1
a
2N
pontú DFT-jét az:
(f (0), (1), ..., f (2N −1 ))
sorozatnak.
hf, W2N (m)i = F012..2N −1 (m) = F024..2N −2 (m) + W2−m N (F135..2N −1 (m)) Periódikus:
F01..2N −1 (m + 2N ) = F (m)
A komplex trigonometrikus függvény miatt:
−(m+2N )
W2N
= −W2−m N
Ezzel így megkapva a rekurziót.
2.3.1. Megjegyzés (N = 8).
Mivel
F01234567 = F0246 + W8−m (F1357 (m))
és
W8−4 = exp(−2πi4/8) = exp(−πi) = −1
F0246 [m + 4] = F0246 (m)
:
F (0) = F0246 (0) + W80 F1357 (0) = F0246 (0) + F1357 (0) F (4) = F0246 (0) + W8−4 F1357 (0) = F0246 (0) − F1367 (0) 17
és
F (1) = F0246 (1) + W8−1 F1357 (1) F (5) = F0246 (1) + W8−5 F1357 (1) = F0246 (1) − W8−1 F1357 (1) F (2) = F0246 (2) + W8−2 F1357 (2) F (6) = F0246 (2) + W8−6 F1357 (2) = F0246 (2) − W8−2 F1357 (2) F (3) = F0246 (3) + W8−3 F1357 (3) F (7) = F0246 (3) + W8−7 F1357 (3) = F0246 (3) − W8−3 F1357 (3) Majd ugyanezt a rekurzió szerint megismételjük a 4 pontú DFT-kre is:
F0246 (m) = F04 (m) + W4−m F26 (m) F1357 (m) = F15 (m) + W4−m F37 (m)
F0246 (0) = F04 (0) + F26 (0) F0246 (2) = F04 (0) − F26 (0) F0246 (1) = F04 (1) + W4−1 F26 (1) F0246 (3) = F04 (1) − W4−1 F26 (1) Hasonlóan
F1357 (m)
is, a két pontú DFT pedig:
F26 (0) = f (2)W20 + f (6)W20 = f (2) + f (6) F26 (1) = f (2)W20 + f (6)w2 −1 = f (2) − f (6)
F04 (0) = f (0) + f (4) F04 (1) = f (0) − f (4) Hasonlóan a maradék 4 darab is.
18
A teljes rekurzió lépéseit az úgynevezett pillangó diagramon szokás ábrázolni. A diagramon az
N
szimbólum a
WNi,j
faktorral való szorzást jelenti,
j
pedig a szorzás el®jelét. A
pirossal jelölt részen látható, hogy a pillangó diagram elnevezés honnan ered. A m¶velet igény is könnyen leolvasható: hiszen
log2 (N )
darab szint lesz, és minden szinten
twiddle factor m¶veletet végzünk. Így az algoritmus sebessége:
N
darab
N log2 (N )
Látható, hogy a vektorok komponensei más sorrendben érkeznek meg a m¶velet végén. Ahhoz, hogy a helyes sorrendet visszaállítsuk egy egyszer¶ bitreverziót kell végezni.
Egy 16 pontú FFT pillangó diagramja
19
2.3.2. A Walsh-rendszerek Megnézzük az FFT algoritmust a Walsh-rendszereken is. Ehhez az egyszer¶ség kedvéért egy saját python implementációt tekintünk. Az algoritmus nagyon hasonló, szinte megegyezik a diszkrét trigonometrikus rendszer esetével. A lényegi különbség az, hogy az
Ln -homogenitás
miatt, a kiemelend® faktor most elt¶nik, így a
+1, −1
marad csak
meg. Ez az implementáción jól látható:
fwht implementáció
A diszkrét trigonometrikus rendszernél, a
+1, −1-nél
még ott lenne a
WN−m
szorzó is.
Ez a függvény ebben a formában a Walsh-Paley rendszer szerint adja meg az együtthatókat. Ahhoz hogy a Hadamard szerinti sorrendben kapjuk meg az eredményt, a bitfordító transzformációt kell használnunk.
Bitfordító transzformáció
A
rev
függvény, a
bin
függvénnyel és úgynevezett
string-slice
m¶veletekkel bináris
stringgé alakítjuk a számot, feltöltjük 0-kal, és azt megfordítjuk, majd egész 10-es számrendszer beli számot készítünk bel®le. A
reverse order 20
függvény elkészíti a permutációs
vektort és ez alapján átrendezi az inputját. Egy példa a futására:
fwht futása
Az els® esetben a Wash-Paley rendszer szerinti sorrendben, másodikban a Hadamard szerinti sorrendben kaptuk meg az eredményt.
21
3. fejezet Alkalmazások
Ebben a fejezetben különböz® alkalmazásokra nézünk példákat. El®ször a párhuzamosítás tesztelését, majd az egyik legelterjedtebb alkalmazást mutatom be, a sz¶rést. Ezt különböz® képeken fogom végre hajtani mégpedig a párhuzamosított Cooley-Tookey helybenjáró algoritmust használva. A párhuzamosítás implementációját
Fixstars
általt
The OpenCL Programming Book cím¶ könyv segítségével tettem meg. A párhuzamosítás hatékonyságának tesztelésére az AMD-SDK 3.0 FFT példaprogamját 2012-ben publikált
használtam.
3.1. Párhuzamosítás tesztelése Az FFT algoritmus legtriviálisabb, párhuzamosítása a következ®képpen néz ki. Az algoritmus minden lépése során, az adott
(n)
fk,l
együtthatók kiszámítása, csak az el®z®
lépésben létrejött együtthatóktól függ, egymástól nem. Tehát minden lépés során az
(n)
fk,l
együtthatók egymástól függetlenül kiszámolhatók, így ezeket párhuzamosan végezhetjük. Ezt a pillangó-diagramon jól láthatjuk. A párhuzamosítás egy R7 260m, 2GB memóriával rendelkez® videókártyán, OpenCL 1.2 keretrendszerrel történt. A GPGPU
sing Units)
(General-Purpose computing on Praphics Proces-
programozásnak hátránya, hogy az adatok írás,olvásasa a kártya memóriája
és a host között lassú folyamat. Ez sokszor annyira lelassítja a folyamatot, hogy meg sem érné CPU helyett GPU-val számolni. A számításkor a komplex vektorban változókat deklaráltam, azaz
f loat
típusú
single precision -nel számoltam, mivel a GPU double típussal
dolgozva, több mint 20-szor lassabb. Sajnos ez a videókártya nem elég nagy teljesítmény¶ ahhoz, hogy az FFT algoritmus-
22
nál látványos eredményt érjen el. Viszont egy ennél körülbelül kétszer er®sebb GPU-val rendelkez® kártya, kétszer gyorsabb eredményt is produkál. (Gyorsabb alatt nagyobb órajellel, vagy több workitemmel, azaz nagyobb párhuzamosíthatósággal rendelkez® kártyát értünk.)
Látható, hogy a vektor méretének növelésével, lassan a GPU kezdi utolérni a CPU
4510U)
(i7
idejét, bár még 25 bites számnál is csak pár milisecundummal gyorsabb.
Teljesítmény diagram
Éppen az említett írás-olvasás lassúsága hátráltatja a GPU-t, amit nagyobb vektoroknál kezd behozni. (Hiszen a CPU ezeknél lassul, de a párhuzamosított algoritmus gyorsabb.)
23
3.2. Sz¶rés Ebben a szakaszban pár képen fogunk éldetektálást végre hajtani. A kép felfogható egy
n × n-es dimenziós mátrixként. Eddig 1 dimenziós adatszerkezeteken használtuk az t
algoritmusokat. Ahhoz, hogy egy képen ezt megtehessük, futtatni kell az t algoritmust a kép minden során, majd minden oszlopán. Ezt úgy szokás megtenni, hogy amikor az oszlopokon futtatnánk, akkor transzponáljuk a képet és azon futtatjuk. Így a számítógép memóriáját jobban kihasználjuk, s®t a párhuzamosítás miatt még inkább, mivel a transzponálás m¶veletét nagyon hatékonyan lehet párhuzamosítani. (A matlab t2 függvénye is transzponálással van megvalósítva).
A képeken a párhuzamosított implementációt hajtuk végre (aminek a sebessége a matlab t2 függvényével
256 × 256 méret¶ pgm képeken körülbelül azonos). Az t algoritmus
a képek jelét, frekvencia tartományba transzformálja, itt sz¶rünk a magas frekvenciájú jelekre (így a képen a hirtelen változásért felel®s alacsony frekvenciájú jeleket (azaz az éleket) nem engedjük át). Ennek eredménye képpen az élek fehéren fognak kiemelkedni. El®ször a standart Léna képet fogjuk vizsgálni, majd egy sajátot.
Eredeti és sz¶rt kép
24
A képre szándékosan éleket rajzoltam, így még jobban látható az algoritmus eredménye.
Eredeti és sz¶rt kép
Az alábbi képeken egy túl szigorú és egy túl enyhe sz¶rés látható.
Szigorú és enyhe sz¶rés
Az implementációban történt változások a Forráskód fejezetben megtekinthet®k.
25
3.3. Polinomok szorzása Tekintsük az alábbi két polinomot:
p(x) := a0 + a1 x + ... + aN −1 xN −1 q(x) := b0 + b1 x + ... + bN −1 xN −1 A két polinom szorzata:
p(x)q(x) = c0 + c1 x + ... + c2N −2 x2N −2 Ekkor minden
1).
ai
együtthatót meg kell szorozni
Tehát legrosszabb esetben
N2
bj
együtthatóval
0 ≤ i, j ≤ N −
darab szorzást kell végezni, (jobb eset ha vannak
együtthatók). Tehát a polinomszorzás m¶veleteigénye:
O(N log(N ))
(∀i, j
O(N 2 ).
0
Ezt az FFT algoritmussal
-re csökkenthetjük. Az implementálást ismét pythonban végeztem.
Polinomok szorzása
Az
INPUT
a
p
és
q
polinomok az együttható vektoruk szerint. A szorzást miatt a
5.
sorban kiszámoljuk mekkorára n®het a két polinom szorzata.
A
9.
4.
és
sorban megnézzük, hogy a nagyobbik hossz kett®hatvány-e. Ehhez kettes számrend-
szerbe váltjuk a számot és megnézzük, hogy tényleg csak az els® helyiértékén szerepel-e 1-es. Ha nem így van, megkeressük a legközelebbi kett®hatványt és eltároljuk a változóban.
26
power
Végül feltöltjtük nullákkal a legközelebbi kett®hatványig, így az t algoritmus gyorsan futhat le. Transzformáljuk
p-t
és
q -t
is, majd pontonként összeszorozzuk ®ket, és végül
visszatranszformáljuk. Ez így 3 db t algoritmust és N darab szorzást igényelt, tehát valóban
O(N log(N ))
lépésszámban kiszámoltuk a szorzatot.
Példa a program futására:
p(x) = x3 + 4x4 + 10x6 q(x) = 5 + 12x − 3x2 + x3 + x7 p(x)q(x) = 5x3 + 32x4 + 45x5 + 39x6 + 124x7 − 30x8 + 10x9 + x10 + 4x11 + 10x13
Program p,q input polinomokkal
A program által adott polinom pedig valóban ez, hiszen a
0-nak
tekinthet®k.
27
10−14
nagyság rend¶ számok
3.4. NIST NIST
(National Institute of Standards and Technology)
Statistical Test Suite egy
statisztikai programcsomag, amely 2001-ben azzal a céllal jött létre, hogy a véletlenszám generátorokat validáljanak. Ez 16 darab tesztet jelent, amelyet akár szoftveres, akár hardware-es bináris sorozatgenerátorra alkalmazhatunk. A tesztek a véltelenszer¶séget próbálják mérni. A tesztek között találhatunk egészen egyszer¶t és egész szosztikáltat is.
NIST Statistical Test Suite
A 6. test a DFT Teszt, másnéven a Spektrális Teszt. Itt a sorozatban szerepl® kiugró értékeket, csúcsokat vizsgáljuk. Ezzel azt próbáljuk mérni, mennyire jelennek meg olyan minták a sorozatban amelyek közel vannak egymáshoz.
3.4.1. Az algoritmus Az
INPUT
egy
a
bináris számsorozat. Ezt el®ször átalakítjük
−1, 1
érték¶vé.
xi := 2ai − 1 Alkalmazzuk DFT-t
x-en: X=x b = DF T (x)
Mivel valósból transzformálunk komplexbe, a transzformált vektornak csak az els®
n/2 da-
rab elemét vizsgáljuk a szimmetria miatt. Megszámoljuk, hogy hány komponens abszolút értéke kisebb mint
h,
legyen ez a szám
h :=
√
N1 ,
2.995732274n,
ahol:
és
28
n
az
a
input hossza.
Számoljuk ki
N0 -át,
ami az elméleti várható értéke (valódi véletlent feltételezve) azon
csúcsok számának, amik kisebbek mint h.
N0 := Számoljuk ki
d-t,
Számoljuk ki a
ami:
0.95n 2
N1 − N0 d= p n(0.95)(0.5)/2
P -értéket:
|d| P := erf c √ 2 Ahol
erf c(x) :=
Ha a kiszámolt
√2 π
R∞ x
2
e−u du
P -érték
kisebb mint
0.01,
akkor a sorozatot nem tekintjük véletlennek.
Egyébként átmegy a teszten és annak tekintjük.
NIST Spectral Test
29
Az implementáció egy stringet vár inputjának, használata az alábbi módon: Ezen az ábrán a
random.org
webservice-t használva kértünk 100 darab valódi véletlen
számot. Ami át is megy a teszten. S®t a kódot ciklusba helyezve, (pár másodperces szünetet tartva a szerver érdekében) 50-b®l 50-szer átment 1000 hosszú stringgel is. A honlap szerint majdnem mindig átmegy ezen a teszten.
Ellenben az alábbi nyilvánvalóban nem véletlen sorozattal, s®t ismeretes hogy mi magunk nehezen tudunk jó véletlen számot generálni. Ezt alá támasztja az is, hogy néhány billenty¶ leütésem nem megy át a teszten. A valóságban persze ajánlott legalább 1000 hosszú stringet, és sok vizsgálat átlagát venni, a megbízhatóbb eredmény érdekében.
30
4. fejezet Forráskód
Az itt leírt kódrészletben csak az implementációtól eltért módosításaimat írtam le. Ezek lényegében csak a kártyáról és a környezetr®l adnak némi információt, így biztosak lehetünk benne, hogy jó eszközön futtatjuk a számításokat. A maradék változtatás csak azért kellett, hogy Visual Studio 2013 compilerrel futtatható legyen a kód. Az
el.cl
és
pgm.h
leok lényegében majdnem változatlanok, de minden kód teljes méretében
megtekinthet® a
https://github.com/LengyelR/FFT-algoritmus
/∗ ... ∗/ #i f d e f __APPLE__ #include
#else #define CL_USE_DEPRECATED_OPENCL_2_0_APIS #include #endif /∗ ... ∗/
int
main ( )
{
/∗ FILE
...
∗/
∗ fp ;
errno_t
kern-
err ;
31
címen.
const char size_t
fileName [ ]
= " kernel . cl " ;
source_size ;
char ∗ s o u r c e _ s t r ; cl_int
i ,
cl_int
n;
j ;
c l _ i n t m;
size_t
gws [ 2 ] ;
size_t
lws [ 2 ] ;
// ( d e v i c e ) char ∗ v e n d o r ; char ∗ device_name ; char ∗ o p e n _ c l _ c _ v e r s i o n ; char ∗ o p e n _ c l _ v e r s i o n ; size_t
//CL_DEVICE_VENDOR //CL_DEVICE_NAME //CL_DEVICE_OPENCL_C_VERSION // CL_DEVICE_VERSION // CL_DEVICE_MAX_WORK_GROUP_SIZE
∗ max_workgroup ;
// ( p l a t f o r m ) char ∗ p r o f i l e = NULL ; char ∗ p l a t f o r m _ v e r s i o n
= NULL ;
/ ∗ Kernel o l v a s a s a ∗ / e r r = f o p e n _ s (& f p ,
if
( err )
fileName ,
"r" );
{ p r i n t f ( "HIBA :
Failed
to
load
k e r n e l . \ n" ) ;
exit (1); } source_str = (
char ∗ ) m a l l o c (MAX_SOURCE_SIZE ) ;
source_size = fread ( source_str ,
1 , MAX_SOURCE_SIZE,
f c l o s e ( fp ) ;
/∗
...
∗/
/ ∗ Device es Platform b e a l l i t a s a ∗ / cl_int
status ;
32
fp ) ;
cl_uint
numPlatforms ;
cl_uint
numDevices ;
size_t
size ;
∗ platforms ;
cl_platform_id cl_device_id
∗ devices ;
clGetPlatformIDs (0 , p r i n t f ( "Number
of
NULL,
&n u m P l a t f o r m s ) ;
P l a t f o r m s : %d \ n \ n " ,
platforms = ( cl_platform_id
∗)
m a l l o c ( numPlatforms ∗
sizeof ( c l _ p l a t f o r m _ i d ) ) ;
c l G e t P l a t f o r m I D s ( numPlatforms ,
clGetDeviceIDs ( platforms [ 1 ] , NULL,
platforms ,
NULL ) ;
CL_DEVICE_TYPE_ALL,
0,
&n u m D e v i c e s ) ;
devices = ( cl_device_id
∗)
m a l l o c ( numDevices ∗
sizeof ( c l _ d e v i c e _ i d ) ) ;
clGetDeviceIDs ( platforms [ 1 ] , numDevices ,
devices ,
CL_DEVICE_TYPE_ALL, NULL ) ;
c o n t e x t = c l C r e a t e C o n t e x t (NULL, NULL,
numPlatforms ) ;
NULL,
numDevices ,
devices ,
&s t a t u s ) ;
cmdQueue = clCreateCommandQueue ( c o n t e x t , devices [1] ,
0 , &s t a t u s ) ;
/ ∗ Platform i n f o ∗ / clGetPlatformInfo ( platforms [ 0 ] , NULL,
profile ,
CL_PLATFORM_PROFILE,
&s i z e ) ;
33
profile
= (
char ∗ ) m a l l o c ( s i z e ) ;
clGetPlatformInfo ( platforms [ 0 ] , size ,
profile ,
CL_PLATFORM_PROFILE,
NULL ) ;
clGetPlatformInfo ( platforms [ 0 ] , NULL,
CL_PLATFORM_VERSION,
platform_version ,
platform_version = (
&s i z e ) ;
char ∗ ) m a l l o c ( s i z e ) ;
clGetPlatformInfo ( platforms [ 0 ] , size ,
CL_PLATFORM_VERSION,
platform_version ,
p r i n t f ( " Platform p r i n t f ( " ( Status
I n f o r m a t i o n s : \ n" ) ; now : %d ) \ n " ,
p r i n t f ( " P r o f l e : %s \ n " , p r i n t f ( " Platform
NULL ) ;
status );
profile );
V e r s i o n : %s \ n " ,
platform_version ) ;
/ ∗ Device i n f o ∗ / clGetDeviceInfo ( devices [ 0 ] , NULL, vendor = (
NULL,
CL_DEVICE_VENDOR,
&s i z e ) ;
char ∗ ) m a l l o c ( sizeof ( char ) ∗ s i z e ) ;
clGetDeviceInfo ( devices [ 0 ] , size ,
vendor ,
NULL ) ;
clGetDeviceInfo ( devices [ 0 ] , NULL,
NULL,
device_name = (
CL_DEVICE_NAME,
&s i z e ) ;
char ∗ ) m a l l o c ( sizeof ( char ) ∗ s i z e ) ;
clGetDeviceInfo ( devices [ 0 ] , device_name ,
NULL,
CL_DEVICE_NAME,
size ,
NULL ) ;
clGetDeviceInfo ( devices [ 0 ] , NULL,
CL_DEVICE_VENDOR,
CL_DEVICE_OPENCL_C_VERSION,
&s i z e ) ;
open_cl_c_version = (
char ∗ ) m a l l o c ( sizeof ( char ) ∗ s i z e ) ;
clGetDeviceInfo ( devices [ 0 ] ,
CL_DEVICE_OPENCL_C_VERSION,
34
size ,
open_cl_c_version ,
clGetDeviceInfo ( devices [ 0 ] , NULL,
NULL ) ;
CL_DEVICE_VERSION,
NULL,
&s i z e ) ;
open_cl_version = (
char ∗ ) m a l l o c ( sizeof ( char ) ∗ s i z e ) ;
clGetDeviceInfo ( devices [ 0 ] , open_cl_version ,
CL_DEVICE_VERSION,
size ,
NULL ) ;
c l G e t D e v i c e I n f o ( d e v i c e s [ 0 ] , CL_DEVICE_MAX_WORK_GROUP_SIZE, NULL,
NULL,
&s i z e ) ;
max_workgroup = ( s i z e _ t
∗ ) m a l l o c ( sizeof ( s i z e _ t ) ∗ s i z e
);
c l G e t D e v i c e I n f o ( d e v i c e s [ 0 ] , CL_DEVICE_MAX_WORK_GROUP_SIZE, size ,
&max_workgroup ,
clGetDeviceInfo ( devices [ 0 ] , NULL,
NULL,
max_cu = ( c l _ u i n t
CL_DEVICE_MAX_COMPUTE_UNITS,
&s i z e ) ;
∗ ) m a l l o c ( sizeof ( c l _ u i n t ) ∗ s i z e
clGetDeviceInfo ( devices [ 0 ] , size ,
NULL ) ;
&max_cu ,
p r i n t f ( " \n\ n D ev i c e
CL_DEVICE_MAX_COMPUTE_UNITS,
NULL ) ;
I n f o r m a t i o n s : \ n" ) ;
p r i n t f ( " Vendor : %s \ n " ,
vendor ) ;
p r i n t f ( " D e v i c e : %s \ n " ,
device_name ) ;
p r i n t f ( " openCL C : %s \ n " , p r i n t f ( " openCL : %s \ n " ,
open_cl_c_version ) ;
open_cl_version ) ;
p r i n t f ( "Max WorkGroup : %d \ n " , p r i n t f ( "Max ComputeUnits :
/∗ return
...
);
max_workgroup ) ;
: %d \ n \ n " ,
∗/
0;
}
35
max_cu ) ;
Irodalomjegyzék
[1] Schipp Ferenc,
Fourier Analízis
[2] Schipp Ferenc,
Párhuzamos FFT Algoritmusok , (vázlat)
[3] Ryoji Tsuchiyama, Takashi Nakamura, Takuro Iizuka, Akihiro Asahara, Jeongdo Son, Satoshi Miki,
The OpenCL Programming Book
[4] Eric W. Hansen,
Fourier Transforms, Principles and Applications
[5] Andrew Rukhin, Juan Soto, James Nechvatal, Miles Smid, Elaine Barker, Stefan Leigh, Mark Levenson, Mark Vangel, David Banks, Alan Heckert, James Dray, San Vo
A Statistical Test Suite for Random Number Generators for Cryptographic Applications [6] Yan-Bin Jia, [7]
Polynomial Multiplication and Fast Fourier Transform
http://developer.amd.com/tools-and-sdks/opencl-zone/ amd-accelerated-parallel-processing-app-sdk
36