Vysoké učení technické v Brně Fakulta strojního inženýrství Ústav fyzikálního inženýrství
DIPLOMOVÁ PRÁCE
Netradiční metody výpočtu difrakčních jevů v optice
Diplomant: Jan Zámečník Vedoucí diplomové práce: Prof. RNDr. Jiří Komrska, CSc. Brno, květen 2004
Prohlášení: Prohlašuji, že jsem diplomovou práci vypracoval samostatně pod vedením prof. RNDr. Jiřího Komrsky, CSc. a že veškerá použitá literatura je v této práci uvedena.
1
Poděkování: Děkuji prof. RNDr. Jiřímu Komrskovi, CSc. za časté konzultace a cenné připomínky během tvorby této diplomové práce. Dále děkuji Ing. Karlu Velechovskému za pomoc při fotografování difrakčních jevů a firmě eMZet, v.o.s. za financování tiskařských prací. V neposlední řadě děkuji rodičům za podporu po celou dobu mého studia.
2
Obsah 1 Úvod
5
2 Difrakce, difrakční zóny a difrakční integrály
6
2.1
Difrakce - experimentální uspořádání . . . . . . . . . . . . . . . . . . . . .
6
2.2
Difrakční zóny . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
7
2.3
Difrakční integrály . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
8
3 Difrakční integrál pro Fraunhoferovu difrakci 3.1
Funkce propustnosti . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
3.2
Fraunhoferova difrakce jako Fourierova transformace vlnové funkce v rovině stínítka . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
8 8
9
4 Difrakční integrál pro Fresnelovy difrakční jevy
10
5 Spojitá Fourierova transformace
11
6 Diskrétní Fourierova transformace (DFT)
12
6.1
Vzorkovací teorém . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 13
6.2
Definice diskrétní Fourierovy transformace . . . . . . . . . . . . . . . . . . 15
6.3
Inverzní diskrétní Fourierova transformace (IDFT) . . . . . . . . . . . . . . 16
6.4
Modifikace diskrétní Fourierovy transformace . . . . . . . . . . . . . . . . . 17
6.5
Vlastnosti diskrétní Fourierovy transformace . . . . . . . . . . . . . . . . . 18
6.6
6.5.1
Linearita diskrétní Fourierovy transformace . . . . . . . . . . . . . 19
6.5.2
Diskrétní Fourierova transformace a posunutí funkce . . . . . . . . 19
6.5.3
Diskrétní Fourierova transformace a posunutí Fourierova obrazu . . 20
6.5.4
Diskrétní Fourierova transformace rozšířené funkce o nulové hodnoty 20
Aplikace diskrétní Fourierovy transformace . . . . . . . . . . . . . . . . . . 21 6.6.1
Interpolace funkčních hodnot pomocí diskrétní Fourierovy transformace . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 21
6.6.2
Diskrétní konvoluce periodických funkcí . . . . . . . . . . . . . . . . 22
6.6.3
Diskrétní konvoluce „neperiodickýchÿ funkcí definovaných na intervalu konečné délky . . . . . . . . . . . . . . . . . . . . . . . . . . . 23
3
6.6.4
Fourierova transformace diskrétní konvoluce a diskrétní Fourierova transformace součinu . . . . . . . . . . . . . . . . . . . . . . . . . . 24
6.7
Příklady výpočtu spojité Fourierovy transformace pomocí diskrétní Fourierovy transformace . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 24 6.7.1 6.7.2 6.7.3
Příklad: h(lT ) = cos(alT ), . . . . . . . . . . . . . . . . . . . . . . . 25 lT Příklad: h(lT ) = rect 2a . . . . . . . . . . . . . . . . . . . . . . . . 26 lT Příklad: rect 2a ∗ rect lT . . . . . . . . . . . . . . . . . . . . . . . 27 2b
7 Rychlá Fourierova transformace (FFT)
28
7.1
Maticové vyjádření diskrétní Fourierovy transformace . . . . . . . . . . . . 28
7.2
Rychlá Fourierova transformace . . . . . . . . . . . . . . . . . . . . . . . . 30
8 Vícerozměrná rychlá Fourierova transformace 8.1
32
Dvourozměrná diskrétní Fourierova transformace . . . . . . . . . . . . . . . 32 8.1.1
Přeskládání kvadrantů diskrétní Fourierovy transformace pro konvenční zobrazení . . . . . . . . . . . . . . . . . . . . . . . . . . . . 33
8.2
Vícerozměrná diskrétní Fourierova transformace . . . . . . . . . . . . . . . 34
9 Výpočet Fraunhoferových a Fresnelových difrakčních jevů pomocí rychlé Fourierovy transformace
34
9.1
Postup výpočtu rychlé Fourierovy transformace . . . . . . . . . . . . . . . 34
9.2
Výpočet Fraunhoferovy difrakce pomocí rychlé Fourierovy transformace . . 35
9.3
Výpočet Fresnelovy difrakce pomocí rychlé Fourierovy transformace . . . . 37
10 Obrazová část
39
10.1 Výpočty . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 39 10.2 Experimenty . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 48 10.2.1 Fraunhoferova difrakce . . . . . . . . . . . . . . . . . . . . . . . . . 48 10.2.2 Fresnelova difrakce . . . . . . . . . . . . . . . . . . . . . . . . . . . 54 11 Závěr
59
4
1
Úvod
V polovině 60. let minulého století vytvořili J.W. Cooley a J.W. Turkey [1] algoritmus, známý pod názvem rychlá Fourierova transformace. Vznik tohoto algoritmu znamenal obrovský posun v rychlosti výpočtu diskrétní Fourierovy transformace. Dokonce rychlá Fourierova transformace se stala pojmem tak významným, že někteří lidé, i z akademických kruhů, o ní hovoří s velikou úctou. Rychlost výpočtu diskrétní Fourierovy transformace pomocí algoritmu rychlé Fourierovy transformace je opravdu hodna obdivu, bohužel mnozí přisuzují rychlé Fourierově transformaci spoustu jí nenáležejících vlastností, což plyne pouze z nepochopení diskrétní Fourierovy transformace. Prvotní chyba tkví v samotném odvození diskrétní Fourierovy transformace (viz např. [2], kap. 6) N −1 n X H h(lT ) exp(−i2πln/N ) , = NT l=0
(1)
kde l = 0, 1, ..., N − 1 a n = 0, 1, ..., N − 1, N je počet vzorků. V odvození diskrétní Fourierovy transformace se vůbec neuvažuje symetrie definičního oboru funkce h(lT ) a navíc absence (resp. odstranění) konstant před sumou (1) nutí autory k přidání speciální kapitoly, zabývající se vztahem mezi diskrétní a spojitou Fourierovou transformací a těmito konstantami. Důsledky těchto chyb (např. uspořádání prvků spočtené Fourierovy transforamce (viz [3], obr. 12.2.2)) jsou potom nesprávně přisuzovány rychlé Fourierově transformaci. V textu diplomové práce se přístupuje k definici diskrétní Fourierovy transformace tak, že je definována přímo jako numerický výpočet spojité Fourierovy transformace 5(1). Hlavní část diplomové práce se zabývá využitím rychlé Fourierovy transformace při výpočtech difrakčních jevů v optice a stanovením měřítka daného difrakčního jevu. V kapitole 10 jsou potom příklady výpočtu některých difrakčních jevů a také ověření shodnosti teorie a experimentu v podkapitole 10.2.
5
2
Difrakce, difrakční zóny a difrakční integrály
Na úvod této kapitoly je vhodné odcitovat zakladatele teorie difrakce F. M. Grimaldiho, který také vytvořil slovo difrakce: ”Světlo se šíří nebo proniká nejen přímo, lomem nebo odrazem, nýbrž někdy ještě také jakýmsi čtvrtým způsobem, difrakcí.” Toto tvrzení dobře vymezuje pojem difrakce. Samotné slovo difrakce je složeno ze dvou latinských slov: disvyjadřující opak, negaci a frangere- s hlavním významem lámat.
2.1
Difrakce - experimentální uspořádání
Mějme optickou soustavu skládající se ze zdroje světla (např. bodového), který je na optické ose (viz obr. 1). Ve vzdálenosti z1 od tohoto zdroje umístíme stínítko. Světlo se nemůže šířit bez omezení, neboť mu v cestě stojí stínítko. Světlo se začne šířit již zmíněným čtvrtým způsobem difrakcí. V rovině pozorování, vzdálené od stínítka o vzdálenost z, není obraz stínítka ostrý. Pozorujeme zde Fresnelův difrakční obrazec. Oblast za stínítkem, kde je obraz ještě ostrý, je oblastí geometrické optiky. Pak následuje oblast Fresnelovy difrakce. Pokud bychom šli dostatečně daleko (resp. do nekonečna) pozorovali bychom Fraunhoferovu difrakci.
Obrázek 1: Experimentální uspořádání pro difrakci světla. Fresnelovy difrakční jevy se objevují při stínové projekci, zvláště pak, použijeme-li malého monochromatického zdroje (resp. koherentního zdroje osvětlení). Fresnelovy difrakční jevy se objevují i v optických soustavách v místech, kde je zobrazení neostré. Rozostřený obraz lze považovat, při koherentním záření, za Fresnelovu difrakci.
6
Fraunhoferova difrakce je pro vysvětlení trochu složitější. Fraunhoferova difrakce by se dala považovat za speciální případ Fresnelovy difrakce, kdy s rovinou pozorování jdeme do nekonečna. Protože bodům roviny pozorování v nekonečnu odpovídají směry, lze chápat tuto situaci tak, že na stínítko dopadá rovnoběžný svazek světla, tj. rovinná vlna, a Fraunhoferova difrakce tedy vyjadřuje směrové rozložení intenzity difraktovaného světla.
2.2
Difrakční zóny
Pro plánování difrakčních experimentů jsou důležité difrakční zóny. Mějme osvětlující vlnu, která dopadá na stínítko S (viz obr. 2). Potom n propuštěných difrakčních zón odpovídá dráhovému rozdílu nλ/2 mezi paprskem jdoucím od středu stínítka M0 a paprskem jdoucím od okraje stínítka (resp. od nejvzdálenějšího bodu stínítka od středu M0 ) v rovině pozorování. Nezáleží přitom na tvaru vlnoplochy osvětlující vlny (např. jestli jde o rovinnou vlnu, kulovou vlnu atd.) Pro Fresnelovu difrakci se běžně používají Fresnelovy zóny (viz např. [4], kap. 5.4), které jsou však v našem případě nevýhodné, neboť pro každý tvar vlnoplochy existuje jiný tvar vztahu pro výpočet počtu propuštěných Fresnelových zón. Z tohoto důvodu, trochu z nouze, jsme zvolili název difrakční zóny. Počet difrakčních zón odpovídá počtu Fresnelových zón pro rovinnou vlnu. Pro poloměr rn n-té difrakční zóny platí √ rn =
r nλz
1+
nλ , 4z
(2)
kde z je vzdálenost roviny pozorování od roviny stínítka a λ je vlnová délka osvětlující vlny.
Obrázek 2: Poloměr rn n-té difrakční zóny. 7
Pro běžné experimentální uspořádání lze vztah (2) aproximovat vztahem √ rn =
nλz.
(3)
rn2 . λz
(4)
Pro počet n difrakčních zón platí n=
Běžným experimentálním uspořádáním je míněno, že vzdálenost z roviny pozorování od stínítka je řádu metrů, vlnová délka λ je řádu stovek nanometrů a poloměr rn je řádu milimetrů. Pomocí počtu n difrakčních zón lze kvantitativně rozdělit oblasti za stínítkem. Jestliže n > 100, jde o oblast geometrické optiky, tj. oblast přímo za stínítkem. Oblast, kde platí n << 1 (resp. n → 0, tj. s rovinou pozorování jdeme do nekonečna), je oblastí Fraunhoferovy difrakce. Oblast mezi je oblastí Fresnelovy difrakce. Nelze přitom přesně stanovit hranice mezi oblastmi.
2.3
Difrakční integrály
Odvozovat difrakční integrály pro Fraunhoferovu a Fresnelovu difrakci by bylo zbytečné, neboť se tato diplomová práce nezabývá odvozením difrakčních integrálů, ale jejich výpočty. Fraunhoferova difrakce může být popsána jako Fourierova transformace funkce propustnosti stínítka (viz [5], kap. 4). Toto bude ukázáno v kapitole 3. Odvození Fresnelova difrakčního integrálu je popsáno např. v [4], kap. 4, 5, 7 a 8. Difrakčním integrálem pro Fresnelovy difrakční jevy se zabývá kapitola 4.
3
Difrakční integrál pro Fraunhoferovu difrakci
Tato kapitola je převzata a upravena z [5], kap 4.
3.1
Funkce propustnosti
Funkci propustnosti t(x, y) lze zavést takto: Předpokládejme, že máme dvojrozměrný objekt (v praxi např. fotografický film), který je alespoň v některých místech transparentní pro nějaké vlnění. Nechť toto vlnění (označme jej ψ0 (x, y)) dopadá na objekt, prochází jím a těsně za ním (tj. v zadní rovině objektu) je charakterizováno funkcí ψ(x, y). 8
Podíl t(x, y) =
ψ(x, y) ψ0 (x, y)
(1)
nazveme funkcí propustnosti, pokud je v rozumných mezích nezávislý na dopadající vlně ψ0 (x, y) (rozumné meze vylučují případy jako téměř tečný dopad rovinné vlny, kulovou vlnu se středem velmi blízko objektu apod.). Funkce propustnosti je obecně komplexní funkcí. V praxi bývají důležité dva případy: (i) Amplitudové objekty. Mají funkci propustnosti ve tvaru t(x, y) = τ (x, y) exp(iε0 ),
(2)
kde τ (x, y) je reálná funkce a ε0 je reálná konstanta. Touto funkcí charakterizujeme např. objekty tvořené otvory v nepropustném stínítku. Potom funkce propustnosti je rovna jedné tam, kde jsou otvory a nule tam, kde otvory nejsou. (ii) Fázové objekty. Mají v propustných částech funkci propustnosti ve tvaru t(x, y) = τ0 exp[iε(x, y)],
(3)
kde τ0 je reálná konstanta a ε(x, y) je reálná funkce polohy. V nepropustných částech je t(x, y) = 0. Tato funkce propustnosti charakterizuje optické elementy, které absorbují světlo jen nepatrně, avšak podstatně ovlivňují fázi vlnění (např.: tenká čočka, fázová mřížka apod.)
3.2
Fraunhoferova difrakce jako Fourierova transformace vlnové funkce v rovině stínítka
Amplitudu a(X, Y ) difraktovaného vlnění ve směru (X, Y,
√
1 − X 2 − Y 2 ) vyjadřuje difrakční
integrál a(X, Y ) =
k 2 2π
Z∞Z ψ(x, y) exp[−ik(Xx + Y y)] dxdy,
(4)
−∞
kde ψ(x, y) je vlnová funkce v rovině stínítka, X a Y jsou směrové kosiny směru, kterým se difraktované světlo šíří, vázané podmínkou X 2 + Y 2 ≤ 1. |a(X, Y )|2 určuje směrové rozložení intenzity za stínítkem. 9
Vlnovnou funkci ψ(x, y) lze vyjádřit pomocí funkce propustnosti takto (viz 3.1(1)) ψ(x, y) = t(x, y) · ψ0 (x, y),
(5)
kde ψ0 (x, y) je vlnová funkce osvětlující vlny v rovině z = 0, tj. v rovině stínítka a t(x, y) je funkce propustnosti stínítka (viz kap. 3.1). Dosazením (5) do (4) dostaneme a(X, Y ) =
k 2 2π
Z∞Z t(x, y)ψ0 (x, y) exp[−ik(Xx + Y y)] dxdy.
(6)
−∞
Využitím tvrzení, že Fourierova transformace součinu dvou funkcí je úměrná konvoluci Fourierových transformací těchto funkcí (viz [5], kap. 9.3), spočteme Fourierovu transformaci jako konvoluci a(X, Y ) = FT{t(x, y)} ∗ FT{ψ0 (x, y)},
(7)
kde FT{f (x, y)} značí Fourierovu transformaci funkce f (x, y). Jestliže osvětlujeme stínítko rovinou vlnou, tj. ψ0 (x, y) = b exp(−ikz), kde b je amplituda vlny, je Fraunhoferova difrakce úměrná Fourierově transformaci funkce propustnosti. To, že Fraunhoferovu difrakci spočteme jako konvoluci Fourierovy transformace funkce propustnosti s Fourierovou transformací vlnové funkce osvětlující vlny, můžeme chápat i tak, že Fraunhoferovu difrakci modulujeme osvětlující vlnou ψ0 (x, y).
4
Difrakční integrál pro Fresnelovy difrakční jevy
Vlnovou funkci v rovině pozorování ve vzdálenosti z (z = konst. > 0) od roviny stínítka vyjadřuje difrakční integrál pro Fresnelovy difrakční jevy (viz [4], 4.4(14)) k exp(ikz) ψ(x, y, z) = −i 2π z
Z∞Z
ik 2 2 ψM (xM , yM ) exp (x − xM ) + (y − yM ) dxM dyM , 2z
−∞
(1) kde ψM (xM , yM ) je vlnová funkce v rovině stínítka. |ψ(x, y, z)|2 určuje rozložení intenzity v difrakčním obrazci. Vlnová funkce ψM (xM , yM ) v rovině stínítka v sobě zahrnuje vlnovou funkci osvětlující vlny a funkci propustnosti stínítka. Vlnovou funkci v rovině stínítka lze nahradit součinem (viz 3.1(1)) 10
ψM (xM , yM ) = t(x, y) · ψ0 (xM , yM ),
(2)
kde t(x, y) je funkce propustnosti stínítka a ψ0 (xM , yM ) je vlnová funkce osvětlující vlny. Rozvinutím druhých mocnin v argumentu fázoru a vytknutím fázoru nezávislého na integračních proměnných před integrál dostaneme k exp(ikz) ik 2 2 ψ(x, y, z) = −i exp (x + y ) × 2π z 2z Z∞Z
ik ik 2 2 (x + yM ) exp − (xxM + yyM ) dxM dyM . (3) ψM (xM , yM ) exp 2z M z
× −∞
Vlnová funkce charakterizující Fresnelovu difrakci je tak vyjádřena Fourierovou trans ik 2 2 formací součinu vlnové funkce ψM (xM , yM ) v rovině difrakčního stínítka a fázoru exp 2z (xM + yM ) . K numerickému výpočtu Fresnelovy difrakce lze tedy použít metod pro výpočet Fourierovy transformace. Postup výpočtu Fraunhoferových a Fresnelových difrakčních jevů bude ukázán v kapitole 10.1.
5
Spojitá Fourierova transformace
Tato diplomová práce se zabývá netradičními metodami výpočtu difrakčních jevů. Analytické výpočet difrakčních integrálů je tradiční metoda výpočtu difrakčních jevů. Fourierova transformace je základem pro výpočty difrakčních jevů, a proto je velmi důležité si ji na tomto místě definovat, neboť z ní vychází definice diskrétní Fourierovy transformace, která je užitečná pro numerické výpočty difrakčních jevů a ve výpočetní technice vůbec. Numerickou metodu, jakožto netradiční metodu výpočtu, použijeme pro výpočty difrakčních jevů, které nelze provést analyticky. Obecně, použijeme ji tehdy, chceme-li počítat difrakční jevy pomocí výpočetní techniky. Definujme tedy Fourierovu transformaci FT{f (~x)} a inverzní Fourierovu transformaci ~ integrály (viz [5], 2.1(1) a 2.1(2)) FT−1 {F (X)}
FT{f (~x)} = AN
∞Z
Z
···
~ · ~x) dN ~x, f (~x) exp(−ik X
(1)
−∞
~ FT−1 {F (X)} = BN
Z
∞Z
··· −∞
11
~ exp(ik X ~ · ~x) dN X. ~ F (X)
(2)
~ jsou absolutně integrovatelné po částech hladké Předpokládá se přitom, že f (~x) a F (X) ~ ∈ EN a že tyto funkce mají nespojitosti komplexní funkce reálných proměnných ~x, X pouze v bodech, které tvoří množinu míry nula v EN . Konstanty A, B, k jsou svázány podmínkou AB =
|k| , 2π
(3)
jež vyplývá z tzv. fundamentální věty o Fourierově transformaci, která říká, že funkce f (~x), pro kterou existuje Fourierův integrál, je v bodech spojitosti f (~x) = FT−1 {FT{f (~x)} = FT FT−1 {f (~x) .
(4)
V bodech nespojitosti funkce f (~x) dává aplikace Fourierovy transformace a inverzní Fourierovy transformace střední hodnotu funkce f (~x) v infinitesimálním okolí bodu nespojitosti (viz [5], 2.2). V různých aplikacích se konstanty A, B, a k definují různě. Konstanty se definují tak, aby ~ rozměrově i jednotkově správný. Např.: jestliže budeme byl vztah mezi souřadnicemi ~x a X uvažovat jednorozměrnou funkci f (x) a souřadnice x bude vyjadřovat čas (jednotka ssekunda), potom při volbě konstant, k = 2π bude souřadnice X vyjadřovat frekvenci (jednotka s−1 = Hz - Hertz). Jelikož se budeme zabývat aplikací Fourierovy transformace v optice, souřadnice x vyjadřuje vzdálenost, souřadnice X je směrový kosinus (kosinus úhlu měřeného od kolmice na směr šíření paprsku) a konstanty zvolíme následovně: A = 1/λ, B = 1, k = 2π/λ. V textu budeme i nadále používat obecné označení A, B, a k pro ~ · ~x a ik X ~ · ~x ve Fourierově transformaci musí snadnější orientaci. Argumenty fázoru −ik X být bezrozměrné.
6
Diskrétní Fourierova transformace (DFT)
V této kapitole se budeme zabývat diskrétní Fourierovou transformací (anglicky - Discrete Fourier Transform, zkr. DFT). Diskrétní Fourierovu transformaci odvodíme poněkud odlišně, než je tomu běžně v literatuře, neboť v literatuře se v závěru odvození musí stanovit vztah mezi diskrétní a spojitou Fourierovou transformací (viz např. [2], kap. 6.4). Definování vztahu mezi diskrétní a spojitou Fourierovou transformací je nelogické, vzhledem k tomu, že definice Fourierovy transformace je jen jedna a diskrétní Fourierova transformace má být pouze nástrojem pro výpočet spojité Fourierovy transformace. Diskrétní 12
Fourierova transformace má sloužit pro výpočet Fourierovy transformace funkcí, který nelze provést analyticky a dále pro výpočty Fourierovy transformace pomocí výpočetní techniky. Omezení tkví v tom, že nelze jednoduše spočítat Fourierovu transformaci pomocí diskrétní Fourierovy transformace v mezích (−∞; ∞), které požaduje definice Fourierovy transformace 5(1) a zároveň nelze spočíst diskrétní Fourierovu transformaci spojitých funkcí. Toto omezení vnáší rozdíl mezi diskrétní a spojitou Fourierovou transformaci.
6.1
Vzorkovací teorém
Odvození diskrétní Fourierovy transformace provedeme pro jednorozměrnou funkci f (x), pro vícerozměrnou funkci je odvození analogické. Odvození diskrétní Fourierovy transformace vychází ze vzorkovacího teorému (anglicky - sampling theorem) definovaného následující větou (převzato z [5], 14 věta 1). Věta 1: Jestliže pro Fourierovu transformaci F (X) funkce f (x) platí F (X) = 0,
pro |X| ≥ X0 > 0,
(1)
je funkce f (x) jednoznačně určena svými hodnotami v bodech xn =
nπ , kX0
n = ..., −2, −1, 0, 1, 2, ... a platí
f (x) =
∞ X n=−∞
f
nπ kX0
sin (kX0 x − nπ) . kX0 x − nπ
(2)
X je proměnná ve Fourierově obraze a X0 je tedy nejvyšší hodnota této proměnné. Libovolnou funkci, jejíž Fourierova transformace je nenulová pouze v konečném intervalu, tedy můžeme vyjádřit hodnotami v bodech vzorkování xn =
nπ . kX0
Absence vysokých hodnot
proměnné X0 zajistí, že funkce f (x) není mezi body vzorkování xn příliš oscilující. Důkaz věty 1 je proveden v [5], str. 119. Body vzorkování xn =
nπ , kX0
funkce f (x), jsou od sebe vzdáleny intervalem (resp. periodou
vzorkování) T =
π . kX0
(3)
V případě optické difrakce vyjadřuje proměnná X prostorovou frekvenci. Potom, dosadímeli do (3) k =
2π , λ
využijeme vztahu
X λ
= Xf , kde Xf je prostorová frekvence (viz [5], kap. 13
5) a označíme
X0 λ
= Xf 0 , tj. podíl maximální hodnoty směrového kosinu X0 a vlnové délky
λ je roven maximální prostorové frekvenci Xf 0 , můžeme vztah (3) přepsat do tvaru T =
1 . 2Xf 0
(4)
Můžeme tedy říci, že maximální perioda vzorkování je rovna převrácené hodnotě dvojnásobku nejvyšší prostorové frekvence Xf 0 . Toto platí i v jiných případech, např. ve vztahu čas-frekvence. Věta 1 o vzorkování funkce f (x) předpokládá znalost Fourierovy transformace (resp. maximální hodnoty směrového kosinu X0 , kdy je Fourierova transformace ještě nenulová) funkce f (x) pro stanovení velikosti vzorkovací periody T . V praxi se však setkáváme spíše s opačnou úlohou, kdy je potřeba určit periodu vzorkování ještě před tím, než známe vlastní funkci f (x). Např. když funkce f (x) má reprezentovat průběh akustického signálu a je potřeba ho zaznamenat. V takovém případě je nutné zvážit jaké hodnoty frekvence (v případě optické difrakce - hodnoty směrového kosinu) jsou pro nás užitečné. Například maximální hodnota směrového kosinu X = 1 (ve skutečnosti se volí maximální hodnota směrového kosinu mnohem menší, neboť by vzorkování bylo příliš jemné). V případě akustického signálu, pokud se jedná o hodnoty frekvencí, které zaznamená lidské ucho, je maximální hodnota frekvence f0 = 22 kHz a perioda vzorkování T = 2, 3 · 10−5 s viz (4). Ve větě 1 o vzorkování funkce f (x) nemusí být rovnost (3) splněna. Jestliže T <
π , kX0
je
funkce f (x) určena více hodnotami, než je potřeba, ale funkce f (x) je určena jednoznačně. V praxi se využívá opačného postupu pro interpolaci funkčních hodnot již navzorkované funkce s periodou vzorkování (3). Ukážeme si to v kapitole 6.6 o aplikacích diskrétní Fourierovy transformace. Jestliže T >
π , kX0
není funkce f (x) určena jednoznačně, tj. není určena dostateč-
ným počtem hodnot. V anglické literatuře se tato skutečnost nazývána ”aliasing effect” (volný překlad efekt jinakosti). Definice ”aliasing” efektu říká, že je to neschopnost rozlišit frekvenci funkce sinus, kterou vzorky vyjadřují (viz [2], Example 5.1). Příklad ”aliasing” efektu je zobrazen v obr. 3. Plnou čarou je zobrazena funkce f (x) a čárkovaně je zobrazena funkce, kterou je také možné proložit danými body vzorkování. Rovnost (3) tedy určuje maximální periodu vzorkování, aby mohla být funkce f (x) vyjádřena přesně v bodech vzorkování xn =
nπ . kX0
V anglické leteratuře je tato perioda
nazývána jako ”Nyquist sampling rate” (Nyquistova vzorkovací perioda). 14
Obrázek 3: Příklad ”aliasing” efektu, adaptováno z [2]. Analogicky k větě 1 o vzorkování funkce f (x), lze vyslovit větu pro vzorkování Fourierovy transformace F (X) funkce f (x) (převzato z [5], 14 věta 2). Využití je např. pro elementární buňku krystalu (mřížku) s mřížkovým parametrem a, kterou lze v jednorozměrném případě definovat funkcí f (x), jež nabývá nenulových hodnot pouze pro 0 < x < a. Věta 2: Jestliže pro funkci f (x) platí pro |x| ≥ x0 > 0,
f (x) = 0,
je funkce F (X) jednoznačně určena svými hodnotami v bodech Xn =
(5) nπ , kx0
n = ..., −2, −1, 0, 1, 2, ... a platí
F (X) =
∞ X n=−∞
F
nπ kx0
sin (kXx0 − nπ) . kXx0 − nπ
(6)
Důkaz věty 2 je proveden v [5], str. 120. Vzorkovaná Fourierova transformace funkce f (x) má analogické vlastnosti jako vzorkovaná funkce f (x).
6.2
Definice diskrétní Fourierovy transformace
Diskrétní Fourierovu transformaci funkce f (x) lze definovat vzhledem ke spojité Fourierově transformaci. V minulé kapitole jsme učinili první krok k definování diskrétní Fourierovy transformace, a tím je navzorkování funkce f (x). Věta 1 o vzorkování funkce f (x) říká, že funkci f (x) lze vyjádřit hodnotami v bodech vzorkování xl =
lπ . kX0
Tohoto faktu
využijeme při numerickém výpočtu jednorozměrného Fourierova integrálu Z
∞
FT{f (x)} = F (X) = A
f (x) exp(−ikXx) dx −∞
15
(1)
(viz 5(1)) obdélníkovou metodou, kde spodní strana obdélníku je rovna periodě vzorkování T (viz 6.1(3)). Pro numerický výpočet integrálu (1) platí ∞ X . F (X) = AT f (lT ) exp(−ikXlT ) ,
(2)
l=−∞
kde T =
π kX0
(viz 6.1(3)).
V následujícím kroku je nutné omezit definiční obor funkce f (x), neboť je obtížné spočíst nekonečnou sumu ve výrazu (2). Definiční obor funkce f (x) omezíme symetricky tak, že ponecháme pouze N (N je sudé číslo, neboť N/2 musí být celé číslo) po sobě jdoucích prvků (resp. vzorků) od indexu −N/2 počínaje. Omezením definičního oboru funkce f (x) se Fourierova transformace F (X) (viz (2)) změní. Pro Fourierovu transformaci funkce f (x), které jsme omezili definiční obor, ponecháme označení F (X) a platí N/2−1
X
F (X) = AT
f (lT ) exp(−ikXlT ) .
(3)
l=−N/2
Aby platila fundamentální věta o Fourierově transformaci (viz 5(4)), je nutné, aby diskrétní Fourierova transformace měla stejný počet N vzorků jako funkce f (x), u které jsme omezili definiční obor. Jelikož funkce f (lT ) splňuje podmínku věty 2 o vzorkování Fourierovy transformace F (X) (viz kap. 6.1), stačí vyjádřit funkci F (X) pouze v bodech vzorkování Xn =
2πn . kN T
Vztah (3) přepíšeme do tvaru F
2πn kN T
N/2−1
= AT
X
f (lT ) exp(−i2πln/N ) ,
(4)
l=−N/2
kde l = −N/2, −N/2+1, ..., N/2−1 a n = −N/2, −N/2+1, ..., N/2−1, Vztah (4) nazveme diskrétní Fourierova transformace funkce f (x) (resp. funkce f (lT ), l = −N/2, −N/2 + 1, ..., N/2 − 1).
6.3
Inverzní diskrétní Fourierova transformace (IDFT)
Obdobně definujeme inverzní diskrétní Fourierovu transformaci 2π f (lT ) = B kN T
N/2−1
X
F
n=−N/2
16
2πn kN T
exp(i2πln/N ),
(1)
kde l = −N/2, −N/2 + 1, ..., N/2 − 1 a n = −N/2, −N/2 + 1, ..., N/2 − 1. Opět se jedná o numerický výpočet jednorozměrného Fourierova integrálu 5(2) obdélníkovou metodou a omezení definičního oboru funkce na konečný počet vzorků N . Důkaz se provede dosazením vztahu (1) do 6.2(4) F
2πn kN T
N/2−1
=
X
AT
l=−N/2
1 N
=
F
m=−N/2
=
B 2π kN T
N/2−1
X
F
2πn kN T
N/n−1
X
F
m=−N/2
2πm kN T
2πm kN T
exp(i2πlm/N ) exp(−i2πln/N )
N/2−1
X
exp(i2πlm/N ) exp(−i2πln/N )
(2)
l=−N/2
.
V důkazu je využito vztahu 5(3) a výpočtu sumy N/2−1
X n=−N/2
N, m = n, exp(i2πlm/N ) exp(−i2πln/N ) = 0, m 6= n.
(3)
To, že suma (3) je rovna nule pro m 6= n, se dokáže vztahem pro součet N členů geometrické řady.
6.4
Modifikace diskrétní Fourierovy transformace
Pro výpočet diskrétní Fourierovy transformace pomocí rychlé Fourierovy transformace, (o rychlé Fourierově transformaci pojednává kapitola 7) je výhodné, aby indexy l a n v definičním vztahu 6.2(4) pro diskrétní Fourierovu transformaci (resp. ve vztahu 6.3(1) pro inverzní diskrétní Fourierovu transformaci) začínaly nulou. Definujme proto funkci h(lT ) takto f (lT ) , h(lT ) = f ((l − N )T ) ,
l = 0, 1, ..., N/2 − 1,
(1)
l = N/2, N/2 + 1, ..., N − 1,
tj. prvky se zápornými indexy l přesuneme za prvky s kladnými indexy. 2πn takto Dále definujme funkci H kN T H
2πn kN T
=
F F
2πn kN T
2π(n−N ) kN T
,
n = 0, 1, ..., N/2 − 1,
,
n = N/2, N/2 + 1, ..., N − 1,
(2)
17
tj. prvky se zápornými indexy n přesuneme za prvky s kladnými indexy. Pro diskrétní Fourierovu transformaci funkce h(lT ) platí H
2πn kN T
= AT
N −1 X
h(lT ) exp(−i2πln/N ) ,
(3)
l=0
kde l = 0, 1, ..., N − 1 a n = 0, 1, ..., N − 1. Pro inverzní diskrétní Fourierovu transformaci funkce H
2πn kN T
platí
N −1 2π X 2πn h (lT ) = B H exp(i2πln/N ), kN T n=0 kN T
(4)
kde l = 0, 1, ..., N − 1 a n = 0, 1, ..., N − 1. Abychom dostali zpět funkci f (lT ) (resp. funkci F
2πn kN T
), je nutné přesunout zpět prvky
odpovídající záporným indexům. V literatuře obecně se neupozorňuje na to, že je nutné přesunout prvky se zápornými indexy za prvky s kladnými indexy a počítá se přímo diskrétní Fourierova transformace nemodifikované funkce f (lT ) pomocí vztahu (3). Nepřesunutí prvků způsobí, že prvky s lichými indexy budou mít stejnou hodnotu, jako kdybychom počítali diskrétní Fourierovu transformaci modifikované funkce h(lT ) pomocí vztahu (3), ovšem se znaménkem mínus. Hodnoty sudých prvků zůstanou nezměněny. Názorně to bude ukázáno na obr. 5 v dalším textu. Přesunutí prvků se zápornými indexy za prvky s kladnými indexy vstupní funkce, tj. funkce, pro kterou chceme spočíst diskrétní Fourierovu transformaci (resp. inverzní diskrétní Fourierovu transformaci), nemá vliv na kvadrát modulu diskrétní Fourierovy transformace (resp. inverzní diskrétní Fourierovy transformace). Dále se budeme zabývat pouze funkcí h(lT ) (resp. funkcí H
6.5
2πn kN T
).
Vlastnosti diskrétní Fourierovy transformace
V kapitole 5 je pro Fourierovu transformaci užito označení FT (resp. FT−1 pro inverzní Fourierovu transformaci). Pro diskrétní Fourierovu transformaci použijeme označení DFT (resp. DFT−1 pro inverzní diskrétní Fourierovu transformaci). Vlastnosti diskrétní Fourierovy transformace jsou analogické vlastnostem spojité Fourierovy transformace (viz [5] nebo [7]).
18
6.5.1
Linearita diskrétní Fourierovy transformace
Linearitu vyjadřuje rovnice analogická rovnici pro linearitu spojité Fourierovy transformace (viz [5], 6.1(1))
DFT
( X
) αj hj (lT )
=
j
X
αj DFT{hj (lT )},
(1)
j
kde αj jsou konstanty (mohou být i komplexní). Důkaz vztahu (1) je založen na záměně pořadí sčítání. 6.5.2
Diskrétní Fourierova transformace a posunutí funkce
Jestliže pro funkci h(lT ) platí DFT{h(lT )} = H
2πn kN T
,
(2)
je DFT{h(lT − sT )} = exp(−i2πns/N )H
2πn kN T
,
(3)
kde sT je posunutí, T je perioda vzorkování a s je celé číslo. Diskrétní Fourierovu transformaci posunutí funkce h(lT ) o jinou hodnotu než je celočíselný násobek s periody vzorkování T , by bylo možné spočítat, ovšem vedlo by to ke komplikacím při aplikaci inverzní diskrétní Fourierovy transformace na diskrétní Fourierovu transformaci posunuté funkce h(lT ). Proto hodnotu jiného posunutí funkce zaokrouhlíme na celočíselný násobek periody vzorkování. To, že velikost posunutí by měla být rovna celočíselnému násobku periody, se v literatuře poněkud obchází, např v [6] se hovoří pouze o posunutí posloupnosti (tj. změně indexu) a její diskrétní Fourierově transformaci. V [2] (kap. 6.5) je napsáno, že body vzorkování lT nahradíme pouze indexem l pro výhodnější značení. Pro důkaz tvrzení (3) se využije definice inverzní diskrétní Fourierovy transformace 6.3(1)
19
h (lT )
h (lT − sT )
=
N −1 2π X 2πn B H exp(i2πln/N ) kN T n=0 kN T
=
N −1 2πn 2π X H exp(i2π(l − s)n/N ) B kN T n=0 kN T
=
N −1 2π X 2πn B H exp(−i2πsn/N ) exp(i2πln/N ). kN T n=0 kN T
= 6.5.3
DFT
−1
(4)
2πn exp(−i2πns/N )H . kN T
Diskrétní Fourierova transformace a posunutí Fourierova obrazu
Jestliže pro funkci h(lT ) platí DFT{h(lT )} = H
2πn kN T
,
(5)
je DFT{h(lT ) exp(i2πsl/N )} = H
2π(n − s) kN T
,
(6)
kde s je celé číslo. Diskrétní Fourierova transformace je posunutá o celočíselný násobek s výrazu
2π , kN T
6.5.4
Diskrétní Fourierova transformace rozšířené funkce o nulové hodnoty
který určuje vzdálenost mezi body vzorkování Fourierovy transformace.
Mějme funkci h(lT ), l = 0, 1, ..., N − 1. Definujme funkci y(sT ), s = 0, 1, ..., M − 1, kde M a N jsou sudá kladná čísla a platí M > N , takto h(sT ), s = 0, 1, ..., N − 1, y(sT ) = 0, s = N, N + 1, ..., M − 1.
(7)
Funkce y(sT ) je tedy rozšířená funkce h(lT ) o M − N nulových hodnot. Pro diskrétní Fourierovu transformaci funkce y(sT ) platí
20
DFT{y(sT )} = Y
2πn kM T
=
AT
M −1 X
y (sT ) exp(−i2πns/M )
s=0
=
AT
N −1 X
h (sT ) exp(−i2πns/M )
(8)
s=0
=
AT
N −1 X
h (sT ) exp(−i2πn(N/M )s/N )
s=0
=
H
2πn(N/M ) kN T
.
V bodech, kde pro index n platí, že nN/M je celé číslo, n = 0, 1, ..., M − 1, je diskrétní Fourierova transformace funkce y(sT ) rovna hodnotám diskrétní Fourierovy transformace funkce h(lT ) v bodech, jejichž indexy l jsou rovny nN/M . O hodnotách v bodech, kde neplatí, že n(N/M ) je celé číslo, pojednává kapitola 6.6.1 o interpolaci funkčních hodnot. Analogicky k výpočtu diskrétní Fourierovy transformace rozšířené funkce h(lT ) se 2πn spočítá i inverzní diskrétní Fourierova transformace funkce H kN rozšířené o nulové T hodnoty.
6.6 6.6.1
Aplikace diskrétní Fourierovy transformace Interpolace funkčních hodnot pomocí diskrétní Fourierovy transformace
V kapitole 6.1 o vzorkovacím teorému jsme se zmínili, že pro periodu vzorkování T nemusí být splněna rovnost 6.1(3). Zmínili jsme se také o využití pro interpolaci funkčních hodnot. 2πn Mějme funkci h(lT ), její diskrétní Fourierovu transformaci H kN a periodu vzorkoT vání T , pro kterou platí T =
π , kX0
(1)
viz 6.1(3), l = 0, 1, ...N − 1 a n = 0, 1, ...N − 1, N je sudé kladné číslo. 2πs Definujme funkci Y kN pro kterou platí T Y
2πs kN T
H = 0,
2πs kN T
,
s = 0, 1, ..., N − 1,
s = N, N + 1, ..., pN − 1,
21
(2)
kde p = 2, 3, 4, ..., tj. funkce Y
2πs kN T
je rovna funkci H
2πn kN T
rozšířené o (p − 1)N
nulových hodnot. Potom pro inverzní diskrétní Fourierovu transformaci y(jT 0 ), kde T 0 = 2πs platí T /p je perioda vzorkování funkce y(jT 0 ), funkce Y kN T 0
y (jT )
−1
pN −1 2πs 2π X 2πs Y = B Y exp(i2πjs/pN ) kN T kN T s=0 kN T
=
DF T
=
N −1 2πs 2π X H exp(i2πjs/pN ) B kN T s=0 kN T
=
h(jT 0 ), j = pl, DF T −1 Y 2πs exp(i2πsm/pN ) , kN T
(3)
j 6= pl,
kde m je celé číslo, pro než platí 1 ≤ m ≤ p − 1. Číslo m určuje počet interpolovaných hodnot mezi jednotlivými body vzorkování funkce h(lT ), které získáme aplikací inverzní 2πs diskrétní Fourierovy transformace na rozšířenou funkci Y kN . T Formálně jsou funkční hodnoty funkce y (jT 0 ) v bodech, pro něž hodnoty indexu nesplňují podmínku j = pl, rovny hodnotám funkce h(lT ) posunuté o m/p násobek periody vzorkování T (resp. o vzdálenost mT 0 ). Analogicky interpolaci funkčních hodnot funkce 2πn h(lT ) lze interpolovat i funkční hodnoty Fourierovy transformace H kN funkce h(lT ). T Rozšíření funkčních hodnot funkce h(lT ) (resp. Fourierovy transformace funkce h(lT ) rozšířené o nulové hodnoty (viz kap. 6.5.4)) nemá vliv na Fourierovu (resp. inverzní Fourierovu) transformaci. Při rozšíření o nulové hodnoty dojde pouze k interpolaci funkčních hodnot. 6.6.2
Diskrétní konvoluce periodických funkcí
Diskrétní konvoluce h(lT ) = h1 (lT ) ∗ h2 (lT ) dvou periodických funkcí h1 (lT ), h2 (lT ) se stejnou periodou N (resp. se stejným počtem N vzorků v jedné periodě), kde l = 0, 1, ..., N − 1, N je sudé kladné číslo, je definována sumou (viz [2], 6.5(6.49)) h(lT ) = h1 (lT ) ∗ h2 (lT ) = T
N −1 X
h1 (jT )h2 ((l − j)T ).
(4)
j=0
Funkce h1 (lT ) a h2 (lT ) musí mít stejný počet vzorků v jedné periodě a musí být navzorkovány stejnou periodou vzorkování T . Vztah (4) je po formální stránce vztahem pro numerický výpočet konvolučního integrálu (viz [5], 9(1)). 22
6.6.3
Diskrétní konvoluce „neperiodickýchÿ funkcí definovaných na intervalu konečné délky
Tato kapitola pojednává o diskrétní konvoluci funkcí h1 (mT1 ) a h2 (nT2 ), které nejsou periodické a jsou nenulové pouze na intervalu konečné délky. Slovo neperiodických v názvu kapitoly je v uvozovkách záměrně, protože se jedná opět o diskrétní konvoluci periodických funkcí, ovšem s tím rozdílem, že nejprve musíme upravit funkce h1 (mT1 ) a h2 (nT2 ) tak, aby byly splněny požadavky pro konvoluci periodických funkcí (viz kap. 6.6.2) a teprve poté použít definičního vztahu 6.6.2(4) pro diskrétní konvoluci periodických funkcí. Úpravu funkcí h1 (mT1 ) a h2 (nT2 ) lze rozdělit do dvou fází. První fázi lze přeskočit, jestliže pro periody vzorkování funkcí h1 (mT1 ) a h2 (nT2 ) platí T1 = T2 = T . Jestliže pro funkce h1 (mT1 ) a h2 (nT2 ) neplatí T1 = T2 , tj. funkce h1 (mT1 ) a h2 (nT2 ) nemají stejnou periodu vzorkování, ale existuje číslo T ∈ (0, ∞) pro které platí T1 = sT, T2 = rT,
(5)
kde s, r jsou celá kladná čísla, lze pomocí interpolace funkčních hodnot (viz kap. 6.6.1), uvést funkce h1 (mT1 ) a h2 (nT2 ) na společnou periodu vzorkování T . Interpolace způsobí změnu počtu vzorků N1 (resp. N2 ). V druhé fázi rozšíříme funkci h1 (mT ) (resp. funkci h2 (nT )) o nulové hodnoty (viz kap. 6.5.4) tak, aby pro nový počet vzorků N funkce h1 (mT ) (resp. funkce h2 (nT )) platilo (viz [2], 7.3(7.6)) N = N1 + N2 − 1.
(6)
Jelikož při konvoluci počítáme plochu pod křivkou součinu dvou funkcí, z nichž jedna je zrcadlově převrácena kolem osy y a přesouvá se přes celý obor hodnot, je nutné počet vzorků N volit tak, aby při aplikaci diskrétní konvoluce existovalo i jednoprvkové překrytí funkcí h1 (mT ) a h2 (nT ). Nyní máme upravené funkce h1 (nT ) a h2 (nT ), které mají stejnou periodu vzorkování a stejný počet vzorků. Každá z těchto funkcí může tvořit jednu periodu periodické funkce, a proto lze použít definičního vztahu (4) pro diskrétní konvoluci periodických funkcí.
23
6.6.4
Fourierova transformace diskrétní konvoluce a diskrétní Fourierova transformace součinu
Nechť DFT{h1 (lT )} = H1
2πn kN T
a DFT{h2 (lT )} = H2
2πn kN T
. Pak
(i) 1 DFT{h1 (lT ) ∗ h2 (lT )} = H1 A
2πn kN T
H2
2πn kN T
(7)
a tedy
DFT
−1
2πn 2πn H1 H2 = A h1 (lT ) ∗ h2 (lT ). kN T kN T
(8)
(ii) DFT{h1 (lT )h2 (lT )} = B H1
2πn kN T
∗ H2
2πn kN T
(9)
a tedy
DFT
−1
H1
2πn kN T
∗ H2
2πn kN T
=
1 h1 (lT )h2 (lT ). B
(10)
Tvrzení i jeho důkaz je analogický spojité Fourierově transformaci konvoluce a spojité Fourierově transformaci součinu (viz [5], kap. 9.3) s rozdílem, že integrace je nahrazena sumací.
6.7
Příklady výpočtu spojité Fourierovy transformace pomocí diskrétní Fourierovy transformace
Diskrétní Fourierova transformace se liší od spojité Fourierovy transformace. Existuje však skupina funkcí, pro které platí, že diskrétní Fourierova transformace je totožná se spojitou Fourierovou transformací. Rovnost mezi diskrétní a spojitou Fourierovou transformací platí pro periodické funkce (resp. funkce definované na konečném intervalu), pro jejichž periodu vzorkování platí vztah 6.1(4) a jsou omezeny intervalem délky N T , která je rovna celočíselnému násobku periody funkce (resp. délka N T je větší nebo rovna délce intervalu, na kterém je funkce definována).
24
6.7.1
Příklad: h(lT ) = cos(alT ),
kde a =
2π NT
2π ). aN
(resp. T =
Pro diskrétní Fourierovu transformaci funkce cos H
2πn kN T
= AT
N −1 X l=0
2πl N
platí
2πl exp(−i2πln/N ) . cos N
(1)
Použitím Eulerovy věty přepíšeme vztah (1) do tvaru H
2πn kN T
N −1 AT X −i2πl i2πl = + exp exp(−i2πln/N ) , exp 2 l=0 N N
(2)
který dále upravíme do tvaru H
2πn kN T
N −1 AT X −i2πl(n − 1) −i2πl(n + 1) = exp + exp . 2 l=0 N N
(3)
Využitím vztahu pro výpočet sumy N členů geometrické řady dostaneme H
2πn kN T
AT = 2
exp[−i2π(n − 1)] − 1 exp[−i2π(n + 1)] − 1 + exp[−i2π(n − 1)/N ] − 1 exp[−i2π(n + 1)/N ] − 1
.
(4)
Vztah (4) je nenulový pouze pro n = 1 (resp. n = N −1). V těchto bodech musíme spočíst limitu lim H n→1
2πn kN T
=
exp[−i2π(n − 1)] − 1 + exp[−i2π(n − 1)/N ] − 1 exp[−i2π(n + 1)] − 1 AT N + , = exp[−i2π(n + 1)/N ] − 1 2 AT lim n→1 2
(5)
a limitu lim
n→N −1
H
2πn kN T
=
AT N . 2
(6)
Úvaha: Řekněme, že bod s indexem n = b diskrétní funkce, který sousedí pouze s nulovými hodnotami, reprezentuje delta funkci δ(n − b). Potom přepíšeme vztah (4), s přihlédnutím k hodnotě limit (5) a (6) do tvaru H
2πn kN T
=
AT N [δ(n − 1) + δ(n − N + 1)] . 2 25
(7)
Přejdeme nyní k funkci F
2πn kN T
přesunutím hodnot, které odpovídají záporným indexům
(viz kap. 6.4) F Substitucí
2πn kN T
2πn kN T
=
AT N [δ(n − 1) + δ(n + 1)] . 2
(8)
= X dostaneme AT N kN T kN T F (X) = δ X −1 +δ X +1 . 2 2π 2π
Využitím podmínky AB =
|k| 2π
(9)
(viz 5(3)) a vztahu pro delta funkci
δ(ax − x0 ) =
1 x0 δ(x − ), |a| a
(10)
(viz [5], 1.3 (7)) upravíme vztah (9) do konečného tvaru 1 2πn 2πn H(X) = δ X− +δ X + . 2B kN T kN T Tento výsledek dostaneme i aplikací spojité Fourierovy transformace 5(1) na funkci cos
(11) 2π x NT
(viz [5], 6.2.1). Tento fakt potvrzuje úvodní odstavec kapitoly 6.7. 6.7.2
Příklad: h(lT ) = rect
Funkce rect
lT 2a
lT 2a
je definována takto 1, lT = rect 0, 2a
Před vlastním výpočtem musíme funkci rect
lT 2a
|lT | < a,
(12)
|lT | > a.
modifikovat tak, že přesuneme prvky se
zápornými indexy za prvky s kladnými indexy (viz kap. 6.4). V obrázku 5 je zobrazen výsledek výpočtu diskrétní Fourierovy transformace pomocí vztahu 6.4(3), kdybychom lT lT nemodifikovali funkci rect 2a . Ihned je zřejmé, proč je nutné modifikovat funkci rect 2a , neboť u prvků s lichými indexy se objeví znaménko mínus. V obrázku 6 je zobrazen lT výsledek výpočtu diskrétní Fourierovy transformace modifikované funkce rect 2a . Počet vzorků N = 512, šířka 2a = 14T .
26
Obrázek 4: Graf funkce rect
2πn kN T
Obrázek 5: Graf funkce H 6.7.3
Příklad: rect
lT 2a
∗ rect
Pro výpočet konvoluce rect
lT 2a
lT 2b
lT 2a
.
spočtené pro nemodifikovanou funkci rect
lT 2a
.
∗ rect
lT 2b
využijeme tvrzení (i) z kapitoly 6.6.4. Výpočet
diskrétní konvoluce provedeme tak, že nejprve spočítáme diskrétní Fourierovu transfor lT maci funkcí rect 2a a rect lT . Spočtené diskrétní Fourierovy transformace obou funkcí 2b mezi sebou vynásobíme a na součin aplikujeme inverzní diskrétní Fourierovu transformaci. lT V obrázku 7 je zobrazena výsledná diskrétní konvoluce rect 2a ∗ rect lT , pro N = 512, 2b 2a = 200T a 2b = 310T .
27
Obrázek 6: Graf funkce H
2πn kN T
spočtené pro modifikovanou funkci rect
Obrázek 7: Graf konvoluce rect
7
lT 2a
∗ rect
lT 2b
lT 2a
.
.
Rychlá Fourierova transformace (FFT)
Existuje mnoho algoritmů pro zrychlení výpočtu diskrétní Fourierovy transformace (viz [2], kap. 8). V této kapitole se budeme zabývat algoritmem pro zrychlení výpočtu diskrétní Fourierovy transformace funkce h(lT ), jejíž počet vzorků je N = 2γ , kde γ je celé kladné číslo.
7.1
Maticové vyjádření diskrétní Fourierovy transformace
Vztah 6.3(3) H
2πn kN T
= AT
N −1 X
h (lT ) exp(−i2πln/N ),
l=0
28
(1)
pro výpočet diskrétní Fourierovy transformace, lze zapsat pomocí matice (viz [2], 8.1(8.4)). Použijeme označení
H
W 2πn kN T h (lT )
=
exp(−i2π/N ),
(2)
=
X(n),
(3)
=
x(l),
(4)
kde X(n) a x(l) chápeme jako posloupnosti o N prvcích. Tyto posloupnosti zapíšeme do ~ a ~x. sloupcových vektorů X Vztah (1) zapíšeme následovně X(0)
= x(0)W 0 +
x(1)W 0
+
x(2)W 0
+ · · · + x(N − 1)W 0
X(1)
= x(0)W 0 +
x(1)W 1
+
x(2)W 2
+ · · · + x(N − 1)W N −1
X(2) .. .
= x(0)W 0 + .. = .
x(1)W 2
+
x(2)W 4
+ · · · + x(N − 1)W 2(N −1)
X(N − 1) = x(0)W 0 + x(1)W N −1 + x(2)W 2(N −1) + · · · + x(N − 1)W (N −1)(N −1) , (5) nebo maticově W0 X(0) W0 X(1) X(2) = W 0 .. .. . . W0 X(N − 1)
W0
W0
1
2
W
W2
W
··· W0 ··· W
N −1
· · · W 2(N −1)
W4
W N −1 W 2(N −1) · · · W (N −1)(N −1)
x(0) x(1) x(2) .. . x(N − 1)
. (6)
Výpočet konstanty W s (s = 0, 1, 2, ..., (N −1)(N −1)) lze zjednodušit, neboť platí (viz [2], 8.2(8.7)) Ws = Ws
mod(N )
,
(7)
kde mod(N ) je zbytek po celočíselném dělení číslem N . Zápis (6) lze zapsat i následovně ~ = W~x, X
(8)
kde W je čtvercová matice N × N , jejíž prvky wi,j tvoří konstanty W ij . Nyní můžeme nahlédnout, kolik matematických operací je potřeba k výpočtu diskrétní Fourierovy transformace. Konstanta W (viz (2)) je komplexní, funkce h(lT ) je obecně komplexní, potom 29
výpočet diskrétní Fourierovy transformace vyžaduje N 2 operací komplexního násobení a N (N − 1) operací komplexního sčítání, tzn. celkem 2N (N − 1) komplexních operací. V následující kapitole si ukážeme, že toto množství operací není potřeba.
7.2
Rychlá Fourierova transformace
V roce 1965 vytvořili J.W. Cooley a J.W. Turkey (viz [1]) algoritmus pro výpočet rychlé Fourierovy transformace (fast Fourier transform - FFT) funkcí, které jsou navzorkovány N = 2γ vzorky, kde γ je celé kladné číslo. Algoritmus využívá metody vybraných posloupností (viz [6], kap. 6.2.1). Rozdělme posloupnost x(l) mající N členů na dvě vybrané posloupnosti o N/2, a to na posloupnost sudých členů a na posloupnost lichých členů. Tento postup bývá nazýván rychlou Fourierovou transformací se základem 2. Pro diskrétní Fourierovu transformaci X(n) posloupnosti x(l) platí
X(n)
=
AT
N −1 X
x(l) exp(−i2πln/N )
l=0 N/2−1
=
AT
X
x(2l) exp(−i2πln/(N/2)) +
(1)
l=0 N/2−1
+ AT
X
x(2l + 1) exp(−i2πln/(N/2)).
l=0
Označíme-li diskrétní Fourierovu transformaci posloupnosti sudých prvků DFT{x(2l)} = X e (n),
n = 0, 1, 2, ..., N/2 − 1,
(2)
a diskrétní Fourierovu transformaci posloupnosti lichých prvků DFT{x(2l + 1)} = X o (n),
n = 0, 1, 2, ..., N/2 − 1,
(3)
pak lze diskrétní Fourierovu transformaci X(n) posloupnosti x(l) vyjádřit dvěma vztahy X(n)
=
X e (n) + W n X o (n),
n = 0, 1, 2, ..., N/2 − 1,
(4)
X(n + N/2)
=
X e (n) − W n X o (n),
n = 0, 1, 2, ..., N/2 − 1,
(5)
30
kde W = exp(−i2π/N ) (viz 7.1(2)). Na rozdělení posloupnosti přišli Danielson a Lanczos (1942) (viz [3]). Rozdělení posloupnosti aplikujeme rekurzivně, tzn. že opět rozdělíme posloupnost X e (n) (resp. posloupnost X o (n)) na dvě posloupnosti výběrem sudých a lichých prvků. Posloupnost dělíme (γ − 1) krát, až máme N/2 dvouprvkových posloupností. Diskrétní Fourierova transformace je na dvouprvkové posloupnosti rovna X(0)
=
x(0) + x(1),
X(1)
=
x(0) − x(1),
viz vztah (4) a (5). Jelikož známe diskrétní Fourierovu tranformaci na N/2 dvouprvkových posloupnostech, použijeme opět vztahu (4) (resp. (5)) a spočteme N/4 diskrétních Fourierových transformací na čtyřprvkových posloupnostech. Toto aplikujeme dále až se dostaneme k N -prvkové diskrétní Fourierově transformaci. Výpočet pomocí rekurze způsobí, že dojde k přeskládání prvků výsledné diskrétní Fourierovy transformace v bitově inverzím pořadí. Bitově invertované pořadí vznikne tak, že indexy prvků vyjádříme binárním číslem a toto binární číslo čteme od konce. Např. číslo k = 3 zapíšeme binárně jako 011. Bitově inverzím pořadí přejde toto číslo na 110, kterému odpovídá index 6. Při výpočtu je výhodné nejprve přeskládat prvky posloupnosti v bitově invertovaném pořadí a teprve poté aplikovat rekurzivně rozdělení diskrétní Fourierovy transformace. Výpočet diskrétní Fourierovy transformace pomocí rekurzivního rozdělení posloupností vyžaduje pouze N log2 N komplexních operací. Diskrétní Fourierova transformace spočtená pomocí algoritmu rychlé Fourierovy transformace se neliší od diskrétní Fourierovy transformace spočtené přímo, tj. pomocí definičním vztahem 6.4(3). Zajímavostí je, že se v literatuře obecně přisuzují rychlé Fourierově transformaci zvláštní vlastnosti, které vnášejí jistou mystičnost do tohoto pojmu. Proto je jistě vhodné zopakovat fakt, že rychlá Fourierova transformace je pouze algoritmem výpočtu diskrétní Fourierovy transfomace a že jakékoli jí přisuzované zvláštní vlastnosti plynou pouze z nepochopení daného problému. Zmiňovat se dále o tomto algoritmu by bylo zbytečné, neboť s velkým rozvojem počítačové techniky v posledních 40 letech, došlo i k velkému rozvoji a propracování algoritmů pro výpočet rychlé Fourierovy transformace. Podrobněji o rychlé Fourierově transformaci pojednává např. [2] (kap. 8), [6] (kap. 6) nebo [3] (kap. 12). Dnes existují nejen pro31
gramové verze algoritmu rychlé Fourierovy transformace, ale i elektronické integrované obvody pro výpočet rychlé Fourierovy transformace. V příloze 1 je uveden zdrojový kód pro algoritmus rychlé Fourierovy transformace pro programovací prostředí Delphi.
8 8.1
Vícerozměrná rychlá Fourierova transformace Dvourozměrná diskrétní Fourierova transformace
Odvození dvourozměrné diskrétní Fourierovy transformace je analogické odvození jednorozměrné diskrétní Fourierovy transformace. Mějme dvourozměrnou komplexní funkci h(l1 T1 , l2 T1 ), kde l1 = 0, 1, 2, ..., N1 − 1, l2 = 0, 1, 2, ..., N2 − 1, kde N1 (resp. N2 ) je počet vzorků v ose x (resp. v ose y). T1 (resp. T2 ) je perioda vzorkování v ose x (resp. v ose y). Potom pro diskrétní Fourierovu transformaci funkce h(l1 T1 , l2 T2 ) platí H
2πn1 2πn2 , kN1 T1 kN2 T2
=
2
A T1 T2
N 1 −1 2 −1 N X X
h(l1 T1 , l2 T2 ) ×
l2 =0 l1 =0
× exp(−i2πl1 n1 /N1 ) exp(−i2πl2 n2 /N2 ),
(1)
kde n1 = 0, 1, 2, ..., N1 − 1 a n2 = 0, 1, 2, ..., N2 − 1. Vztah (1) upravíme do tvaru H
2πn1 2πn2 , kN1 T1 kN2 T2
=
2
A T1 T2
N 2 −1 X
exp(−i2πl2 n2 /N2 ) ×
l2 =0
×
N 1 −1 X
h(l1 T1 , l2 T2 ) exp(−i2πl1 n1 /N1 ),
(2)
l1 =0
neboť fázor exp(−i2πl2 n2 /N2 ) nezávisí na vnitřní sumě. Použijeme-li označení
H
h(l1 T1 , l2 T2 ) 2πn1 2πn2 , kN1 T1 kN2 T2
=
x(l1 , l2 ),
(3)
=
X(n1 , n2 ),
(4)
kde x(l1 , l2 ) (resp. X(n1 , n2 )) je matice x (resp. matice X), můžeme chápat výpočet dvourozměrné diskrétní Fourierovy transformace tak, že nejprve spočteme diskrétní Fourierovu 32
transformaci na všech řádcích matice x a na vzniklé matici spočteme diskrétní Fourierovu transformaci na všech sloupcích a tím dostaneme matici X. Pro výpočet diskrétní Fourierovy transformace použijeme algoritmu rychlé Fourierovy transformace. Pro výpočty pomocí rychlé Fourierovy transformace je výhodné, aby platilo T1 = T2 = T a N1 = N2 = N , tj. aby perioda vzorkování v ose x a ose y byla stejná a počet vzorků v ose x a ose y byl stejný. 8.1.1
Přeskládání kvadrantů diskrétní Fourierovy transformace pro konvenční zobrazení
V kapitole 6.4 jsme původní funkci f (lT ) modifikovali na funkci h(lT ). U dvourozměrné funkce h(l1 T, l2 T ) je pro zpětné zobrazení funkce f (l1 T, l2 T ) nutné zaměnit kvadranty I. za III., a II. za IV. Přeskládání kvadrantů je zobrazeno na příkladu matice 4 × 4 v obr. 8.
Obrázek 8: Přeskládání kvadrantů matice pro konvenční zobrazení výpočtu dvourozměrné diskrétní Fourierovy transformace. a) vypočtená matice, b) přeskládaná matice. Převzato z [2], str. 245.
33
8.2
Vícerozměrná diskrétní Fourierova transformace
Mějme p-rozměrnou komplexní funkci h(l1 T1 , T2 l2 , ..., lp Tp ), kde l1 = 0, 1, 2, ..., N1 − 1, l2 = 0, 1, 2, ..., N2 − 1, . . ., lp = 0, 1, 2, ..., Np − 1, kde N1 , N2 , . . ., Np jsou počty vzorků v osách x1 , x2 , . . ., xp . T1 , T2 , . . ., Tp jsou periody vzorkování v osách x1 , x2 , . . ., xp . Potom pro diskrétní Fourierovu transformaci funkce h(l1 T1 , T2 l2 , ..., lp Tp ) platí H
2πn1 2πn2 2πnp , , ..., kN1 T1 kN2 T2 kNp Tp
Np −1 p
= A T1 T2 · · · Tp
X
...
lp =0
N 2 −1 N 1 −1 X X
h(l1 T1 , T2 l2 , ..., lp Tp ) ×
l2 =0 l1 =0
× exp(−i2πl1 n1 /N1 ) exp(−i2πl2 n2 /N2 ) · ... · exp(−i2πlp np /Np ), (5) kde n1 = 0, 1, 2, ..., N1 − 1, n2 = 0, 1, 2, ..., N2 − 1, . . ., np = 0, 1, 2, ..., Np − 1. Vztah (5) je uveden pro úplnost. V této práci se budeme zabývat pouze dvourozměrnou rychlou Fourierovou transformací. Dvourozměrná (resp. vícerozměrná) diskrétní Fourierova transformace má analogické vlastnosti (resp. aplikace) jako jednorozměrná diskrétní Fourierova transformace (viz kap. 6.5 a kap. 6.6). Pokud bude v dalším textu řeč o rychlé Fourierově transformaci, bude tím míněna diskrétní Fourierova transformace spočtená pomocí algoritmu rychlé Fourierovy transformace.
9
Výpočet Fraunhoferových a Fresnelových difrakčních jevů pomocí rychlé Fourierovy transformace
9.1
Postup výpočtu rychlé Fourierovy transformace
Na začátku je nejprve nutné vysvětlit postup výpočtu rychlé Fourierovy transformace pro dvourozměrné funkce. Výpočet rozdělíme do několika kroků: 1. Zadání funkce f (x, y). 2. Funkci f (x, y) navzorkujeme v ose x a y stejnou periodou vzorkování T a vyjádříme ji pouze v bodech vzorkování (x, y) = (l1 T, l2 T ). 3. Definiční obor funkce f (l1 T, l2 T ) omezíme na dvourozměrný konečný interval h−N T /2; (N/2 − 1)T i × h−N T /2; (N/2 − 1)T i, počet vzorků v jedné ose N = 2γ (γ je kladné celé číslo). 34
4. Funkci f (l1 T, l2 T ) modifikujeme na funkci h(m1 T, m2 T ), kde m1 = 0, 1, ..., N − 1 a m2 = 0, 1, ..., N − 1, tj. zaměníme kvadranty I. za III. a II. za IV. 5. Body funkce h(m1 T, m2 T ) ztotožníme s prvky matice X, tj. xm1 ,m2 = h(m1 T, m2 T ). 6. Spočteme jednorozměrnou rychlou Fourierovu transformaci na všech řádcích matice X, tím vznikne matice X1 . 7. Spočteme jednorozměrnou rychlou Fourierovu transformaci na všech sloupcích matice X1 , tím vznikne matice Y. 8. Prvky ys1 ,s2 matice Y ztotožníme s prvky funkce H kde s1 = 0, 1, ..., N − 1 a s2 = 0, 1, ..., N − 1, tj. H
2πs1 2πs2 , kN T kN T 2πs1 , kN T kN T
2πs2
,
= ys1 ,s2 .
2πs1 2πs2 9. zaměníme kvadranty I. za III. a II. za IV. funkce H kN , tím vznikne funkce , T kN T 1 2πn2 F 2πn , , kde n1 = 0, 1, ..., N − 1 a n2 = 0, 1, ..., N − 1, která je rychlou FoukN T kN T rierovou transformací funkce f (x, y). Jestliže se jedná o výpočty difrakčních jevů, potřebujeme navíc spočíst intenzitu světla. 2 1 2πn2 1 2πn2 , , . Na závěr tedy spočteme kvadrát modulu funkce F 2πn , tj. F 2πn kN T kN T kN T kN T
9.2
Výpočet Fraunhoferovy difrakce pomocí rychlé Fourierovy transformace
V kapitole 3 bylo řečeno, že Fraunhoferovu difrakci lze spočíst jako Fourierovu transformaci funkce propustnosti. Difrakční integrál pro Fraunhoferovu difrakci je tvaru (viz 3(6)) Z∞Z a(X, Y ) =
t(x, y)ψ0 (x, y) exp[−ik(Xx + Y y)]dxdy,
(1)
−∞
kde t(x, y) je funkce propustnosti a ψ0 (x, y) je vlnová funkce osvětlující vlny. Při výpočtech se budeme zabývat pouze případy kdy stínítko osvětlujeme rovinou vlnou b exp(−ikz). Jelikož nás zajímá intenzita difraktovaného světla, stačí spočíst pouze Fourierovu transformaci funkce propustnosti, neboť | exp(−ikz)|2 = 1. Konkrétní hodnotu amplitudy b osvětlující vlny nemusíme uvažovat, neboť v závěru výpočtu budeme zobrazovat pouze relativní intenzitu difraktovaného světla vůči maximální hodnotě intenzity difraktovaného světla. 35
Difrakční integrál pro Fraunhoferovu difrakci spočteme tak, že spočítáme rychlou Fourierovu transformaci funkce propustnosti t(x, y). Není ovšem možné jednoznačně říci, jak omezit definiční obor funkce propustnosti t(x, y). Vzhledem k reciproké vlastnosti Fourierovy transformace je vhodné omezit definiční obor funkce propustnosti t(x, y) tak, že šířka stínítka tvoří asi 1/20 až 1/40 rozměru N T (N T je šířka omezeného definičního oboru). Jestliže bychom příliš omezili funkci t(x, y), tzn. šířka stínítka by byla srovnatelná s rozměrem N T , potom by středová část difrakčního obrazce byla příliš malá k dobrému rozlišení detailů. Vzdálenost dvou sousedních bodů dX ve výsledném obrazci je
2π . kN T
Jde o směrový
kosinus, tzn. kosinus úhlu měřeného od kolmice na směr šíření paprsku. Směrový kosinus je bezrozměrný. Pro praktické účely je však výhodnější znát úhlovou vzdálenost dvou sousedních bodů. Úhlovou vzdálenost dα dvou sousedních bodů spočítáme podle vzorce 2π π π dα = − arccos(dX) = − arccos , 2 2 kN T
(2)
kde jednotky dα jsou radiány (rad). Maximální hodnota úhlu, který difrakční obrazec zobrazuje N N · dα = αm = 2 2
2π π − arccos , 2 kN T
(3)
kde N je počet vzorků. Celý difrakční obrazec potom zobrazuje Fraunhoferovu difrakci v rozmezí h−αm ; αm i.
(4)
Pro běžné rozměry stínítka (tj. N T = (10 až 100) mm) a vlnovou délku λ = 632, 8 nm je úhlová vzdálenost dvou sousedních bodů dα = 0, 063 · 10−3 rad. Pro počet vzorků N = 1024 a rozměr N T = 100 mm by spočtený difrakční obrazec zobrazoval směrové rozložení intenzity difraktovaného světla v rozmezí h−3, 08 · 10−3 ; 3, 08 · 10−3 i rad.
36
9.3
Výpočet Fresnelovy difrakce pomocí rychlé Fourierovy transformace
Difrakční integrál pro výpočet Fresnelovy difrakce pomocí rychlé Fourierovy transformace je tvaru (viz 4(3)) k exp(ikz) ik 2 2 ψ(x, y, z) = −i exp (x + y ) × 2π z 2z Z∞Z ×
ik ik 2 2 (x + yM ) exp − (xxM + yyM ) dxM dyM . (1) ψM (xM , yM ) exp 2z M z
−∞
kde ψM (xM , yM ) je vlnová funkce v rovině stínítka. Fázor x2 + y 2 exp ik z + 2z
(2)
v integrálu (1) nemá význam pro výpočet intenzity difrakčního obrazce, a proto jej položíme roven jedné. Přepíšeme tedy integrál (1) do tvaru k ψ(x, y, z) = −i 2πz
Z∞Z
ik 2 ik 2 (x + yM ) exp − (xxM + yyM ) dxM dyM . ψM (xM , yM ) exp 2z M z
−∞
(3) V integrálu (3) lze závislost na z a λ postihnout jediným parametrem (viz 2.2(4)) n=
ρ20 kρ2 = 0, λz 2πz
(4)
udávajícím počet propuštěných difrakčních zón stínítkem, kde ρ0 je nejvzdálenější bod stínítka od středu stínítka M0 (viz obr. 2) (resp. poloměr kružnice opsané kolem stínítka se středem v bodě M0 ). S pomocí parametru (4) přepíšeme integrál (3) do tvaru kρ2 ψ x, y, 0 2πn
in =− 2 ρ0
Z∞Z −∞
= ψ1 (x, y, n) = x2 + y 2 ψM (xM , yM ) exp iπn M 2 M ρ0
37
xxM + yyM exp −i2πn dxM dyM . (5) ρ20
Tento integrál použijeme pro výpočet Fresnelovy difrakce pomocí rychlé Fourierovy transformace. Konstantu − ρin2 před integrálem položíme rovnu jedné, neboť budeme zobrazovat 0
pouze relativní intenzitu difrakčního obrazce. Vlastní výpočet provedeme tak, že spočteme rychlou Fourierovu transformaci funkce x2 +y 2 ψM (xM , yM ) exp iπn Mρ2 M a zobrazíme její kvadrát modulu. Opět omezíme definiční 0
obor této funkce tak, aby šířka stínítka tvořila asi 1/20 až 1/40 rozměru N T . Stěžejní část výpočtu je určení měřítka difrakčního obrazce. Měřítko určíme z porovnání argumentu fázoru rychlé Fourierovy transformace −i2π(l1 n1 + l2 n2 )/N,
(6)
a argumentu fázoru Fourierovy transformace (5) −i2πn
xxM + yyM , ρ20
(7)
pro osu x. V ose y je měřítko stejné, neboť požadujeme, aby byla stejná perioda vzorkování T v ose x a ose y a zároveň byl stejný počet vzorků N v jednotlivých osách. Abychom zjistili vzdálenost dx dvou sousedních bodů v difrakčním obrazci, dosadíme do fázoru rychlé Fourierovy transformace změnu indexu dl1 = 1 (resp. indexu dn1 = 1) a do fázoru Fourierovy transformace (5) dosadíme změnu dxM = T . Dosazením a úpravou dostaneme vztah −i2πl1 n1 /N = −i2πn −i2πdl1 dn1 /N = −i2πn 1/N = n
xxM ρ20 xdxM ρ20
(8)
dxT . ρ20
Pro vzdálenost dx dvou sousedních bodů difrakčního obrazce tedy platí dx =
ρ20 . nN T
(9)
V praxi je potřeba zobrazit v difrakčním obrazci úsečku, která udává určitou délku v (např. v = 5 mm). Pro počet bodů s, které je potřeba pro zobrazení úsečky v platí s=
v nN T = · v. dx ρ20
Tuto usečku vyneseme do výsledného difrakčního obrazce. 38
(10)
10
Obrazová část
V této kapitole jsou ukázány výsledky výpočtů některých difrakčních jevů. Všechny výpočty difrakčních jevů jsou pomocí rychlé Fourierovy transformace (viz kap. 9.2 a kap. 9.3). Dále u některých difrakčních jevů je ukázána shodnost teorie a experimentu. Fotografie difrakčních jevů byly pořízeny na optické lavici na Ústavu fyzikálního inženýrství VUT v Brně. Vzdálenosti na optické lavici byly měřeny pomocí pásma s přesností 1 mm, rozměry stínítek byly měřeny měřící lupou s přesností 0, 05 mm. Fotografický materiál byl černobílý film FOMAPAN Classic 100. Vyvolaný negativ byl naskenován do počítače pomocí skeneru HP ScanJet 5370 a dále digitálně zpracován. U výsledků výpočtů, které srovnáváme s experimentem, jsou barvy ve stupni šedi (hodnoty 0 ÷ 255) přetransformovány podle křivky, odpovídající svým průběhem přibližně křivce zčernání filmu (viz obr. 9).
Obrázek 9: Křivka zčernání. Obrázek převzat a upraven z programu Adobe Photoshop 7.0 CE.
10.1
Výpočty
V této části jsou uvedeny pouze příklady výpočtu difrakčních jevů.
39
(a)
(b)
(c)
(d)
(e)
(f)
(g)
(h)
(i)
Obrázek 10: Výpočet Fraunhoferovy a Fresnelovy difrakce na čtverci o straně a = 5 mm, vlnová délka difraktovaného světla λ = 632, 8 nm. a) zadání, b) Fraunhoferova difrakce, Fresnelova difrakce v rovinné vlně pro různý počet n propuštěných difrakčních zón (vodorovná úsečka vyznačuje 5 mm), c) n = 1, d) n = 2, e) n = 3, f) n = 4, g) n = 5, h) n = 6, i) n = 7.
40
(a) zadání
(b) Fraunhoferova difrakce
(c) n = 1
41
(d) n = 2
(e) n = 3
(f) n = 4
42
(g) n = 5
(h) n = 6
(i) n = 7 Obrázek 11: Výpočet Fraunhoferovy a Fresnelovy difrakce na kruhu (vlevo) o poloměru 25 mm, OK-filtru (kvadrantní fázová destička omezená kruhovým otvorem) o poloměru 25 mm osvětleném rovinnou vlnou (uprostřed) a OK-filtru o poloměru 25 mm osvětleném gaussovským svazkem (vpravo) (poměrný pokles intenzity na okraji stínítka je 1/e2 ). a) zadání, b) Fraunhoferova difrakce, Fresnelova difrakce pro různý počet n propuštěných difrakčních zón (vodorovná úsečka vyznačuje 50 mm), c) n = 1, d) n = 2, e) n = 3, f) n = 4, g) n = 5, h) n = 6, i) n = 7. OK-filtru se dá výužít pro vytyčování přímek (viz [8]).
43
(a)
(b)
Obrázek 12: Výpočet Fraunhoferovy difrakce na dvojčetné Siemensově hvězdici. Poloměr a = 2, 5 mm, vlnová délka difraktovaného světla λ = 632, 8 nm. a) zadání, b) Fraunhoferova difrakce.
(a)
(b)
Obrázek 13: Výpočet Fraunhoferovy difrakce na trojčetné Siemensově hvězdici. Poloměr a = 2, 5 mm, vlnová délka difraktovaného světla λ = 632, 8 nm. a) zadání, b) Fraunhoferova difrakce.
44
(a)
(b)
Obrázek 14: Výpočet Fraunhoferovy difrakce na tenkém mezikruží o poloměru a = 2, 5 mm, vlnová délka difraktovaného světla λ = 632, 8 nm. Profil řezu skrz střed difrakčního obrazce odpovídá kvadrátu Besselovy funkce J0 . a) zadání, b) Fraunhoferova difrakce.
(a)
(b)
Obrázek 15: Výpočet Fraunhoferovy difrakce na pěti kruhových otvorech rovnoměrně rozmístněných podél kružnice. Vnější poloměr a = 2, 5 mm, poloměr kruhového otvoru b = a/20, vlnová délka difraktovaného světla λ = 632, 8 nm. a) zadání, b) Fraunhoferova difrakce. 45
(a)
(b)
Obrázek 16: Výpočet Fraunhoferovy difrakce na sedmi kruhových otvorech rovnoměrně rozmístněných podél kružnice. Vnější poloměr a = 2, 5 mm, poloměr kruhového otvoru b = a/20, vlnová délka difraktovaného světla λ = 632, 8 nm. a) zadání, b) Fraunhoferova difrakce.
(a)
(b)
Obrázek 17: Výpočet Fraunhoferovy difrakce na osmi kruhových otvorech rovnoměrně rozmístněných podél kružnice. Vnější poloměr a = 2, 5 mm, poloměr kruhového otvoru b = a/20, vlnová délka difraktovaného světla λ = 632, 8 nm. a) zadání, b) Fraunhoferova difrakce. 46
(a)
(b)
Obrázek 18: a) Jean Baptiste Joseph Fourier, b) rychlá Fourierova transformace z portrétní perokresby.
47
10.2
Experimenty
10.2.1
Fraunhoferova difrakce
Abychom mohli ukázat shodu mezi výpočtem a experimentem u Fraunhoferovy difrakce, musíme nejprve určit měřítko (resp. zvětšení) pro dané experimentální uspořádání, neboť výsledky výpočtu Fraunhoferovy difrakce jsou v jednotkách radiánů, kdežto experimenty jsou v jednotká milimetrů. Schéma experimentálního uspořádání je na obrázku 19.
Obrázek 19: Schéma experimentálního uspořádání pro Fraunhoferovu difrakci, adaptováno z [5]. V kapitole 2.1 bylo řečeno, že Fraunoferovu difrakci bychom pozorovali v nekonečnu. Z praktického hlediska je to dosti složité, a proto využijeme spojné čočky, neboť čočka fokusuje rovnoběžné svazky do své ohniskové roviny. Obrazová ohnisková rovina F 0 čočky je obrazem roviny v nekonečnu a body této roviny odpovídají směrům předmětového prostoru čočky (viz [4], kap. 3.2 nebo [9], kap 8.3.1). Podle experimentálního uspořádání (viz obr. 19) je tedy Fraunhoferova difrakce v ohniskové rovině F 0 čočky C2 (obrazová ohnisková vzdálenost f 0 čočky C2 je f 0 = 1200 mm). Na stínítko S dopadá rovinná vlna. Fraunhoferovu difrakci pozorujeme objektivem Obj. (parametry objektivu jsou z hlediska stanovení měřítka nepodstatné) v rovině pozorování P ve vzdálenosti d = 480 mm od obrazové ohniskové roviny F 0 čočky C2 , neboť Fraunhoferova difrakce v ohnisku čočky je příliš malá. Pro stanovení měřítka jsme vyfotografovali Fraunhoferovu difrakci na osmi kruhových otvorech ležících na jedné přímce. Poloměr každého kruhového otvoru je r = 0, 25 mm a vzdálenost středů jednotlivých otvorů je a = 1 mm. Pro dané stínítko je velice jednodu48
ché vypočítat analytický tvar Fraunhoferovy difrakce (To je učiněno např. v [5], kap. 11). Z výsledku výpočtu vyplývá, že vzdálenost mezi sousedními hlavními maximy je
π . ka
Pro
úhlovou vzdálenost α mezi sousedními hlavními maximy potom platí α= kde k =
2π , λ
π 2
− arccos
π rad, ka
(11)
λ je vlnová délka difraktovaného světla.
Vzdálenost y dvou sousedních maxim vyfotografovaného difrakčního obrazce změříme. Z úhlové vzdálenosti α a vzdálenosti y mezi sousedními hlavními maximy určíme měřítko tak, že 1 · 10−3 rad u spočteného difrakčního obrazce odpovídá
y α
mm ([y] = mm) u vyfoto-
grafovaného difrakčního obrazce, a naopak 1 mm u vyfotografovaného difrakčního obrazce odpovídá
α y
· 10−3 rad ([α] = 1 · 10−3 rad) u spočteného difrakčního obrazce.
V našem případě platí α = 0, 316 · 10−3 rad a y = 3, 13 mm. Platí tedy, že 1 · 10−3 rad u spočteného difrakčního obrazce odpovídá 9, 905 mm u vyfotografovaného difrakčního obrazce, a naopak 1 mm u vyfotografovaného difrakčního obrazce odpovídá 0, 101·10−3 rad u spočteného difrakčního obrazce. Pro výpočet Frauhoferovy difrakce musíme zvolit vhodnou periodu vzorkování T a stanovit rozměr N T (viz kap. 9.2). Vypočítaný difrakční obrazec potom zobrazuje úhlové rozložení intenzity difraktovaného světla v rozmezí h−αm ; αm i rad, kde αm = 2π N π − arccos kN (viz 9.2 (3)). Teoretická šířka (resp. výška) vyfotografovaného difrakč2 2 T ního obrazce by potom měla velikost 2m αy . Na obrázku 20 je fotografie Fraunhoferovy difrakce pro výše zmiňované stínítko a výpočet Fraunhoferovy difrakce na tomto stínítku. Nad měřítkem je vždy uveden rozměr, který měřítko vyjadřuje v rovině pozorování (platí i pro Fresnelovu difrakci).
49
(a)
(b)
(c) Obrázek 20: Fraunhoferova difrakce na osmi kruhových otvorech ležících na jedné přímce. Poloměr každého kruhového otvoru je r = 0, 25 mm a vzdálenost středů jednotlivých otvorů je a = 1 mm, vlnová délka difraktovaného světla λ = 632, 8 nm. a) zadání, b) spočtená Fraunhoferova difrakce, c) fotografie Fraunhoferovy difrakce.
50
(a)
(b)
(c) Obrázek 21: Fraunhoferova difrakce na osmi kruhových otvorech ležících rovnoměrně podél kružnice o poloměru r0 = 5 mm. Poloměr každého kruhového otvoru je r = 0, 25 mm, vlnová délka difraktovaného světla λ = 632, 8 nm. a) zadání, b) spočtená Fraunhoferova difrakce, c) fotografie Fraunhoferovy difrakce.
51
(a)
(b)
(c) Obrázek 22: Fraunhoferova difrakce na kosodélníkovém otvoru. Delší strany mají délku 3, 4 mm, kratší strany mají délku 2, 3 mm a vzdálenost dvou delších stran je 2 mm, vlnová délka difraktovaného světla λ = 632, 8 nm. a) zadání, b) spočtená Fraunhoferova difrakce, c) fotografie Fraunhoferovy difrakce.
52
(a)
(b)
(c) Obrázek 23: Fraunhoferova difrakce na elipsovém otvoru. Hlavní poloosa a = 0, 9 mm a vedlejší poloosa b = 0, 65 mm, vlnová délka difraktovaného světla λ = 632, 8 nm. a) zadání, b) spočtená Fraunhoferova difrakce, c) fotografie Fraunhoferovy difrakce.
53
10.2.2
Fresnelova difrakce
Experimentální uspořádání pro Fresnelovu difrakci bylo ukázáno již v obrázku 1 (viz kap. 2.1). Stanovení měřítka pro Fresnelovu difrakci bylo ukázáno v kapitole 9.3.
(a)
(b)
Obrázek 24: Fresnelova difrakce na kruhovém otvoru o poloměru r = 1, 15 mm. Vzdálenost stínítka od bodového zdroje z1 = 70, 1 cm, vzdálenost roviny pozorování od stínítka z = 205, 4 cm, vlnová délka difraktovaného světla λ = 632, 8 nm. a) fotografie Fresnelovy difrakce, b) výpočet Fresnelovy difrakce.
54
(a)
(b) Obrázek 25: Fresnelova difrakce na OK-filtru o poloměru r = 1, 95 mm, vlnová délka difraktovaného světla λ = 632, 8 nm. Vzdálenost stínítka od bodového zdroje z1 = 86, 0 cm, vzdálenost roviny pozorování od stínítka z = 294, 8 cm. a) fotografie Fresnelovy difrakce, b) výpočet Fresnelovy difrakce.
55
(a)
(b)
(c)
(d)
Obrázek 26: Fresnelova difrakce na dvou kruhových otvorech o poloměru r = 1 mm. Vzdálenost středů kruhových otvorů a = 2, 6 mm, vlnová délka difraktovaného světla λ = 632, 8 nm. Vzdálenost stínítka od bodového zdroje z1 = 70, 1 cm. a) fotografie Fresnelovy difrakce, b) výpočet Fresnelovy difrakce pro vzdálenost roviny pozorování od stínítka z = 79, 1 cm. c) fotografie Fresnelovy difrakce, d) výpočet Fresnelovy difrakce pro vzdálenost roviny pozorování od stínítka z = 162, 4 cm.
56
(a)
(b)
(c)
(d)
Obrázek 27: Fresnelova difrakce na dvou kruhových otvorech o poloměru r = 1 mm. Vzdálenost středů kruhových otvorů a = 2, 03 mm, vlnová délka difraktovaného světla λ = 632, 8 nm. Vzdálenost stínítka od bodového zdroje z1 = 70, 1 cm. a) fotografie Fresnelovy difrakce, b) výpočet Fresnelovy difrakce pro vzdálenost roviny pozorování od stínítka z = 79, 1 cm. c) fotografie Fresnelovy difrakce, d) výpočet Fresnelovy difrakce pro vzdálenost roviny pozorování od stínítka z = 162, 4 cm.
57
(a)
(b)
Obrázek 28: Fresnelova difrakce na pravoúhlém rohu. Zajímavé na této difrakci je, že ji pozoroval již Grimaldi v polovině 17. století (viz [10] (podle [4])). Pozoroval tzv. laloky, které se objevují v geometrickém stínu za stínítkem. Výpočet Fresnelovy difrakce na tomto stínítku je poněkud obtížný, neboť hranice pravoúhlého rohu sahají do nekonečna, proto lze spočíst difrakci pouze na velkém obdélníku. Z tohoto důvodu také není uvedeno měřítko u vypočteného obrázku. Výpočtem difrakce na tomto stínítku se zabýval např. Komrska (viz [11], str. 191-208). Vlnová délka difraktovaného světla λ = 632, 8 nm. Vzdálenost stínítka od bodového zdroje z1 = 115, 4 cm, vzdálenost roviny pozorování od stínítka z = 65, 4 cm. a) fotografie Fresnelovy difrakce, b) výpočet Fresnelovy difrakce.
58
11
Závěr
V diplomové práci je zvolen přístup k definici diskrétní Fourierovy transformace, kdy definice diskrétní Fourierovy transformace vychází z numerického výpočtu spojité Fourierovy transformace. Běžné definice (resp. odvození) diskrétní Fourierovy transformace (viz např. [2], kap. 6) mohou vést čtenáře k mylné domněnce, že existují dvě Fourierovy transformace. Je tedy vhodné připomenout, že definice Fourierovy transformace je jen jedna a že diskrétní Fourierova transformace je pouze nástrojem pro výpočet spojité Fourierovy transformace. Diskrétní Fourierova transformace má sloužit pro výpočet Fourierovy transformace funkcí, který nelze provést analyticky a dále pro výpočty Fourierovy transformace pomocí výpočetní techniky. Výpočet difrakčních jevů v optice pomocí rychlé Fourierovy transformace je pravděpodobně nejefektivnější způsob výpočtu. Přináší ovšem některá úskalí, zvláště při výpočtech difrakčních jevů na stínítcích sahajících svými rozměry do nekonečna. Navíc vlastní zadání tvaru stínítka není vždy přesné, neboť např. zadání kruhového otvoru do pravoúhlé sítě bodů činí problémy. Tuto skutečnost dobře ilustruje obrázek 23b (viz kap. 10.2), kdy hlavní minima nemají kruhový tvar. I přes to je shodnost mezi výpočtem a experimentem velice dobrá (viz kap. 10.2) a rychlé Fourierovy transformace lze tedy s výhodou využívat pro výpočty difrakčních jevů. Na závěr je vhodné upozornit na to, že rychlé Fourierovy transformace, jako i jiných prostředků moderní výpočetní techniky, je nutné užívat s rozvahou a pochopením, neboť praxe mnohých ukazuje, že neučiníme-li tak, dojdeme většinou k špatnému výsledku a zbytečné ztrátě času.
59
Reference [1] Cooley, J., Tuckey J.: An Algorithm for the Machine Calculation of Complex Fourier Series. Mathematics of Computation 19, 1965, 297-301. [2] Brigham, E. O.: The Fast Fourier Transform and its Applications. 2nd ed., PrenticeHall, Inc. Engelwood Clifs, New Jersey 1988. [3] Press H. W. et al.: Numerical recipes in Fortran 77: The art of scientific computing. Cambridge University Press, Cambridge 1992. [4] Komrska, J.: Difrakce světla. VUTIUM, Brno 2001. [5] Komrska, J.: Fourierovské metody v teorii difrakce a ve strukturní analýze. VUTIUM, Brno 2001. [6] Čížek, V.: Diskrétní Fourierova transformace a její použití. Nakladatelství technické literatury, Praha 1981. [7] Bracewell, R. M.: The Fourier Transform and its Applications. McGraw-Hill, Inc. New York 1965. [8] Velechovský, K.: Vytyčování přímek využitím Fresnelovy difrakce. Disertační práce, FSI Vysoké učení technické v Brně, Brno 2003. [9] Goodman, J. W.: Introduction to Fourier Optics. 2nd ed., McGraw-Hill, New York 1996. [10] Grimaldi F. M.: Physico-mathesis de lumine, coloribus, et iride. Bonaniae 1665. [11] Komrska, J.: Scalar Difraction Theory in Electron Optics. Advances in Electronics and Electron Physics. (L. Marton, ed.), Vol. 30. Academic Press, Inc., New York and London 1971.
60
