Adatfeldolgozás - egyenesillesztés

Egy életből ellesett példán keresztül azt próbálom megmutatni, hogy a Linux szerszámosláda segítségével milyen kényelmesen megoldható a mérési adatok kiértékelése és az eredmények megjelenítése - megközelítőleg zérus programfejlesztéssel.

Számos mintán és energián végeztünk mérési sorozatokat. Minden mérési sorozat megkezdése előtt a vákuumban in situ mintatisztítást végeztünk. Ezután a mintáról rugalmasan visszaszórt elektronok intenzitását mértük (mert valaki perverz módon éppen erre volt kíváncsi...) a mintatisztítástól eltelt idő függvényében. A mérést többször megismételve, szisztematikus intenzitáscsökkenést tapasztaltunk - ez a minta ultavákuumban történű elszennyeződésével(!) magyarázható. Ezen adatsorokat egy egyszerű (pl. lineáris) függvénnyel megillesztve, a mintatisztítás időpontjára (a legtisztább felületi állapotra) extrapolálhatunk.

1. lista Az adatokat ehhez hasonló eredményfájlok formájában kaptam:

#Ni 1B       
#1500 eV
#Measurement time: 1995. dec. 1.
#Notebook: ATOMKI ESA-31 No6. pages: 1-3.
Evaluated on 11 Dec, 95 at 10:51
from file rug389
time= 4.
Region # 1 between 
 1496.99 eV and  1501.47 eV   with markers at  1496.99 eV and  1501.47 eV
                                           raw        error 
            energy=                      1499.9+/-   6E-4      eV
              area=                      24024+/-   21.51 Counts/s*eV
           height =                      26682+/-  103.73 Counts/s
             FWHM =                     .79834+/-  .00428 eV

Evaluated on 11 Dec, 95 at 10:52
from file rug392
time= 23.
Region # 1 between 
 1496.99 eV and  1501.47 eV   with markers at  1496.99 eV and  1501.47 eV
                                           raw        error 
            energy=                      1499.9+/-   7E-4      eV
              area=                      23757+/-  21.402 Counts/s*eV
           height =                      26361+/-  103.11 Counts/s
             FWHM =                     .80078+/-  .00415 eV
...  És így tovább ...

Gyakorlatilag csak a from file, a time és az area adatok kellenek, így az alábbi AWK program segítségével a filenév.res fájlokból filenév.dat fájlokat készítünk. A #-jellel kezdődő megjegyzés sorokat változatlanul átvisszük a kimenő fájlba, s kiirjuk az összetartozó time és area értékeket (a spektrumazonosító csak emlékeztetőnek kell).

2. lista A res.awk program listája

NR == 1        { split(FILENAME,t,"."); fname=t[1] ".dat"
                 print "# " FILENAME > fname
               }
/^#/           {print $0 >> fname
                next
               }
/^from.file/   { spname=$NF }
/^time=/       { time= $NF  }
/area=/        { split($2,t,"+"); area=t[1]
                 print spname,time,area >> fname
               }

Néhány kézmozdulat (awk -f res.awk ni386425.res ) egy villanás, és máris készen áll az egyenesillesztés kiinduló adatsora. Mint látható, a kiindulási adatfájl nevét is "megörökítettem" az első megjegyzés sorában).

3. lista A NI386425.dat fájl listája, melyet a fenti program készített

# NI386425.RES
#Ni 1B
#1500 eV
#Measurement time: 1995. dec. 1.
#Notebook: ATOMKI ESA-31 No6. pages: 1-3.
rug389 4. 24024
rug392 23. 23757
rug395 42. 23633
rug398 63. 23410
rug401 81. 23295
rug404 103. 23181
rug407 123. 23270
rug413 123. 22563
rug414 236. 22318
rug417 260. 22307
rug420 281. 21923
rug423 302. 21811

Ezek után jöhet az egyenesillesztés! Kényelmi okokból egy meglevő FORTRAN eljárást vettem át, melyen csak apróbb módosításokat kellett végezni (indexes változóknál szögletes zárójelre cseréltem az eredetileg kerek zárójeleket, átírtam a programciklusok fejlécét - és nagy élvezettel kiirtottam a COMMON tömbökkel meg fájlnyitogatással kapcsolatos szerencsétlenkedéseket).

Az új "főprogram" változatlan formában átmásolja a megjegyzés sorokat, az x és y tömbökbe begyűjti az adatokat, s az adatsor végén (fájl vége vagy elválasztó üres sor észlelésekor) meghívja az illesztőrutint, mely egyúttal az eredmények kiírását is elvégzi.

A program "túlteljesíti" feladatát: a NF == 0 {if (n > 2) {lsq(n,x,y); n=0} } betoldás és az m számláló révén azt is kezelni tudja, ha egy bemenő adatfájlban több adatsor van - ezeket egyenként megilleszti, s az eredményfájlba is elválasztó sorokat tesz.

4. lista Az egyenesillesztő program (lsq.awk) listája

