Meß­da­ten­ana­ly­se mit LibreOffice

In der Ent­wick­lungs­pha­se des noch zu beschrei­ben­den Fre­quenz­zäh­lers habe ich eini­ge Meß­rei­hen erstellt. Eine davon war der Fre­quenz­gang eines 8 MHz Quar­zes auf dem CPU-Board, die ande­re der Fre­quenz­gang eines 100 MHz Quarz­os­zil­la­tors auf dem Fre­quenz­zäh­ler selbst. Bei­de Fre­quenz­gän­ge wur­den über die Tem­pe­ra­tur gemes­sen, die mit Käl­te­spray (Butan­gas, nicht Rau­chen, kein offe­nes Feu­er aber offe­ne Fen­ster!) bzw. einem Föhn geän­dert wurde.

In bei­den Fäl­len wur­de die Tem­pe­ra­tur mit dem Tem­pe­ra­tur­sen­sor auf dem CPU-Board gemes­sen, einem TMP275. Zwi­schen ‑40 °C und +125 °C hat er eine Tole­ranz von +/- 1 °C. Die tat­säch­li­che Tem­pe­ra­tur des Quar­zes und beson­ders des Quarz­os­zil­la­tors dürf­ten um eini­ge °C von der gemes­se­nen Tem­pe­ra­tur abwei­chen, weil sie etwas ent­fernt davon mon­tiert sind. Es soll hier aber nicht um die Prä­zi­si­on der Mes­sung gehen, son­dern ein­zig und allein um die Aus­wer­tung mit Libre­of­fice Calc. Der Fre­quenz­zäh­ler spei­chert die Meß­da­ten in einem CSV-File ab, das mit Libre­Of­fice direkt ein­ge­le­sen wer­den kann. Die bei­den Open­Of­fice Calc Datei­en ste­hen als Refe­renz und für eige­ne Expe­ri­men­te zum Down­load bereit:

In der ersten Tabel­le der jewei­li­gen Datei ste­hen die ori­gi­nal Meß­da­ten. Die wei­te­ren Tabel­len sind nach­fol­gend beschrieben.

Aus­wer­tung

Zur Aus­wer­tung emp­fiehlt es sich, zunächst eine Pivot-Tabel­le („Daten->Pivot-Tabelle->Einfügen) zu erstel­len. Als Zei­len­fel­der wählt man die Tem­pe­ra­tur und als Daten­fel­der die gemes­se­ne Fre­quenz Fcheck. Für die Daten­fel­der inter­es­siert uns aber nicht die Sum­me, son­dern der Mit­tel­wert. Das muß vor dem Erzeu­gen der Pivot-Tabel­le ein­ge­stellt wer­den. Nun sucht Calc alle in der Tabel­le vor­kom­men­den Tem­pe­ra­tu­ren, sor­tiert sie auf­stei­gend und ord­net jeder Tem­pe­ra­tur den Mit­tel­wert aller für die­se Tem­pe­ra­tur gefun­de­nen Meß­wer­te zu. So erhal­ten wir eine Tabel­le der gemit­tel­ten Meß­wer­te von etwa ‑30 °C bis +70 °C (Tabel­le: „Pivot-Tabel­le“).

Da wir einen 8 MHz bzw. einen 100 MHz Quarz­os­zil­la­tor gemes­sen haben, lie­gen alle Meß­wer­te ziem­lich nahe bei 8.000.000 Hz oder 100.000.000 Hz. Ohne ver­nünf­ti­ge Ska­lie­rung sieht man da nur eine hori­zon­ta­le Linie über der Tem­pe­ra­tur. Damit die Meß­wer­te sinn­voll dar­ge­stellt wer­den kön­nen, erstellt man in die­sem Fall am besten eine Spal­te mit der Dif­fe­renz der Meß­wer­te zum Nomi­nal­wert, denn alle Wer­te wei­chen nur sehr wenig davon ab. Die­ses Del­ta läßt sich dann gra­fisch über der Tem­pe­ra­tur darstellen.

Die lin­ken Dia­gram­me zei­gen die unver­än­der­ten Ori­gi­nal­da­ten. Man erkennt deut­lich den Tem­pe­ra­tur­gang, aber die vie­len Aus­rei­ßer sind erschreckend. Die Ori­gi­nal­da­ten ent­hal­ten die sehr schnel­le Abkühl­pha­se mit Käl­te­spray und die fast genau­so schnel­le Auf­heiz­pha­se mit dem Föhn. Dar­an schließt sich jeweils eine Pha­se lang­sa­men Auf­hei­zens bzw. lang­sa­men Abküh­lens auf Zim­mer­tem­pe­ra­tur an.

