Astronomy Answers
Frequency Determination

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

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

This page explains some relatively quick methods for estimating characteristics (such as the frequency) of periodic components of a set of measurements.

1. From Measurements at Regular Intervals

Fig. 1: Measurements Diagram

Fig. 2: Meetreeksgrafiek met proefreeks

1.1. Use a Fourier Spectrum

A fast method to get some insight is to determine the Fourier spectrum of the data. One can use an FFT algorithm for this, if the time between successive measurements is constant. (FFT stands for Fast Fourier Transform.) If a programming language has a function for this, then it is usually called something like POWER or FFT. The Fourier spectrum shows how much "energy" there is in the data for each whole frequency (i.e., for sines and cosines that exactly fit a whole number of times in the data). That energy is equal to the square of the amplitude. Mathematics tells us that it is always possible to decompose an arbitrary data set into harmonic components in this way. It is much better for the calculations if the data set has an average value equal to 0, so if this is not yet the case, then we calculate the average value and subtract it from the data set first.

Fig. 3: Frequency Spectrum

It is clear that the greatest energy is at a frequency of 10, but there is also activity at other frequencies. There is some sort of a peak and a valley near frequency 80, too. We saw before that a frequency of 10 does not fit well enough. A Fourier spectrum only looks at harmonic components that exactly fit the data set a whole number of times, and it would be quite a coincidence if the data set consists of just one harmonic component (a combination of a sine and a cosine with the same frequency) that just happens to exactly fit the data set a whole number of times, so even if there really is only a single harmonic component in the data set, chances are that a Fourier spectrum of that data set will show activity at every frequency.

The advantage of a Fourier spectrum is that it can be calculated relatively quickly. The disadvantages are that a harmonic component is only confined to a single frequency in a Fourier spectrum if that component happens to fit the data set exactly a whole number of times, and that different harmonic components usually distort each other in such a Fourier spectrum, like the peak near frequency 80 that looks rather ugly with a dip in it, because around that frequency there is also a significant contribution from the peak at frequency 10.

1.2. Simple Estimates

1.2.1. Simplest Estimate: Method 1

The simplest estimate of the frequency of a harmonic component from a Fourier spectrum is to locate all peaks in it (where there is greater energy than at neighboring frequencies) and to measure the frequency and amplitude (the square root of the energy) of those peaks.

For our data set, this yields

# 

so there is also a peak at frequency 32. This peak is far more difficult to see in the diagram than the peak at frequency 78, yet the estimate for the amplitude at frequency 32 is a lot higher than that at frequency 78. That is because of the "slope" from the peak at frequency 10.

1.2.2. Better Estimate: Method 2

If you calculate the Fourier spectrum of a single harmonic component with positive frequency \( f_0 \) that need not fit the data set a whole number of times (i.e., also if \( f_0 \) is not an integer number), then you find that the energy \( P_f \) at frequency \( f \) is equal to about

\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}

where \( β \) is the phase of the harmonic component at the beginning of the data set. The sinc function has its greatest value when its argument is equal to 0, so the terms with \( f - f_0 \) in it are greatest near \( f = f_0 \) and the terms with \( f + f_0 \) in it are greatest near \( f = -f_0 \). So, there are two equally high peaks, symmetrically around frequency 0: one peak at positive frequency and one peak at negative frequency. The terms from the negative peak (the peak at negative frequency) are smaller at positive frequencies than the terms from the positive peak are. If we ignore the terms from the negative peak (with \( f + f_0 \) in it, which have little influence for \( f \gg 0 \)), then the previous equation leads to

\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}

This no longer depends on \( β \). If \( f_x \) is the integer frequency at which the spectrum has its highest peak, and \( f_0 = f_x + k \), then \( k \) is between −1 and +1. With

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

the preceding two numbered equations lead to

\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}

With this one can find the frequency and amplitude of all peaks from the frequency spectrum. The procedure is as follows:

1. Find the integer frequency \( f_0 \) where the frequency spectrum shows greater energy than for neighboring frequencies
2. Use equation \ref{eq:k1} to estimate the frequency difference \( k \). The frequency \( f \) of the component can then be estimated at \( f_0 + k \). (If \( |k| \) is greater than 1, then the assumptions upon which the current method is based do not hold for the current data set. See below.)
3. Use equation \ref{eq:A1} to estimate the amplitude of the component.

