Wavelet transformace v anal´ yze a ˇ c´ıslicov´ em zpracov´ an´ı biomedic´ınsk´ ych sign´ al˚ u a obraz˚ u

(1)

VYSOK ´ A ˇ SKOLA CHEMICKO - TECHNOLOGICK ´ A V PRAZE FAKULTA CHEMICKO - INˇ ZEN ´ YRSK ´ A

USTAV PO ˇ ´ CÍTA ˇ COV´ E A ˇ RÍDICÍ TECHNIKY

DIPLOMOV ´ A PR ´ ACE

Wavelet transformace v anal´ yze a ˇ c´ıslicov´ em zpracov´ an´ı biomedic´ınsk´ ych sign´ al˚ u a obraz˚ u

2003/2004

Vypracovala:

Andrea Gavlasov´a

Vedouc´ı pr´ace:

Prof. Ing. Aleˇs Proch´azka, CSc.

Praha Kvˇ eten 2004

(2)

Prohlaˇsuji, ˇze jsem diplomovou práci vypracovala samostatnˇe a pouˇzila jen pramen˚u, které uvád´ım v soupisu bibliografických citac´ı.

(3)

Summary

The goal of diploma thesis is in the study of wavelet transform use for analysis and digital signal and image processing. First part of this work includes the study of signal and image analysis using 1-D and 2-D Fourier transform and simple digital filtering methods with application of window functions to reduce disturbing signal and undesirable image components. Further stu- dies are devoted to possibilities of signal decomposition and reconstruction by wavelet transform using different wavelet functions, general global and local thresholding, modification of selected wavelet coefficients by different methods based upon threshold limits estimation to achieve mini- mal error value, and application of hard and soft thresholding. All these methods were used for modification of coefficients after image decomposition by a selected wavelet function and results for each function were discussed. Comparison of results includes also rejection of additive components by digital filtration and by application of wavelet transform. The main part of diploma thesis is devoted to application of these methods for rejection of image noise to achieve image components enhancement and to provide information for image components feature extraction and classification.

Key words - discrete Fourier transform, digital filtering, discrete wavelet transform, Wavelet decomposition and reconstruction, MR image, global and local thresholding, error value mini- malisation, biomedical image processing

(4)

Souhrn

Zamˇeˇren´ım této diplomové práce je vyuˇzit´ı wavelet transformace v analýze a ˇc´ıslicovém zpra- cován´ı signál˚u a obraz˚u. Souˇcást´ı je také analýza obraz˚u a signál˚u s vyuˇzit´ım 1-D a 2-D diskrétn´ı Fourierovy transformace a metody jednoduché ˇc´ıslicové filtrace s aplikac´ı okénkové funkce k potlaˇcován´ı ruˇsivých sloˇzek signálu a obrazu. Dále práce zahrnuje moˇznosti dekompozice a rekonstrukce pomoc´ı wavelet transformace s pouˇzit´ım r˚uzných wavelet funkc´ı, modifikace d´ılˇc´ıch koeficient˚u rozkladu pomoc´ı r˚uzných metod prahován´ı, k nimˇz patˇr´ı univerzáln´ı globáln´ı prahován´ı, univerzáln´ı lokáln´ı prahován´ı, prahován´ı na základˇe minimalizace chyby, a jejich vari- ant hard a soft volby prahu. Vˇsechny metody byly pouˇzity pro modifikaci koeficient˚u po rozkladu obrazu vˇsemi wavelet funkcemi a výsledky metod pro jednotlivé funkce srovnány. Srovnáno je také potlaˇcován´ı aditivn´ıch sloˇzek metodou ˇc´ıslicové filtrace a wavelet transformac´ı. Hlavn´ı ˇcást práce je tedy vˇenována aplikac´ım tˇechto metod pro potlaˇcován´ı ruˇsivých sloˇzek a snahou je dosaˇzen´ı zvýraznˇen´ı obrazu pro potˇreby lepˇs´ı klasifikace obrazových informac´ı.

Kl´ıˇcová slova - diskrétn´ı Fourierova transformace, ˇc´ıslicová filtrace, diskrétn´ı Wavelet transformace, Wavelet dekompozice a rekonstrukce, MR obraz, globáln´ı a lokáln´ı prahován´ı, minimalizace chyby

(5)

Podˇ ekov´ an´ı

Nejprve bych ráda podˇekovala Prof. Ing. Aleˇsovi Procházkovi, CSc., mému vedouc´ımu, za jeho veden´ı, cenné rady, pˇripom´ınky a podnˇety nejen pˇri výzkumu a absolvován´ı studia.

Dále bych velmi ráda podˇekovala mým rodiˇc˚um za jejich stálou podporu, trpˇelivost a lásku. Bez nich by tato práce nemohla nikdy vzniknout.

Nakonec bych chtˇela podˇekovat Vladim´ırovi za jeho trpˇelivost, l´asku a podporu.

Praha Andrea Gavlasov´a

10. kvˇetna 2004

(6)

Obsah

Abstract i

Souhrn ii

Podˇekov´an´ı iii

Obsah iv

Seznam obr´azk˚u v

Seznam tabulek vi

1 Uvod´ 1

2 Fourierova a wavelet transformace 3

2.1 Fourierova transformace . . . 3

2.2 Diskr´etn´ı Fourierova transformace . . . 4

2.3 Ok´enkov´a Fourierova transformace . . . 6

2.4 Dvourozmˇern´a Fourierova transformace . . . 7

2.5 Princip wavelet transformace . . . 8

2.5.1 Podobnosti Fourierovy a wavelet transformace . . . 8

2.5.2 Rozd´ıly mezi Fourierovou a wavelet transformac´ı . . . 9

(7)

2.7 Haarova wavelet funkce . . . 11

2.8 Dekompozice a rekonstrukce sign´alu . . . 11

2.9 Dekompozice a rekonstrukce obrazu . . . 14

2.10 Daubechies D4 wavelet funkce . . . 18

3 Poˇrizov´an´ı MR obraz˚u 20 3.1 Princip magnetick´e rezonance . . . 21

4 Potlaˇcován´ı ruˇsivých sloˇzek signál˚u a obraz˚u 23 4.1 Potlaˇcován´ı pomoc´ı Fourierovy transformace . . . 23

4.1.1 Výsledky ˇc´ıslicové filtrace simulovaných signál˚u a reálného MR obrazu . . 25

4.2 Potlaˇcov´an´ı pomoc´ı wavelet transformace . . . 27

4.2.1 Zpracov´an´ı sign´al˚u . . . 27

4.2.2 Zpracov´an´ı obraz˚u . . . 31

4.2.3 Testovac´ı matice . . . 32

4.2.4 Re´aln´y MR obraz . . . 38

5 Závˇer 43 Literatura 46 A Kódován´ı obraz˚u 47 A.1 Barva . . . 47

A.2 Reprezentace obrazu . . . 49

A.3 Barevn´a reprezentace obrazu . . . 50

A.3.1 Indexovan´y obraz . . . 50

A.3.2 Intenzitn´ı obraz . . . 51

A.3.3 Cernob´ıl´ˇ y obraz . . . 51

(8)

A.3.4 RGB obraz . . . 51

A.4 Vizualizace obrazu . . . 51

A.5 Tvorba videosekvenc´ı . . . 58

A.6 Montage . . . 58

A.7 Immovie . . . 59

A.8 Movie . . . 59

A.9 Movie2avi . . . 60

A.10 Getframe . . . 60

B Zdrojov´e k´ody Matlabu 61 B.1 Seznam soubor˚u . . . 61

B.2 V´ypis program˚u . . . 63

(9)

Seznam obr´ azk˚ u

2.1 Zobrazen´ı (a) posloupnosti x(n) = sin(2πf0n) pro n = 0, 1, . . . , N − 1 a detekce

frekvenˇcn´ıch sloˇzek z (b) frekvenˇcn´ıho spektra f0 = 0.2 Hz. . . 5

2.2 Zobrazen´ı (a) posloupnosti x(n), která vznikla souˇctem dvou sinusových funkc´ı s r˚uznými frekvencemi a (b) detekce frekvenˇcn´ıch sloˇzek f0 = 0.2 a f1= 0.4 Hz. 5 2.3 Zobrazen´ı (a) simulovaného signálu s lineárnˇe rostouc´ı frekvenc´ı a (b) STFT pˇri dˇelen´ı signálu na 12 výbˇerových oken. Velikost kaˇzdého okna je 16s. . . 7

2.4 Zobrazen´ı (a) posloupnosti 2D sign´alu a (b) detekce frekvenˇcn´ıch sloˇzek. . . 8

