Modellering av ultraljudsbilder med tv ˚adimensionella AR-modeller

(1)

2005:136 CIV

E X A M E N S A R B E T E

Modellering av ultraljudsbilder med tvådimensionella AR-modeller

Johan Svensson

Luleå tekniska universitet Civilingenjörsprogrammet

Teknisk fysik

Institutionen för Systemteknik EISLAB

2005:136 CIV - ISSN: 1402-1617 - ISRN: LTU-EX--05/136--SE

(2)

Modellering av ultraljudsbilder med tv ˚adimensionella AR-modeller

Johan Svensson

Lule˚ a tekniska universitet Institutionen f¨ or Systemteknik

29 april 2005

(3)

(4)

S AMMANFATTNING

M˚alet med det här examensarbetet var att p˚a datorsimulerat vis placera ut ett antal spridare och göra ultraljudsavbildning p˚a dessa, och därefter ur ultraljudsbilden försöka skatta hur m˚anga spridarna var. Detta är en början p˚a ett forskningsarbete som har som m˚al att kunna mäta porösitet med hjälp av ultraljud. Tanken är att spridarna ska vara en förenklad modell av porösitet.

För att förhoppningsvis kunna skatta egenskaperna i bilden, skattas bildens autoregres- siva (AR) motsvarighet. AR-modellen skattas i sin tur med hjälp av minsta kvadratmetoden.

För att testa om det gick att skatta antalet spridare gjordes här en stor mängd ult- raljudsavbildningar med olika antal spridare P˚a var och en av avbildningarna görs en AR-analys. Men det visade sig att alla AR-modeller är s˚a gott som lika oberoende av hur m˚anga spridare som använts. Det skulle kunna vara p˚a grund av att n˚agot är gjort p˚a fel sätt men det troligaste är att AR-modellen inte lämpar sig för detta ändam˚al.

En annan fr˚aga som kommit fram under arbetets g˚ang är hur porösitet egentligen ska modelleras. Att använda sig av en massa spridare som är slumpmässigt spridda över ett omr˚ade är förmodligen inte en representativ modell.

Arbetet har trots detta resulterat i en modell hur simulering och ultraljudsavbildning med efterf¨oljande analys ska g˚a till.

iii

(5)

(6)

A BSTRACT

The goal of this master’s thesis was to make a computer simulation with a number of ultrasound scatterers, and then perform ultrasound imaging of these, and finally, from the ultrasound image decide how many scatterers that were used. This is the beginning of a scientific work where the goal is to be able to measure the porosity of materials with ultrasound. The main idea is that scatterers should be a simplified model of the porosity.

To be able to estimate the properties of the image, the image’s autoregressive (AR) counterpart is first estimated. The AR-model is calculated with help of the least-squares method.

To be able to test if it is possible to estimate the number of scatterers, a large number of ultrasound images was simulated with a varying number of scatterers. An AR model was estimated for each image. But, as it turned out, every AR-model were nearly identical, independent of the number of scatterers. The reason for this could be that something was made in the wrong way, but most probably the AR-model was the wrong solution for this problem.

Another question that has arisen during the work was how the porosity really should be modeled. The use of a large number of scatterers in a random position inside an area is probably the wrong way to describe the porosity.

Despite these problems this work has resulted in a model for making the ultrasound simulation and how to make the following analysis.

v

(7)

(8)

F ¨ ORORD

M˚alet med detta projekt var fr˚an början att verkligen att kunna mäta porösiteten i ett simulerat material och dessutom möjligen att kanske göra det p˚a riktigt. Men sedan visade sig problemet vara betydligt sv˚arare än väntat, och arbetet begränsades till att testa om det gick att skatta antalet spridare fr˚an en simulerad ultraljudsbild med hjälp av en AR-modell.

Jag vill passa p˚a att tacka Johan Carlson som är en utmärkt handledare, lätt att prata med och bra p˚a att ge tips. Jag kommer speciellt att minnas en dag sommaren 2004 d˚a jag och Johan tog en tur ut till Magnus Lundbergs sommarstuga och hade en lektion om tv˚adimensionella processer.

Jag vill ocks˚a tacka min bror Kennet Svensson för det ekonomiska stödet de sista m˚anad- erna, annars hade det kanske inte varit möjligt att avsluta arbetet.

vii

(9)

(10)

I ^NNEH ALL ˚

Kapitel 1: Inledning 1

1.1 Bakgrund . . . 1

1.2 M˚al . . . 1

1.3 Avgr¨ansning av arbetet . . . 2

1.4 Uppl¨agg . . . 2

Kapitel 2: Teori 3 2.1 Syntetisk aperturavbildning . . . 3

2.2 Modellering av specklebild . . . 12

2.3 J¨amf¨orelse av AR-parametrar . . . 22

Kapitel 3: Simuleringar och andra beräkningar 23 3.1 Första körningen . . . 23

3.2 Andra k¨orningen . . . 24

3.3 Tredje k¨orningen . . . 25

3.4 Fj¨arde k¨orningen . . . 26

3.5 Femte k¨orningen . . . 26

Kapitel 4: Resultat 29 4.1 Syntetisk aperturavbildning . . . 29

4.2 Parameterskattning och j¨amf¨orelser av AR-parametrar . . . 31

4.3 Diskussion . . . 36

Kapitel 5: Slutsatser 41

(11)

(12)

K APITEL 1 Inledning

1.1 Bakgrund

Inom medicinen utvecklas keramiska material som ska fylla ut där ben saknas och samtidigt stimulera ˚aterväxt. Dessa material behöver allts˚a vara porösa. Men porösa material

är ocks˚a svagare, och bencementen f˚ar därför inte vara för porösa.

Porösiteten mäts idag med hjälp av Archimedes princip som är b˚ade förstörande och extremt tidskrävande. Det har visat sig att det g˚ar att mäta vissa mekaniska egenskaper med ultraljud. Fr˚agan är om det g˚ar att mäta porösiteten?

En annan tillämpning som porösitetsmätning med ultraljud skulle kunna användas till

¨ar att hitta sprickor i betong. Sprickorna ter sig d˚a som avvikelser i por¨ositeten.

1.2 M ˚al

Syftet med det här examensarbetet är att ta reda p˚a om en tv˚adimensionell autoregressiv (AR) modell är lämplig för att beskriva ”materialegenskaper” som att skatta antalet punktformiga spridare som använts i en simulering. Detta ska i viss mening motsvara porösiteten, men hur riktig porösitet ska modelleras tas inte upp här. Vidare är syftet med arbetet att f˚a fram en modell, i form av MATLAB-funktioner, för hur datorsimuleringar

1

(13)

2 Inledning av dessa problem ska göras med toolboxen Field II, inkl. efterföljande lobformning och analys. Syftet med en lobformning är att det ska vara möjligt att göra en lokaliserad bedömning av den egenskap som ska skattas.

1.3 Avgr ¨ansning av arbetet

Att i MATLAB programmera nödvändiga verktyg, utföra simulerad ultraljudsavbildning och skatta ultraljudbildens AR-process. Därefter försöka se om AR-parametrarna g˚ar att använda till att skatta antalet spridare, eller andra egenskaper.

1.4 Uppl ¨agg

Arbetet bestod av matematiskt h¨arleda eller begripa matematiska principer, men framf¨or allt att programmera i MATLAB.

(14)

K ^APITEL 2 Teori

I detta kapitel tas den teoretiska bakgrunden till arbetet upp.