If we apply this to the data set then we find

# 

The amplitude of the largest peak is now estimated to be quite a bit higher, and also the amplitude of the peak at frequency 33.

1.2.3. Even Better Estimate: Method 3

The strange shape (in the preceding diagram) of the peak at frequency 78 and the near-invisibility of the peak at frequency 33 derive from the broad "slopes" that come from the tallest peak, which disturb the lower peaks. So, we must try to make the slopes of the tallest peak more steep so that they don't distort the other peaks so much. We do that by multiplying the data set with a mask or apodization function that goes to zero at the beginning and end of the data set. The apodization function we'll use here is a cosine that has values between 0 and 2, has exactly one period in the data set, and goes to zero at the edges. (This is a form of the Hanning mask or Hanning function.) When such a mask is applied, it is important that the average of the data set be equal to zero.

Fig. 4: Frequentiespectrum

Note that the vertical axis of this diagram has a much greater range than the previous diagram. In the previous spectrum, the energy ratio between the top and bottom of the diagram was about a thousand million (10⁹), but in this diagram it is ten thousand million times greater still (or about 10¹⁹ in total). The peaks are slightly fatter, but the slopes are far steeper, so the smaller peaks stand out much more here than in the previous diagram.

In the same way as before one can calculate the Fourier spectrum of a single harmonic component of frequency \( f_0 \) and amplitude \( A \). Just like before, we ignore the contributions from the negative peak. We then find

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

This leads (in a similar way as before) to the estimate

\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}

(If \( |k| \) is greater than 1, then the assumptions upon which the current method is based do not hold for the current data set. See below.)

Applied to the data set, this yields

# 

Now the peak at frequency 78 is estimated to be higher than the peak at frequency 32.

1.2.4. Method 4: General Quadratic

The two preceding methods assume that the components in the Fourier spectrum are pure harmonics, which means that they are assumed to have a fixed frequency and amplitude, but this is quite rare in practice. Usually the frequency and amplitude of all components vary a bit with time, which makes the peaks in the Fourier spectrum wider. If the variation gets too large, then those methods might not give any estimate for how much the "real" frequency differs from the Fourier peak. In such cases you can make the simpler assumption that the relationship between energy and frequency near the peaks in the Fourier spectrum is quadratic. Then we find (in a similar way as before) that

\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}

You can apply a mask here, too, to make the slopes of the peaks more steep. Applied to our data set this method yields

# 

1.2.5. Remarks

1. The resolving power of a frequency spectrum is limited. If the measurements last for a time \( T \), then the frequencies are multiples of \( 1/T \). Frequency differences smaller than \( 1/T \) cannot be measured from such a frequency spectrum. Two harmonic components of which the frequencies differ from each other by less than that amount cannot be told apart in the frequency spectrum.
2. Fig. 5: Frequency Spectrum Showing Aliasing

   
   A harmonic component with a frequency that is greater than the greatest frequency from the spectrum does show up in the spectrum, but at a frequency that is smaller than its own frequency. If the data set lasts a time \( T \) and contains \( N \) data values, then the greatest frequency in the spectrum is about equal to \( \frac{N}{2T} \). If the frequency \( f \) of a harmonic component is greater than that greatest frequency, then that component ends up in the spectrum at a frequency equal to \( \left| f - N \left[ \frac{f}{N} \right] \right| \). This is called aliasing. An example is shown in Figure 5, which shows three different harmonic components that all intersect in the measurement points (indicated by small squares). One cannot tell from measurements at those point which of the components it is.
3. If the data set also contains a component that is not periodic (or perhaps is periodic at a period that is far greater than the length of the data set), then the energy of that component shows up at the lowest frequencies, where it distorts any "true" peaks and boosts them to undeservedly greater amplitudes. To suppress this effect, one can subtract a running average from the data set, but that suppresses not just non-periodic components but also periodic components with periods that are longer than the width of the running average.
4. Methods 2 and 3 ignore the influence from negative peaks, which is strongest at the smallest frequencies (still smaller than the influence of the positive peaks, but not a whole lot smaller). Because of this, the results will tend to be less accurate for the smallest frequencies than for the greater frequencies.
5. The steepening of the slopes of the peaks that method 3 yields has a price, which is that the peaks get fatter at their top than they were for method 2. If two peaks are very close together, then method 3 can have more trouble to separate them than method 2 has. It seems that method 2 is sometimes more accurate than method 3 at finding the characteristics of the greatest component, but method 3 is almost always better at finding the smaller components.