Die rech­ten Dia­gram­me zei­gen nur die lang­sa­men Pha­sen und die Aus­rei­ßer sind hier kom­plett ver­schwun­den (Tabel­len: „Aus­wahl“ und „Pivot-Tabel­le_1“). Die Kur­ven­form bleibt unver­än­dert. Offen­sicht­lich tun in bei­den Fäl­len die schnel­len Tem­pe­ra­tur­än­de­run­gen nicht gut, die Schwin­gung reißt ver­mut­lich kurz­zei­tig ab.

Erzeu­gen einer Ausgleichskurve

Alle wei­te­ren Aus­wer­tun­gen erfol­gen nun mit den bei­den „lang­sa­men Kur­ven“, denn sie spie­geln offen­sicht­lich die wah­ren Ver­hält­nis­se. Das Ziel ist nun, die jewei­li­ge Kur­ve so zu glät­ten, daß die geglät­te­te Kur­ve mög­lichst genau der gemes­se­nen Kur­ve ent­spricht, ohne deren Meß­ab­wei­chun­gen abzu­bil­den. Die Stan­dard­me­tho­de zum Erzeu­gen einer soge­nann­ten Aus­gleichs­kur­ve ist die Aus­gleichs­rech­nung. Für ein mathe­ma­ti­schen Lösungs­ver­fah­ren muß die­ser Aus­gleichs­rech­nung ein Modell zugrun­de lie­gen. Im all­ge­mei­nen Fall ist das ein Poly­nom der Form

y = a0 + a1 * x + a2 * x2 + … + an * xn

In der Pra­xis muß man die­se Rei­he natür­lich früh­zei­tig irgend­wo abbre­chen, denn sonst wird der Rechen­auf­wand immer grö­ßer. Es sei vor­ab ver­ra­ten, daß bei den gezeig­ten Kur­ven min­de­stens ein Poly­nom drit­ten Gra­des nötig ist, also n = 3. Das liegt an den zwei loka­len Mini­ma bzw. Maxi­ma (jeweils grob bei 0 °C und bei 60 °C), bei denen die erste Ablei­tung des Poly­noms null sein muß. Um zwei Null­stel­len zu haben, muß die erste Ablei­tung aber min­de­stens einen qua­dra­ti­schen Term haben.

Pro­bie­ren wir den­noch zum Spaß und zu Übungs­zwecken ein­fach mal eine Gera­de aus, bre­chen also nach n = 1 ab (Tabel­le „Aus­gleichs­ge­ra­de“):

y = a0 + a1 * x

  • Wir tra­gen wenn mög­lich gut erra­te­ne Start­wer­te für die Para­me­ter a0 und a1 in jeweils eine Zel­le der Tabel­le ein (hier C2 und C3). In eine Spal­te (hier: Spal­te E) neben den Meß­wer­ten tra­gen wir dann die For­mel für die Gera­de ein, mit Bezug auf die Para­me­ter ($C$2 und $C$3) und fül­len die­se Spal­te nach unten auf (obe­re Zel­le und alle Zel­len bis zum Ende der Daten aus­wäh­len und Strg‑D drücken).
  • In einer wei­te­ren Spal­te berech­nen wir dann das Qua­drat des Feh­lers, den der so errech­ne­te Punkt zum tat­säch­lich gemes­se­nen Wert hat (hier: Spal­te F). Auch die­se Spal­te fül­len wir nach unten auf und errech­nen in einer wei­te­ren Zel­le (G2) die Sum­me all die­ser Feh­ler­qua­dra­te. Das Ziel ist nun, a0 und a1 (C2 und C3) so zu bestim­men, daß die Sum­me der Feh­ler­qua­dra­te mini­mal wird.

Von Hand ist das aus­ge­spro­chen müh­sam, aber zum Glück bie­tet Calc dafür eine auto­ma­ti­sche Lösung, den „Sol­ver“. Man star­tet Extras->Solver und gibt die Zel­le mit der Sum­me der Feh­ler­qua­dra­te (hier G2) als Ziel­zel­le an. Der Ziel­wert soll mini­mal wer­den und als ver­än­der­li­che Zel­len geben wir C2 und C3 in der Form $C$2:$C$3 an. Der Sol­ver rech­net nun eini­ge Sekun­den, denn er muss die Wer­te in C2 und C3 solan­ge anpas­sen, bis er einen Mini­mal­wert findet.

