Filtry a konvoluce: vyhlazení, hrany a odstranění šumu

Nechť x = (x[n])_{n ∈ ℤ} je posloupnost reálných čísel (signál) a M ∈ ℕ je liché, M = 2m + 1 pro nějaké m ∈ {0, 1, 2, …}. Centrovaný klouzavý průměr délky M signálu x je posloupnost y = (y[n])_{n ∈ ℤ} daná pro každé n ∈ ℤ předpisem

y[n] = 1/M ∑_j=−m^m x[n + j].

Hodnota y[n] je tedy průměr vzorku x[n] a m sousedů na každé straně. Centrování zajistí, že výstup není proti vstupu zpožděný. Má-li signál jen konečně mnoho vzorků x[0], …, x[L − 1], chybějí na okrajích sousedé; tady je nahrazujeme krajní hodnotou (x[n] = x[0] pro n < 0 a x[n] = x[L − 1] pro n ≥ L) a čísla v tabulce počítáme jen ve vnitřní části (bez 40 krajních vzorků z každé strany).

Nechť M ∈ ℕ, M ≥ 1, a e₁, …, e_M jsou nezávislé náhodné veličiny, které mají všechny tentýž rozptyl σ² (σ ≥ 0). Pak

var(1/M ∑_j=1^M e_j) = σ²/M,   a tedy směrodatná odchylka průměru je σ/√M.

Důkaz. Rozptyl součtu nezávislých veličin je součet rozptylů a konstanta vstupuje do rozptylu na druhou, var(cX) = c² var(X). Proto var(1/M ∑ e_j) = 1/M² · M σ² = σ²/M. ∎

Pro σ = 0,4: M = 5 dá 0,4/√5 = 0,179, M = 21 dá 0,087, M = 61 dá 0,051. V tabulce je o něco víc (0,195; 0,099; 0,065), protože je to jediná konkrétní realizace šumu; průměr přes 300 různých realizací vychází 0,178; 0,087 a 0,051.

Barvy podle vzorce výše; chyba = √(šum)² + 1/2 · (zkreslení)².

M = 21: šum 0,087, zkreslení 1 − 0,929 = 0,071:

chyba ≈ √0,087² + 1/2 · 0,071² = √0,0076 + 0,0025 = 0,101.

M = 61: šum 0,051, zkreslení 1 − 0,491 = 0,509:

chyba ≈ √0,051² + 1/2 · 0,509² = √0,0026 + 0,130 = √0,132 = 0,364 — vládne zkreslení.

M = 5: šum 0,179, zkreslení jen 0,004, chyba ≈ 0,179 — vládne šum.

Projdeme-li všechna lichá M, vzorec má minimum při M = 19 a 21 (obě 0,101). Na naší jediné realizaci šumu vyšlo 0,124 — o něco víc, ale se stejným vítězem.

Nechť x = (x[n])_{n ∈ ℤ} je posloupnost reálných čísel (signál) a h = (h[k])_{k ∈ ℤ} posloupnost reálných čísel, která má jen konečně mnoho nenulových členů (jádro, též filtr). Signál konečné délky chápeme jako posloupnost doplněnou nulami. Konvoluce signálu x s jádrem h je posloupnost x * h daná pro každé n ∈ ℤ předpisem

(x * h)[n] = ∑_{k ∈ ℤ} h[k] · x[n − k]

(součet má jen konečně mnoho nenulových sčítanců). Pro takto definovanou konvoluci platí:

1. linearita: (αx + βz) * h = α(x * h) + β(z * h) pro všechny signály x, z a všechna α, β ∈ ℝ;
2. nezávislost na posunu: označíme-li (S_dx)[n] = x[n − d] posun signálu o d ∈ ℤ vzorků, pak (S_dx) * h = S_d(x * h) — posunutý vstup dá stejně posunutý výstup;
3. impulzní odezva: pro jednotkový impulz δ (δ[0] = 1, δ[n] = 0 pro n ≠ 0) je δ * h = h.

Jádro je tedy zároveň impulzní odezva filtru — to, co filtr udělá z jediného jedničkového vzorku.

(1, 2, 3) * (1, 1), indexy od nuly; výsledek má délku 3 + 2 − 1 = 4:

x = (1, 2, 3), h = (1, 1), y = x * h; barvy podle vzorce výše:

y[0] = h[0] · x[0] = 1 · 1 = 1,   y[1] = h[0] · x[1] + h[1] · x[0] = 1 · 2 + 1 · 1 = 3,

