Dictaat introductie MatLab

Document Sample
Dictaat introductie MatLab Powered By Docstoc
					              MatLab
            Een korte inleiding
            Nol Bendermacher
         Rinske de Graaff Stoffers
                April 2005




MatLab                 1             2-3-12
1        Oriëntatie op het cursusonderdeel
Voor veel richtingen in de psychologie is signaalverwerking een noodzakelijk en
belangrijk onderwerp. Het gebruik van software om gemakkelijk geavanceerde
bewerkingen op signalen uit te voeren is inmiddels gemeengoed geworden. Daarom is
deze cursus (2 studiepunten) opgezet. Na 80 uur, die voornamelijk uit zelfstudie zullen
bestaan, heb je een goed idee van wat er met MatLab allemaal mogelijk is. Om MatLab
daarna zelfstandig te kunnen gebruiken zal echter veel meer tijd en inspanning nodig zijn.
Deze cursus is dus niet meer dan een inleiding. Na voltooiing van dit onderdeel moet je
in staat zijn in MatLab:

    1.   de basisfuncties te gebruiken
    2.   de juiste taal te gebruiken om MatLab aan te sturen
    3.   signaaldata op te vragen uit een bestand
    4.   een Fourier-analyse uit te voeren
    5.   signalen te filteren
    6.   je resultaten weer te geven in een overzichtelijke plot

Het rooster van de cursus is te vinden op blackboard.ru.nl (onder Course Documents).
Daarop zie je ook welke stof in welke week besproken wordt en waar je terecht kunt voor
het maken van de computeropdrachten.




MatLab                                          2                                       2-3-12
2        Introductie MatLab
MatLab is een afkorting van matrix laboratory. Het is een programmeertaal voor
technische berekeningen. Het integreert rekenen, visualiseren en programmeren in een
gemakkelijke omgeving waarbij problemen en oplossingen uitgedrukt worden in
gebruikelijke mathematische notatie. Door de feedback van veel gebruikers heeft MatLab
zich ontwikkeld tot een programma, dat op veel universiteiten wordt gebruikt, zowel in
de opleiding van studenten als in onderzoek. MatLab wordt ingezet voor data-acquisitie,
data-exploratie, data-analyse, visualisatie, ontwikkeling van algoritmen, modelvorming,
simulatie en applicatie-ontwikkeling.

Het MatLab-systeem bestaat uit vijf hoofdonderdelen:

   De taal: de programmeertaal waarmee je MatLab bestuurt.

   De werkomgeving: Dit is de omgeving waarin je als MatLab-gebruiker werkt om je
    gegevens binnen te halen, te bewerken en op te slaan.

   Het grafische systeem: Dit geeft je de mogelijkheid twee- en driedimensionale
    afbeeldingen te maken en deze tot in detail aan te passen.

   De mathematische functies: Een enorme collectie rekenalgoritmes, van zeer simpele
    functies, zoals optellen en aftrekken, tot zeer complexe, zoals fast Fourier transforms.

   Het Application Program Interface (API): Dit geeft je de mogelijkheid om MatLab te
    koppelen aan programma's in C en in Fortran.

Naast de basisomgeving van MatLab bestaan er een aantal gereedschapskisten
('toolboxes') die ieder een bepaald deelterrein bestrijken. Ze bestaan uit een collectie
MatLab-functies (M-files). Wanneer je in de basisomgeving een commando uit een
bepaald deelterrein geeft, wordt dit commando opgeroepen in de betreffende toolbox en
vervolgens uitgevoerd. Toolboxes die voor psychologen relevant zijn, zijn bijvoorbeeld
de optimization toolbox (non-lineaire curve fitting en regressie), de signal processing
toolbox (signaalverwerking, filtering en frequentie-analyse), de neural network toolbox
(ontwerp, implementatie, visualisatie en simulatie van neurale netwerken) en de statistics
toolbox (geavanceerde statistische analyses, non-lineaire curve fitting en regressie). In dit
cursusonderdeel zullen we vooral gebruik maken van de signal processing toolbox. Een
onderdeel waar we geen gebruik van zullen maken, maar dat wel het noemen waard is, is
“Simulink”. Hiermee kun je blokdiagrammen ontwerpen om systemen weer te geven,
vervolgens hun gedrag te simuleren en hun prestatie te evalueren, om zodoende tot een
blokdiagram te komen dat het systeem zo goed mogelijk benadert. Ten slotte biedt
MatLab ook een Compiler, die koppelingen met programma’s in andere programmeer-
talen (zoals C en C++) mogelijk maakt.




MatLab                                          3                                          2-3-12
Mocht je geïnteresseerd zijn in meer informatie over de verschillende onderdelen van
MatLab, dan kun je terecht op de homepage van Mathworks (de fabrikant):
www.mathworks.com

3        Kennismaking met MATLAB
Alvorens we op een meer inhoudelijk niveau aan de slag kunnen, moeten we eerst
kennismaken met de omgeving waarin we zullen gaan werken. In dit hoofdstuk zitten een
aantal opdrachten die je daarbij moeten helpen. Het is belangrijk, dat je deze opdrachten
daadwerkelijk uitvoert. Je zult zien dat je op deze manier al heel snel vertrouwd raakt met
de algemene structuur en de commando’s die binnen MatLab gangbaar zijn. Hierdoor
kun je je in latere hoofdstukken meer en meer aan de inhoudelijke uitdaging overgeven,
zonder nog met elementaire problemen te stoeien.

3.1      Opstarten
Wanneer je MatLab opstart, verschijnt er een sober en eenvoudig venster voor je (zie
figuur 1). Je zult al meteen zien dat het dezelfde opbouw heeft als bij ieder ander
programma onder Windows.




Figuur 1: Het beginscherm




MatLab                                         4                                         2-3-12
Bovenaan in je scherm zie je een balk met commando’s als file, edit, view, window, help.
Tot dusver is er dus niets nieuws onder de zon. Onder deze balk zit nog een balk (de
'toolbar') met daarop een aantal knoppen waarmee je snel een actie kunt uitvoeren. In het
vervolg van dit hoofdstuk zullen de belangrijkste onderdelen van deze balken
geïntroduceerd worden. Voordat we dat gaan doen, zullen we eerst eens kijken wat voor
soort commando’s je kunt intypen.

3.2       Commando’s en interpunctie
Zoals je ondertussen wel weet, is MatLab een wiskundig programma. Het is dus logisch
dat een groot deel van de commando’s die MatLab kent, wiskundig van aard zijn. Sterker
nog, de meeste commando’s komen overeen met de gebruikelijke wiskundige notatie
(een overzicht van deze commando’s krijg je door help gevolgd door een spatie en een
slash (help /) in te typen).

Opdracht 3a:

Typ in het Command Window van MatLab het volgende in:
1+2 <enter>


Typ nu in:
1          +      2   <enter>


MatLab leest je opdracht van links naar rechts, zonder hierbij te letten op de afstand
tussen de verschillende commando's. Wanneer je op enter drukt, begint MatLab met het
verwerken van het commando.

Typ in:
1 + 2; <enter>


Je zult tot de conclusie komen dat MatLab je het antwoord nu schuldig blijft, althans, het
laat het antwoord nu niet zien. Om toch het antwoord in beeld te krijgen typ je:
ans <enter>


Het lijkt op dit moment onzinnig dat MatLab je na een puntkomma geen uitkomst laat
zien, maar wanneer er later grote tussenberekeningen nodig zijn, zul je nog dankbaar
gebruik maken van deze mogelijkheid.
Je kunt ook meerdere commando's op één regel zetten; achter elk commando behalve het
laatste moet dan een komma staan of een puntkomma: een komma als je het resultaat wilt
zien en een puntkomma als je dat niet wilt.

Probeer nu zelf een aantal eenvoudige berekeningen uit te voeren (denk aan: optellen,
aftrekken, delen, vermenigvuldigen, machtsverheffen, gebruik van haakjes etc.).

Typ nu in:



MatLab                                        5                                         2-3-12
 1 a 2 <enter>


Je krijgt dan een foutmelding van MatLab. In dit geval is de foutmelding standaard
(missing operator, comma, or semi-colon). We zullen nu eens proberen een wat
duidelijkere foutmelding uit te lokken.

Typ in:
 1 + (2+3     <enter>


Je ziet dat de foutmelding je nu zeer concrete informatie geeft, die je helpt de fout te
verbeteren. In de praktijk zul je zien dat er meer fout dan goed gaat (althans bij eerste
pogingen) en de foutmelding van MatLab is altijd een belangrijk aanknopingspunt om je
fout te vinden. Dit geldt zowel voor eenvoudige interpunctiefouten zoals hierboven, als
voor fouten ergens in een grotere en complexere reeks commando's (oftewel: een
programma).

Tot slot van deze opdracht nog een aantal tips.