Wie die bei­den Dia­gram­me zei­gen, fin­det der Sol­ver nach eini­gen Sekun­den einen plau­si­blen Wert. Die oran­ge dar­ge­stell­ten Gera­den nähern die gemes­se­ne Kur­ve best­mög­lich an. Bei nur zwei Varia­blen und den linea­ren Eigen­schaf­ten einer Gera­den kann man auch sicher sein, ein ech­tes Mini­mum gefun­den zu haben.

Aus­gleichs­kur­ve mit Polynomen

Die Aus­gleichs­ge­ra­de aus dem vori­gen Kapi­tel ist bes­ser als nichts, aber sie befrie­digt nicht wirk­lich. Für bes­se­re rech­ne­ri­sche Annä­he­rung an die Meß­kur­ve kommt man nicht umhin, höher­gra­di­ge Poly­no­me zu ver­wen­den. Der Rechen­weg unter­schei­det sich dabei nicht grund­le­gend von dem der Aus­gleichs­ge­ra­den, nur die For­mel wird etwas kom­pli­zier­ter. Wir neh­men ein­fach wei­te­re Ter­me der Form anxn dazu.

Je grö­ßer n wird, je höher also die Ord­nung des Poly­noms wird, umso genau­er kön­nen wir die Meß­kur­ve annä­hern. Pro­bie­ren wir es der Rei­he nach durch, zunächst mit einem Poly­nom zwei­ten Grades:

Aus den oben schon genann­ten Grün­den ver­wun­dert es nicht, daß das die Meß­kur­ve noch nicht trifft. Wie erwar­tet gibt es nur einen loka­len Extrem­wert, ein Mini­mum bei etwa 20 °C. Ver­su­chen wir es nun mit einem Poly­nom drit­ten Grades:

Damit kom­men wir der Form der gemes­se­nen Kur­ve schon sehr nahe. Gleich­zei­tig wird die Sum­me der Feh­ler­qua­dra­te mit stei­gen­dem Grad des Poly­noms immer kleiner.

Um es auf die Spit­ze zu trei­ben, kann man für den 100 MHz Quarz­os­zil­la­tor auch ein Poly­nom vier­ten und fünf­ten Gra­des aus­pro­bie­ren, das wei­te­re klei­ne Ver­bes­se­run­gen bringt:

Tem­pe­ra­tur­kom­pen­sa­ti­on per Software

Das sind net­te Spie­le­rei­en mit einem Spreadsheet-Pro­gramm, aber was bringt das gan­ze nun? Ganz ein­fach, das so gefun­de­ne Poly­nom gestat­tet uns eine soft­ware­ge­steu­er­te Tem­pe­ra­tur­kom­pen­sa­ti­on. Wir kön­nen die mut­maß­li­che tat­säch­li­che Fre­quenz des Oszil­la­tors bei einer bestimm­ten Tem­pe­ra­tur aus sei­ner Nomi­nal­fre­quenz annähern.

Das zeigt die Tabel­le Poly_n5 für bei­de Quarz­os­zil­la­to­ren. In der Spal­te H wur­de die nach der Kom­pen­sa­ti­on ver­blei­ben­de Abwei­chung zur gemes­se­nen Fre­quenz errech­net und in Spal­te I in ppm (parts per mil­li­on) umge­rech­net. Die ver­blei­ben­de Abwei­chung liegt nun außer an den Tem­pe­ra­tur­gren­zen bei weni­ger als 1 ppm. Die unkom­pen­sier­te Abwei­chung ist zum Ver­gleich in Spal­te J dar­ge­stellt und deren Maxi­mum liegt bei knapp 40 ppm bzw. gut 16 ppm, was für einen mit +/- 25 ppm spe­zi­fi­zier­ten Quarz auch schon recht gut ist.

Man beach­te die unter­schied­li­che Ska­lie­rung der y‑Achsen. Von den extre­men Tem­pe­ra­tur­be­rei­chen abge­se­hen, hat die Tem­pe­ra­tur­kom­pen­sa­ti­on die Fre­quenz­sta­bi­li­tät also min­de­stens um den Fak­tor 10 ver­bes­sert. Die deut­li­chen Abwei­chun­gen bei etwa 25 °C lie­gen an der Meß­me­tho­de, denn die Zim­mer­tem­pe­ra­tur war jeweils der End­punkt der Mes­sun­gen. Hier gab es mut­maß­lich Meß­feh­ler durch den abge­setz­ten Tem­pe­ra­tur­sen­sor auf dem CPU Board, der durch die Eigen­erwär­mung eini­ge Grad mehr anzeigt, als beim Quarz wirk­lich herrschen.

