2. Komputerowa obróbka danych

2.7. Filtracje częstotliwościowe

Alternatywnym do przestrzennego sposobem filtracji są metody częstotliwościowe, co wynika przede wszystkim ze skutecznej formy analizy sygnałów oraz wygodnych metod projektowania filtrów za pomocą reprezentacji częstotliwościowej. Wykorzystuje się tutaj fakt, że liniowej filtracji splotowej w przestrzeni obrazu jest  równoważna filtracja realizowana w fourierowskiej dziedzinie częstotliwościowej jako iloczyn transformat Fouriera sygnału i odpowiedzi impulsowej filtru.

W podrozdziale 2.2, zdefiniowano trygonometryczne szeregi fourierowskie (2.1) w postaci trygonometrycznej oraz warunki dokładnego opisu sygnałów (warunki Dirichleta) za pomocą takiej reprezentacji. 

Mając sygnał  T okresowy  f(t) (tj.  f(t)=f(t+T)) dla t \in (-\infty,\infty),  można go przedstawić korzystając z syntetycznej, zespolonej postaci wykładniczej szeregów Fouriera (przypomnienie (2.2)

f(t)= \sum_{k \in \mathbb{Z}} b_k e^{j i \omega_0 t}

(4.44) 

gdzie b_i=1/T\int\limits_{-T/2}^{T/2}f(t)e^{-j k \omega_0 t}dt.

Transformacja Fouriera jest rozszerzeniem koncepcji wykładniczych szeregów Fouriera na klasę wszystkich sygnałów f \in \mathrm{L}^2(\mathbb{R}), w tym nieokresowych. Jest więc to przekształcenie \mathcal{F}: \mathrm{L}^2(\mathbb{R}) \rightarrow \mathrm{L}^2(\mathbb{R}) zakładające szczególne rozumienie sygnałów nieokresowych. 

Jeśli f(t) jest nieokresowy i rozciąga się w całej dziedzinie liczb rzeczywistych, wtedy traktowany jest jako granicę skończonych sygnałów okresowych z T\rightarrow \infty, gdzie \omega_0 \rightarrow 0, co prowadzi do modelowego uciąglenia częstotliwości k \omega_0 \rightarrow \omega. Stąd zamiast sumy w (4.44) należy zapisać całkę, co prowadzi do ciągłego zbioru współczynników F(\omega), takich że 

F(\omega)=\int_{-\infty}^{\infty} f(t)cos \omega t dt + j \int_{-\infty}^{\infty} f(t)sin \omega t dt = \int_{-\infty}^{\infty} f(t)e^{-j\omega t} dt

(4.45) 

Wtedy sygnał reprezentowany jest jako  

f(t)=\frac{1}{2 \pi}\int_{-\infty}^{\infty} F(\omega)e^{j\omega t} d\omega

(4.46) 

Wyrażenie (4.45) jest transformacją Fouriera dająca ciągłe widmo częstotliwościowe sygnału, zaś (4.46) jest odwrotnym przekształceniem Fouriera.

Ponieważ funkcje rozwinięcia fourierowskiego, w tym przypadku funkcje bazowe transformacji Fouriera są zespolone, mamy więc transformatę F(\omega) \in \mathbb{C}, czyli zespolone widmo fourierowskie z częścią rzeczywistą \Re (F(\omega)) i urojoną  \Im (F(\omega)), amplitudą |F(\omega)| i fazą \angle F(\omega):

F(\omega)=\Re (F(\omega)) + \Im (F(\omega))=|F(\omega)|e^{j\angle F(\omega)}

(4.47) 

gdzie  |F(\omega)|=\sqrt{\Re^2(F(\omega)) + \Im^2 (F(\omega))} oraz \angle F(\omega) = \arctan \frac{\Im (F(\omega))}{\Re (F(\omega))}.
 
W przypadku sygnałów nieokresowych, określonych w skończonym przedziale -[T/2,T/2] przyjmuje się okresowe (cykliczne) rozszerzenie (powielenie) tego sygnału na przedział (-\infty,\infty) z okresem T. Wtedy, podobnie jak dla sygnałów okresowych rozciągniętych na całą dziedzinę \mathbb{R}, mając klasę sygnałów f(t) \in \mathrm{L}^2_T transformacja przyjmuje postać (analogicznie do (4.50)

F(u)=1/T\int\limits_{-T/2}^{T/2}f(t)e^{-j u \frac{2\pi}{T} t}dt

(4.48) 

gdzie F(u) \in \mathrm{l}^2(\mathbb{Z}), zaś odwrotna

f(t)= \sum_{u \in \mathbb{Z}} F(u) e^{j u \frac{2\pi}{T} t}

(4.49) 

W przypadku sygnałów dyskretnych f(k) \in \mathrm{l}^2(\mathbb{Z}), rozumianych analogicznie jako okresowe w sensie przejścia granicznego, transformacja jest iloczynem skalarnym sygnału z bazą funkcji sinusoidalnych 

F(\omega)=\sum_{k \in \mathbb{Z}} f(k) e^{-j k \omega}

(4.50) 

gdzie F({\omega}) \in \mathrm{L}^2([-\pi,\pi]), zaś odwrotna

f(k)= \frac{1}{2 \pi} \int_{- \pi}^{\pi} F(\omega) e^{j k \omega} d\omega

(4.51) 

Najbardziej praktyczny, wykorzystywany w realiach systemów cyfrowych, sposób liczenia transformacji Fouriera dla dyskretnych sygnałów opisanych skończonym ciągiem wartości (f(k))_{k=0}^{K-1} \in \mathbb{R}^K (z możliwością reprezentacji sygnałów zespolonych (f(k))_{k=0}^{K-1} \in \mathbb{C}^K) określają zależności formy prostej

F(u)= \sum_{k=0}^{K-1} f(k) e^{-\frac{j 2\pi k u}{K}},\quad u=0,\ldots, K-1 

(4.52) 

i odwrotnej

f(k) = \frac{1}{K} \sum_{u=0}^{K-1} F(u) e^{\frac{j 2\pi k u}{K}},\quad k=0,\ldots, K-1 

(4.53) 

Zakładamy przy tym cykliczne rozszerzenie f(k) z okresem K. Wyrażenia (4.52) oraz  (4.53) definiują tzw. dyskretną transformację Fouriera, w skrócie DFT (ang. Discrete Fourier Transform).

Dwuwymiarowa dyskretna transformacja Fouriera - 2D DFT, użyteczna w analizie częstotliwościowej obrazów postaci f(k,l), przyjmuje separowalną postać poprzez sekwencyjne złożenie przekształcenia wzdłuż jednego kierunku (opisanego współrzędną k danych obrazowych), a później drugiego (opisanego współrzędną l). Poprzez analogię do (4.52) przyjmuje ona postać  

F(u,v)= \sum_{k=0}^{K-1} \sum_{l=0}^{L-1} f(k,l) e^{-\frac{j 2\pi k u}{K}} e^{-\frac{j 2\pi l v}{L}}  

(4.54) 

gdzie u=0,\ldots, K-1,\ v=0,\ldots, L-1.
Ze względu na brak informacji lokalnych w fourierowskim widmie obrazu, w wielu zastosowaniach przekształcenie 2D DFT (lub zwykle jego szybsza wersja 2D FFT - Fast Fourier Transform) realizowane jest przy blokowym podziale dziedziny. Transformata obliczana jest w każdym z bloków niezależnie, co pozwala zachować lokalny (o lokalności decyduje wielkość elementarnego bloku z pokrycia dziedziny) charakter analizy częstotliwościowych właściwości obrazu. Nawiązuje to do koncepcji okienkowej transformacji Fouriera. W zależności od lokalnych cech widma sygnału można stosować dobrane metody filtracji.

Operacja filtracji obrazów w dziedzinie częstotliwości zakłada wykorzystanie filtrów cyfrowych o dobranych, selektywnych właściwościach, np. redukcji składowych wysokoczęstotliwościowych (dolnoprzepustowy filtr wygładzający) czy zachowania jedynie lokalnej informacji o charakterze szybkozmiennym (górnoprzepustowy filtr pasmowy). Stosowane są także filtry pasmowo-przepustowe, zachowujące cechy sygnału reprezentowane przez określony zakres widma (tj. częstotliwości pomiędzy dolną i górną częstotliwością graniczną lub wokół wybranej częstotliwości środkowej, z określoną szerokością pasma przepustowego). Filtry o dokładnie odwrotnych właściwościach, tłumiące sygnał w ustalonym zakresie częstotliwości nazywane są pasmowo-zaporowymi. W przypadku operacji lokalnych są to filtry o skończonej odpowiedzi impulsowej - SOI. 

Charakterystyka częstotliwościowa filtru dyskretnego dokonywana jest za pomocą transformacji Fouriera jego skończonej odpowiedzi impulsowej \{h_n\}_{n=0}^{N-1}, tak że

H(u)= \sum_{n=0}^{N-1} h_n e^{-\frac{j 2\pi n u}{N}} 

(4.55) 

gdzie H(u) nazywana jest funkcją przenoszenia (inaczej transmitancją) filtru. Analogicznie do (4.54) oraz (4.55) określana jest dwuwymiarowa transmitancja H(u,v) filtrów stosowanych do przetwarzania obrazów.

Liniowej filtracji splotowej g(k,l)=f*h(k,l), realizowanej w przestrzeni obrazu za pomocą filtru h_{m,n} według (4.18), odpowiada (jest równoważna - zgodnie z twierdzeniem o splocie: \mathcal{F}\{f*h\}=\mathcal{F}\{f\}\cdot \mathcal{F}\{h\}  filtracja realizowana w fourierowskiej dziedzinie częstotliwościowej za pomocą iloczynu transformaty Fouriera sygnału i transmitancji filtru 

G(u,v)=F(u,v) \cdot H(u,v)

(4.56) 

Stąd obraz przetworzony

g(k,l) =\mathcal{F}^{-1} \{F(u,v) \cdot H(u,v)\} 

(4.57) 

Dyskretna postać funkcji przenoszenia filtru ma duże znaczenie implementacyjne, jednak w wielu przypadkach wygodniej jest analizować właściwości filtrów (o potencjalnie nieograniczonej odpowiedzi impulsowej - NOI) za pomocą ich ciągłej charakterystyki częstotliwościowej (analogicznie do (4.50)

H(\omega)=\sum_{n \in \mathbb{Z}} h_n e^{-j n \omega}

(4.58) 

a w przypadku dwuwymiarowym

H(\omega_x,\omega_y)=\sum_{m \in \mathbb{Z}} \sum_{n \in \mathbb{Z}} h_{m,n} e^{-j m \omega_x} e^{-j n \omega_y}

(4.59) 

Idealny filtr dolnoprzepustowy o pulsacji granicznej (odcięcia) \omega_0 jest opisany przez moduł transmitancji (charakterystyka amplitudowa) postaci 

|H(\omega)|= \left \{ \begin{array}{ll} 1 & \textrm{dla $\omega \leq \omega_0$}\\ 0 & \textrm{w p.p.} \end{array} \right.

(4.60) 

Daje on efekt przepuszczenia bez żądnych zmian wszystkich składowych sygnału o częstotliwościach z użytecznego zakresu widma (pasmo przepustowe), przy całkowitym usunięciu  pozostałych (z pasma zaporowego), bez pasma przejściowego (zakres częstotliwości częściowo tłumionych, leżących na granicy zakresów przepustowego i zaporowego). 

Transmitancja dolnoprzepustowego filtru idealnego w zakresie pasma przepustowego ma więc postać H(\omega)=e^{-j 2\Pi \frac{\omega}{\omega_0}} (w paśmie zaporowym jest zerowa), przy liniowej fazie \angle H(\omega)=-  \frac{2\Pi}{\omega_0}\omega

Niestety, takich filtrów nie sposób wykorzystać w praktyce (jest nierealizowalny fizycznie) -- jego skończone widmo (tj. przyjmujące niezerowe wartości jedynie w domkniętym obszarze dziedziny) implikuje NOI o postaci funkcji \textrm{sinc} x=\sin x/x (nie jest zachowany warunek przyczynowości, bo odpowiedź impulsowa filtru nie ma końca). Obcinanie odpowiedzi impulsowej do postaci SOI powoduje odkształcenie zmierzonych według  (4.61) i (4.62) postaci H(d), pojawią się niekorzystne efekty zniekształceń fazowych. 

Większość uniwersalnych filtrów obrazowych ma symetryczną charakterystykę częstotliwościową, co upraszcza definiowanie ich transmitancji jako H(d), gdzie odległość od początku układu częstotliwościowych współrzędnych d=\sqrt{\omega_x^2+\omega_y^2}. Nie wyklucza to oczywiście możliwości dobierania filtrów o anizotropowej charakterystyce przepustowej H(\omega_x, \omega_y), dostosowanej do specyfiki przetwarzanych obrazów.

W przypadku idealnego filtru obrazowego, dolnoprzepustowe o pulsacji granicznej d_0 mamy więc analogiczną charakterystyka amplitudową postaci 

|H_i(d)|= \left \{ \begin{array}{ll} 1 & \textrm{dla $d \leq d_0$}\\ 0 & \textrm{w p.p.} \end{array} \right.

(4.61) 

Odwrotne zadanie przepuszczania wysokich częstotliwości przy eliminacji niskich można  uzyskać za pomocą analogicznego, idealnego filtru górnoprzepustowego:

|H_i(d)|= \left \{ \begin{array}{ll} 1 & \textrm{dla $d \geq d_0$}\\ 0 & \textrm{w p.p.} \end{array} \right.

(4.62) 

Stosowanie filtrów obrazowych o charakterystyce zbliżonej do gwałtownego skoku filtru idealnego jest źródłem - np. w filtracji dolnoprzepustowej - tzw. efektów pierścieniowych (ringing effect) powodowanych zjawiskiem Gibbsa (brakuje wyciętych z widma wysokich częstotliwości do reprezentacji krawędziowych skoków funkcji jasności), z powielaniem krawędzi i słabą redukcją szumu impulsowego (rys. 4.19). Takie właściwości próbuje się aproksymować za pomocą realnych filtrów o łagodniejszych charakterystykach  w zakresie pasma przejściowego, najlepiej o znikomych oscylacjach w pasmach przepustowym i zaporowym.

Rys. 4.19 Przykład filtracji w dziedzinie częstotliwościowej za pomocą dolnoprzepustowego filtru ''idealnego'' (aproksymowanego), z widocznym efektem pierścieniowym; kolejno w szeregu u góry obraz źródłowy, obraz po filtracji z wykorzystaniem transmitancji (4.61) przy rozmiarze d_0 równym 20% dziedziny obrazu oraz obraz po analogicznej filtracji ze zmniejszeniem względnego rozmiaru d_0 do 10%; w dolnym rzędzie zamieszczono   częstotliwościowe widma amplitudowe odpowiednich obrazów, z jednoznacznie selektywną rolą filtru ''idealnego''.

Przykładowe postacie charakterystyk częstotliwościowych filtrów stosowanych w przetwarzaniu obrazów są następujące:

  • filtry Butterwortha, o maksymalnie płaskiej charakterystyce amplitudowej w paśmie przepustowym, bez oscylacji (zafalowań) w pasmach przepustowym i zaporowym, o charakterystyce fazowej najbardziej zbliżonej do liniowej, jednak przy stosunkowo wolno opadającym widmie pasma przejściowego; 
    • uśredniający (dolnoprzepustowy) 
|H_B|(d)=\dfrac{1}{1+(d/d_0)^{2n}} 

(4.63) 

gdzie n \in \mathbb{Z} jest rzędem filtru -- wraz ze wzrostem n rośnie nachylenie charakterystyki w paśmie przejściowym, kosztem jednak coraz bardziej nieliniowej fazy; przy odpowiednio dużym n pojawiają się efekty Gibbsa; 

  • górnoprzepustowy  
|G_B(d)|=\dfrac{1}{1+(d_0/d)^{2n}} 

(4.64) 

  • filtry gaussowskie 
    • uśredniający 
|H_g(d)|=e^{-(-1/2(d/d_0)^2)} 

(4.65) 

  • górnoprzepustowy 
|G_g(d)|=1-e^{-1/2(d/d_0)^2} 

(4.66) 

  • filtry wykładnicze rzędu n (uogólnienie filtrów gaussowskich)
    • uśredniający 
|H_e(d)|=e^{-(d/d_0)^n} 

(4.67) 

  • górnoprzepustowy 
|G_e(d)|=1-e^{-(d/d_0)^n} 

(4.68) 

W praktyce przy przetwarzaniu obrazów stosowane są zwykle rzeczywitoliczbowe, przybliżające postacie transmitancji opisanych powyżej filtrów częstotliwościowych, tak że H(d)\approx |H_i(d)|. Zróżnicowanie charakterystyk amplitudowych tych filtrów ukazano na rys. 4.20, zaś efekty filtracji na rys. 4.21.

Rys. 4.20 Charakterystyki amplitudowe wybranych filtrów -- kolejno od lewej do prawej, z góry na dół pary filtrów dolno- i górnoprzepustowych: idealny, Butterwortha rzędu n=1, n=4, n=10, gaussowski oraz  wykładniczy rzędu n=1, n=4, n=10.


  
Rys. 4.21 Efekty filtracji w dziedzinie częstotliwościowej; w dwóch rzędach zamieszczono przekształcone obrazy wraz z ich widmami amplitudowymi (pod spodem), odpowiednio dla filtrów dolnoprzepustowych Butterwortha z n=1 i n=10, wykładniczego z n=1, Gaussa, a także górnoprzepustowych -- Butterwortha z n=1 oraz ''idealnego'' przy względnym rozmiarze d_0 równym 20\%.