Astronomy Answers
Frequency Determination

\( \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. Introduction

Fig. 1: Measurements Diagram

Fig. 2: Meetreeksgrafiek met proefreeks

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

3. Simple Estimates

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

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

where \( \sinc_{π}(x) \) is a function equal to \( \frac{\sin(π x)}{π x} \) if \( x \) is unequal to 0, and 1 if \( x \) is equal to 0, and \( β \) 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.

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

3.4. Method 4: General Quadratic

The two preceding methods assume that the components in the Fourier spectrum are pure haromics, 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

# 

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

4. The Accuracy

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

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

languages: [en] [nl]

http://aa.quae.nl/en/reken/frequentie.html;
Last updated: 2016−02−07