Touchstone Dateien, üblicherweise mit der Endung .s1p oder .s2p sind lesbare Textdateien, die Netzwerk-Parameter enthalten. Diese Dateien sind zu einem Standard geworden. Sie können von vielen Programmen der HF-Meß- und Simulationstechnik gelesen und geschrieben werden. Mitunter ist es praktisch, sie auch mit einem Spreadsheetprogramm zu bearbeiten. Hier soll an einem einfachen Beispiel, den mit einem VNWA gemessenen s11 Parametern, die prinzipielle Vorgehensweise beschrieben werden.
Ein reales Beispiel
Hier ist das Ergebnis der s11-Messung einer kernlosen Zylinderspule. Die VNWA Betriebssoftware zeigt folgendes an:
s11-Messung an einer kernlosen Zylinderspule
Die Meßergebnisse können als Touchstone-Datei exportiert werden:
Hier wurde als Format Real und Imaginärteil gewählt. Andere Formate (Magnitude und Winkel bzw. dB und Phase) sind auch wählbar und können genauso gut weiterverarbeitet werden. Die ersten Zeilen sehen folgendermaßen aus:
! ListType=Lin
# MHz S RI R 50
5.0000000 0.0273505 0.9908113
5.0237559 0.0325082 0.9905086
Die erste Zeile startet mit einem „!“ und ist ein beliebiger Kommentar. Das „#“ in der zweiten Zeile kennzeichnet die Options-Zeile. In diesem Fall besagt sie, daß Frequenzen (erste Spalte) in MHz angegeben sind. Das „S“ bedeutet, daß s‑Parameter folgen und zwar im Format RI, also Realteil (zweite Spalte) und Imaginärteil (dritte Spalte). R kennzeichnet den Referenzwiderstand in Ohm, in diesem und den meisten anderen Fällen 50 Ohm. Dann folgen beliebig viele Zeilen mit s‑Parametern im beschriebenen Format. Andere Formate sollen hier nicht besprochen werden. Sie können der oben genannten Spezifikation entnommen werden.
Bevor man diese Datei nun mit einem Spreadsheet-Programm, wie z.B. LibreOffice, weiterverarbeiten kann, muß man sie in ein importfähiges Format umwandeln. Das geht am einfachsten, indem man sie mit einem beliebigen Text-Editor in ein CSV-Format umwandelt. CSV erwartet im einfachsten Fall eine Beschreibung der nachfolgenden Spalten in der ersten Zeile, gefolgt von den Daten. Spaltenelemente werden am besten durch „;“ getrennt, andere Trennzeichen sind aber auch möglich. Dann müssen noch auf den meisten europäischen PCs die Dezimaltrennzeichen von Punkt auf Komma geändert werden. Die importierbare CSV-Datei sieht dann so aus:
In LibreOffice importiert sieht das dann folgendermaßen aus:
s1p-File in LibreOffice importiert
Damit kann man arbeiten! Da LibreOffice mit einigen eingebauten Funktionen auch komplexe Zahlen bearbeiten kann, ist es hilfreich, die s11-Parameter zunächst in eine komplexe Zahl umzuwandeln. Das geht mit der Funktion Komplexe(real;imag;„j“) in Spalte D, wie hier für die erste Zeile gezeigt:
=KOMPLEXE(B2;C2;"j")
„real“ ist der Realteil und „imag“ der Imaginärteil der zu generierenden komplexen Zahl. Der dritte Parameter gibt an, wie die imaginäre Einheit genannt werden soll. In der Elektrotechnik wird normalerweise ein „j“ gewählt.
Da alle weiteren Berechnungen auf der komplexen Impedanz Z beruhen, sollte diese als nächstes berechnet werden. Das geht über folgende Formel:
1 + s11
Z = Z0 * ---------
1 - s11
In Libre Office wird die Formel dann folgendermaßen in Spalte E codiert:
=IMPRODUKT(50;IMDIV(IMSUMME(1;D2);IMSUB(1;D2)))
Die Bezeichnungen der Funktionen sind eigentlich selbsterklärend: IMSUMME() und IMSUB() berechnen die Summe bzw. die Differenz zweier komplexer Zahlen, IMDIV() den Quotienten und IMPRODUKT() das Produkt. Die Arbeit mit komplexen Zahlen wir damit zum Kinderspiel. Das Spreadsheet sieht nun folgendermaßen aus:
s1p-File mit s11 und Z als komplexen Zahlen
Aus der Impedanz und deren Komponenten X und R lassen sich nun wie in diesem Beitrag zusammengefasst weitere Parameter berechnen, z.B. die Induktivität L = X/ω und die Güte Q = X/R. Die dazu verwendeten LibreOffice-Funktionen sind:
Induktivität L [µH]: =IMAGINÄRTEIL(E2)/(2*PI()*A2)
Der Term „2∗PI()∗A2“ im Nenner entspricht dabei „2∗PI()∗f“, also ω. Da die Frequenz im MHz angegeben ist, wird die Induktivität ohne weitere Umrechnung in µH ausgegeben.
Güte Q: =ABS(IMAGINÄRTEIL(E2)/IMREALTEIL(E2))
Damit die Güte auch im kapazitiven Bereich positiv bleibt, wurde hier noch die ABS() Funktion verwendet.
Die Berechnung weiterer Parameter und die Weiterverarbeitung beispielsweise zum Glätten der Kurven und zur graphischen Darstellung sei dem geneigten Leser überlassen. Hier ist das Libre Office File zum Experimentieren:
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.
Originaldaten des 8 MHz Quarzoszillators
nur langsame Phase
Originaldaten des 100 MHz Quarzoszillators
nur langsame Phase
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.
8 MHz Quarzoszillator mit Ausgleichsgerade
100 MHz Quarzoszillator mit Ausgleichsgerade
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:
8 MHz Quarzoszillator mit Polynom zweiten Grades
100 MHz Quarzoszillator mit 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:
8 MHz Quarzoszillator mit Polynom dritten Grades
100 MHz Quarzoszillator mit 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:
8 MHz Quarzoszillator mit Polynom vierten Grades
100 MHz Quarzoszillator mit Polynom vierten Grades
8 MHz Quarzoszillator mit Polynom fünften Grades
100 MHz Quarzoszillator mit Polynom fünften Grades
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.
8 MHz Quarzoszillator unkompensiert
8 MHz Quarzoszillator temperaturkompensiert
100 MHz Quarzoszillator unkompensiert
100 MHz Quarzoszillator temperaturkompensiert
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:
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:
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):
100 MHz Quarzoszillator mit Polynom dritten Grades aus der Kurvendiskussion
100 MHz Quarzoszillator mit Polynom dritten Grades
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.