Metoden som används här, har i stora drag tre delar. Först ultraljudssimulering, sedan beamforming/lobforming och sist AR-parameterskattning. Först simuleras ultraljudsdata fram fr˚an en beskrivning av hur ett material ser ut. Ultraljudspulser sänds ut och tas emot igen som i en radar eller en sonar. Därefter lobformas detta, dvs en karta som visar varifr˚an ultraljudsekona kommer, skapas. Denna karta visar b˚ade riktning och avst˚and.

Denna karta analyseras av en AR-parameterskattare, som försöker f˚a fram de stokastiska egenskaperna i form av parametrar till en AR-process. Utifr˚an dessa parametrar är det sedan fr˚agan om det g˚ar att skatta n˚agra materialegenskaper.

2.1 Syntetisk aperturavbildning

Ultraljudsmätningar tillsammans med lobformning kallas för ultraljudsavbildning. I ult- raljudssammanhang används vanligen en uppsättning (p˚a engelska: array) av omvandla- relement (p˚a engelska: transducer). Elementen omvandlar elektriska signaler till ultraljud och tillbaka igen till en elektrisk signal. Det finns element bara kan omvandla i en riktning och de som klarar dubbelvägsomvandling, dvs de som är b˚ade ”mikrofon” och ”högtalare”

samtidigt. I denna rapport förutsätts alla element vara dubbelvägselement om inget annat sägs. Med de multipla elementen kan källan för ett eko bestämmas med hjälp av lobformning. Dessutom kan pulsen som sänds göras starkare i bestämda riktningar (dvs lobformas vid utsändningen).

Syntetisk apertur kan egentligen betyda lite olika saker. Dels kan det vara att bara ett 3

(15)

4 Teori

(a) (b)

Elementets rörelse

z

x

Figur 2.1: Enkelelement och multielement. H¨ar s¨ander alla elementen p˚a multielements varianten samtidigt vilket ger en str˚ale riktad rakt fram. Spridarna sprider egentligen lika ˚at alla h˚all, bilden visar bara v˚agorna som studsar tillbaka mot ultraljudutrustningen.

element sänder och alla tar emot (det omvända vore bara slöseri med resurser). Eller s˚a används bara ett element för b˚ade sändning och mottagning, och d˚a behövs ingen uppsättning med element. Aperturen byggs upp efterhand med hjälp av flera enskilda mätningar med elementet i olika positioner i sidled. I detta arbete används enkelelements- metoden. Metoden har fördelen att aperturen lätt kan göras stor och kan ändra form relativt lätt, dessutom är den mycket billigare. Nackdelen är att det som mäts inte f˚ar röra p˚a sig mellan n˚agon mätning, men eftersom materialet här inte rör p˚a sig är inte detta n˚agot problem. Se figur 2.1.

2.1.1 Datainsamling

De simulerade ultraljudsmätningarna är gjorda med ett simuleringspaket som heter Field II, som simulerar ultraljudutrustning inkl objekt. Field II är utvecklat Jørgen Arendt Jensen p˚a Danmarks Tekniske Universitet, och är en gratis verktygsl˚ada (toolbox) till MATLAB. Version 3.0 fr˚an 17 April, 2002, är det som har används.

Field II g˚ar att hitta p˚a http://www.es.oersted.dtu.dk/staff/jaj/field/. Den som vill ha en n¨armare introduktion om Field kan ocks˚a l¨asa [1].

Field använder sig av relativt komplexa beräkningar för att f˚a fram hur ljudv˚agorna fortplantar sig. Dessa sker i tre dimensioner även om objekten här bara är modellerade i

(16)

2.1. Syntetisk aperturavbildning 5

tv˚a dimensioner, sidled, x, och djupled, z (se figur 2.1).

Det finns ett antal standardelement och flerelementsupps¨attningar i vilka geometrin kan

ändras. Här används en enkelelements ”piston”, vilket är ett enkelt element i form av en cylinder. Elementet är inte flyttbart i Field, utan är alltid positionerat i origo, rik- tat i z-riktningen. Istället för elementet f˚ar objekten flyttas. Objekten best˚ar enbart av punktformiga spridare som sprider v˚agen sfäriskt. Spridarna har inställningsbar sprid- ningsstyrka. Ett stort antal s˚adana används s˚a att det inte g˚ar att se en enskild spridare.

Vidare s˚a ska den elektriska puls som sänds till elementet och överföringsfunktionen hos elementet definieras. Även samplingsfrekvens och ljudhastighet ska ställas in.

När allt är rätt inställt och i rätt position ges ett beräkningskommando och en signal med ekon f˚as som svar. Detta repeteras för alla positioner och en datamängd som kallas för RF-data byggs upp, där RF är en förkortning av radiofrequency. Förmodligen för att det oftast handlar om höga frekvenser och förmodligen ocks˚a att mycket kommer ifr˚an radarteknologin. RF-data där elementet rört sig längs en linje kallas ocks˚a ofta för B-Scan. (Där elementet st˚att still kallas det för A-Scan, och där elementet rört sig fram och tillbaka över en yta s˚a kallas det för C-Scan.)

Ett praktiskt problem är att Field inte börjar svarssignalen förrän signalen skiljer sig fr˚an noll (den hör inte sig själv vid utsändningen). och slutar inte förrän den ˚ater klingat av mot noll (med viss noggrannhet). I svaret fr˚an Field f˚as, förutom signalen, en tidangivelse som säger när signalen började. Problemet är att beroende p˚a hur spridarna befinner sig i förh˚allande till elementet s˚a kommer signalerna att bli olika l˚anga och datastrukturen för RF-data är ordnad som en matris. Dvs, som det är gjort här är varje kolumn en mätning och varje rad en sampling som ska ha samma tidsoffset i alla mätningar. Alla kolumnerna m˚aste ha samma storlek. Problemet löses genom att ha tv˚a ”dummyspridare” som b˚ada stannar kvar mitt framför elementet. En som är närmast av alla punkter, men inte s˚a nära att det blir beräkningsfel. Och en som är längre bort än alla andra spridare i alla mätningar (hänsyn m˚aste tas till diagonalen). Se figur 2.2.

2.1.2 Lobformning

M˚alet med lobformningen är att utifr˚an RF-data göra en rekonstruktion av var ifr˚an ekona kommer, i form av en bild eller en matris i MATLAB. Denna rekonstruktion kallas för ultraljudsbild eller specklebild. En spridare som ger upphov till ett eko som kommer vara olika l˚angt borta fr˚an elementet beroende vilken x-position elementet kommer vara i. I RF-data ser en spridare ut som en v˚agfront i form av en hyperbel.

(17)

6 Teori

Elementets rörelse

Dummyspridarna följer med

Området som ska lobformas

Dummyspridarna följer med z

x

Figur 2.2: Field ger signaler med variabel längd och start som beror p˚a den första och sista samplingen som är skild fr˚an noll. För att f˚a signaler med samma start och längd används vad som kan kallas ”dummyspridare”. Dessa följer med elementet, men h˚aller samma avst˚and fr˚an denna. Tv˚a stycken används, en som är närmare och en som är längre bort än alla datapunkter, oberoende av vilken x-position elementet befinner sig i (inom ett omr˚ade). Cirkelsektorn i figuren visar att den längsta diagonalen bör beaktas när den bortre positionen ska bestämmas.

Bilden beräknas genom att summera längs varje motsvarande hyperbel för varje bildpunkt i specklebilden. Se figur 2.3 och 2.4.

För att kunna göra en s˚adan summering s˚a m˚aste tiderna till hyperbeln för varje mätning beräknas. Tiden för en mätning kommer att vara ∆t = ∆t1+ ∆t2, där ∆t1 är tiden för pulsen att n˚a spridaren och ∆t₂ är tiden fr˚an spridaren till elementet. Används syntetisk apertur är vägen fram och tillbaka densamma och allts˚a är ∆t₁ = ∆t₂. Tiden enkel väg

