Fourierova transformace: rozklad signálu na frekvence

Nechť T > 0 a x : ℝ → ℝ je periodická funkce s periodou T, tj. x(t + T) = x(t) pro všechna t ∈ ℝ. Číslo 1/T je její základní frekvence. Nechť a₀, a₁, a₂, … a b₁, b₂, … ∈ ℝ jsou Fourierovy koeficienty funkce x (další definice). Fourierova řada funkce x je řada

a₀ + ∑_k=1^∞ [a_k cos(2πkt/T) + b_k sin(2πkt/T)];

její K-tý částečný součet S_K(t) obsahuje a₀ a členy k = 1, …, K. Sinusovky s frekvencí k/T, k ≥ 1, se nazývají harmonické (k-tá harmonická).

„Mírné předpoklady“, za kterých řada funkci skutečně vyjádří (Dirichletova věta): je-li x na periodě po částech spojitě diferencovatelná, pak řada v každém bodě t konverguje k x(t−) + x(t+)/2, kde x(t−) a x(t+) jsou jednostranné limity. V bodech spojitosti tedy k x(t), ve skoku k průměru obou stran.

Nechť T > 0 a x : ℝ → ℝ je periodická s periodou T a integrovatelná na [0, T]. Její Fourierovy koeficienty jsou čísla

a₀ = 1/T ∫₀^T x(t) dt,

a_k = 2/T ∫₀^T x(t) cos(2πkt/T) dt,   b_k = 2/T ∫₀^T x(t) sin(2πkt/T) dt   pro k = 1, 2, …

Funguje to díky ortogonalitě: pro celá k, l ≥ 1 je ∫₀^T sin(2πkt/T) cos(2πlt/T) dt = 0 a pro k ≠ l je nulový i integrál součinu dvou sinů a dvou kosinů, kdežto ∫₀^T sin²(2πkt/T) dt = T/2. Vynásobení signálu k-tou sinusovkou a integrace proto z celé řady „vyloví“ jediný člen, b_k · T/2; činitel 2/T to T/2 vyrovná.

Obdélník s periodou T = 1: x(t) = 1 pro 0 < t < ½ a −1 pro ½ < t < 1.

a₀ = 0, protože kladná a záporná půlka mají stejnou plochu. a_k = 0, protože obdélník je lichý (x(−t) = −x(t)) a kosinus sudý — součin je lichý a jeho integrál přes periodu nulový.

∫₀^½ sin(2πkt) dt = 1 − cos πk/2πk a druhá půlka dá totéž (záporný signál krát záporný sinus), takže b_k = 2 · 2 · 1 − cos πk/2πk = 2(1 − cos πk)/πk.

Pro liché k je cos πk = −1, tedy b_k = 4/πk: b₁ = 1,273, b₃ = 0,424, b₅ = 0,255. Pro sudé k je cos πk = 1 a b_k = 0.

Kontrola uprostřed kladné půlky, t = 1/4, pro K = 3 (barvy podle vzorce níže): 4/π · (sin(π/2)/1 + sin(3π/2)/3 + sin(5π/2)/5) = 1,273 · (1 − 0,333 + 0,2) = 1,103; jeden člen 1,273. Součet se k hodnotě 1 blíží zleva i zprava, jak se přidávají členy s opačnými znaménky.

Nechť x je obdélník s periodou 1 z příkladu výše a pro K ≥ 1 nechť

S_K(t) = 4/π ∑_k=1^K sin(2π(2k − 1)t)/2k − 1

je jeho částečný součet. Pak

lim_{K → ∞} max_{t ∈ (0, 1/2)} S_K(t) = 2/π ∫₀^π sin u/u du ≈ 1,1790.

Překmit nad hodnotu 1 je tedy 0,179, a protože skok z −1 na 1 má výšku 2, je to 0,179/2 ≈ 9 % výšky skoku — pro velká K nezávisle na K. Jev popsal Henry Wilbraham (1848) a znovu J. Willard Gibbs (1899).