Typ eens een lang commando in met daarin een eenvoudige interpunctiefout,
bijvoorbeeld
 1+5/3*5+(3*4)+34/3*(3*2 <enter>


Als je nu op grond van de foutmelding het commando wilt verbeteren, kun je dat doen
door de gehele formule opnieuw in te tikken, maar er is ook een tweede, snellere
methode. Druk op de pijltjestoets 'omhoog'. Je zult zien dat bij iedere druk op deze toets
een voorgaand commando verschijnt. Blijkbaar heeft MatLab alles wat je hebt ingetypt
automatisch bewaard. Je kunt dus heen en weer scrollen naar vorige en volgende
commando's. Zo'n commando kun je dan aanpassen en opnieuw uitvoeren (via <enter>).
Ook handig om te weten, is dat namen van variabelen altijd met een letter dienen te
beginnen, gevolgd door een willekeurig aantal letters, cijfers of underscores (let op: voor
MatLab zijn hoofdletters en kleine letters niet gelijk!). Tenslotte dien je in getallen een
decimale punt in plaats van een decimale komma te gebruiken.

3.3       Command Window, Workspace en M-Files
Het venster waar we tot nu toe in gewerkt hebben, wordt het Command Window
genoemd. De commando’s die je hierin geeft, worden direct uitgevoerd in de workspace.
In het Command Window heeft een <enter> dus een handeling tot gevolg. Het schrijven
van programma's doe je daarom niet hier, maar daarover straks meer.

Opdracht 3b:

We gaan nu een klein programmaatje maken.
Typ in:




MatLab                                            6                                      2-3-12
 a=2; <enter>
 b=3; <enter>
 c=a+b; <enter>


Wanneer je nu c intypt (en een <enter>), krijg je de uitkomst van a+b. Maar als je op
deze manier de som van twee andere getallen zou willen berekenen, zou je weer eerst a
en b een waarde moeten geven en dan c uitrekenen. Hier komt het voordeel van de M-file
voor de dag.

We gaan nu naar een ander venster in MatLab, het zogenaamde M-file-venster, waar we
zullen zien dat het rijtje instructies uit opdracht 3b (licht aangepast) een echt programma
wordt.

Opdracht 3c:

Ga met de muis naar de bovenste balk op je scherm en klik File aan. Achter File zitten
veel nieuwe mogelijkheden, waarvan we er nu slechts één nodig hebben en wel de
mogelijkheid New. Kies M-file. Zodra je deze optie hebt aangeklikt, verschijnt er een
nieuw venster. Maak dit voor het gemak even paginagroot (zie rechtsboven in beeld).

Je ziet nu bovenaan het venster dat het automatisch de naam ‘untitled’ heeft gekregen.
We gaan deze M-file meteen opslaan onder een andere naam. Klik weer op File en
vervolgens op Save As. MatLab slaat M-files standaard op in de directory ‘work’.
Aangezien meerdere mensen dezelfde opdrachten maken op de onderwijs-PC’s, is het
verstandig eerst een eigen map aan te maken binnen deze directory. Geef deze map je
eigen naam. Bewaar de M-file vervolgens in je eigen map onder een logische naam (deze
moet beginnen met een letter), zodat je hem later nog terug kunt vinden (in dit geval
opdracht3c). Gebruik nooit punten of komma’s in namen van files: dat geeft problemen
in de workspace.

Typ nu in het venster van de M-file opdracht3c het volgende in:
a=input (‘a=’); <enter>
b=input (‘b=’); <enter>
c=a+b


Figuur 2 laat deze situatie zien.




MatLab                                         7                                         2-3-12
Figuur 2:      Een opdracht in de M-file editor

N.B.: Mocht je een commando in willen typen dat langer is dan één regel, geef dan aan
      het eind van het eerste deel drie puntjes (…) gevolgd door een enter. Als je nu je
      tekst op de volgende regel typt, begrijpt MatLab dat het om één commando gaat.

Save dit programmaatje en zet het M-file venster onderaan op je scherm neer (zie
rechtsboven). Teruggekomen in het Command Window kun je het programmaatje nu
testen. Om te zorgen dat MatLab niet automatisch in de standaarddirectory ‘work’ gaat
zoeken, maar in jouw eigen subdirectory, typ je eerst Addpath(‘c:\matlab6p1\work\…’).
Op de plaats van de … voer je uiteraard de naam in, zoals jij die gegeven hebt aan je
eigen subdirectory. Typ opdracht3c (en <enter>) in en het programma wordt uitgevoerd.

N.B. Elke keer als je MatLab opnieuw start en bestanden uit je eigen subdirectory wilt
     gebruiken, moet je eenmalig dit Addpath-commando geven.

 Maak nu de M-file weer groot en bekijk het programma. Allereerst valt je het commando
input op. Het zorgt ervoor, dat je programma om input vraagt. Ten tweede viel je op dat
de tekst 'a=' in beeld verscheen toen je het programma draaide. Dat is de tekst, die in het
programma achter het woord input is gegeven. Tot slot staan er achter de eerste twee
regels wel puntkomma's, maar achter de derde niet. De reden hiervoor is eenvoudig: je
wilt niet meer uitkomsten zien dan nodig is, en in dit geval geeft alleen de derde regel een
interessante uitkomst.

Opdracht 3d:

Maak nu zelf een soortgelijk M-file-programma en controleer of het werkt. Het mag
uiteraard uit meer dan 3 regels bestaan.

Het belangrijkste verschil tussen een Command Window (met daaraan gekoppeld de
workspace) en een M-file zal nu duidelijk zijn. Tevens is helder geworden waarom je ze
beide nodig hebt (resp. uitvoeren en schrijven van een programma). Om de volgende keer
dat je MatLab opstart je M-file weer terug te vinden, dien je op File  Open te klikken.
Dan vind je o.a. de bewaarde M-files.




MatLab                                         8                                          2-3-12
We keren tenslotte nog even terug naar het Command Window en we beschrijven hier de
mogelijkheden om je workspace op te slaan, danwel leeg te maken en om de variabelen
in je workspace weer te geven. Onder File “Save Workspace as” kun je de variabelen die
je tot dan toe in de workspace hebt gecreëerd bewaren, zodat je een volgende keer weer
verder kunt gaan. Je haalt de workspace weer op door op File  Open te klikken. Let op:
de mogelijkheid om met de pijltjestoetsen te scrollen begint in deze situatie met een
schone lei: de oude commando's gaan verloren. Je moet dus uit je hoofd weten wat je aan
instructies had gegeven. Dit nadeel is er de oorzaak van dat je deze functie in de praktijk
niet zo vaak zult gebruiken. Om het Command Window ‘op te schonen’, kan je het
commando ‘clc’ geven. Hiermee krijg je een leeg Command Window, maar blijven de
variabelen in je workspace bewaard (uiteraard vervalt ook nu de scrollmogelijkheid). Wil
je ook de bewaarde variabelen in je workspace verwijderen, kies dan onder Edit voor de
optie Clear Workspace. Een overzicht van de variabelen in je workspace krijg je tenslotte
met de commando's 'who' (alleen de namen) en 'whos' (gedetailleerde informatie).

Opdracht 3e:

Definieer een variabele (bijv. c: “a=2; <enter> b=3; <enter> c=a+b; <enter>”). Maak je
Command Window schoon en vraag de waarde van c op. Maak nogmaals je Command
Window schoon en verwijder opgeslagen variabelen uit je workspace om met een schone
lei verder te kunnen gaan…

3.4 Vectoren en matrices
In de voorbeelden tot nu toe hebben we gewerkt met enkelvoudige getallen, hetzij als
constanten zoals 3 of 4, hetzij als symbolen zoals a of b. MatLab maakt echter ook
uitgebreid gebruik van matrices en vectoren. (MatLab = matrix laboratory). Een matrix
is een rechthoekige ordening van getallen. Ze bevat een bepaald aantal rijen en een
bepaald aantal kolommen. Bijvoorbeeld:
 x = [1,     2, 3,      4,   5;
      1,     4, 9,     16, 25;
      1,     8, 27,    64, 125]

In dit voorbeeld is x een 3 bij 5 matrix. Bij het definiëren worden de getallen ingetikt
tussen rechte haakjes (engels: brackets). Rijen worden van elkaar gescheiden door
puntkomma's en de getallen binnen een rij door komma's of spaties. We hadden ook
mogen schrijven:
x = [1,2,3,4,5;1,4,9,16,25;1,8,18,32,125].

De eerste rij kan nog korter, nl. 1:5, wat je kunt lezen als 1 t/m 5.

Wanneer een matrix maar één rij bevat, spreken we ook wel van een rijvector; wanneer
een matrix maar één kolom bevat spreken we van een kolomvector. MatLab kent heel
veel bewerkingen op matrices. Eén van de meest simpele is transpositie: bij het
transponeren van een matrix worden rijen en kolommen verwisseld (de matrix wordt dus
als het ware gekanteld).


MatLab                                          9                                        2-3-12
In ons voorbeeld levert het commando:
         y = transpose(x)
of nog korter:
         y = x'
de volgende waarde voor y op:
  [ 1,      1,   1;
    2,      4,   8;
    3,      9, 27;
    4,     16, 64;
    5,     25, 125 ]

In hoofdstuk 4 gaan we nog wat nader in op de taal die MatLab gebruikt. Gezien de
beperkte omvang van deze cursus, gaat het te ver om uitgebreid in te gaan op de principes
van de matrixalgebra. Wie echter intensief met MatLab gaat werken, zal er niet aan
ontkomen om zich de matrixalgebra eigen te maken, aangezien vrijwel alle operaties die
je in MatLab uitvoert hierop gebaseerd zijn.

3.5      De help-functie van MatLab
We zijn nu aangekomen bij een handige functie van MatLab, namelijk de help-functie.
Hierin staan alle commando’s en functies die MatLab kent genoemd en uitgelegd. In het
begin zal deze uitleg erg technisch overkomen, maar na goed en aandachtig lezen valt er
veel nuttige informatie uit af te leiden. We gaan eens een kijkje nemen in de
helpbestanden. Klik rechtsboven Help aan. Er verschijnen nu een aantal mogelijkheden.

De bovenste mogelijkheid is Full Product Family Help. Dat is een functie waar je veel
mee te maken zult krijgen. Onder het tabblad Contents zijn schematisch alle onderdelen
van MatLab te vinden. Je kunt door te dubbelklikken op het door jou gekozen onderdeel
verder in de schematische structuur terechtkomen, net zolang tot je bij dat onderdeel
aankomt waar je vragen over hebt. Dit wandelen door de structuur van MatLab is een vrij
omslachtige manier van zoeken. Efficiënter is het om via het tabblad Index op trefwoord
te zoeken.

Opdracht 3f:

Ga naar het Help Window en zoek ‘fft’. Zoek eerst via het tabblad Contents (fft staat
voor ‘fast fourier transform' en dus zal het wel te vinden zijn bij de signal processing
toolbox). Zoek vervolgens via het Tabblad Index. Veronderstel ook eens dat je de term
‘fft’ niet kent, maar op zoek bent naar een functie die een Fourier-analyse uitvoert.

De tweede, derde en vierde mogelijkheid onder het menu Help, zijn MATLAB Help,
Using the Desktop en Using the Command Window. Door één van deze mogelijkheden
aan te klikken, kom je in hetzelfde venster als wanneer je voor de optie Full Product
Family Help kiest. Het verschil is echter dat je nu in de schematische structuur direct bij
respectievelijk MATLAB Help, Using the Desktop en Using the Command Window
terechtkomt. De laatste relevante mogelijkheid onder Help is Demos.




MatLab                                        10                                         2-3-12
Opdracht 3g:

Probeer, om een indruk te krijgen van de mogelijkheden van MatLab, een aantal demo’s
uit.

Gebruik maken van het menu-item Help is één manier om de helpbestanden te bereiken.
Een druk op de “?”-knop op de toolbar is een tweede manier en leidt tot hetzelfde
resultaat als de keuze voor de optie Full Product Family Help.

Een derde manier om hulp te krijgen, is via commando’s in het Command Window. Met
de commando’s help en helpwin gevolgd door een onderwerp dat in de helpbestanden
van MatLab bekend is, roep je respectievelijk een helpbestand op in het Command
Window of in een apart venster. De voorkeur gaat uiteraard uit naar de laatste van de
twee mogelijkheden, want je wilt je Command Window zo overzichtelijk mogelijk
houden.

Opdracht 3h:

Probeer nu beide mogelijkheden eens uit. Gebruik hiervoor weer 'fft' als het onderwerp
waarover je meer wilt weten.

De commando’s help en helpwin geven een korte beschrijving van de functie en een korte
uitleg over de te gebruiken syntax. Wil je echter meer gedetailleerde informatie over een
functie, dan kun je het commando doc gevolgd door een onderwerp gebruiken.

Opdracht 3i:

Vraag om hulp voor de functie ‘plot’. Doe dit zowel met het commando helpwin als met
het commando doc en vergelijk de twee helpbestanden.

Met de commando’s help, helpwin en doc kun je alleen zoeken op onderwerpen waarvan
je de naam (zoals MatLab die hanteert) kent. In sommige gevallen ken je deze naam
echter niet en voor die gevallen biedt het commando lookfor uitkomst. Met dit commando
gevolgd door een trefwoord, laat je MatLab in de eerste regel van elk helpbestand zoeken
naar het door jou opgegeven trefwoord. Voor iedere treffer toont MatLab de eerste regel
van het helpbestand. M.b.v. het help of doc commando kun je dan verder zoeken.
Wanneer je het commando lookfor -all (gevolgd door een onderwerp) geeft, doorzoekt
MatLab niet alleen de eerste regel van de helpbestanden, maar de gehele inhoud van de
helpbestanden.

Opdracht 3j:

Zoek naar informatie over Fourier-analyses, met als trefwoord ‘fourier’.

Tip: Hoewel de helpfunctie van MatLab veel informatie geeft, is het soms misschien
handig een elementaire handleiding voor het werken met MatLab bij de hand te hebben.
Deze is te vinden op:
http://www.mathworks.com/academia/student_center/tutorials/launchpad.html


MatLab                                        11                                        2-3-12
4        Taal-elementen
In het voorgaande hebben we kennisgemaakt met de MatLab-taal. In dit hoofdstuk zetten
we de zaken wat meer systematisch op een rij.

4.1      Objecten
De belangrijkste objecten waarmee MatLab werkt zijn:

Matrix:        Een rechthoek met getallen
Vector:        Een kolom met getallen; een n bij 1 matrix dus
Rij-vector:    Een rij met getallen; een 1 bij n matrix dus
Scalar:        Eén getal; een 1 bij 1 matrix dus

De namen van deze objecten moeten beginnen met een letter en de eerste 31 tekens
moeten uniek zijn. Let erop, dat MatLab onderscheid maakt tussen hoofdletters en kleine
letters. Objecten worden automatisch aangemaakt. Na het commando 'a = [1:10]'
bijvoorbeeld, is het object a een rijvector met 10 elementen. Het object blijft in deze
vorm bestaan totdat het opnieuw gedefinieerd wordt of opgeruimd met het commando
'clear a'. Bij de definitie scheiden spaties of komma's de elementen van een rij en
scheiden puntkomma's de rijen: X = [1,2,3;4,5,6] definieert deze matrix:        1 2 3
                                                4 5 6

Getallen worden opgeslagen met een precisie van ongeveer 16 decimale cijfers en een
bereik van ongeveer 10-308 t/m 10+308

Tot nu toe hebben we al regelmatig gebruik gemaakt van de dubbele punt om een
rijvector aan te geven. De dubbele punt speelt echter nog meer rollen. De volgende
voorbeelden maken het duidelijk:

Notatie             Betekenis
1:4                 de rijvector [1,2,3,4]
1:2:7               de rijvector [1,3,5,7]; de 2 geeft een stapgrootte aan
100:-7:50           de rijvector [100,93,86,79,72,65,58,51]
0:pi/4:pi           0, 0.7854, 1.5708, 2.3562, 3.1416
A(1:k,j)            de eerste k elementen van de j-de kolom van A
A(:,j)              alle elementen van de j-de kolom van A

Opdracht 4a:

Probeer eens uit wat de volgende reeks commando's oplevert.
X =   [1:3; 1:2:5]
V =   X(:,2)
Q =   X(:,end)
W =   V'
Num   = W(1)




MatLab                                          12                                   2-3-12
4.2      Speciale functies om matrices te vullen
Een functie is een bewerking die wordt aangegeven met een naam, meestal gevolgd door
één of meer argumenten tussen haakjes. Denk aan vertrouwde (?) vormen als sin(x),
cos(x), max(x,y) enz. Je kunt het resultaat meteen toewijzen aan een object of het
gebruiken in een expressie, bijvoorbeeld: z = cos(x)*sin(y).

Er zijn een aantal speciale functies, die het definiëren en vullen van matrices en vectoren
vergemakkelijken:
Z = zeros(2,4)                levert een 2 bij 4 matrix met nullen
W = ones(3,3)                 geeft een 3 bij 3 matrix met enen
F = 5*ones(3,3)               F is natuurlijk een 3 bij 3 matrix met in elke cel het getal 5
R = rand(1,10)                Rij met 10 trekkingen uit een uniforme verdeling (0  x < 1)
R = 6*rand(1,10))             Rij met 10 trekkingen uit een uniforme verdeling (0  x < 6)
R = fix(6*rand(1,10))         Rij met 10 trekkingen uit een uniforme verdeling,
                              afgekapt tot hele getallen (0  x  5)
