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 wiki a MathWorld nebo si projděte interaktivní prezentaci.

• 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 CT_brain.jpg a CT_lungs.jpg . Vytvořenou funkci odevzdejte a do zprávy vložte několik ilustračních obrázků pro rozptyly filtru <latex>\sigma = (10, 50, 100)</latex>. [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 CT_brain.jpg a 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.

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.

Nechť <latex>f(t) \overset{\mathcal{F}}{\leftrightarrow} F(\omega) \, , \; g(t) \overset{\mathcal{F}}{\leftrightarrow} G(\omega)</latex>, kde <latex>\mathcal{F}</latex> je Fourierova transformace. Pak platí:

• Linearita: <latex>a \cdot f(t) + b \cdot g(t) \overset{\mathcal{F}}{\leftrightarrow} a \cdot F(\omega )+b \cdot G(\omega ) </latex> .
• Symetrie: <latex>F(\omega) \overset{\mathcal{F}}{\leftrightarrow} F^*(-\omega)</latex>, <latex>F^*</latex> je funkce komplexně sdružená k <latex>F </latex> .
• Posun: <latex>f(t-t_0) \overset{\mathcal{F}}{\leftrightarrow} e^{-j\omega t_0}F(\omega ) </latex> .
• Změna měřítka: <latex>f(t/a) \overset{\mathcal{F}}{\leftrightarrow} |a| \cdot F(a \omega ) </latex> .
• Derivace: <latex>\frac{\partial ^n}{\partial t^n}f(t) \overset{\mathcal{F}}{\leftrightarrow} (j\omega)^nF(\omega ) </latex> .
• Integrál: <latex>\int_{-\infty}^t f(t)\,{\rm d}t \overset{\mathcal{F}}{\leftrightarrow} = \frac {F(\omega)}{j\omega}+2\pi\delta(\omega)c</latex> kde <latex>\int_{-\infty}^{\infty} f(t)-c \,{\rm d}t =0 </latex>
• Parsevalova věta (norma): <latex>\int_{-\infty}^{\infty}{|f(t)|^2\partial t}=\int_{-\infty}^{\infty}{|\mathcal{F}\big(F(\omega )|^2\partial \omega}</latex> .
• Parsevalova věta (skalární souèin): <latex>\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}</latex> .
• Konvoluce: <latex>f(t)\ast g(t) \overset{\mathcal{F}}{\leftrightarrow} F(\omega ) \cdot G(\omega ) </latex>
• Dualita: je-li <latex>f(t) \overset{\mathcal{F}}{\leftrightarrow} F(\omega)</latex> pak <latex>F(t) \overset{\mathcal{F}}{\leftrightarrow} 2\pi f(-\omega)</latex>.

Ještě budete potřebovat FT Gaussiánu: <latex>e^{-a t^2} \overset{\mathcal{F}}{\leftrightarrow} \sqrt{\frac{\pi}{a}} \cdot e^{-\pi^2 \frac{\omega^2}{a}}</latex>

Příklad 1: Vypočítejte FT pro Diracův impulz <latex>\delta (t) = \left\{\begin{matrix} \infty & t=0
0 & t \neq 0 \end{matrix}\right.</latex> , <latex>\int_{-\infty}^{\infty} \delta (t) \cdot f(t) \; dt = f(0)</latex>.

Příklad 2: Vypočítejte FT pro jednotkový skok <latex>H (t) = \left\{\begin{matrix} 0 & t<0
1 & t \geq 0 \end{matrix}\right.</latex>

Poznámka: funkce <latex>H(t)</latex> 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 <latex>\delta (t-t_0)</latex>.

Příklad 4: Vypočítejte FT pro jednotkový obdélníkový impuls <latex>I (t) = \left\{\begin{matrix} 1 & |t|<\frac{1}2{}
0 & otherwise \end{matrix}\right.</latex> .

Příklad 5: Vypočítejte FT pro <latex>\sin (t)</latex>

Příklad 6: Vypočítejte FT pro <latex>\cos (t)</latex>

Vektorový zápis rovnice s použitím vektoru souřadnic <latex>\mathbf{x}</latex> a vektoru úhlových frekvencí <latex>\boldsymbol{\omega}</latex>: <latex> F(\boldsymbol{\omega})=\int_{\mathbb{R}^n}{f(t)e^{-j\boldsymbol{\omega}\cdot\mathbf{x}}{\rm d} \mathbf{x}}\, </latex> . Rozmyslete rovnice pro frekvence <latex>2\pi\mathbf{f}=\boldsymbol{\omega}</latex>.

Rovnice platí pro libovolný počet dimenzí. Inverze je <latex> 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}}\, </latex>, kde d je počet dimenzí.

• Separabilita: pokud <latex>f(x,y)=f(x)f(y)</latex> pak <latex>\mathcal{F}\big(f(x,y)\big)=\mathcal{F}\big(f(x)\big)\mathcal{F}\big(f(y)\big)</latex>. 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 <latex>\phi</latex> je 1D Fourierovou transformací projekce téhož obrazu pod úhlem <latex>\phi</latex>.

Vlastnosti uvedené výše pro 1D FT platí analogicky.

Příklad 1: Vypočítejte FT pro Diracův impuls <latex>\delta (x,y) = \left\{\begin{matrix} \infty & x=y=0
0 & otherwise \end{matrix}\right.</latex> , <latex>\int_{-\infty}^{\infty} \int_{-\infty}^{\infty} f(x,y) \delta (x,y) \; dx \, dy = f(0,0)</latex>.

Příklad 2: Vypočítejte FT pro komplexní exponencielu: <latex>f(x,y) = e^{j(u_0 x + v_0 y)}</latex>.

Příklad 3: Vypočítejte FT pro harmonickou funkci: <latex>f(x,y) = \sin(u_0x+v_0y)</latex>

Příklad 4 *: Vypočítejte FT pro harmonickou funkci: <latex>f(x,y) = \cos(u_0x+v_0y)</latex>

Příklad 5 *: Vypočítejte FT pro 2D obdélníkovou funkci: <latex>rect_{2D}(x,y) = \left\{\begin{matrix} 1 & |x|<0.5, \; |y|<0.5, \;
0 & otherwise \end{matrix}\right.</latex>.

Pro diskrétní signál <latex>f_d(x,y)\colon(0,1,\ldots N_x-1)\times (0,1,\ldots N_y-1)\to\mathbb{R}</latex> definujeme diskrétní Fourierovu transformaci jako <latex> 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})}}} </latex> a zpětnou transformaci jako <latex> 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})}}}</latex> .

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.