Such­stra­te­gien

Bei der hier gezeig­ten Meß­da­ten­ana­ly­se kommt der Calc-Sol­ver schnell an sei­ne Gren­zen. Am Ende steht ein Poly­nom fünf­ter Ord­nung mit sechs unbe­kann­ten Para­me­tern a0 .. a5. Der Sol­ver muß im Prin­zip alle mög­li­chen Para­me­ter-Wer­te durch­pro­bie­ren, die Qua­drat­sum­me bil­den und ver­su­chen, die­se Sum­me zu mini­mie­ren. Ohne halb­wegs pas­send gewähl­te Start­wer­te ist das ein Unter­fan­gen, das auf heu­ti­gen Rech­nern nicht lös­bar ist. Der Such­raum ist unend­lich groß und der Sol­ver muß ihn auf irgend­ei­ne Wei­se ein­schrän­ken. Sei­ne Suche gilt als been­det, wenn die Qua­drat­sum­me bei jeder klei­nen Ände­rung eines Para­me­ters wie­der grö­ßer wird. Dabei besteht aber immer die Gefahr, ein loka­les Mini­mum zu fin­den, das weit weg vom tat­säch­li­chen Mini­mum liegt. Es kann immer sein, daß hin­ter den Ber­gen bei den sie­ben Zwer­gen ein Mini­mum liegt, das tau­send­mal schö­ner als das bis­her gefun­de­ne ist. Man soll­te dem Sol­ver also durch geeig­ne­te Anfangs­wer­te hel­fen. Wie macht man das?

Suk­zes­si­ve Erhö­hung des Polynom-Grades

Die ein­fach­ste Metho­de, die hier auch ange­wen­det wur­de, ist die suk­zes­si­ve Erhö­hung des Poly­nom Gra­des. Man star­tet mit nied­ri­ger Ord­nung und lässt den Sol­ver dazu pas­sen­de Para­me­ter suchen. Geht man schritt­wei­se zu höhe­ren Ord­nun­gen, dann sind die­se Para­me­ter zwar nicht mehr ganz rich­tig, aber auch nicht ganz falsch. Sie sind ein guter Start­wert. Man kann auch zunächst mal nur den besten neu­en Para­me­ter suchen und erst dann eine wei­te­re Suche über alle bis­he­ri­gen Para­me­ter star­ten. So ver­bes­sert sich das Ergeb­nis nach und nach.

Kur­ven­dis­kus­si­on

Eine Kur­ven­dis­kus­si­on kann beim Erra­ten gün­sti­ger Anfangs­wer­te für die Para­me­ter hel­fen. Zur Kur­ven­dis­kus­si­on benö­tigt man die Ablei­tun­gen der Funk­ti­on. Gehen wir zunächst von einem Poly­nom drit­ten Gra­des aus, so erge­ben sich fol­gen­de Ableitungen:

y = a0 + a1 * x + a2 * x2 + a3 * x3

y‘ = a1 + 2 * a2 * x + 3 * a3 * x2

y“ = 2 * a2 + 6 * a3 * x

Aus der ersten Glei­chung erkennt man sofort, daß a0 = y(0) ist, denn wenn x = 0 ist, dann fal­len alle ande­ren Ter­me weg. a0 ist beim 8 MHz Oszil­la­tor also etwa 212 und beim 100 MHz Oszil­la­tor etwa 355.

Die erste Ablei­tung einer Funk­ti­on ist an loka­len Extrem­wer­ten gleich null. Beim 100 MHz Oszil­la­tor fin­den wir ein loka­les Mini­mum etwa bei x=3 und ein loka­les Maxi­mum bei etwa x=60. Mit der ersten Ablei­tung erhal­ten wir dar­aus zwei Gleichungen:

y‘ = a1 + 2 * a2 * 3 + 3 * a3 * 32 = 0

y‘ = a1 + 2 * a2 * 60 + 3 * a3 * 602 = 0