y[2] = h[0] · x[2] + h[1] · x[1] = 1 · 3 + 1 · 2 = 5,   y[3] = h[1] · x[2] = 1 · 3 = 3.

Na nesymetrickém jádru je vidět, proč je tam minus: (1, 2, 3) * (1, −1) = (1, 2 − 1, 3 − 2, −3) = (1, 1, 1, −3) — každé číslo je „současný mínus předchozí“, tedy diference. Bez otočení (tzv. korelace, ∑ h[k] x[n + k]) by vyšlo totéž s opačným znaménkem. Konvoluční vrstvy neuronových sítí ve skutečnosti počítají korelaci; protože se jádro učí, na tom nezáleží.

Nechť x, z : ℤ → ℝ jsou signály, g, h : ℤ → ℝ jádra (posloupnosti s konečně mnoha nenulovými členy) a α, β ∈ ℝ. Pak platí:

1. komutativita: x * h = h * x, kde (h * x)[n] = ∑_k x[k] · h[n − k] (důkaz substitucí k → n − k);
2. asociativita: (x * g) * h = x * (g * h);
3. linearita: (αx + βz) * h = α(x * h) + β(z * h).

Z asociativity plyne praktický důsledek: dva filtry za sebou (nejdřív g, pak h) se dají nahradit jediným, s jádrem g * h. Tak vznikne derivace gaussovky v oddílu o hranách. Má-li signál x konečnou délku |x| a jádro h délku |h|, stojí přímý výpočet |x| · |h| násobení.

Nechť f_s > 0 je vzorkovací frekvence (vzorků za sekundu) a x : ℤ → ℝ signál, pro který ∑_n |x[n]| < ∞. Jeho Fourierova transformace v diskrétním čase je funkce 𝓕(x) : ℝ → ℂ,

𝓕(x)(f) = ∑_{n ∈ ℤ} x[n] · exp(−2πi fn/f_s),   f ∈ ℝ

(spojitá obdoba diskrétní Fourierovy transformace X_k konečného signálu x₀, …, x_N−1 z dílu Fourierova transformace; f je frekvence v Hz). Nechť h je jádro (konečně mnoho nenulových členů). Frekvenční odezva filtru s jádrem h je funkce H = 𝓕(h) : ℝ → ℂ,

H(f) = ∑_{k ∈ ℤ} h[k] · exp(−2πi fk/f_s).

Číslo |H(f)| je zesílení a arg H(f) fázový posun, které filtr způsobí sinusovce s frekvencí f (věta níže).

Konvoluční věta: za uvedených předpokladů platí pro každé f ∈ ℝ

𝓕(x * h)(f) = 𝓕(x)(f) · H(f),

tj. Fourierova transformace konvoluce je součin Fourierových transformací.

Filtr se nazývá dolní propust, když propouští frekvence pod jistou mezní frekvencí f_c (|H(f)| je blízko 1 pro |f| < f_c) a tlumí frekvence nad ní (|H(f)| je blízko 0 pro f_c < |f| ≤ f_s/2); jak přesně „blízko“, je věcí návrhu filtru, ne definice. Obráceně (tlumí nízké, propouští vysoké) je to horní propust.

Nechť f_s > 0, h je jádro s frekvenční odezvou H, f ∈ ℝ a x je komplexní sinusovka s frekvencí f,

x[n] = exp(2πi fn/f_s),   n ∈ ℤ.

Pak y = x * h splňuje y[n] = H(f) · x[n] pro všechna n ∈ ℤ.

Důkaz.

y[n] = ∑_k h[k] exp(2πi f(n − k)/f_s) = exp(2πi fn/f_s) · ∑_k h[k] exp(−2πi fk/f_s) = H(f) · x[n]. ∎

Výstup je tatáž sinusovka vynásobená komplexním číslem H(f): |H(f)| ji zesílí nebo zeslabí, arg H(f) posune. Protože se každý signál dá rozložit na sinusovky (díl o Fourierově transformaci) a konvoluce je lineární, filtr prostě vynásobí každou frekvenci svým H(f) — to je konvoluční věta.

f_s = 500 Hz, M = 25, f = 50 Hz. Do okna 25/500 s = 0,05 s se vejde 50 · 0,05 = 2,5 periody — necelý počet, takže nula to není. Barvy podle vzorce výše:

Čitatel: |sin(π · 50 · 25/500)| = |sin(2,5π)| = 1. Jmenovatel: 25 · sin(π · 50/500) = 25 · sin(0,1π) = 25 · 0,309 = 7,73.

|H(50)| = 1/7,73 = 0,13