2.5 Porovnán´ı Fourierovy a wavelet analýzy, (a)-horˇs´ı ˇcasové rozliˇsen´ı, lepˇs´ı frekv. rozliˇsen´ı, (b)-lepˇs´ı ˇcasové rozliˇsen´ı, horˇs´ı frekv. rozliˇsen´ı, (c)-standardn´ı DWT, která má nejˇsirˇs´ı okénko pˇri n´ızkých frekvenc´ıch, kde se nacház´ı nejd˚uleˇzitˇejˇs´ı detaily. . . 9

2.6 Shannonova funkce a jej´ı spektr´aln´ı odhad. . . 10

2.7 Mallatovo sch´ema dekompozice a rekonstrukce sign´alu. . . 12

2.8 Zobrazen´ı p˚uvodn´ıho simulovan´eho sign´alu. . . 13

2.9 Perfektn´ı rekonstrukce p˚uvodn´ıho sign´alu. . . 14

2.10 Zobrazen´ı wavelet koeficient˚u a jejich poˇctu v z´avislosti na ´urovni dekompozice. . 14

2.11 Mallatovo schéma dekompozice a rekonstrukce obrazu a jednorozmˇerného pˇr´ıpadu signálu. . . 15

2.12 Zobrazen´ı p˚uvodn´ı testovac´ı matice a dekompozice do dvou matic D, A a jejich n´asledn´a dekompozice do ˇctyˇr matic AA, AD, DA, DD. . . 17

2.13 Rekonstrukce matice na základˇe zachovaných dekompoziˇcn´ıch koeficient˚u. Koneˇcný obraz je stejný jako p˚uvodn´ı . . . 17

(10)

2.14 Rekonstrukce matice z modifikovaných koeficient˚u. Koeficienty, které byly modi- fikovány jednoduˇse rozpoznáme z grafu wavelet koeficient˚u. . . 18

4.1 C´ıslicov´ˇ a filtrace simulovaného 1D signálu znázorˇnuj´ıc´ı (a) zaˇsumˇenou posloupnost, (b) filtrovanou posloupnost, (c) spektrum p˚uvodn´ıho signálu, (d) okénkovou funkci odstraˇnuj´ıc´ı vysoké frekvence a druhou neˇzádouc´ı frekvenci, (e) spektrum filtrované posloupnosti. . . 25 4.2 C´ıslicov´ˇ a filtrace simulovaného 2D signálu zobrazuje (a) signál s neˇzádouc´ı frek-

venc´ı, (b) signál s odfiltrovanou frekvenc´ı, (c) spektrum p˚uvodn´ıho signálu, (d) okénkovou funkci potlaˇcuj´ıc´ı druhou frekvenci v p˚uvodn´ım signálu, (e) spektrum filtrovaného signálu. . . 25 4.3 C´ıslicov´ˇ a filtrace reálného MR obrazu. Znázorˇnuje (a) p˚uvodn´ı MR obraz, (b)

MR obraz s pˇridaným ˇsumem, (c) filtrovaný MR obraz,(d) spektrum zaˇsumˇeného obrazu s vysokými frekvencemi, (e) 2D okénkovou funkci, která potlaˇcuje vysoké frekvence a zachovává pouze d˚uleˇzité n´ızko-frekvenˇcn´ı sloˇzky, (f ) spektrum fil- trovaného obrazu. . . 26 4.4 Filtrovaný MR obraz. . . 26 4.5 Zobrazen´ı p˚uvodn´ıch koeficient˚u (modˇre) a koeficient˚u po hard prahován´ı (ˇcervenˇe). 28 4.6 Zobrazen´ı p˚uvodn´ıch koeficient˚u (modˇre) a koeficient˚u po soft prahován´ı (ˇcervenˇe).

28

4.7 Zobrazen´ı wavelet koeficient˚u po dekompozici do páté úrovnˇe. Dlouhá okénka v 5.

´

urovni poskytuj´ı nejlepˇs´ı frekvenˇcn´ı rozliˇsen´ı. Úzká okénka v 1. úrovni maj´ı n´ızké frekvenˇcn´ı rozliˇsen´ı, ale dobré rozliˇsen´ı ˇcasové. . . 29 4.8 Analýza signálu zahrnuj´ıc´ı (a) p˚uvodn´ı signál, (b) rekonstruovaný signál, (c)

scalogram, (d) scaling a wavelet koeficienty, p˚uvodn´ı po dekompozici (fialová), modifikované po globáln´ım hard prahován´ı pro rekonstrukci (tyrkysová). . . 30 4.9 Zobrazen´ı p˚uvodn´ı testovac´ı matice. . . 34 4.10 Testovac´ı matice s pˇridaným ˇsumem. . . 34 4.11 Grafické srovnán´ı chyb RMSE pˇri hard modifikaci koeficient˚u uvedenými meto-

dami a pˇri wavelet transformaci do 1. stupnˇe. . . 34 4.12 Nejlepˇs´ı filtrovan´a matice pˇri transformaci funkc´ı Daubechies Db8 do 1.´urovnˇe a

aplikac´ı glob´aln´ı hard varianty metody minim´aln´ı chyby. . . 34

(11)

4.13 Grafické srovnán´ı chyb RMSE pˇri hard modifikaci koeficient˚u uvedenými meto-

dami a pˇri wavelet transformaci do 2. stupnˇe. . . 35

4.14 Nejlepˇs´ı filtrovaná matice pˇri transformaci funkc´ı Daubechies Db8 do 2.úrovnˇe a aplikac´ı globáln´ı hard varianty metody minimáln´ı chyby. . . 35

4.15 Grafické srovnán´ı chyb RMSE pˇri hard modifikaci koeficient˚u uvedenými metodami a pˇri wavelet transformaci do 3. stupnˇe. . . 35

4.16 Nejlepˇs´ı filtrovaná matice pˇri transformaci funkc´ı Daubechies Db8 do 3.úrovnˇe a aplikac´ı globáln´ı hard varianty metody minimáln´ı chyby. . . 35

4.17 Grafické srovnán´ı chyb RMSE pˇri soft modifikaci koeficient˚u uvedenými metodami a pˇri wavelet transformaci do 1. stupnˇe. . . 36

4.18 Nejlepˇs´ı filtrovaná matice pˇri transformaci funkc´ı Daubechies Db8 do 1.úrovnˇe a aplikac´ı globáln´ı soft varianty metody minimáln´ı chyby. . . 36

4.23 Nejlepˇs´ı v´ysledek wavelet transformace testovac´ı matice. . . 37

4.24 Zobrazen´ı p˚uvodn´ı re´aln´eho MR obrazu. . . 39

4.25 Reálný MR obraz s pˇridaným ˇsumem. . . 39

4.27 Nejlepˇs´ı filtrovaný MR obraz pˇri transformaci funkc´ı Daubechies Db8 do 1.úrovnˇe a aplikac´ı globáln´ı hard varianty metody minimáln´ı chyby. . . 39

(12)

4.30 Grafické srovnán´ı chyb RMSE pˇri hard modifikaci koeficient˚u uvedenými meto-

dami a pˇri wavelet transformaci do 3. stupnˇe. . . 40

4.33 Nejlepˇs´ı filtrovaný MR obraz pˇri transformaci funkc´ı Daubechies Db8 do 1.úrovnˇe a aplikac´ı globáln´ı soft varianty metody minimáln´ı chyby. . . 41

4.38 Nejlepˇs´ı v´ysledek wavelet transformace MR obrazu. . . 42

5.1 Zobrazen´ı výˇrez˚u reálného MR obrazu (a) po ˇc´ıslicové filtraci 2D obdéln´ıkovým okénkem, (b) po wavelet transformaci do 2. stupnˇe funkc´ı Daubechies Db4, (c) p˚uvodn´ıho (bez ˇsumu). . . 43

A.1 Rozdˇelen´ı elektromagnetick´eho vlnˇen´ı podle jeho frekvenc´ı [M.03a]. . . 47

A.2 Zobrazen´ı hladk´e kˇrivky na zaˇr´ızen´ı s dan´ym rozliˇsen´ım [M.03b]. . . 50

A.3 Matice obrazu a barevn´a paleta [M.03a]. . . 50

A.4 Obr´azky z´ıskan´e pomoc´ı pˇr´ıkazu imshow. . . 53

A.5 Intenzitn´ı obrazy z´ıskan´e pomoc´ı pˇr´ıkazu imshow. . . 53

A.6 Intenzitn´ı a ˇcernob´ıl´y obraz z´ıskan´y pomoc´ı pˇr´ıkazu imshow. . . 53