är ∆t₁ = d(m, x, z)/c, där d(m, x, z) är avst˚andet mellan elementet och spridaren som en funktion av spridarens x och z-koordinat samt elementets position m. Avst˚andet delas med ljudhastigheten c. Positionen m motsvaras, om mätpunkterna är lika fördelade längs x-axeln, av x_s = m∆x. Avst˚andet d(m, x, z)/c beräknas med Pythagoras sats:

d²(m, x, z) = z² + (x − m∆x)² ⇒

⇒ d(m, x, z) =p

z²+ (x − m∆x)². (2.1)

(18)

P

RF-data Speckle-bild

Tid Avstånd

(a) (b)

Elementets rörelse Elementets rörelse

Figur 2.3: Summering l¨angs en hyperbel i RF-data (a) till en punkt i specklebilden (b). Bilden

¨ar gjord p˚a material fr˚an [2].

Spridare

Figur 2.4: Principen för hur lobforming fungerar. De gr˚aa förtjockningarna är fördröjningar.

Om ett eko kommer fr˚an punkten som fördröjningsoperationen fokuserar p˚a, kommer dessa att bli jämsides. S˚a som det är gjort i lobformaren i detta arbete, summeras dock bara ett sampel fr˚an varje mätning. Denna bild har sitt ursprung i [3].

(19)

8 Teori

Position i sidled (mm)

Avstånd från elementet (mm)

Med kompensering

−20 −10 0 10 20

20 22 24 26 28 30 32 34 36 38 40

Utan kompensering

−20 −10 0 10 20

20 22 24 26 28 30 32 34 36 38 40

Figur 2.5: Bilder med och utan pulslängdskompensation. Cirklarna visar positionen där spridarna är placerade enligt bildernas koordinater. Skalan till vänster visar djupet som lobformaren

”tror” att den fokuserar p˚a. I den högra bilden stämmer inte detta överens med verkligheten – fokuseringen sker närmare än vad som egentligen är meningen. En l˚ang pulslängd har använts här s˚a att fenomenet blir tydligt.

Kombinerar man ∆t(m, x, z) = ∆t₁+ ∆t₂ med ekvation 2.1:

∆t(m, x, z) = ∆t₁+ ∆t₂ = d(m, x, z)

c +d(m, x, z)

c = 2d(m, x, z)

c .

Fr˚an Field f˚as en starttid som anger när den första samplingen togs emot, mätt fr˚an starten av utsändningen. Om den tiden kallas för t0, s˚a ska tiden för ekot vara ∆t − t0 in i en viss mätning. Denna tid mäter dock tiden till ekots början – inte mitten – vilket kan ge fokusproblem. Oftast är pulsen s˚a kort att det inte spelar s˚a stor roll men är pulsen l˚ang hamnar fokus för nära och spridarna ser ut att ligga längre bort än vad de egentligen gör, samt att de kan se lite ofokuserade ut, se figur 2.5. För att korrigera tiden s˚a m˚aste den halva mottagna pulsens längd, t_puls/2, adderas till ∆t:

∆t(m, x, z) − t₀+ t_puls/2.

Denna puls är faltad med elementets sändöverföringsfunktion och igen med mottag- ningsöverföringsfunktionen.

Den fördröjning som d˚a ska göras är:

∆t_{ef f}(m, x, z) = ∆t(m, x, z) − t₀+ t_puls/2.

(20)

2.1. Syntetisk aperturavbildning 9 Men det m˚aste göras om till samplingsindex, n, om lobformningen ska kunna göras p˚a RF- data i MATLAB. Om fsär samplingsfrekvensen s˚a är förh˚allandet mellan samplingsindex och tid:

∆t_{ef f}(m, x, z) = n f_s.

F¨ordr¨ojningen i samplingar blir d˚a n = f_s∆t(m, x, z) − f_st₀+ f_st_puls/2. Eftersom ett index m˚aste vara ett heltal s˚a avrundas uttrycket: n = bf_s∆t(m, x, z) − f_st₀+ f_st_puls/2+ 1/2c.

Klamrarna bxc betyder avrundning ner˚at – ”floor” p˚a engelska – och bx + 1/2c blir d˚a

”vanlig” avrundning.

Om man sl˚ar ihop allt ovanf¨or i ett uttryck f¨or n f˚as:

n =

$

2fspz²+ (x − m∆x)²

c − f_st₀+ f_st_puls/2+ 1/2

% ,

alternativt

n =



2 s

f_s cz

2

+

(x − m∆x)f_s c

2

− f_st₀+ f_st_puls/2+ 1/2



.

Hur bilderna p˚a RF-data är ritade i denna rapport förklaras i figur 2.6 och 2.7. Hur specklebilderna är gjorda förklaras i figur 2.8.

I m˚anga av figurerna med RF-data och specklebilder, som figur 2.7 och 2.8, beräknas enveloppen av signalerna. Detta görs p˚a följande sätt:

x_envelopp[n] = |x[n] + jH {x[n]}|,

där H är Hilberttransformen och j den imaginära enheten. Den som vill veta mer om Hilberttransformen kan exempelvis läsa [4].

Om ∆x är större än halva v˚aglängden,_2f^c , uppst˚ar spatiell vikning. Dvs ”rumslig” vikning – i motsats till tidsmässig vikning. Vikningen förklaras bättre i figur 2.9.

2.1.3 D ¨ampningskompensation

Eftersom ultraljudsekon som kommer fr˚an spridare som ligger l¨angre bort ¨ar svagare

än de som är närmare, s˚a m˚aste RF-data korrigeras om det ska g˚a att f˚a n˚agorlunda stationärt stokastiskt data att analysera. Detta m˚aste göras innan lobformningen d˚a den förstör en del information.

(21)

10 Teori

Tid från utsändning (µs)

−40 −30 −20 −10 0 10 20 30 40

30

40

50

60

70

80

90

100

Figur 2.6: RF-data. Bilden visar enveloppen av signalerna. Värdena har sedan logaritmerats s˚a att det g˚ar att se alla detaljer i bilden. Det g˚ar att se effekten av dummy-punkterna upptill och nertill i bilden. Omr˚adet innanför den streckade linjen är där lobformaren kommer att fokusera.

Här har tv˚a metoder använts, den ena mer vetenskaplig än den andra, men de funkar

¨

and˚a n˚agorlunda bra.

Den ena metoden är en ”quick’n’dirty” variansnormalisering där standardavvikelsen över alla mätningar för given tidpunkt beräknas. Dvs standardavvikelsen beräknas som en (diskret) funktion av tiden. Sedan normaliseras RF-data s˚a att standardavvikelsen blir den samma över alla tidpunkter.

Den andra metoden för dämpningskompensation placerar en spridare i olika positioner och mäter hur mycket av signalen som kommer tillbaka. Punkten varieras b˚ade i z-led och x-led, p˚a s˚a vis f˚as en bild av ljudstr˚alen som sänds ut av elementet. bilden summeras i x-led, och en reflektionskoefficient f˚as som en funktion av avst˚andet. Se figur 2.10. En signal som sänds ut och dämpas likformigt kontinuerligt av mediet den färdas i, f˚ar en exponentiell funktion, e^−αz.

(22)

(a)

−15 −10 −5 0 5 10 15

30

35

40

45

50

(b)

−15 −10 −5 0 5 10 15

30

35

40

45

50

(c)

−15 −10 −5 0 5 10 15

30

