Astronomy Answers: Frequency Determination

Astronomy Answers
Frequency Determination


[Astronomy Answers] [Dictionary] [AnswerBook] [Universe Family Tree] [Science] [Starry Sky] [Planet Positions] [Calculate] [Colophon]

1. From Measurements at Regular Intervals ... 1.1. Use a Fourier Spectrum ... 1.2. Simple Estimates ... 1.2.1. Simplest Estimate: Method 1 ... 1.2.2. Better Estimate: Method 2 ... 1.2.3. Even Better Estimate: Method 3 ... 1.2.4. Method 4: General Quadratic ... 1.2.5. Remarks ... 1.3. The Accuracy ... 1.3.1. Test with the Data Set ... 1.3.2. Testing with the VSOP Model ... 2. Uit tijdstippen van gebeurtenissen ... 2.1. De Methode ... 2.2. Voorbeelden ... 2.3. De Uitleg

\(\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. 1: Measurements Diagram
In astronomy many things have periodical components. Every celestial object with a closed (elliptical) orbit regularly returns to about the same position in its orbit. Some stars get periodically brighter and dimmer again. The magnetical activity of the Sun increases and decreases again. Suppose that you have a set of measurements of some characteristic of which you believe it may contains such periods. Take, for example, the measurements displayed in Figure 1. During 1024 seconds, one measurement was taken each second, and the results are shown in the diagram (with time along the horizontal axis and the measured values along the vertical axis). How can you find the periods that contribute to this series of measurements?

Fig. 2: Meetreeksgrafiek met proefreeks
Fig. 2: Meetreeksgrafiek met proefreeks
One method is to guess at the period, and to compare the measurements to what they should be if your guess is correct. Figure 2 shows the same data again, but also some test data for a guessed harmonic component with a frequency of 10 (i.e., with exactly 10 full periods in the data set) and an amplitude (greatest deviation from the average) of 1. The test data is reasonably close, but not good enough. It is not tall enough, and it still runs slowly out of step with the measurements. It is a lot of work to keep adjusting your guess and try again.

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
Fig. 3: Frequency Spectrum
Figure 3 shows the frequency spectrum of our data set. The horizontal axis shows frequency \( f \), which says how many times the harmonic component fits in the observations. The vertical axis shows the energy \( P \), which is the square of the amplitude.

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

# 
\({f_0}\) \({A}\)
1 10 0.915
2 32 0.0065
3 78 0.0024

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

# 
\({f_0}\) \({A}\)
1 10.30033 1.066
2 32.66 0.016
3 78.26 0.0027

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
Fig. 4: Frequentiespectrum
The frequency spectrum of the apodized data set is shown in Figure 4. For convenience I again show the previous frequency spectrum (without the data apodization), too, but now as a dashed line.

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 (109), but in this diagram it is ten thousand million times greater still (or about 1019 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

# 
\({f_0}\) \({A}\)
1 10.29955 1.04978
2 77.6024 0.0019970
3 31.60 0.0012

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

# 
\({f_0}\) \({A}\)
1 10.18609 1.001697
2 77.6994 0.0018485
3 31.25 0.00095

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
    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:

# 
\({f_0}\) \({A}\) \({β}\)
1 10.3 1.05 36°
2 77.6 0.002
3 31.4 0.001

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

\({Δf_0}\) \({\frac{ΔA}{A}}\)
# 
Method 1 Method 2 Method 3 Method 4 Method 1 Method 2 Method 3 Method 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

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−8 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.

\({λ}\)
 # 
\({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

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.

\({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

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.

\({β}\)
 # 
\({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

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

\({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

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−8 AU.

\({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

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

\({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

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

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.

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.

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 \).

Fig. 6: Fase voor verschillende frequenties, één frequentie in de data
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
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
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
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.

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 \).

  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.



[AA]

languages: [en] [nl]

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