Het gebruik van de Gabor transformatie bij de analyse van het geluid van carillonklokken: een inleiding in de Gabor transformatie

Het gebruik van de Gabor transformatie bij de analyse van het

geluid van carillonklokken

Citation for published version (APA):

Cornelissen, A. H. P. J. (1993). Het gebruik van de Gabor transformatie bij de analyse van het geluid van

carillonklokken: een inleiding in de Gabor transformatie. (DCT rapporten; Vol. 1993.039). Technische Universiteit Eindhoven.

Document status and date: Gepubliceerd: 01/01/1993 Document Version:

Uitgevers PDF, ook bekend als Version of Record Please check the document version of this publication:

• A submitted manuscript is the version of the article upon submission and before peer-review. There can be important differences between the submitted version and the official published version of record. People interested in the research are advised to contact the author for the final version of the publication, or visit the DOI to the publisher's website.

• The final author version and the galley proof are versions of the publication after peer review.

• The final published version features the final layout of the paper including the volume, issue and page numbers.

Link to publication

General rights

Copyright and moral rights for the publications made accessible in the public portal are retained by the authors and/or other copyright owners and it is a condition of accessing publications that users recognise and abide by the legal requirements associated with these rights. • Users may download and print one copy of any publication from the public portal for the purpose of private study or research. • You may not further distribute the material or use it for any profit-making activity or commercial gain

• You may freely distribute the URL identifying the publication in the public portal.

If the publication is distributed under the terms of Article 25fa of the Dutch Copyright Act, indicated by the “Taverne” license above, please follow below link for the End User Agreement:

www.tue.nl/taverne

Take down policy

If you believe that this document breaches copyright please contact us at:

openaccess@tue.nl

Door: A.€IPJ. &i.~eEssen id. nr.: 300692

WFW-rapport 93.065

Stage begel&derc: ir. A. de Mlaker %r.

R.

Kodde

Inhoudsopgave

1 Inleiding 3

2 Beschrijving van de Gabor transformatie 4

2.1 Continue gabor transformatie

. . .

4

2.2 Discrete Gabor transformatie

. . .

4

2.3 Discrete Gabor transformatie voor periodieke signalen

. . . .

6

2.4 De bi-orthogonale functie

. . .

7

2.5 Oversampelen

. . .

8

2.6 De overeenkomsten tussen Gabor en Fouiier transfainziiie

. .

10

3 Implementaties 12 3.1 De bi-orthogonale functie

. . .

12

3.2 Gabor transformatie

. . .

13

3.3 Inverse Gabor transformatie

. . .

13

3.4 Oversampelen

. . .

13

3.5 Fourier Transformatie

. . .

14

4 Testen 14 4.1 De bi-orthogonale functie

. . .

14

4.2 Gabor transformatie

. . .

16

4.3 Inverse Gabor transformatie

. . .

19

4.4 Oversampelen

. . .

19

4.5 Gabor met Fourier Algoritme

. . .

21

4.7 Kloksignaal

. . .

22

. . .

4.6 Ruis 21 5 Conclusies 25 6 Aanbevelingen 26 A Uitwerking H-matrix 26 B Gabor programma’s 27 B . l MaakGamC

. . .

27 B.2 GaborC

. . .

29 B.3 GabInvC

. . .

30 B.4 MaakGamO

. . .

32 B.5 GaborO

. . .

35

B.6 GabInvO

. . .

36 B.7 Gabor

. . .

37 C Test programma’s 39 41.1 RiOrthog i i i i i i i i i i i : : : i : i i : i i 39 C.2 RechtHk

. . .

40 C.3 GaussVen

. . .

41 C.4 TestAOO

. . .

42 C.5 GabRecht

. . .

43 C.6 GabGauss

. . .

44 C.7 Vensters

. . .

45 C.8 BandBr

. . .

48 C.9 InvAOO

. . .

49 C.10 TestInvG

. . .

50 C.11GeEjk

. . .

51 C.12 BiOrOver

. . .

52 C.13 GamQver

. . .

54 C.14 OverSam

. . .

55 C.15FFTGab

. . .

56 C.16Ruis

. . .

57 C.17 Invloed

. . .

58 C.18 KlokSim

. . .

60 C.19 MlokGel

. . .

63 D Hulp programma’s 64 D.l Effect

. . .

64 D.2 FlipIt

. . .

65 D.3 Gauss

. . .

66 D.4 NaNs

. . .

66 D.5 Normaal

. . .

67 D.6 PlotGab

. . .

67 D.7 SinSweep

. . .

68

Samenvatting

De Gabor transformatie zou een mogelijkheid kunnen zijn om het inslinger gedrag van carillonklokken te kunnen analyseren. Dit is het gedrag van een carïïonMok direct na de aansiag met de kiepei, de klok is dan nog niet overaj in trilling. Hieraan is het akoesitisch inslingergedrag gekoppeld. Deze Gabor transformatie wordt beschreven en geïmplementeerd in een aantal MATLAB programma’s.

Met deze programma’s worden een aantal tests uitgevoerd, waaruit blijkt dat de Gabor transformatie iets meer rekenwerk vereist dan de Short Time Fourier Transformation, hierbij wordt het signaal in verschillende stukken opgedeeld, waarna van deze stukken de Fourier transformatie bepaald wordt. Ook is het zo dat de Gabor transformatie een iets betere resolutie geeft dan deze STFT.

Het blijkt dat de Gabor traasformatie met de door mij gekozen basis- functies te weinig resolutie heeft om de geluidssignalen voldoende nauwkeu- rig te bekijken. Een betere resolutie zou verkregen kunnen worden door het gebruik van andere basisfuncties of gebruik te maken van een soort “zoom- Gabor”.

1

Inleiding

In 1991 is door Van Klinken [5] onderzoek gedaan naar het inslingeren van carillonklokken. Bij dit onderzoek zijn metingen gedaan van het trillings- gedrag van verschillende klokken. Hierbij werd geprobeerd door gebruik t e maken van bandfilters de verschillende eigenfrequenties afzonderlijk zicht- baar t e maken. Dit leverde problemen op omdat de aanlooptijd van het filter van dezelfde orde is als de inslingertijd van de afzonderlijke frequentie van de klok.

Van Klingen beveelde aan om met behulp van Short Time Fourier Trans- formation, eventueel met overlap, naar het signaal in zowel tijd- als frequen- tiedomein t e bekijken. De

STFT

bleek niet voldoende resolutie t e hebben om het inslinger gedrag t e laten zien.

Almgren [l] beweert dat de Gabor transformatie [4] een veel betere reso- lutie (zowel voor tijd als frequentie) geeft dan de STFT. In mijn stage wordt de Gabor transformatie in MATLAB geïmplementeerd om t e onderzoeken of het inslingergedrag wel te analyseren is met de Gabor transformatie.

Eerst wordt een korte beschrij- Het verslag is als volgt opgebouwd.

ving gegeven van de Gabor transformatie. Daarna wordt beschreven hoe deze t:a-r,shrmatie in MATLATY geïmplementeerd wordt. Waarna deze pro- grarnaa's op verschillende manieren getest worden. Tot slot wordt ook nog een gemeten geluid geanalyseerd.

2

Beschrijving van

de

Gabor transformatie

2.1 Continue

gabor

transformatie

De Gabor transformatie (zie [ 4 ] ) is een methode voor tijd-frequentie analyse. Hierbij wordt een signaal opgebouwd verondersteld uit de sommatie van in tijd en frequentie verschoven basisfuncties:

z ( t ) = a,,,h(t - mT)einnt (1)

m n

Waarbij de volgende relatie geldt:

R T = 2T (2)

Dit betekent dat het signaal wordt bescheven met een discrete ves.zarne- Zing van basis functies die verschoven zijn over discrete afstanden mT en gemoduleerd met discrete frequenties

no.

Bastiaans [2] heeft bewezen dat de Gabor coëfficienten eenvoudig bepaald kunnen worden als we een functie y(t) kunnen vinden die bi-orthogonaal is

ten opzichte van de verzameling basis functies volgens:

J

h(t)y*(t - rnT)e-in"tdt = SmSn

am,n =

s

z ( t) y *( t

- mT)e-in"tdt

(3) Hierbij staat de asterix (*) voor complex conjungeren en SI, voor de Kron- ecker delta (SO = 1 en SI, = O voor k

#

O).