Náčrt důkazu. První maximum za skokem leží v t = 1/4K. Označme u_k = π(2k − 1)/2K; body jsou od sebe vzdálené Δu = π/K a vyplní interval (0, π). Pak sin(2π(2k − 1)t)/2k − 1 = sin u_k/u_k · Δu/2 a částečný součet je 4/π · 1/2 · ∑ sin u_k/u_k Δu — Riemannův součet integrálu výše. Náš výpočet: 1,273 pro jeden člen, 1,180 pro deset, 1,179 pro dvě stě.

Nechť N ≥ 1 a x[0], …, x[N − 1] ∈ ℂ (typicky reálné vzorky signálu vzorkovaného frekvencí f_s > 0). Nechť i je imaginární jednotka, i² = −1, a exp(iθ) = e^iθ komplexní exponenciála (další definice). Diskrétní Fourierova transformace (DFT) vzorků je posloupnost X[0], …, X[N − 1] ∈ ℂ,

X[k] = ∑_n=0^N−1 x[n] · exp(−2πikn/N),   k = 0, 1, …, N − 1.

Index k se nazývá koš a odpovídá mu frekvence f_k = k · f_s/N (v Hz). Amplitudové spektrum je posloupnost absolutních hodnot |X[k]| vynesená proti f_k; arg X[k] je fáze příslušné složky. Hodnota |X[k]| tedy říká, jak silně je frekvence f_k v signálu zastoupená.

Nechť i je imaginární jednotka, i² = −1, a ℂ = {a + ib : a, b ∈ ℝ} jsou komplexní čísla. Pro θ ∈ ℝ je komplexní exponenciála dána Eulerovým vzorcem

e^iθ = exp(iθ) = cos θ + i sin θ

(shoduje se s dosazením iθ do Taylorovy řady exponenciály). Je to bod jednotkové kružnice v komplexní rovině, otočený od 1 o úhel θ proti směru hodinových ručiček; e^−iθ = cos θ − i sin θ je otočení na druhou stranu. Nechť z = a + ib, a, b ∈ ℝ.

1. Absolutní hodnota (vzdálenost od nuly): |z| = √a² + b².
2. Komplexně sdružené číslo: z^* = a − ib.
3. Pro z ≠ 0 je argument (fáze) arg z to jediné číslo z intervalu (−π, π], pro které z = |z| e^{i arg z} — úhel, který z svírá s kladnou reálnou osou.

Signál x = (1, 2, 1, 0) je stejnosměrná složka 1 plus sinusovka s amplitudou 1, která se za čtyři vzorky otočí jednou: x[n] = 1 + sin(πn/2).

Pro N = 4 je exp(−2πikn/4) = (−i)^kn, takže stačí mocniny −i: 1, −i, −1, i (barvy podle vzorce výše).

Čtení: stejnosměrná složka je X[0]/N = 4/4 = 1, amplituda sinusovky 2 · |X[1]|/N = 2 · 2/4 = 1. Fáze arg(−2i) = −90° říká „sinus, ne kosinus“. Koš k = 2 je prázdný a X[3] je komplexně sdružené k X[1] — viz další větu.

Nechť N ≥ 1, x[0], …, x[N − 1] ∈ ℝ a X je jejich DFT. Pak pro každé k = 1, …, N − 1 platí

X[N − k] = X[k]^*

(komplexně sdružené číslo). Důkaz: pro celé n je exp(−2πi(N − k)n/N) = exp(−2πin) · exp(2πikn/N) = 1 · exp(−2πikn/N)^*, sdružení součtu je součet sdružení a reálné x[n] se sdružením nemění. ∎

Důsledky: nová informace je jen v koších 0 až N/2 (pro 8 192 vzorků jich je 4 097) a energie sinusovky se dělí napůl mezi k a N − k. Proto se amplituda počítá jako 2 · |X[k]|/N — dvojka vrací polovinu ze zrcadlového koše, dělení N odstraní to, že se sčítalo N vzorků.