1.3. The Accuracy

1.3.1. Test with the Data Set

The harmonic components that compose the data set are as follows:

# 

Now we can determine how accurate the various estimates were. The next table shows the difference between the estimate and the true value for the frequency (\( Δf_0 \)) and the relative difference for the amplitude (\( \frac{ΔA}{A} \)).

# 

Method 3 seems best in general (if the data set contains only pure harmonic components). This is confirmed by a large number of tests with various randomly constructed data sets (with an exponential distribution of amplitudes). The accuracy of the frequency appears to be roughly constant for method 1, and, for methods 2 and 3, to increase with about the 2.7-th power of the average spacing (in the spectrum) between the harmonic peaks, but with an advantage of about a factor of 10 for method 3 compared with method 2. Method 4 is much better than method 1, but often worse than methods 2 and 3 if the data set contains only pure harmonic components. If the components are not pure harmonics, then methods 2 and 3 might not work at all, and then method 4 is best.

The relative accuracy of the amplitudes seems roughly constant for method 1, to increase with about the square of the spacing (in the spectrum) between the harmonic peaks for method 2, and with about the third power of the average spacing for method 3.

1.3.2. Testing with the VSOP Model

Below, we test the methods on the VSOP model of the motion of the planets. One can calculate the positions of the planets from this model for a few thousand years into the future and into the past. The model consists of a large number of harmonic components, so we can test how well we can find those components again in the positions of the planets.

Below is a table that shows the ten largest harmonic components in the ecliptic heliocentric longitude \( λ \) of the Earth relative to the equinox and ecliptic of the date according to the VSOP model (\( f_\text{vsop} \) is the frequency in radians per Julian millennium of 365250 days, \( A_\text{vsop} \) is the amplitude in 10⁻⁸ radians), as given by [Meeus], and also the ten largest harmonic components calculated from \( λ \) (calculated for each 10 days during a period of 10,000 years around the year 2000) using methods 1-4 (the other \( f \) and \( A \), in the same units as for the VSOP model). To suppress non-periodic components at small frequencies I removed a running average over 68 years from the data sets before trying to detect harmonic components in them.

 # 

Among the 10 greatest components found through method 1 there are 6 of the 10 greatest components from the VSOP list. For method 2 this number is 5, and for methods 3 and 4 it is 9. The next two tables display the absolute (for the frequency \( f \)) or relative (for the amplitude \( A \)) error for the first three VSOP components (\( Δ_1 \), \( Δ_2 \), \( Δ_3 \)), and also the minimum (\( \min \)), average (\( μ \)), maximum (\( \max \)), and standard error (\( σ \)) of that same kind of error for the first ten VSOP components, for methods 1-4 (M1 - M4). To determine such an error for each of the VSOP components, I looked for the estimated component (from methods 1-4) that was closest to it.

For example, method 2 underestimates the frequency of the greatest VSOP component by 0.0565 units, and underestimates the amplitude by a fraction of 0.0044 (i.e., by 0.44%). Method 3 is usually the best one (with the smallest \( σ \) that is yet much greater than the \( μ \)) for both the frequency and the amplitude.

Next is a similar table for the ecliptic heliocentric latitude of the Earth. Meeus provides only 5 components for this case, but I did, when calculating the positions, also take the smaller components into account that Meeus does not mention.

 # 

All four methods find the five greatest VSOP components, and in the correct order, too.

Here, too, method 3 is usually best, and method 2 is almost as good. Methods 1 and 4 systematically underestimate the amplitude.

Next is a similar table for the distance \( r \) between the Earth and the Sun, measured in units of 10⁻⁸ AU.

 # 

Among the ten greatest components found by method 1 there are 5 of the ten greatest VSOP components. For method 2 the number is 3, and for methods 3 and 4 it is 10 (and all in the correct order, too).

Here, too, method 3 is usually best. For the frequency, method 4 is second best, but for the amplitude it is method 2. For the frequency, method 2 is usually even worse than method 1. Methods 1 and 4 systematically underestimate the amplitude.

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.

languages: [en] [nl]

//aa.quae.nl/en/reken/frequentie.html;
Last updated: 2021-07-19