Aus zwei Glei­chun­gen mit drei Unbe­kann­ten kann man durch Mul­ti­pli­ka­ti­on mit einem geeig­ne­ten Fak­tor und anschlie­ßen­der Addi­ti­on bei­der Glei­chun­gen eine Glei­chung mit zwei Unbe­kann­ten machen. Eli­mi­nie­ren wir zunächst a2 durch Mul­ti­pli­ka­ti­on der ersten Glei­chung mit ‑20 und anschlie­ßen­der Addition:

‑20 * (a1 + 2 * a2 * 3 + 3 * a3 * 32) + a1 + 2 * a2 * 60 + 3 * a3 * 602 = 0

ergibt nach Auflösung:

a1 = 540 * a3

Durch Mul­ti­pli­ka­ti­on mit ‑1 und anschlie­ßen­der Addi­ti­on lässt sich alter­na­tiv a1 eliminieren:

a1 + 2 * a2 * 60 + 3 * a3 * 602 – (a1 + 2 * a2 * 3 + 3 * a3 * 32) = 0

und nach Ausmultiplizieren

a2 = ‑94,5 * a3

Damit haben wir eine direk­te Abhän­gig­keit zwi­schen a1 und a3 bzw. a2 und a3 gefun­den. Die­se bei­den For­meln kann man direkt in die Zel­len für a1 und a2 ein­ge­ben und braucht sie dann nicht mehr in die Liste der Unbe­kann­ten für den Sol­ver auf­zu­neh­men. Der Such­raum ist also wei­ter deut­lich eingeschränkt.

Kön­nen wir noch mehr tun? Ja klar, wir ken­nen nun Nähe­rungs­wer­te für a0, a1 und a2 und kön­nen die in die Glei­chung des Poly­noms einsetzen:

y = 355 + 540 * a3 * x – 94,5 * a3 * x2 + a3 * x3

Ohne zu ver­ges­sen, daß alles sowie­so nur Nähe­rungs­wer­te sind, suchen wir aus der Meß­rei­he einen Wert x und den zuge­hö­ri­gen Wert y her­aus. Dabei ver­su­chen wir einen reprä­sen­ta­ti­ven Wert zu fin­den, also Son­der­fäl­le an den Tem­pe­ra­tur­gren­zen und auch die loka­len Extrem­wer­te zu ver­mei­den. Wie wär’s, mehr oder weni­ger will­kür­lich, mit x=40 und y=773? Mit die­sen Wer­ten kön­nen wir nun einen plau­si­blen Start­wert für a3 ausrechnen:

773 = 355 + 540 * a3 * 40 – 94,5 * a3 * 402 + a3 * 403

Aus die­ser Glei­chung ergibt sich a3 = ‑0,006372 und durch Ein­set­zen die­ses Wer­tes in die oben schon gefun­de­nen Glei­chun­gen a1 = ‑3,44088 und a2 = 0,602154.

Damit haben wir dem Sol­ver das Leben leicht­ge­macht. Wir haben plau­si­ble Start­wer­te für alle Para­me­ter des Poly­noms drit­ter Ord­nung. Ohne den Sol­ver auch nur zu bemü­hen, deckt sich die­se Kur­ve schon erstaun­lich gut mit der Meß­kur­ve (links):

Rechts ist noch­mal zum Ver­gleich das oben schon gezeig­te Dia­gramm des mit dem Sol­ver gefun­de­nen Poly­noms der­sel­ben Ord­nung gezeigt. Über gro­ße Tei­le des Dia­gramms ist die Deckung des von Hand gefun­de­nen Poly­noms sogar bes­ser. Aller­dings ist die Sum­me der Feh­ler­qua­dra­te den­noch höher, denn um den Dia­gramm­teil bei hohen Tem­pe­ra­tu­ren kor­rekt zu model­lie­ren feh­len die höher­gra­di­gen Terme.

Ganz am Schluß kann man dann a0 bis a3 noch­mal durch den Sol­ver ver­bes­sern las­sen. Der fin­det dabei tat­säch­lich eine bes­se­re Lösung, als bei der Suche aufs Gera­te­wohl im ersten Ansatz. Die gering­ste Sum­me der Feh­ler­qua­dra­te ergibt sich bei a0 = 339,27, a1 = ‑3,023, a2 = 0,703 und a3 = ‑0,00858.

Mit dem­sel­ben Vor­ge­hen kann man natür­lich auch den Fall des 8 MHz Oszil­la­tors noch­mal durch­spie­len, was aber hier kei­ne neu­en Erkennt­nis­se mehr bringt. Das sei daher dem geneig­ten Leser überlassen.