Vstup: amplitudové spektrum S[0], …, S[M] ≥ 0 (koše 0 až M = N/2 z DFT délky N), vzorkovací frekvence f_s > 0 a práh p ≥ 0. Výstup: seznam upřesněných frekvencí vrcholů v Hz.

1. Projdi koše k = 1, …, M − 1 a vyber ty, pro které S[k] > p, S[k] > S[k − 1] a S[k] > S[k + 1] (ostrá lokální maxima nad prahem).
2. Pro každý vybraný koš polož a = S[k − 1], b = S[k], c = S[k + 1].
3. Spočti posun vrcholu paraboly proložené body (−1, a), (0, b), (1, c): d = 1/2 · a − c/a − 2b + c. Protože b je ostré maximum, je jmenovatel záporný a −1/2 < d < 1/2.
4. Přidej do výstupu frekvenci f = (k + d) · f_s/N.

Algoritmus skončí, když projde všechny koše. Proč je potřeba: koše jsou od sebe 8 000/8 192 = 0,977 Hz a 440 Hz leží v koši 440 · 8 192/8 000 = 450,56 — mezi koši. Vrchol se proto upřesní parabolou proloženou třemi sousedními koši:

pro každý koš k, který je vyšší než práh a než oba sousedé:
    a, b, c ← S[k − 1], S[k], S[k + 1]
    d ← ½ · (a − c) / (a − 2b + c)      # posun vrcholu paraboly, −½ … ½ koše
    f ← (k + d) · fs / N                # upřesněná frekvence v Hz

U tónu a: nejvyšší je koš 451 (440,43 Hz) s hodnotou b = 0,8790, sousedé a = 0,8058 (koš 450) a c = 0,2023 (koš 452). d = 1/2 · 0,6035/0,8058 − 1,7580 + 0,2023 = 0,3018/−0,7499 = −0,402, takže f = 450,598 · 8 000/8 192 = 440,04 Hz. Skutečná frekvence 440 Hz, chyba 0,04 Hz — dvacetkrát méně než šířka koše.

Nechť f_s > 0 je vzorkovací frekvence a N ≥ 1 počet vzorků. Délka záznamu je T = N/f_s (v sekundách). Frekvenční rozlišení DFT je vzdálenost frekvencí sousedních košů

Δf = f_k+1 − f_k = f_s/N = 1/T.

Dva tóny, jejichž frekvence se liší o méně než zhruba Δf (s oknem o něco víc), proto ve spektru splynou v jeden vrchol.

Rozdíl je 10 Hz, takže potřebujeme Δf < 10 Hz, tedy T > 1/10 s a N = f_s · T > 8 000 · 0,1 = 800 vzorků (barvy podle vzorce výše).

Nechť N ≥ 2, f_s > 0 a x[n] = A sin(2πf · n/f_s + φ), n = 0, …, N − 1, jsou vzorky sinusovky s frekvencí 0 < f < f_s/2. Frekvence f leží na koši, právě když f = k · f_s/N pro nějaké celé k; pak je X[k′] = 0 pro všechny koše k′ ∉ {k, N − k}. Prosakování (spectral leakage) je jev, kdy f na koši neleží: pak je X[k′] ≠ 0 i v koších daleko od f — energie sinusovky se rozlije do celého spektra.

Okenní funkce je posloupnost w[0], …, w[N − 1] ∈ [0, 1]; spektrum s oknem je DFT posloupnosti w[n] · x[n]. Okno se volí tak, aby se signál na okrajích úseku plynule ztlumil k nule. Hannovo okno je

w[n] = 1/2 − 1/2 cos(2πn/N) = sin²(πn/N),   n = 0, …, N − 1.

Nechť N ≥ 2, x[0], …, x[N − 1] ∈ ℂ, X je jejich DFT, w Hannovo okno délky N a X_w DFT posloupnosti w[n] · x[n]. Indexy X počítejme dokola modulo N (X[−1] = X[N − 1], X[N] = X[0]). Pak pro každé k = 0, …, N − 1

X_w[k] = 1/2 X[k] − 1/4 X[k − 1] − 1/4 X[k + 1].

