AstronomieAntwoorden
Frequentiebepaling

\(\def\|{&}\DeclareMathOperator{\D}{\bigtriangleup\!} \DeclareMathOperator{\d}{\text{d}\!}\)

\( \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. Uit metingen op regelmatige tijden

Fig. 1: Meetreeksgrafiek

Fig. 2: Meetreeksgrafiek met proefreeks

1.1. 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.

1.2. Eenvoudige schattingen

1.2.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.

1.2.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}

\begin{equation} \sinc_{π}(x) ≡ \left\{ \begin{split} \dfrac{\sin(πx)}{πx} \| \qquad (x ≠ 0) \\ 1 \| \qquad (x = 0) \end{split} \right. \end{equation}

waarin \( β \) de fase is 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.

1.2.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.

1.2.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

# 

1.2.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.

1.3. De nauwkeurigheid

1.3.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.

1.3.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.

2. Uit tijdstippen van gebeurtenissen

Stel, je hebt een reeks van \( n \) oplopende tijdstippen \( t_i \) (zonder dubbele), waarop een interessante gebeurtenis voorkwam. Hoe kun je perioden \( P \) in die reeks vinden?

Hieronder staat een methode beschreven die ik bedacht om dit probleem op te lossen. Het blijkt beter te passen om die beschrijving te doen aan de hand van frequentie \( ν \) in plaats van periode \( P \). Ze zijn verbonden door \( νP = 1 \), dus \( ν = 1/P \) en \( P = 1/ν \).

2.1. De Methode

1. Bereken de tijdsduur \( T \) van de meetreeks, het verschil tussen de laatste en eerste meettijd:

   \begin{equation} T = t_{\text{max}} − t_{\text{min}} \end{equation}

2. Bepaal het kleinste verschil \( \D t_{\text{min}} \) tussen de opvolgende tijdstippen:

   \begin{equation} \D t_{\text{min}} = \min_i(t_i − t_{i−1}) \end{equation}

3. Doe nu voor alle frequenties \( f_m \) tussen \( 2/T \) en \( 1/\D t_{\text{min}} \) in stappen van \( 1/(2T) \) het volgende:

   1. Bereken de fase \( φ_i \) van elk tijdstip \( t_i \):

      \begin{equation} φ_i = f_mt_i \bmod 1 = f_mt_i − ⌊f_mt_i⌋ \end{equation}

      Hier betekent \( \bmod 1 \) dat je de rest bepaalt bij deling door 1, en betekent \( ⌊\cdot⌋ \) dat je afrond naar beneden naar het dichtstbijzijnde gehele getal. Er geldt \( 0 \le φ_i \lt 1 \).

   2. Sorteer de fases in oplopende volgorde. Die gesorteerde fases noemen we \( φ'_k \), met \( k \) lopend van 1 tot \( n \).

   3. De rest volgt later...

4. Als de hoogste waarde van \( s \) duidelijk hoger is dan \( 1/n \) dan is dat een goede aanwijzing dat de bijbehorende frequentie \( ν_m \) in de buurt ligt van een frequentie \( ν_0 \) die bij de metingen past.

   Als \( s_m = 1 \) dan is \( ν_0 = ν_m \).

   Als \( s_m \lt 1 \) dan kan dat zijn omdat \( ν_m \) niet precies gelijk is aan \( ν_0 \) (omdat \( ν_0 \) niet in het rijtje uitgeprobeerde frequenties zit), of omdat er tijdstippen van meer dan één periode in de meetreeks zitten, of omdat de periode van het verschijnsel een beetje met de tijd varieert, of omdat er metingen ontbreken, of omdat er ruis op de metingen zit.

   Je kunt nu het frequentieinterval van \( ν_{m−1} \) tot \( ν_{m+1} \) verdelen in (bijvoorbeeld) 10 gelijke delen en dan bovenstaand recept herhalen voor die nieuwe lijst van frequenties en op die manier steeds dichter bij \( ν_0 \) te komen.

   Ik ben nog op zoek naar manieren om sneller bij \( ν_0 \) te komen die ook werken als de meetreeks op diverse manieren niet perfect is.

2.2. Voorbeelden

2.3. De Uitleg

Stel, je vermoedt dat alleen frequentie \( ν \) voorkomt in deze meetreeks, dus dat alle waarnemingen voldoen aan \( t_i = t_0 + iP = t_i = t_0 + i/ν \) dus \( νt_i = νt_0 + i \). Dan is \( νt_i = νt_0 \pmod{1} \) voor elke meting. Dit is de basis voor onze methode.

We definiëren de fase \( φ \):

\begin{equation} φ = νt \bmod 1 \end{equation}

De fase \( φ \) is tussen 0 en net onder 1, voor elke frequentie. Bereken de fase \( φ \) voor alle tijdstippen \( t_i \).

Fig. 6: Fase voor verschillende frequenties, één frequentie in de data

Figuur 6 toont de fase \( φ \) als een functie van tijd \( t \) voor onze voorbeeldmeetreeks met alleen \( ν_1 \). Voor frequentie 0.087 is de fase steeds hetzelfde. Voor de andere frequenties is de fase steeds anders.

---

Fig. 7: Fase voor één frequentie, twee frequenties in de data

Figuur 7 toont de fase \( φ \) als functie van tijd \( t \) voor de meetreeks met beide frequenties erin, met \( ν_1 \) als de testfrequentie waarmee de fase berekend is. We zien veel meetpunten met dezelfde fase 0.6351 − die horen bij \( ν_1 \). De andere meetpunten horen bij \( ν_2 \) en tonen fases die ongeveer gelijk verdeeld zijn tussen 0 en 1.

---

Sorteer de fase \( φ \) in oplopende volgorde. \( φ_k \) is de \( k \)-de fase in oplopende volgorde van fase (met \( k \) lopend van 0 tot en met \( n − 1 \)). Onthoud de \( i \) voor elke \( k \). Voor het gemak definiëren we

\begin{equation} x = \dfrac{k}{n} \end{equation}

dus \( x \) geeft aan welke fractie van de gesorteerde meetpunten we al gehad hebben. Dan lopen zowel \( φ \) als \( x \) van 0 tot net onder 1, ongeacht het aantal meetpunten.

Fig. 8: Gesorteerde fase voor verschillende perioden

Figuur 8 toont de gesorteerde fase \( φ \), voor dezelfde frequenties als in figuur 6. Voor frequentie 0.087 die echt in de meetreeks voorkomt is er een horizontaal plateau op fase 0.6351 − dat zijn de meetpunten die bij die frequentie horen. Voor frequentie 0.200 die ver van de juiste frequentie ligt zijn de fases redelijk gelijkmatig verdeeld tussen 0 en 1 en is de grafiek ongeveer een rechte lijn van linksonder naar rechtsboven.

---

Als er meetpunten bij frequentie \( ν_0 \) horen, dan kun je met op de fase gebaseerde detectiemethoden ook veelvouden en omgekeerde veelvouden van die frequentie detecteren. Als er meetpunten bij frequentie \( ν_0 \) horen dan is de bijbehorende fase \( ν_0t \bmod 1 = φ_ν \) constant voor die meetpunten. De fase voor een veelvoud \( mν_0 \) van die frequentie (met \( m \) een heel getal groter dan 0) is dan

\begin{equation} mν_0t_i = mν_0 \left( t_0 + \dfrac{i}{ν_0} \right) = mφ_ν + mi = mφ_ν \pmod{1} \end{equation}

dus ook dan geven al die meetpunten een constante fase, die een veelvoud (mod 1) is van \( φ_ν \). De fase van een omgekeerd veelvoud \( ν_0/m \) van de beginfrequentie is

\begin{equation} \dfrac{ν_0t_i}{m} = \dfrac{ν_0}{m} \left( t_0 + \dfrac{i}{ν_0} \right) = \dfrac{φ_ν + i}{m} \pmod{1} \end{equation}

dus zijn de meetpunten dan verdeeld over \( m \) klassen met in elke klasse een vaste fase die een veelvoud van \( 1/m \) verschoven is ten opzichte van \( φ_ν \).

Fig. 9: Gesorteerde fase voor (omgekeerde) veelvouden van de frequentie

Figuur 9 toont de gesorteerde fase \( φ \) voor veelvouden en omgekeerde veelvouden van frequentie \( ν = ν_0 = 0.087 \) waarbij een deel van de meetpunten horen.

Het horizontale plateau bij \( φ = φ_ν = 0.635 \) wordt gevormd door de meetpunten die bij \( ν_0 \) horen.

Voor een veelvoud van die basisfrequentie \( ν_0 \) verschuift het plateau naar hetzelfde veelvoud (modulus 1) van \( φ_ν \). Dus voor \( 2ν_0 \) verschuift het plateau naar \( 2φ_ν \bmod 1 = 2×0.635 \bmod 1 = 0.270 \), en voor \( 3ν_0 \) naar \( 3φ_ν \bmod 1 = 0.905 \).

Voor een omgekeerd veelvoud van \( ν_0 \) (dus \( ν_0 \) gedeeld door een heel getal) wordt het plateau opgesplitst in kortere plateau's, elk bij een andere fase. Voor frequentie 0.087/2 verschuift de helft van het plateau naar fase \( φ_ν/2 = 0.318 \) en de andere helft naar \( (φ_ν + 1)/2 = 0.818 \). Voor frequentie 0.087/3 wordt het plateau net zo verdeeld in drieën.

---

\( \D φ \) is het verschil tussen opeenvolgende gesorteerde fasen:

\begin{equation} \D φ_k = \left\{ \begin{split} φ_k − φ_{k−1} \| \qquad (1 \le k \lt n) \\ φ_0 + 1 − φ_{n−1} \| \qquad (k = 0) \end{split} \right. \end{equation}

De reden voor de aanwezigheid van \( \D φ_0 \) is dat zonder die regel de resultaten af kunnen hangen van de keuze van het nulpunt van de tijd (\( t = 0 \)), en dat is ongewenst. Nu is de som van de gesorteerdefaseverschillen altijd gelijk aan 1.

Het zou rekentechnisch fijn zijn als we als maat voor de kwaliteit van de periode konden gebruiken:

\begin{equation} s_1 = \sum_{k=0}^{n−1} f(\D φ_k) \end{equation}

waarin \( f(\D φ) \) een nog nader te bepalen functie is. Helaas blijkt dit in de praktijk alleen te kunnen werken voor meetreeksen waarin nagenoeg alle meetpunten bij één periode horen. We moeten dus niet alleen de losse faseverschillen gebruiken, maar ook naar hun volgorde kijken, of de fases zelf gebruiken.

Bekijk figuur 8 nog eens. Voor een verkeerde periode is de grafiek ongeveer een rechte lijn van linksonder naar rechtsboven. Die lijn noemen we hier de diagonaal. Voor de goede periode komen sommige punten van de grafiek ver van de diagonaal, en dat kunnen we gebruiken als maat voor de kwaliteit van de periode. Het is daarbij een probleem dat hoever de meetpunten voor de goede periode van de diagonaal zijn afhangt van het gekozen nulpunt van de tijd.

Als we een kwaliteitsmaat willen baseren op de afstand van de meetpunten van de diagonaal, dan is één manier om met die afhankelijkheid om te gaan om de kwaliteitsmaat uit te rekenen voor

\begin{equation} φ' = (φ + c) \bmod 1 \end{equation}

voor veel faseverschuivingen \( c \) regelmatig verdeeld tussen 0 en 1, en dan de uitkomsten op een toepasselijke manier samen te voegen. In de praktijk doen kwaliteitsmaten gebaseerd op de afstand tot de diagonaal het ook niet zo goed.

De beste kwaliteitsmaat die ik heb gevonden is gebaseerd op het ontdekken van opeenvolgende lage \( \D φ \), dus van de plateau's in grafieken van gesorteerde \( φ \) versus \( x \), zoals figuren 8 en 9.

Stel, de frequentie \( ν \) die we proberen ligt heel dicht bij frequentie \( ν_0 \) die in de data voorkomt, dus

\begin{equation} \D ν = ν − ν_0 \end{equation}

is klein. \( m \) van de \( n \) meetpunten horen bij frequentie \( ν_0 \), en de rest niet. Dan zal er in een grafiek zoals figuur 8 een plateau zijn dat een fractie \( m/n \) breed is in \( x \). De tijdstippen van de \( m \) meetpunten zijn

\begin{equation} t_i = t_0 + \dfrac{i}{ν_0} \end{equation}

voor hele getallen \( i \), en de overeenkomstige fases uitgerekend met frequentie \( ν \) zijn

\begin{equation} φ_i = νt_i = νt_0 + \dfrac{iν}{ν_0} = φ_0 + \dfrac{i\D ν}{ν_0} \pmod{1} \end{equation}

Als de \( i \) varieert van \( i_{\text{min}} \) tot \( i_{\text{max}} \) (met \( i_{\text{max}} − i_{\text{min}} + 1 \ge m \)) dan is het gesorteerde faseverloop

\begin{equation} Φ = |φ_{i_{\text{max}}} − φ_{i_{\text{min}}}| = (i_{\text{max}} − i_{\text{min}}) \dfrac{|\D ν|}{ν_0} ≈ T \left| \D ν \right|\end{equation}

Alleen als \( |Φ| \ll 1 \) kunnen we detecteren dat we in de buurt van de juiste frequentie zijn. We kiezen \( Φ_0 = 1/10 \ll 1 \) als toegestane bovengrens van \( Φ \). We willen de frequentie ook kunnen detecteren als alle \( n \) meetpunten bij de frequentie horen, dus de gemiddelde fasestap op het plateau moet voldoen aan \( \D φ \le Φ_0/n \).

1. Begin bij \( k = −1 \).
2. Zoek de eerstvolgende \( k \) waarvoor \( \D φ_k \le Φ_0/n \). Als die er niet is (vóór \( k = n \)) dan ben je klaar. Als die er wel is dan heb je het begin van een kandidaat-plateau gevonden. Die \( k \) noemen we \( k_1 \).
3. Ga verder zolang de gemiddelde fasestap op het tot dan toe gevonden kandidaat-plateau niet groter is dan \( Φ_0/n \). Als je bij \( k = n − 1 \) komt en ook daarvoor aan die eis is voldaan, ga dan verder met \( k = 0 \). De laatste \( k \) waarvoor aan de eis is voldaan noemen we \( k_2 \). Als \( k_2 \lt k_1 \) en er is overlap met het eerstgevonden plateau, gooi dan dat eerstgevonden plateau weg.
4. Als \( k_2 = k_1 \) dan is er maar 1 faseverschil in het kandidaat-plateau, en dus maar hooguit 2 bijbehorende meetpunten. Dan wijzen we het kandidaat-plateau af en gaan terug naar punt 2.
5. De kwaliteitsmaat \( s_2 \) voor de frequentie \( ν \) is

   \begin{equation} s_2 = (k_2 − k_1 + 1) \bmod n \end{equation}

   Met deze kwaliteitsmaat kun je zien dat je in de buurt bent van een gezochte frequentie, en ongeveer hoeveel meetpunten er bij die frequentie horen, maar niet zo goed hoe ver je van de gezochte frequentie bent.

6. Bekijk de verzameling van fasestappen van \( \D φ_{k_1} \) tot en met \( \D φ_{k_2} \). Als ze allemaal bij dezelfde frequentie horen dan zouden ze allemaal veelvouden van dezelfde waarde moeten zijn (namelijk van \( |\D ν|/ν_0 \)). Bepaal de deelverzamelingen die allemaal (praktisch) veelvouden zijn van dezelfde waarde. Als \( \D φ_a \bmod \D φ_b \lt \min(\D φ_a, \D φ_b)/100 \) dan vinden we dat \( \D φ_a \) en \( \D φ_b \) bij dezelfde deelverzameling horen. Alle deelverzamelingen met minder dan 3 leden wijzen we af. Als er geen deelverzamelingen overblijven, ga dan terug naar punt 2.
7. Bereken voor elke overblijvende deelverzameling het gemiddelde \( μ \) van de fasestappen \( \D φ \) in die deelverzameling.

   Bekijk voor de eerste twee \( k \) in elke deelverzameling de bijbehorende \( i \). Als de \( i \) afnemen, dan nam de fase af met de tijd. Vermenigvuldig dan het gemiddelde \( μ \) met −1.

   Een schatting voor de gezochte frequentie \( ν_0 \) is dan \( ν/(1 + μ) \).

   Het aantal leden van de deelverzameling is een schatting voor hoeveel meetpunten er bij de gezochte frequentie horen.

   Als er meerdere deelverzamelingen naar dezelfde \( ν_0 \) wijzen dan tellen de meetpunten op.

talen: [en] [nl]

//aa.quae.nl/nl/reken/frequentie.html;
Laatst vernieuwd: 2021-07-19