AstronomieAntwoorden
Frequentiebepaling

\( \DeclareMathOperator{\sinc}{sinc} \DeclareMathOperator{\sgn}{sgn} \)

Deze bladzijde legt een paar relatief snelle manieren uit om de eigenschappen (zoals de frequentie) te schatten van periodieke componenten in een meetreeks.

1. Inleiding

Fig. 1: Meetreeksgrafiek

Fig. 2: Meetreeksgrafiek met proefreeks

2. Gebruik een fourierspectrum

Een snelle methode om inzicht te krijgen is om het fourierspectrum van de meetreeks te bepalen. Hiervoor kan een FFT-algoritme gebruikt worden, als er telkens even veel tijd tussen opeenvolgende metingen zat. (FFT betekent Fast Fourier Transform = Snelle Fouriertransformatie.) Als een programmeertaal daar een functie voor heeft dan heeft die meestal een naam zoals POWER of FFT. Het fourierspectrum geeft aan hoeveel "energie" er in de data zit voor elke hele frequentie (dwz. voor sinussen en cosinussen die precies een heel aantal keren in de meetperiode passen). Die energie is gelijk aan het kwadraat van de amplitude. Volgens de wiskunde is het mogelijk om elke willekeurige meetreeks op die manier in componenten te verdelen. Het is voor een frequentiespectrum rekentechnisch een stuk beter als de meetreeks een gemiddelde waarde gelijk aan 0 heeft, dus als dat niet zo is dan berekenen we eerst de gemiddelde waarde en trekken die dan van de meetreeks af.

Fig. 3: Frequentiespectrum

Het is wel duidelijk dat de grootste energie zit bij een frequentie van 10, maar dat er ook bij andere frequenties wat te doen is. Rond een frequentie van 80 is ook nog iets met een piek en een dal te zien. We hadden boven al gezien dat een frequentie van precies 10 niet goed genoeg past. Een fourierspectrum kijkt alleen naar harmonische componenten die een heel aantal keren in de meetreeks passen, en het zou wel erg toevallig zijn als je meetreeks precies één harmonische component (een combinatie van een sinus en een cosinus met dezelfde frequentie) bevat die ook nog eens precies een heel aantal keren in de meetreeks past, dus zelfs als er eigenlijk maar één periode in de meetreeks zit is de kans groot dat een fourierspectrum bij elke frequentie wat laat zien.

Het voordeel van een fourierspectrum is dat het relatief snel te berekenen is. De nadelen zijn dat een fourierspectrum een harmonische component alleen bij precies één frequentie aangeeft als die component toevallig precies een heel aantal keren in de waarneemperiode past, en dat verschillende harmonische componenten elkaar meestal verstoren in zo'n fourierspectrum, zoals de piek bij frequentie 80 die er nogal lelijk uit ziet met ook een stuk "deuk" erbij, omdat er rond die frequentie ook nog een flinke bijdrage is van de piek bij frequentie 10.

3. Eenvoudige schattingen

3.1. Eenvoudigste schatting: Methode 1

De eenvoudigste schatting van de frequentie van een harmonische component met behulp van een fourierspectrum is om alle pieken te vinden (waar de energie groter is dan bij de buren) en daarvan de frequentie en amplitude (de vierkantswortel van de energie) te bepalen.

Voor onze meetreeks levert dit

# 

dus er is ook nog een piek bij frequentie 32. Deze piek is in de grafiek veel moeilijker te zien dan de piek bij frequentie 78, maar toch is de schatting voor de amplitude van de piek bij 32 een stuk hoger dan die bij 78. Dat komt natuurlijk vanwege de "helling" van de piek van frequentie 10.

3.2. Betere schatting: Methode 2

Als je uitrekent wat het fourierspectrum van een enkele harmonische component met positieve frequentie \( f_0 \) is die misschien niet een heel aantal keren in de meetperiode past (dus ook als \( f_0 \) geen geheel getal is) dan vind je dat de energie \( P_f \) bij frequentie \( f \) ongeveer gelijk is aan