R = 1+fix(6*rand(1,10)) Rij met 10 uitkomsten van een nagebootste dobbelsteen
R = randn(4,4)                R is een 4 bij 4 matrix met random trekkingen uit een
                              standaard normale verdeling. Handig voor simulaties van
                              normaal verdeelde variabelen!

Het verwijderen van een rij, kolom of cel gebeurt door er een lege verzameling aan toe te
kennen: 'Z(:,2) = []' verwijdert de tweede kolom uit Z.

Je kunt ook objecten aan elkaar knopen (concatenate): 'RR = [R,R*2;R*2,R*4]'.

Tip:
Als je in een simulatie iets wil laten gebeuren met kans p kun je dat als volgt doen:
if rand < p
  .... laat het gebeuren
end




MatLab                                         13                                         2-3-12
4.3      Operatoren en functies
Waarden worden aan objecten toegekend via expressies (zeg maar formules):
object = expressie. Deze expressies gebruiken operatoren en functies. De objecten die
binnen de expressie genoemd worden noemt men operands. De meest voorkomende
operatoren zijn:
+     Optellen van scalars, matrices of beiden; de som van twee matrices wordt berekend
      door corresponderende elementen op te tellen; ze moeten natuurlijk dezelfde maten
      hebben. De som van een matrix en een scalar wordt verkregen door bij elk element
      van de matrix de scalar op te tellen.
-     Aftrekken van scalars, matrices of beiden; deze operator werkt analoog aan de
      vorige.
*     Vermenigvuldiging van scalars, matrices of beiden. Als beide operands een matrix
      of een vector zijn, wordt een matrix-vermenigvuldiging uitgevoerd als gedefinieerd
      door de lineaire algebra. Dat is heel iets anders dan een element bij element
      vermenigvuldiging.
.*    Dit is een element bij element-vermenigvuldiging. De operands moeten dezelfde
      maten hebben, tenzij er één een scalar is.
/     Rechter matrix-deling. A/B komt grofweg neer op A*INV(B)
./    Element bij element deling. A./B betekent, dat elk element van A wordt gedeeld
      door het corresponderende element van B. A en B moeten dezelfde maten hebben of
      B moet een scalar zijn.
\     Linker matrix-deling. A\B komt grofweg neer op INV(A)*B. A en B moeten dezelfde
      maten hebben of A moet een scalar zijn. Als A een m bij n matrix is met m  n en B
      een kolomvector met m elementen, dan geeft X = A\B de kleinste kwadraten-
      oplossing voor de set vergelijkingen AX = B.
.\    De element bij element linker deling: A.\B betekent, dat ieder element van B wordt
      gedeeld door het corresponderende element van A. Beide operands moeten dezelfde
      maat hebben, tenzij de eerste scalar is.
^     Machtsverheffing van matrices of scalaren. A^p met een positief geheel getal p
      betekent, dat A p-1 keer met zichzelf wordt vermenigvuldigd. A moet dan natuurlijk
      vierkant zijn: A^2 is hetzelfde als A*A. Als p een negatief geheel getal is wordt de
      inverse van A een aantal keren met zichzelf vermenigvuldigd: A^-2 =
      INV(A)*INV(A). Als p een niet-geheel getal is, wordt het nog ingewikkelder.
.^    Element bij element machtsverheffing. Elk element in de eerste operand wordt tot
      de macht verheven, die in het corresponderende element van de tweede operand
      wordt aangegeven. Beide operands moeten dezelfde maten hebben of de tweede
      moet scalar zijn.
'     Transpositie: A' levert een gekantelde matrix op: rijen zijn kolommen geworden en
      omgekeerd. Als bijvoorbeeld A = [1 2 3; 4 5 6] dan is A' = [1 4; 2 5; 3
      6].




MatLab                                        14                                        2-3-12
Veel commando's kun je op twee manieren uitvoeren:

1 als commando, dat wil zeggen in de vorm operatie operand(s),
   bijvoorbeeld:                    Load MijnWerk
2 als functieaanroep, bijvoorbeeld: bestandsnaam = input('welk        bestand: ');
                                     Load('bestandsnaam');

In die laatste vorm wordt de waarde gebruikt van het object, dat je als argument
meegeeft, niet zijn naam! Het voordeel is, dat je nu niet tijdens het schrijven van het
programma al hoeft te weten hoe het bestand heet. Met de opdracht bestandsnaam =
input('welk bestand: '); vraag je tijdens het uitvoeren van het programma om een
naam. Die naam wordt de waarde van bestandsnaam. Het programma geeft daarna het
argument 'bestandsnaam' door aan de functie. De conversatie verloopt dus als volgt:
welk bestand: 'mijnspul.mat'
waarna alle variabelen uit mijnspul.mat zijn binnengehaald in de workspace.

4.4      Eigen functies
Naast de functies, die standaard in MatLab zitten, kun je ook nog je eigen functies
definiëren. Ze moeten dan in een bestand komen met de extensie '.m'. De naam van dat
bestand moet overeenkomen met de naam van de functie. De eerste regel van het bestand
moet aangeven hoe de functie aangeroepen moet worden; daarna volgen één of meer
commentaarregels die een korte beschrijving geven. Deze beschrijving komt voor de dag
als iemand het commando 'help functie' intikt. Die uitleg eindigt bij de eerstvolgende
regel, die niet met een percentageteken begint.

Let op: de percentagetekens moeten in de eerste positie van de regel staan zonder
        spaties ervoor.

Een voorbeeld: het bestand meansd.m zou er zo uit kunnen zien:

function [mean,stdev] = meansd(x)
% [m,s]= meansd(x): m and s give mean and stand. deviation of vector x.
 n = size(x,1)
 mean = sum(x)/n;
 stdev = sqrt(sum((x - mean).^2)/n);


De namen die tussen het woord function en het =teken staan (hier mean en stdev) zijn de
lokale namen voor de waarden die de functie berekent. Bij het gebruik van de functie
worden ze vervangen door de namen in de aanroep, bijvoorbeeld:
         [gem,std] = meansd(data);
De    naam achter het =teken (hier meansd)is de naam die gebruikt moet worden om de
functie op te roepen. Daarachter staan tussen haakjes de argumenten (hier x), dat wil
zeggen de waarden die de functie nodig heeft voor haar berekening.




MatLab                                       15                                      2-3-12
De variabelen binnen deze functie zijn lokaal: ze hebben dus niets te maken met de
variabelen in de workspace van waaruit de functie wordt aangeroepen. Ook de
argumenten (x in het voorbeeld) zijn lokale namen voor de doorgegeven objecten. Zoals
je ziet, kan de functiewaarde bestaan uit een hele verzameling objecten. In de aanroep
moet je dan natuurlijk ook het juiste aantal argumenten meegeven.

INLEVEROPDRACHT 4.b:

Schrijf een functie voor het berekenen van de t-statistic voor het vergelijken van de
gemiddelden van twee groepen. Argumenten zijn twee vectoren x1 en x2 met de scores
van de proefpersonen in elk van de groepen. De formule voor de t-statistic is:
                  x1  x 2
        t=
                   2
                  s1 s 2
                       2
                  n1 n 2
waarbij si de standaarddeviatie van xi is en ni het aantal elementen in groep i.

4.5      Keuzes en herhalingen
In een m-file kun je de volgorde van uit te voeren bewerkingen regelen met een aantal
speciale commando's of, zoals dat in programmeertalen heet: statements. Ze maken
gebruik van speciale sleutelwoorden (keywords). We zullen ze hier kort bespreken.

if ... elseif ... else ... end

Met het if-statement geef je aan, dat een groep van statements wel of niet moet worden
uitgevoerd afhankelijk van een conditie. De conditie wordt aangegeven met behulp van
een logische expressie, d.w.z. een uitdrukking die ja ( niet nul) of nee (nul) oplevert. Met
de keywords elseif en else kun je alternatieve groepen aangeven. Ieder blok dat begint
met if moet worden afgesloten met het keyword end. Hieronder een eenvoudig
voorbeeld:
function maximum = grootste(a,b,c)
% grootste(a,b,c) geeft de grootste van de waarden in a, b en c
if a >= b & a >= c maximum = a;
elseif b >= c maximum = b;
else maximum = c;
end


In de logische expressies worden vergelijkingsoperatoren gebruikt:
>       groter dan
>=      groter dan of gelijk aan
<       kleiner dan
<=      kleiner dan of gelijk aan
==      gelijk aan
~=      ongelijk aan




MatLab                                         16                                         2-3-12
Ook worden enkele logische operatoren gebruikt:
&      én
|      of
~      niet
en natuurlijk allerlei functies die logische waarden teruggeven.

Je moet er goed aan denken dat MatLab een matrix-taal is: als a en b scalars zijn levert
a == b gewoon 1 of 0 op, maar als het matrices zijn is het resultaat een matrix van enen
en nullen. In het statement "If a == b doeiets", zal de opdracht doeiets alleen
worden uitgevoerd als alle elementen van a gelijk zijn aan de corresponderende
elementen van b. Als je wil weten of de matrices gelijk zijn moet je dit schrijven:
isequal(a,b). Nog een paar handige functies zijn in dit verband:
any(x)        1 als minstens één element van x niet nul is
all(x)        1 als geen enkel element van x nul is
isempty(x) 1 als x nul elementen bevat

switch en case

