opgaver:Uge4

Ligesom i sidste uge skal vi simulere nogle støjede data ud fra en kendt model. Start med at lave dine model-data, og tilføj så støj ved at lægge et tal til hvert af punkternes tælletal-værdi. Støjen skal være et tal der er Gaussisk fordelt omkring nul, med en bredde der er usikkerheden på punktet (og usikkerheden på punktet er jo givet ved tællestatistik!). Husk at tælletal er positive heltal, så du er nødt til at runde støjen op eller ned til heltal.

Det kan være en god idé at læse lidt i MATLAB hjælpen om funktionerne randn, round og abs.

Det er vigtigt at holde styr på forskellen mellem bredden af spektrogram-peaken som du laver data over, og bredden på den støj du lægger til dine simulerede datapunkter.

Først defineres parametrene for den Gaussiske linieprofil-model. Her defineres signalets amplitude a, baggrundens amplitude d, centrum af peaken x0 og bredden af peaken sig

a = 100;

b = 10;

x0 = 5000;

sig = 50;

videre skal akserne for eksperimentet bestemmes; start og slut for hvor der skal måles, x_min og x_max, samt step-størrelsen x_step.

x_min = 4500;

x_max = 5500;

x_step = 10;

Med disse definitioner er det nu muligt automatisk at lave sin x-akse med

x_exp = x_min : x_step : x_max;

og til sidst kan man også finde hvor lang tid der måles per datapunkt,

t = 1;

Nu kan man udregne sit model-datasæt ved at indsætte alle sine informationer i en Gaussisk funktion,

y_model = t * a * exp( - (( x_exp - x0 ).^2) ./ (2*sig^2) ) + t * b;

Nu er det tid til at lægge støj på punkterne, så man får et simuleret datasæt. Først og fremmest findes en vektor der indeholder usikkerheden på hvert punkt af sin model. Modellen er for tælletal, så denne usikkerhed er bare $\sqrt{N}$, og

err_model = sqrt(y_model);

Det nye datasæt kan så laves ved at lægge en Gaussisk støj med bredden givet element-vist i err_model, runde af til nærmeste heltal, og til sidst sørge for at tallet er positivt:

y_exp = abs(round( y_model + err_model .* randn(size(err_model)) ));

Usikkerheden på de simulerede datapunkter er dog ikke den usikkerhed man fandt for modellen. Hvis man målte disse punkter i den virkelige verden ville man i stedet bruge

err_exp = sqrt(y_exp);

Man kan nu lave et fit til de simulerede data, ved hjælp af MATLAB's indbyggede fit funktion. Denne gang skal man fitte til en funktion man selv skriver ind:

f1 = fit(x_exp',y_exp','a*exp(-((x-b).^2)/(2*c^2))+d',...

'Startpoint', [a*t x0 sig b*t]);

Her er Startpoint option'en de første gæt til de fire parametre, altså amplitude, centrum, bredde og baggrund. Her er de bare givet som modellens værdier - men i den virkelige verden kender man selvfølgelig ikke modellen perfekt, og man må bare give nogle cirka-gæt.

De fundne fit-parametre kan man få fra f1 objektet. Selve fit funktionen kalder dem for a, b, c og d - efter hvad der blev skrevet som fit-funktion. De numeriske værdier af de bedste parameter-værdier fås så ved f.eks. f1.a. De fire interessante fittede parametre er altså

a_fit = f1.a/t;

x0_fit = f1.b;

sig_fit = f1.c;

b_fit = f1.d/t;

De simulerede datapunkter for en Gaussisk linieprofil, sammen med et Gaussisk fit.

og med disse kan man udregne fit-linien

x_fit = x_min:0.05:x_max;

y_fit = t*a_fit(j)*exp(-((x_fit-x0_fit(j)).^2)./(2*sig_fit(j)^2)) + t*b_fit(j);

Til sidst kan man plotte de simulerede datapunkter sammen med fit-linien,

figure

errorbar(x_exp,y_exp,err_exp,'b.','MarkerFaceColor', [0.4 0.4 0.8],...

'MarkerEdgeColor',[0.4 0.4 0.8])

hold on

plot(x_fit,y_fit,'g-')

axis([x_min-50 x_max+50 0 (a+b)*1.2*t])

xlabel('Bølgelængde [Å]')

ylabel('Intensitet [tællinger/s]')

Figuren kan ses her til højre.