35

40

45

50

Figur 2.7: RF-dataomr˚adet innanför den streckade linjen. Bilderna visar olika sätt att visa RF-data. Bild (a) är r˚a obehandlad data. Det är dessa som beräkningarna sker p˚a. Bild (b) är enveloppen av samma data. Bild (c) är samma som (b) fast bildpunkterna är logaritmerade.

(23)

12 Teori

(a)

−15 −10 −5 0 5 10 15

20 22 24 26 28 30 32 34 36 38 40

(b)

−15 −10 −5 0 5 10 15

20 22 24 26 28 30 32 34 36 38 40

Figur 2.8: Specklebilder. Bild (a) är r˚a obehandlad data. Bild (b) är enveloppen av (a). De suddiga strecken vid varje punkt är en mild effekt av spatiell vikning d˚a en kort puls har stor frekvensspridning.

Förlusten i mottagen signal är dock större än vad dämpningen orsakar, pga av att str˚alen dessutom sprids. Men det bör inte vara n˚agot problem om inte ultraljudspulsen har allt för l˚ag frekvens relativt hur stort elementet är, för d˚a sprids str˚alen mera sfäriskt.

2.2 Modellering av specklebild

Det som ska undersökas är om en specklebild kan betraktas som en tv˚adimensionell stokastisk autoregressiv process. För att kunna jämföra tv˚a realiseringar av specklebilderna s˚a ska parametrarna till motsvarande AR-process skattas. Förhoppningsvis g˚ar det att ur parametrarna säga n˚agonting om vilka villkor specklebilderna kommer ur, exempelvis hur m˚anga spridare som bilderna är gjorda med.

2.2.1 Endimensionell AR-modell

En AR-process är en stokastisk process som kan beskrivas som en viktad summa av sina tidigare värden plus ett vitt (gaussiskt) matningsbrus, ”felet”, hos processen. Processen har pga ˚aterkopplingen ett oändligt impulssvar, och är därför det samma som ett IIR-

(24)

2.2. Modellering av specklebild 13

−20 −15 −10 −5 0 5 10 15 20

20

22

24

26

28

30

32

34

36

38

40

Figur 2.9: Spatiell vikning. Tv˚a spridare finns i bilden, en i (-15,0 22,0) mm – nästan osynlig, och den andra i (15,0 38,0) mm. Speciellt den andra spridaren uppvisar stora sidlober som suddiga band förmodligen i form en hyperbel. Bandet är fläckvis starkare där ett maximum för en sidlob finns.

filter (Infinite Impulse Response) som inte har nollställen i z-transformdomänen, dvs täljaren i filtrets z-transform motsvarighet är en konstant skild fr˚an noll.

Matematiskt ¨ar definitionen f¨or en endimensionell AR-process y(n) enligt [5]:

y(n) =

p

X

k=1

a_k· y(n − k) + ε(n) ∀n, (2.2)

där a_k är AR-parametrarna, ε(n) det vita matningsbruset och p är ordningen p˚a processen. Väntevärdena ska vara

E[ε(n)] = 0, E[ε(n)²] = σ², E[ε(n₁)y(n₂)] = 0, n₂ < n₁,

där σ är standardavvikelsen. AR-processen använder de p senaste utvärdena och nuva- rande värde p˚a ε(n) för att generera ett nytt utvärde.

(25)

14 Teori

(a)

−20 −10 0 10 20

20 22 24 26 28 30 32 34 36 38

(b)

−20 −10 0 10 20

20 22 24 26 28 30 32 34 36 38

20 25 30 35

1 1.2 1.4 1.6 1.8 2 2.2 2.4

x 10⁻¹⁸ (c)

Styrka

Anpassning Mätning

20 25 30 35

0.6 0.7 0.8 0.9 1 1.1 1.2 1.3 1.4 1.5

(d)

Styrka

Figur 2.10: Kompensation för dämpning. Bild (a) visar hur mycket eko som f˚as tillbaka om en spridare placeras i olika positioner framför elementet. Det g˚ar att se effekten av b˚ade str˚alens spridning och ren dämpning i bilden. De ljusa ränderna tros bero p˚a interferens i elementet, annars är orsaken oklar. Kurvan (c) är styrkan i (a) summerad radvis. Till den har en expo- nentialfunktion, c · e^−α·z, anpassats. Den funktionen används sedan för att jämna ut styrkan.

Bild (b) och (d) visar effekten av detta.

(26)

D +

+ +

D

y[n − 3]

y[n − 1]

a3

a2

ε[n]

a1

y[n − 2]

D

y[n]

Figur 2.11: AR-3 process (p=3). Insignalen ε[n] är vitt gaussiskt brus. Blocken med ”D” är fördröjningsblock. Där det är ett a ovanför, multipliceras signalen med denna. Figuren är hämtad ur [4].

Figur 2.11 visar ett blockschema av en AR-3 process. En AR-process kan vara instabil, vilket den ¨ar om processens z-transforms poler inte befinner sig innanf¨or enhetscirkeln.

Dvs absolutvärdet av rötterna till ekvationen, z^p− a₁z^p−1. . . − a_p−2z²− a_p−1z − a_p = 0, ska vara mindre än ett.

2.2.2 Tv ˚adimensionell AR-modell

En tv˚adimensionell AR-process fungerar p˚a ungefär samma sätt som endimensionell men AR-parametrarna är nu uppställda i tv˚a dimensioner ungefär som en matris, men som inte behöver vara kvadratisk. Kausaliteten gör att det f˚ar inte finnas parametrar för

”framtida” värden. Tv˚a vanliga uppställningar som är kausala är NSHP (Nonsymmetrical Half Plane) och QP (Quarter Plane). Se figur 2.12

Den matematiska definitionen av en tv˚adimensionell AR-process ¨ar enligt [6]:

y(n, m) = X

k1,k2∈S_qp

a_k₁_,k₂ · y(n − k₁, m − k₂) + ε(n, m), ∀n, m. (2.3)

(27)

16 Teori







a4,3 a4,2 a4,1 a4,0 a4,−1 a4,−2 a4,−3

a3,3 a3,2 a3,1 a3,0 a3,−1 a3,−2 a3,−3

a_2,3 a_2,2 a_2,1 a_2,0 a_2,−1 a_2,−2 a_2,−3 a_1,3 a_1,2 a_1,1 a_1,0 a_1,−1 a_1,−2 a_1,−3 a0,3 a0,2 a0,1 •













a4,3 a4,2 a4,1 a4,0

a3,3 a3,2 a3,1 a3,0

a_2,3 a_2,2 a_2,1 a_2,0 a_1,3 a_1,2 a_1,1 a_1,0 a0,3 a0,2 a0,1 •







NSHP QP

Figur 2.12: Exempel p˚a NSHP resp. QP AR-parameteruppställningar. NSHP är den mest ge- nerella kausala beskrivningen, men QP är vanligare och mer teori finns beskriven. Observera att i en del figurer i denna rapport, är indexen negerade.

Om NSHP anv¨ands (variant av) g¨aller,

S_nshp = {k₁, k₂| k₁ = 0, 1 ≤ k₂ ≤ p; 1 ≤ k₁ ≤ q, −p ≤ k₂ ≤ p}, och om QP anv¨ands,

S_qp = {k₁, k₂| k₁ = 0, 1 ≤ k₂ ≤ p − 1; 1 ≤ k₁ ≤ q, 0 ≤ k₂ ≤ p}.

Ar p eller q noll s˚¨ a ¨ar processen endimensionell och ekvivalent med 2.2.

2.2.3 Tv ˚adimensionell stabilitet