De Gabor coëfficienten kunnen nu met de volgende formule bepaalt wor- den:

( 4 ) 2.2 Discrete

Gabor

transformatie

Het in de tijd continue signaal wordt gediscretiseerd met een intervallengte A. Voor discrete signalen luidt de Gabor transformatie (conform vergelij- king i) als volgt:

De Gabor coëfficienten kunnen nu bepaald worden met:

(6) am,n = x[[t]y*[t - r n ~ l e - ~ ~ "

t

Nu worden de discrete vergelijkingen zo omgeschreven dat slechts gehele getallen tussen de rechte haken staan. Dit gebeurt voor vergelijking 5 door in deze vergelijking de tijd t te substitueren door:

Dan volgt de volgende vergelijking:

x [ A k ] = a,,,h[Ak - mT]einQAk

m n

Waarna de termen tussen de rechte hakken door

A

gedeeld worden, zodat overblijft :

Voor de eenvoud wordt nu de constante W gedefinieerd:

(10)

w

= ,iRA

Als deze constante ingevuld wordt in vergelijking 9 ontstaat tenslotte:

mT

.[I;:] = am,nh -

á]

Wrik

m n

Dit kan ook voor vergelijking 6 gedaan worden. Substitutie van verge- lijking 7 levert:

(12)

am+ = ~[Ab]y*[Ak - mT]e-in"Ak

k

Tussen de rechte haken delen door A leidt tot:

Figuur 1: eindige signalen worden als periodiek beschouwd

2.3

Om de hiervoor besproken vergelijkingen in een computer programma te kunnen implementeren worden nu slechts periodieke of eindige signalen be- schouwd. Eindige signalen worden periodiek beschouwd zie figuur l. Voor een signaal met een periode tijd t p zijn nu

P

sampels nodig:

Discrete Gabor transformatie voor periodieke signalen

t p = P A (15)

Nu worden een aantal nieuwe parameters gedefinieerd. Te weten: Het aantal tijdpunten ( M ) en het aantal frequentiepunten ( N ) . De parameters

fi

en T

kunnen dan als volgt gesubstitueerd worden:

Met vergelijking 2 volgt hieruit:

T N = -

A

Zodat vergelijking 11 als volgt kan worden geschreven:

en vergelijking 14 op de volgende manier:

Een alternatieve defenitie voor W wordt dan:

(21)

127r

W = e N

Door het nemen van een periodiek signaal, worden de coëfficienten ook pe- riodiek:

am,n = am+M,n+N (22)

De Gabor transformatie wordt dan:

M - 1 N-1

$4

= i [ k - mN]Wnk

m=O n=O Zie ook [8]. Waarbij

h

gedefinieerd is als:

h [ k ] = h[k

+

IP] (24)

1

Voor de Gabor coëfficienten geldt dan:

P-1

= .[k]?*[k - mN]W-"k

k=l

2.4

De

bi-orthogonale functie

Het probleem voor het vinden van de Gabor coëfficienten is omgezet in het probleem voor het vinden van de bi-orthogonale functie y. Deze bi- orthogonaliteitseis luidt voor discrete periodieke functies:

P-1

(hik

+

~ N I W - ~ ~ ) ?*[/cl = S,S, Wat ook geschreven kan worden als:

(P[k

+

~ N I W ~ ~ )

?*[-I

=

S

,

S

,

v

O

5

m

5

M - i A O

5

n

5

N - 1 (26)

v

O

5

m

5

M - i

A

O

5

n

5

N - i (27) k=O P-1 k=O

Deze vergelijking kan als volgt in matrix vorm geschreven worden:

"

=

H * y = g - (28)

E e r i n is de matrix E! eeri.

P

x

P matrix en wiin Y en I/ kolom-m-en met lengte

B.

De matrix

W

is op de volgende manier onder t e verdelen in M x M

submatrices: --J-- L --- - b [ Z N ] O **e

-wo

wo

wo

...

wo

-

wo

w1

w2

...

WN-I

wo

w2

w4

**a W N - ~

1:

h * [ Z N t ' ] ~

wo

WN-1 WN-2

...

WO O k [ l N

+

N - 11 -

Voor deze N

x

N submatrices geldt:

2.5 Oversampelen

Bij de hiervoor uitgevoerde Gabor transformatie worden evenveel onafhan- kelijke Gabor coëfficienten gevonden als er sampels genomen zijn. Door de restrictie R T = 27r t e laten vallen is het mogelijk om meer coëfficienten te berekenen. Deze coëfficienten zijn dan niet onafhankelijk van elkaar. We gaan bij oversampelen uit van het geval:

Dit betekent:

Nu worden twee nieuwe variabelen gedefinieerd:

P

N i =

P

Mi = - N

De Gabor representatie luidt dan als volgt:

M-1 N-I

(33)

(34)

m=O n=O

De Gabor coëfficienten kunnen dan op de volgende manier bepaald worden:

P-1

= 2[k]Y*[k - mN1]W-nk (36)

k=O

De bi-orthogonaliteits voorwaarde waaraan

Y

moet voldoen luidt als volgt:

P

-6

6 ~ O < m < M 1 - 1 A O < n < N I - 1 M N

(h[k

- m N ] W F n k )

7 * [ k ]

= P-1 k=O (37) W1 = e N i (38) Waarbij geldt: i2lr - Vergelijking 37 kan worden uitgeschreven tot:

P

P-1

(hast[k

+

m N ] W T k ) Y [ k ] = =SmSnVO

<

m

<

M l - 1 A O

<

n

5

N I - i

k=O

(39) Ook deze vergelijking is als matrix vergelijking t e schrijven:

In dit geval is H een MiN1 x P matrix die onderverdeeld kan worden in

Ml x M sub-matrices met afmetingen Ni x N l . De kolom y met lengte P

bestaat ook nu weer uit de waarden van de bi-orthogonale functie van h. De

kolom (lengte M I N I ) bestaat, op het eerste element na, ook weer geheel uit nullen. Het eerste element is in tegenstelling bij het kritisch sampelen niet 1 , maar

m.

P

Deze matrix vergelijking is onderbepaald, er zijn dus meerdere oplos- singen mogelijk. De door Weder en Raz [8] gesuggereerde gegeneraliseerde inverse oplossing:

y = H T ( M . H ) T -1

-9

geeft een y met minimum energie. Voor het bewijs zie Qian [7] paragraaf 2.

Ook

is het mogelijk om een aantal (maximaal P-IMlNl) waarden van y op t e leggen. Als aan de rand van het venster een aantal nullen opgelegd worden, zal het venster een in de tijd nauwkeuriger spectrum geven. Indien niet het maximale aantal waarden wordt opgelegd zal ook weer een onderbepaald stelsel ontstaan, wat weer met de gegeneraliseerde inverse oplossen bepaald kan worden.

2.6 De overeenkomsten tussen Gabor en Fourier transfor- matie

De Fourier transformatie (zie bijvoorbeeld [6] paragraaf 2.3.2):

-

x

= F ( g ) (42)

of per element uitgeschreven:

P-1

X [ n ] = Z[k]WF"k

k=O

waarbij WF is gedefinieerd als:

lijkt erg veel op de manier bepaald worden:

i2T

W F = e 7

(43)

(44) waarop de Gabor coëfficienten van een m-rij

k=O

Deze vergelijking kan als het product van ~ [ k ] met y*[- - mNl] worden samengevoegd tot:

geschreven worden als:

Y[kI

=

.[-Ir*[-

- m&I (46)

P - 1

a,[n] = y[k]w-"k ₍₄₇₎

Om

W

en

WF

gelijk aan elkaar t e krijgen wordt de sommatie opgesplitst: M-1 l(N1+1)-1

Of:

M-1 Ni -1

I=O k=O

a&] =

Y[-

+

~NlIW-"]" (49)

Tot slot worden de sommatie tekens gewisseld:

Zodat een rn-rij van de coëfficienten matrix als volgt bepaald kan worden:

Waarbij y gedefinieerd is als:

-1

M-1

k 0

Als N I een macht van twee is kan voor de Fourier transformatie het Fast Fourier Transform algoritme gebruikt worden. De benodigde rekenoperaties (en daardoor de ook de benodigde tijd) voor het bepalen van de Gabor coëfficienten nemen daardoor aanzienlijk af.

Figuur 2 laat een grafische weergave van het nu afgeleide FFT-algoritme voor de Gabor transformatie zien. Linksboven staat het t e transformeren signaal, rechtsboven de reeds bepaalde bi-orthogonale functie. Deze func- ties worden met elkaar vermenigvuldigd, wat de functie linksonder oplevert. Deze functie wordt in een aantal

(6

= N I , in dit geval 4) stukken gehakt, waarna deze stukken bij elkaar opgeteld worden. Het resultaat is rechts- onder te zien. Als op deze functie de F F T wordt uitgevoerd, ontstaan de Gabor coëfficienten voor deze tijdstap. Door de bi-orthogonale functie t e verschuiven en de hiervoor beschreven procedure t e herhalen kunnen de Ga- bor coëfficienten op de andere tijdstippen bepaald worden.

"O 5 10 15 u)

Figuur 2: Grafische weergave van het FFT-algoritme voor de Gabor trans-

formatie -

3

Implementaties

3.1

De bi-orthogonale functie

De voorgaande methode om de bi-orthogonale functie t e bepalen wordt nu geïmplementeerd in een matlabprogramma. Dit matlabprogramma staat af- gedrukt in bijlage B.1. Als parameters voor deze functie moeten worden opgegeven: de basisfunctie (h), het aantal tijdstappen (M) en het aantal fre- quentiestappen

(N).

De uitvoer van deze functie is de bi-orthogonale functie van h (gam).

Na het uitrekenen van de bi-orthogonale functie y verschijnt een getal op het beeldscherm. Dit getal geeft de maximale waarde aan van het imaginaire deel. Dit getal zou gelijk aan nul moeten zijn. Dit is echter niet het geval, het getal wat verschijnt ligt (op mijn computer) in de ordegrootte van dit getal wordt bepaald door de nauwkeurigheid van de computer.

3.2

Gabor transformatie

De gabor transformatie wordt geïmplementeerd in het programma GaborC, t e vinden in bijlage B.2. De parameters voor deze functie zijn: het te ana- lyseren signaal (x) i d e bi-orthogonale functie gamma, het aantal tijdstappen

(M), het aantal frequentiestappen(N) en de tijd tussen twee sampels (dt). De uitvoer van deze functie is: de Gabor coëfficienten (a), een tijdvector (t) en een frequentie vector (f).

Het programma rekent de Gabor coëfficienten een voor een uit. Als het programma een coëfficient heeft uitgerekend verschijnen van deze coëfficient de m en n op het scherm met de waarde van deze coëfficient

.

In het programma wordt een index f[i] aangegeven met f

Ci+ll

omdat matlab zijn indices laat beginnen met 1 in tegenstelling tot de algemene literatuur waar de indices met O beginnen. Dit betekent dus dat wordt aangegeven met a Cm+ 1

,

n+

ll.

3.3

Inverse Gabor transformatie

De inverse Gabor transformatie wordt geïmplementeerd in het programma GabInvC, t e vinden in bijlage B.3. De parameters voor deze functie zijn: de coëfficienten matrix (a) en de basis functie (h). De uitvoer van deze functie is het gereconstrueerde signaal (x).

Het programma rekent de functie punt voor punt uit. Als het programma een punt heeft uitgerekend laat het programma het nummer van dat punt zien plus de waarde die het programma daarbij heeft uitgerekend. Na het uitrekenen van de gehele functie verschijnt de maximale imaginaire waarde in beeld. Omdat er slechts met reële functies gewerkt wordt zou dit getal

nul moeten zijn. Door de rekenonnauwkeurigheid van mijn computer ligt

deze waarde in de ordegrootte van Om de verdere verwerking van de resultaten gemakkelijker te maken wordt slechts het reële gedeelte van het resultaat naar buiten gebracht.

3.4

Oversampelen

Door de programma’s voor het kritisch sampelen aan t e passen zoals in pa- ragraaf 2.5 beschreven staat, ontstaan de programma’s: MaakGamO, GaborO en GabinvO, deze programma’s staan afgedrukt in respectievelijk bijlage B.4, B.5 en B.6.

De invoer en uitvoer van deze programma’s is hetzelfde als bij kritisch sampelen, met uitzondering van het programma MaakGamO. Bij dit pro-

gramma kan (het is dus niet noodzakelijk) een rij gamin mee gegeven worden in deze vector kunnen de opgelegde waarden van y geplaatst worden, als een waarde door de computer berekend moet worden moet daar de MATLAB constante NaN (Not a number) worden ingevuld.

3.5

Fourier Transformatie

Het programma GaborO wordt op de manier die in paragraaf 2.6 beschreven staat aangepast, zodat gebruik gemaakt kan worden van de Fourier transfor- matie. Het aangepaste programma heet Gabor en is t e vinden in bijlage B.7. De in- en uitvoer van Gabor zijn hetzelfde als bij GaborC. Omdat de coëffienten nu niet meer stuk voor stuk, maar per rij bepaald worden is het niet meer nuttig om deze coëfficienten stuk voor stuk op het scherm t e laten zien, in plaats daarvan wordt het nummer van de uitgerekende rij op het beeldscherm gezet.

De inverse Gabor transformatie is niet zo eenvoudig om te schrijven tot een algortme waarbij gebruik gemaakt kan worden van de FFT. Dit komt doordat ook in de basisfunctie de index

k

zit.

De bepaling van de bi-ortogonale functie zou sneller bepaald kunnen worden door gebruik te maken van de speciale structuur van de

H

matrix. Als vaak van de Gabor transformatie gebruik gemaakt zal gaan worden, is het nuttig om een verzameling van bi-orthogonale functies op te slaan die telkens opnieuw gebruikt kunnen worden.

4

Testen

De testen worden uitgevoerd op een IBM kloon met een 80386 C P U en 80387 co-processor draaiend op 33 Mhz. Ik maak gebruik van 386-MATLAB7 versie 3.5g. De snelheid van mijn computer is (volgens de benchmarks van 386-MATLAB) 26,0022 keer zo snel als een PC/XT, wat overeenkomt met 461 Kflops/s.

4.1

De bi-orthogonale functie

Met het programma MaakGamC worden een aantal tests uitgevoerd. In de eerste test wordt gekeken of als van een functie de bi-orthogonale functie wordt uitgerekend de bi-orthogonale functie van deze bi-orthogonale functie weer de oorspronkelijke functie is. Dit wordt gedaan met het programma: BiOrthog dat afgedrukt staat in bijlage C.1. Dit programma genereert

een aantal willekeurige functies, bepaalt hier de bi-orthogonale functie van, waarna van deze bi-orthogonale functie weer de bi-orthogonale functie wordt bepaald. Hiervan wordt het absolute verschil met de oorspronkelijke functie bepaald. Tot slot wordt het gemiddelde van al deze verschillen bepaald. Op mijn computer gaf dat een verschil in de orde grootte van lû-15. Dit is terug te voeren op de nauwkeurigheid van de computer.

De tweede test is de bi-orthogonale functie bepalen van een rechthoe- kig venster. Dit gebeurt met het programma RechtHk zie bijlage C.2. De breedte van het rechthoekige venster is gelijk aan de stapgrootte. Dit bete- kent dat de bi-orthogonale functie van dit venster gelijk moet zijn aan het venster zelf. Dit blijkt (behoudens enkele rekenonnauwkeurigheiden) ook te

kloppen.

Als laatste wordt het door Gabor gesuggereerde Gauss venster getest met het programma GaussVen, t e vinden in bijlage (2.3. Hieruit blijkt dat de gevonden bi-orthogonale functie visueel lijkt op de functie die door Bastiaans [3] analytisch is afgeleid, zie figuur 3. Volgens Wexler en Raz [8] komt dit

Venster met bi-orthogde ñmctie

o * 6 *

I

O 10 20 30 40 50 60 70

Figuur 3: Een Gauss venster (onderbroken lijnen) met zijn bi-orthogonale functie (ononderbroken lijnen)

4.2

Gabor

transformatie

De eerste test die gedaan wordt is kijken of de Gabor transformatie van de basisfunctie voor elke Gabor coëfficient O oplevert, behalve voor a o , ~ waar de s m d e 1 meet zijc. Deze test wordt iiitgevoerd met het programma TestAOO, afgedrukt in bijlage C.4. Hieruit blijkt inderdaad dat alle co-

efficienten gelijk aan nul zijn, behalve ao,o dat de waarde l heeft. Als iets nauwkeuriger gekeken wordt zien we dat de maximale absolute waarde waar- mee de berekende coëfficienten van de theoretische afwijken in de orde van ligt. Dit is weer terug te voeren op de reken onnauwkeurigheid van mijn computer.

De tweede test is de transformatie met een rechthoekig venster. Deze transformatie komt overeen met de short- time fourier transformatie (STFT).

In het programma GabRecht, afgedrukt in bijlage C.5, wordt de transfor- matie uitgevoerd met een sinus signaal dat precies in het venster past. De verwachte resultaten (Een rechte lijn ter hoogte van de frequentie van de sinus) verschijnen op het scherm. Ook is de de vouwfrequentie goed t e zien.

In dit programma worden verschillende signalen met steeds hogere frequen- tie getransformeerd. Hierdoor is het effect van aliassing goed te zien. Met het programma GabGauss wordt een sinus signaal getransformeerd dat niet precies in het venster past. Eerst wordt dat gedaan met een rechthoekige basisfunctie, de signaallek is dan duidelijk zichtbaar. Als in plaats van een rechthoekige basisfunctie, de functie van Gauss als basisfunctie genomen wordt komen de Gabor coëfficienten meer rond de frequentie van het signaal t e liggen.

Bij de volgende test wordt de sine sweep met een aantal verschillende basisfuncties getransformeerd. Dit gebeurt met het programma vensters, afgedrukt in bijlage ap:vensters. De verschillende basisfuncties waarmee de Gabor transformatie uitgevoerd worden zijn:

e Rechthoek e Bartlett

e Gauss

e Hanning

e Hamming

Alle basisfuncties hebben een effectieve breedte (gedefinieerd als de som van alle punten gedeeld door het maximum van de functie) die 16 zo dicht moge- lijk benadert. Exact 16 is niet in alle gevallen mogelijk omdat de basisfunc- tie een discreet aantal punten moet bevatten. Ook worden de bi-orthgonale functie van deze functies als basisfunctie gebruikt. De hoogtelijnen van

de absolute waarden van de Gabor coëfficienten zijn afgebeeld in figuur 4. Omdat de Gabor coëfficienten bij reëele functies symetrisch zijn is alleen het

Tul lil

r

c

Figuur 4: De hoogtelijnen van de absolute waarden van de Gabor coëfficien- ten bij verschillende basisfuncties

onderste gedeelte afgebeeld.

Bij de rechthoekige basisfunctie (STFT) wordt een redelijke represen- tatie verkregen. De band die op het beeldscherm verschijnt is recht. Bij de “Gewone basisfuncties” (Bartlett, Gauss, Hamming en Hanning) zijn de hoogtelijnen van de coëfficienten matrix niet meer recht, de coëfficienten op de “diagonaal” zijn afwisselend hoog en laag. Het resultaat van de Bartlett functie is zelfs zeer slecht. Bij de zogenaamde bi-orthogonale basisfuncties is het resultaat iets beter dan bij de rechthoekige basisfunctie, de hoogtelijnen van de coëfficienten matrix zijn weer recht, en de band die op het scherm verschijnt is kleiner dan bij de rechthoekige basisfunctie, vooral bij de ba- sisfunctie die bi-orthogonaal ten opzichte van Gauss is, en in iets mindere

mate bij de basisfuncties die bi-orthogonaal zijn ten opzichte van Hamming en Hanning.

Het effect van de bandbreedte (gedefinieerd als de som van alle punten gedeeld door de maximale waarde van het venster) wordt getest met het programma B a n d B r (zie bijlage C.8). Hierbij wordt Gabor transformatie uitgevoerd met de bi-orthogonale functie van de functie van Gauss. Hier- bij wordt de effectieve breedte van de bi-orthogononale functie steeds opge- voerd. Hierbij zien we dan dat de band die op het beeldscherm smaller wordt naarmate de effectieve bandbreedte toeneemt. Als de effectieve bandbreedte 32 (De bandbreedte bij STFT bij deze afmetingen van de coëfficienten ma- trix) overschrijdt, neemt de band op het scherm weer langzaam toe. Om

dit duidelijker t e maken zijn in figuur 5 de absolute waarden van de Gabor coëfficienten van de achtste rij in een grafiek uitgezet. Door dit “spelen”

Gabor coefficienten bij veischiuende bandbreedtes

Frequentie

Figuur 5: De absolute waarde van de Gabor coëfficienten van de achtste rij in de coëfficienten matrix bij verschillende effectieve breedtes (van boven naar beneden: 64,48,40, 32, 24 en 8)

met de effectieve bandbreedte laat Almgren [i] de Gabor transformatie er een stuk beter uit zien dan de STFT. Almgren legt een Hamming venster over het rechthoekige gedeelte, waardoor de effectieve bandbreedte afneemt

,

zodat de resolutie slechter wordt.

4.3 Inverse Gabor transformatie

De eerste test die wordt uitgevoerd is Ejken of e a c~ëfE:laenten matrix waarvan alle waarden O zijn, behalve ao,o dat gelijk aan 1 is bij inverse Gabor transformatie de basisfunctie opleverd. Dit wordt getest met het programma InvAOO dat staat afgedrukt in bijlage C.9. Dit blijkt exact zo het geval te zijn.

De tweede test wordt uitgevoerd met het programma TestInvG (zie bij- lage C.10. Met dit programma wordt een willekeurig signaal gegenereerd. Daarna wordt met een willekeurige basisfunctie de Gabor transformatie uit- gevoerd. Hierna wordt met de gevonden Gabor coëfficienten het signaal weer gereconstrueerd. Dit gereconstrueerde signaal wordt vergeleken met het oorspronkelijke signaal. Het verschil tussen deze signalen zou theore- tisch nul moeten zijn, maar is door afrondfouten van mijn computer van de orde

4.4

Oversampelen

Het programma Gelijk, te vinden in bijlage C.11, test of met de pro- gramma’s voor oversampelen dezelfde resultaten verkregen worden als met de programma’s voor kritisch sampelen. Voor het bepalen van de bi- orthogonale functie ligt het verschil in de orde van het bepalen van de Gabor coëfficienten en het reconstrueren van het signaal levert voor beide programma’s exact hetzelfde resultaat.

Het programma M a a k G a m O wordt met het programma BiOrOver, af- gedrukt in bijlage C.12, op soortgelijke wijze getest als MaakGamC met BiOrthog. Het absolute verschil tussen de oorspronkelijke en de bi- orthogonale functie van de bi-orthogonale functie hiervan ligt in de orde van 10-l~.

Met het programma GamOver, dat t e vinden is in bijlage C.13, wordt het bepalen van een bi-orthogonale functie waarbij aan de rand van deze functie nullen opgedwongen worden. De resultaten van dit programma zijn t e zien in figuur 6. In deze figuur is te zien dat de bi-orthogonale functies grilliger worden als het aantal gedwongen nullen toeneemt, op de grafiek is duidelijk t e zien dat de energie dan ook toeneemt.

Tot slot wordt het effect van oversampelen getest. Dit gebeurt met het programma OverSam, t e vinden in bijlage C.14. Hierin wordt van een sine

Aamal5d*.mpm*uipuntui

Figuur 6: De verschillende gammafuncties (links) en de energie van deze verschillende functies (rechts)

sweep die bestaat uit 128 punten de Gabor coëfficienten bepaald als vier keer wordt oversampeld (Twee keer in de tijd en twee keer in de frequentie). Deze coëfficienten worden vergeleken met de coëfficienten die onstaan als een sine sweep van 512 punten met kritisch sampelen getransformeerd wordt. Uit

deze test blijkt dat de coëfficienten die met oversampelen bepaald zijn een minder scherp beeld opleveren, zie ook figuur 7.

Figuur 7: Het verschil in Gabor coëfficienten bij kritisch samepelen (links) en oversampelen (rechts)

4.5

Gabor

met

Fourier

Algoritme

Met het programma FFTGab, afgedrukt in bijlage (3.15 worden GaborO (zon- der FFT-algoritme) en Gabor (met FFT-algoritme) met elkaar vergeleken, wa,t hetreft benodigde tijd en rekenoperaties (Flops). De resultaten van deze vergelijking zijn afgedrukt in tabel 1. Het maximale absolute verschil tussen

I

Tijd

Is1

Flops

I

, -

_ _

1

GaborO

I

17,63 154012

1

I

Gabor

I

0,33 2301

I

Tabel 1: Vergelijking tussen GaborO en Gabor

de coëfficienten die met en zonder FFT-algoritme bepaald zijn is van de orde dit is t e wijten aan afrond fouten tijdens het rekenproces.

4.6

Ruis

In laatste serie testen die worden uitgevoerd wordt gekeken naar de invloed die ruis heeft op de Gabor transformatie.

In de eerste test, die wordt uitgevoerd met het programma Ruis (zie bijlage C.16). Met dit programma wordt het effect bekeken dat een steeds sterker wordend ruissignaal op de Gabor transformatie heeft. Dit levert de resultaten figuur 8.

Bij de tweede test wordt gekeken naar de invloed van het venster op het effect van ruis. Om dit t e bekijken wordt met het programma Invloed, afgedrukt in bijlage C.17. In dit programma worden met verschillende basis- functies de Gabor coëfficienten bepaald van een ongestoord signaal, hiervan worden de Gabor coëfficienten afgetrokken die bepaald zijn van het signaal met ruis. Dit verschil laat bij iedere basisfunctie nog kleine toppen zien, dat betekent dat ruis de extremen van de Gabor coëfficienten wat afzwakt. Hierna worden de verschillen tussen de Gabor coëfficienten van het signaal met en zonder ruis bij de verschillende basisfuncties vergeleken met het ver- schil tussen de Gabor coëfficienten van het signaal met en zonder ruis bij een rechthoekig venster. Hieruit blijkt dat er weinig verschil is tussen de verschillende basisfuncties met betrekking tot het effect van ruis. Als heel nauwkeurig wordt gekeken blijkt dat, bij het gebruik van basisfuncties die bi-orthogonaal zijn ten opzichte van de functies van Gauss en Hanning, de

Figuur 8: De invloed van het ruisniveau op de Gabor transformatie toppen van de Gabor coëfficienten iets minder afgezwakt worden dan bij de andere basisfuncties.

4.7

Kloksignaal

In de hiervoor beschreven testvoorbeelden is gebruik gemaakt van de sine sweep. Het voordeel van deze functie is dat het ideale resultaat (een diago- nale lijn) en de afwijkingen daarvan, gemakkelijk bepaald kunnen worden.

Bij carillonklokken is het waarschijnlijker dat frequenties opkomen en af- nemen. In het programma KlokSim (zie bijlage (3.18) wordt dit gesimuleerd. Hierbij wordt een signaal opgebouwd uit drie signalen met verschillende fre- quenties waarvan in de tijd de amplidudes gemoduleerd worden. Bij deze test worden de coëfficienten bepaald met:

o Short Time Fourier Transformation

o STFT met Hamming

o Gabor transformatie

De basisfunctie van de Gabor transformatie is de bi-orthogonale functie van de functie van Gauss. Bij de STFT leidde dit tot een “ruwe” verzameling van coëfficienten. Als het Hamming venster over de rechthoek wordt ge-

zet, wordt de mesh een stuk gladder, maar dit gaat wel ten koste van de resolutie. De Gabor transformatie geeft tenslotte een “gladde” verzameling coëfficienten met ongeveer dezelfde resolutie als de

STFT.

Tot slot wordt een signaal dat door Van Klinken gemeten is geanalyseerd. De hoogtefijnen van de absolute waardes van de Gabor coZ%cieriten zijn afgebeeld in figuur 9. Deze coëfficienten zijn bepaald met als bi-orthogoaak

Ge115

I

I

~~

“O 0.01 0.M 0.03 0.04 0.05 0.06 0.07

Tijd [SI

Figuur 9: De hoogtelijnen van de absolute waarden van de Gabor coëfficien- ten van het signaal “Ge115”

functie de Gausse functie met een effectieve breedte van 256,

M

= 8 en

N = 256. We zien dan dat het signaal slechts frequenties tot ongeveer 30000 Hz bevat. Om deze reden worden slechts de onderste 16 coëfficienten (tot 2929,7 Hz) afgebeeld. In deze afbeelding zijn de eigenfrequenties 585,9, 976,6 en 1367,2 Hz duidelijk zichtbaar. Ook nog t e herkennen, maar niet

zo duidelijk zijn de eigenfrequenties 1757,8 en 2343,8. Dit geeft de volgende

verhoudingen reeks: i

n

L I

1 - 1- - 2- - 3 - 4

3 3 (53)

Deze reeks komt niet overeen met de voor mij bekende frequentieverhoudin- gen van klokken. De mij bekende klokken die er het meest op lijken zijn de

gme klok met verhoudingen: 1 2

1 - 1- - 2 - 3 - 4 En de grote terts klok met verhoudingen:

1

1 - 2 - 2- - 3 - 4

2

(54)

( 5 5 ) Dat de gevonden freqenties niet exact overeen komen met de theoretische frequenties komt omdat bij de Gabor transformatie slechts naar discrete frequentisch gekeken kan worden.

Om wat nauwkeuriger naar de frequenties van deze klok te kunnen kijken wordt de Fourier transformatie van dit signaal bepaald. Deze Fourier coëffi- ciienten staan afgebeeld in figuur 10. In deze figuur staan ook de gemiddelde

Figuur 10: De Fourier coefficienten van het signaal “Ge115”, ook is de ge- middelde waarde van de Gabor coëfficienten tegen de frequentie uitgezet waarden van de Gabor coëfficienten afgebeeld. Deze figuur laat de volgende eigenfrequenties zien:

o 231,9 Hz o 463,9 Hz o 585,9 Hz

o 927,7 Hz

o 1367,2 Hz

De verhouding tussen deze frequenties is:

1 - 2 - 2,5 - 4 - 5 , 9 ₍₅₆₎ Ge115 komt dus van een grote terts klok, maar met een lagere grondtoon dan uit de Gabor transformatie gehaald werd. Hieruit blijkt dat de resolutie van de Gabor transformatie t e grof is om de eigenfrequenties goed t e kunnen bepalen. De grondtoon wordt ondergesneeuwd door de gecombineerde punt van de priem en de terts.

Omdat Ge115 slechts de enige opname is die teruggevonden is van het werk van Van Klinken is het niet mogelijk om toch nog te kijken of er verschil t e zien is tussen klokken die goed en slecht inslingeren. Ik verwacht d a t ook hier de resolutie te kort schiet om het inslingergedrag t e kunnen bekijken.

5

Conclusies

o De Gabor transformatie kost iets meer rekenwerk dan STFS.

o De Gabor transformatie geeft

,

indien een goede basisfunctie gekozen

worden, een iets betere resolutie dan STFT.

o Goede basisfuncties zijn functies die bi-orthogonaal zijn ten opzichte

van de functies van Gauss, Hamming en Hanning.

o De effectieve bandbreedte van de gebruikte functie heeft invloed op de resolutie.

o De beste resolutie wordt verkregen als de effectieve bandbreedte gelijk aan het totaal aantal punten gedeeld door het aantal tijdstappen

(5

= f l l ) wordt genomen.

o De Gabor transformatie geeft, bij de door mij gekozen basisfunctie, onvoldoende resolutie om de afzonderlijke eigenfrequenties t e bekijken. o De Gabor transformatie geeft zeer waarschijnlijk onvoldoende resolutie

6

Aanbevelingen

e De door mij geïmplementeerde MATLAB progra.mma’s zijn rechtoe rechttaan geprogrammeerd. De snelheid van de programma’s is waar- schijnlijk op t e voeren door de programma‘s zo om t e schrijven dat de speciseke mogelijkheden van MATLAB optimaal gebruikt worden.

e Ik heb voor de basisfuncties gebruik gemaakt van een aantal functies die al standaard in MATLAB geïmplementeerd waren. Andere func- ties, zoals bijvoorbeeld de Pseudo Wigner-Ville, of de optimale functies zoals die door Qian e.a. [7] zijn afgeleid, zouden een beter resultaat kunnen geven.

o Bij de transformatie van het opgenomen carillongeluid is t e zien dat

de energie in het lage frequentiegebied zit. De coëfficienten in het hoge gebied worden eigenlijk voor niets bepaald. Uit het signaal zou meer informatie gehaald kunnen worden als een soort ‘zoom-Gabor’, vergelijkbaar met zoom FFT, wordt gebruikt.

A Uitwerking H-matrix

In deze bijlage wordt de matrix uit vergelijking 28 gegeven voor het geval

M = 4 en N = 4. Hieruit volgt dat:

P = M N = 2 * 4 = 8 ( 5 7 )

hi =

h*[i]

( 5 8 )

Voor het uitschrijven van vergelijking 27 wordt een nieuwe schrijfwijze in- gevoerd :

De H-matrix wordt dan: mn O0 o 1 02 03 10 11 12 13 hoWO ho W o ho

W o

ho W o h4 W o h4 W o h4 W o h4 W o hl W o hl W 1 hl W 2 h5 W o h5 W 1 h5 W 2 h5 W 3 hl W 3 h3 W o h3W3 h3W6 h3 W 9 h7 W o h7 W 3 h7 W 6 h7 W 9 h4 W o h4 W 4 h4 W s h4W12 ho W o ho W s howl2 how4 h5 W o h5 W 5 h5Wl0 h5W15 hl W o hl WIO hl W 5 hlW15 h7W0 h7W7 h7 WI4 h7W21 h3 W o h3 W 7 h3 W14 h3 W21 (59

Voor de constante W geldt:

(60)

WN+n = W N

IE

dit g e ~ d is N geEjk aan

4:

Bierlit volgt voor de voorgaande H-matrix:

mn

O0

o1 02 o3 10 11 12 13 We zien ho W o ho W o ho W o ho W o h* W o h4 W o h4 W o h4

W o

hl W o hl W 1 hl W 2 hi W 3 h5 W o h5 W 1 h5 W 2 h5

W 3

h2 W o h2 W 2 h2 W o hz W 2 h6W0 h6W2 hpj W o h6W2 h3 W o h3 W 3 h3 W 2 h3 W 1 h7 W o h7 W 3 h7 W 2 h7W1 h4 W o h4 W o h4 W o h4 W o ho W o ho W o ho W o ho W o h5 W o h5 W 1 h5 W 2 hl W o hl W 1 hl W 2 h5

w3

hl W 3 h6

W o

h6W2 h6W0 h6

w 2

h2 W o h2 W 2 h2 W o h2

W 2

h7 W o h7

W 3

h7W2 hTW1 h3W0 h3 W 3 h3W2 h3W1

iier nu inderdaad 2 x 2 submatrices met de afmeting 4 x 4.

I3

Gabos programma’s

B . l MaakGamC

function gam=maakgamc(h,M,N) ;

x

X

gam=maakgam(h,M,N) ;

x

x

1

1 N

frequentiestappen.

x

X

X

Ì! h : Basis functie

x

M

: Aantal tijdstappen

x

N

: Aantal frequentiestappen

x

x

x

Deze functie bepaalt de bi-orthogonale reeks voor de functie h voor de onderverdeling in

M

tijdstappen en Deze functie is slechts geschikt voor kritisch sampelen Omdat dit programma slechts geschreven is voor reele functies wordt de conjungatie

van

h niet uitgevoerd.

x

gam : Bi-orthogonale functie van h

P=max ( s i z e (h) ;

i f P == M*N

% Bepaal de l e n g t e van de f u n c t i e .

X

Controle of de l e n g t e w e l goed i s . U=exp(i*Bpi/N); Maak de constante W.

f o r n=l:N-1

end

f o r n=l:N-1

WW(n)=W^(-n); % Maak de e e r s t e r i j van de 'WY-matrix

WW(n,:)=WW(l,:).^n;

X

De andere r i j e n van de

X

'

W ' -matrix maken. end

WW= [ones (1, N)

X

De O-de r i j en kolom v u l l e n ones (N- 1

,

1 )

WWl

;

X

m e t W ^ O = l .

H=[1;

f o r l = O : M - 1

Hl=WW*diag(h(l:W));

X

De Hl-matrix berekenen. H=[H

H 1 1

;

X

Deze H1 aan de e e r s t e r i j

%

h=[h(N+l:P) h(l:N)];

X

De h-vector N p o s i t i e s

X

verschuiven zodat de

X

volgende h ' s goed s t a a n . van de H matrix toevoegen.

end

f o r 1=0:M-2 % De volgende r i j e n van de H m a t r i x d i e H=CH;

X

ten opzichte van e l k a a r verschoven z i j n

H ( l * N + l : (l+l)*N

,

N + l : P

...

X

toevoegen. H ( l * N + l : ( l + l ) * N

,

1 :

N I ] ;

end

Xmesh(abs (Hl ;

%pause

%contour (abc (HI ;

%pause

gam=(H\nu)’; % De gamma vector uitrekenen

MaxIm=max(imag(gam)) % Geef de maximale waarde van

X

net imaginaire deei. gam=real(gam); % in gevai van reele functies is de

% bi-orthgonale functie ook reeel.

1

Het imaginaire gedeelte is het gevolg

%

van rekenonnauwkeurigheden en heeft

% vertragingen in het rekenproces tot

% gevolg daarom wordt dat weggehaald end

E.2 GahssC

%

x

% % x : Te analyseren signaal % % M : Aantal tijdstappen % N : Aantal frequentiestappen

% dt : Tijd tussen 2 sampels

1

% a : De coefficienten matrix % t : Tijdvector

X

f

: Frequentievector

function [a, t ,fl =gaborc(x ,gamma,M ,N, dt) ;

%

[a, t ,f] =gaborc (x ,gamma,

M,

N,

dt) ;

Bepaalt de Gabor transformatie bij kritisch sampelen

gamma : Bi-ortogonale functie van de basisfunctie

%

% Definitie van enkele constanten P=max(size(x));

if M*N==P

% Bepaal de lengte van het signaal

% Controle of de lengte wel goed is W=exp(i*2*pi/N);

%

Maak de constante

W

a=zeros(M,N); % Alle coefficienten op O zetten for m=O:M-1 % Voor iedere m

for n=O:N-1 % Voor iedere

n

for k=O:P-1 3i, Som over k

c=k-m*N;

1

Index voor de functie 'gamma' cc=c;

while cc<O % Gamma is een periodieke functie cc=cc+P ; % met periode P. Zet de index zo

% dat deze binnen CO. .P-1 valt] end;

gammamn=gamma (cc+ 1) *W^ (-n*k) ;

%

Bepaal de zogenaamde gammamn voor dit punt

%

Tel het product van deze gammamn met xCkl a(m+l,n+l)=a(m+l,n+l)+x(k+l)*gammamn;

%

% g a m m Ckl *xCkl

.

op bij de hiervoor uitgerekende end

disp( Cm n a(m+l ,n+l>l) ;

% Laat zien hoever het programma is plus het

%

%

resultaat zodat de gebruiker niet naar een zwart scherm hoeft te kijken.

end end

t=O:dt*N:dt*P-dt; % Bepaal de tijdvector

f=O:l/(dt*N):(N-l)/(dt*N); % Bepaal de frequentievector end

B.3 GabInvC

function x=gabinvc(a,h) ;

%

x=gabinvc (a

,

h) %

X

%

%

Reconstrueert het signaal weer uit de Gabor coefficienten. Deze functie is slechts geschikt voor kritisch gesampelde

% Gabor transformatie

%

% a : De Gabor coeficienten

% h : De basisfunctie

h x : Het gereconstrueerde signaal

i/i

%

O'

[M

,

NI

=size (a) ;

% M

: Het aantal tijdstappen

% N : Het aantal frequentiestappen P=max(size(h));

if M*N==P

x

Bepaal de lengte van het signaal

% Controle of de lengte wel goed is W=exp(i*2*pi/N); % Maak de constante

W

x=zeros(l,P);

%

Alle waarden op O zetten

for k=O:P-1 % Voor iedere k

for m=O:M-1 % Som over m for n=O:N-1 % Som over n

c=k-m*N ; % Index voor de functie 'h' cc=c;

while cc<O

end ;

X

binnen

[O.

.p-11

% h is een periodieke functie cc=cc+P;

%

met periode P. Zet de index

hmn=h(cc+l) *WA (n*k) ;

%

Bepaal de zogenaamde

hmn

voor

1

deze m en n

X

Tel het product van deze hmn met xckl

%

% hmnCkl*x[kl

x(k+l)=x(k+l)+a(m+l,n+l)*hmn;

op bij de hiervoor uitgerekende end

disp( Ck x(k+l)l)

% laat zien hoever het programma is plus het

% resultaat zodat e gebruiker niet naar een

% zwart scherm hoeft te kijken end

#axImag=max(imag(x)) x=real(x);

end

B.4

MaakGamO

function gam=maakgamo (h,

M

,

N

,gamin) ;

% % gam=maakgamo(h,M,N,gamin)

x

%

1

%

N

frequentiestappen.

%

%

% % h : Basisfunctie % M : Aantal tijdstappen % N : Aantal frequentiestappen

% gamin : De opgelegde waarde van gamma (deze parameter

% hoeft niet meegegeven te worden, als een waarde

% niet opgelegd wordt, moet hier de constante NaN

x

opgegeven worden.

%

X

gam : Di-orthogonale functie

van

h

%

Deze functie bepaalt de bi-ortogonale reeks voor de functie h voor de onderverdeling in M tijdstappen en Deze functie is ook geschikt voor oversampelen. Omdat dit programma slechts geschreven is voor reele functies wordt de conjungatie

van

h niet uitgevoerd.

%

P=max(size(h)); % Bepaal de lengte van de functie.

Ml=P/N;

Nl=P/M;

X

Definieer de nieuwe parameters

f o r n=l:Nl-1

WWí(n)=WiA(-n) ;

1

Maak de e e r s t e r i j van de

i/,

’

W’ -submatrix end f o r n=2:Nl-1 WWl(n,:)=WWl(l,:).^n;

X

De andere r i j e n van de

X

’

W ’ -submatrix maken. end

W W l = Cones ( i ,N i )

X

De O-de r i j en kolom v u l l e n ones(Mi-i,l)

W W i 1 ;

Y, m e t W A O = í .

ww=

c1;

f o r m=l:M end W W l = W W ;

ww=

cww

W W i l

;

X

De e e r s t e r i j van de W-matrix maken

f o r m=2:M1 end

ww=cww;

W W i l

;

X

De andere r i j e n van de W-matrix

HH=

Cl

; f o r n=l:Nl end

HH=

CHH; hl ;

X

De e e r s t e r i j H-submatrices f o r 1=0:M1-2

X

De volgende r i j e n van de H m a t r i x d i e t e n H H ( l * N l + l : ( l + l ) * N l

,

N + l : P )

...

1

toevoegen H H ( l * N l + l : ( l + l ) * N l

,

1 :

N I ] ;

HH=[HH; opzichte van e l k a a r verschoven z i j n

end

H = HH

.*

WW;

X

S t e l de H-matrix samen u i t

nu=[P/(M*N); zeros((Ml*Nl)-1,l)l ;

%

De nu vector

if nargin == 4 % Als een gamma-in vector wordt opgegeven

% Terug tellen omdat data gewist wordt

X

en bij vooruittellen de nummering

X

niet meer klopt for n=P:-l:l if gamin(n) *= NaN nu=nu-gamin(n)*H(:,n); H=[H(: ,l:n-1) H(: ,n+l:max(size(H~~)I; end end end

%

De gamma vector uitrekenen

%

De gamma vector uitrekenen

X

Hier wordt de gegeneraliseerde inverse oplossing genomen gam=(H’*inv(H*H’)*nu)’;

%

De gamma’s van gamma-in tussenvoegen if nargin == 4 for n=l:P if gamin(n) --= NaN gam=Cgam( :

,

1 :n-l) gaminhl

. .

gam(: ,n:max(size(gam>>>l ; end end end MaxIm=max(imag(gam)) gam=real(gam);

1

Geef de maximale waarde van

1

het imaginaire deel. In geval van reele functies is de

1

bi-orthgonale functie ook reeel.

% Het imaginaire gedeelte is het gevolg van

X

rekenonnauwkeurigheden en heeft

X

X daarom wordt dat weggehaald

B.5 GaborO

function [a,t ,f]=gaboro(x,gamma,M,N,dt);

%

11.

%

X

Bepaalt de Gabor transformatie

%

% x : Het te analyseren signaal

%

gamma : Het window

%

M

: Het aantal tijdstappen

% N : Het aantal frequentiestappen

% dt : Tijd tussen 2 sampels

% % a : De coefficienten matrix % t : Tijdvector % f : Frequentievector

x

Ta $ t f

I

=gaboro

P=max(size(x)); % Bepaal de lengte van het signaal Ml=P/N;

N

l=P/M

;

1

Definieer de nieuwe parameters

W=exp(i*2*pi/N);

X

Maak de constante

W

a=zeros(M,N); % Alle coefficienten op O zetten for m=O:M-1

1

Voor iedere m

%

Voor iedere n

%

Som over k for n=O:N-1

for k=O:P-1

c=k-m*N1; % Index voor de functie 'gamma' cc=c;

while ccC0 % Gamma is een periodieke functie cc=cc+P;

1

met periode P. Zet de index zo

dat deze binnen [O. .P-1 valt]

%

end end t=O: dt*N1 :i

X

X

t*P-I gammamn=gamma(cc+l) *W^(-n*k) ;

a(m+l ,n+l) =a(m+l ,n+l)

+x

(k+l) *gammamn;

1

Bepaal de zogenaamde gammamn voor dit punt

X

Tel het product van deze gammarm met xCkl

X

op bij de hiervoor uitgerekende

1

gammamn Ckl *x Ckl

.

end

disp( [m n a(m+l ,n+l)l) ;

X

Laat zien hoever het programma is plus het resultaat zodat de gebruiker niet naar een zwart scherm hoeft te kijken.

t;

X

Bepaal de tijdvector df=l/(P*dt) ;

f=O:df*Ml:df*P-df;

1

Bepaal de frequentievector

B.6 GabInvO

function x=gabinvo(a,h) ;

X

X

x=gabinvo(a,h)

x

X

X

X

a : De Gabor coefficienten

X

h : De basisfunctie

'/o

Reconstrueert het signaal weer uit de Gabor coefficienten.

X

X

x

: Het gereconstrueerde signaal

P=max(size(h));

X

Bepaal de lengte van het signaal [M,N]=size(a) ;

X M : Het aantal tijdstappen

X N

: Het aantal frequentiestappen Ml=P/N;

N 1 =P/M ; W=exp(i*2*pi/N);

X

Maak de c o n s t a n t e W

..

x = z e r o s ( i , P j ; A Â i i e waarden op O z e t t e n f o r k=O:P-1

X

Voor i e d e r e k

X

Som over m f o r n=O:N-1

X

Som over n f o r m=O:M-1

c=k-m*Nl ;

X

Index voor de f u n c t i e 'h'

cc=c ;

while cc<O

end ;

X

binnen [O. .p-11 hmn=h(cc+i)*W^(n*k) ;

X

Bepaal de zogenaamde hmn voor

X

deze m en n

X

Tel h e t product van deze hmn met xCkl

X

X hmnCkl*xCkl

h i s een p e r i o d i e k e f u n c t i e cc=cc+P;

X

met p e r i o d e P. Zet de index

x (k+ 1) =x (ka 1) +a (ma 1

,

n+ 1 ) *hmn ; op b i j de h i e r v o o r u i t g e r e k e n d e end end d i s p ( Ck x ( k + l ) l )

X

l a a t z i e n hoever h e t programma i s p l u s h e t

X

r e s u l t a a t zodat e g e b r u i k e r n i e t n a a r een

X

zwart scherm h o e f t t e k i j k e n end

MaxImag=max (imag (XI x=real (x) ; end

B.7 Gabor

f u n c t i o n Ca

,

t

,

f

1

=gabor (x

,

gamma

,

M

,N , d t ) ;

X

X

X

Ì! Ì! Ì! Ì! ti,

x

Ì! Ì! Ì! Ì! Ì! Ì! Ì!

x

Ì!

Bepaalt de Gabor transformatie.

Hierbij wordt, indien mogelijk, gebruik gemaakt van het FFT algoritme.

X : Te analyseren signaal

gamma : Bi-orthogonale functie van de basisfunctie

M

: Aantal tijdstappen

N

: Aantal frequentiestappen dt : Tijd tussen 2 sampels a : De coefficienten matrix t : Tijdvector

f : Frequentievector

P=max(size(x)); Ì! Bepaal de lengte van het signaal Ml=P/N;

Nl=P/M;

1

Definieer de nieuwe parameters

W=exp(i*a*pi/N); Ì! Maak de constante

w

a=

[I

; Ì! De coeff icienten leeg maken

for m=O:M-1

X

Voor iedere m y=x

.*

gamma;

yl=zeros (1 ,N) ;

for l=O:Ml-I Ì! Sommeren over 1 end

a=Ca;

am1

;

yl=yl+y(l*N+l: (l+I)*N) ;

am=f f t

(y

i) ;

gamma= [gamma(

(M-

i) *Nl+l : P) gamma( 1 : (M-I) *Ni)] ;

disp (m) ;

Ì! De gamma voor de volgende cyclus goed zetten end

t=O:dt*Nl:dt*P-dt;

X

Bepaal de tijdvector df = 1

/

(P*dt ) ;

f=O:df*Ml:df*P-df;

X

Bepaal de frequentievector

C

Test programma’s

(2.1

BiOrthog

disp(’Bi0rthog’);

disp(’Dit programma test of de bi-orthogonale functie’); disp(’van de bi-orthogonale functie

van

een functie’); disp(’weer dezelfde functie opleverd’);

disp(’

’

;

disg(’Druk op een toets’); pause

h=(rand(1,64))-0.5; h=normaal (h) ;

bar (h) ;

title(’Een willekeurige functie’); pause

disp(’Berekening van de bi-orthogonale functie’); gam=maakgamc(h,8,8);

bar

(gam)

;

title(’De bi-orthogonale functie hiervan’); pause

disp(’Berekening van de bi-orthogonale functie’); hl=maakgamc (gam, 8,8) ;

bar (hi) ;

title(’De bi-orthogonale functie de vorige functie’); pause

bar(h-hl) ;

title(’Het verschil van de eerste en de laatste functie’); dif =abs (h-hi) ;

disp(’Dit wordt voor N willekeurige functie bepaald’); N=50

dispi’ï~ruk op een toets’;;

pause for n=2:N disp(n) ; h=(rand(l ,641 -0.5; h=normaal(h); gam=rnaakgamc(h,8,8); hl=maakgamc(gamY8,8) ;

dif =dif +abc (h-hi) ;

end

bar(dif

/NI

;

title(’Gemidde1de verschil’); C.2

RechtHk

disp (

’

RechtHk

’

;

disp(’Dit programma test een aantal eigenschappen van’); disp(’de door MaakGamC bepaalde bi-orthogonale functie’); disp(”);

disp(’Druk op een toets’); pause

h= [zeros (1,28) ones (i, 8) zeros

(I,

2811 ;

h=nomaal (h) ;

bar (h) ;

title( ’Een rechthoekig venster’) ;

pause

disp( ’Berekening van de bi-orthogonale functie’) ;

gam=maakgamc (h, 8,8) ;

bar (gam) ;

pause b a r (h-gam) ; t i t l e 0 H e t bewijs: h

-

gam’)

C . 3

Gaeissven

d i s p (

’

GaussVen’ d i s p ( > L a a t de bi-orthogonale f u n c t i e van h e t ’ ) d i s p ( ’Gausse v e n s t e r z i e n b i j v e r s c h i l l e n d e ’ ) d i s p ( waarden van s

’

; d i s p (

’ ’

; disp(’Druk op een t o e t s ’ ) ; s=[.7 .6 .5 .4 .3

.11;

n = l e n g t h ( s ) ; N=64 ; x3= [i NI ; y3=[0 01 ; f o r i = l : n pause x = l i n s p a c e ( - l , l , N ) ; x = g a u s s ( x , s ( i ) ) ; x=normaal

(XI

; b a r ( x ) ; t i t i e O V e n s t e r ’ ) ; pause d i s p ( ’Bereken bi-orthogonale f u n c t i e ’ ) ; d i s p ( s ( i ) ) gam=maakgamc (x ,8,8) ; [xl,yll=bar(gam) ; [x2, y21 =bar (x) ; t i t l e o v e n s t e r met bi-orthogonale f u n c t i e ’ ) ; p l o t ( x i , y 1

,

x2 ,y2, x3 ,y31 ;

end

C.4

TestAOO

disp(?TestAûû’)

disp( it prograiia test of de Gabor trzmsfcrmatie win’) ;

disp( ’een willekeurig signaal met als basisfunctie dat

’1

;

disp(’ze1fde signaal de matrix:’); disp(’ 1 O O

..

O’); disp(’ O O O’); disp(’

.

_{. ’ I ;}

disp(’

.

. ’ I ;

disp(> O

..

O’); disp (

’

opleverd.

’

; disp(’

’1;

disp(’Druk op een toets’); pause

x=rand(l,64) -0.5; x=normaal (x) ;

title(>Een willekeurig testsignaal (x)

’1

;

pause; plot (XI ;

disp(’Bereken bi-orthgonale functie van x’); gam-aakgamc (x ,8,8) ;

disp( ’Bekeken Gabor coeff icienten’) ;

[a,t ,f

1

=gaborc (x ,gam, 8,8,1/50) ;

plotgab(abs(a),t,f,’punt m=O,n=O is 1, de rest is O’) pause;

a(1, i)=a(i, 11-1;

disp(’De maximale absolute waarde van de coefficienten’); max(max(abs(a)))

C.5

GabRecht

disp (

’

GabRecht

’

;

disp ( ’Dit programma voert een Gabor transformatie van een’ ) ;

disp

<

sirdsvomlg s l g z â â l U i t m e t eP,n rechth^ekig venster

’

1

;

diSF(”);

disp(’Druk op een toetsy);

% Rechthoekig venster

h=[zeros(l,28) ones(l,8) zeros(l,28)1; h=nomaaï (h) ;

% Maak hulp tijd voor een mooie sinus th=linspace (O, 2*pi, 65) ;

th=th( 1 : 64) ;

% Maak de echte tijd dt=1/50 ;

t=O:dt:63*dt; for c=1:5

pause

X

Maak sinusvormig testsignaal x=sin(th*c*8) ;

title(’Het ingangssignaal’); pause

plot (t 9x1

disp(>Het berekenen

van

de Gabor coefficienten’);

[amn, time ,f reql =gaborc(x ,h,8,8 ,dt) ;

X

h=gamma plotgab (abs (amn) ,time

,

f req,

. .

>Gabor transformatie van een sinus signaal’); end

C.6 GabGauss

disp (

’

GabGauss

’

;

disp( ’Dit programma voert een Gabor transformatie van een’) ;

disp(>sin~sï~ri.,ig signaal dat n i e t precies in het’) ;

disp(>(rechthoekig) venster past.’);

disp( ’Hierna wordt datzelf de signaal getransformeerd met

’

) ;

disp(’een Gausse basisfunctie’); disp(’

’1;

disp(’Druk op een toets’); pause

X

Rechthoekig venster

h=[zeros(l,28) ones(l,8) zeros(l,28)1; h=nomaal (h) ;

h=f lipit (h) ;

X

Maak de tijd dt = 1/50 ;

t=O:dt:63*dt;

X

Maak sinusvormig testsignaal x=sin(t*50) ;

title( ’Het ingangssignaal’) ;

pause plot (t ,XI

disp(’Het berekenen van de Gabor coefficienten’); [amnr ,time ,freq] =gaborc(x ,h, 8,8, dt) ; % h=gamma

plotgab(abs(amnr) ,time,freq,.

.

pause h=linspace(-1,1,64); h=gauss (h, .35) ; h=normaal (h) ; plot (h) title(’De basisfunctie’) ; pause

disp( ’Bereken bi-orthogonale functie’) ;

gam=maakgamc (h

,

8,8) ;

pause plot (gam)

gam=f

lipit (gam) ;

disp(’Het berekenen van de Gabor coefficienten’); [amng

,

time

,

f req] =gaborc (x ,gam, 8,8

,

dt) ;

plotgab(abs(amng),time,freq,.. pause

’Gabor transformatie met Gauss’);

plotgab(abs(amng)-abs(amnr),time,freq,.. ’Coefficienten Gauss

-

Rechthoek’);

(3.7

Vensters

disp(’Vensters’);

disp(’Met dit programma wordt een sine sweep met’); disp(’verschil1ende basisfuncties getransformeerd’); disp(’

’1;

disp( ’Druk op een toets’ ;

pause M=4 ;

N=M*S;

P=M*N; dt=1/50 ; x=sinsweep(P) ; plot(O:dt:(P-l)*dt,x); title(’De sine sweep’); pause

hrech=normaal([zeros(l, (P-N)/2) ones(1,N) zeros(1, (P-N)/2)]); plot (hrech) ;