====== Fourierova transformace ====== V tomto cvičení si zopakujeme Fourierovu transformaci (FT) a její vlastnosti. Začneme s příklady na 1D a 2D FT, potom přejdeme do diskrétního prostoru. FT si vyzkoušíme v MATLABu na úloze filtrace obrázků. Pro zopakování Fourierovy transformace si přečtěte [[http://cs.wikipedia.org/wiki/Fourierova_transformace|wiki]] a [[http://mathworld.wolfram.com/FourierTransform.html|MathWorld]] nebo si projděte [[http://betterexplained.com/articles/an-interactive-guide-to-the-fourier-transform/|interaktivní prezentaci]]. ===== Zadání ===== * Spočítejte příklady označené červeně v části [[#Dvourozměrná Fourierova transformace]], kde každý bude hodnocen 1 bodem. [2b] * [[#Filtrace ve frekvenční oblasti]] - Vytvořte v MATLABu funkci pro filtraci dolnopropustním filtrem pomocí násobení ve frekvenčním spektru a aplikujte ji na obrázky {{courses:a6m33zsl:ct_brain.jpg?linkonly|CT_brain.jpg}} a {{courses:a6m33zsl:ct_lungs.jpg?linkonly|CT_lungs.jpg}} . Vytvořenou funkci odevzdejte a do zprávy vložte několik ilustračních obrázků pro rozptyly filtru \sigma = (10, 50, 100). [2b] * [[#Filtrace pomocí konvoluce]] - Vytvořte v MATLABu funkci která filtruje dolnopropustním filtrem obrázek pomocí konvoluce v obrazovém prostoru a aplikujte ji na obrázky {{courses:a6m33zsl:ct_brain.jpg?linkonly|CT_brain.jpg}} a {{courses:a6m33zsl:ct_lungs.jpg?linkonly|CT_lungs.jpg}}. Výsledek by se měl co nejvíce přibližoval k filtraci ve frekvenčním prostoru (z předešlé úlohy); určete vztah parametrů. Vytvořenou funkci odevzdejte a do zprávy vložte několik ilustračních obrázků pro různé rozptyly filtru. [1b] Do zprávy vždy stručně napište postup, vložte výsledné obrázky a k jednotlivým obrázkům pak napište zvolené parametry. ===== Jednorozměrná Fourierova transformace (Opakování) ===== Fourierova transformace je vyjádření časově závislého signálu pomocí harmonických signálů, tj. funkcí sin a cos. Slouží pro převod signálů z časové oblasti do oblasti frekvenční. Tato část slouží pro opakování, řešení příkladů do zprávy nepište. ==== Vlastnosti ==== Nechť f(t) \overset{\mathcal{F}}{\leftrightarrow} F(\omega) \, , \; g(t) \overset{\mathcal{F}}{\leftrightarrow} G(\omega), kde \mathcal{F} je Fourierova transformace. Pak platí: * Linearita: a \cdot f(t) + b \cdot g(t) \overset{\mathcal{F}}{\leftrightarrow} a \cdot F(\omega )+b \cdot G(\omega ) . * Symetrie: F(\omega) \overset{\mathcal{F}}{\leftrightarrow} F^*(-\omega), F^* je funkce komplexně sdružená k F . * Posun: f(t-t_0) \overset{\mathcal{F}}{\leftrightarrow} e^{-j\omega t_0}F(\omega ) . * Změna měřítka: f(t/a) \overset{\mathcal{F}}{\leftrightarrow} |a| \cdot F(a \omega ) . * Derivace: \frac{\partial ^n}{\partial t^n}f(t) \overset{\mathcal{F}}{\leftrightarrow} (j\omega)^nF(\omega ) . * Integrál: \int_{-\infty}^t f(t)\,{\rm d}t \overset{\mathcal{F}}{\leftrightarrow} = \frac {F(\omega)}{j\omega}+2\pi\delta(\omega)c kde \int_{-\infty}^{\infty} f(t)-c \,{\rm d}t =0 * Parsevalova věta (norma): \int_{-\infty}^{\infty}{|f(t)|^2\partial t}=\int_{-\infty}^{\infty}{|\mathcal{F}\big(F(\omega )|^2\partial \omega} . * Parsevalova věta (skalární souèin): \int_{-\infty}^{\infty}{f(t)g^*(t)\partial t}=\int_{-\infty}^{\infty}{\mathcal{F}\big(F(\omega )\mathcal{F}\big(g^*(t)\big)\partial \omega} . * Konvoluce: f(t)\ast g(t) \overset{\mathcal{F}}{\leftrightarrow} F(\omega ) \cdot G(\omega ) * Dualita: je-li f(t) \overset{\mathcal{F}}{\leftrightarrow} F(\omega) pak F(t) \overset{\mathcal{F}}{\leftrightarrow} 2\pi f(-\omega). Ještě budete potřebovat FT Gaussiánu: e^{-a t^2} \overset{\mathcal{F}}{\leftrightarrow} \sqrt{\frac{\pi}{a}} \cdot e^{-\pi^2 \frac{\omega^2}{a}} ==== Vypočítejte ==== **Příklad 1:** Vypočítejte FT pro Diracův impulz \delta (t) = \left\{\begin{matrix} \infty & t=0 \\ 0 & t \neq 0 \end{matrix}\right. , /* když víte že platí */ \int_{-\infty}^{\infty} \delta (t) \cdot f(t) \; dt = f(0). /* Řešení: \mathcal{F}\left (\delta (t) \right )(\omega) = \int_{-\infty}^{\infty} \delta (t) \cdot e^{j\omega t} \; dt = e^{j\omega \cdot 0} = 1 */ **Příklad 2:** Vypočítejte FT pro jednotkový skok H (t) = \left\{\begin{matrix} 0 & t<0 \\ 1 & t \geq 0 \end{matrix}\right. /* , když víte že platí H(t) = \int_{-\infty}^t \delta(\tau) \; d\tau. */ /* Řešení: \mathcal{F}\left (H(t) \right )(\omega) = \mathcal{F}\left (\int_{-\infty}^t \delta(\tau) \; d\tau \right )(\omega) = \frac{1}{j\omega} \cdot \mathcal{F}(\delta)(\omega) + 2 \pi / 2 \delta(\omega) = \frac{1}{j\omega} \cdot 1 + \pi \delta(\omega) */ //Poznámka:// funkce H(t) ani její obraz nejsou integrovatelné, výpočet s využitím výše zmíněných symbolických pravidel však je použitelný. **Příklad 3:** Vypočítejte FT pro \delta (t-t_0). **Příklad 4:** Vypočítejte FT pro jednotkový obdélníkový impuls I (t) = \left\{\begin{matrix} 1 & |t|<\frac{1}2{} \\ 0 & otherwise \end{matrix}\right. . **Příklad 5:** Vypočítejte FT pro \sin (t) /* = \frac{1}{2j}e^{jt} - \frac{1}{2j}e^{-jt}. */ **Příklad 6:** Vypočítejte FT pro \cos (t) /* = \frac{1}{2j}e^{jt} + \frac{1}{2j}e^{-jt}. */ ===== Dvourozměrná Fourierova transformace ===== /* Nechť f(t)\colon\mathbb{R}^2\to\mathbb{R}, pak F(u,v)=\int_{-\infty}^{\infty}\int_{-\infty}^{\infty}{f(x,y)e^{-2\pi j(xu+yv)}{\rm d} x{\rm d}y je Fourierovým obrazem funkce f(x,y). Zpětná transformace: f(x,y)=\int_{-\infty}^{\infty}\int_{-\infty}^{\infty}{F(u,v)e^{2\pi j(xu+yv)}{\rm d} u{\rm d}y. . */ Vektorový zápis rovnice s použitím vektoru souřadnic \mathbf{x} a vektoru úhlových frekvencí \boldsymbol{\omega}: F(\boldsymbol{\omega})=\int_{\mathbb{R}^n}{f(t)e^{-j\boldsymbol{\omega}\cdot\mathbf{x}}{\rm d} \mathbf{x}}\, . Rozmyslete rovnice pro frekvence 2\pi\mathbf{f}=\boldsymbol{\omega}. Rovnice platí pro libovolný počet dimenzí. Inverze je f(\mathbf{x})=\frac{1}{(2\pi)^d}\int_{\mathbb{R}^n}{f(t)e^{j\boldsymbol{\omega}\cdot\mathbf{x}}{\rm d}\boldsymbol{\omega}}\, , kde //d// je počet dimenzí. ==== Vlastnosti ==== * Separabilita: pokud f(x,y)=f(x)f(y) pak \mathcal{F}\big(f(x,y)\big)=\mathcal{F}\big(f(x)\big)\mathcal{F}\big(f(y)\big). Této vlastnosti se využívá pro urychlení výpočtu FT - rozmyslete.\\ * Věta o centrálním řezu: řez 2D Fourierovou transformací obrazu pod úhlem \phi je 1D Fourierovou transformací projekce téhož obrazu pod úhlem \phi. /* Rovnicí se tento fakt dá zapsat: \hat{P}_\phi(u,v)=\mathcal{F}\big(o(\xi,\nu)\big) , kde \hat{P}_\phi(u,v) je řez 2D FT pod úhlem \phi a o(\xi,\nu) je projekce obrazu o. */ Vlastnosti uvedené výše pro 1D FT platí analogicky. ==== Vypočítejte ==== **Příklad 1:** Vypočítejte FT pro Diracův impuls \delta (x,y) = \left\{\begin{matrix} \infty & x=y=0 \\ 0 & otherwise \end{matrix}\right. , \int_{-\infty}^{\infty} \int_{-\infty}^{\infty} f(x,y) \delta (x,y) \; dx \, dy = f(0,0). /* \mathcal{F}_{2D}(\delta)(u,v) = \mathcal{F}(\delta)(u) \cdot \mathcal{F}(\delta)(v) = 1 \cdot 1 = 1 */ **Příklad 2:** Vypočítejte FT pro komplexní exponencielu: f(x,y) = e^{j(u_0 x + v_0 y)}. /* \mathcal{F}_{2D}(f)(u,v) = \mathcal{F}_{2D}(e^{j(u_0 x + v_0 y)})(u,v) = \mathcal{F}_{2D}(e^{ju_0 x} \cdot e^{jv_0 y})(u,v) \\ \;\;\;\;\; = \mathcal{F}(e^{ju_0 x})(u) \cdot \mathcal{F}(e^{jv_0 y})(v) = 2\pi \delta (u-u_0) \cdot 2\pi \delta (v-v_0) = 4\pi^2 \cdot \delta (u-u_0) \cdot \delta (v-v_0) */ **Příklad 3:** Vypočítejte FT pro harmonickou funkci: f(x,y) = \sin(u_0x+v_0y) /* , když víte že: \sin(u_0x+v_0y) = \frac{1}{2j} e^{j(u_0 x + v_0 y)} - \frac{1}{2j} e^{-j(u_0 x + v_0 y)}. */ **Příklad 4 *:** Vypočítejte FT pro harmonickou funkci: f(x,y) = \cos(u_0x+v_0y) /*, když víte že: \sin(u_0x+v_0y) = \frac{1}{2j} e^{j(u_0 x + v_0 y)} + \frac{1}{2j} e^{-j(u_0 x + v_0 y)}. */ **Příklad 5 *:** Vypočítejte FT pro 2D obdélníkovou funkci: rect_{2D}(x,y) = \left\{\begin{matrix} 1 & |x|<0.5, \; |y|<0.5, \; \\ 0 & otherwise \end{matrix}\right.. ===== Diskrétní Fourierova transformace - DFT ===== /* Diskrétní signál vzniká obyčejně jako měření nějaké (spojité) veličiny. Pro zpracování obrazu (nejen) je důležité, aby byla dodržena vzorkovací věta, tedy měřená veličina musí být frekvenčně omezená s maximálním kmitočtem menším nebo rovným polovině vzorkovací frekvence. f_{max}\leq \frac{f_{vz}}{2} . Pokud vzorkovací věta není splněna, nelze měřený signál rekonstruovat z měření. Tomuto fenoménu se říká aliasing. V praxi se používají antialiasingové filtry na vstupech měřicích členů, které frekvenčně omezí měřenou veličinu. */ Pro diskrétní signál f_d(x,y)\colon(0,1,\ldots N_x-1)\times (0,1,\ldots N_y-1)\to\mathbb{R} definujeme diskrétní Fourierovu transformaci jako F_d(u,v)=\sum_{x=0}^{N_x-1}{\sum_{y=0}^{N_y-1}{f_d(x,y)e^{-2\pi j(\frac{xu}{N_x}+\frac{yu}{N_y})}}} a zpětnou transformaci jako f_d(u,v)=\frac{1}{N_x N_y}\sum_{u=0}^{N_x-1}{\sum_{v=0}^{N_y-1}{F_d(x,y)e^{2\pi j(\frac{xu}{N_x}+\frac{yu}{N_y})}}} . //Opakování:// DFT předpokládá, že signál je periodický s periodou //N//. Je-li diskrétní sekvence získána vzorkováním, nesmí mít původní spojitý signál frekvence větší než polovina frekvence vzorkovací, jinak dojde k aliasingu. //Doplňující úloha:// Porovnejte analytické výsledky výše uvedených příkladů s numerickými výsledky získanými v Matlabu pomocí FFT. /* ==== Vztah DFT ke spojité FT ==== Abychom mohli reprezentovat transformovanou funkci konečným počtem frekvencí (bází, harmonických funkcí) musí být vzor periodický (vzpomeňte na Fourierovy řady). DFT tedy implicitně předpokládá, že je měřený signál periodický je, tedy že po posledním vzorku přijde opět vzorek první, a před prvním vzorkem předcházel poslední vzorek naší (měřené) posloupnosti. Je to tedy frekvenční omezenost a periodicita funkce které nám umožňují výpočet DFT diskrétního signálu. Implicitní periodicita signálu má ale i negativní vlastnost, pokud signál není svým vlastním periodickým prodloužením (např. \sin(x) na intervalu \langle 0,2\pi+1\rangle vyskytnou se ve FT signálu frekvence, které bychom nečekali, říká se tomu rozmývání spektra. */ ===== Filtrace pomocí FFT a konvoluce ===== V MATLABu je FFT implementována ve funkcích: //fft//, //ifft//, //fft2//, //ifft2// (přečtěte si "help [[http://www.mathworks.com/help/matlab/ref/fft2.html|fft2]]"). Stáhněte si obrázky {{courses:a6m33zsl:ct_brain.jpg?linkonly|CT_brain.jpg}} a {{courses:a6m33zsl:ct_lungs.jpg?linkonly|CT_lungs.jpg}} a šablonu požadovaných funkcí //{{courses:a6m33zsl:imfilterfft_lowpass.m?linkonly|imFilterFFT_lowpass}}// a //{{courses:a6m33zsl:imfilterconv_lowpass.m?linkonly|imFilterConv_lowpass}}//. Načtěte si každý z obrázků do MATLABu (pomocí [[http://www.mathworks.com/help/matlab/ref/imread.html|imread]]), konvertujte na černobílý (je-li potřeba) a do typu double, vypočítejte jeho Fourierovo spektrum a vykreslete jej. /* Fourierův obraz se skládá z komplexních čísel, které se dají složit z magnitudy a fáze: F(u,v) = |F(u,v)| \cdot e^{j \angle F(u,v) }, kde |F(u,v)| = \sqrt{ F_R(u,v)^2 + F_I(u,v)^2}, \angle F(u,v) = \tan^{-1}\frac{F_I(u,v)}{F_R(u,v)}. */ Při vizualizaci spektra nás zajímá magnituda, které nese informaci o síle frekvence (u,v). /* Skript který použijeme kresli výkonové spektrum */ Na vizualizaci můžete používat funkci {{courses:a6m33zsl:draw_fft.m?linkonly|draw_fft.m}}, která nakreslí obrázek s výkonovým Fourierovým spektrem |F(u,v)|^2 a aplikuje na něj logaritmus pro lepší využití rozsahu šedé škály. image_fft = fft2(image1); fig = draw_fft(image_fft); //__Poznámka:__ Při práci z frekvenčním spektru v MATLABe myslete na to, že není centrované, tj. nízké frekvence se nacházejí v rozích matice a vysoké frekvence uprostřed obrázku. Na centrovaní můžete použít příkaz [[http://www.mathworks.com/help/matlab/ref/fftshift.html|fftshift]], a na zpětnou transformaci [[http://www.mathworks.com/help/matlab/ref/ifftshift.html|ifftshift]]. Funkce fftshift pouze centruje spektrum (přesouvá jednotlivé kvadranty obrázku), neprovádí Fourierovu transformaci.// Obrázek CT_brain.jpg a jeho centrované výkonové spektrum: {{ :courses:a6m33zsl:fft_brain.png?700 |CT mozku a jeho spektrum}} ==== Filtrace ve frekvenční oblasti ==== Filtraci, neboli potlačení vysokých frekvencí můžeme provést ve Fourierově spektru vynásobením F(u,v) \cdot G_\sigma(u,v) kde //G// je dolnopropustní filter G_\sigma(u,v)= e^{\frac{-2\pi(u^2+v^2)}{\sigma^2}} o velikosti \sigma. Pro jednoduchost jsme vynechali normalizační koeficient. //Poznámka:// velikost filtru v prostorové oblasti je nepřímo úměrná \sigma, viz Fourierova transformace Gaussiánu výše. **Úloha:** Vytvořte v MATLABu funkci //{{courses:a6m33zsl:imfilterfft_lowpass.m?linkonly|[img] = imFilterFFT_lowpass( image, sigma )}}// pro filtraci Gaussovským filtrem pomocí násobení ve frekvenčním spektru a aplikujte ji na vstupní obrázky {{courses:a6m33zsl:ct_brain.jpg?linkonly|CT_brain.jpg}} a {{courses:a6m33zsl:ct_lungs.jpg?linkonly|CT_lungs.jpg}}. Vyzkoušejte několik rozptylů \sigma = (10, 50, 100), subjektivně vyberte nejlepší parametr tak, aby se odstranil sum z obrázků, ale zachovaly se hranice objektů. Jak parametr jste vybrali? Vložte filtrované obrázky. {{ :courses:a6m33zsl:fft_dolnipropust.png?800 |Dolní propusť pomocí FFT}} **Návod:** - Vypočítejte fourierovu transformaci //image_fft// vstupního obrázku. - Vytvořte matici //filter_fft// o stejné velikosti jako vypočítaná matice vstupního obrázku filter_fft = zeros(size(image_fft)); - Do matice uložte hodnoty filtru G_\sigma. Nulovou frekvenci na začátku zvolíme uprostřed a krok 1 (z definice DFT). * Můžete si vytvořit pomocné matice, které budou obsahovat frekvence pro každou položku: [sz1 sz2] = size(image_fft); fft_range1 = (-sz1/2):(sz1-sz1/2-1); fft_range2 = (-sz2/2):(sz2-sz2/2-1); [ids1,ids2] = ndgrid(fft_range1,fft_range2); dist2 = ids1.^2 + ids2.^2; * Vypočítáme hodnoty filtrupodle parametru //sigma//. filter_fft(:,:) = (sigma*sqrt(2*pi))^(-1) * exp(-dist2/(2*sigma^2)); /* filter_fft(:,:) = exp(-2*pi*dist2/GAUSS_SIGMA^2); */ * Posuneme spektrum z centrovaného formátu nazpět. (Podívejte se jak to vypadá.) filter_fft = ifftshift(filter_fft); - Uděláme filtraci násobením: image_filt_fft = image_fft .* filter_fft; draw_fft(image_filt_fft,3); - Vypočítáme výsledný obraz pomocí IFFT. Předpokládáme, že je reálný. image_filt_rec = ifft2(image_filt_fft); image_filt_rec = real(image_filt_rec); ==== Filtrace pomocí konvoluce ==== /* Nechť F(u,v) a G(u,v) jsou Fourierovy transformace spojitých obrazových funkci f(x,y) a g(x,y), potom */ Fourierova transformace konvoluce /* g(u,v) = f(u,v) * g(u,v) = \int_{-\infty}^{\infty} \int_{-\infty}^{\infty} f(\alpha,\beta) \cdot g(u-\alpha,v-\beta) \; d\alpha \, d\beta */ je součin Fourierových transformací: \mathcal{F}_{2D}(f*g)(u,v) = F(u,v) \cdot G(u,v) Tuto vlastnost využijeme na filtraci pomocí konvoluce s maskou g(x,y), ale v diskrétním prostoru. **Úloha:** Vytvořte v MATLABu funkci //{{courses:a6m33zsl:imfilterconv_lowpass.m?linkonly|[img] = imFilterConv_lowpass( image, sigma )}}// pro filtraci Gaussovským filtrem a aplikujte na obrázky{{courses:a6m33zsl:ct_brain.jpg?linkonly|CT_brain.jpg}} a {{courses:a6m33zsl:ct_lungs.jpg?linkonly|CT_lungs.jpg}}. Vytvořte konvoluční masku g_{\sigma}(x,y) pro filter G_{\sigma}(u,v). Proveďte konvoluci v obrazovém prostoru pomocou MATLAB-prikazu: [[http://www.mathworks.es/es/help/matlab/ref/conv2.html|conv2]]. Vypočítejte \sigma pro filtraci v obrazovém prostoru tak, aby se výsledek co nejvíce přibližoval k filtraci ve frekvenčním prostoru. Velikost filtru zvolte tak, aby velikost aproximační chyby byla přijatelné. Vložte do zprávy filtrované obrázky. {{ :courses:a6m33zsl:fft_konvoluce.png?800 |Konvoluce}} **Návod:** - Vytvořte si Gaussovskou masku //G// v obrazovém prostoru, tak aby odpovídala dříve (předešlá úloha) použitým maskám ve Foourierově spektru. Velikost masky zvolte tak, aby zanedbaná část Gaussiánu byla dostatečně malá. Střed Gaussiánu je uprostřed matice. - Filtr normalizujte na soucet 1. - Volitelně vykreslete //G// jako 2D povrch příkazem //mesh//. - Konvoluci aplikujeme pomocí funkce //conv2// v MATLABu. image_filt_conv = conv2(img,G,'same'); //Poznámka:// Pro filtraci v MATLABu se nejčastěji používají funkce[[http://www.mathworks.com/help/images/ref/fspecial.html|fspecial]] a [[http://www.mathworks.com/help/matlab/ref/filter2.html|filter2]]. {{tag>medical exercise Fourier FT filter}}