In der Entwicklungsphase des noch zu beschreibenden Frequenzzählers habe ich einige Meßreihen erstellt. Eine davon war der Frequenzgang eines 8 MHz Quarzes auf dem CPU-Board, die andere der Frequenzgang eines 100 MHz Quarzoszillators auf dem Frequenzzähler selbst. Beide Frequenzgänge wurden über die Temperatur gemessen, die mit Kältespray (Butangas, nicht Rauchen, kein offenes Feuer aber offene Fenster!) bzw. einem Föhn geändert wurde.
In beiden Fällen wurde die Temperatur mit dem Temperatursensor auf dem CPU-Board gemessen, einem TMP275. Zwischen ‑40 °C und +125 °C hat er eine Toleranz von +/- 1 °C. Die tatsächliche Temperatur des Quarzes und besonders des Quarzoszillators dürften um einige °C von der gemessenen Temperatur abweichen, weil sie etwas entfernt davon montiert sind. Es soll hier aber nicht um die Präzision der Messung gehen, sondern einzig und allein um die Auswertung mit Libreoffice Calc. Der Frequenzzähler speichert die Meßdaten in einem CSV-File ab, das mit LibreOffice direkt eingelesen werden kann. Die beiden OpenOffice Calc Dateien stehen als Referenz und für eigene Experimente zum Download bereit:
In der ersten Tabelle der jeweiligen Datei stehen die original Meßdaten. Die weiteren Tabellen sind nachfolgend beschrieben.
Auswertung
Zur Auswertung empfiehlt es sich, zunächst eine Pivot-Tabelle („Daten->Pivot-Tabelle->Einfügen) zu erstellen. Als Zeilenfelder wählt man die Temperatur und als Datenfelder die gemessene Frequenz Fcheck. Für die Datenfelder interessiert uns aber nicht die Summe, sondern der Mittelwert. Das muß vor dem Erzeugen der Pivot-Tabelle eingestellt werden. Nun sucht Calc alle in der Tabelle vorkommenden Temperaturen, sortiert sie aufsteigend und ordnet jeder Temperatur den Mittelwert aller für diese Temperatur gefundenen Meßwerte zu. So erhalten wir eine Tabelle der gemittelten Meßwerte von etwa ‑30 °C bis +70 °C (Tabelle: „Pivot-Tabelle“).
Da wir einen 8 MHz bzw. einen 100 MHz Quarzoszillator gemessen haben, liegen alle Meßwerte ziemlich nahe bei 8.000.000 Hz oder 100.000.000 Hz. Ohne vernünftige Skalierung sieht man da nur eine horizontale Linie über der Temperatur. Damit die Meßwerte sinnvoll dargestellt werden können, erstellt man in diesem Fall am besten eine Spalte mit der Differenz der Meßwerte zum Nominalwert, denn alle Werte weichen nur sehr wenig davon ab. Dieses Delta läßt sich dann grafisch über der Temperatur darstellen.
Die linken Diagramme zeigen die unveränderten Originaldaten. Man erkennt deutlich den Temperaturgang, aber die vielen Ausreißer sind erschreckend. Die Originaldaten enthalten die sehr schnelle Abkühlphase mit Kältespray und die fast genauso schnelle Aufheizphase mit dem Föhn. Daran schließt sich jeweils eine Phase langsamen Aufheizens bzw. langsamen Abkühlens auf Zimmertemperatur an.
Die rechten Diagramme zeigen nur die langsamen Phasen und die Ausreißer sind hier komplett verschwunden (Tabellen: „Auswahl“ und „Pivot-Tabelle_1“). Die Kurvenform bleibt unverändert. Offensichtlich tun in beiden Fällen die schnellen Temperaturänderungen nicht gut, die Schwingung reißt vermutlich kurzzeitig ab.
Erzeugen einer Ausgleichskurve
Alle weiteren Auswertungen erfolgen nun mit den beiden „langsamen Kurven“, denn sie spiegeln offensichtlich die wahren Verhältnisse. Das Ziel ist nun, die jeweilige Kurve so zu glätten, daß die geglättete Kurve möglichst genau der gemessenen Kurve entspricht, ohne deren Meßabweichungen abzubilden. Die Standardmethode zum Erzeugen einer sogenannten Ausgleichskurve ist die Ausgleichsrechnung. Für ein mathematischen Lösungsverfahren muß dieser Ausgleichsrechnung ein Modell zugrunde liegen. Im allgemeinen Fall ist das ein Polynom der Form
y = a0 + a1 * x + a2 * x2 + … + an * xn
In der Praxis muß man diese Reihe natürlich frühzeitig irgendwo abbrechen, denn sonst wird der Rechenaufwand immer größer. Es sei vorab verraten, daß bei den gezeigten Kurven mindestens ein Polynom dritten Grades nötig ist, also n = 3. Das liegt an den zwei lokalen Minima bzw. Maxima (jeweils grob bei 0 °C und bei 60 °C), bei denen die erste Ableitung des Polynoms null sein muß. Um zwei Nullstellen zu haben, muß die erste Ableitung aber mindestens einen quadratischen Term haben.
Probieren wir dennoch zum Spaß und zu Übungszwecken einfach mal eine Gerade aus, brechen also nach n = 1 ab (Tabelle „Ausgleichsgerade“):
y = a0 + a1 * x
- Wir tragen wenn möglich gut erratene Startwerte für die Parameter a0 und a1 in jeweils eine Zelle der Tabelle ein (hier C2 und C3). In eine Spalte (hier: Spalte E) neben den Meßwerten tragen wir dann die Formel für die Gerade ein, mit Bezug auf die Parameter ($C$2 und $C$3) und füllen diese Spalte nach unten auf (obere Zelle und alle Zellen bis zum Ende der Daten auswählen und Strg‑D drücken).
- In einer weiteren Spalte berechnen wir dann das Quadrat des Fehlers, den der so errechnete Punkt zum tatsächlich gemessenen Wert hat (hier: Spalte F). Auch diese Spalte füllen wir nach unten auf und errechnen in einer weiteren Zelle (G2) die Summe all dieser Fehlerquadrate. Das Ziel ist nun, a0 und a1 (C2 und C3) so zu bestimmen, daß die Summe der Fehlerquadrate minimal wird.
Von Hand ist das ausgesprochen mühsam, aber zum Glück bietet Calc dafür eine automatische Lösung, den „Solver“. Man startet Extras->Solver und gibt die Zelle mit der Summe der Fehlerquadrate (hier G2) als Zielzelle an. Der Zielwert soll minimal werden und als veränderliche Zellen geben wir C2 und C3 in der Form $C$2:$C$3 an. Der Solver rechnet nun einige Sekunden, denn er muss die Werte in C2 und C3 solange anpassen, bis er einen Minimalwert findet.
Wie die beiden Diagramme zeigen, findet der Solver nach einigen Sekunden einen plausiblen Wert. Die orange dargestellten Geraden nähern die gemessene Kurve bestmöglich an. Bei nur zwei Variablen und den linearen Eigenschaften einer Geraden kann man auch sicher sein, ein echtes Minimum gefunden zu haben.
Ausgleichskurve mit Polynomen
Die Ausgleichsgerade aus dem vorigen Kapitel ist besser als nichts, aber sie befriedigt nicht wirklich. Für bessere rechnerische Annäherung an die Meßkurve kommt man nicht umhin, höhergradige Polynome zu verwenden. Der Rechenweg unterscheidet sich dabei nicht grundlegend von dem der Ausgleichsgeraden, nur die Formel wird etwas komplizierter. Wir nehmen einfach weitere Terme der Form anxn dazu.
Je größer n wird, je höher also die Ordnung des Polynoms wird, umso genauer können wir die Meßkurve annähern. Probieren wir es der Reihe nach durch, zunächst mit einem Polynom zweiten Grades:
Aus den oben schon genannten Gründen verwundert es nicht, daß das die Meßkurve noch nicht trifft. Wie erwartet gibt es nur einen lokalen Extremwert, ein Minimum bei etwa 20 °C. Versuchen wir es nun mit einem Polynom dritten Grades:
Damit kommen wir der Form der gemessenen Kurve schon sehr nahe. Gleichzeitig wird die Summe der Fehlerquadrate mit steigendem Grad des Polynoms immer kleiner.
Um es auf die Spitze zu treiben, kann man für den 100 MHz Quarzoszillator auch ein Polynom vierten und fünften Grades ausprobieren, das weitere kleine Verbesserungen bringt:
Temperaturkompensation per Software
Das sind nette Spielereien mit einem Spreadsheet-Programm, aber was bringt das ganze nun? Ganz einfach, das so gefundene Polynom gestattet uns eine softwaregesteuerte Temperaturkompensation. Wir können die mutmaßliche tatsächliche Frequenz des Oszillators bei einer bestimmten Temperatur aus seiner Nominalfrequenz annähern.
Das zeigt die Tabelle Poly_n5 für beide Quarzoszillatoren. In der Spalte H wurde die nach der Kompensation verbleibende Abweichung zur gemessenen Frequenz errechnet und in Spalte I in ppm (parts per million) umgerechnet. Die verbleibende Abweichung liegt nun außer an den Temperaturgrenzen bei weniger als 1 ppm. Die unkompensierte Abweichung ist zum Vergleich in Spalte J dargestellt und deren Maximum liegt bei knapp 40 ppm bzw. gut 16 ppm, was für einen mit +/- 25 ppm spezifizierten Quarz auch schon recht gut ist.
Man beachte die unterschiedliche Skalierung der y‑Achsen. Von den extremen Temperaturbereichen abgesehen, hat die Temperaturkompensation die Frequenzstabilität also mindestens um den Faktor 10 verbessert. Die deutlichen Abweichungen bei etwa 25 °C liegen an der Meßmethode, denn die Zimmertemperatur war jeweils der Endpunkt der Messungen. Hier gab es mutmaßlich Meßfehler durch den abgesetzten Temperatursensor auf dem CPU Board, der durch die Eigenerwärmung einige Grad mehr anzeigt, als beim Quarz wirklich herrschen.
Suchstrategien
Bei der hier gezeigten Meßdatenanalyse kommt der Calc-Solver schnell an seine Grenzen. Am Ende steht ein Polynom fünfter Ordnung mit sechs unbekannten Parametern a0 .. a5. Der Solver muß im Prinzip alle möglichen Parameter-Werte durchprobieren, die Quadratsumme bilden und versuchen, diese Summe zu minimieren. Ohne halbwegs passend gewählte Startwerte ist das ein Unterfangen, das auf heutigen Rechnern nicht lösbar ist. Der Suchraum ist unendlich groß und der Solver muß ihn auf irgendeine Weise einschränken. Seine Suche gilt als beendet, wenn die Quadratsumme bei jeder kleinen Änderung eines Parameters wieder größer wird. Dabei besteht aber immer die Gefahr, ein lokales Minimum zu finden, das weit weg vom tatsächlichen Minimum liegt. Es kann immer sein, daß hinter den Bergen bei den sieben Zwergen ein Minimum liegt, das tausendmal schöner als das bisher gefundene ist. Man sollte dem Solver also durch geeignete Anfangswerte helfen. Wie macht man das?
Sukzessive Erhöhung des Polynom-Grades
Die einfachste Methode, die hier auch angewendet wurde, ist die sukzessive Erhöhung des Polynom Grades. Man startet mit niedriger Ordnung und lässt den Solver dazu passende Parameter suchen. Geht man schrittweise zu höheren Ordnungen, dann sind diese Parameter zwar nicht mehr ganz richtig, aber auch nicht ganz falsch. Sie sind ein guter Startwert. Man kann auch zunächst mal nur den besten neuen Parameter suchen und erst dann eine weitere Suche über alle bisherigen Parameter starten. So verbessert sich das Ergebnis nach und nach.
Kurvendiskussion
Eine Kurvendiskussion kann beim Erraten günstiger Anfangswerte für die Parameter helfen. Zur Kurvendiskussion benötigt man die Ableitungen der Funktion. Gehen wir zunächst von einem Polynom dritten Grades aus, so ergeben sich folgende 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 Gleichung erkennt man sofort, daß a0 = y(0) ist, denn wenn x = 0 ist, dann fallen alle anderen Terme weg. a0 ist beim 8 MHz Oszillator also etwa 212 und beim 100 MHz Oszillator etwa 355.
Die erste Ableitung einer Funktion ist an lokalen Extremwerten gleich null. Beim 100 MHz Oszillator finden wir ein lokales Minimum etwa bei x=3 und ein lokales Maximum bei etwa x=60. Mit der ersten Ableitung erhalten wir daraus zwei Gleichungen:
y‘ = a1 + 2 * a2 * 3 + 3 * a3 * 32 = 0
y‘ = a1 + 2 * a2 * 60 + 3 * a3 * 602 = 0
Aus zwei Gleichungen mit drei Unbekannten kann man durch Multiplikation mit einem geeigneten Faktor und anschließender Addition beider Gleichungen eine Gleichung mit zwei Unbekannten machen. Eliminieren wir zunächst a2 durch Multiplikation der ersten Gleichung mit ‑20 und anschließender 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 Multiplikation mit ‑1 und anschließender Addition lässt sich alternativ 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 direkte Abhängigkeit zwischen a1 und a3 bzw. a2 und a3 gefunden. Diese beiden Formeln kann man direkt in die Zellen für a1 und a2 eingeben und braucht sie dann nicht mehr in die Liste der Unbekannten für den Solver aufzunehmen. Der Suchraum ist also weiter deutlich eingeschränkt.
Können wir noch mehr tun? Ja klar, wir kennen nun Näherungswerte für a0, a1 und a2 und können die in die Gleichung des Polynoms einsetzen:
y = 355 + 540 * a3 * x – 94,5 * a3 * x2 + a3 * x3
Ohne zu vergessen, daß alles sowieso nur Näherungswerte sind, suchen wir aus der Meßreihe einen Wert x und den zugehörigen Wert y heraus. Dabei versuchen wir einen repräsentativen Wert zu finden, also Sonderfälle an den Temperaturgrenzen und auch die lokalen Extremwerte zu vermeiden. Wie wär’s, mehr oder weniger willkürlich, mit x=40 und y=773? Mit diesen Werten können wir nun einen plausiblen Startwert für a3 ausrechnen:
773 = 355 + 540 * a3 * 40 – 94,5 * a3 * 402 + a3 * 403
Aus dieser Gleichung ergibt sich a3 = ‑0,006372 und durch Einsetzen dieses Wertes in die oben schon gefundenen Gleichungen a1 = ‑3,44088 und a2 = 0,602154.
Damit haben wir dem Solver das Leben leichtgemacht. Wir haben plausible Startwerte für alle Parameter des Polynoms dritter Ordnung. Ohne den Solver auch nur zu bemühen, deckt sich diese Kurve schon erstaunlich gut mit der Meßkurve (links):
Rechts ist nochmal zum Vergleich das oben schon gezeigte Diagramm des mit dem Solver gefundenen Polynoms derselben Ordnung gezeigt. Über große Teile des Diagramms ist die Deckung des von Hand gefundenen Polynoms sogar besser. Allerdings ist die Summe der Fehlerquadrate dennoch höher, denn um den Diagrammteil bei hohen Temperaturen korrekt zu modellieren fehlen die höhergradigen Terme.
Ganz am Schluß kann man dann a0 bis a3 nochmal durch den Solver verbessern lassen. Der findet dabei tatsächlich eine bessere Lösung, als bei der Suche aufs Geratewohl im ersten Ansatz. Die geringste Summe der Fehlerquadrate ergibt sich bei a0 = 339,27, a1 = ‑3,023, a2 = 0,703 und a3 = ‑0,00858.
Mit demselben Vorgehen kann man natürlich auch den Fall des 8 MHz Oszillators nochmal durchspielen, was aber hier keine neuen Erkenntnisse mehr bringt. Das sei daher dem geneigten Leser überlassen.