Důkaz: podle Eulerova vzorce je cos θ = e^iθ + e^−iθ/2, takže w[n] = 1/2 − 1/4 exp(2πin/N) − 1/4 exp(−2πin/N). Vynásobení signálu exp(2πin/N) posune celé spektrum o jeden koš: ∑_n x[n] exp(2πin/N) exp(−2πikn/N) = ∑_n x[n] exp(−2πi(k − 1)n/N) = X[k − 1]; obdobně exp(−2πin/N) dá X[k + 1]. Zbytek je linearita DFT. ∎

Okno tedy ve frekvencích jen zprůměruje každý koš se sousedy vahami (−¼, ½, −¼). Daleko od tónu se prosakování mění z koše na koš jen pomalu (nenormovaná DFT v koších 149, 150, 151 je −2,815; −2,750; −2,687, imaginární části nulové), takže ½ X[k] a ¼ (X[k − 1] + X[k + 1]) jsou skoro stejná čísla a téměř se odečtou: ½ · (−2,750) − ¼ · (−2,815 − 2,687) = −1,375 + 1,3755 = 0,0005. U tónu, kde se koše mění prudce, k odečtení nedojde — vrchol jen rozšíří na tři koše. Na vzorku 785,2 Hz to platí do zaokrouhlovací chyby (rozdíl pod 10⁻¹²). Je to tatáž konvoluce, o které je díl Filtry a konvoluce.

Sinusovka vzdálená δ košů od koše k v něm bez okna zanechá zhruba podíl |sin(πδ)/πδ| své amplitudy. Náš tón leží přesně mezi koši, δ = 1/2: sin(π/2)/π/2 = 2/π ≈ 0,637, tedy −3,9 dB. S Hannovým oknem je tatáž ztráta jen −1,4 dB, tedy 0,85. Obě hodnoty uvádí jako „scalloping loss“ přehled okenních funkcí F. J. Harrise (1978).

Ve vzdálenosti δ = 49,5 koše (koš 150, asi 390 Hz od tónu) dá odhad 1/πδ = 1/155,5 = 0,0064 amplitudy, vůči vrcholu 0,637 tedy 0,0101, neboli −39,9 dB. Přesný výpočet dá −41,5 dB; rozdíl dělá zrcadlový obraz tónu na záporné frekvenci, který se do koše 150 přičte s opačnou fází. Úroveň bez okna klesá jen jako 1/δ (6 dB na oktávu), s Hannovým oknem jako 1/δ³ (18 dB na oktávu, opět Harris 1978) — odtud −110 dB.

Nechť N ≥ 1 a x[0], …, x[N − 1] ∈ ℝ jsou vzorky signálu. Výkon signálu je průměr druhých mocnin

P = 1/N ∑_n=0^N−1 x[n]².

Sinusovka s amplitudou A (přes celé periody) má výkon A²/2, protože průměr sin² přes periodu je 1/2; šum s nulovou střední hodnotou a směrodatnou odchylkou σ má (očekávaný) výkon σ². Nechť P_signál, P_šum > 0 jsou výkony užitečného signálu a šumu. Poměr signálu k šumu a jeho hodnota v decibelech jsou

SNR = P_signál/P_šum,   SNR_dB = 10 · log₁₀ SNR.

Pro poměr amplitud A₁/A₂ se píše 20 · log₁₀ A₁/A₂, protože výkon roste s druhou mocninou amplitudy.

Tady: P_tón = 0,1²/2 = 0,005, P_šum = 1² = 1, SNR = 0,005 = 1/200, tedy 10 · log₁₀ 0,005 = −23,0 dB.

Zisk pro N = 8 192: 10 · log₁₀ 8 192/2 = 10 · log₁₀ 4 096 = 36,1 dB, takže SNR_koš = −23,0 + 36,1 = 13,1 dB (barvy podle vzorce výše), poměr výkonů 20,5 a poměr amplitud √20,5 ≈ 4,5.