— přesně hodnota z textu výše. Pro 40 Hz se vejdou právě 2 periody a čitatel je sin(2π) = 0.

Nechť L ∈ ℕ, L ≥ 1, a e[1], …, e[L] ∈ ℝ. Efektivní hodnota (RMS, root mean square) této posloupnosti je číslo

RMS(e) = √1/L ∑_n=1^L e[n]² ≥ 0,

tedy odmocnina z průměru druhých mocnin. Je to „typická velikost“ kmitání, které má průměr nula — obyčejný průměr by se vynuloval. U šumu s nulovou střední hodnotou je to (výběrová) směrodatná odchylka; všechny „odchylky“ a „zbylý šum“ v tomto dílu jsou efektivní hodnoty. Sinusovka s amplitudou A má (přes celé periody) efektivní hodnotu A/√2, protože průměr sin² přes periodu je 1/2.

Brum před filtrem: 0,3/√2 = 0,212.

Brum po filtru: M = 10 při 500 Hz je okno 0,02 s a do něj se vejde 50 · 0,02 = 1 perioda brumu — přesně celá, čitatel sin(π) = 0 a brum zmizí.

Hrot QRS: složka 25 Hz se do okna vejde jen napůl (0,5 periody); barvy podle vzorce |H(f)| v oddílu o frekvencích:

|H(25)| = |sin π/2|/10 · sin π/20 = 1/10 · 0,156 = 0,64.

Proto hrot klesne — z 0,98 na 0,81.

Nechť x : ℤ → ℝ je signál. Jeho diference je posloupnost y daná pro každé n ∈ ℤ předpisem

y[n] = x[n] − x[n − 1].

Je to konvoluce y = x * d s jádrem d[0] = 1, d[1] = −1 a d[k] = 0 pro ostatní k, krátce d = (1, −1). Odhaduje derivaci signálu (změnu na jeden krok vzorkování) a je to horní propust (věta níže).

Nechť s > 0 (šířka ve vzorcích) a r = ⌈4s⌉. Derivace gaussovky je jádro

h[k] = −k · exp(−k²/2s²)/C   pro k = −r, …, r,   h[k] = 0 jinak,

kde C = ∑_k=1^r k · exp(−k²/2s²) > 0 je normovací konstanta. Je to (až na kladný násobek) navzorkovaná derivace gaussovky g(t) = exp(−t²/2s²), protože g′(t) = −t/s² g(t). Až na normování se přibližně rovná jádru d * g — diferenci zkonvolvované s navzorkovanou gaussovkou — takže signál se zároveň vyhladí i derivuje.

Nechť (e[n])_{n ∈ ℤ} jsou nezávislé náhodné veličiny, které mají všechny tentýž rozptyl σ², a h je jádro (konečně mnoho nenulových členů). Pak pro každé n ∈ ℤ

var((e * h)[n]) = σ² ∑_k h[k]²,

tedy jádro změní směrodatnou odchylku nezávislého šumu ze σ na σ · √∑_k h[k]². Speciálně diference d = (1, −1) dá rozptyl σ² + σ² = 2σ².

Důkaz. (e * h)[n] = ∑_k h[k] e[n − k] je součet nezávislých veličin; rozptyl součtu nezávislých veličin je součet rozptylů (i u rozdílu je to součet) a var(cX) = c² var(X). ∎

Diference: σ√2 = 0,2 · 1,414 = 0,28 (změřeno 0,27), a hrana výšky 1 dá vrchol 1 — vrchol je jen 3,5krát nad šumem. Derivace gaussovky níže: √∑_k h[k]² = 0,47, šum 0,2 · 0,47 = 0,095 (změřeno 0,10), vrchol hrany 1 — přes 10krát nad šumem.

Ve frekvencích: diference má odezvu |H(f)| = 2|sin(π f/f_s)|, nula pro f = 0 a 2 na Nyquistově frekvenci f_s/2. Konstantu zahodí a nejrychlejší kmitání zdvojnásobí — proto horní propust.

Nechť M ∈ ℕ je liché, M = 2m + 1 s m ∈ {0, 1, 2, …}. Pro reálná čísla a₁, …, a_M označme a₍₁₎ ≤ a₍₂₎ ≤ … ≤ a_(M) tatáž čísla seřazená vzestupně (se zachováním opakování). Medián je prostřední z nich,

med(a₁, …, a_M) = a_(m+1).

Nechť x : ℤ → ℝ je signál. Mediánový filtr délky M z něj vytvoří signál y daný pro každé n ∈ ℤ předpisem

