FFT command line application?
Estimating Multiple Frequencies

The Goertzel algorithm is more efficient than the DFT for a small number of frequency bins and can be easily implemented in open source software such as Octave or Python. 

More info is available on Wikipedia https://en.wikipedia.org/wiki/Goertzel_algorithm including the handy rule of thumb as to when it is more efficient:

$$M \le \frac{5N_2}{6N}log_2(N_2)$$

Estimating A Single Tone in Low SNR Conditions

For estimation of a single tone in low SNR conditions see this paper by Rim Elasmi_Ksibi, Hichem Bessbes, Roberto Lopez-Lacarace and Sofiane Cherife: https://www.researchgate.net/publication/220228033_Frequency_estimation_of_real-valued_single-tone_in_colored_noise_using_multiple_autocorrelation_lags which extracts the estimate for $cos(\omega_0)$ from the samples of the autocorrelation, where $\omega_0 \in [0, \pi] $.  Which derives the estimate for $\omega_0$ as: 

$$\cos(\hat{\omega_0}) = \frac{\sum_{k=p}^q \hat{r_k}(\hat{r}_{k-1}+\hat{r}_{k+1})}{2\sum_{k=p}^q \hat{r}_k^2}$$

Where $\hat{r}_k$ is the unbiased autocorrelation for the observed samples y given as:

$$\hat{r}_k = \frac{1}{N-k}\sum_{n=k+1}^N y_n y_{n-k}$$

and p and q are any integers q>p large enough such that the noise in the samples compared are independent.  The larger the range of q-p the more processing required but the lower the noise in the estimate (so you can make that trade). If you choose a p that is less than the lag for noise independence, then you will have no advantage from additional processing for p+n until p+n is at the lag where the noise is independent. You can access this from the autocorrelation of the noise process alone to determine the lag at which the autocorrelation is 0. For example with white noise the autocorrelation = 0 for any offset, meaning all noise samples are indpendent in which case p can be as low as 1.  

Intuitive Explanation of Autocorrelation As a Frequency Discriminator

The above gives the actual frequency estimate with all scaling accounted for, and you can trade the computational complexity with the noise of the estimate approaching the CRLB as detailed in the referenced paper. What follows is to provide and intuitive understanding of how this works. This works on the simple principle that the product of a sinusoid and a phase shifted version of the same sinusoid is scaled by the cosine of the phase as given by the following trigonometric identity:

$$cos(\alpha)cos(\beta) = cos(\alpha+\beta) + cos(\alpha-\beta)$$

So when the frequency is the same and only the phase is different, the product is:

$$cos(\omega_c t+\phi)cos(\omega_c t) = cos(\phi) + cos(2\omega t + \phi)$$

If we average (low pass filter) the above, the $cos(2\omega t+ \phi)$ term goes to zero and we are left with $\cos(\phi)$. This shows how the product is a phase detector. When we delay and multiply (as done in the autocorrelation!), the delay produces a signal with the same frequency but a phase shift. The resulting phase measured by the phase detector (product) is the change in phase over that delay which by definition is frequency! (Frequency is the derivative of phase).

A commmon frequency discriminator topology is to delay and multiply (a frequency discriminator produces an output value that is proportional to the frequency of the input): 



Each sample of the Autocorrelation Function is a delay and multiply with a different delay value for each. The above referenced paper is simply scaling each result back to be $cos(\omega)$ and averaging to minimize the noise contribution and improve the estimate. In the plot below the vertical axis is crossing the horizontal axis at $-\pi/2$ to be at the point of maximum slope (operating point when used as a discriminator): 



This all applies to complex signals as well, in which case a complex conjugate multiply is done as shown with phase detector topologies below comparing real signals to complex signals. This suggests for a single complex tone the use of either the real or imaginary output of the complex conjugate multiplication to get a similar $cos(\omega_0)$ (real out, I) or $sin(\omega_0)$ (imag our, Q) result but with further processing a direct result of $\omega_0$ is obtained using:

$$\omega_0=atan2(Q, I)$$

Where atan2 is the 2-argument arctangent with I and Q are the real and imaginary results of the complex conjugate multiplication, suggesting how the referenced approach for a single real sinusoid can also be extended to the case for a single complex tone.

And for a single complex tone in high SNR conditions the estimate is trivial since the normalized angular frequency is the phase change from one complex sample to next, which is readily extracted from complex conjugate multiplication of the two samples:

$$Ae^{j\omega_0} = y[n-1]y[n]^*$$

With $\omega_0$ extracted using the atan2 function on the real (I) and imaginary (Q) result of the product as $atan2(Q,I)$. This ends up with the following in terms of $y[n]=I[n]+jQ[n]$:

$$\omega_0=atan2 \bigg( \frac{I[n]I[n-1]+Q[n]Q[n-1]}{I[n]Q[n-1]-Q[n]I[n-1]}\bigg)$$

(And there are numerous efficient estimators for the atan2 process that can be used to further simplify this, including the iterative CORDIC rotator when cycle times to iterate are  more available than multipliers and look up tables.)



What is useful and interesting from this is the imaginary portion of the autocorrelation function for any waveform will be proportional to the frequency offset of that waveform, which is useful for carrier recovery implementations for radio receivers! This is demonstrated below in the result for autocorrelation of a complex additive white Gaussian noise signal with a frequency offset in one direction ($e^{j\omega_o t}$) as plotted on a complex plane showing the real and imaginary terms of the autocorrelation.