Šum v jednom koši: komplexní koeficient X[k] má střední kvadrát Nσ², po přepočtu na amplitudu 2/N · |X[k]| je efektivní hodnota 2σ/√N = 2/90,5 = 0,0221. Absolutní hodnota takového šumu má Rayleighovo rozdělení s mediánem 0,0221 · √ln 2 = 0,0184 — to je ten „typický koš šumu“ 0,019 z odstavce výše. Tón s amplitudou 0,1 je nad ním 0,1/0,0221 = 4,5krát, přesně podle zisku výše.

Nejvyšší ze zhruba 4 090 košů šumu bude kolem 0,0221 · √ln 4 090 ≈ 0,064; změřili jsme 0,061. Proto je dobré, aby tón vyčníval aspoň několikanásobně: jen nad typickým košem nestačí, musí převýšit i ten nejvyšší.

FFT (rychlá Fourierova transformace) je každý algoritmus, který spočítá DFT délky N za řádově N log₂ N operací místo N² u výpočtu přímo podle definice. Nejznámější je varianta Cooleyho a Tukeyho (1965):

Vstup: N = 2^m vzorků x[0], …, x[N − 1] ∈ ℂ, m ≥ 0. Výstup: jejich DFT X[0], …, X[N − 1].

1. Je-li N = 1, vrať X[0] = x[0] (podmínka zastavení rekurze).
2. Jinak rozděl vzorky na sudé x[0], x[2], …, x[N − 2] a liché x[1], x[3], …, x[N − 1].
3. Rekurzivně spočti E = FFT(sudé) a O = FFT(liché), obě délky N/2.
4. Pro k = 0, …, N/2 − 1 polož t = W^k O[k], kde W = exp(−2πi/N), a X[k] = E[k] + t, X[k + N/2] = E[k] − t.
5. Vrať X.

Sudé vzorky (1, 1): E = (1 + 1, 1 − 1) = (2, 0). Liché vzorky (2, 0): O = (2 + 0, 2 − 0) = (2, 2).

Otočné činitele pro N = 4: W⁰ = 1, W¹ = exp(−2πi/4) = exp(−πi/2) = −i.

Totéž jako přímý výpočet v oddílu o DFT, jen se čtyřmi násobeními (z nichž jediné netriviální, krát −i) místo šestnácti.

Nechť N = 2^m, m ≥ 0, a M(N) je počet komplexních násobení, které provede FFT výše (včetně násobení jedničkou). Pak

M(N) = N/2 · log₂ N.

Důkaz. Z algoritmu M(N) = 2 M(N/2) + N/2 pro N ≥ 2 a M(1) = 0. Na úrovni rekurze j = 0, …, m − 1 je 2^j podúloh, každá po N/2^j+1 násobeních, tedy dohromady N/2 násobení na úroveň; úrovní je m = log₂ N. ∎

Pro N = 8 192 = 2¹³: 4 096 · 13 = 53 248, proti N² = 67 108 864 u přímé DFT, poměr 1 260. Násobení jedničkou se tu pro jednoduchost počítají také.

Nechť d ≥ 2 je sudé (rozměr vektorů modelu) a pos ∈ {0, 1, 2, …} pozice slova ve vstupu. Sinusové kódování pozice podle článku Vaswani a kol., Attention Is All You Need (2017), je vektor PE(pos) ∈ ℝ^d, jehož složky s indexy 0, …, d − 1 jsou pro j = 0, …, d/2 − 1

PE(pos, 2j) = sin(pos/10 000^2j/d),   PE(pos, 2j + 1) = cos(pos/10 000^2j/d).

Vektor PE(pos) se přičte k vektoru slova na pozici pos.

Každá dvojice složek je jedna „harmonická“: vlnové délky rostou geometricky od 2π (pro j = 0) po 2π · 10 000^(d−2)/d, při d = 512 asi 2π · 9 647 ≈ 60 600 pozic. Rychlé složky rozliší sousední slova, pomalé vzdálená — podobně jako ručičky hodin. Dvojice sin a cos se volí proto, že posun o pevný počet pozic je pak pro každou dvojici jen otočení o pevný úhel, tedy lineární operace, kterou se síť snadno naučí.