Efficient calculation of exact frequency from a broad frequency range (1 Hz - 1000 Hz)

Question

I have some input signal sampled at 1-10 MS/s. Let's assume that it contains one strong frequency component (with frequency anywhere in the range from 1 Hz to 1000 Hz or even higher, changing over time) plus a number of higher harmonics (with quite a bit lower amplitude than the main frequency) and some noise.

The frequency may rise from zero to 1 kHz in a few seconds. The idea is to calculate the result approximately every 10-20 periods of the signal.

I need an algorithm which is as fast as possible to calculate the precise main frequency in real time.

I'm currently testing a couple of methods (like this or this) which can already give a pretty accurate result (under 0.01% when there are not too many higher harmonics) if I know an approximate frequency upfront and feed the FFT with cca. 10 periods (if the signal is around 100 Hz, I would use 100 ms slices; in case of 1 Hz, I'm ok with a lower accuracy, so taking just 1 - 2 periods is fine in that case). All these algorithms use just two or three bins from the FFT (around the highest amplitude) to calculate the exact frequency.

The problem is that in order to cover a wide frequency range without any idea what the base frequency is, I would first need to calculate FFT a lot of times. For example, first calculate all chunks of 10 ms just in case there is some 1 kHz component present. If it's not, maybe switch to 20 ms, then 40 ms ... and so on until we find the first signal length which adequately covers cca. 10 periods of the signal. That might slow down the algorithm 10-fold compared to the situation when frequency is approximately known upfront.

My question is: is there any way to first calculate 10 ms FFTs for each chunk of the signal, and if the frequency turns out to be too low, somehow combine results from two neighbouring 10 ms + 10 ms FFTs with a method that would require less computation resources than calculating FFT of the full 20 ms from scratch (given that I probably only need one or two FFT bins for a 20 ms signal anyway, and in this case I should know which ones).

Or any other idea of the fastest way to determine an approximate frequency quickly? (Some even use zero-crossing to do that which works fine as long as you are looking at nearly perfect sine curve, but that's not useful when higher frequencies and noise are present.)

Answer 1

Thanks for citing my blog article on frequency calculation. Based on your emphasis on efficiency I am going to recommend a different approach which should get you equivalent or superior results with way fewer calculations.

The first thing you should do is smooth your signal heavily with exponential smoothing. This will squash your harmonics and mitigate noise. You can carry your results forward buffer by buffer so it is extremely efficient. While you are doing the smoothing, you can also check for zero crossings. The signal need not be a "nearly perfect sine curve" for zero crossings to give you a good estimate of the frequency. After heavy smoothing, the zero crossing should correspond to your fundamental frequency.

If your signal is sufficiently sine like, you can calculate the frequency quite accurately by using the three point (simplest case) formula with the appropriate spacing that you can find here: Exact Near Instantaneous Frequency Formulas Best at Peaks (Part 1). You will want to select the center point at the peak (or trough) and the other two points at the nearest zero crossings. This formula will give you as good results as if you did interpolation at the zero crossings and calculated the frequency from there. You can then calculate the frequency at every peak and trough, that is every half cycle, and generate a new sequence of values. You can then apply a best polynomial fit on intervals of this sequence to give you very good calculations of the frequency and how it is varying.

A more accurate, but slightly more computationally expensive approach would be to calculate two bins of the DFT and use the formula found in my article: Two Bin Exact Frequency Formulas for a Pure Real Tone in a DFT. For the best result the frequency should fall halfway between the bins so you want to select a frame interval of a whole number of cycles and a half. Two and a half would be a good compromise between getting an accurate reading with a varying frequency and getting a better reading with multiple cycles. Point spacing is key for saving calculations. Four points per cycle will give the best results, so you want to be close to that. For two and a half cycles, that means ten points, which is not many. Since you are calculating the bin values separately, there is no advantage to a power of two. You can build lookup tables for the ten points for the two bin and three bin ahead of time to make it even more efficient. Pick a uniform spacing parameter that gets you as close to two and a half cycles as you can. Your calculations result in a cycles per frame value, then use your spacing value to convert to cycles per whatever your measure is.

Feel free to add any followup questions in your question and I will respond with an edit on my own.

---

I'm glad to hear you find my blog articles useful.

"How do you "heavy smooth" second harmonic when the base frequency could be 100x higher? "

Exponential smoothing has the property that the higher the frequency the more the tone will be attenuated. If you are still getting multiple zero crossing then you still have higher harmonics present with significant amplitude. Increase your smoothing parameter even more. If you are working in floating point, you don't have to worry about destroying your fundamental frequency. It will always be less dampened than the higher frequencies.

Since you have such a high sampling rate, you may also want to downsample your signal some as well. There are many ways to do this, but since this is not an area of expertise for me, I will hold off on recommendations. I can say the lowest calculation method is to just sample every nth point, but this does nothing to mitigate alias frequencies. Another simple method is average intervals to get a new sequence.

The nice thing about picking a whole plus a half number of cycles for your frame is that the second harmonic (twice the fundamental frequency), and all even harmonics above it, are a whole (or very close to) number of cycles in the frame so they will be orthogonal to your basis vectors. If there are a lot of harmonics present, you may want to increase the number of points above four per cycle to push the Nyquist frequency up for your frame length.

You can find my email address on my profile page. Feel free to contact me there.

P.S. Thanks to the upvoters.

Answer 2

So, since 1 MS/s is more than plenty than observe a 1 kHz signal, we'll stick with that. It's already highly oversampling!

The frequency may rise from zero to 1 kHz in a few seconds.

"A few seconds" is a couple million samples; this is an extremely slow change.

The idea is to calculate the result approximately every 10-20 periods of the signal.

So you're aiming for a robust rather than a very fast estimator; 10 periods is "pretty much", depending on your estimator.

The rest of your question is about FFT, which isn't the frequency estimator of choice here. You've got a very parametric signal model – a single tone that rises in frequency – so a parametric estimator will work better than an unparametric one. Furthermore, in your vastly oversampling approach, a DFT is a bad choice because your frequency resolution is fixed to sampling_rate/len_fft, and that's bad news for someone observing changing frequencies. So, don't use the DFT (FFT implementation of that).

There's a lot of parametric estimators – you're looking for a spectral estimator for line spectra; among these, you actually want a tone estimator; so, either something extremely simple like the $f_{est.}=\arctan\left(\frac{x[n]}{x[n_1]}\right)$ with $x$ being the analytic signal to your observed real-valued signal, or just a complex baseband version of that; you can use that as an input to e.g. a Kalman filter if you want to use the knowledge you have about the frequency trajectory.

A simple PLL / costas loop with a frequency error output would do, as well. In fact, there's a lot of beauty in the simplicity of that.

If you've got more than one tone, you could go for things like ESPRIT, which will give you an arbitrary fine frequency resolution with good variance properties. (ESPRIT could lend itself very well to detecting frequency change, but it would require modification in the rotational step)

Answer 3

Estimating "one strong frequency component" should work well using LPC analysis.