y[n] = med(x[n − m], …, x[n + m]).

Každý vzorek se tedy nahradí mediánem hodnot v okně kolem něj (okraje konečného signálu doplňujeme krajní hodnotou jako u průměru). Není to konvoluce — je nelineární (věta níže).

Okno pěti vzorků s jedním výstřelem: (0,1; 0,2; 3,1; 0,3; 0,2). Průměr = 3,9/5 = 0,78 — výstřel ho vytáhl čtyřikrát výš, než leží ostatní vzorky. Seřazeno: (0,1; 0,2; 0,2; 0,3; 3,1), medián je prostřední, 0,2. Výstřel by mohl být 3 nebo 3 000 a medián by se nezměnil.

Podle definice je výstup (m + 1). nejmenší hodnota okna. Dokud je v okně nejvýš m výstřelů, leží výstup mezi hodnotami „slušných“ vzorků; přebít medián dokážou až tři výstřely v okně pěti vzorků.

Nechť M ≥ 3 je liché. Pak mediánový filtr délky M není lineární zobrazení: existují signály x, z, pro které medián součtu x + z není součtem mediánů. Proto neexistuje jádro h, pro které by mediánový filtr byl konvolucí s h.

Důkaz. Stačí jedno okno. Nechť okno signálu x obsahuje m jedniček a m + 1 nul a okno signálu z jedinou jedničku na místě, kde má x nulu, jinak nuly. Nul je v obou oknech většina, takže oba mediány jsou 0; okno x + z má ale m + 1 jedniček, tedy většinu, a medián 1 ≠ 0 + 0. Pro M = 3: med(0, 0, 1) = 0 a med(1, 0, 0) = 0, ale med((0, 0, 1) + (1, 0, 0)) = med(1, 0, 1) = 1. Každá konvoluce je lineární, takže medián konvolucí být nemůže — a nemá ani frekvenční odezvu. ∎

Souvislost se ztrátami: pro čísla x₁, …, x_M je průměr 1/M ∑ x_j to jediné c ∈ ℝ, které minimalizuje ∑_j (x_j − c)², a medián (pro liché M) to jediné c, které minimalizuje ∑_j |x_j − c|. Kvadrát velké odchylky je obří, takže jeden výstřel průměr strhne; absolutní hodnota roste jen lineárně a výstřel „táhne“ stejnou silou, ať je jakkoli velký.

Barvy podle vzorce výše. Dvě kostky:

P(X + Y = 7) = ∑_k=1⁶ 1/6 · 1/6 = 6/36

— dvojice (1, 6), (2, 5), (3, 4), (4, 3), (5, 2), (6, 1). Rozdělení dvou kostek v šestatřicetinách: 1, 2, 3, 4, 5, 6, 5, 4, 3, 2, 1 pro součty 2 až 12.

Tři kostky: P(10) = ∑_k=1⁶ P(třetí = k) · P(první dvě dají 10 − k). Součty 9, 8, 7, 6, 5, 4 mají na dvou kostkách 4, 5, 6, 5, 4 a 3 případy z 36, takže

P(10) = 1/6 · 4 + 5 + 6 + 5 + 4 + 3/36 = 27/216.

Nechť I : ℤ² → ℝ je obrázek (I[i, j] je jas pixelu v řádku i a sloupci j; obrázek konečné velikosti doplníme mimo jeho okraj) a nechť r ∈ {0, 1, 2, …} a K : {−r, …, r}² → ℝ je jádro velikosti (2r + 1) × (2r + 1). Dvourozměrná konvoluce obrázku I s jádrem K je obrázek I * K daný pro všechna (i, j) ∈ ℤ² předpisem

(I * K)[i, j] = ∑_a=−r^r ∑_b=−r^r K[a, b] · I[i − a, j − b].

Je to totéž co v jednom rozměru, jen se jádro posouvá po řádcích i sloupcích. Například jádro 3 × 3 s řádky (1, 0, −1), (2, 0, −2), (1, 0, −1) (Sobelův operátor) je diference ve vodorovném směru spojená s vyhlazením (1, 2, 1) ve svislém — dvourozměrná obdoba derivace gaussovky.

Cena: přímý výpočet stojí na každý pixel (2r + 1)² násobení. U jednorozměrného signálu délky N a jádra délky L je to N · L násobení, přes FFT řádově N log₂ N — vyplatí se to, když je jádro dlouhé (ozvěna sálu má desítky tisíc vzorků), u krátkých jader neuronových sítí (3 × 3) ne.