Aven en tv˚¨ adimensionell AR-process kan vara instabil, men h¨ar ¨ar teorin mycket sv˚arare.

H¨ar ges bara en fingervisning av problemet.

Det finns en tv˚adimensionell motsvarighet av z-transformen. Eftersom den transformen inneh˚aller koordinater hos tv˚a stycken komplexa variabler, z₁ och z₂, s˚a blir motsvarig- heten för enhetscirkeln, en fyrdimensionell enhetshypersfär. Vidare s˚a blir polynomen som bestämmer var nollställena är, tv˚adimensionella, och teorin att lösa dessa är inte lika välutvecklad eller enkel att först˚a som för endimensionella polynom. Nollställena till s˚adana polynom blir generellt ytor. Exempelvis uttrycket 1−bz₁⁻¹z₂⁻¹är noll d˚a z₂ = b/z₁. Avst˚andet fr˚an origo till n˚agot nollställe hos uttrycket

1 − X

k1,k2∈S_qp

a_k₁_,k₂z₁^−k¹z^−k₂ ² (2.4) f˚ar inte vara st¨orre eller lika med ett n˚agonstans. Dvs det f˚ar inte finnas n˚agot punkt, k(z₁, z₂)k ≥ 1, d¨ar uttrycket 2.4 blir noll.

Eftersom det är utanför ramen för detta arbete att beskriva om ett tv˚adimensionellt filter

¨ar stabilt eller inte s˚a h¨anvisas den intresserade till [7].

(28)

2.2.4 Endimensionell parameterskattning

En AR-skattning av en stokastisk process är ekvivalent med en spektralskattning eftersom en AR-process är ett linjärt filter som matas med vitt brus vilket beskrevs tidigare. AR- parametrarna är filtrets koefficienter. Här används minsta kvadratmetoden för att skatta AR-parametrarna. AR-parameterskattning beskrivs utförligt i [8], men även i [4].

Ekvation 2.2 kan skrivas p˚a matrisform.

y[n] = y[n − 1] y[n − 2] . . . y[n − p]

| {z }

ϕ^T(n)





 a₁ a₂ ... a_p





 + ε[n]

Ar detta gjort f¨¨ or ett antal v¨arden, y[n₁], y[n₂], . . . y[n_N], kan ett ekvationssystem byggas upp:





 y[n₁] y[n2]

... y[n_N]







| {z }

Y

=







y[n₁− 1] y[n₁− 2] . . . y[n₁− p]

y[n2− 1] y[n2− 2] . . . y[n2− p]

... ... . .. ... y[n_N − 1] y[n_N − 2] . . . y[n_N − p]







| {z }

Φ^T





 a₁ a2

... a_p







| {z }

θ

+





 ε[n₁] ε[n2]

... ε[n_N]







| {z }

E

(2.5)

Uttrycket kan skrivas kortare som

Y = Φ^Tθ + E. (2.6)

Detta är ett ekvationssystem som adderas med en störning, E. Det är parametrarna i θ som ska skattas ur Φ och Y . Det görs genom att lösa ekvationssystemet Y = Φ^Tθ.

Men eftersom antalet rader, N , i ekvation 2.5 bör vara mer än antalet parametrar, p, är systemet överbestämt, och lösningen finns d˚a generellt inte. Istället används det θ som ger minst fel i minsta kvadratmening (dvs felet mäts som vanlig norm eller tv˚anorm som det ocks˚a kallas, beräknas med pythagoras sats). Att antalet rader är mer än antalet parametrar är bara bra eftersom känsligheten för störningen E d˚a minskar.

Minsta kvadratl¨osningen f˚ar man genom att multiplicera ekvationen Y = Φ^Tθ med trans-ˆ ponatet av systemmatrisen Φ^T,

ΦY = ΦΦ^Tθ,ˆ (2.7)

vilket är en projektion p˚a ett rum av lägre dimension. Därefter löses systemet och som svar f˚as:

θ = (ΦΦˆ ^T)⁻¹ΦY.

(29)

18 Teori Observera att ekvationen Y = Φ^Tθ i allm¨ˆ anhet är falsk! Inte förrän efter projektionen, ekv. 2.7, gäller ett likhetstecken.

Matrisen som f˚as ur uttrycket (ΦΦ^T)⁻¹Φ kallas för pseudoinversen till Φ^T, förutsatt att ΦΦ^T har full rang, vilket den har om Φ^T har full kolumnrang. Detta betyder att systemet inte är underbestämt. Är systemet underbestämt funkar inte minsta kvadratmetoden. Är systemet exakt bestämt ger minsta kvadratmetoden den exakta lösningen.

Men här är snarare problemet istället enorma och mer än väl överbestämda dataserier som tar för mycket tid om beräkningen skulle ske p˚a all data. Därför används ibland bara en delmängd av raderna hos Y och Φ^T i ekvation 2.5.

När det gäller minnes˚atg˚angen, g˚ar det att använda sig av ett trick. Istället för att bygga upp hela Φ^T p˚a en och samma g˚ang, s˚a byggs bara en eller ett antal rader hos Φ^T. Kalla den matrisen för Φ^T_i , där i är numret p˚a delen av Φ^T som använts. D˚a kan ett mellanresultat beräknas som R_i = Φ_iΦ^T_i . Denna del adderas till en ackumulator R_ack som när hela Φ^T g˚atts igenom blir det samma som ΦΦ^T. Det samma gäller ΦY som delas upp i ri = ΦiYi, där Yi d˚a är ett eller flera värden motsvarande Φi. Kolumnmatrisen riadderas p˚a samma sätt till en ackumulator, r_ack, som till slut blir det samma som ΦY . Om Φ^T är tex 2133139 × 99 stor (faktisk storlek p˚a en matris som använts), och varje element är ett flyttal av typen ”double precision”, ˚atta bytes, s˚a tar den upp 2133139·99·8 = 1689446088 bytes eller ca 1, 7 GB. Matrisen Ri tar upp 99 · 99 · 8 = 78408 bytes, ca 78 kB. Figur 2.13 visar ett exempel.

Den som vill l¨asa mer om minsta kvadratmetoden h¨anvisas till [9] eller [10].

2.2.5 Tv ˚adimensionell parameterskattning

En tv˚adimensionell skattning sker p˚a samma sätt som en endimensionell skattning. Det finns ingen anledning att bygga p˚a en dimension till ekvation 2.5. Parametermatrisen rullas ut rad för rad eller kolumn för kolumn s˚a att den blir endimensionell och motsvarar θ i ekvation 2.5 och 2.6. I MATLAB kan det göras med reshape. Samma sak görs för Y , Φ^T och E. I figur 2.14 ges ett exempel.

Figur 2.15, 2.16 och 2.17 visar ett exempel p˚a en realisering och p˚af¨oljande analys av en AR-process som gjordes i MATLAB f¨or att testa om AR-skattningen fungerade.

(30)

= =

(b) (a)

Φ^T_i

Φi Ri Φi Yi ri

Figur 2.13: Minnessn˚al parameterskattning. Istället för att bygga upp hela Φ p˚a samma g˚ang, s˚a kan en rad eller ett antal rader byggas upp. Och sedan göra en multiplikation med transponatet av denna, (a), och p˚a s˚a sätt f˚a ett delresultat Ri av matrisen R = ΦΦ^T som oftast är betydligt mindre än Φ. Detta delresultat adderas efterhand som Φ g˚as igenom, till en matris som blir R när all data g˚atts igenom. D˚a bör samma sak ske för r där Φ multipliceras med Y , (b).

Matriserna i figuren bygger p˚a matriserna i figur 2.14.