\begin{equation} P_f = A^2 (\sinc^2_{π}(f - f_0) + \sinc^2_{π}(f + f_0) + 2\sinc_{π}(f - f_0)\sinc_{π}(f + f_0)\cos(2β)) \end{equation}

waarin \( \sinc_{π}(x) \) een functie is gelijk aan \( \frac{\sin(π x)}{π x} \) als \( x \) ongelijk is aan 0, en 1 als \( x \) gelijk is aan 0, en \( β \) de fase van de harmonische component is aan het begin van de waarnemingen. De sinc-functie heeft de grootste waarde als haar argument gelijk is aan 0, dus zijn de termen met \( f - f_0 \) erin het grootste in de buurt van \( f = f_0 \) en de termen met \( f + f_0 \) erin het grootste in de buurt van \( f = -f_0 \). Er zijn dus twee even hoge pieken, symmetrisch rondom frequentie 0: één piek bij een positieve frequentie, en één piek bij een negatieve frequentie. De termen van de negatieve piek (de piek bij negatieve frequentie) zijn bij positieve frequenties kleiner dan die van de positieve piek. Als je de termen van de negatieve piek (met \( f + f_0 \) erin) negeert (die voor \( f \gg 0 \) weinig invloed hebben), dan volgt uit de vorige vergelijking dat

\begin{align} P_f & = A^2 \sinc^2_{π}(f - f_0) \\ P_{f+1} & = P_f \frac{f^2}{(f + 1)^2} \\ P_{f-1} & = P_f \frac{f^2}{(f - 1)^2} \end{align}

Dit hangt nu niet meer af van \( β \). Als \( f_x \) de hele frequentie is waarbij in het frequentiespectrum de hoogste piek zit, en \( f_0 = f_x + k \), dan is \( k \) tussen −1 en +1. Met

\begin{equation} q = \frac{P_{f_x}^2}{P_{f_{x+1}}P_{f_{x-1}}} \end{equation}

volgt uit de twee voorgaande genummerde vergelijkingen dan

\begin{align} k & = \frac{\sgn(P_{f_{x+1}} - P_{f_{x-1}})}{\sqrt{1 - \sqrt{q}}} \label{eq:k1} \\ A = \frac{\sqrt{P_{f_x}}}{\sinc_{π}(k)} \label{eq:A1} \end{align}

Hiermee kun je dan de frequentie en amplitude van alle pieken uit het frequentiespectrum schatten. De procedure is dan als volgt:

1. Vind de hele frequentie \( f_0 \) waar het frequentiespectrum meer energie vertoont dan voor de buurfrequenties.
2. Gebruik vergelijking \ref{eq:k1} om de frequentieafwijking \( k \) te schatten. De frequentie \( f \) van de component is dan geschat op \( f_0 + k \). (Als \( |k| \) groter is dan 1, dan kloppen de aannames waarop de huidige methode gebaseerd is niet voor deze meetreeks. Zie dan verderop.)
3. Gebruik vergelijking \ref{eq:A1} om de amplitude van de component te schatten.

Als we dit toepassen op de meetreeks van hierboven dat vinden we

# 

De amplitude van de grootste piek wordt nu een stuk hoger geschat, en ook die van de piek bij frequentie 33.

3.3. Nog betere schatting: Methode 3

De rare vorm (in de vorige grafiek) van de piek bij frequentie 78 en de bijna-onzichtbaarheid van de piek bij frequentie 33 komen natuurlijk van de brede "hellingen" die van de hoogste piek komen, die de lagere pieken verstoren. We moeten dus proberen om de hellingen van elke piek steiler te laten lopen zodat ze andere pieken niet zoveel vervormen. Dat doen we door de meetreeks te vermenigvuldigen met een masker dat de randen van de meetreeks dempt. Het masker dat we hier zullen gebruiken is gelijk aan een cosinus die waarden tussen 0 en 2 heeft, precies eenmaal in de meetperiode past, en aan de randen de waarde 0 heeft. (Dit is een vorm van het hanningmasker of de hanningfunctie.) Bij toepassing van zo'n masker is het heel belangrijk dat het gemiddelde van de meetreeks gelijk is aan nul.

