Bestimmung von Frequenzänderungen
                      ▀▀▀▀▀▀▀▀▀▀▀▀▀▀▀▀▀▀▀▀▀▀▀▀▀▀▀▀▀▀▀▀▀
                         periodischer Funktionen aus
                         ▀▀▀▀▀▀▀▀▀▀▀▀▀▀▀▀▀▀▀▀▀▀▀▀▀▀▀
                          Messungen phasenbekannter
                          ▀▀▀▀▀▀▀▀▀▀▀▀▀▀▀▀▀▀▀▀▀▀▀▀▀
                                  Zeitpunkte
                                  ▀▀▀▀▀▀▀▀▀▀

                                               Rolf Schröder, 1999-06-07

      Vielfältig treten in der Natur schwingende bzw.  rotierende
      Systeme auf, die uns durch die Exaktheit, mit der sie scheinbar
      ihre Frequenz beibehalten, überraschen.  Denken wir nur an die
      Genauigkeit, mit der die Jahreslänge oder gar die Tageslänge
      eingehalten wird!  Dieses Verhalten hat zum Teil in Form von
      Erhaltungssätzen, wie etwa dem der Drehimpulserhaltung, seinen
      Niederschlag gefunden.

      Aber auch andere Vorgänge laufen mit außerordentlich hoher
      Präzision ab, wie uns verschiedenste Uhrenarten demonstrieren:
      angefangen mit der Pendeluhr über Uhren mit Unruhe bis zur
      Quarzuhr.  Die schärfsten uns bekannten Resonanzen aber treten im
      Bereich der Atom- und Kernphysik auf.

      All die hier beispielhaft aufgeführten Systeme zeigen trotz ihrer
      vom Prinzip her diskret definierten Eigenfrequenz kleine, aber
      dennoch meßbare Frequenzänderungen.  Diese haben ihre Ursache
      einerseits in äußeren Störungen, wie z. B.  die Störung der
      Erdbahn durch Jupiter, die andererseits aber auch auf eine
      Änderung der inneren Eigenschaften zurückgehen, wie z. B.  die
      Veränderungen im Kristallgitter bei Quarzuhren.  Ja selbst bei
      Atomuhren, bei denen die Genauigkeit letztlich durch die
      Heisenbergsche Unschärferelation gegeben ist, wird die theoretisch
      mögliche relative Genauigkeit von ca.  3*10e-25 (Wasserstoffmaser)
      nicht erreicht, da es im Labor unmöglich ist, äußere Störeinflüsse
      auf das Wasserstoffatom völlig auszuschalten.  Praktisch erreicht
      man heute aber eine relative Genauigkeit von etwa 10e-15.

      Die Frequenzen schwingender Systeme sind durch gewisse Eigen-
      schaften, wie Geometrie, Massenverteilung, innere Kräfte
      (Zustandsgleichung!) usw., bestimmt.  Man kann die zeitlichen
      Veränderungen solcher Größen durch Frequenzmessungen mit großer
      Genauigkeit ermitteln.  Noch höhere Genauigkeiten werden erzielt,
      wenn statt einfacher Frequenzmessungen die Zeitpunkte des
      Phasenablaufs der Schwingung registriert werden.

      Diese Methode wollen wir im folgenden ausführlich darstellen.  Wir
      beschränken uns dabei auf die Messung ungedämpft schwingender
      Systeme, deren Schwingungsablauf sich durch periodische Funktionen
      beschreiben läßt.

      Gegeben sei eine über dem Einheitsintervall periodische Funktion
      h, so daß

          h(E) = h(E+1)  für beliebige reelle E ,                    (1)

      mit E als Phase.  Wir wollen den Begriff der periodischen Funktion
      auch für Funktionen verwenden, die im strengen Sinne nicht die
      Bedingung (1) erfüllen.  Wir nennen eine zeitabhängige Funktion
      f(t) periodisch in der Phase, wenn es eine Funktion h nach (1)
      gibt und dazu eine streng monotone stetige Phasenfunktion E(t)
      gibt mit dE/dt ≥ k > 0, k eine Konstante, so daß die Funktion f(t)
      beschrieben werden kann durch:

          f(t) = h(E(t)) .                                           (2)

      Haben wir z.B.  eine in der Zeit periodische Funktion mit einer
      konstanten Frequenz Oo , so gilt für die Phase:

          E(t) = Oo∙(t-to) + Eo .                                    (3)

      Ist die Frequenz zeitabhängig, O = O(t), so können wir leider
      nicht in (3) Oo∙(t-to) durch O(t) ersetzen.  Es würde nämlich
      bedeuten, daß die Phase E zu einem bestimmten Zeitpunkt t1 nur von
      der Frequenz zu diesem Zeitpunkt abhinge.  Das kann aber nicht
      sein, da Änderungen der Frequenz vor diesem Zeitpunkt ebenfalls
      die Phase zur Zeit t1 beeinflussen müssen!  Richtig ist, daß

          dE = O(t)∙dt  ,                                            (4)

      und daß sich daher E durch Integration über die Zeit ergibt:

                 t
                 ⌠
          E(t) = │ O(t')dt' + Eo                                     (5)
                 ⌡
                 to

      mit der Phase Eo = E(to) als Integrationskonstante zum
      (willkürlich vorgegebenen) Zeitpunkt to.

      Unsere Aufgabe besteht nun darin, aus N Zeitpunktmessungen ti
      (z. B.  die Zeiten bestimmter Kurvenpunkte, wie Maxima, Minima
      usw.) der Funktion f(t) zu bekannten Phasen Ei die
      Frequenzfunktion O(t) zu bestimmen.

      Dazu wird O(t) zweckmäßig parametrisiert.  Die Art der Parametri-
      sierung wird wesentlich bestimmt durch die physikalischen Eigen-
      schaften des selbsterregt schwingenden Systems.  Wir wollen hier
      Systeme behandeln, bei denen mit zwei Arten von Frequenzänderun-
      gen zu rechnen ist:  Erstens erwarten wir, daß das System langsame
      Änderungen (wie Entwicklungseffekte oder Ermüdungserscheinungen
      usw.) aufweist, einen Trend.  Zweitens vermuten wir periodische
      Einflüsse im Sinne einer Frequenzmodulation (wie z. B.  die
      Modulation der Erdrotation mit der Jahreslänge).

      Den Trend beschreiben wir durch ein Polynom bis zum Grad NP, die
      Frequenzmodulation beschreiben wir durch eine Fourierreihe bis zur
      Ordnung NF.  Für O(t) erhalten wir daher:

                    (n)      n
                 NP Oo∙(t-to)
          O(t) = Σ  ──────────
                n=o     n!
                                                                     (6)
                   NF
                 + Σ  Am∙cos(2π∙m(t-to)/Po)+Bm∙sin(2π∙m(t-to)/Po) .
                  m=1

      Als Parameter haben wir die Zeitableitungen des Trends im
      Zeitpunkt to von nullter bis NP-ter Ordnung, die Fourierkoeffizi-
      enten für die Glieder erster bis NF-ter Ordnung und die Grund-
      periode Po der Modulation.  Mit (6) erhält man so aus (5):

                    (n)      (n+1)
                 NP Oo∙(t-to)
          E(t) = Σ ────────────── + Eo
                n=o    (n+1)!

                   NF
                 + Σ  Am∙Po/(2π∙m)∙sin(2π∙m(t-to)/Po)                (7)
                  m=1

                 + Bm∙Po/(2π∙m)∙(1-cos(2π∙m(t-to)/Po))  .

      Setzen wir die N Messungen ti zu den Phasen Ei in (7) ein, so
      haben wir ein Gleichungssystem mit N Gleichungen und
      M = NP + 3 + 2∙NF Unbekannten.  Anzustreben ist eine Lösung nach
      der Methode der kleinsten Quadrate, da diese biasfreie Schätzungen
      der Parameter mit minimalen Varianzen liefert.

      Zur direkten Lösung des Gleichungssystems (7) müßten die Messungen
      ti explizit auftreten, und alle Unbekannten müßten linear im
      System vorkommen.  Leider sind die Meßwerte ti implizit in (7)
      vorhanden, und die Unbekannte Po ist nichtlinear in den Gleichun-
      gen enthalten.  Man könnte das System zu lösen versuchen, indem
      man es nach den Meßwerten und den Unbekannten linearisiert, um
      durch Iteration die Lösung zu finden.  Dann hat man das Problem,
      geeignete Anfangswerte zu finden.

      Glücklicherweise ist es aber gar nicht nötig, das System (7) nach
      den Meßwerten ti zu linearisieren.  Wir können nämlich für unsere
      Fälle von der allgemein gut erfüllten Annahme ausgehen, daß die
      Frequenz O(t) relativ kleine Schwankungen um eine feste mittlere
      Frequenz Oo ausführt, die dem Betrage nach eine gewisse Schranke
      C, die wesentlich kleiner als eins (und größer Null) ist, nicht
      übersteigt:

          │O(t)-Oo│
          ───────── < C << 1 .                                       (8)
             Oo

      Aufgrund der daher nahezu linearen Beziehung zwischen den ti und
      den Ei = E(ti) können statt der eigentlichen Meßwerte ti die
      Phasen Ei als Meßwerte angesehen werden, und die ti sind als fest
      gegeben anzusehen.  Mit Hilfe von (4) können wir aus den
      Meßfehlern σ(ti) der ti dann auch einen entsprechenden Fehler
      σ(Ei) der Ei angeben:

          σ(Ei) = Oo∙σ(ti) .                                         (9)

      Auch die Gewichte Pi zur Lösung des Gleichungssystems (7) nach der
      Methode der kleinsten Quadrate mit den Phasen Ei als Messungen
      sind dann bekannt, sie sind proportional zu 1/σ(ti)², da

                      2               2
          Pi = 1/σ(Ei)  = 1/(Oo∙σ(ti))     .                        (10)

      Wäre die Grundperiode Po der Modulation fest gegeben, so könnte
      man (7) jetzt direkt lösen.  Tatsächlich ist aber Po unbekannt.
      Eine bezüglich Po iterative Lösung ist nur möglich, wenn bereits
      ein hinreichend guter Näherungswert bekannt ist.  Wir diskutieren
      daher zunächst andere Lösungswege, und sei es, um wenigsten einen
      guten Anfangswert Poo für eine anschließende Iteration nach Po zu
      finden.  Eine Möglichkeit der Lösung des Gleichungssystems (7)
      besteht darin, daß wir die Anzahl der Fourierkoeffizienten soweit
      erhöhen, bis die Zahl der insgesamt zu bestimmenden Unbekannten
      gleich der Zahl der Gleichungen ist.  Im Falle äquidistanter
      Zeitpunkte ti hätten wir dann eine klassische Fourieranalyse (hier
      jedoch mit zusätzlichen Polynomgliedern).  Die Periode Po wäre
      dann durch

          Po = 2∙NF∙Dt   (Dt: Abstand der Zeitpunkte ti)            (11)

      festgelegt und wir hätten nicht die Qual der Wahl von Po.  Dieses
      Vorgehen würde sogar für nichtäquidistante ti funktionieren (mit
      Dt = (tN-t1)/(N-1) als mittlerem Zeitabstand), solange nur die in
      (7) auftretenden Funktionen bezüglich der Verteilung der ti linear
      unabhängig blieben.

      Dieses Verfahren wäre zweckmäßig, wenn in einer größeren
      Datenmenge unbekannte Frequenzen gesucht werden und keine
      Signifikanzaussage, insbesondere über die Polynomglieder, benötigt
      wird.

      Ist man jedoch an einer physikalisch relevanten Aussage der
      Polynomglieder interessiert, so müssen die Varianzen der
      Ei = E(ti), i=1,N, bekannt sein, so daß mit Hilfe der über (7)
      durch Linearisierung nach den Parametern gegebenen Matrix die
      Varianzen der Parameter (und damit deren Signifikanzen) bestimmt
      werden können.

      Sind jedoch die Varianzen der Ei unbekannt, so können diese
      dennoch mit Hilfe einer Ausgleichsrechnung nach der Methode der
      kleinsten Quadrate geschätzt werden, allerdings muß dazu die Zahl
      der zu bestimmenden Parameter beträchtlich verringert werden, da
      sonst die Signifikanzaussage (in Gestalt der mittleren Fehler)
      wenig verläßlich ist.  Das bedeutet aber auch eine beträchtliche
      Verringerung der Zahl der Fourierparameter, so daß nicht annähernd
      alle über das Meßzeitintervall N*Dt möglichen orthogonalen
      Fourierkomponenten zugleich bestimmbar sind.  Besonders schwer-
      wiegend wird dieses Problem, wenn, wie leider häufig der Fall, die
      Zahl der Meßwerte nicht sehr groß ist.

      Im allgemeinen werden wir daher anders vorgehen müssen.  Wir haben
      zwei durch die Länge des Meßzeitraumes gegebene grundsätzliche
      Schranken für den Bereich der überhaupt bestimmbaren
      Modulationsfrequenzen:  Einerseits können schwerlich Frequenzen
      ermittelt werden, deren Perioden größer als N*Dt sind, da dann die
      Korrelation zwischen Polynom- und Fourierkomponenten gegen eins
      strebt.  Andererseits besagt das Abtasttheorem, daß Perioden
      kleiner als etwa 2*Dt nicht bestimmbar sind.  Wir brauchen also
      nur nach Perioden Po zu suchen, die die folgende Bedingung (12)
      erfüllen:

          N∙Dt ≥ Po ≥ 2∙Dt∙NF  .                                    (12)

      Wir lösen nun unter Vorgabe von Polynomgrad NP und Fourierordnung
      NF das Gleichungssystem (7) mit der Grundperiode Po als unabhängig
      gegebener Variablen.  Wir erhalten so die Schätzung der Varianz
      (bzw.  des mittleren Fehlers) des Meßwertes mit Einheitsgewicht
      als Funktion von Po.  Beim absoluten Minimum dieser Funktion haben
      wir dann die wahrscheinlichste Lösung gefunden.  Ob die Modulation
      mit der Grundperiode Po tatsächlich signifikant ist, läßt sich mit
      einer anschließenden Ausgleichsrechnung unter Einbeziehung von Po
      als Unbekannter aus den Varianzen der Fourierkomponenten
      ermitteln.

                                    * * *

      An dieser Stelle muß der gewissenhafte Statistiker warnend den
      Zeigefinger erheben:  Indem für eine quasi dichte Periodenfolge
      Ausgleichungen durchgeführt werden, zeigt die damit berechnete
      Varianz als Funktion von Po einen um die tatsächliche Varianz
      schwankenden Verlauf, und zwar auch im Falle nicht vorhandener
      periodischer Glieder in den Daten!  D.h., der Erwartungswert
      dieser Funktion ist zwar gleich der Varianz der Meßwerte, das
      Minimum dieser Funktion ist jedoch systematisch kleiner als die
      Varianz der Meßwerte.  Daher werden Perioden vorgetäuscht.  Der
      Elimination solcher nicht signifikanter Perioden werden wir uns
      noch später zuwenden.

                                    * * *

      Praktisch führen wir für eine Folge äquidistanter Modulations-
      grundfrequenzen Oo = 1/Po innerhalb der durch (12) definierten
      Grenzen Ausgleichslösungen nach der Methode der kleinsten Quadrate
      durch, und entnehmen der Lösung mit kleinster Varianz der Meßwerte
      einen Anfangswert Poo für die anschließende iterative Ausgleichs-
      lösung nach allen Parametern (für die (7) nur nach Po linearisiert
      wurde).  Das Ganze muß für eine Reihe von Polynomgraden NP und
      Fourierordnungen NF durchgeführt werden.  Man wählt diejenige
      Lösung aus, für die einerseits die Zahl M der Parameter möglichst
      klein ist, und für die andererseits keine signifikante
      Verkleinerung der Varianz des Meßwertes mit Einheitsgewicht bei
      Vergrößerung von NP bzw.  NF mehr zu erzielen ist.  Wir können
      dann annehmen, das die Frequenz O(t) im Rahmen der Meßfehler nach
      (6) mit entsprechenden NP,NF beschrieben wird.

      Da wir die Zahl M der Unbekannten immer merklich kleiner als die
      Zahl N der Gleichungen haben müssen, muß insbesondere auch NF
      relativ klein sein, d.h., die Modulation muß sich mit wenigen
      Fouriergliedern darstellen lassen.

      Tatsächlich ist das praktisch immer der Fall.  Hätte man z.B.
      eine rechteckförmige Modulation, so müßte diese in (7) durch eine
      unendliche Fourierreihe einer Sägezahnfunktion dargestellt werden.
      Die Amplitudenbeträge der Fourierreihe in (7) nehmen mit
      wachsendem Index m aber erheblich schneller ab als die in (6), da
      (7) durch Integration von (6) erzeugt wurde.  Es bedeutet für uns,
      daß wir im Rahmen der durch die Meßfehler gegebenen Signifikanz-
      grenzen praktisch immer mit wenigen Fouriergliedern auskommen.
      Pessimistisch formuliert bedeutet es aber auch, daß sich die
      Frequenzmodulation grundsätzlich nur bis auf wenige Fourierglieder
      signifikant bestimmen läßt (und damit übrigens auch die zuvor
      erwähnte "große" Fourieranalyse unpraktikabel ist).

      Solange in den Daten keine Modulation außerhalb der Bedingung (12)
      enthalten ist, und solange die Grundkomponente der Fourierreihe
      nicht zu stark mit dem Polynom korreliert ist, sind keine Probleme
      bei der Reduktion zu erwarten.

      Treten jedoch Modulationsperioden in den Daten größer als Po auf,
      so können sie durch das Polynom dargestellt werden, insbesondere,
      wenn sie nur wenige Fourierkomponenten enthalten.  Damit stellen
      aber die Polynomglieder nicht mehr den Trend im obigen Sinne
      unverfälscht dar, eine Trennung von Modulation und Trend ist dann
      prinzipiell unmöglich.  Das ist übrigens auch der Fall, wenn
      Fourierkomponenten durch die Wahl eines zu hohen Polynomgrades NP
      immer besser approximierbar (bzw.  mit Polynomkomponenten immer
      stärker korreliert) sind.

      Treten andererseits Modulationsanteile mit P < 2∙Dt in den Daten
      auf, so wirken diese für die Ausgleichung wie eine zusätzliche
      Streuung der Meßwerte, solange die Meßwerte ti über die Phasen der
      Modulationsanteile statistisch gleichverteilt sind.  Dann ist auch
      eine systematische Verfälschung des Trends nicht zu erwarten.

                                   -------

      Für die Beobachtung weiterer phasenbekannter Zeitpunkte ist es von
      Interesse, diese mit möglichst großer Sicherheit im voraus zu
      berechnen.  Wir wollen also für bekannte Parameter A der
      Frequenzfunktion O(t) und der zum Zeitpunkt to bekannten Phase Eo
      für eine beliebige Phase E den zugehörigen Zeitpunkt t = t(E) mit
      Varianz σ(t) bestimmen.  Aufgrund der Eigenschaft (8) können wir
      uns eine erste Näherung t1 nach

          t1 = to + (E-Eo)/Oo                                       (13)

      verschaffen.  Sodann linearisieren wir E(t) an der Stelle t1,

          E = E(t1) + O(t1)∙(t-t1)  +  Rest                         (14)

      mit O(t1) = dE/dt(t1), so daß man für die nächst bessere Näherung
      t2 erhält:

          t2 = t1 + (E-E(t1))/O(t1)  .                              (15)

      Wir haben so eine Iterationsvorschrift erhalten, und da wegen (8)
      die Ableitung der Funktion E(t) nach der Zeit nahezu konstant ist,
      ist die Lösung t nach wenigen Schritten im Rahmen der
      Rechengenauigkeit erreicht.

      Es bleibt noch die Varianz von t zu bestimmen.  Da uns von der
      Ausgleichslösung des Gleichungssystems (7) nach der Methode der
      kleinsten Quadrate auch die Kovarianzmatrix CA der A und von Eo
      bekannt ist, kann man über (15) die Varianz σ(t) der Zeit t
      erhalten.  Dazu betrachten wir im Ausdruck (15) die Auswirkung
      einer Variation DA aller Parameter A (im folgenden denken wir uns
      Eo mit in den A enthalten) auf t2.  Wir können die daraus
      resultierende Variation Dt2 von t2 durch die Matrizengleichung
      (16) beschreiben:

          Dt2 = DTDA˜͂∙DA   .                                        (16)

      Dabei ist DTDA˜͂ der transponierte Spaltenvektor der partiellen
      Ableitungen von t2 nach den Parametern A, und DA ist der
      Spaltenvektor der Variation von A.  Dann erhält man nach bekannten
      Regeln der Statistik die Kovarianzmatrix Ct von t2 (die in unserem
      Fall ein Skalar ist) aus (16) nach:

          Ct = DTDA˜͂∙CA∙DTDA   .                                    (17)

      Da im Konvergenzfall t2 in t übergeht, erhält man für DTDA:

                        2
          DTDA = -1/O(t) ∙DEDA(t)                                   (18)

      mit DEDA(t) als Spaltenvektor der partiellen Ableitungen von E(t)
      nach den Parametern A zur Zeit t.

      Anmerkung:

      Hat man bereits to so gewählt, daß E(to) = Eo in der Nähe von Null
      liegt, so kann man mit sehr hoher Genauigkeit nach (13) den
      Zeitpunkt t(E = 0) angeben:  t(E = 0) = -Eo/Oo.  Da in der
      Umgebung von E = 0 t(E) praktisch nur von den fehlerbehafteten
      Größen Eo,Oo nach (13) abhängt, ist zur Berechnung der Varianz
      Ct(E = 0) nur die Kovarianzmatrix CEO der beiden Größen Eo,Oo
      nötig.  Mit (dt/dEo,dt/dOo) = DtDEO˜͂ (von (13)) gilt daher
      ebenfalls mit sehr großer Genauigkeit:
      Ct(E = 0) = DtDEO˜͂∙CEO∙DTDEO.

      Zur Beschreibung der Änderung des physikalischen Systems haben wir
      bisher die Frequenzfunktion O(t) benutzt.  Für manche Zwecke - zum
      Teil auch aus Gewohnheit - ist es wünschenswert, statt der
      Frequenz die Periode P(t) zu benutzen.  Speziell am Trend- bzw.
      Polynomanteil der Periodendarstellung ist man interessiert.  Wir
      können eine umkehrbar eindeutige Transformation der einander
      entsprechenden Parameter von O(t) und P(t) einschließlich
      Kovarianzmatrix exakt angeben, d.h., sind von O(t) G Komponenten
      der Taylorreihe bekannt, so können umkehrbar eindeutig G
      Komponenten der Taylorreihe der Funktion P(t) = 1/O(t) bestimmt
      werden; allerdings muß das konstante Glied der Taylorreihe von
      O(t) ungleich Null sein.

      Wer nicht an rechentechnischen Einzelheiten interessiert ist, der
      kann den folgenden Abschnitt auch übergehen.

      Wir gehen von der Grundbeziehung zwischen O(t) und P(t) aus:

          O(t)∙P(t) = 1  .                                          (19)

      Für die n-te Ableitung von (19) erhalten wir für einen Polynomteil
      in O(t) mit G Komponenten, d.h.  mit nullter bis (G-1)-ter
      Zeitableitung:

              (n)       n        (n-i)   (i)
          [O(t)∙P(t)] = Σ (n÷i) ∙ O(t) ∙ P(t) = 0  für n=1, ..,G-1 .(20)
                       i=0
                             ( (n÷i) lies: n über i )

      Wir definieren rechnerkonform:

              (i-1)          (i-1)
          Pi = Po   und  Oi = Oo    ; i=1,... ,G  ,                 (21)

      der Index "o" an den Größen P bzw. O bedeutet, daß sie für den
      Zeitpunkt to gelten.  Die Ausdrücke (19) und (20) lauten damit:

          O1 ∙ P1 = 1 und

           n                                                        (22)
           Σ ((n-1)÷(i-1)) ∙ On-i+1 ∙ Pi = 0  für n=2,... ,G  .
          i=1

      Nun können wir (22) rekursiv lösen:

          P1 = 1/O1   und

                        n-1                                         (23)
          Pn = -(1/O1) ∙ Σ ((n-1)÷(i-1)) ∙ On-i+1 ∙ Pi  für n=2,... ,G .
                        i=1

      Wir wollen (23) nochmals aufschreiben, jedoch mit Betonung der
      funktionellen Verknüpfungen:

          P1 = P1(O1)  = 1/O1               für i=1  und
                                                                    (24)
          Pi = Pi(O1,... ,Oi;P1,... ,Pi-1)  für i>1  .

      Zur Berechnung der Kovarianzmatrix der Pi benötigen wir die
      totalen Ableitungen der Pi nach den Oj, wobei die totale bzw.
      partielle Ableitung mit D bzw.  d geschrieben wird:

                             i-1
          DPi/DOj = dPi/dOj + Σ dPi/dPr ∙ dPr/dOj
                             r=1                                    (25)
                                         0
          für  j <= i  und  0 sonst und  Σ ... = 0  .
                                        r=1

      Die Fehlertransformation erfolgt über:

                i
          DPi = Σ DPi/DOj∙DOj
               j=1

                i            i-1
              = Σ { dPi/dOj + Σ (dPi/dPr)∙(dPr/dOj) }∙DOj           (26)
               j=1           r=1

                        für i = 1,... ,G  ,

      in Matrizenschreibweise (man beachte, daß dPr/dOj = 0 für r < j):

          DP = T∙DO        mit

          Tij = DPi/DOj

                         i-1
              = dPi/dOj + Σ (dPi/dPr)∙(dPr/dOj)  für  j <= i        (27)
                         r=j

                              i-1
              und 0 sonst und  Σ ... = 0  .
                              r=i

      T ist eine untere Dreiecksmatrix.  Da nun aus (23) folgt:

          dPi/dOj = -P1∙((i-1)÷(i-j))∙Pi-j+1  für  j=1,... ,i  und
                                                                    (28)
          dPi/dPr = -P1∙((i-1)÷(r-1))∙Oi-r+1  für  r < i

      können wir (27) auch so schreiben:

                                           i-1
        Tij = -P1∙((i-1)÷(j-1))∙Pi-j+1 - P1∙Σ ((i-1)÷(r-1))∙Oi-r+1∙Trj
                                           r=j
                                                                    (29)
                          0
        für  j <= i  und  Σ... = 0  (und ((i-1)÷(i-j))=((i-1)÷(j-1)) ).
                          0

      Mit der Formel (29) läßt sich T rekursiv aufbauen.  Da von der
      Ausgleichslösung auch die Kovarianzmatrix CO der Oi, i=1,...  ,G
      bekannt ist, kann nun die Kovarianzmatrix der Pj, j=1,...  ,G
      berechnet werden:

          CP = T∙CO∙T˜͂  .                                           (30)

      Man kommt bei geschickter Programmierung mit nur einem Vektor der
      Länge G und einer GxG-Matrix aus, um aus dem O-Vektor und der
      Kovarianzmatrix CO den P-Vektor und die Kovarianzmatrix CP zu
      berechnen.  Dabei gehen allerdings die ursprünglichen O, CO
      verloren.

                                    * * *

      In diesem Abschnitt wollen wir zunächst einen Aspekt bei der
      Untersuchung von Periodenänderungen bei δ-Cep-Sternen aufgreifen:
      wir wollen uns mit der Ursache des periodischen Gliedes in (6)
      befassen, und zwar speziell mit der Entscheidung, ob dieses durch
      die Bewegung des Cepheiden als Mitglied eines Doppelsternes
      verursacht sein kann.

      Zur Abschätzung des Problems nehmen wir eine Kreisbahn an und
      werten entsprechend auch nur das erste Glied der Fourierreihe von
      (6) aus, mit AO bezeichnen wir die Amplitude (AO = √(A1²+B1²)) in
      der Frequenzdarstellung.  Im folgenden bedeute der Index 1 stets
      Masse, Geschwindigkeit usw.  des δ-Cep-Sternes.

      Der Dopplereffekt führt für die wahre Bahngeschwindigkeit v1 des
      Sternes 1 (M1) relativ zum Schwerpunkt zu:

          v1/c = AO∙P/sin(i),  (P=1/O), i=Bahnneigung.              (31)

      Andererseits ist der Bahnradius a1 mit v1 verknüpft:

          v1 = 2π∙a1/U,  so daß                                     (32)

          a1 = AO∙P∙U∙c/(2π∙sin(i)),                                (33)

      mit U als Umlaufszeit des Doppelsterns (gleich Grundperiode Po der
      Fourierkomponente in (6)).  Mit dem Schwerpunktsatz

          a1∙M1 = a2∙M2                                             (34)

      und den Definitionen

          a = a1 + a2    (Abstand Stern 1 - Stern 2)
                                                                    (35)
          M = M1 + M2    (Gesamtmasse des Systems)

      folgt aber:

      a1 = a∙(1-M1/M).                                              (36)

      Es sei darauf hingewiesen, daß a1 sich nach (33) aus der
      Beobachtung ergibt, es muß jedoch für sin(i) ein plausibler Wert
      angenommen werden, da i in unseren Fällen nicht bestimmbar ist.

      Wir wollen nun folgende Maßeinheiten verwenden:

          Zeiteinheit   = siderisches Jahr as = 365.256366 Tage,    (37)
          Längeneinheit = Astronomische Einheit AU = 1.4959787066E11 m,
          Masseneinheit = Sonnenmasse MS = 1.98892E30 kg,

      so daß das 3.  Keplersche Gesetz in diesen Einheiten lautet:

                 3
           2    a
          U  = ───     .                                            (38)
                M

      Im folgenden betrachten wir die Masse M1 als gegeben.  Da U und a1
      als aus der Beobachtung gegebene Parameter anzusehen sind (33),
      können wir aus (36) und (38) eine kubische Gleichung zur
      Bestimmung der Gesamtmasse M des Doppelsternes aufstellen:

                3           2
          (M-M1)  - Mmin ∙ M  =  0   , mit                          (39)

                    3
                  a1                         3
          Mmin = ───── = (AO∙P∙c/(2π∙sin(i))) ∙ U .                 (40)
                   U²
          (Nebenbemerkung:  AO∙P ist dimensionslos;
          c = 2.99792458E08 m/s; c/2π = 10065.30545 AU/as )

      Da stets Mmin > 0 ist, hat (39) nur eine reelle Nullstelle M > 0.
      Bis auf sin(i) sind alle Größen in dem Ausdruck Mmin beobachtbar
      und daher bekannt.  Setzen wir sin(i) = 1, so hat Mmin für
      gegebene AO,P,U seinen minimalen Wert.  Offenbar erhalten wir dann
      aus (39) als Lösung auch das kleinstmögliche M, da für sin(i)
      ungleich 1 der Wert von Mmin größer wird und damit auch ein
      größeres M als Lösung resultiert.

      Zusammenfassend stellen wir fest, daß mit sin(i) = 1 aus (39) eine
      untere Schranke für die Gesamtmasse des Doppelsternsystems folgt.
      Ist diese merklich größer als die Cepheidenmasse M1, so ist die
      Masse nicht mit der Hypothese eines Doppelsternes verträglich, da
      die zweite Komponente dann sichtbar wäre, es sei denn, wir hätten
      es mit einem unsichtbaren Schwarzen Loch zu tun!

      Übrigens können wir sogar von der (unsinnigen) Annahme ausgehen,
      daß M1=0 ist.  Als Lösung von (39) erhalten wir dann als
      Minimalmasse gerade Mmin (40), die kleiner ist als die Lösung M
      bei Annahme M1 > 0.  Ist also bereits Mmin zu groß für einen
      plausiblen Doppelstern, so braucht man (39) nicht mehr zu lösen.

                                    * * *

      An dieser Stelle sollen einige Abschätzungen zur Genauigkeit der
      Frequenz- bzw.  Periodenänderungsmessungen angegeben werden.

      Fall 1:
 
      Nach (7) betrachten wir zunächst mit NP = 1 und NF = 0 den Fall:

                       .
          Ei = Oo∙ti + Oo∙ti²/2 + Eo + vi                           (41)

      mit i=1,Nm Beobachtungen; vi ist der scheinbare Fehler.  Mit den
      Varianzen der Ei:  V(Ei) = σ(Ei)², der Einheitsmatrix I und der
      Annahme, daß alle Varianzen V(Ei) gleich VE sind, ist die
      Kovarianzmatrix der Ei:

          CE = VE∙I   .                                             (42)

      Nach (41) wird der Zusammenhang zwischen dem Meßvektor E und dem                             .
      Parametervektor A =(Oo,Oo,Eo)˜͂ durch die folgende Matrix F
      beschrieben:

              ┌                  ┐
              │ t1   0.5∙t1²  1  │
              │ .        .    .  │
              │ .        .    .  │
          F = │ .        .    .  │                                  (43)
              │ .        .    .  │
              │ .        .    .  │
              │ tNm  0.5∙tNm² 1  │       .
              └                  ┘

      Die Kovarianzmatrix CA der Parameter A ist bei einer Lösung nach
      der Methode der kleinsten Quadrate gegeben durch

          CA = i(F˜͂∙PE∙F)                                           (44)

      (Vorzeichen "i" bedeutet "Inverse von"), dabei ist PE = iCE die
      Gewichtsmatrix der Ei .

      Wir erhalten:

                        ┌      2             3               ┐
                        │   Σti       0.5∙Σti         Σti    │
                        │                                    │
                        │        3           4           2   │
          iCA = 1/VE  ∙ │ 0.5∙Σti     .25∙Σti     0.5∙Σti    │      (45)
                        │                                    │
                        │                    2               │
                        │     Σti     0.5∙Σti         Nm     │    .
                        └                                    ┘

      Wählt man die ti symmetrisch zum Nullpunkt und äquidistant
      verteilt, so sind die Summen über ungerade Potenzen der ti Null.
      Bezeichnen wir den Beobachtungszeitraum mit T = tNm - t1, und
      schreibt man genähert Σti² = Nm∙T²/12, Σ(ti²)² = Nm∙(T²)²/80, so
      erhält man für CA :

                            ┌     2                      ┐
                            │ 12∙T        0         0    │
                            │                            │
                       4    │                        2   │
          CA = VE/(Nm∙T ) ∙ │    0       720    -30∙T    │          (46)
                            │                            │
                            │              2       4     │
                            │    0    -30∙T     9∙T /4   │    .
                            └                            ┘

      Da mit (8) und (9) σ(Ei) = Oo∙σ(ti) ist, erhalten wir schließlich
      für die Varianz bzw.  Streuung der Frequenzänderung:

            .
          V(Oo) = 720∙Oo²∙V(ti)/(Nm∙(T²)²)  bzw.                    (47)

            .
          σ(Oo) = Oo∙σ(ti)/T²∙√(720/Nm)  .                          (48)

          .          .
      Da dPo = -Po²∙dOo erhalten wir auch die Streuung der Perioden-
      änderung:

            .
          σ(Po) = (Po∙σ(ti)/T²)∙√(720/Nm)  ;                        (49)

                       ∙
          (genähert: σ(Po) = 27∙(Po∙σ(ti)/T²)/√Nm)   .             (49a)

      Fall 2:
 
      Zur Anwendung der Formel (49) ist es nützlich, aus der
      Meßgenauigkeit σ(fi) der Funktionswerte der periodischen Funktion
      sowie aus Amplitude Af, Periode Po und der Anzahl Ne von Messungen
      die Genauigkeit σ(ti) eines einzelnen Phasenpunktes abzuschätzen.
      Dabei sei die Änderung der Periode über eine oder wenige Phasen
      (über die gemessen wird) so klein, daß man P(t) = Po setzen darf.
      Als periodische Funktion nehmen wir eine einfache Sinusfunktion,
      so daß wir die folgenden Meßgleichungen haben:

          fi = a1 + a2∙cos(2π∙ti/Po) + a3∙sin(2π∙ti/Po) + vi  .     (50)

      Die Gewichtsmatrix Pf der Messungen fi ist gleich

          Pf = (1/V(fi))∙I   ,                                      (51)

      die Bedeutung von V und I ist dieselbe wie im Fall 1 .  Ganz
      analog zu Fall 1 haben wir die Matrix F :

              ┌                                         ┐
              │  1    cos(2π∙t1/Po)     sin(2π∙t1/Po)   │
              │  .         .                 .          │
              │  .         .                 .          │
          F = │  .         .                 .          │           (52)
              │  .         .                 .          │
              │  .         .                 .          │
              │  1    cos(2π∙tNe/Po)    sin(2π∙tNe/Po)  │    .
              └                                         ┘

      Wir nehmen an, daß die Meßorte äquidistant gerade so über den
      Meßzeitraum der Länge g*Po (mit g > 0 und ganzzahlig) verteilt
      sind, daß

          tNe-t1 = g∙Po∙(Ne-1)/Ne  mit  Ne > 2g   .                 (53)

      Die Inverse der Kovarianzmatrix Ca der Parameter a1,a2,a3 lautet
      dann:

                            ┌                 ┐
                            │  1     0     0  │
                            │                 │
          iCa = Ne/V(fi) ∙  │  0    1/2    0  │                     (54)
                            │                 │
                            │  0     0    1/2 │     ,
                            └                 ┘

      und damit erhalten wir für Ca :

                           ┌                 ┐
                           │  1     0     0  │
                           │                 │
          Ca = V(fi)/Ne ∙  │  0     2     0  │                      (55)
                           │                 │
                           │  0     0     2  │     .
                           └                 ┘

      Die Lösung a2,a3 kann in der Form

          a2 = Af∙cos(2π∙te/Po)                                     (56)
          a3 = Af∙sin(2π∙te/Po)

      geschrieben werden, so daß sich nach (50) auch schreiben läßt:

          fi = a1 + Af∙cos(2π∙te/Po-2π∙ti/Po) + vi  .               (57)

      te beschreibt gerade den Zeitpunkt, an dem die Funktion ein
      Maximum hat.  Wir können ohne Einschränkung der Allgemeinheit aus
      (56) die Streuung σ(te) bestimmen, die identisch ist mit der
      gesuchten Streuung σ(ti) in (49),(49a).  Es folgt aus (56) :

          te = Po/2π∙arctg(a3/a2)    .                              (58)

      Da a2 und a3 nicht korreliert sind, kann die Streuung σ(te) über
      das einfache Fehlerfortpflanzungsgesetz berechnet werden.  Man
      erhält so mit den aus der Matrix (55) zu entnehmenden Streuungen

          σ(a2) = σ(a3) = σ(fi)∙√(2/Ne)                             (59)

      schließlich:

          σ(te) = (σ(fi)/Af)∙(Po/π)/√(2Ne)    .                     (60)

      Ersetzen wir in (49) σ(ti) durch (60), so erhalten wir:

            .
          σ(Po) = (σ(fi)/Af)∙(Po/T)²∙(√360/π)/√(Nm∙Ne) ;            (61)
                       .
          (genähert: σ(Po) = 6∙(σ(fi)/Af)∙(Po/T)²/√(Nm∙Ne)).       (61a)