Met het switch-statement kies je één groep statements uit een aantal groepen en wel op
basis van de waarde van een variabele of expressie. De groepen worden aangegeven met
de woorden case en otherwise. De hele constructie eindigt met een end-statement. Een
voorbeeld:
function resultaat = reken(a,op,b)
% reken(a,op,b) voert operatie op ('*','/','+' of '-') uit
% met operands a en b
switch op
case '*'; resultaat = a .* b
case '/'; resultaat = a ./ b
case '+'; resultaat = a + b
case '-'; resultaat = a - b
otherwise; disp(['ongeldige bewerking: ',op])
end


for

De for-constructie wordt gebruikt om een groep statements te herhalen, waarbij één
bepaalde scalar steeds een andere waarde aanneemt. De groep statements wordt
afgesloten met het keyword end. De volgende functie bijvoorbeeld gebruikt twee geneste
for-loops om een matrix element voor element na te lopen.
function zoveel = count(mat,x)
% count(x,getal) telt in matrix mat het aantal elementen met waarde x
zoveel = 0;
[nrow,ncol] = size(mat);
for i = 1:nrow;
  for j = 1:ncol;
    if mat(i,j) == x zoveel = zoveel + 1;
    end
  end
end



MatLab                                         17                                     2-3-12
while

De while-constructie geeft na het woordje while een logische expressie. Wanneer deze
expressie waar (lees niet nul) oplevert wordt de hele reeks reeks statements na de
expressie en vóór het bijbehorende end-statement uitgevoerd. Daarna wordt de expressie
opnieuw ge-evalueerd; wanneer de expressie nog steeds waar oplevert, wordt de reeks
statements opnieuw uitgevoerd enzovoort, totdat de expressie onwaar (lees nul) oplevert.
x = 0
while x < 1 | x > 2;
  x = input(' geef het aantal dimensies (1 of 2): ');
end


continue en break

Met het continue-statement binnen een for- of while-loop geef je aan, dat de lopende
herhaling niet hoeft te worden afgemaakt, maar dat meteen met de volgende moet worden
begonnen.
Met het break-statement geef je ook aan, dat de lopende herhaling moet worden
afgebroken, maar nu moet meteen de loop worden verlaten. Bij geneste loops wordt
alleen de binnenste verlaten. Dit statement gebruik je bijvoorbeeld in een iteratief proces
wanneer je het doel van de iteratie bereikt hebt.

INLEVEROPDRACHT 4c:

Analyseer onderstaande m-file en vind het antwoord op de volgende vragen:
   1. Wat krijg je te zien als je het volgende intikt: help pieken
   2. Wat is de rol van de variabele stapgrens
   3. Het for-commando vormt het begin van een iteratieproces; bij welke regel houdt
      het te herhalen deel op.
   4. Wat doen al die end-statements op het einde?




MatLab                                        18                                         2-3-12
function [toppen,waarden,ntop] = pieken(signaal,drempel);
% function [toppen,waarden,ntop] = pieken(signaal,drempel);
%
% pieken(signaal,drempel)
% berekent plaatsen en waarden van de pieken in signaal en telt ze
% Input : signaal (1:n)
% Drempel: drempelwaarde in stand. dev. boven gemiddelde:
%          waarden daaronder zijn nooit toppen
% Output : toppen(1:ntop)= indices van extreme waarden
%               waarden(1:ntop) = bijbehorende extreme waarden
%          ntop = aantal pieken
%
s = signaal;
n = length(s);
grens = mean(s)+ drempel.*std(s);
stapgrens = 2;
volgnummer = 0;
for i = 2:n-1,
  if (s(i) >= s(i-1)) & (s(i) >s (i+1)) & (s(i) > grens)
    index = i;
    if volgnummer == 0
      volgnummer = 1;
      toppen(volgnummer) = index;
      waarden(volgnummer) = s(index);
    end;
    if (volgnummer >= 1) & (index - toppen(volgnummer) > stapgrens)
      volgnummer = volgnummer + 1;
      toppen(volgnummer) = index;
      waarden(volgnummer) = s(index);
    elseif (volgnummer > 1) & (index - toppen(volgnummer) <= stapgrens)
      if s(index) > waarden(volgnummer)|
         toppen(volgnummer) = index;
         waarden(volgnummer)= s(index);
      end;
    end;
  end
end
ntop = volgnummer;


4.6      Opslaan en opvragen van bestanden
De inhoud van je workspace kun je in een file opbergen met het commando 'save naam'.
Voor naam moet je dan een file-naam geven. Als die geen extensie bevat zet MatLab er
standaard '.mat' achter. Als je niet alle variabelen wilt opslaan, kun je na de naam nog een
lijst met variabelen-namen opgeven. Alleen de genoemde variabelen worden dan
opgeslagen. Als de file-naam in een (string-)variabele staat, bijvoorbeeld omdat je hem
hebt opgevraagd, kun je de functie-achtige vorm van het save-commando gebruiken,
bijvoorbeeld: save('bestandsnaam','a','b') om de variabelen a en b op te slaan in
de file waarvan de naam in bestandsnaam staat.




MatLab                                         19                                         2-3-12
De opgeslagen informatie wordt door MatLab in binaire vorm (d.w.z. in zijn interne
representatie) opgeslagen. Editors en tekstverwerkers weten daar doorgaans geen raad
mee. Als je de informatie meer leesbaar wilt opslaan kun je aan het einde van het
commando nog '-ASCII' of '-ASCII -DOUBLE' toevoegen. Getallen worden dan in ASCII
formaat (8 of 16 decimale cijfers) opgeslagen1. Daarmee kun je natuurlijk wel enige
precisie kwijt raken.

Een eenmaal opgeslagen bestand kun je opvragen met het commando 'load naam'. De
naam kan natuurlijk het volledige pad aangeven. Als dat niet zo is, zoekt MatLab in de
current directory en in zijn eigen search path.

De current directory kun je instellen via de balk boven het command window. De
standaard-instelling kun je aanpassen door (buiten MatLab om) rechts te klikken op de
snelkoppeling naar MatLab en dan te klikken op de knop Eigenschappen. Achter de tekst
'Beginnen in' vind (of zet) je de current directory.

Het search path is een lijst van directories. Als een bestand niet met zijn volledige pad
wordt aangegeven, zal MatLab alleen zoeken in de directories, die in dit pad staan.
Subdirectories zijn niet automatisch inbegrepen. Je kunt het pad opgeven door in het file-
menu 'Set Path' te kiezen of door het commando 'PathTool' in te tikken.

Je kunt ook nog een tussenweg bewandelen tussen volledig pad en alleen maar de naam.
Je geeft dan alleen het laatste deel van het pad. Je moet dan wel in plaats van backslashes
slashes gebruiken, bijvoorbeeld: 'TestBestanden/test1.mat'. Als de naam spaties bevat
moet je de functionele vorm gebruiken: load('naam').

Na de bestandsnaam mag je nog variabelen noemen, bij voorbeeld: load mijnwerk x y.
In dat geval zal MatLab alleen de opgegeven variabelen (x en y) ophalen. Die moeten dan
natuurlijk wel in het bestand zijn opgeslagen. Je kunt in die lijst van variabelen ook de
joker '*' gebruiken: bijvoorbeeld: item*b voor alle variabelen waarvan de naam begint
met item en eindigt met b.

Je kunt de inhoud van het bestand (of de opgegeven selectie) ook toekennen aan een
object, bijvoorbeeld S = load('MijnArray').

Als in het load-commando een bestandsnaam geen extensie bevat zoekt MatLab naar een
file met de opgegeven naam en, als dat niets oplevert, naar een file met de naam
naam.mat. In beide gevallen wordt aangenomen, dat het om een bestand in MatLabs
eigen binaire vorm gaat. Als de naam een andere extensie dan .mat bevat wordt
aangenomen, dat het om een ASCII-bestand gaat (dus om makkelijk leesbare tekst).

Met 'load -mat naam ' of 'load -ascii naam' geef je zelf het type bestand aan; de
extensie doet er dan niet meer toe.



1
 Dit soort vlaggen, beginnend met een min-teken, zul je nog wel meer tegenkomen. Ze passen typisch in
de UNIX-stijl.


MatLab                                             20                                              2-3-12
Hoe je m-files bewaart is al bekend. Nog even een tip: bewaar je eigen m-files nooit in
$matlabroot/toolbox/matlab. Ze kunnen worden overschreven als je een nieuwe versie
van MatLab installeert. Bovendien werken wijzigingen in die directory vertraagd door:
om snelheid te winnen laadt MatLab hun adressen al bij het opstarten in het
werkgeheugen.

Als je niet het load-commando gebruikt, maar zomaar een naam intikt (bijvoorbeeld
'Blaf') verloopt het zoekproces als volgt:
1. MatLab zoekt naar een variabele met de naam 'Blaf'.
2. Als die er niet is zoekt MatLab naar een 'built-in functie met die naam.
3. Als die ook niet wordt gevonden zoekt MatLab in de current directory naar het
     bestand 'Blaf.m'.
4. Als dat niets oplevert zoekt MatLab in de directories van het search path naar een file
     'Blaf.m' en kiest het eerste gevonden bestand.
Als je op meerdere plaatsen bestanden hebt staan met dezelfde naam kun je met het
commando 'Which naam' achterhalen welk door MatLab gekozen is.




MatLab                                        21                                        2-3-12
5        Signaalopbouw
In dit hoofdstuk zullen we een sinusvormig signaal opbouwen. In het verlengde hiervan
zul je meteen kennismaken met de grafische mogelijkheden van MatLab. In het begin zal
er nog veel voorgekauwd worden, maar naarmate we verder komen in de stof zal meer en
meer het initiatief en de vindingrijkheid van jou moeten komen. Daarbij kun je gebruik
maken van je kennis over de structuur en functies (denk aan Help) van MatLab en van je
inhoudelijke kennis over signaalverwerking (zoals behandeld in het theoretisch gedeelte
van de cursus).




Figuur 3:     Een signaal wordt gesampled