2.2.6 Residual och val av ordning

Vid val av antalet parametrar, undersöks AR-parametrarna själva och n˚agot som kallas för residualen av skattningen.

Om datamängden som analyserats är representerad i Y och Φ och de skattade parametrarna är ˆθ, s˚a är bY = Φ^Tθ prediktionen av Y . I det endimensionella fallet kan det sesˆ som en skattning av vad nästkommande värde kommer att bli med hjälp av föreg˚aende värden och AR-parametrarna. Felet i denna skattning kallas för residualen och är det samma som F = Y − bY . Dvs

F = Y − Φ^Tθ.ˆ (2.8)

Om ˆθ är lika med θ s˚a är felet F det samma som den realisering av det vita drivbruset E, som Y och Φ^T är skapad ur. Detta gäller förutsatt att den antagna modellen är korrekt.

Om residualen inte är helt vit och uppvisar mönster, bör antalet parametrar ökas. Men parametrarnas värden avtar mot noll med avst˚andet fr˚an origo, och är oftast obetydliga efter bara n˚agra f˚a parametrar. Om alla parametrar av signifikant värde är inkluderade utan att residualen vit, tyder det p˚a att en AR-process inte är den rätta beskrivningen.

(31)

20 Teori

(a) QP

(b) θ

(c) Data

(d) Y

(e) Φ^T

Figur 2.14: Exempel p˚a hur utrullning (QP) kan g˚a till. Parametermatrisen (a), som ska skattas, rullas ut till kolumnmatrisen θ. Datamängden som ska analyseras är (c). Kolumnmatrisen Y , f˚ar sina värden fr˚an det rektangulära omr˚adet nere till höger under den tunna streckade linjen.

Matrisen Φ, f˚ar sina värden ur hela datamängden. Värdena i Φ plockas ur datamängden efter parametermatrisens form runt det aktuella värdet (dvs värdet som kompletterar (a) s˚a den blir rektangulär). Det aktuella värdet ska vara samma som värdet p˚a samma rad i Y . Ett värde i datamängden har f˚att en gr˚a ruta och ett annat en liten ring. Detta för att det ska g˚a att se var dessa värden kan ˚aterfinnas i Y och Φ. Samma gäller för cirkeln och krysset i (a) och (b).





−0.4950 0.4950 −0.4950 0.4950 −0.4950 0.4950

−0.4950 0.4950 0





−2.5 −2 −1.5 −1 −0.5 0 0.5

−2.5

−2

−1.5

−1

−0.5

0

0.5

−0.4

−0.3

−0.2

−0.1 0 0.1 0.2 0.3 0.4

Figur 2.15: Tv˚a representationer av AR-parametrar. Eftersom en del parametermatriser är s˚a stora, är det lättare att representera dessa som bildpunkter i en bild. Värdet längst ner till höger som här är noll, saknar betydelse för AR-processen. Men om matrisen ska stabilitetstestas används alla värden och d˚a bör värdet längst ner till höger vara minus ett beroende p˚a att ekvation 2.2 se ut som den gör.

(32)

" 0.0123 0.0060 −0.0061 0.0082 −0.0117 0.0069 −0.0111 −0.0246 −0.0110 0.0064

−0.0027 −0.0137 −0.4805 0.4734 −0.5228 0.0010 0.0016 0.4455 −0.4645 0.4804 0.0021 −0.0087 −0.4859 0.4657 0

#

−4 −3 −2 −1 0

−4

−3

−2

−1

0 −0.4

−0.2 0 0.2 0.4

Figur 2.16: Skattade AR-parametrar, ur realiseringen i figur 2.17, som är gjord med parametrarna i figur 2.15. Här har fem g˚anger fem parametrar skattas, men bara tre g˚anger tre har använts för att generera bruset. De överflödiga blir nära noll. Realiseringen är 64 g˚anger 64 stor.

Vitt matningsbrus

20 40 60

10 20 30 40 50

60 −3

−2

−1 0 1 2 3

AR−filtrerat

20 40 60

10 20 30 40 50 60

−5 0 5

Residual efter analys

20 40 60

10 20 30 40 50 60

−3

−2

−1 0 1 2 3

Figur 2.17: Bild (a) visar en realisering av vitt brus som sedan i bild (b) filtrerats med parametrarna i figur 2.15. Därefter skattas parametrarna som i figur 2.16 och residualen (c) beräknas ur dessa parametrar. Residualen blir lite mindre, eftersom den behöver ett ”insvängnings”-omr˚ade motsvarande AR-matrisen. Realiseringen är 64 g˚anger 64 stor. Residualen är 60 g˚anger 60.

Om spatiella fenomen uppträder tyder det p˚a att processen inte är stationär, vilket är ett krav för AR-processer.

Notera att vid skattningen spelar det ingen roll i vilken ordning raderna kommer till Y och Φ^T. Men när residualen beräknas m˚aste E kunna omvandlas till en bild där bildpunkterna kommer i rätt ordning. Vidare s˚a g˚ar det inte att hoppa över rader i Y och Φ^T.

Figur 2.17 visar en ideal residual.

(33)

22 Teori

2.3 J ¨amf ¨ orelse av AR-parametrar

Jämförelsen av AR-parametrarna fr˚an olika specklebilder sker genom att θ, dvs den utrul- lade varianten av parametrarna, transponeras och läggs rad för rad i en jämförelsematris.

Förutom att den kontrolleras rent visuellt s˚a görs en singulärvärdesuppdelning.

Singulärvärdesuppdelning p˚aminner om egenvärdesuppdelning. Singulärvärdesuppdel- ningen ser ut som följer:

A = U ΣV^T,

där A är en m × n matris, U är en m × m ortogonal matris med normerade kolumnvektorer, V är en n × n ortogonal matris med normerade kolumnvektorer, och Σ är en m × n diagonalmatris. Värdena i diagonalen hos Σ kallas för singulärvärden. Singulärvärdena

är vanligen ordnande med största först, och sedan i sjunkande ordning. Kolumnerna i U och V kallas för singulärvektorer. Skillnaden mellan egenvärdesuppdelning och sin- gulärvärdesuppdelning är att i singulärvärdesuppdelningen behöver inte matrisen vara kvadratisk, och singulärvärdena är alltid icke-negativa reella tal även om A är komplex.

Sedan ¨ar oftast de tv˚a vektormatriserna olika varandra b˚ade i storlek och inneh˚all.

Likheten med egenvärden är att singulärvärdena är kvadratroten av egenvärdena till A^TA eller AA^T. Singulärvektormatrisen U är egenvektorsmatrisen till AA^T, och V är egenvektorsmatrisen till A^TA.

Om det bara finns ett signifikant singulärvärde hos jämförelsematrisen, betyder detta att alla rader i det närmaste (beroende p˚a de andra singulärvärdena) är en multiplikation av den första raden. Har matrisen U inte n˚agot beroende mellan rader för m˚anga eller f˚a spridare, tyder detta p˚a att AR-parametrarna inte beskriver förh˚allande för detta.

I [10] st˚ar det lite om singulärvärdesuppdelning och hur den ska beräknas med numeriskt robusta metoder.

(34)

K ^APITEL 3 Simuleringar och andra

ber ¨akningar

Detta kapitel tar upp vilka beräkningar som gjordes. Den största tiden av detta arbete utgjordes av programmering av MATLAB-funktioner och p˚a dessa gjordes ett antal da- torkörningar. Eftersom det är tv˚adimensionella och ibland tredimensionella datamängder som hanteras s˚a blir beräkningstiden snabbt väldigt stor.