A.7 Colormap editor - prvku v matici s hodnotou 0,75 je pˇriˇrazena barva z colormapy s indexem 49. Barva je vyj´adˇrena v RGB nebo HSV hodnot´ach. . . 54

(13)

A.9 Colormap editor pro testovac´ı matici zobrazenou pomoc´ı image (pˇr. prvky matice s hodnotou 32 maj´ı svˇetle zelenou barvu, oznaˇcenou v colormapˇe indexem 33). . 55 A.10 Zobrazen´ı testovac´ı matice pomoc´ı imagesc. . . 56 A.11 Colormap editor pro testovac´ı matici zobrazenou pomoc´ı imagesc . . . 56 A.12 Zobrazen´ı testovac´ı matice pomoc´ı imagesc s pouˇzit´ım clims. . . 57 A.13 Colormap editor pro testovac´ı matici zobrazenou pomoc´ı imagesc z pouˇzit´ım clims . 57 A.14 Zobrazen´ı v´ıcesn´ımkového obrazového pole, které vzniklo pouˇzit´ım pˇr´ıkazu montage. 59 A.15 Zobrazen´ı jednotlivých sn´ımk˚u z´ıskaných pomoc´ı getframe. . . 60

(14)

Seznam tabulek

4.1 Chyba RMSE pˇri filtraci n´ızko-frekvenˇcn´ım filtrem. . . 26 4.2 Srovn´an´ı chyb RMSE pˇri hard a soft modifikaci uveden´ymi metodami a pˇri wavelet

transformaci do 1. stupnˇe vˇsemi wavelet funkcemi. . . 33 4.3 Srovn´an´ı chyb RMSE pˇri hard a soft modifikaci uveden´ymi metodami a pˇri wavelet

transformaci do 3. stupnˇe vˇsemi wavelet funkcemi. . . 38

(15)

Kapitola 1

Uvod ´

Oblast digitáln´ıho zpracováván´ı je oblast zpracován´ı obraz˚u a signál˚u. Odstranˇen´ı ˇsumu a následná rekonstrukce obraz˚u tvoˇr´ı základn´ı ˇcást metod zpracován´ı obraz˚u s mnoha aplikacemi vˇcetnˇe zpracován´ı inˇzenýrských obraz˚u a zlepˇsován´ı biomedic´ınských struktur. K nejvýnamnˇejˇs´ım metodám pro zpracován´ı obraz˚u patˇr´ı wavelet transformace. Ta m˚uˇze být pouˇzita pro dekompozici obrazu a odstranˇen´ı neˇzádouc´ıch sloˇzek. U znehodnocených signál˚u m˚uˇzeme pˇri pouˇzit´ı filtru zlepˇsit kvalitu signálu, z´ıskat d˚uleˇzité informace. S t´ım je pˇredevˇs´ım spojována diskrétn´ı Fourierova transformace. Obˇe transformace tvoˇr´ı ned´ılnou souˇcást oblasti DSP, která nahrazuje analogovou techniku.

Práce je rozdˇelena do 6 ˇcást´ı, které zahrnuj´ı

• matematický popis Fourierovy transformace, základn´ı frekvenˇcn´ı analýzu signálu pomoc´ı Fourierovy transformace, základy wavelet transformace, dekompozici signálu a obrazu pomoc´ı wavelet transformace a následnou rekonstrukci z hlediska ovˇeˇren´ı wavelet transformace a posouzen´ı vzájemné shody a rozd´ıly mezi transformacemi.

• techniku poˇrizován´ı biomedic´ınských obraz˚u metodou magnetické rezonance

• aplikaˇcn´ı ˇcást diskrétn´ı Fourierovy transformace - ˇc´ıslicová filtrace simulovaných dat a posléze MR obrazu a diskrétn´ı wavelet transformace - dekompozice, modifikace metodami prahován´ı a rekonstrukce k odstranˇen´ı ruˇsivých sloˇzek vˇcetnˇe ovˇeˇren´ı navrˇzených postup˚u na simulovaných datech

• uˇzit´ı výsledných algoritm˚u pro analýzu a zpracován´ı biomedic´ınských dat

• závˇer, ve kterém je vyhodnoceno srovnán´ı metod a obou transformac´ı

• kódován´ı obraz˚u a vizualizaci v pˇr´ıloze, kde jsou shrnuty základy reprezentace obraz˚u a tvorby videa

• zdrojov´e k´ody na CD

(16)