Als een signaal wordt gemeten, gebeurt dat door op regelmatige tijden het onderliggende
fysieke kenmerk te meten (bijvoorbeeld de luchtdruk op een microfoon). Het aantal
metingen (samples) per seconde heet de sampling frequentie (in Hz: lees Hertz). De
reeks getallen die dat oplevert vormt een beschrijving van het signaal en bij voldoende
hoge sampling frequentie is die reeks voldoende om het signaal te reproduceren. Volgens
dat principe werken alle digitale media. Het signaal uit figuur 3 zou dus kunnen worden
opgeslagen als: [0 0.16 0.32 0.46 0.60 0.72 0.82 0.90 0.96 1.0 0.98 .....
enz.

5.1      Een sinusfunctie
We beginnen bij het begin. Een sinusvormig signaal is een signaal, waarvan de sterkte
(de amplitude, bijvoorbeeld de luchtdruk op je oor) beschreven kan worden als een
sinusvormige functie van de tijd.




MatLab                                      22                                        2-3-12
Allereerst zullen we zien hoe we zo'n sinusvormig signaal kunnen programmeren en
vervolgens hoe we het grafisch kunnen weergeven. Hiervoor zijn meteen een tweetal
MatLab-functies nodig. De eerste is linspace. Je gebruikt linspace om een fatsoenlijke
tijdsas te maken: een rijvector met een indeling van de tijd. De drie parameters achter
linspace staan achtereenvolgens voor het begin van de tijdsas (in seconden), het einde van
de tijdsas (eveneens in seconden) en het aantal punten waarin je de tijdsas onderverdeelt.

De tweede MatLab-functie is de zogenaamde plot-functie. Die heb je nodig om een
afbeelding van het signaal te maken. De parameters die je achter plot aangeeft zijn de x-
as, de y-as, en eventueel het teken (bijv ‘o’) waarmee de data-punten in de figuur worden
weergegeven. Voer je deze laatste parameter niet in dan maakt MatLab automatisch een
lijn.

Opdracht 5a:

Typ in een M-File (noem hem bijvoorbeeld sinus1) het volgende in en draai het
programma vervolgens.
t= linspace (0,1,101); % tijdsas
a= sin(t);             % a = bijbehorende amplitudes
plot(t,a)              % plaatje van het signaal


Kun je ontdekken wat er aan de hand is met de grafiek. Waarom is er geen sinus
zichtbaar zoals je die kent?

Wanneer je goed kijkt naar de functie sin(t) zoals we die hebben gedefinieerd, dan is het
logisch dat we slechts een deel van de sinus in beeld kregen. De periode van een
sinusgolf is nu eenmaal 2*pi.

Opdracht 5b:

Typ in een M-File (denk aan de logische naam) het volgende in en controleer of er nu wel
een hele sinusgolf tevoorschijn komt.
t= linspace (0,1,101); % t = tijds-as
a= sin(2*pi*t);        % a = amplitudes
plot(t,a)              % plaatje


Dit sinusvormig signaal heeft een frequentie van 1Hz (één complete golf in één seconde).
We kunnen nu weer een stapje verder gaan en een andere frequentie invoeren. Met deze
toegevoegde frequentie-parameter geven we aan hoe vaak de sinusgolf in het door ons
reeds gedefinieerde tijdsbestek van één seconde moet voorkomen.

Opdracht 5c:

Typ in een M-File het volgende in en kijk wat er gebeurt als je het programma draait.
d= 1;                         % d = duur in seconden




MatLab                                        23                                         2-3-12
sf = 100;                     %   sf = sampling frequentie
t= linspace (0,d,d*sf);       %   t = tijds-as
f= 2;                         %   f = frequentie in Hertz
a= sin(2*pi*f*t);             %   a = amplitude
plot(t,a)                     %   plaatje


De variabele sf in opdracht 5c definieert de sampling frequentie. Deze geeft aan
hoeveel metingen (getallen) er worden gebruikt om één seconde te beschrijven.

Probeer nu eens een andere waarde voor f uit en kijk vervolgens wat er gebeurt.

We hebben tot dusver nog nauwelijks gekeken naar de amplitude. We kunnen het signaal
sterker of zwakker maken via een nieuwe parameter v (volume). Laten we een nieuwe
sinusvorm programmeren met deze nieuwe parameter.

Opdracht 5d:

Typ het volgende in en kijk wat er gebeurt. Sla het uiteraard op in een nieuwe M-file.
d= 1;                         %   d = duur in seconden
sf = 100;                     %   sf = sampling frequentie
t= linspace (0,d,d*sf);       %   t = tijds-as
f= 2;                         %   f = frequentie in Hertz
a= sin(2*pi*f*t);             %   a = amplitude
v= 0.3;                       %   v = volume
a2= v*sin(2*pi*f*t);          %   a2 = kleinere amplitude
plot(t,a,t,a2)                %   plaatje van de twee functies


Zoals je ziet, is het heel goed mogelijk om in één plot twee verschillende signalen weer te
geven.

We gaan nu kijken wat er gebeurt als we twee sinusvormigen bij elkaar optellen. Om
twee plaatjes samen weer te geven moeten we gebruik maken van de functie subplot.
Deze functie verdeelt het figuurscherm in een tabel met een aantal cellen. Achter het
subplot-commando specificeer je daarom de grootte van de tabel (uitgedrukt in een
parameter voor het aantal rijen, gevolgd door een parameter voor het aantal kolommen)
en de plaats waar het betreffende plaatje in deze tabel moet staan (uitgedrukt in een
“plaats-parameter”, waarbij de cellen van links naar rechts en vervolgens van boven naar
beneden geteld worden, de cel linksboven heeft dus plaatsnummer 1).

Opdracht 5e:

Typ in een nieuwe M-file het volgende in. Probeer voordat je het resultaat bekijkt, het
eerst eens te voorspellen.




MatLab                                        24                                         2-3-12
d= 1;                    % d = duur in seconden
sf = 100;                % sf = sampling frequentie
t= linspace (0,d,d*sf); % t = tijds-as
f= 2;                    % f = frequentie eerste sinus
a1= sin(2*pi*f*t);       % a = amplitude eerste sinus
v= 0.3;                  % volume
f2= 10;                  % frequentie tweede sinus
a2= v*sin(2*pi*f2*t);    % a2 = amplitude tweede sinus
a3= a1+a2;               % a3 = de som van de twee sinussen
subplot(2,1,1), plot(t,a,t,a2);
subplot(2,1,2), plot(t,a3)


Kwam je voorspelling uit?

We zijn nu op een punt aangekomen waarop alle relevante parameters de revue zijn
gepasseerd.

   Of toch niet? Eigenlijk hebben we een parameter vergeten, nl. de fase. Twee sinusvormige signalen
   kunnen dezelfde frequentie en dezelfde amplitude hebben en toch verschillen, doordat ze ten
   opzichte van elkaar in de tijd verschoven zijn. De fase wordt in de beschrijving opgenomen door bij
   het argument van de functie sin een waarde op te tellen. We zullen daar in deze inleiding verder
   geen aandacht aan besteden.

Laten we eens kijken of we naast de grafische representatie van het door ons opgebouwde
signaal ook een auditieve representatie kunnen maken. Hiervoor hebben we weer een
nieuw commando nodig, te weten sound. Sound geeft de computer de opdracht om
gegevens naar de geluidskaart te sturen. Je moet aan sound twee parameters doorgeven.
Allereerst moet je aangeven welke vector het signaal beschrijft en ten tweede wat de
sampling frequentie is.

Opdracht 5f:

Schrijf het volgende programma in een M-File en draai het. Zoals je ziet, zit er in de regel
die het geluid moet produceren tweemaal een pause-commando. Het getal dat tussen
haakjes achter het commando staat (de parameter) geeft aan hoeveel seconden er
gewacht moet worden voordat het volgende commando wordt uitgevoerd. Als je die
parameter weglaat (en ook de haakjes) blijft MatLab wachten totdat je een toets indrukt.
De drie frequenties die hier zijn gegeven beschrijven de tonen a, b, en cis.
d= 2;                                %   duur 2 seconden
sf= 20000;                           %   redelijke sampling frequentie voor geluid
t= linspace (0,d,d*sf);              %   tijdsas
fa= 440;                             %   frequentie van a
a= sin(2*pi*fa*t);                   %   signaal van a
fb= 493.883;                         %   frequentie van b
a2= sin(pi*2*fb*t);                  %   signaal van b
fcis= 554.365;                       %   frequentie van cis
a3= sin(pi*2*fcis*t);                %   signaal van cis
sound(a,20000);pause(2);             %   een a
sound(a2,20000);pause(2);            %   een b
sound(a3,20000);                     %   een cis




MatLab                                               25                                                  2-3-12
INLEVEROPDRACHT 5g:

Pas het bovenstaande programma zo aan dat je het begin van Vader Jacob kunt afspelen
met zuivere sinusvormigen. Houd er rekening mee, dat de tijdsduur per toon korter moet
worden (en dat dus ook de parameter bij het pause-commando verandert).

Het is nu tijd om dit hoofdstuk af te sluiten met een wat realistischer input dan een
zuivere sinusvorm. We gaan een stukje muziek binnenhalen om dit vervolgens te laten
horen en zien. We halen het signaal binnen met het commando wavread, waarmee een
wav bestand door MatLab kan worden gelezen. Let wel, hoewel dit nieuwe
muzieksignaal er heel anders uitziet dan de voorgaande sinusvormigen, is en blijft het te
beschrijven als een som van sinusvormigen. Dit is namelijk de crux van de Fourier-
analyse, die we verderop zullen bespreken.

Opdracht 5h:

Zoek op je PC naar een wav-file (bijvoorbeeld het openingsdeuntje van Windows). Typ
het volgende in een M-file en draai het programma vervolgens (uiteraard voer je bij c:\…
de naam in (inclusief het juiste pad) van de door jou gekozen file).
[y,sf] = wavread('c:\...');       %   het signaal
d = length(y)/sf;                 %   duur = aantal samples / sampling freq.
t=linspace (0,d,length(y));       %   tijds-as
sound (y,sf)                      %   geluid
plot(t,y)                         %   beeld


Probeer nu eens zelf een ander wav-bestand van het internet te halen en in je programma
te stoppen. Let erop dat het geluidsfragment niet te lang mag duren.

We zijn aan het einde gekomen van dit hoofdstuk. We gaan er in het vervolg van uit dat
de opbouw van een sinusvormig signaal en de daarmee samenhangende parameters (zoals
frequentie, volume en sampling frequentie) helder zijn.




MatLab                                       26                                        2-3-12
6        Fouten opsporen in M-files
Waarschijnlijk heb je in de vorige opdrachten al gemerkt, dat er bij het programmeren
geregeld een foutje in je M-file sluipt, waardoor de opdrachten in de file niet meer
uitvoerbaar zijn. Soms is zo’n foutje gemakkelijk op te sporen, maar vaak ook niet en in
die gevallen is het handig dat MatLab een debugger ter beschikking heeft. Deze
debugger helpt je om je M-file systematisch te doorzoeken op fouten en deze te
herstellen. In dit hoofdstuk zullen we eerst een indeling geven van het type fouten dat je
kunt maken en vervolgens voor elk van die fouten aangeven hoe je ze kunt opsporen en
herstellen.

6.1      Soorten fouten in M-files en hun opsporing
Grofweg kun je twee soorten fouten onderscheiden:

       Syntax fouten:     deze fouten kun je vergelijken met spelfouten (bijvoorbeeld
                           het verkeerd spellen van een functie of het vergeten van een
                           haakje/aanhalingsteken). Deze fouten zorgen ervoor dat een
                           M-file niet uitgevoerd kan worden.

       Run-time fouten: deze fouten zijn meestal algoritmisch van aard (bijvoorbeeld
                         het gebruiken van de verkeerde variabele of het geven van een
                         onjuiste formule voor een berekening). Deze fouten worden
                         zichtbaar, doordat MatLab onverwachte resultaten produceert.

Syntax fouten zijn meestal vrij eenvoudig op te sporen. Wanneer je in het Command
Window de opdracht geeft om een M-file te runnen en er zit een syntaxfout in deze M-
file, dan verschijnt er een foutmelding. De foutmelding beschrijft de fout en geeft ook de
regel aan waarin de fout zich bevindt. Door te klikken op het onderstreepte deel van de
foutmelding, wordt de betreffende M-file geopend. De cursor wordt dan automatisch
geplaatst op de regel waar de fout ontdekt is. Dat hoeft niet echt de plek van de fout te
zijn; als je bijvoorbeeld een openhaakje vergeet, kun je tegen een foutmelding aanlopen
bij het (niet) corresponderende sluithaakje.

Opdracht 6a

Typ in een nieuwe M-file het volgende programmaatje, spoor de fout die erin zit op en
herstel deze.
t= linspace (0,1,101); % t = tijds-as
a= sin(2*pi*t);        % a = amplitudes
plot(t;a)              % plaatje


Zoals je wellicht gemerkt hebt, geeft MatLab wel een aanwijzing waar de fout ontdekt is,
maar geeft het programma niet altijd aan hoe je de fout moet herstellen. Je zult dus in
veel gevallen alsnog de helpfunctie moeten gebruiken om de juiste syntax te vinden.



MatLab                                        27                                        2-3-12
Run-time fouten zijn moeilijker om op te sporen. Er zijn een aantal manieren om dit te
doen: handmatig of met behulp van de debugger. In de volgende paragraaf zullen we
ingaan op het gebruik van de debugger; in deze paragraaf eerst nog een tip om fouten
handmatig op te sporen (dit kan in sommige gevallen sneller gaan dan het doorlopen van
het gehele debugproces). Zoals gezegd, herken je run-time fouten aan onverwachte
resultaten. Om te onderzoeken hoe vroeg in de berekeningen er iets misgaat, kun je de
“;”-tekens die je gegeven hebt aan het einde van elke regel in je M-file verwijderen. De
tussenresultaten worden nu niet onderdrukt, maar weergegeven als je de file runt. De fout
zit uiteraard vóór het eerste onverwachte tussenresultaat. Mocht deze simpele
opsporingsmethode er niet toe leiden dat je de fout vindt, dan kun je de debugger
inschakelen.

6.2      Het gebruik van de debugger
De debugger kan op twee manieren aangestuurd worden: in het venster van de M-file
(met de opdrachten die te vinden zijn onder het kopje “debug”) en in het Command
Window (hier moet je de commando’s zelf intypen). De methoden kunnen door elkaar
gebruikt worden, maar in deze paragraaf zal alleen het debuggen in de M-file aan de orde
komen.

Om de werking van de debugger verder uit te leggen, zullen we gaan werken aan de hand
van een voorbeeld. In dit voorbeeld wordt gebruik gemaakt van de wortel-functie.

      Deze functie berekent op een primitieve manier de wortel uit een getal (mits dat niet
      al te groot of te klein is). De berekening is gebaseerd op de volgende constatering:

          1           =   12
          1+3         =   22
          1+3+5       =   32
          1+3+5+7     =   42
          enzovoort

      Als je nu de wortel uit 16 wil berekenen hoef je alleen maar te kijken hoeveel van
      deze termen 1, 3, 5, 7, ... je nodig hebt om bij 16 te komen.

      Deze berekening is wel erg grof. Je kunt ze verbeteren door het getal eerst met
      bijvoorbeeld één miljoen te vermenigvuldigen, dan de wortel te berekenen op
      bovenstaande manier en vervolgens het resultaat door duizend te delen. De functie
      wortel vermenigvuldigt het getal met één miljoen, roept dan de functie wortel1 voor
      de berekening via bovenstaande methode en deelt de uitkomst door 1000.

Helaas (?) is er een foutje geslopen in de m-files waarin deze functies zijn gerealiseerd.
De m-files zien er nu als volgt uit:

function [wortel] = wortel(x)
% y = wortel(x) trekt de wortel uit x mits x niet al te klein of groot
wortel = wortel1(x*1000000)/1000;




MatLab                                         28                                         2-3-12
function [wortel1] = wortel1(x)
% y = wortel1(x) trekt ruwweg de wortel uit x                    (maar nu even niet)
som = 0;
term = 1;
teller = 0;
while term < x
    teller = teller + 1;
    som = som + term;
    term = term + 2;
end
wortel1 = teller - 0.5;


Opdracht 6b

Copy en paste deze functies in twee aparte M-files (noem deze resp. wortel.m en
wortel1.m). Om het debuggen te kunnen illustreren, is er in één van de M-files een fout
gestopt. Run wortel.m met als argument bijvoorbeeld 16. Typ dus in het Command
Window: wortel(16). Je zult ontdekken dat het antwoord niet klopt.

De fout die zich in één van de M-files bevindt, kun je gaan opzoeken m.b.v. de debugger.
De eerste stap in het debug-proces is uiteraard het openen van een M-file. Mocht je nog
wijzigingen aanbrengen in de file, dan moet je deze altijd opslaan vóórdat je begint met
debuggen (niet opgeslagen files kun je niet debuggen). Ter voorbereiding op het
debuggen, plaats je vervolgens zogenaamde “breakpoints”. Deze plaats je op plekken
waarop de uitvoering van de M-file onderbroken moet worden, zodat je de waarden van
de variabelen (waarvan je vermoedt dat ze problemen geven) kunt bekijken. Breakpoints
kun je plaatsen door op de breakpoint-kolom (de kolom rechts van het regelnummer) te
klikken.

Opdracht 6c

Plaats in de file wortel.m een breakpoint.

Een breakpoint wordt als het goed is weergegeven met een rood bolletje. Verschijnt er
een grijs bolletje, dan is het breakpoint niet geldig. Dit kan twee oorzaken hebben: de file
is niet opgeslagen of de file bevat nog een syntax fout. Uiteraard moeten beide problemen
verholpen zijn voordat met debuggen begonnen kan worden. Tijdens het debuggen kun je
overigens niet alleen de M-file debuggen waarin je werkt, maar ook de M-files waarvan
gebruik wordt gemaakt in de M-file die je debugt. In dit geval is het daarom dus handig
de file wortel.m te gebruiken als file die je gaat debuggen, aangezien daarin de file
wortel1.m ook gebruikt wordt en je die dus tegelijkertijd kunt debuggen.




MatLab                                         29                                         2-3-12
N.B. Onder het “Breakpoints-menu” vind je ook een aantal opties die beginnen met “stop
     if...”. In plaats van het plaatsen van breakpoints, kun je ook deze opties aanvinken.
     Wanneer je een M-file runt (onder het “Debug-menu”), zal MatLab automatisch
     stoppen op de plek waar het programma een run time fout tegen komt en hier een
     melding van maken. De verschillende opties geven aan voor welke fouten MatLab
     moet stoppen; voor een omschrijving van deze opties kun je de helpfunctie
     raadplegen. (Uiteraard spoort MatLab geen fouten op die het gevolg zijn van het
     gebruik van bijvoorbeeld een niet-bedoelde, maar wel correcte functie).

Als de breakpoints geplaatst zijn, kan het debuggen beginnen. Je start het debuggen door
in je Command Window de opdracht te geven de M-file uit te voeren.

Opdracht 6d

Run de M-file nogmaals met als argument 16. Typ dus in het Command Window:
wortel(16).

De M-file wordt nu gerund tot de regel waar het eerste breakpoint geplaatst is. De
debugger biedt nu een aantal hulpmiddelen om de fout op te sporen. Ten eerste kun je
tijdens het debuggen de waarden van de variabelen opvragen. Dit doe je eenvoudig door
in de M-file met je cursor op een variabelenaam te gaan staan. Wanneer MatLab een
onverwachte waarde laat zien, weet je dat er op dat punt of eerder in de file een fout in de
file geslopen is. Het is van belang je te realiseren dat je de waarde te zien krijgt die de
variabele heeft op het moment dat de uitvoering van de M-file tijdelijk wordt
onderbroken.

Er zijn er een aantal mogelijkheden die allemaal te vinden zijn onder “debug”. Een korte
beschrijving van deze mogelijkheden:
     Step:       voer de regel waarop de cursor staat uit en stap daarna naar de volgende
                  regel.
     Step in: voer de regel waarop de cursor staat uit en stap in de M-file die
                  aangeroepen wordt.
     Step out: voer de aangeroepen M-file uit, stap uit de aangeroepen M-file, keer
                  terug naar de M-file die je aan het debuggen was en pauzeer op het punt
                  waarop je terugkomt.
     Continue: voer de commando’s tot het volgende breakpoint uit en pauzeer daar.

Opdracht 6e

Stap door de beide M-files met de boven beschreven opties en houd daarbij de volgorde
aan die MatLab zelf ook zou gebruiken bij het runnen van deze files. (Houd er daarbij
rekening mee dat sommige commando’s meerdere malen uitgevoerd moeten worden).
Bekijk in elke stap de waarden van de variabelen en spoor op waarom het met deze twee
files niet lukt om de juiste wortel uit te rekenen.




MatLab                                         30                                         2-3-12
Als je de fout hebt opgespoord, kun je deze gaan corrigeren. Dit doe je in een aantal
stappen. Ten eerste stop je het debug-proces, met de optie “Exit Debug Mode”, onder
“debug”. Breng de wijzigingen aan. Sla de M-file op en run de M-file opnieuw. Mocht
deze nog steeds fouten bevatten, doorloop het debug-proces dan nogmaals. Krijg je bij
het runnen geen onverwachte resultaten meer, dan kun je tevreden zijn!

INLEVEROPDRACHT 6f

Herstel de fout, sla de juiste versie van de M-files op en lever deze in als eindproduct van
de opdrachten 6a t/m 6f.

7.       Fourier-analyse
In dit hoofdstuk zullen we Fourier-analyses gaan uitvoeren. We beginnen zoals
gewoonlijk eenvoudig en werken langzaam naar wat gecompliceerdere toepassingen.
Zoals je weet is de Fourier-analyse een manier om een complex signaal te beschrijven als
de som van een aantal sinusvormige signalen. Die beschrijving geeft je een veel
duidelijker beeld van de structuur van het oorspronkelijke signaal. En vanuit die structuur
is het mogelijk gericht bepaalde (storende) frequenties weg te filteren. Dat filteren van
een signaal komt in het volgende hoofdstuk aan bod.

7.1 Een eenvoudig signaal opbouwen
Om een goed beeld te krijgen van wat een Fourier-analyse nu in concreto doet, gaan we
eerst zelf een eenvoudig signaal opbouwen, waarvan we precies weten welke frequenties
erin zitten. Vervolgens voeren we een Fourier-analyse uit om te controleren of deze het
signaal weer kan ontleden in de oorspronkelijke frequenties.

Opdracht 7.a:

Schrijf een programma dat vraagt (input) om twee frequenties en vervolgens twee
sinusvormige signalen met deze frequenties optelt en de som laat zien. Noem het
uiteindelijke signaal voor het gemak x en definieer de duur in een variabele met de naam
d.

De frequenties zullen verkeerd worden weergegeven als ze groter zijn dan de helft van de
sampling frequentie (zie de wet van Nyqist). Een frequentie, die hoger is dan de helft van
de samplingfrequentie wordt als het ware gespiegeld ten opzichte van deze halve
samplingfrequentie. Wanneer bijvoorbeeld de samplingfrequentie 200Hz is en het signaal
een frequentie 130Hz bevat, dan zie je die na Fourier-analyse in het frequentiespectrum
terug als een frequentie van 70 (30 onder de halve samplingfrequentie in plaats van 30
erboven). Houd daarmee rekening bij het uitproberen.

Verfraai de grafische weergave door er correcte titels aan toe te voegen. Raadpleeg
hiervoor één van de mogelijke hulpmethoden.




MatLab                                         31                                         2-3-12
Voer het programma uit en bekijk het resultaat. We zijn nu klaar voor een Fourier-
analyse op dit signaal.

Opdracht 7.b:

Voeg aan het programma van opdracht 7.a de volgende regels toe.
a= abs(fft(x));                    %   Fourier-analyse: correspondeert met
                                   %   freq. 0 t/m sf (sf = sampling frequentie)
a2 = a(1:d*sf/2);                  %   freq. 1 t/m (sf/2)
f=(1/d)*(0:length(a2)-1);          %   frequentie-as op maat
subplot(2,1,2), plot(f,a2);


Let er hierbij op dat de vector x staat voor het signaal dat je in opgave 7.a gemaakt hebt;
d staat voor de duur in seconden.

Vanwege de wet van Nyquist laten we alleen de linker helft van het frequentiespectrum
zien; die helft wordt geplaatst in a2. Het commando fft(x) geeft zijn resultaten alsof het
signaal maar één seconde duurt. Als het signaal d seconden duurt, zal een frequentie van
bijvoorbeeld 1000Hz dus worden gezien als d*1000Hz. Om in de plot alle frequenties tot
aan de Nyquist-frequentie (sf / 2) te kunnen zien moeten we dus de eerste d * sf / 2
waarden bekijken. Omdat we op de x-as de echte frequenties willen hebben corrigeren we
f door te delen door d. Als je alleen geinteresseerd bent in frequenties beneden een
bepaalde waarde, bijvoorbeeld 1500Hz, kun je natuurlijk a2 inkorten: a2 =
a(1:d*1500);

Wanneer je dit aangepaste programma draait zul je zien dat de frequenties die je in het
signaal hebt gestopt door de Fourier-analyse netjes worden teruggevonden. Het overzicht
van de frequenties en hun amplitudes heet het frequentiespectrum van het signaal.

We hebben nu twee beschrijvingen gezien van hetzelfde signaal. De eerste beschrijft de
amplitude als een functie van de tijd. De tweede, het resultaat van de Fourier-analyse,
beschrijft de amplitude als een functie van de frequentie. Beide beschrijvingen, die in het
tijdsdomein en die in het frequentiedomein, zijn bruikbaar. In de ene situatie is het
tijdsdomein handig, in de andere het frequentiedomein. De stap van tijdsdomein naar
frequentiedomein wordt gezet wanneer een Fourier-analyse wordt uitgevoerd op het
signaal in de tijd. De stap van frequentiedomein naar tijdsdomein wordt gezet wanneer
op basis van het frequentiespectrum het signaal wordt opgebouwd als een som van
sinusvormige signalen.

7.2      Een signaal binnenhalen
In plaats van een eenvoudig zelfgemaakt signaal heb je in de praktijk meestal te maken
met veel complexere signalen, waar je zelf niet of slechts ten dele de hand in hebt gehad.
We gaan in de volgende opdracht een Fourier-analyse uitvoeren op een binnengehaald
signaal (gelezen via het commando wavread).




MatLab                                        32                                          2-3-12
INLEVEROPDRACHT 7.c:

Schrijf een programma dat een signaal binnenhaalt en er vervolgens een Fourier-analyse
op uitvoert. Noem het binnengehaalde signaal y en bereken de duur als d. Gebruik voor je
Fourier-analyse de volgende opzet:
a= abs(fft(y));                       %   Fourier-analyse: correspondeert met
                                      %   frequenties 0 t/m sf
a2 = a(1:ceil(length(a)/2));          %   frequenties 0 t/m (sf/2)
f=(1/d)*(0:length(a2)-1);             %   frequentie-overzicht


De functie ceil(x) levert het kleinste hele getal, dat groter is dan x of eraan gelijk. Het
geeft dus een afronding naar boven. (engels: ceiling = plafond)

Plot zowel het signaal in het tijdsdomein als het signaal in het frequentiedomein.

8.       Filteren
Zoals we in het vorige hoofdstuk hebben gezien, is het uitvoeren van een Fourier-analyse
een krachtige manier om een op het oog zeer complex signaal te beschrijven als de som
van een aantal eenvoudige sinusvormige signalen. Een dergelijke analyse maakt het
mogelijk bepaalde (storende) frequenties uit het signaal te verwijderen. Om dat te laten
zien, beginnen we weer met een eenvoudig signaal en werken we langzaam naar meer
realistische en dus ook complexere signalen toe.

8.1      Filteren van een eenvoudig signaal
Bij het filteren wordt een signaal zodanig veranderd, dat bepaalde frequenties onderdrukt
worden. Dat gebeurt door elke binnenkomende waarde te vervangen door een gewogen
som van de laatste zoveel voorafgaande binnenkomende waarden en (bij sommige filters)
de laatste zoveel aangepaste waarden. (Omdat je voor het filteren van elk punt een aantal
voorgaande punten nodig hebt, zul je zien dat in het gefilterde signaal de eerste zoveel
punten wegvallen. De lengte van dat weggevallen gedeelte is gelijk aan de halve lengte
van het filter.) In het volgende voorbeeld gebruiken we een FIR-filter (FIR = Finite
Impulse Response). Een FIR-filter gebruikt alleen de voorafgaande binnenkomende
waarden.

De kunst is natuurlijk om de goede gewichten te vinden. Dat kan (o.a.) via een
commando als:

b = fir2(n,f,m);


De parameters n, f en m hebben daarbij de volgende betekenis:

f is een vector met uit te filteren frequenties. Die frequenties moeten echter zijn gedeeld
door de Nyquist-frequentie (dat is dus de halve sampling frequentie). De waarden moeten
bovendien in oplopende volgorde worden opgegeven: van de laagste frequentie naar de
hoogste.


MatLab                                         33                                        2-3-12
m  is een vector met de corresponderende amplitudes zoals men die uit het filter wil
krijgen: een nul betekent complete onderdrukking, een 1 betekent volledig doorlaten.

n is de lengte van het filter, d.w.z. hoeveel voorgaande waarden meedoen in de gewogen
som. Hoe groter n, hoe beter het filter zal werken, maar ook hoe meer tijd het filteren zal
kosten. De lengte van het filter moet altijd even zijn.

Zoals je ziet, is het filter onafhankelijk van het te filteren signaal. Je kunt het één keer
uitrekenen en vervolgens op allerlei signalen toepassen. Zolang ze dezelfde sampling
frequentie hebben, zal het filter altijd dezelfde frequenties wegfilteren of afzwakken.

Ter verduidelijking een voorbeeld:
We hebben een signaal gemeten met sampling frequentie 500. De Nyquist-frequentie is
dus 250. Het gewenste filter heeft de volgende karakteristiek:

Frequentie     f               m
20             0.08            1.0
40             0.16            0.0
60             0.24            1.0

We willen dus de frequentie 40Hz onderdrukken. De vraag is natuurlijk wat er moet
gebeuren met frequenties dicht bij die 40. Het is daarom verstandig om nog een paar
waarden toe te voegen rond de 40 Hz. Neem ook de nul en de Nyquist-frequentie op.
Bijvoorbeeld:

Frequentie     f               m
0              0.080           1.0
20             0.080           1.0
38             0.152           1.0
40             0.160           0.0
42             0.168           1.0
60             0.240           1.0
250            1.000           1.0

Als het oorspronkelijke signaal in de vector x is opgeslagen, kunnen we het nu als volgt
filteren:
sf = 500;                               %   sampling frequentie;
ny = sf/2;                              %   Nyquist-frequentie
n = 100;                                %   lengte van het filter (moet even zijn)
freq = [0,20,38,40,42,60,250];          %   de frequenties van 0 t/m ny
f = freq / ny;                          %   frequenties genormeerd
m = [1, 1, 1, 0, 1, 1, 1];              %   gewenste relatieve sterkte
b = fir2(n,f,m);                        %   het filter
y = fftfilt(b,x);                       %   het gefilterde signaal


Het gefilterde signaal staat nu in de vector y.




MatLab                                            34                                      2-3-12
Opdracht 8.a:

Maak eerst een programma dat om twee frequenties vraagt en een signaal berekent dat is
opgebouwd uit twee sinusvormigen met de opgegeven frequenties. Ook moet het
programma op dit signaal een Fourier-analyse uitvoeren en de resultaten overzichtelijk
plotten.

Voeg nu aan dit programma een aantal opdrachten toe, die ervoor zorgen, dat een
frequentie naar keuze kan worden uitgefilterd. Als de uit te filteren frequentie in de
variabele Weg staat, kun je het filteren als volgt voorbereiden:
r    =   Weg/(0.5*sf); % Frequentie gedeeld door        halve sampling frequentie
n    =   100;                                           % Filter-lengte
f    =   [0.0, 0.9*r, r, min(1.1*r,0.999), 1.0];        % frequenties genormeerd
m    =   [1, 0.1, 0, 0.1, 1];                           % gewenste doorlaatsterkte


Als we de vierde waarde in f zouden definieren als 1.1*r zou zij groter kunnen zijn dan 1
en dat is niet toegestaan (het zou een frequentie voorstellen boven de Nyquist-frequentie).
Daarom is de waarde gedefinieerd als min(1.1*r,0.999), het minimum van 1.1*r en
0.999.

Vervolmaak nu het programma door een Fourier-analyse op het gefilterde signaal uit te
voeren en alle vier de resultaten op een overzichtelijke manier te plotten.

8.2       Filteren van een complex signaal
Nu een wat realistischer voorbeeld. We gaan weer gebruik maken van een wav- bestand.
Daarbij gaan we niet één frequentie uitfilteren, maar alle hoge frequenties. Het is
hiervoor raadzaam om nog eens te kijken naar het filter in het eenvoudige voorbeeld.

INLEVEROPDRACHT 8.b:

Schrijf een programma dat de volgende taken uitvoert:

1.    Lees een wav-bestand.
2.    Laat het gelezen signaal horen.
3.    Voer een Fourier-analyse uit op het gelezen signaal.
4.    Laat het signaal zien in het tijdsdomein en in het frequentiedomein.
5.    Laat vragen om een bovengrens voor de frequenties in het signaal.
6.    Filter alle frequenties boven de bovengrens uit (een kopie van) het signaal.
7.    Laat het gefilterde signaal horen.
8.    Voer een Fourier-analyse uit op het gefilterde signaal.
9.    Laat via een viervoudige plot het oorspronkelijke signaal en het gefilterde signaal
      zien in het tijdsdomein en in het frequentiedomein.




MatLab                                        35                                         2-3-12
9        Plotten
MatLab bevat uitgebreide mogelijkheden voor het maken van grafieken en plots.

9.1      Het commando en de functie Plot
Plot(Y)
Plot geeft een 2-dimensionale weergave van een reeks getallenparen. In zijn eenvoudigste
vorm geef je als enig argument één vector mee: Plot([2 4 6 8 10]) geeft een rechte te
zien van het punt (1,2) linksonder naar het punt (5,10) rechtsboven. Op de x-as komen de
getallen 1, 2, ... op de y-as de elementen van de vector.
Plot(X,Y)
Met deze vorm bevat X de punten op de X-as en Y die op de Y-as. X en Y mogen zowel
rij- als kolomvectoren zijn.
Plot(X,Y,X2,Y2)
Op deze manier worden er twee lijnen getekend in verschillende kleuren: de ene wordt
bepaald door X en Y, de andere door X2 en Y2.
Hold on
Zorgt ervoor, dat de huidige plot wordt vastgehouden, d.w.z. dat een volgend plot-
commando zijn resultaat daaraan toevoegt.
Hold off
Zorgt ervoor, dat een volgend plot-commando het laatste plaatje wist en alle as-
instellingen weer op hun default-waarde zet
Plot(X,Y,':ob',X2,Y2)
Op deze manier worden er twee lijnen getekend in verschillende kleuren: de ene wordt
bepaald door X en Y, de andere door X2 en Y2. Het derde argument geeft een paar
stijlkenmerken voor de eerste plot: de dubbele punt voor het lijntype (puntjes), de letter o
voor de data-punten en de letter b voor de kleur (blauw). De specificatiestring kan de
volgende elementen bevatten:
-       doorlopende lijn
--      onderbroken lijn
:       lijn bestaat uit puntjes
_.      lijnen en puntjes om en om
+       data-punten zijn plusjes
o       data-punten zijn rondjes
*       data-punten zijn sterretjes
.       data-punten zijn punten
x       data-punten zijn kruisjes
s       data-punten zijn vierkantjes (squares)
d       data-punten zijn ruitjes (diamonds)
^       data-punten zijn driehoekjes met de punt omhoog
v       data-punten zijn driehoekjes met de punt omlaag


MatLab                                         36                                         2-3-12
>        data-punten zijn driehoekjes met de punt naar rechts
<        data-punten zijn driehoekjes met de punt naar links
p        data-punten zijn vijfpuntige sterren (pentagrams)
h        data-punten zijn zespuntige sterren (hexagrams)
r        kleur is rood
g        kleur is groen
b        kleur is blauw
c        kleur is cyaam
m        kleur is magenta
y        kleur is geel (yellow)
k        kleur is zwart (black; b is voor blauw)
w        kleur is white

Als je geen specificatie-string geeft, geeft Plot de verschillende lijnen kleuren in de
volgorde die is aangegeven in ColorOrder. Dat is een zoveel bij 3 matrix met RGB-
waarden tussen nul en één: ze geven voor iedere kleur de sterkte van de rood-, groen- en
blauw-component.      De      inhoud    van    ColorOrder     kun     je   zien     met:
get(gca,'ColorOrder').          Je     kunt     de     waarden      veranderen      via:
set(0,'DefaultAxesColorOrder',kleuren), waarbij kleuren een zoveel bij 3 matrix
is.

De defaultwaarden zijn:
0.00     0.00   1.00   (blauw)
0.00     0.50   0.00   (groen)
1.00     0.00   0.00   (rood )
0.00     0.75   0.75   (magenta)
0.75     0.00   0.75   (cyaan)
0.75     0.75   0.00   (geel)
0.25     0.25   0.25   (grijs)

Als alle kleuren aan de beurt zijn geweest, wordt de reeks herhaald met de volgende
waarde van LineStyleOrder enzovoort. MatLab accepteert maximaal 4 lijnstijlen, die je
kunt opgeven in één string met rechte strepen tussen de stijlen, bijvoorbeeld:
set(0,'DefaultAxesLineStyleOrder','-|--|:'). Je kunt beide specificaties
combineren in één set-statement. Probeer eens het volgende:
set(0,'DefaultAxesColorOrder',[1 0 0;0 1 0;0 0 1],               ← hier geen Enter
      'DefaultAxesLineStyleOrder','-|--|:')
t = 0:pi/20:2*pi;
a = ones(length(t),9);
for i = 1:9
    a(:,i) = sin(t-i/5)';
end
plot(t,a)




MatLab                                         37                                     2-3-12
9.2      Subplot: een matrix van plots
Het commando subplot stelt je in staat een matrix te maken met in elke cel een plot. Dit
commando verdeelt namelijk het figuurscherm in een aantal rechthoeken. Achter het
subplot-commando specificeer je daarom de grootte van de tabel (uitgedrukt in een
parameter voor het aantal rijen, gevolgd door een parameter voor het aantal kolommen)
en de plaats waar het eerstvolgende plaatje in deze matrix moet staan (uitgedrukt in een
“plaats-parameter”, waarbij de cellen van links naar rechts en vervolgens van boven naar
beneden geteld worden, de cel linksboven heeft dus plaats 1). Een voorbeeld:
 subplot(2,2,1); plot(t,x)
 subplot(2,2,2); plot(t,y)
 subplot(2,2,3); plot(t,z)


In dit voorbeeld wordt het plaatje verdeeld in vier rechthoeken. Linksboven komt de plot
van t tegen x, rechtsboven die van t tegen y en links onder die van t en z. De rechthoek
rechtsonder blijft leeg.

Als je een plot op een andere plaats wil hebben kun je dat zo doen:
subplot('position',[links onder breedte hoogte])
Als de zo gekozen plot over een al bestaande subplot heen valt, wordt die laatste
verwijderd. De argumenten links, onder, breedte en hoogte geven de afstand van de
linkerrand, de afstand van de onderrand, de breedte en de hoogte aan, dit alles als
proporties van de totale figuur, dus met getallen tussen 0 en 1.

9.3      Plot: nog meer vormgeving
Als een plot gemaakt is, kun je er nog een heleboel aan veranderen met behulp van de
plot-editor: klik op het pijltje links omhoog om de editor te activeren. Dubbelklik daarna
bijvoorbeeld op een van de assen. Alle opties voor het bewerken van de assen verschijnen
dan in een venster.

Maar je kunt ook vanuit het command window (of vanuit een m-file) werken. In het plot-
commando kunnen, na de eerdere argumenten nog een hele reeks van kenmerken worden
opgegeven. Dat gebeurt ook weer via paren van parameters. De eerste is een string, die
aangeeft om welk kenmerk het gaat, de tweede is de gewenste waarde van het kenmerk:

'MarkerEdgeColor',kleur       : kleur van de rand van de merktekens (data-punten)
'MarkerFaceColor',kleur       : vulkleur voor de merktekens
'LineWidth',dikte                                              1
                              : lijndikte in punten (1 punt is 72 inch)
'MarkerSize',grootte          : grootte van de merktekens in punten




MatLab                                        38                                        2-3-12
De kleuren kunnen op drie manieren worden aangegeven: via hun RGB-waarden (de
sterkte van rood, groen en blauw, tussen 0 en 1), via een letter of via een woord:

RGB-waarde      Letter   Woord
[1 1 0]         y        yellow
[1 0 1]         m        magenta
[0 1 1]         c        cyan
[1 0 0]         r        red
[0 1 0]         g        green
[0 0 1]         b        blue
[1 1 1]         w        white
[0 0 0]         k        black

De plaats van de maatstreepjes en hun labels kun je aanpassen door na het plot-
commando set-commando te gebruiken. De volgende twee commando's
set(gca,'XTick',0:5:20)
set(gca,'XTickLabel',{'0','|','10','|','20'})
zetten op de x-as streepjes bij de waarden 0, 5, 10,15,20 en labelen ze als aangegeven.

Je kunt aan een plot nog meer opmaak toevoegen, zoals titel-regels, as-labels en legenda.
De volgende commando's doen dat:
xlabel 'aantal spruitjes in 2002'
ylabel 'aantal spruitjes in 2003'
title 'spruitjesproductie'
text(-pi/4,sin(-pi/4),'\leftarrow sin(-\pi\div4)')


9.4      Plot: de assen
Hoe de x- en de y-as van een plot eruit moeten zien valt te regelen met de commando's
axis, xlim en ylim. Geef deze commando's nadat de plot gemaakt is. De minimale en
maximale waarden zet je alsvolgt:
axis([minx,maxx,miny,maxy])           % grenswaarden voor x en y-as
xlim([minx,maxx])                     % grenswaarden voor x-as
ylim([miny,maxy])                     % grenswaarden voor y-as


Andere mogelijkheden zijn:
axis   auto     %   gebruik   de default instellingen     voor beide assen
xlim   auto     %   gebruik   de default instellingen     voor de x-as
ylim   auto     %   gebruik   de default instellingen     voor de y-as
axis   manual   %   houd de   geldende schaal vast
xlim   manual   %   houd de   geldende schaal voor de     x-as vast




MatLab                                        39                                          2-3-12
ylim manual % houd de geldende schaal voor de y-as vast
axis tight % Zet de grenzen gelijk aan die van de data
axis fill   % Schaal zo, dat het hele plaatje wordt gebruikt
            % (alleen als PlotBoxAspectRatioMode of DataAspectRatioMode op manual staan)
axis ij     % zet de oorsprong links boven
axis xy     % zet de oorsprong links onder
axis equal % gelijke intervallen op beide assen even groot weergeven
axis image % idem, maar leg de randen strak om de data
axis square % beide assen even groot
axis normal % maak 'axis square' en 'axis equal' ongedaan
axis vis3d % bevriest aspect ratio om rotatie mogelijk te maken
axis off    % geen labels, maatstreepjes of achtergrond
axis on     % maak axis off ongedaan


Onderstaande functie circle gebruikt het commando axis square om te zorgen, dat er
echt een cirkel zichtbaar wordt en niet een of andere ellips.

function circle = circle(x,y,r);
% circle(x,y,r) plots a circle with center (x,y) and radius r
t=[0:2:360];
degtorad=pi/180; % van graden naar radialen
xc=x+r.*cos(t.*degtorad);
yc=y+r.*sin(t.*degtorad);
plot(xc,yc,'g');
axis square


9.5      Plot: de legenda
De legenda (het leesvoorschrift) bij een plot kun je regelen met het commando legend.
Je moet dat commando geven nadat de plot gemaakt is. De volgende opties zijn
beschikbaar:
legend('string1','string2',...)           %   de omschrijving per lijn
legend(matrix)                            %   idem: teksten in rijen van matrix
legend('off')                             %   geen legenda
legend('hide')                            %   maak legenda onzichtbaar
legend('show')                            %   maak legenda zichtbaar
legend('boxoff')                          %   geen rechthoek om de legenda
legend('boxon')                           %   teken rechthoek om de legenda

legend(...,-1)                            %   legenda rechts buiten de assen
legend(...,0)                             %   legenda binnen de assen
                                          %   (zoek meest lege plek)
legend(...,1)                             %   legenda rechts boven binnen de assen
                                          %   (dit is de default)
legend(...,2)                             %   legenda links boven binnen de assen
legend(...,3)                             %   legenda links onder binnen de assen
legend(...,4)                             %   legenda rechts onder binnen de assen


Als een plot in beeld is, kun je de legenda verslepen. Als je op een tekst dubbelklikt kun
je hem aanpassen.




MatLab                                         40                                          2-3-12
INLEVEROPDRACHT 9

Genereer een normaalverdeling met de onderstaande M-file:

y = randn(1,1000000);
freq = zeros(1,100);
for i = 1:1000000;
    j = round((y(i)+5) * 10);
    if j > 0 & j < 101         % laat uitschieters buiten beschouwing
        freq(j) = freq(j) + 1;
    end
end

Leg stap voor stap uit wat er gebeurt in deze M-file

Plot vervolgens de verdeling (gebruik freq op de y-as) in een magenta stippellijn, en
voorzie de y-as en de totale figuur van een titel. Plot dan in dezelfde figuur (maar in een
grijze ononderbroken lijn) de verdeling die je krijgt als je 500000 getallen random laat
selecteren uit een normaalverdeling.

10         Afgeleiden
Wanneer je onderzoek doet op het terrein van motoriek (schrijven, oogbewegingen, ...)
zul je al gauw te maken krijgen met begrippen als snelheid en versnelling. Zo is er
bijvoorbeeld een direct verband tussen versnelling en het gebruik van kracht. Nu is
snelheid de afgeleide van plaats en versnelling de afgeleide van snelheid. Daarom een
kleine uitleg over het berekenen van afgeleiden uit gesamplede gegevens.

De afgeleide van een functie geeft aan hoeveel de functiewaarde verandert als haar
argument verandert. Een eenvoudig voorbeeld is positie als functie van tijd. Een auto rijdt
van Nijmegen naar Arnhem. Om drie uur is hij op de Waalbrug, om twee over drie is hij
6 km. verder. Als de functie f(t) de afstand van de Waalbrug (in meters) aangeeft op
tijdstip t (in seconden na 3 uur), dan kunnen we dus zeggen: f(0) = 0 en f(120) = 6000.
De gemiddelde snelheid gedurende die ene minuut was dus 6000/2 = 3000 meter per
seconde (180 km. per uur!). Verkeersagenten zijn goed in dit soort wiskunde. De snelheid
op een bepaald moment bepalen we door de duur tussen de tijdstippen heel klein te
nemen : dus in feite door de gemiddelde afstand te berekenen gedurende een heel kleine
periode. Deze snelheid v (als functie van het tijdstip t) is de afgeleide van de positie:

          f (t  h)  f (t  h)
vt                           voor heel kleine h
                  2.h

De afgeleide van de snelheid (dus de verandering van de snelheid in de loop van de tijd)
noemen we de versnelling2. De berekening zou als volgt kunnen:



2
    En uit versnelling en massa valt weer de uitgeoefende kracht af te leiden


MatLab                                                    41                             2-3-12
 xydata = load('xydata.txt')
 sf = 100;                        % sampling frequentie
 snelheidxy = afgeleid(xydata,sf) % snelheid in x- en in y-richting
 % snelheid overall (Pythagoras):
   snelheid = [snelheidxy(:,1).^2 + snelheidxy(:,2).^2].^0.5
   versnelling = afgeleid(snelheid,sf) % versnelling
 % snelheid x, snelheid y, snelheid overall en versnelling:
   snelheid = [snelheidxy,snelheid,versnelling]


De onderstaande functie gebruikt een iets andere formule voor de berekening van de
afgeleide. Zij berekent de afgeleide op tijdstip (lees: samplenummer) i uit de twee
voorgaande en de twee volgende waarden, maar de dichtsbijzijnde hebben een zwaarder
gewicht:
function [afgeleide] = afgeleid(y,sf)
% x = afgeleid(y,sf): derivatives for the columns of y;
% sf = sampling frequency (in Hz);
% y = the matrix: rows represent subsequent samples
% (met dank aan Ron Jacobs en Ruud Meulenbroek)
[m,n] = size(y); % m rijen en n kolommen
for i=3:m-2,
  afgeleide(i,:) = (-y(i-2,:)*1.5-y(i-1,:)*4+y(i+1,:)*4+y(i+2,:)*1.5)*(sf/14);
end
afgeleide(1,:) = afgeleide(3:m); afgeleide(2,:) = afgeleide(3:m);
afgeleide(m,:) = afgeleide(m-2:m); afgeleide(m-1,:) = afgeleide(m-2:m);




MatLab                                    42                                      2-3-12

				
DOCUMENT INFO
Shared By:
Categories:
Tags:
Stats:
views:49
posted:3/2/2012
language:Dutch
pages:42