Fig. 4: Frequentiespectrum

Merk op dat de vertikale as van deze grafiek een veel groter bereik heeft dan in de vorige grafiek. In het vorige spectrum was de energieverhouding tussen de bovenkant en de onderkant van de grafiek ongeveer een miljard (10⁹), maar in deze grafiek is het nog eens tien miljard maal zo groot (in totaal ongeveer 10¹⁹). De pieken zijn iets dikker geworden, maar de hellingen lopen een stuk steiler en daarom valt de kleinere piek hier veel beter op dan in de vorige grafiek.

Op dezelfde wijze als voorheen kun je het fourierspectrum van een enkele harmonische component met frequentie \( f_0 \) en amplitude \( A \) uitrekenen. Net als voorheen negeren we de bijdragen van de negatieve piek. We vinden dan

\begin{equation} P_f = A^2 \left( \frac{\sinc_{π}(f - f_0)}{(f - f_0)^2 - 1} \right)^2 \end{equation}

Hieruit volgt (op vergelijkbare wijze als voorheen) de schatting

\begin{align} q & = \frac{P_{f_x}}{\sqrt{P_{f_{x+1}} P_{f_{x-1}}}} \\ k & = \sqrt{\frac{q - 4}{q - 1}} \sgn(P_{f_{x+1}} - P_{f_{x-1}}) \label{eq:k3} \\ A & = (1 - k^2) \frac{\sqrt{P_{f_x}}}{\sinc_{π}(k)} \end{align}

(Als \( |k| \) groter is dan 1, dan kloppen de aannames waarop de huidige methode gebaseerd is niet voor deze meetreeks. Zie dan verderop.)

Toegepast op de meetreeks geeft dat

# 

Nu wordt de piek bij frequentie 78 hoger ingeschat dan die bij frequentie 32.

3.4. Methode 4: Algemeen Kwadratisch

De twee voorgaande methoden gaan er van uit dat de componenten in het fourierspectrum puur harmonisch zijn, dus met een enkele, welbepaalde frequentie en amplitude, maar dit komt in de praktijk niet zo vaak voor. Meestal variëren de frequentie en amplitude van alle componenten een beetje met de tijd, waardoor de pieken in het fourierspectrum dikker worden. Als de variatie te groot wordt dan geven die methoden misschien wel helemaal geen schatting voor hoeveel de "echte" frequentie van de fourierpiekfrequentie afwijkt. In zulke gevallen kun je de eenvoudigere aanname doen dat het verband tussen energie en frequentie bij de pieken in het fourierspectrum kwadratisch is. Dan vinden we (op vergelijkbare wijze als voorheen)

\begin{align} k & = -\frac{1}{2} \frac{P_{f_{x+1}} - P_{f_{x-1}}}{P_{f_{x+1}} + P_{f_{x-1}} - 2P_{f_x}} \\ A & = \sqrt{P_{f_x} + \frac{1}{4} (P_{f_{x+1}} - P_{f_{x-1}}) k} \end{align}

Je kunt dan ook weer een masker toepassen om de hellingen van de pieken steiler te maken. Toegepast op onze meetreeks geeft deze methode

# 

3.5. Opmerkingen

