0
$\begingroup$

Are there any command line programs for Windows, preferably free and stand-alone, which can report the peak/strongest frequency within a given range of frequencies? I need something like this to automate finding the frequency of a calibration signal which slowly drifts.

$\endgroup$
  • $\begingroup$ Hi! Sorry, you're asking for a program that fulfills your specifications – these questions are best-case borderline off-topic. You could, however, write a program yourself that does that. Can you tell us more about your signal, and the tools you're using? (also, the title stands somewhat in conflict with your question's body, so this really calls for more background on what signal you're dealing with, and how much precision with how much observation you need) $\endgroup$ – Marcus Müller Jun 22 '19 at 18:18
  • $\begingroup$ I need the strongest frequency from within a range with white noise and 1 strong calibration sine signal. The added calibration signal shifts as the receiver has no TCXO, so this allows me to determine its characteristics and calibrate the other signals. From this I can determine the calibration signals' amplitude and relate to its known power. I want to automate this as there are 576 audio files to process. 'FFT 1.0' by Lionel Loudet sidstation.loudet.org/fft-en.xhtml does FFT from the command line. I coded a script to extract the info but it takes a lot of processing time. $\endgroup$ – Petoetje59 Jun 23 '19 at 14:18
  • $\begingroup$ why is the script slow? I'm not sure a standalone program would be faster, since in scripting languages, you'll typically just call a FFT function from a library (that is very fast), so I'd assume the overhead is negligible. $\endgroup$ – Marcus Müller Jun 25 '19 at 6:57
  • $\begingroup$ In this script the fft.exe is the culprit - it gobbles up most of the processing time. In all it takes 732 secs to process just a 5 minute file, so visually determining the calibration frequency by just watching it in SpecLab is way faster, but still very time consuming to do. $\endgroup$ – Petoetje59 Jun 26 '19 at 12:02
  • $\begingroup$ yes, but having a different program do the same will not be faster. You need to have a better way of estimating the frequency than to use the FFT, not a different program to encapsulate the FFT. $\endgroup$ – Marcus Müller Jun 26 '19 at 12:13
1
$\begingroup$

The calibration tone is a pure sine wave. I solved the problem without resorting to FFT - by using the "Sine fitting algorithms" (4-parameter method) described in an annex to IEEE-STD-1057/1241.

| improve this answer | |
$\endgroup$
1
$\begingroup$

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

Frequency Discriminator

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

Autocorrelation as a Freq Disc

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

phase detectors

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.

Autocorrelation for complex frequency offset

| improve this answer | |
$\endgroup$
  • $\begingroup$ Hey, I fixed the formula for you. There was a "+' missing which made the numerator third order and the denominator second order. If you double the amplitude of the signal, the result has to stay the same. This formula looks like a boat load of calculations, but I am going to compare it to mine. Might take a while to get to though. $\endgroup$ – Cedron Dawg Dec 28 '19 at 21:17
  • $\begingroup$ Thanks for reading- this is for low SNR where you have the option to use a boatload of computations to the extent you want to reduce noise (you select). In high SNR you can just do the multiplication of two samples delayed from each other to get $Acos(\omega)$ and the ideal delay to use in this case is close to where the phase would be $\pi/2$ given the result has the highest slope, hence sensitivity. For complex tones with high SNR it is even easier since you can just do the complex conjugate multiplication between any two samples within $2\pi$ rotation and then use atan2/N to get $\omega$! $\endgroup$ – Dan Boschen Dec 28 '19 at 21:27
  • $\begingroup$ It's even easier than that with high SNR in the time domain. $$ \frac{\hat{r}_{k-1}+\hat{r}_{k+1} }{ \hat{r}_k } $$ has the same form as $$ \frac{ y[k-1]+y[k+1] }{ y[k] } $$ which is known as Turner's three point formula. I've generalized this into various families of formulas in three blog articles: dsprelated.com/showarticle/1051.php, dsprelated.com/showarticle/1056.php, and dsprelated.com/showarticle/1074.php. Yes, I'm curious to compare. Boat loads of calculations don't cost anywhere near what they used to. ;-) $\endgroup$ – Cedron Dawg Dec 28 '19 at 22:21
  • $\begingroup$ Check out my followup. That is with my original 3-bin formula. My 2-bin formula is even better except when the frequency is really close to a bin, then the 3-bin formula is a little more robust. Both are exact in the noiseless case. BTW, thanks for the paper reference! $\endgroup$ – Cedron Dawg Dec 28 '19 at 23:05
  • $\begingroup$ Cool - thanks Cedron. Would be cool to compare apples to apples including processing metrics as I have a use for efficient tone estimators (in fact I have an approach that is done with just adds and compares)— your approach is on DFT bin’s correct? Could it work with a 2 point DFT (would be ideal) or do your exclude the DC bin because of divide by zero issues? When you say “first of it’s kind” for an exact estimate with no noise, isn’t the cosine(angle) from the product exact? $\endgroup$ – Dan Boschen Dec 28 '19 at 23:21
0
$\begingroup$

If you have a high SNR signal and can bookend your signal with two DFT bins, I don't think you can do better than https://www.dsprelated.com/showarticle/1095.php with an implementation shown in https://www.dsprelated.com/showarticle/1284.php if accuracy is important. For a noiseless, non-integer frequency (cycles per frame), the formula works with any two bins, as it is mathematically exact. For frequencies near the middle between two bins, it is the most robust by my testing.


Here is the graphic comparison which is sparking my curiosity:

From the Elasmi-Ksibi, et. al., paper cited by D.B.:

Fig 3 MSE vs SNR

This shows the autocorrelation formula's results.

From: http://www.tsdconseil.fr/log/scriptscilab/festim/index-en.html

MSE vs SNR

There is a full write up available as a PDF in the upper right corner under the link "Resources: Comparison of different frequency estimation algorithms (pdf)"

This shows the comparisons of several formulas. Notice that the green line, which is my original 3-bin real valued formula, is the only one that hugs the line into the high SNR territory. (It is exact in the noiseless case, first of its kind.)

What is more interesting though, is that it also seems to do better in the low SNR range. This may not be a valid comparison based on noise type, not sure.

| improve this answer | |
$\endgroup$
  • 1
    $\begingroup$ (continued from under D.B.'s, curse the comments to chat policy) They were really clever in the autocorrelation formulas by multiplying eq (3) by $r_k$ before summing, ensuring the denominator will be a sum of squares and thus always non-zero. $$ $$ The near instantaneous article formulas work on either real or complex tones. Quite a difference to the DFT ones. The main point of these formulas is to be able to calculate the frequencies in a short duration, much shorter than a cycle. For a rapidly changing tone, this is an advantage. $\endgroup$ – Cedron Dawg Dec 29 '19 at 0:18

Not the answer you're looking for? Browse other questions tagged or ask your own question.