Add frequency estimators by Jacobsen and Quinn #503
Conversation
mbaz
commented
Jun 28, 2023
mbaz
commented
This commit adds two new frequency estimation algorithms:
* `jacobsen`, a very fast algorithm that performs a three-point
curve fit around the DFT peak. Not super accurate.
* `quinn`, a very accurate iterative algorithm based on ARMA(2,2).
It requires an initial frequency estimate (which may be
provided by `jacobsen`).
Both algorithms work with real and complex signals.
The `quinn` functions now only return a tuple with two values; the first is the estimated frequency, and the other is a boolean which is true if the algorithm ran for the maximum allowed number of iterations. The `converged` return value is somewhat redundant (and unavailable in the complex case).
|
Is this good to merge? |
Codecov ReportAll modified and coverable lines are covered by tests ✅
Additional details and impacted files@@ Coverage Diff @@
## master #503 +/- ##
==========================================
+ Coverage 97.58% 97.63% +0.05%
==========================================
Files 18 18
Lines 3147 3216 +69
==========================================
+ Hits 3071 3140 +69
Misses 76 76 ☔ View full report in Codecov by Sentry. |
|
Added a couple of tests to include the two lines that were missing coverage. |
|
Perhaps good to wait for a bit for @martinholters to take a look. |
Co-authored-by: Martin Holters <[email protected]>
Co-authored-by: Martin Holters <[email protected]>
Co-authored-by: Martin Holters <[email protected]>
Co-authored-by: Martin Holters <[email protected]>
Co-authored-by: Martin Holters <[email protected]>
Co-authored-by: Martin Holters <[email protected]>
Co-authored-by: Martin Holters <[email protected]>
|
@martinholters Thank you for the review! |
Before this commit, the `jacobsen` function did not perform any interpolation when the DFT peak occurred at frequencies 0 or Nyquist (`fs/2`). Using the symmetry of the DFT, the estimation formula can be both improved (by interpolating in all cases, with no exceptions), and simplified (since the DFT values around the peak are (often) complex conjugates). This is slightly tricky because, in order to save computation effort, we use the `rfft` instead of the `fft`. The formula can be simplified when the peak occurs at 0, and when it occurs at Nyquist and `N` is even. The result is overall improved error performance near zero and near Nyquist, even though in some frequency ranges the error increases.
|
@mbaz In preparation for a 0.8 release, would it be possible for you to get this PR ready? |
|
@ViralBShah I'm still actively working on it. Two things happened: (1) I realized applying Jacobsen to real signals requires a bit of care, and (2) my work load peaked the past two weeks. My plan at the moment is to have a single Jacobsen estimator based on the DFT. If/when I work out the issues with an RFFT-based estimator, I'll make another PR. How does that sound? What is the timeline for 0.8? I will try to finish up this PR in the next couple of days. |
|
Thank you. There is no tight deadline, and we can always add new features in new releases. I am just thinking out aloud, but if we can get as much of the backlog of PRs in as we can in the next 2 weeks or so - that will make for a nice 0.8 release. @martinholters I would of course defer to you on that. |
|
Agree on the 0.8 release. Would be nice to get of bunch of things done for it which are (and some of them have been for while) in an almost-ready state. OTOH, I don't want to release 0.8 too hastily in case we discover something else that requires a breaking change. My gut feeling is that 2 to 4 weeks from now might be a good time. |
|
@mbaz Just a gentle nudge to see if possible to finish this PR. 0.8 may be in a week or two. But of course, if you are busy, there's always 0.9! |
Jacobsen is sub-optimal for real signals. While the algorithm's complexity can be reduced, compared to the complex case, by taking advantage of the `rfft` transform, it is not really worth it given the large estimation error involved.
|
The test failure seems to be related to codecov. |
* Add explanation/warning about using Jacobsen with real signals to Jacobsen's docstring. * Also document that, if input is real, the frequency estimate will be positive. * Add tests for Jacobsen with real signals. The signal's frequency in these tests is close to Fs/4, where the estimation error is at its lowest.
When the frequency peak occurs at one of the spectrum extremes, use the fact that the DFT is periodic to obtain the correct spectrum values.
|
Thank you too, @martinholters -- the code is much better now than it was in the original PR. |