Vˇsechny programy jsem provedeny v Matlabu Version 6.5.0.180913a Release13, seznam zdro- jových kód˚u je uveden v pˇr´ıloze a programové kódy jsou uloˇzeny na CD. Souhrn a vybrané výsledky práce jsou rovnˇeˇz zveˇrejnˇeny na WWW stranˇe autorky (http://dsp.vscht.cz).

(17)

Kapitola 2

Fourierova a wavelet transformace

Pouˇz´ıvaj´ı se jako základn´ı matematické metody pro zpracován´ı signál˚u a obraz˚u.

2.1 Fourierova transformace

Fourierova transformace je názornou metodou teorie signál˚u, která reprezentuje signál pomoc´ı jeho frekvenˇcn´ı charakteristiky. Stejnou metodu lze také vyuˇz´ıt pˇri zpracován´ı obrazu, který je také moˇzno vyjádˇrit jako superpozici sinusových funkc´ı r˚uzných fáz´ı a amplitud. V pˇr´ıpadˇe zpra- cován´ı obrazu se nepracuje se spojitým, ale s diskrétn´ım obrazem a t´ım se tedy nˇekteré postupy zjednoduˇsuj´ı [KA94]. Pro reprezentaci funkce jej´ı frekvenˇcn´ı charakteristikou máme k dispo- zici dvˇe reprezentace. V prvn´ım pˇr´ıpadˇe, kdy pracujeme s obrazovými prvky, se pohybujeme v ˇ

casové oblasti, v druhém pˇr´ıpadˇe hovoˇr´ıme o oblasti frekvenˇcn´ı. Pˇrechod mezi tˇemito oblastmi nám umoˇzˇnuje Fourierova transformace.

Definice spojit´e Fourierovy transformace

Pˇr´ım´a transformace - pˇrevede sign´al do frekvenˇcn´ı oblasti

X(jf ) = Z ∞

−∞

x(t)e^{−j2πf t}dt (2.1)

Zpˇetná (inverzn´ı) transformace - pˇrevede signál z frekvenˇcn´ı do ˇcasové oblasti

x(t) = Z ∞

−∞

X(f )e^{j2πf t}df (2.2)

(18)

2.2 Diskr´ etn´ı Fourierova transformace

Pˇri zpracován´ı signál˚u pomoc´ı poˇc´ıtaˇcové techniky se pracuje s koneˇcným poˇctem vzork˚u a pro spojité funkce to tedy znamená, ˇze lze pracovat pouze s urˇcitým poˇctem vzork˚u tˇechto funkc´ı.

Pracuje se tedy s diskrétn´ımi pr˚ubˇehy i ve frekvenˇcn´ı oblasti. Signály v ˇcasové i frekvenˇcn´ı oblasti maj´ı koneˇcný poˇcet hodnot N a pˇri výpoˇctech se povaˇzuj´ı za periodické. Transformace, která umoˇzˇnuje pˇrechody mezi ˇcasovou oblast´ı, kde nezávisle promˇennou budeme znaˇcit n, a frek- venˇcn´ı oblast´ı, kde nezávisle promˇennou budeme znaˇcit k, je tzv. finitn´ı Fourierova transformace.

Nazývá se diskrétn´ı Fourierova transformace(DFT) [KA94]. Analyzuje signál nalezen´ım jeho am- plitudy a fázového spektra. Napˇr. DFT jedné ˇcisté sinusovky nalezne jeden p´ık ve frekvenˇcn´ım spektru a tˇri p´ıky pro signál sloˇzený ze tˇr´ı sinusovek.

Pro výpoˇcet této transformace byly vypracovány algoritmy, kterým se ˇr´ıká rychlá Fourierova transformace(Fast Fourier Transform - FFT). T´ımto se DFT stala výkonným nástrojem pro frekvenˇcn´ı analýzu a ˇc´ıslicovou filtraci.

Definice DFT

Pˇr´ım´a transformace - pˇrevede sign´al do frekvenˇcn´ı oblasti

X(k) =

N −1

X

n=0

x(n)e^−j^2π^N^kn, k = 0, 1, 2, . . . , N − 1 (2.3)

Zpˇetná (inverzn´ı) transformace - pˇrevede signál z frekvenˇcn´ı do ˇcasové oblasti

x(n) = 1 N

N −1

X

k=0

X(k)e^j^2π^N^kn, n = 0, 1, 2, . . . , N − 1 (2.4)

DFT umoˇzˇnuje detekovat frekvenˇcn´ı sloˇzky sign´alu, coˇz je zn´azornˇeno na obr. 2.1, 2.2.

(19)

0 10 20 30 40 50 60 70

−1

−0.5 0 0.5 1

(a) Prubeh posloupnosti x(n) s frekvenci f 0

n

x(n)

0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1

0 5 10 15 20 25 30

(b) Frekvencni spektrum

fk

X(k)

Obr´azek 2.1: Zobrazen´ı (a) posloupnosti x(n) = sin(2πf0n) pro n = 0, 1, . . . , N − 1 a detekce frekvenˇcn´ıch sloˇzek z (b) frekvenˇcn´ıho spektra f₀ = 0.2 Hz.

0 10 20 30 40 50 60 70

−2

−1 0 1 2

(a) Prubeh dvou posloupnosti x(n) s frekvencemi f 0 a f

1

n

x(n)

0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1

0 5 10 15 20 25 30

(b) Frekvencni spektrum

fk

X(k)

Obrázek 2.2: Zobrazen´ı (a) posloupnosti x(n), která vznikla souˇctem dvou sinusových funkc´ı s r˚uznými frekvencemi a (b) detekce frekvenˇcn´ıch sloˇzek f₀= 0.2 a f₁= 0.4 Hz.

(20)

2.3 Ok´ enkov´ a Fourierova transformace

V roce 1946 upravil Denis Gabor FT pro vyhodnocován´ı malých ˇcást´ı signál˚u najednou a tato technika dostala název zpracován´ı signálu v okénc´ıch. Gaborova úprava nazvaná krátkodobá Fourierova tranformace (STFT) zakresluje signál jako dvojrozmˇernou funkci ˇcasu a frekvence.

STFT reprezentuje jakýsi druh kompromisu mezi ˇcasovým a frekvenˇcn´ım pohledem na signál.

Poskytuje informaci o oboj´ım, kdy a pˇri jaké frekvenci nastane zmˇena signálu. Avˇsak m˚uˇzeme z´ıskat pouze informaci s omezeným stupnˇem pˇresnosti, která je dána velikost´ı okénka. Zat´ımco tento kompromis m˚uˇze být uˇziteˇcný, stinnou stránkou je, ˇze si vybereme jednu velikost okénka, která je pak stejná pro vˇsechny frekvence. Hodnˇe signál˚u vˇsak vyˇzaduje flexibilnˇejˇs´ı pˇr´ıstup, který umoˇzˇnuje mˇenit velikost okénka pro pˇresnˇejˇs´ı urˇcen´ı ˇcasu nebo frekvence [TM97].

Základn´ım principem je rozdˇelen´ı signálu na malé ˇcásti, u nichˇz m˚uˇzeme pˇredpokládat do- stateˇcnou stacionaritu. To provedeme multiplikac´ı okénkové funkce a signálu. Okénko se po- souvá v ˇcase a kaˇzdý výˇrez signálu okénkem je uloˇzen do jednoho sloupce nové matice. Na kaˇzdý sloupec pak aplikujeme FFT a vzniklou ˇcasovˇe-frekvenˇcn´ı závislost zobraz´ıme.

Definice spojit´e STFT

ST F T (t, ω) = Z ∞

−∞

x(τ )w^∗(τ − t)e^−jωτdτ (2.5)

w . . . ok´enkov´a funkce

∗ . . . komplexn´ı konjunkce t . . . ˇcasové posunut´ı okénka x(τ ) . . . signál

Definice diskr´etn´ı STFT

ST F T (n, ω) =

∞

X

m=−∞

x(n + m)w(m)e^−jωm (2.6)

Velikost okének urˇcuje, zda bude signál v kaˇzdém oknˇe stacionárn´ı. Napˇr. nˇekolik krátkých okének zvyˇsuje pravdˇepodobnost, ˇze bude signál v kaˇzdém oknˇe stacionárn´ı. Dokud kratˇs´ı okénka mohou zlepˇsovat schopnost DFT analyzovat nestacionárn´ı signály, tak jsou to dobré výsledky.

Ale pokud nezmˇen´ıme rozliˇsen´ı, tak bude v kratˇs´ım okénku ménˇe bod˚u, které znamenaj´ı ménˇe detail˚u ve spektru signálu. Kratˇs´ı ˇcasové okénko poskytuje dobré ˇcasové rozliˇsen´ı, ale horˇs´ı frekvenˇcn´ı rozliˇsen´ı, protoˇze je ˇcas sn´ımán´ı signálu pˇr´ıliˇs krátký.

Dobré ˇcasové a frekvenˇcn´ı rozliˇsen´ı nen´ı moˇzné dosáhnout ve stejný ˇcas. Je obt´ıˇzné vybrat pˇrimˇeˇrenou délku okna a vhodnou pozici hranice okna [dV02].

(21)

Casov´ˇ e a frekvenˇcn´ı rozliˇsen´ı v z´avislosti na d´elce okna

—————————————————————————————–

Krátké okno Dobré ˇcasové rozliˇsen´ı Horˇs´ı frekvenˇcn´ı rozliˇsen´ı Sirok´ˇ e okno Horˇs´ı ˇcasové rozliˇsen´ı Dobré frekvenˇcn´ı rozliˇsen´ı

—————————————————————————————–

Obrázek 2.3: Zobrazen´ı (a) simulovaného signálu s lineárnˇe rostouc´ı frekvenc´ı a (b) STFT pˇri dˇelen´ı signálu na 12 výbˇerových oken. Velikost kaˇzdého okna je 16s.

2.4 Dvourozmˇ ern´ a Fourierova transformace

2D Fourierova transformace m´a analogickou vˇetˇsinu vlastnost´ı s jednorozmˇernou Fourierovou transformac´ı.

Pˇr´ım´a 2D DFT

X(k, l) =

M −1

X

m=0 N −1

X

n=0

x(m, n)e^−j2π(^km^M⁺^ln^N⁾ (2.7)

(22)

Inverzn´ı 2D DFT

x(m, n) = 1 M N

M −1

X

k=0 N −1

X

l=0

X(k, l)e^j2π(^km^M ⁺^ln^N⁾ (2.8)

Obr´azek 2.4: Zobrazen´ı (a) posloupnosti 2D sign´alu a (b) detekce frekvenˇcn´ıch sloˇzek.

2.5 Princip wavelet transformace

Wavelet transformace pˇredstavuje alternativu ke krátkodobé Fourierovˇe transformaci a je ˇcást´ı hlavn´ıho smˇeru DSP. Wavelet analýza pˇredstavuje dalˇs´ı logický krok: okénkovou techniku

s promˇenlivou velikost´ı okénka. Základem transformace je ˇcasovˇe omezená funkce h(t).

2.5.1 Podobnosti Fourierovy a wavelet transformace

Rychlá Fourierova transformace (FFT) a diskrétn´ı wavelet transformace (DWT) jsou lineárn´ı operace, které vytváˇr´ı datovou strukturu, která obsahuje log2n ˇcást´ı r˚uzných délek obvykle doplnˇené a transformované do r˚uzných vektor˚u délky 2ⁿ. Matematické vlastnosti matic, které se pouˇz´ıvaj´ı pˇri transformac´ıch jsou rovnˇeˇz podobné. Inverzn´ı matice jak pro FFT, tak i pro DWT je odvozena ze základn´ı matice. Obˇe transformace pracuj´ı s okénky.

(23)

2.5.2 Rozd´ıly mezi Fourierovou a wavelet transformac´ı

Diskrétn´ı Fourierova transformace i wavelet transformace poskytuj´ı informace o frekvenˇcn´ım obsahu signálu. Diskrétn´ı Fourierovy transformace mohou poskytnout pˇresnou frekvenˇcn´ı lokalizaci, ale ne pˇresnou ˇcasovou lokalizaci uvnitˇr okénka. Jak jiˇz bylo ˇreˇceno, ˇcasové rozliˇsen´ı se zlepˇsuje, kdyˇz se zuˇzuje okénko.

Napˇr. jeden signál s 32 vzorky, m˚uˇze být rozdˇelen pomoc´ı okénkové funkce na osm signál˚u po ˇctyˇrech vzorc´ıch nebo na ˇctyˇri signály po osmi vzorc´ıch. Lepˇs´ı ˇcasové rozliˇsen´ı znamená horˇs´ı frekvenˇcn´ı rozliˇsen´ı. At’ jsou frekvence rozdˇeleny jakkoliv, rozliˇsen´ı je stejné pro vˇsechna okénka.

Zvládnut´ı frekvenˇcn´ıho rozliˇsen´ı je významným pˇr´ınosem diskrétn´ı wavelet transformace (DWT).

Pˇri vysokých frekvenc´ıch poskytuje DWT horˇs´ı frekvenˇcn´ı rozliˇsen´ı, ale dobré ˇcasové rozliˇsen´ı.

Pˇri n´ızkých frekvenc´ıch je horˇs´ı ˇcasové rozliˇsen´ı, ale lepˇs´ı frekvenˇcn´ı rozliˇsen´ı. DWT pˇresnˇe urˇcuje frekvence z n´ızkofrekvenˇcn´ıch signál˚u pouˇzit´ım ˇsirokého okénka a pˇresnˇe urˇcuje ˇcas z vy- sokofrekvenˇcn´ıch signál˚u pouˇzit´ım úzkých okének. Na druhé stranˇe, ˇcasové rozliˇsen´ı je horˇs´ı právˇe pˇri niˇzˇs´ıch frekvenc´ıch protoˇze okénka jsou ˇsiroká a frekvenˇcn´ı rozliˇsen´ı je horˇs´ı pˇri vy- sokých frekvenc´ıch, protoˇze okénka jsou pˇr´ıliˇs úzká k pˇresné identifikaci frekvenc´ı. DWT pokryt´ı ˇ

casovˇe-mˇeˇr´ıtkového pole je lepˇs´ı neˇz pokryt´ı pomoc´ı DFT. Mˇeˇr´ıtko je nepˇr´ımo úmˇerné frekvenci. Kromˇe toho jsou wavelet transformace pˇrizp˚usobivé a umoˇzˇnuj´ı rovnˇeˇz jiné moˇznosti.

Standardn´ı DWT má nejuˇzˇs´ı filtry pˇri n´ızkých frekvenc´ıch, kde se obvykle nacház´ı nejd˚uleˇzitˇejˇs´ı detaily signálu. Jiná DWT m˚uˇze um´ıstit nejuˇzˇs´ı filtry do jiných oblast´ı frekvenˇcn´ıho spektra, vˇsude, kde poˇzadujeme v´ıce podrobnost´ı o signálu [dV02]. Výhodou wavelet transformace tedy

0 4 8

0.5 0.4375 0.375 0.3125 0.25 0.1875 0.125 0.0625 0

cas (s)

frekvence

(a) STFT ANALYZA 2/8

0 2 4 6 8

cas (s) (b) STFT ANALYZA 4/4

0 1 2 3 4 5 6 7 8

1 2 3 4

cas (s)

meritko

(c) WAVELET DEKOMPOZICE

Obrázek 2.5: Porovnán´ı Fourierovy a wavelet analýzy, (a)-horˇs´ı ˇcasové rozliˇsen´ı, lepˇs´ı frekv.

rozliˇsen´ı, (b)-lepˇs´ı ˇcasové rozliˇsen´ı, horˇs´ı frekv. rozliˇsen´ı, (c)-standardn´ı DWT, která má nejˇsirˇs´ı okénko pˇri n´ızkých frekvenc´ıch, kde se nacház´ı nejd˚uleˇzitˇejˇs´ı detaily.

je, ˇze se okénka mˇen´ı. Naproti tomu diskrétn´ı Fourierova transformace pouˇz´ıvá goniometrické funkce sinus a kosinus stejné délky. Zapamatujme si, ˇze wavelet transformace nemá jednotný soubor základn´ıch funkc´ı jako Fourierova transformace, která pouˇz´ıvá jen funkce sinus a kosinus.

(24)

2.6 Analytick´ a wavelet funkce

Wavelet funkci nazýváme ”vlnkovou” funkc´ı a to na základˇe podobnosti pr˚ubˇehu tˇechto funkc´ı s ”vlnkou” definovanou na daném intervalu hodnot nezávisle promˇenné. Jej´ı dilatace umoˇzˇnuje analýzu daného signálu s r˚uzným rozliˇsen´ım v návaznosti na kompresi spektra této funkce ve frekvenˇcn´ı oblasti.

Pˇr´ıkladem analyticky definovan´e wavelet funkce je napˇr. Shannonova wavelet funkce h(t) = sin(π ·₂^t)

π ·₂^t · cos(3 · π · t

2) (2.9)

Souvislost mezi dilatac´ı t´eto funkce a odhadem spektra je uvedena na obr. 2.6.

0 100 200 300 400

−1 0 1

a = 1

Shannonova f. h(t)

0.1 0.2 0.3 0.4 0.5

0.005 0.01 0.015 0.02 0.025

spektrum Shannonovy f.

0 100 200 300 400

−1 0 1

a = 2

komprese h(t)

0.1 0.2 0.3 0.4 0.5

0.005 0.01 0.015 0.02 0.025

dilatace spektra

0 100 200 300 400

−1 0 1

a = 0.75

dilatace h(t)

0.1 0.2 0.3 0.4 0.5

0.005 0.01 0.015 0.02 0.025

komprese spektra

0 100 200 300 400

−1 0 1

a = 0.5

dilatace h(t)

t

h(t)

0.1 0.2 0.3 0.4 0.5

0.005 0.01 0.015 0.02 0.025

komprese spektra

f

H(t)

Wavelet funkce Spektrální odhad

Obr´azek 2.6: Shannonova funkce a jej´ı spektr´aln´ı odhad.

Délka wavelet funkce nám umoˇzˇnuje volit interval dané ˇcasové ˇrady a t´ım i lepˇs´ı urˇcen´ı oblasti zmˇen v dané ˇradˇe. Výbˇer tvaru wavelet funkce v daném intervalu, umoˇzˇnuje volit typ deteko- vané nestacionarity. D´ıky tˇemto vlastnostem z´ıskáváme rozsáhlé moˇznosti analýzy a zpracován´ı daného signálu s moˇznost´ı volby délky ˇcasového intervalu, ve kterém se nestacionarita signálu vyskytla.

(25)

2.7 Haarova wavelet funkce

Z´akladn´ı Wavelet funkc´ı je Haarova funkce.

Máme dánu posloupnost {s_n}^N_n=1 = {s₁, s₂, . . . , s_N} . Pak pˇri kaˇzdé úrovni rozkladu zadané posloupnosti pomoc´ı Haarovy funkce z´ıskáme N/2 koeficient˚u v obou vˇetv´ıch dekompoziˇcn´ıho schematu, pˇriˇcemˇz

• prvn´ı vˇetev slouˇz´ı k detekci rychle promˇenných sloˇzek signálu a pˇr´ısluˇsné wavelet koeficienty se urˇcuj´ı ze vztahu [L.01]

dn= s_n− s_n+1

2 , n = 1, 2, . . . , N (2.10)

• druhá vˇetev slouˇz´ı k detekci pomalu promˇenných sloˇzek signálu a pˇr´ısluˇsné koeficienty (”scaling coefficients”) se urˇcuj´ı pomoc´ı vztahu [L.01]

a_n= s_n+ s_n+1

2 , n = 1, 2, . . . , N (2.11)

V pˇr´ıpadˇe diskrétn´ı wavelet transformace je moˇzné podvzorkován´ı dvˇema (downsampling), tzn.

ˇ

ze p˚uvodn´ı signál s_ndélky N je opˇet reprezentován N koeficienty. Z´ıskané koeficienty a_nse stávaj´ı vstupem do dalˇs´ıho kroku wavelet transformace. V pˇr´ıpadˇe, ˇze N = 2^L je tedy moˇzný rozklad do L úrovn´ı. P˚uvodn´ı posloupnost N prvk˚u je nahrazena vypoˇctenou ˇradou koeficient˚u v tomto poˇrad´ı: posledn´ı vektor a_L, následován ˇradou wavelet koeficient˚u jejichˇz velikost vektor˚u je vzr˚ustaj´ıc´ı mocninou dvou (2⁰, 2¹, 2², . . . , N/2).

Ve wavelet terminologii z mˇeˇr´ıtkové funkce (”scaling function”) vypoˇc´ıtáme koeficienty a_na z wavelet funkce vypoˇc´ıtáme wavelet koeficienty dn.

Rovnice (2.11) - mˇeˇr´ıtková funkce odpov´ıdá n´ızko-frekvenˇcn´ımu filtru (LP-filter) a rovnice (2.10) - wavelet funkce odpov´ıdá vysoko-frekvenˇcn´ımu filtru (HP-filter).

2.8 Dekompozice a rekonstrukce sign´ alu

V blokové podobˇe je rozklad daného signálu na dvˇe ˇcásti a jeho následná rekonstrukce uvedena na obr. 2.7 , který pˇredstavuje Mallatovo dekompoziˇcn´ı schema. Rozklad lze pˇritom realizovat pomoc´ı rovnic (2.10) a (2.11). Výstupn´ı signál z(n) = s(n), pokud jsou koeficienty rozkladu nezmˇenˇené.

(26)

↓

↑

↑ FAZE DEKOMPOZICE FAZE REKONSTRUKCE

Konvoluce signálu

Konvoluce signálu a

vzorkování dolu

a

vzorkování nahoru

D R

Puvodní signál

[s(n)]

Konecný signál

[z(n)]

d1(n)

a1(n)

d2(n)

a2(n)

d1(n)

a1(n) ↓

↓

↑

Konvoluce signálu

vzorkování dolu

a

vzorkování nahoru

D R

Puvodní signál

[s(n)]

Konecný signál

[z(n)]

d1(n)

a1(n)

d2(n)

a2(n)

d1(n)

a1(n) ↓

↓

↑

Konvoluce signálu

vzorkování dolu

a

vzorkování nahoru

D R

Puvodní signál

[s(n)]

Konecný signál

[z(n)]

d1(n)

a1(n)

d2(n)

a2(n)

d1(n)

a1(n) ↓

↓

↑

Konvoluce signálu

vzorkování dolu

a

vzorkování nahoru

D R

Puvodní signál

[s(n)]

Konecný signál

[z(n)]

d1(n)

a1(n)

d2(n)

a2(n)

d1(n)

a1(n) ↓

↓

↑

Konvoluce signálu

vzorkování dolu

a

vzorkování nahoru

D R

Puvodní signál

[s(n)]

Konecný signál

[z(n)]

d1(n)

a1(n)

d2(n)

a2(n)

d1(n)

a1(n) ↓

↓

↑

Konvoluce signálu

vzorkování dolu

a

vzorkování nahoru

D R

Puvodní signál

[s(n)]

Konecný signál

[z(n)]

d1(n)

a1(n)

d2(n)

a2(n)

d1(n)

a1(n)

Obrázek 2.7: Mallatovo schéma dekompozice a rekonstrukce signálu.

Pro následnou rekonstrukci po dekompozici do prvé úrovnˇe lze z tˇechto rovnic z´ıskat p˚uvodn´ı signál [L.01] pomoc´ı vztah˚u:

s_n= a_n+ d_n (2.12)

s_n+1= a_n− d_n (2.13)

Zavedeme-li koeficienty mˇeˇr´ıtkov´e funkce [h₀, h₁] = [1, 1] a koeficienty wavelet funkce [g0, g1] = [1, −1] lze vztahy (2.12) a (2.13) vyj´adˇrit ve tvaru

sn= an· h₀+ dn· h₁, n = 1, 3, 5, 7 (2.14) sn+1= an· g₀+ dn· g₁, n = 1, 3, 5, 7 (2.15)

Ten lze formulovat pomoc´ı konvolutorn´ıho souˇcinu

s_n=

L−1

X

k=0

h_k· v_n+k, n = 1, 3, 5, 7 (2.16)

sn+1=

L−1

X

k=0

g_k· v_n+k, n = 1, 3, 5, 7 (2.17)

kde L znaˇc´ı d´elku filtru (v naˇsem pˇr´ıpadˇe L = 2) pro vektor v ve tvaru

v = [a1, d1, a3, d3, a5, d5, a7, d7]⁰ (2.18)

(27)

V maticov´e formˇe lze tuto operaci popsat ve tvaru

s = H · v (2.19)

kde

s = [s1, s2, . . . , sN]⁰ (2.20)

a H je konvolutorn´ı matice

H =







h₀ h₁ 0 0 0 0 0 0

g₀ g₁ 0 0 0 0 0 0

0 0 h0 h1 0 0 0 0

0 0 g₀ g₁ 0 0 0 0

0 0 0 0 h0 h1 0 0

0 0 0 0 g₀ g₁ 0 0

0 0 0 0 0 0 h0 h1

0 0 0 0 0 0 g0 g1







V pˇr´ıpadˇe modifikace koeficient˚u rozkladu dané funkce lze potlaˇcit ˇcásti signálu (ˇc´ıslicová filtrace). Na obr. 2.8 a 2.9 je uveden pˇr´ıklad dekompozice a rekonstrukce dané posloupnosti.

1 2 3 4 5 6 7 8

−1

−0.8

−0.6

−0.4

−0.2 0 0.2 0.4 0.6 0.8 1

N

f(x)

Pùvodní signál

Obrázek 2.8: Zobrazen´ı p˚uvodn´ıho simulovaného signálu.

(28)

1 2 3 4 5 6 7 8

−1

−0.8

−0.6

−0.4

−0.2 0 0.2 0.4 0.6 0.8 1

N

s(N)

Signál po Wavelet rekonstrukci

Obr´azek 2.9: Perfektn´ı rekonstrukce p˚uvodn´ıho sign´alu.

1 2 3 4 5 6 7 8

−1

−0.5 0 0.5 1

Signál

1 2 3 4 5 6 7 8

−1

−0.5 0 0.5 1

WT koeficienty

0.5 1 1.5 2 2.5 3 3.5 4 4.5

3

2

1

20 40 60

Obrázek 2.10: Zobrazen´ı wavelet koeficient˚u a jejich poˇctu v závislosti na úrovni dekompozice.

2.9 Dekompozice a rekonstrukce obrazu

Základn´ı princip dekompozice a rekonstrukce signálu pomoc´ı wavelet transformace jsme si ukázali na obr. 2.7.

Na obr. 2.11 si ukáˇzeme princip dekompozice a rekonstrukce obrazu. Stejný pˇr´ıstup, který je

(29)

↓

↑

Sloupcova konvoluce

Radkova konvoluce

Sloupcova konvoluce a

podvzorkovani

a podvzorkovani

a prevzorkovani

D.1 D.2 R.1 R.2

D

A

DD

DA

AD

AA

D

A Puvodni

signal nebo obraz [X(n,m)]

Konecny signal nebo obraz [Z(n,m)]

Obrázek 2.11: Mallatovo schéma dekompozice a rekonstrukce obrazu a jednorozmˇerného pˇr´ıpadu signálu.

Výsledná matice [Z(n, m)] = [X(n, m)], pokud z˚ustanou matice rozkladu DD, DA, AD, AA pro rekonstrukci nezmˇenˇené.

DEKOMPOZICE

V prvn´ım stupni dekompozice (L = 1) docház´ı nejprve ke zpracován´ı a podvzorkován´ı matice X o rozmˇerech [M/2^L−1, N/2^L−1] po sloupc´ıch, kde

M . . .poˇcet ˇr´adk˚u matice N . . .poˇcet sloupc˚u matice

L . . .aktu´aln´ı ´uroveˇn (stupeˇn) dekompozice pˇriˇcemˇz

• v horn´ı vˇetvi aplikujeme vysoko-frekvenˇcn´ı filtr a z´ıskáme matici s rozmˇery [M/2^L, N/2^L−1] tzn. (sniˇzujeme poˇcet ˇrádk˚u na polovinu, poˇcet sloupc˚u z˚ustává), která obsahuje rychle promˇenné sloˇzky

D((i + 1)/2, j) = (X(i, j) − X(i + 1, j))/2, i = 1, 3, . . . , M j = 1, 2, . . . , N (2.21)

• ve spodn´ı vˇetvi aplikujeme n´ızko-frekvenˇcn´ı filtr a z´ıskáme matici s rozmˇery [M/2^L, N/2^L−1] tzn. (sniˇzujeme poˇcet ˇrádk˚u na polovinu, poˇcet sloupc˚u z˚ustává), která zachovává d˚uleˇzité pomalu promˇenné sloˇzky

A((i + 1)/2, j) = (X(i, j) + X(i + 1, j))/2, i = 1, 3, . . . , M j = 1, 2, . . . , N (2.22)

(30)

Z´ıskané matice A, D dále podvzorkováváme po ˇrádc´ıch jak je naznaˇceno v Mallatovˇe schematu.

• na matici D, která vznikla v horn´ı vˇetvi aplikujeme jak vysoko-frekvenˇcn´ı, tak n´ızko- frekvenˇcn´ı filtr a z této matice z´ıskáme dvˇe matice o rozmˇerech [M/2^L, N/2^L] tzn. (sniˇzujeme poˇcet sloupc˚u na polovinu, poloviˇcn´ı poˇcet ˇrádk˚u z˚ustává)

DD(i, (j + 1)/2) = (D(i, j) − D(i, j + 1))/2, j = 1, 3, . . . , N i = 1, 2, . . . , M/2^L (2.23) DA(i, (j+1)/2) = (D(i, j)+D(i, j+1))/2, j = 1, 3, . . . , N i = 1, 2, . . . , M/2^L (2.24)

• na matici A, která vznikla ve spodn´ı vˇetvi aplikujeme také vysoko-frekvenˇcn´ı a n´ızko- frekvenˇcn´ı filtr a z této matice z´ıskáme opˇet dvˇe matice o rozmˇerech [M/2^L, N/2^L] tzn.

(sniˇzujeme poˇcet sloupc˚u na polovinu, poloviˇcn´ı poˇcet ˇrádk˚u z˚ustává)

AD(i, (j +1)/2) = (A(i, j)−A(i, j +1))/2, j = 1, 3, . . . , N i = 1, 2, . . . , M/2^L (2.25) AA(i, (j +1)/2) = (A(i, j)+A(i, j +1))/2, j = 1, 3, . . . , N i = 1, 2, . . . , M/2^L (2.26)

Pˇri dalˇs´ım stupni dekompozice (L = 2) se celý postup opakuje a vstupn´ı matic´ı pro dekompozici se stává matice AA, která vznikla z p˚uvodn´ı matice X podvzorkován´ım n´ızko-frekvenˇcn´ımi filtry po sloupc´ıch a ˇrádc´ıch.

REKONSTRUKCE

Pro následnou rekonstrukci m˚uˇzeme pouˇz´ıt p˚uvodn´ı matice DD, DA, AD, AA a z´ıskáme p˚uvodn´ı obrazovou matici X o rozmˇerech [M, N ] nebo, pokud modifikujeme matice DD, DA, AD z´ıskáme modifikovanou obrazovou matici Xm s rozmˇery [M, N ].

• ˇrádková rekonstrukce matice D o rozmˇeru [M/2^L, N/2^L−1] tzn. (zvyˇsujeme poˇcet sloupc˚u o polovinu, poloviˇcn´ı poˇcet ˇrádk˚u z˚ustává) z matic DD a DA

D(i, 2 ∗ j − 1) = DA(i, j) + DD(i, j), i = 1, 2, . . . , M/2^L, j = 1, 2, . . . , N/2^L (2.27) D(i, 2 ∗ j) = DA(i, j) − DD(i, j), i = 1, 2, . . . , M/2^L, j = 1, 2, . . . , N/2^L (2.28)

• ˇrádková rekonstrukce matice A o rozmˇeru [M/2^L, N/2^L−1] tzn. (zvyˇsujeme poˇcet sloupc˚u o polovinu, poloviˇcn´ı poˇcet ˇrádk˚u z˚ustává) z matic AD a AA

A(i, 2 ∗ j − 1) = AA(i, j) + AD(i, j), i = 1, 2, . . . , M/2^L, j = 1, 2, . . . , N/2^L (2.29) A(i, 2 ∗ j) = AA(i, j) − AD(i, j), i = 1, 2, . . . , M/2^L, j = 1, 2, . . . , N/2^L (2.30)

• sloupcová rekonstrukce matice X o rozmˇeru [M/2^L−1, N/2^L−1] tzn. (zvyˇsujeme poˇcet ˇrádk˚u o polovinu, poˇcet sloupc˚u z˚ustává) z matic A a D

X(2 ∗ i − 1, j) = A(i, j) + D(i, j), i = 1, 2, . . . , M/2^L, j = 1, 2, . . . , N (2.31)

(31)

Na následuj´ıc´ıch obrázc´ıch si ukáˇzeme dekompozici testovac´ı matice a jej´ı následnou rekonstrukci s pouˇzit´ım p˚uvodn´ıch a modifikovaných matic DD, DA, AD. Tyto matice AA, AD, DA a DD pak zjednoduˇsenˇe pˇrevád´ıme do vektoru a zobrazujeme je jako wavelet koeficienty. Z grafu pak m˚uˇzeme jednoduˇse rozpoznat, které koeficienty jsme modifikovaly a které jsme zachovaly.

PUVODNI OBRAZ

1. Dve matice [M/2,N] 2. Ctyri matice [M/2,N/2]

Obrázek 2.12: Zobrazen´ı p˚uvodn´ı testovac´ı matice a dekompozice do dvou matic D, A a jejich následná dekompozice do ˇctyˇr matic AA, AD, DA, DD.

50 100 150 200 250

−1

−0.5 0 0.5 1

WAVELET KOEFICIENTY

KONECNY OBRAZ−CHYBA: 4.1885e−017

Obrázek 2.13: Rekonstrukce matice na základˇe zachovaných dekompoziˇcn´ıch koeficient˚u.

Koneˇcn´y obraz je stejn´y jako p˚uvodn´ı

(32)

50 100 150 200 250

−0.4

−0.2 0 0.2 0.4

WAVELET KOEFICIENTY

KONECNY OBRAZ−CHYBA: 0.44527

Obrázek 2.14: Rekonstrukce matice z modifikovaných koeficient˚u. Koeficienty, které byly modi- fikovány jednoduˇse rozpoznáme z grafu wavelet koeficient˚u.

2.10 Daubechies D4 wavelet funkce

Pro srovnán´ı s Haarovou wavelet funkc´ı se zde jen krátce zm´ın´ım o dalˇs´ı z rozsáhlé skupiny wavelet funkc´ı a to o Daubechies D4. Jak uˇz název napov´ıdá jsou scaling a wavelet koeficienty z Daubechies D4 wavelet transformace poˇc´ıtány s pouˇzit´ım ˇctyˇr koeficient˚u h0, h1, h2, h3 a g₀, g₁, g₂, g₃. Daubechies mˇeˇr´ıtková a wavelet funkce posouvá koeficienty o dvˇe m´ısta doleva pˇri kaˇzdém kroku wavelet transformace.

Jestliˇze porovn´ame Haarovu a Daubechies D4 konvolutorn´ı matici, tak u Haarovy matice nen´ı ˇ

zádný pˇrekryv mezi pˇredchoz´ımi a následnými páry wavelet a scaling koeficient˚u, jak je tomu u Daubechies D4.

Daubechies matice





h0 h1 h2 h3 0 0 . . . g₀ g₁ g₂ g₃ 0 0 . . . 0 0 h0 h1 h2 h3 . . .





(33)

Haarova matice







h₀ h₁ 0 0 0 0 . . . g0 g1 0 0 0 0 . . . 0 0 h₀ h₁ 0 0 . . . 0 0 g0 g1 0 0 . . . 0 0 0 0 h0 h1 . . . 0 0 0 0 g₀ g₁ . . . ... ... ... ... ... ... ...







Po vysoko-frekvenˇcn´ı filtraci Haarovou funkc´ı dostaneme výsledek, který odráˇz´ı rozd´ıl mezi sudými a lichými prvky. Rozd´ıl mezi lichým prvkem a jeho sudým následn´ıkem se neodráˇz´ı v pásmu koeficient˚u vypoˇc´ıtaném jedn´ım krokem. Naproti tomu, existuje pˇrekryv mezi následovn´ıky vysoko-frekvenˇcn´ıch Daubechies filtr˚u, takˇze zmˇena mezi dvˇema prvky se odráˇz´ı ve výsledku.

Daubechies D4 je pˇresnˇejˇs´ı, dokud se zmˇena na vstupu odráˇz´ı ve výsledc´ıch v kaˇzdém kroku.

Pouˇzit´ı Daubechies algoritm˚u pˇrináˇs´ı vyˇsˇs´ı nároky na výpoˇcty a komplikovanˇejˇs´ı algoritmus.

Z´avis´ı na typu aplikace, zdali se v˚uˇci cenˇe vyplat´ı vyˇsˇs´ı pˇresnost Daubechies algoritmu.

(34)

Kapitola 3

Poˇ rizov´ an´ı MR obraz˚ u

Témˇeˇr pˇred 2500 lety zformuloval Hippokrates zásady, na kterých buduje lékaˇrská vˇeda dodnes.

Kromˇe zásad etických zavedl také pojmy diagnostiky a správnˇe rozpoznal, ˇze úspˇech terapie závis´ı na správném rozpoznán´ı pˇr´ıˇciny nemoci, tedy na co moˇzná nejpˇresnˇejˇs´ı diagnostice.

Modern´ı ´era zobrazovac´ı techniky zaˇcala roku 1895, objeven´ım paprsk˚u X W. C. Roentgenem.

Paprsky X, známé jako rentgenové záˇren´ı, vznikaj´ı pˇri interakci rychlých elektron˚u s hmotou a d´ıky velmi krátké vlnové délce jsou schopny prozáˇrit lidské tˇelo. Standardn´ı rentgenová dia- gnostika proˇsla v pr˚ubˇehu dvacátého stolet´ı velkým vývojem, avˇsak ˇc´ım dál v´ıce zˇrejmé byly i jej´ı nedostatky.

Prvn´ı nedostatek byl ten, ˇze orgány, kterými paprsek postupnˇe proˇsel, byly znázornˇeny sumárnˇe.

Jejich obrazy se pˇrekrývaly. Druhý, v souˇcasné dobˇe stále v´ıce obecnému povˇedom´ı zˇrejmý nedostatek, je vedlejˇs´ı úˇcinek rentgenového záˇren´ı, který pacienta ohroˇzuje a proto je snaha jeho omezen´ı na nejmenˇs´ı moˇznou m´ıru.

Tyto nedostatky se snaˇz´ı ˇreˇsit modern´ı tomografick´e metody.

Metoda poˇc´ıtaˇcové tomografie umoˇzˇnuje pˇresné zobrazen´ı vyˇsetˇrovaných lidských orgán˚u a do- konalé vzájemné rozliˇsen´ı jednotlivých typ˚u lidských tkán´ı. Jej´ı nevýhodou vˇsak z˚ustává uˇzit´ı rentgenového záˇren´ı. To m˚uˇze modifikovat nebo zniˇcit genový kód ozáˇrených bunˇek a zp˚usobit tak vznik zhoubného rakovinného bujen´ı.

Metoda magnetické rezonanˇcn´ı tomografie nemá na rozd´ıl od poˇc´ıtaˇcové tomografie za následek prakticky ˇzádné neˇzádouc´ı úˇcinky na vyˇsetˇrovaného pacienta. Vytváˇr´ı vysoce kvalitn´ı sn´ımky lidského tˇela. Nevýhodou této metody je jej´ı stále jeˇstˇe velká nákladnost, která vede mnohdy pouze ke komplementárn´ımu nasazen´ı této metody [ZI00].

(35)

3.1 Princip magnetick´ e rezonance

Kaˇzdý atom se skládá z jádra, tvoˇreného protony a neutrony, a z elektron˚u, které krouˇz´ı po elektronových drahách, Vˇsechny protony i neutrony rotuj´ı s vysokou frekvenc´ı kolem své osy.

Obsahuje-li atomové jádro lichý poˇcet proton˚u a neutron˚u, je jeho rotaˇcn´ı impuls, který je roven souˇctu rotaˇcn´ıch impuls˚u vˇsech jeho ˇcást´ı, nenulový, coˇz má za následek rotaci celého atomového jádra. Tento rotaˇcn´ı impuls oznaˇc´ıme I a nazveme jaderný spin. Jaderný spin m˚uˇze nabývat pouze diskrétn´ıch hodnot, daných vztahem:

|I| = hpI_q(I_q+ 1) (3.1)

kde h je Planckova konstanta, h = 1, 0545 · 10⁻³⁴Js, a Iq je kvantové ˇc´ıslo jaderného spinu, které m˚uˇze nabývat pouze diskrétn´ıch hodnot celoˇc´ıselného násobku jedné poloviny (tedy 1/2, 1, 3/2 atd.).

Kaˇzdý jaderný spin vykazuje magnetický moment µ.

µ = γI (3.2)

kde γ je gyromagnetická konstanta. Nen´ı-li atomové jádro vystaveno p˚usoben´ı vnˇejˇs´ıho pole, jsou vˇsechny orientace jeho magnetického momentu energeticky rovnocenné. Vystav´ıme-li je vˇsak p˚usoben´ı statického homogenn´ıho pole o intenzitˇe B0, které je orientováno ve smˇeru osy z, zvýˇs´ı se t´ım energie jádra o pˇr´ıdavnou potenciáln´ı energii o velikosti

E = −µzB0. (3.3)

Zmˇena energie prob´ıhá nespojitˇe, pˇreskokem mezi sousedn´ımi energetickými hladinami. V izolo- vaném homogenn´ım magnetickém poli jsou pˇreskoky mezi energ. hladinami z d˚uvodu platnosti zákona o zachován´ı energie zakázány. Mohou být indukovány pouze interakc´ı, tzv. magnetic- kou rezonanc´ı, atomového jádra s vysokofrekvenˇcn´ım elektromagneickým signálem o frekvenci rovné frekvenci rotace magnetického momentu atomového jádra [ZI00].

Z hlediska klasické fyziky je atomové jádro s lichým poˇctem proton˚u a neutron˚u elementárn´ı koule rotuj´ıc´ı kolem své vlastn´ı osy, která má kladný elektrický náboj.

Lidské tˇelo je sloˇzeno ze ˇctyˇr základn´ıch stavebn´ıch kamen˚u, a to z vod´ıku, uhl´ıku, kysl´ıku a dus´ıku, pˇriˇcemˇz atom vod´ıku je v lidském tˇele zastoupen nejˇcastˇeji.

Mnoˇzina vˇsech proton˚u vod´ıku v lidském tˇele nevykazuje ˇzádný magnetický moment, nebot’

magnetické momenty jednotlivých proton˚u se d´ıky zcela náhodnému rozdˇelen´ı smˇer˚u os rotace navzájem ruˇs´ı. Tento stav se zmˇen´ı, vystav´ıme-li lidské tˇelo p˚usoben´ı silného homogenn´ıho mag- netického pole. Magnetické momenty jednotlivých spin˚u se nasmˇeruj´ı rovnobˇeˇznˇe k silokˇrivkám p˚usob´ıc´ıho pole, a to jak paralelnˇe, tedy do smˇeru vektoru p˚usob´ıc´ıho pole, tak i antiparalelnˇe, tedy do smˇeru opaˇcného.

(36)

Tyto dva stavy maj´ı rozd´ılnou energetickou ´uroveˇn a jsou v tˇele nerovnomˇernˇe zastoupeny.

Antiparaleln´ı uspoˇrádán´ı je energeticky nároˇcnˇejˇs´ım stavem. Pˇri intenzitˇe vnˇejˇs´ıho magnetického pole rovné 1, 5 Tesla odpov´ıdá 1000000 antiparalelnˇe uspoˇrádaných proton˚u zhruba 1000007 proton˚u uspoˇrádaných paralelnˇe.

Uváˇz´ıme-li, ˇze v jednom krychlovém centimetru lidské tkánˇe je obsaˇzeno zhruba 10¹⁶ proton˚u, tak se kaˇzdý objemový element pacientova tˇela stane d´ıky takto vzniklé magnetické nerovnováze magnetickým dipólem s mˇeˇritelným magnetickým momentem.

Jedinou cestou, jak tento moment zmˇeˇrit, je zmˇena jeho orientace.

Protony rotuj´ı s precesn´ı frekvenc´ı ω. Vyˇsleme-li kolmo ke smˇeru vnˇejˇs´ıho magnetického pole impuls o kmitoˇctu rovném právˇe ω, dojde v d˚usledku rezonance k pˇrenosu ˇcásti energie na protony vykonávaj´ıc´ı pohyb. To se projev´ı pˇreskokem urˇcitého poˇctu proton˚u z niˇzˇs´ı energetické hladiny na hladinu vyˇsˇs´ı. Vyslán´ı impulsu má za následek jeˇstˇe dalˇs´ı zmˇenu pohybu proton˚u.

Výsledný magnetický moment, jehoˇz výchoz´ı orientace je ve smˇeru osy z se vlivem rezonance s elektromagnetickým impulsem rozloˇz´ı na dvˇe sloˇzky. Pomˇer obou sloˇzek je dán intenzitou vyslaného impulsu. Tuto intenzitu je moˇzno volit tak, aby vektorová sloˇzka ve smˇeru osy z zcela vymizela.

Um´ıst´ıme-li nyn´ı v bl´ızkosti rotuj´ıc´ıho magnetického pole indukˇcn´ı c´ıvku, indukuje se v této c´ıvce stˇr´ıdavý elektrický proud úmˇerný intenzitˇe tohoto pole. Zmˇeˇren´ım intenzity indukovaného proudu je pak moˇzno hledanou intenzitu magnetického pole stanovit.

Po odeznˇen´ı elektromagnetického impulsu se rotuj´ıc´ı protony nacházej´ı v energeticky nevyváˇzeném stavu a zaˇcnou se vracet do stavu energetické rovnováhy, do stavu pˇred vyslán´ım impulsu [ZI00].