This problem relates to digital signal processing (DSP). For serious acquaintance with the theme can be recommended the book of Blahut
Methods for this problem solving are divided into linear (based on FFT) and non-linear (based on autoregressive models). Typical for linear methods are using of FFT algorythm with weighted Hamming window. Linear model presents signal as sum of known fixed harmonics with the constant step, and it's easy to realization.
The only drawback of the method is the limited accuracy of the frequency measurement.
But if you want to discover the frequency beyond the scale of the FFT, then the best choice becomes Durbin method based on AR model. I already wrote about this method on the Russian site, where you can find texts of programs. Here I give the translation of the article.
One of the best forecasting methods - autoregressive. In any case, the periodic bursts catch it.
Signal Model
The method - a model of autoregression (AR) of order $k$ for sample $i=0,1,\dots, n-1:$
$$x_i + f_0x_{i-1} + f_1x_{i-2} + f_2x_{i-3} + \dots + f_{k-1}x_{i-k} = 0,$$
where $i = k, k+1,\dots, n-1.$
The desired frequencies are the roots of the characteristic equation
$$1+f_0\omega+f_1\omega^2+f_2\omega^3+\dot+f_{k-1}\omega^{k}=0.$$
The order of the model should correspond approximately to the complexity of signal.
Toeplitz symmetry
According to the method of least squares, the discrepancy
$$\varepsilon = \left<(x_i + f_0x_{i-1} + f_1x_{i-2} + f_2x_{i-3} + ... + f_kx_{i-k})^2\right>_{i=k\dots n-1}$$
should be minimized.
For this, should equate to zero the partial derivatives by the AR coefficients $f_s$ for $s = 0\dots k-1$.
This leads to the following system of equations special:
\begin{cases}
a_0f_0 + a_1f_1 + a_2f_2 + ... + a_{k-2}f_{k-2} + a_{k-1}f_{k-1} = -a_1,\\
a_1f_0 + a_0f_1 + a_1f_2 + ... + a_{k-3}f_{k-2} + a_{k-2}f_{k-1} = -a_2,\\
a_2f_0 + a_1f_1 + a_0f_2 + ... + a_{k-4}f_{k-2} + a_{k-3}f_{k-1} = -a_3,\\
\dots\\
a_{k-2}f_0 + a_{k-3}f_1 + a_{k-4}f_2 + ... + a_0f_{k-2} + a_1f_{k-1} = -a_{k-1},\\
a_{k-1}f_0 + a_{k-2}f_1 + a_{k-3}f_2 + ... + a_1f_{k-2} + a_0f_{k-1} = -a_k,
\end{cases}
where $a$ is the autocorrelation function (ACF).
The first feature is the form of LHS matrix with the same coefficients on the main diagonal and that all the coefficients on the diagonal, equidistant from the main. Matrix with the symmetry of this type are called Toeplitz, and Levinson algorithm is applicable to them.
The second feature - that the free column members formed of the same elements of the matrix $a$, and the algorithm of Levinson-Durbin is applicable to such a system.
The third feature is in fact that considering matrix is Toeplitz only approximate. For example, for a second-order model, the exact form of the equations is:
$$\begin{cases}
\left<x_{i-1}^2\right>f_0 + \left<x_{i-1}x_{i-2}\right>f_1 = \left<x_ix_{i-1}\right>\\
\left<x_{i-1}x_{i-2}\right>f_0 + \left<x_{i-2}^2\right>f_1 = \left<x_ix_{i-2}\right>,
\end{cases}$$
where $i$ is the index of the amount of each type of $\left<u_i\right>$ runs from $2$ to the maximum.
It is not difficult to note that:
$1)$ the items of the main diagonal correspond to element zero autocorrelation function samples for computation are shifted by one element.
$2)$ the items on the secondary diagonal and the first item in the column of free terms correspond to the first element of the autocorrelation function, within diagonal elements are equal, and sample for free member moved toward them.
This can lead to unexpected effects for sequences with a strong increasing (decreasing) trend.
So, for the sample of the first $10$ members of the Fibonacci sequence
$x_i = \{1, 1, 2, 3, 5, 8, 13, 21, 34, 55\}$
we obtain the system of equations in form
$$\begin{cases}
1869f_0 + 1155f_1 = - 3024,\\
1155f_0 + 714f_1 = - 1869
\end{cases}$$
with the correct solution $f_0 = -1,\ f_1 = -1,$ but the expected Toeplitz symmetry of the problem doesn't exist.
The use of the ACF system gives
$$\begin{cases}
4893f_0 + 3024f_1 = - 3024,\\
3024f_0 + 4893f_1 = - 1869,
\end{cases}$$
using of which in this case is inappropriate. The main method of dealing with these effects - centering of the sample. But the main moment is checking the possibility of a serious simplification of the model for using of Durbin algorithm.
Durbin algorithm
The Durbin idea is to search for a solution as the matrix of dimension $n+1$ in the form of:
$$\left(f'_0, f'_1,\dots, f'_{n-1}, f'_n\right) = \left(f_0, f_1, ..., f_{n-1}, 0\right) + \beta\left(f_{n-1},f_{n-2},\dots,f_0,1\right).$$
Indeed, the substitution of this solution in a matrix of dimension $n+1,$ taking into account the dimension $n,$ gives the system:
$$\begin{cases}
-a_1 - \beta a_n + a_n\beta = -a_1,\\
-a_2 - \beta a_{n-1} + a_{n-1}\beta = -a_2,\\
\dots\\
-a_{n-1} - \beta a_1 + a_1\beta = -a_n,\\
a_nf_0 + a_{n-1}f_1 + a_{n-2}f_2 + ... + a_2f_{n-2} + a_1f_{n-1}\\ + \beta (a_nf_{n-1} + a_{n-1}f_{n-2} + a_{n-2}f_{n-3} + ... + a_2f_1 + a_1f_0 + a_0) = -a_{n+1}.
\end{cases}$$
$N$ the first equations of the system are satisfied for any value of $\beta$, the latter - at
$$\beta = - \dfrac{a_{n+1} + a_nf_0 + a_{n-1}f_1 + a_{n-2}f_2 +\dots+ a_2f_{n-2} + a_1f_{n-1}}{a_nf_{n-1} + a_{n-1}f_{n-2} + a_{n-2}f_{n-3} + \dots + a_2f_1 + a_1f_0 + a_0}.$$
durbin() function returns an array of AR vectors for a total of $n$ range.
For $n=1\ f = -\dfrac{a_1}{a_0},$ the following vectors of coefficients are calculated recursively (first $\beta$, and then $f$).
Software implementation
The program in PHP has the following features:
Centering array $center()$.
The scalar product of the vectors of the second vector shifted $scalar\_prod()$.
The output of the array with a text commentary $print\_array()$.
The output of the system of linear equations of the second order and its solutions $print\_s()$.
The calculation of the autocorrelation function (ACF) $acf().$
Compare accurate and Toeplitz second order systems and solutions for a given sample $compare\_s().$
Durbin algorithm for solving the system of equations of a special form $durbin()$.
Test Durbin algorithm $test\_durbin()$.
Results
The program showed high efficiency in a limited (sine) sample with a large volume of raw data.