\(\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.
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.
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.
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
# | \({f_0}\) | \({A}\) |
---|---|---|
1 | 10 | 0,915 |
2 | 32 | 0,0065 |
3 | 78 | 0,0024 |
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.
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:
Als we dit toepassen op de meetreeks van hierboven dat vinden we
# | \({f_0}\) | \({A}\) |
---|---|---|
1 | 10,30033 | 1,066 |
2 | 32,66 | 0,016 |
3 | 78,26 | 0,0027 |
De amplitude van de grootste piek wordt nu een stuk hoger geschat, en ook die van de piek bij frequentie 33.
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.
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 (109), maar in deze grafiek is het nog eens tien miljard maal zo groot (in totaal ongeveer 1019). 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
# | \({f_0}\) | \({A}\) |
---|---|---|
1 | 10,29955 | 1,04978 |
2 | 77,6024 | 0,0019970 |
3 | 31,60 | 0,0012 |
Nu wordt de piek bij frequentie 78 hoger ingeschat dan die bij frequentie 32.
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
# | \({f_0}\) | \({A}\) |
---|---|---|
1 | 10,18609 | 1,001697 |
2 | 77,6994 | 0,0018485 |
3 | 31,25 | 0,00095 |
De harmonische componenten waaruit de meetreeks gebouwd werd zijn als volgt:
# | \({f_0}\) | \({A}\) | \({β}\) |
---|---|---|---|
1 | 10,3 | 1,05 | 36° |
2 | 77,6 | 0,002 | 0° |
3 | 31,4 | 0,001 | 0° |
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} \)).
\({Δf_0}\) | \({\frac{ΔA}{A}}\) | |||||||
---|---|---|---|---|---|---|---|---|
# | Methode 1 | Methode 2 | Methode 3 | Methode 4 | Methode 1 | Methode 2 | Methode 3 | Methode 4 |
1 | −0,3 | +0,00033 | −0,00045 | −0,11 | −0,12 | +0,015 | −0,00021 | −0,046 |
2 | +0,4 | +0,66 | +0,0024 | +0,099 | +2,3 | +0,36 | −0,0015 | −0,076 |
3 | +0,6 | +1,26 | +0,20 | −0,15 | +1,4 | +14,7 | +0,18 | −0,050 |
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.
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−8 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.
\({λ}\) | ||||||||||
---|---|---|---|---|---|---|---|---|---|---|
# | \({f_\text{vsop}}\) | \({A_\text{vsop}}\) | \({f_1}\) | \({A_1}\) | \({f_2}\) | \({A_2}\) | \({f_3}\) | \({A_3}\) | \({f_4}\) | \({A_4}\) |
1 | 6283,0757 | 3341656 | 6283,1853 | 2957169 | 6283,0191 | 3326893 | 6283,0366 | 3304973 | 6283,0874 | 3213528 |
2 | 12566,1572 | 34894 | 12565,7423 | 23781,5 | 12566,0368 | 35193,11 | 12566,0147 | 34050,8 | 12566,0137 | 31518,94 |
3 | 5753,3848 | 3497 | 5753,5128 | 3159,722 | 6069,1218 | 6556,893 | 5753,3852 | 3496,726 | 5753,4425 | 3419,802 |
4 | 3,5231 | 3418 | 77713,5775 | 2666,237 | 5753,2840 | 3970,837 | 77713,7709 | 3134,605 | 77713,6994 | 2983,824 |
5 | 77713,7734 | 3136 | 7860,2648 | 2597,325 | 5880,2491 | 3630,401 | 7860,4192 | 2676,357 | 7860,3537 | 2591,483 |
6 | 7860,4194 | 2676 | 6069,5570 | 2477,555 | 5855,4863 | 3434,766 | 3930,2098 | 2342,756 | 3930,1724 | 2323,368 |
7 | 3930,2097 | 2343 | 3930,1324 | 2328,372 | 77713,7717 | 3136,410 | 11506,7698 | 1324,200 | 529,6594 | 1270,505 |
8 | 11506,7695 | 1324 | 5885,4597 | 1470,287 | 6681,9815 | 3059,166 | 529,6725 | 1270,298 | 11506,8305 | 1221,383 |
9 | 529,6910 | 1273 | 5223,8403 | 1312,458 | 7860,1163 | 2852,389 | 1577,4341 | 1196,691 | 1577,4788 | 1095,007 |
10 | 1577,3435 | 1199 | 5879,8048 | 1300,653 | 5885,8305 | 2839,734 | 5884,9271 | 990,4211 | 5884,8816 | 977,9457 |
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.
\({f}\) | |||||||
---|---|---|---|---|---|---|---|
\({Δ_1}\) | \({Δ_2}\) | \({Δ_3}\) | \({\min}\) | \({μ}\) | \({\max}\) | \({σ}\) | |
M1 | 0,1096 | −0,4149 | 0,1280 | −0,4149 | −0,0385 | 0,2560 | 0,2237 |
M2 | −0,0565 | −0,1204 | −0,1008 | −0,3031 | 0,0073 | 0,6496 | 0,2492 |
M3 | −0,0391 | −0,1425 | 0,0004 | −0,1433 | −0,0255 | 0,0906 | 0,0703 |
M4 | 0,0117 | −0,1436 | 0,0577 | −0,1695 | −0,0256 | 0,1353 | 0,0946 |
\({A}\) | |||||||
---|---|---|---|---|---|---|---|
\({Δ_1}\) | \({Δ_2}\) | \({Δ_3}\) | \({\min}\) | \({μ}\) | \({\max}\) | \({σ}\) | |
M1 | −0,1151 | −0,3185 | −0,0964 | −0,3185 | −0,0557 | 0,0379 | 0,1058 |
M2 | −0,0044 | 0,0086 | −0,0178 | −0,0178 | 0,0087 | 0,0612 | 0,0230 |
M3 | −0,0110 | −0,0242 | −0,0001 | −0,0242 | −0,0016 | 0,0230 | 0,0116 |
M4 | −0,0383 | −0,0967 | −0,0221 | −0,0967 | −0,0269 | 0,0187 | 0,0325 |
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.
\({β}\) | ||||||||||
---|---|---|---|---|---|---|---|---|---|---|
# | \({f_\text{vsop}}\) | \({A_\text{vsop}}\) | \({f_1}\) | \({A_1}\) | \({f_2}\) | \({A_2}\) | \({f_3}\) | \({A_3}\) | \({f_4}\) | \({A_4}\) |
1 | 84334,6641 | 280 | 84334,7981 | 258,1 | 84334,6603 | 279,8 | 84334,6646 | 279,2 | 84334,7219 | 272,5 |
2 | 5507,5532 | 102 | 5507,2119 | 76,91 | 5507,4655 | 102,2 | 5507,4691 | 102,2 | 5507,4022 | 94,04 |
3 | 5223,6899 | 80 | 5223,8403 | 66,76 | 5223,6270 | 81,32 | 5223,6410 | 79,94 | 5223,7014 | 76,04 |
4 | 2352,8701 | 44 | 2353,0529 | 42,58 | 2352,9514 | 44,46 | 2352,9227 | 44,55 | 2353,0000 | 43,44 |
5 | 1577,3400 | 32 | 1577,0795 | 28,1 | 1577,2574 | 32,18 | 1577,2653 | 32,21 | 1577,1857 | 30,72 |
6 | 1048,0353 | 15,54 | 1047,7396 | 23,08 | 1047,7791 | 22,03 | 1047,7835 | 20,58 | ||
7 | 5856,5570 | 15,01 | 6283,4136 | 17,84 | 5856,3454 | 17,06 | 6283,2864 | 16,03 | ||
8 | 6283,1853 | 14,21 | 5856,3852 | 17,03 | 6283,1853 | 15,90 | 5856,4586 | 15,98 | ||
9 | 10213,3177 | 13,88 | 9437,6669 | 15,25 | 9437,6385 | 15,20 | 10213,3463 | 14,14 | ||
10 | 14143,4501 | 10,36 | 10213,3738 | 14,06 | 10213,3177 | 14,13 | 9437,7010 | 13,21 |
Alle vier de methoden vinden de vijf grootste VSOP-componenten, en ook nog in dezelfde volgorde.
\({f}\) | |||||||
---|---|---|---|---|---|---|---|
\({Δ_1}\) | \({Δ_2}\) | \({Δ_3}\) | \({\min}\) | \({μ}\) | \({\max}\) | \({σ}\) | |
M1 | 0,1341 | −0,3413 | 0,1503 | −0,3413 | −0,0269 | 0,1828 | 0,2523 |
M2 | −0,0038 | −0,0877 | −0,0629 | −0,0877 | −0,0311 | 0,0813 | 0,0711 |
M3 | 0,0005 | −0,0841 | −0,0489 | −0,0841 | −0,0309 | 0,0526 | 0,0571 |
M4 | 0,0579 | −0,1510 | 0,0114 | −0,1543 | −0,0212 | 0,1299 | 0,1272 |
\({A}\) | |||||||
---|---|---|---|---|---|---|---|
\({Δ_1}\) | \({Δ_2}\) | \({Δ_3}\) | \({\min}\) | \({μ}\) | \({\max}\) | \({σ}\) | |
M1 | −0,0781 | −0,2460 | −0,0387 | −0,2460 | −0,1034 | −0,0323 | 0,0874 |
M2 | −0,0008 | 0,0015 | 0,0165 | −0,0008 | 0,0066 | 0,0165 | 0,0070 |
M3 | −0,0030 | 0,0019 | −0,0007 | −0,0030 | 0,0034 | 0,0125 | 0,0062 |
M4 | −0,0267 | −0,0780 | −0,0494 | −0,0780 | −0,0413 | −0,0127 | 0,0248 |
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−8 AE.
\({r}\) | ||||||||||
---|---|---|---|---|---|---|---|---|---|---|
# | \({f_\text{vsop}}\) | \({A_\text{vsop}}\) | \({f_1}\) | \({A_1}\) | \({f_2}\) | \({A_2}\) | \({f_3}\) | \({A_3}\) | \({f_4}\) | \({A_4}\) |
1 | 6283,0757 | 1670700 | 6283,1853 | 1478578. | 6283,0191 | 1663435. | 6283,0366 | 1652461. | 6283,0874 | 1606753. |
2 | 12566,1514 | 13956 | 12565,7423 | 9518,163 | 12566,0367 | 14080,45 | 12566,0130 | 13582,83 | 12566,0136 | 12592,51 |
3 | 77713,7734 | 3084 | 77713,5775 | 2621,098 | 6069,1215 | 3268,973 | 77713,7704 | 3080,909 | 77713,6995 | 2933,586 |
4 | 5753,3848 | 1628 | 7860,2648 | 1515,749 | 77713,7717 | 3083,366 | 5753,3858 | 1628,202 | 5753,4425 | 1592,821 |
5 | 7860,4194 | 1576 | 5753,5128 | 1475,316 | 5959,7856 | 2244,096 | 7860,4191 | 1575,710 | 7860,3537 | 1525,868 |
6 | 11506,7695 | 925 | 6069,5570 | 1233,091 | 5956,6438 | 2226,028 | 11506,7709 | 923,9234 | 11506,8306 | 852,9342 |
7 | 3930,2100 | 542 | 5960,2296 | 805,2115 | 5954,3899 | 2199,643 | 3930,2094 | 542,3906 | 3930,1726 | 537,9689 |
8 | 5884,9268 | 472 | 5957,0880 | 797,9427 | 5952,2461 | 2195,110 | 5884,9287 | 472,2971 | 5884,8817 | 466,1149 |
9 | 5507,5532 | 346 | 5953,9464 | 791,5154 | 6609,7268 | 1926,811 | 5507,6453 | 344,5373 | 5507,7109 | 328,0761 |
10 | 5223,6938 | 329 | 5952,6898 | 789,1468 | 6611,6113 | 1914,480 | 5223,7757 | 328,2695 | 5223,8040 | 326,4411 |
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).
\({f}\) | |||||||
---|---|---|---|---|---|---|---|
\({Δ_1}\) | \({Δ_2}\) | \({Δ_3}\) | \({\min}\) | \({μ}\) | \({\max}\) | \({σ}\) | |
M1 | 0,1096 | −0,4091 | −0,1960 | −0,4091 | 0,0623 | 0,5329 | 0,2746 |
M2 | −0,0565 | −0,1146 | −0,0017 | −0,3045 | 0,1103 | 0,9054 | 0,3990 |
M3 | −0,0391 | −0,1384 | −0,0030 | −0,1384 | −0,0003 | 0,0921 | 0,0631 |
M4 | 0,0117 | −0,1378 | −0,0740 | −0,1378 | 0,0038 | 0,1577 | 0,0920 |
\({A}\) | |||||||
---|---|---|---|---|---|---|---|
\({Δ_1}\) | \({Δ_2}\) | \({Δ_3}\) | \({\min}\) | \({μ}\) | \({\max}\) | \({σ}\) | |
M1 | −0,1150 | −0,3180 | −0,1501 | −0,3180 | −0,0863 | 0,0140 | 0,0994 |
M2 | −0,0043 | 0,0089 | −0,0002 | −0,0240 | 0,0004 | 0,0243 | 0,0118 |
M3 | −0,0109 | −0,0267 | −0,0010 | −0,0267 | −0,0045 | 0,0007 | 0,0086 |
M4 | −0,0383 | −0,0977 | −0,0488 | −0,0977 | −0,0348 | 0,0107 | 0,0347 |
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.
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/ν \).
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}
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}
Doe nu voor alle frequenties \( f_m \) tussen \( 2/T \) en \( 1/\D t_{\text{min}} \) in stappen van \( 1/(2T) \) het volgende:
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 \).
Sorteer de fases in oplopende volgorde. Die gesorteerde fases noemen we \( φ'_k \), met \( k \) lopend van 1 tot \( n \).
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.
Als voorbeeld nemen we een meetreeks waarin twee frequenties zitten, \( ν_1 = 0.087 \) en \( ν_2 = 0.136 \). Dat betekent dat er in de eerste reeks elke \( 1/0.087 ≈ 11.49 \) seconden iets gebeurt, en in de tweede reeks elke \( 1/0.136 ≈ 7.35 \) seconden. De eerste reeks loopt van \( t = 7.3 \) tot \( t = 145.23 \) (13 punten), en de tweede reeks loopt van \( t = 2.3 \) tot \( t = 142.01 \) (20 punten).
We kijken hoe onze kwaliteitsmaten het doen apart voor de meetreeks van elk van de twee frequenties, en ook voor de meetreeks waar de metingen van beide frequenties vermengd zijn.
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.
Voor onze meetreeks met alleen \( ν_1 \) vinden we dat \( ν_1t_i \bmod 1 = 0.6351 \) voor alle gebeurtenissen, en voor onze meetreeks met alleen \( ν_2 \) vinden we dat \( ν_2t_i \bmod 1 = 0.3125 \).
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 \).
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.
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 \( φ_ν \).
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.
Dit is te zien in figuur 9 voor de grafieken die horen bij de basisfrequentie en veelvouden (2× en 3×) daarvan. De plateau's zijn dan even lang maar liggen bij verschillende fases. Merk op dat het plateau grotendeels rechts van de diagnonaal ligt als het plateau een lage fase heeft, en grotendeels links van de diagnonaal als het plateau een hoge fase heeft. Als het plateau een fase rond 0.5 heeft dan ligt het ongeveer evenveel links en rechts van de diagnonaal. Daardoor hangt een maat zoals de gemiddelde afstand van de meetpunten tot de diagonaal af van de fase van het plateau.
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 \).
\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.
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.
//aa.quae.nl/nl/reken/frequentie.html;
Laatst vernieuwd: 2021-07-19