opgaver:Uge4

I denne opgave skal I, lige som i sidste uge, fabrikere jeres egne data - denne gang en Gaussisk linjeprofil med en bestemt bredde og placering, med Poisson statistik på tælletallene.

Opgave 1 - Spektrometer

Vi forestiller os et spektrometer, der tæller fotoner som funktion af bølgelængde; én bølgelængde ad gangen. Vi skal måle en enkelt spektrallinje. Linjens amplitude er 100 tællinger per sekund, linjens centrum 5000 Å, og den har bredde $\sigma$ 50 Å. Der er en baggrund på 10 tællinger per sekund. Usikkerheden på tællinger er givet som $\sqrt{N}$ tællestatistik, og antages gaussisk.

Spørgsmål 1 - Simulation af data

Simulér en måleserie fra 4500 Å til 5500 Å med step af 10 Å, ét sekund per punkt. Brug MATLAB fit funktionen til at fitte en Gaussisk linjeprofil til dataen. Vis de simulerede data og det bedste fit i samme plot.

Hint 

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.

Løsning 

Først defineres parametrene for den Gaussiske linjeprofil-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 linjeprofil, sammen med et Gaussisk fit.

og med disse kan man udregne fit-linjen

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-linjen,

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.

Spørgsmål 2 - Flere måleserier

Kør 100 måleserier og lav histogrammer over fitværdierne. Brug histogrammerne til at bestemme præcisionen af eksperimentet.

Hint 

Du skal bare bruge den samme kode som du brugte til spørgsmål 1, men nu skal du "pakke" den ind i et loop, hvor du yderligere gemmer dine fit-parametre i arrays, så de kan plottes bagefter.

Løsning 

For at gøre dette, skal man lave et loop over adskillige data-simulationer og fits. Som det første skal man initialisere arrays der skal gemme fit-parametrene:

N = 100;

a_fit = zeros(1,N);

b_fit = zeros(1,N);

x0_fit = zeros(1,N);

sig_fit = zeros(1,N);

Herefter kan man starte sit loop lige før man udregner sine eksperiment-data, altså

for j = 1:N

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

...

Fit parametrene gemmes, og så kan loopet stoppes igen,

...

a_fit(j) = f1.a/t;

x0_fit(j) = f1.b;

sig_fit(j) = f1.c;

b_fit(j) = f1.d/t;

end

Herefter kan man lave histogrammer over parametrene ved hjælp af MATLAB's hist funktion, og se på formen af dem. F.eks. kan man bestemme bredden af distributionerne, for at se hvor nemt det var for fittet at finde den rigtige værdi.

Spørgsmål 3 - Planlægning af måling

(SVÆRT) Hvis der er 100 sekunder til rådighed, hvordan skal målingen så planlægges for at bestemme linjens midtpunkt bedst? Og hvad med linjens bredde?

Hint 

Tænk over hvordan du ville gøre hvis du faktisk sad ved et eksperiment og intet vidste på forhånd. Ville du bare gå i gang fra en ende af og måle punkt efter punkt, eller ville du hoppe rundt og håbe på at du var heldig? Tænk over forskellige fordele og ulemper.

Når du har tænkt over de forskellige måder du kunne finde på at måle, så prøv at lave målingerne ved hjælp af MATLAB, og kør et loop så du kan sammenligne distributionerne for de fundne fitte-parametre.

Løsning 

Som en første iteration kan man se hvad der sker med statistikken hvis man måler på en ligeligt opdelt linje, med samme tælletid for hvert punkt. Start med at omdefinere din værdi af t til at blive udregnet med

t_meas = 100;

t = t_meas / length(x_exp);

hvor t_meas er den totale tid der skal måles i. Nu kan du køre dit loop, og se på distributionen af dine fitte-parametre. Prøv evt. at ændre på step-størrelsen af dine målinger, så du får færre eller flere datapunkter, mens du stadig tæller den samme tid totalt. Ændrer dette noget ved distributionerne?

Du kan også prøve at lave en dataserie hvor du først måler meget kort tid på mange ligeligt fordelte punkter, og så bagefter måler i lidt længere tid nogle af de steder hvor det ser ud til at der er en peak.

Hvis du gerne vil kode lidt mere avanceret, kan du jo forsøge at kode en søger-mekanisme, som altid vil forsøge at finde det højeste punkt af en peak.

Der er rigtig mange muligheder for at måle på det givne signal, og mange af dem er lige gode. Det vigtige med denne opgave er at I tænker over hvordan I ville optimere den tid I er givet til et sådant eksperiment.

Spørgsmål 4 - Stor baggrund

(SVÆRT) Lad baggrunden være 100 tællinger per sekund og amplituden være 10 tællinger per sekund. Hvor længe skal der tælles for at det kan afgøres, om der er en peak (amplitude signifikant større end nul)?

Løsning 

For at udregne svaret, løses det hvornår toppunktet for peaken ligger $n$ gange sin usikkerhed plus $n$ gange usikkerheden på baggrunden, væk fra baggrunden. Dette løses med

$\begin{align} (a+b)t - n\sqrt{(a+b)t} &= bt + n \sqrt{bt} \\ at - n\sqrt{(a+b)t} &= n \sqrt{bt} \\ at &= n (\sqrt{bt} + \sqrt{(a+b)t}) \\ a^2 t^2 &= n^2 (bt + (a+b)t + 2 \sqrt{b(a+b)}t) \\ a^2 t &= n^2 (2b + a + 2\sqrt{b(a+b)}) \\ t &= \dfrac{n^2 (2b + a + 2\sqrt{b(a+b)})}{a^2} \end{align}$

så hvis f.eks. $a=10$, $b=100$ og $n=1$, så er $t=4.2\,\text{s}$ for et enkelt målepunkt i toppen af peaken. Hvis man vil have bedre sikkerhed, så f.eks. $n=2$ (også kaldt et signifikant signal), så er $t=16.8\,\text{s}$.

Samlede løsninger

• Et samlet dokument med MATLAB kode til at løse alle opgaverne ovenfor kan hentes her: opgaver_uge4_spektrometer.m