Ljudhastigheten i alla simuleringar är satt till 1540 m/s. Detta är en typisk ljudhastighet i mänsklig vävnad. Elementets radie är alltid 1 mm med matematiska element p˚a en storlek av 0,05 mm. I figur 3.1 visas hur elementet ser ut.

3.1 F ¨ orsta k ¨ orningen

Första simuleringen hade 20 000 spridare och dimensioner enligt figur 3.2. Det var den största enskilda simuleringen och tog ca 34 timmar p˚a en 1,7 GHz Pentium 4 och 768 MiB¹ arbetsminne. Men det skulle ta alldeles för l˚ang tid att göra beräkningar p˚a denna konfi- guration när det behövs flera simuleringar för att kunna göra en statistisk undersökning.

1Enheten ”MiB” st˚ar för mebibyte, 1048576 eller 2²⁰bytes. Detta för att göra skillnad fr˚an SI prefixet M dvs mega eller 10⁶. Vidare finns kiB, kibibyte eller 1024 och GiB, gibibyte eller 2³⁰. Detta är en IEC standard.

23

(35)

24 Simuleringar och andra ber¨akningar

−1

−0.5 0

0.5 1

−1

−0.5 0 0.5 1

−1

−0.5 0 0.5 1

y (mm)

x (mm)

z (mm)

Figur 3.1: Elementet som användes till alla simuleringar. Det som visas är cirkelformade mem- branet. Rutorna är de matematiska elementen som använts för att beräkna hur elementets membran rört p˚a sig.

3.2 Andra k ¨ orningen

Ovriga simuleringar hade mindre uppl¨¨ osning men gjordes i grupp om minst 25 stycken med olika variationer av spridare. De hade 10 000 till 50 000 spridare, i steg om 10 000.

För att testa om det utgjorde n˚agon skillnad s˚a slumpades den reflektiva styrkan p˚a spridarna olika mellan simuleringarna, alla med ett väntevärde p˚a ett och 98% inom ett intervall som dels var 0,2-1,8 sedan 0,4-1,6, 0,6-1,4, 0,8-1,2, och till sist en med konstant styrka. Allts˚a fem olika antal spridare och fem olika variationer av amplituder, totalt 25 st. Figur 3.3 visar ett exempel p˚a RF-data fr˚an andra körningen.

Andra till femte k¨orningen gjordes p˚a en dubbelprocessors 2,8 GHz Xeon-maskin med 4 GiB arbetsminne. Dessa simuleringar tog, trots en snabb maskin, dagar till veckor att slutf¨ora.

(36)

3.3. Tredje k¨orningen 25

−20 −15 −10 −5 0 5 10 15 20

30

40

50

60

70

80

Figur 3.2: RF-data med 20 000 spridare, alla med samma styrka. Figuren är av samma typ som figur 2.6, men omr˚adet som ska lobformas g˚ar ända ut i vänster och höger kant. En ultraljudspuls

är sänd och mottagen i 2078 mätpunkter. Samplingsfrekvensen är 40 MHz. Signalen som sändes var p˚a 10 MHz.

3.3 Tredje k ¨ orningen

Tredje försöket är det samma som det andra men där de 25 enskilda simuleringarna har repeterats ytterligare fyra g˚anger. Detta s˚a att det ska g˚a att se hur en enskild konfi- guration varierar. Det blev inte s˚a mycket efterarbete p˚a denna eftersom kanteffekterna upptäcktes vara för stora i b˚ade specklebilden och residualen. Kanteffekten kommer sig av att omr˚adet som ultraljudselementet söker av är lika stor som omr˚adet (i x-led) där spridarna finns. Spridare som är p˚a gränsen av omr˚adet kommer efter lobforming att ha sämre fokus i x-led och se utsmetade ut. Specklebilden blir d˚a spatiellt beroende i x-led.

Att AR-modellen f˚ar problem med detta syns i residualen i figur 4.5. Residualen beskrivs i avsnitt 2.2.6.

(37)

26 Simuleringar och andra ber¨akningar

−10 −8 −6 −4 −2 0 2 4 6 8 10

30

35

40

45

50

55

Figur 3.3: Exempel fr˚an simuleringsomg˚ang tv˚a. RF-data med 20 000 spridare, alla med samma styrka. Figuren är av samma typ som figur 2.6, men omr˚adet som ska lobformas g˚ar ända ut i vänster och höger kant. En ultraljudspuls är sänd och mottagen i 260 mätpunkter. Samplings- frekvensen är 50 MHz. Signalen som sändes var p˚a 5 MHz.

3.4 Fj ¨arde k ¨ orningen

I den fjärde ändrades simuleringarna s˚a att elementet g˚ar en bit utanför spridarnas omr˚ade. Detta gör att kanteffekterna, som beskrevs tidigare, blir mindre. Här repete- rades de 25 simuleringarna fyra g˚anger. RF-data blev som i figur 3.4.

3.5 Femte k ¨ orningen

L˚angt efter att den fjärde körningen var klar upptäcktes att i vissa simuleringar var spridarna placerade p˚a samma koordinater. Det visade sig bero p˚a att MATLAB inte gör en ”randomize”-operation vid starten. MATLAB liksom nästan alla andra datorprogram använder pseudoslumpgeneratorer i form av en algoritm som utifr˚an ett inre tillst˚and

(38)

3.5. Femte k¨orningen 27

−20 −15 −10 −5 0 5 10 15 20

30

35

40

45

50

55

60

65

Figur 3.4: Exempel fr˚an simuleringsomg˚ang fyra. RF-data med 40 000 spridare, alla med samma styrka. Figuren är av samma typ som figur 2.6. En ultraljudspuls är sänd och mottagen i 520 mätpunkter varav 260 utanför spridarnas x-intervall. Samplingsfrekvensen är 50 MHz. Signalen som sändes var p˚a 5 MHz.

genererar ett v¨arde som ska vara sv˚art att f¨orutse utan kunskap om detta tillst˚and.

Men när MATLAB startas utg˚ar tydligen programmet fr˚an samma tillst˚and. Tillst˚andet m˚aste slumpas manuellt i skripten som körs och det m˚aste göras för rand respektive randn var för sig, eftersom de har var sitt tillst˚and. Detta kallas för att ”randomisera”

tillst˚and. Vanligen används klockan för att randomisera. I MATLAB görs detta med rand(’state’,sum(100*clock)) och precis samma argument för randn.

När detta var ˚atgärdat var det bara att göra en ny körning.

(39)

(40)

K ^APITEL 4 Resultat

Det h¨ar kapitlet beskriver resultaten, dels med skattningen av antalet spridare och resultaten p˚a v¨agen dit.

4.1 Syntetisk aperturavbildning

Eftersom parameterskattningen kan vara känslig för hur specklebilden ser ut är det viktigt att f˚a denna del att fungera s˚a bra som möjligt. Denna del var ocks˚a den mest arbetsamma när det gäller programmering och felsökning.

4.1.1 Simulering av ultraljudsutrustning

Det största problemet med att f˚a simuleringarna att fungera, var att f˚a Field att ge datasekvenser (signaler) som är lika l˚anga som börjar i samma tidpunkt. Detta löstes med dummy-spridare som ”spänner upp” omr˚adet som ska bli RF-data.

4.1.2 D ¨ampningskompensation

Det tog ett tag komma fram till hur dämpningen skulle mätas. Resultatet blev tv˚a metoder, en som är ganska beräkningskrävande och en snabb men inte lika noggrann. Hur

29

(41)

30 Resultat

Utan kompensation

−20 −10 0 10 20

20 22 24 26 28 30 32 34 36 38