V MATLABu je FFT implementována ve funkcích: fft, ifft, fft2, ifft2 (přečtěte si “help fft2”).

Stáhněte si obrázky CT_brain.jpg a CT_lungs.jpg a šablonu požadovaných funkcí imFilterFFT_lowpass a imFilterConv_lowpass. Načtěte si každý z obrázků do MATLABu (pomocí imread), konvertujte na černobílý (je-li potřeba) a do typu double, vypočítejte jeho Fourierovo spektrum a vykreslete jej. Při vizualizaci spektra nás zajímá magnituda, které nese informaci o síle frekvence <latex>(u,v)</latex>. Na vizualizaci můžete používat funkci draw_fft.m, která nakreslí obrázek s výkonovým Fourierovým spektrem <latex>|F(u,v)|^2</latex> 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 fftshift, a na zpětnou transformaci 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:

Filtraci, neboli potlačení vysokých frekvencí můžeme provést ve Fourierově spektru vynásobením <latex>F(u,v) \cdot G_\sigma(u,v)</latex> kde G je dolnopropustní filter <latex>G_\sigma(u,v)= e^{\frac{-2\pi(u^2+v^2)}{\sigma^2}}</latex> o velikosti <latex>\sigma</latex>. Pro jednoduchost jsme vynechali normalizační koeficient.

Poznámka: velikost filtru v prostorové oblasti je nepřímo úměrná <latex>\sigma</latex>, viz Fourierova transformace Gaussiánu výše.

Úloha: Vytvořte v MATLABu funkci [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 CT_brain.jpg a CT_lungs.jpg. Vyzkoušejte několik rozptylů <latex>\sigma = (10, 50, 100)</latex>, 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.

Návod:

1. Vypočítejte fourierovu transformaci image_fft vstupního obrázku.
2. Vytvořte matici filter_fft o stejné velikosti jako vypočítaná matice vstupního obrázku

   filter_fft = zeros(size(image_fft));

3. Do matice uložte hodnoty filtru <latex>G_\sigma</latex>. 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));

   • Posuneme spektrum z centrovaného formátu nazpět. (Podívejte se jak to vypadá.)

     filter_fft = ifftshift(filter_fft);

4. Uděláme filtraci násobením:

   image_filt_fft = image_fft .* filter_fft;
   draw_fft(image_filt_fft,3);

5. 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);

Fourierova transformace konvoluce je součin Fourierových transformací: <latex>\mathcal{F}_{2D}(f*g)(u,v) = F(u,v) \cdot G(u,v)</latex>

Tuto vlastnost využijeme na filtraci pomocí konvoluce s maskou <latex>g(x,y)</latex>, ale v diskrétním prostoru.

Úloha: Vytvořte v MATLABu funkci [img] = imFilterConv_lowpass( image, sigma ) pro filtraci Gaussovským filtrem a aplikujte na obrázkyCT_brain.jpg a CT_lungs.jpg. Vytvořte konvoluční masku <latex>g_{\sigma}(x,y)</latex> pro filter <latex>G_{\sigma}(u,v)</latex>. Proveďte konvoluci v obrazovém prostoru pomocou MATLAB-prikazu: conv2. Vypočítejte <latex>\sigma</latex> 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.

Návod:

1. 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.
2. Filtr normalizujte na soucet 1.
3. Volitelně vykreslete G jako 2D povrch příkazem mesh.
4. 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í funkcefspecial a filter2.