BEGIN   { n=0 ; m=0 }
NR == 1 { split(FILENAME,t,"."); fname=t[1] ".fit" }
/^#/    { print $0 >> fname; next }
NF == 3 { n++ ; x[n]=$2 ; y[n]=$3 }
NF == 0 { if (n > 2) { lsq(n,x,y); n=0 } }
END     { if (n > 2) { lsq(n,x,y) } }

function lsq(n,x,y) {
#
#     A legkisebb negyzetek modszerenek alkalmazasa
#     y=a+bx egyenletu egyenes illesztese a megadott pontokra
#
#     n       az adatpontok szama
#     x       a fuggetlen valtozo adatainak tombje
#     y       a fuggo valtozo adatainak tombje
#     a       a tengelymetszet erteke
#     b       iranytangens
#     ycalc   a szamitott ertekek tombje
#     stda    a tengelymetszet korrigalt empirikus szorasa
#     stdb    az iranytangens korrigalt empirikus szorasa
#
      sx=0.0
      sy=0.0
      sxx=0.0
      syy=0.0
      sxy=0.0
      sx2=0.0
      seps2=0.0
#
      for (i=1; i<=n; i++) {
          sx=sx+x[i]
          sy=sy+y[i]
          sxx=sxx+x[i]*x[i]
          sxy=sxy+x[i]*y[i]
          syy=syy+y[i]*y[i]
                           }
#
      a=(sx*sxy-sxx*sy)/(sx*sx-n*sxx)
      b=(sx*sy-n*sxy)/(sx*sx-n*sxx)
#
      for (i=1; i<=n; i++) {
          ycalc[i]=a+b*x[i]
          diff[i]=y[i]-ycalc[i]
          seps2=seps2+diff[i]*diff[i]
                           }
#
      koresz2=seps2/(n-2)
      szaml=n*sxx-sx*sx
      stda=sqrt(koresz2*sxx/szaml)
      stdb=sqrt(koresz2*n/szaml)

# A rajzoloprogram szamara elvalaszto sort irunk, 
# ha mar nem az elso adatsornal tartunk
 if (m != 0) { print "&" >> fname }

# ne felejtsuk el leptetni a szamlalot!
 m++
 
# Kiirjuk az eredmenyt (x,y,ucalc), de az y ertekeket ugy normaljuk,
# hogy a tengelymetszet eppen 100 % legyen
 print "# Intercept = " a " +/- " stda "   Slope = " b " +/- " stdb >> fname
 for (i=1; i<=n; i++) {
 printf("%8g %8g %8g \n", x[i],y[i]*100/a,ycalc[i]*100/a) >> fname
                       }
}

Lássuk a medvét... akarom mondani az eredményt! Az illesztés eredménye (a megjegyzés soroktol eltekintve) három oszlopból áll, mely az idő, a mért intenzitások és az illesztési paraméterekből számolt intenzitások adatait tartalmazza. Az intenzitásokat úgy normáljuk, hogy a tengelymetszet (a time = 0 értekhez tartozó számolt intenzitás) éppen 100 legyen, a relatív csökkenés mértéke így majd az ábráról könnyen leolvasható.

5. lista Az egyenesillesztés eredményének (NI386425.fit) listája

# NI386425.RES
#Ni 1B
#1500 eV
#Measurement time: 1995. dec. 1.
#Notebook: ATOMKI ESA-31 No6. pages: 1-3.
# Intercept = 23890.6 +/- 91.9209   Slope = -6.8219 +/- 0.540425
       4  100.559  99.8858 
      23  99.4409  99.3432 
      42  98.9219  98.8007 
      63  97.9885   98.201 
      81  97.5071  97.6871 
     103    97.03  97.0589 
     123  97.4025  96.4878 
     123  94.4432  96.4878 
     236  93.4176  93.2611 
     260  93.3716  92.5758 
     281  91.7643  91.9761 
     302  91.2955  91.3765 

Nézze meg az ember!

Az illesztés eredményét (akár egyszerre többet is) megjeleníthetjük az xvgr program segítségével:

xvgr -autoscale xy -nxy fájlnév1.fit fájlnév2.fit...

Ha tetszik az ábra, akár ki is nyomtathatjuk (vagy PostScript fájlt készíthetünk). Persze előtte nem árt egy kis szépítgetés. A kisérleti adatsort lehet pontokkal (filled circle) az illesztett egyenes pedig folytonos vonallal ábrázolni, stb.

Nem lehetne ezt a sok bíblődéssel járó beállítást automatizálni valahogyan? De igen! Csak egyszer kell megcsinálni a beállítást, s a paramétereket elmenthetjük pl. a graph.sty nevű fájlba. Az újabb eredmények megjelenitésénél ezt a -param graph.sty opcióval használhatjuk fel. Például:

xvgr -autoscale xy -nxy -param graph.sty fájlnév1.fit fájlnév2.fit...

Ha egyszer kényelmesnek születettünk...

Sok adatsorunk van, ugye nem fogjuk egyenként beírogatni a fájlok nevét? De nem ám! Készítsünk rá valami egyszerű shell scriptet!

"Hogyan is megvan esztet mondani matyarul? Á, becs fájl! Capisco!"


URL of this page: http://esca.atomki.hu
last modification:
Page maintained by Istvan Cserny   <cserny@atomki.hu>