Med kompensation

−20 −10 0 10 20

20 22 24 26 28 30 32 34 36 38

20 25 30 35 40

0.8 1 1.2 1.4 1.6 1.8 2 2.2

x 10⁻¹⁷ Utan kompensation

Intensitet (envelopp)

20 25 30 35 40

1200 1400 1600 1800 2000 2200 2400 2600 2800 3000 3200

Med kompensation

Intensitet (envelopp)

Figur 4.1: Utan och med dämpningskompensering. De nedre bilderna är en intensitetssum- mering längs varje rad hos de övre bilderna, s˚a att skillnaden framg˚ar tydligare. Skillnaden i medelintensitet som kan ses i de nedre figurerna har ingen betydelse för AR-skattningen.

metoderna fungerar beskrivs i avsnitt 2.1.3. Resultatet blev som i figur 4.1.

Dämpningskompensationen enligt den första metoden görs fr˚an mätningar av reflektionen i 64 g˚anger 64 positioner. I figur 4.2 visas hur reflektionerna s˚ag ut för de tv˚a typer av element som använts i simuleringarna.

När RF-data som blivit kompenserad med den snabba metoden blev envelopdetekte- rad uppträdde en tendens till ränder i bilden som inte finns i motsvarande grad p˚a

(42)

4.2. Parameterskattning och j¨amf¨orelser av

AR-parametrar 31

Första

−20 −10 0 10 20

20 22 24 26 28 30 32 34 36 38

Andra till femte

−10 −5 0 5 10

25 26 27 28 29 30 31 32 33 34

Figur 4.2: Reflektionerna fr˚an olika koordinater. Den vänstra bilden visar reflektionerna fr˚an den första simuleringen, den högra visar hur de övriga s˚ag ut.

icke kompenserade bilder och bilder kompenserad med den l˚angsamma metoden. AR- parametrarna fick sm˚a störningar som s˚ag ut som lite mer brus/variationer p˚a raderna ovanför parametrarna med större värden. Metoden har ocks˚a andra nackdelar, därför används den första metoden här.

4.1.3 Lobformning

Lobformningen inleddes med att göra laboration 4 i kursen medicinsk signalbehandling, sms046. Vissa problem med spatiell vikning tog ett tag att lösa eftersom det tog ett tag att först˚a att det var vikning.

4.2 Parameterskattning och j ¨amf ¨ orelser av AR-parametrar

Denna del var sv˚arast rent teoretiskt, speciellt att först˚a hur specklebilderna skulle analyseras d˚a det inte finns n˚agon teori bakom det, utan mest en massa gissningar. Det var till en början sv˚art att först˚a hur AR-skattningen skulle överföras fr˚an en dimension till flera.

(43)

32 Resultat

4.2.1 AR-skattning

Det största problemet här var att matrisen Φ^T tog mycket minne under beräkningen.

L¨osningen var dels den metoden som inte bygger upp alla rader hos Φ^T och Y , och dels

”Divide and Conquer”-metoden som bygger upp Φ^T del för del och utför beräkningarna p˚a dessa delar. Detta beskrivs i avsnitt 2.2.4. Den första metoden ger bra resultat även om bara en br˚akdel av raderna i matriserna byggs upp, den gör ocks˚a att beräkningen g˚ar mycket snabbare.

Ett annat problem med AR-skattningen som har med lobformningen att göra, var att de i första simuleringarna sökte elementet bara av x-koordinater där det fanns spridare.

Detta ledde till att eko-hyperblar blev trunkerade ganska nära dess vertex om de l˚ag nära kanten. D˚a blir det sämre fokus p˚a de respektive spridarna till dessa hyperblar, de blir utsmetade i x-led. Denna kanteffekt g˚ar att ana i figur 4.3.

Kanteffekten blir tydligare hos AR-parametrarna i figur 4.4, och i residualen (beskrivs i avsnitt 2.2.6), i figur 4.5, för de tidigare simuleringarna syns ränder, möjligen med hy- perbolisk form. Detta tyder starkt p˚a att dessa simuleringar inte är att betrakta som sta- tionärt brus. Men residualen i 4.5 för nyare simuleringar visar ocks˚a att AR-skattningen inte lyckas f˚anga upp all dynamik, eftersom residualen inte är vitt brus.

AR-parametermatriser med formatet av en 10 × 10 QP-matris visade sig vara fullt tillräckligt. Större matriser ger inget bättre resultat. QP-matrisen beskrivs i fig 2.12.

Figur 4.6 visar hur en realisering av de högra AR-parametrarna i figur 4.4 kan se ut, jämför detta med den högra specklebilden i figur 4.3. Notera att styrkan i figur 4.6 är mycket större än i figur 4.3. Detta har inget samband med parametrarna utan bestäms enbart av matningsbruset.

4.2.2 J ¨amf ¨ orelse av AR-parametrar

Enligt beskrivningen i kap. 2.3, gjordes en jämförelsematris med alla AR-parametrar ord- nade som radvektorer. Denna inneh˚aller parametrar fr˚an körning fyra och fem. (Körning fyra anses trots allt tillräckligt bra för analys trots misstaget med lika slumptalssekven- ser som beskrivs i avsnitt 3.5). Jämförelsematrisen byggdes upp fr˚an AR-parametrar i formatet av en 10 × 10 QP-matris, dvs det blir 99 enskilda parametrar p˚a varje rad i jämförelsematrisen (ett värde i matrisen används inte, vilket beskrivs i figur 2.12). P˚a denna gjordes en singulärvärdesuppdelning.

Modellering av ultraljudsbilder med tv ˚adimensionella AR-modeller

E X A M E N S A R B E T E

Modellering av ultraljudsbilder med tvådimensionella AR-modeller

Johan Svensson

Modellering av ultraljudsbilder med tv ˚adimensionella AR-modeller

Johan Svensson

Lule˚ a tekniska universitet Institutionen f¨ or Systemteknik

29 april 2005

S AMMANFATTNING

A BSTRACT

F ¨ ORORD

I NNEH ALL ˚

K APITEL 1 Inledning

1.1 Bakgrund

1.2 M ˚al

1.3 Avgr ¨ansning av arbetet

1.4 Uppl ¨agg

K APITEL 2 Teori

2.1 Syntetisk aperturavbildning

2.1.1 Datainsamling

2.1.2 Lobformning

2.1.3 D ¨ampningskompensation

2.2 Modellering av specklebild

2.2.1 Endimensionell AR-modell

2.2.2 Tv ˚adimensionell AR-modell

2.2.3 Tv ˚adimensionell stabilitet

2.2.4 Endimensionell parameterskattning

2.2.5 Tv ˚adimensionell parameterskattning

2.2.6 Residual och val av ordning

2.3 J ¨amf ¨ orelse av AR-parametrar

K APITEL 3 Simuleringar och andra

ber ¨akningar

3.1 F ¨ orsta k ¨ orningen

3.2 Andra k ¨ orningen

3.3 Tredje k ¨ orningen

3.4 Fj ¨arde k ¨ orningen

3.5 Femte k ¨ orningen

K APITEL 4 Resultat

4.1 Syntetisk aperturavbildning

4.1.1 Simulering av ultraljudsutrustning

4.1.2 D ¨ampningskompensation

4.1.3 Lobformning

4.2 Parameterskattning och j ¨amf ¨ orelser av AR-parametrar

4.2.1 AR-skattning

4.2.2 J ¨amf ¨ orelse av AR-parametrar

I ^NNEH ALL ˚

K ^APITEL 2 Teori

K ^APITEL 3 Simuleringar och andra

K ^APITEL 4 Resultat