1. Het onderscheidend vermogen van een frequentiespectrum is beperkt. Als de meetreeks een tijd \( T \) duurt, dan zijn de frequenties veelvouden van \( 1/T \). Frequentieverschillen kleiner dan \( 1/T \) kunnen uit zo'n frequentiespectrum niet gemeten worden. Twee harmonische componenten waarvan de frequenties minder van elkaar verschillen dan dat kun je in het frequentiespectrum zeker niet uit elkaar houden.
2. Fig. 5: Frequentiespectrum met aliassing
   
   Een harmonische component die een frequentie heeft die hoger is dan de grootste frequentie die in het frequentiespectrum staat komt wel in het spectrum terecht, maar op een frequentie die kleiner is dan zijn eigen frequentie. Als de meetreeks een tijd \( T \) duurt en \( N \) meetwaarden omvat, dan is de hoogste frequentie in het spectrum ongeveer gelijk aan \( \frac{N}{2T} \). Als de frequentie \( f \) van een harmonische component groter is dan die grootste frequentie, dan komt die component in het spectrum terecht bij frequentie \( \left| f - N \left[ \frac{f}{N} \right] \right| \). Een voorbeeld is te zien in figuur 5, die drie verschillende harmonische componenten toont die elkaar snijden in de meetpunten (aangegeven door vierkantjes). Uit de metingen op die meetpunten kun je dus niet afleiden welke component het nu precies was.
3. Als de meetreeks ook een component bevat die niet periodiek is (of misschien periodiek met een periode die veel groter is dan de lengte van de meetreeks), dan komt de energie van die component in het frequentiespectrum terecht bij de kleinste frequenties, waar het eventuele pieken verstoort en onterecht een hoge amplitude bezorgt. Om dit tegen te gaan kun je een lopend gemiddelde van de meetreeks aftrekken, maar daarmee onderdruk je niet alleen niet-periodieke componenten maar ook wel-periodieke componenten met perioden die langer zijn dan de breedte van het lopende gemiddelde.
4. Methoden 2 en 3 negeren de invloed van negatieve pieken, die het sterkst is bij de kleinste frequenties (nog steeds kleiner dan de invloed van de positieve pieken, maar niet heel veel kleiner). Daarom zullen de resultaten bij de kleinste frequenties minder nauwkeurig zijn dan bij grotere frequenties.
5. Het steiler laten lopen van de hellingen van de pieken door methode 3 heeft ook een prijs, namelijk dat de pieken aan de bovenkant wat breder worden dan bij methode 2. Als twee pieken heel dicht bij elkaar staan, dan kan methode 3 meer moeite hebben om ze uit elkaar te houden dan methode 2. Het lijkt erop dat methode 2 soms nauwkeuriger is dan methode 3 om de grootste component te vinden, maar methode 3 is bijna altijd beter om kleinere componenten te vinden.

4. De nauwkeurigheid

4.1. Testen met de proefmeetreeks

De harmonische componenten waaruit de meetreeks gebouwd werd zijn als volgt:

# 

Nu kunnen we bepalen hoe nauwkeurig de verschillende schattingen waren. De volgende tabel toont het verschil tussen de schatting en de echte waarde voor de frequentie (\( ∆f_0 \)) en het relatieve verschil voor de amplitude (\( \frac{∆A}{A} \)).

# 

Methode 3 lijkt in het algemeen het beste (als de meetserie alleen zuivere harmonische componenten bevat). Dit wordt bevestigd door een groot aantal tests met diverse willekeurige samengestelde meetseries (met een exponentiële verdeling van amplituden). De nauwkeurigheid van de frequenties lijkt voor methode 1 ongeveer constant te zijn, en voor methoden 2 en 3 toe te nemen met ongeveer de 2,7-de macht van de gemiddelde ruimte (in het spectrum) tussen de harmonische pieken, maar wel met een gemiddeld voordeel van ongeveer een factor 10 voor methode 3 ten opzichte van methode 2. Methode 4 is een stuk beter dan methode 1, maar vaak slechter dan methoden 2 en 3 als de meetserie alleen zuivere harmonische componenten bevat. Maar als de componenten in de meetserie niet zuiver harmonisch zijn, dan werken methoden 2 en 3 soms helemaal niet, en dan is methode 4 het beste.

De relatieve nauwkeurigheid van de amplituden lijkt voor methode 1 ongeveer constant te zijn, voor methode 2 toe te nemen met ongeveer het kwadraat van de gemiddelde ruimte (in het spectrum) tussen de harmonische pieken, en voor methode 3 met ongeveer de derde macht van de gemiddelde ruimte tussen de harmonische pieken.

4.2. Testen met het VSOP-model

Hieronder testen we de methoden op het VSOP-model van de beweging van de planeten. Uit dit model kan men de posities van de planeten uitrekenen voor een paar duizend jaar in de toekomst en in het verleden. Het model is opgebouwd uit een groot aantal harmonische componenten, dus we kunnen testen hoe goed we die componenten weer uit de posities van de planeten kunnen halen.

Hieronder staat een tabel die de tien grootste harmonische componenten uit de eclipticale heliocentrische lengtegraad \( λ \) van de Aarde ten opzichte van de equinox en ecliptica van de datum volgens het VSOP-model laat zien (\( f_\text{vsop} \) is de frequentie in radialen per juliaans millennium van 365250 dagen, \( A_\text{vsop} \) is de amplitude in 10⁻⁸ radialen), zoals gegeven door [Meeus], en ook de tien grootste harmonische componenten uitgerekend uit \( λ \) (berekend voor elke 10 dagen in een periode van 10.000 jaar rond het jaar 2000) volgens methoden 1−4 (\( f \) en \( A \), in de zelfde eenheden als voor het VSOP-model). Om niet-periodieke componenten bij kleine frequenties te onderdrukken verwijderde ik eerst een lopend gemiddelde over 68 jaar uit de meetreeksen voor dat ik er harmonische componenten in probeerde te vinden.

 # 

Tussen de 10 grootste componenten gevonden door methode 1 staan 6 van de 10 grootste componenten uit het VSOP-lijstje. Voor methode 2 is dat 5, en voor methoden 3 en 4 is dat 9. De volgende twee tabellen tonen de absolute (voor de frequentie \( f \)) of relatieve (voor de amplitude \( A \)) fout voor de eerste drie VSOP-componenten (\( ∆_1 \), \( ∆_2 \), \( ∆_3 \)) en ook het minimum (\( \min \)), gemiddelde (\( μ \)), maximum (\( \max \)) en standaardafwijking (\( σ \)) van diezelfde soort fouten voor de eerste tien VSOP-componenten, voor methoden 1−4 (M1 - M4). Om zo'n fout te bepalen zocht ik voor elk van de VSOP-componenten naar de met methode 1−4 geschatte component die daar het dichtste bij lag.

Bijvoorbeeld, methode 2 schat de frequentie van de grootste VSOP-component 0,0565 eenheden te klein en de amplitude een fractie 0,0044 (dus 0,44%) te klein. Methode 3 doet het meestal het beste (met de kleinste \( σ \) die toch flink groter is dan de \( μ \)) voor zowel de frequentie als de amplitude.

Hieronder staat een vergelijkbare tabel voor de eclipticale heliocentrische breedtegraad \( β \) van de Aarde. Meeus geeft voor dit geval maar 5 componenten, maar in de berekeningen van de posities heb ik ook alle kleinere componenten meegenomen die Meeus niet noemt.

 # 

Alle vier de methoden vinden de vijf grootste VSOP-componenten, en ook nog in dezelfde volgorde.

Ook hier doet methode 3 het meestal het beste, en methode 2 is bijna even goed. Methoden 1 en 4 onderschatten de amplitude systematisch.

Hieronder staat een vergelijkbare tabel voor de afstand \( r \) tussen de Aarde en de Zon, gemeten in eenheden van 10⁻⁸ AE.

 # 

In de tien grootste componenten gevonden door methode 1 staan 5 van de grootste tien VSOP-componenten. Voor methode 2 is dat 3, en voor methoden 3 en 4 is het 10 (en nog allemaal in de goede volgorde ook).

Ook hier is methode 3 meestal het beste. Voor de frequentie volgt methode 4, maar voor de amplitude methode 2. Voor de frequentie is methode 2 hier meestal nog slechter dan methode 1. Methoden 1 en 4 onderschatten de amplitude systematisch.

talen: [en] [nl]

http://aa.quae.nl/nl/reken/frequentie.html;
Laatst vernieuwd: 2016−02−07