# New frequency estimator from two DFT bins

Started by February 6, 2017
```>
>This subject is so basic and pure that I don't understand why it
>hasn't been treated in DSP litterature.
>
>Michael

Isn't that about the same thing that Galileo asked the Catholic Church?

Ced
---------------------------------------
Posted through http://www.DSPRelated.com
```
```On Wed, 08 Feb 2017 14:08:06 +0100, Michael Plet <me@home.com> wrote:

>On Tue, 07 Feb 2017 22:55:37 -0600, "martinvicanek" <99986@DSPRelated>
>wrote:
>
>>>Hi Group
>>>
>>>I have derived this estimator. The limited tests I have done shows
>>>high accuracy.
>>>
>>>
>>>Let k be the index of the DFT bin with the largest magnitude.
>>>Let j be the index of the DFT bin with the largest magnitude
>>>neigboring bin k.
>>>That is j=k-1 or j=k+1.
>>>
>>>Now let P=Tan(PI*k/N) and Q=Tan(PI*j/N)
>>>
>>>Then the normalized frequency is estimated by:
>>>
>>>Freq=(N/PI)*Arctan(Sqr(P*Q*(Im[k]*P-Im[j]*Q)/(Im[k]*Q-Im[j]*P)))
>>>
>>>where Sqr is Square root and Im[] is the imaginary part of X[].
>>>
>>>
>>>I was wondering if anyone has a test suite, so it would be possible to
>>>compare with other estimators?
>>>
>>>
>>>Regards,
>>>Michael Plet
>>
>>Hi Michael,
>>
>>your estimator seems to be equivalent to mine
>>(http://vicanek.de/articles/FreqFromTwoBins.pdf) if you choose a = b = 1.
>>I have uploaded an extended version of my writeup to facilitate
>>translation between my centered DFT notation and the standard DFT. Perhaps
>>you can get some inspiration from there for a formula based on Re[k] (and
>>some combination of formulas to optimize SNR)?
>>
>>Best
>>Martin
>>---------------------------------------
>>Posted through http://www.DSPRelated.com
>
>Hi Martin
>
>If a = b = 1, wouldn't my formula include Re[k]?
>
>I have tried to insert some of my test bin values into your formula.
>The result was not what I expected. I think I need to rotate the DFT
>bins before applying your formula. I will look into that.
>
>Thank you for the link to the updated document.
>I am sure I will find inspiration in that.
>
>This subject is so basic and pure that I don't understand why it
>hasn't been treated in DSP litterature.
>
>Michael

Hi Martin

Cedron ran a test showing that you are right about a = b = 1. Well
spotted Martin!
So my (independently discovered) formula is a special case of your
formula.

Michael
```
```On Wed, 08 Feb 2017 14:08:06 +0100, Michael Plet <me@home.com> wrote:

>On Tue, 07 Feb 2017 22:55:37 -0600, "martinvicanek" <99986@DSPRelated>
>wrote:
>
>>>Hi Group
>>>
>>>I have derived this estimator. The limited tests I have done shows
>>>high accuracy.
>>>
>>>
>>>Let k be the index of the DFT bin with the largest magnitude.
>>>Let j be the index of the DFT bin with the largest magnitude
>>>neigboring bin k.
>>>That is j=k-1 or j=k+1.
>>>
>>>Now let P=Tan(PI*k/N) and Q=Tan(PI*j/N)
>>>
>>>Then the normalized frequency is estimated by:
>>>
>>>Freq=(N/PI)*Arctan(Sqr(P*Q*(Im[k]*P-Im[j]*Q)/(Im[k]*Q-Im[j]*P)))
>>>
>>>where Sqr is Square root and Im[] is the imaginary part of X[].
>>>
>>>
>>>I was wondering if anyone has a test suite, so it would be possible to
>>>compare with other estimators?
>>>
>>>
>>>Regards,
>>>Michael Plet
>>
>>Hi Michael,
>>
>>your estimator seems to be equivalent to mine
>>(http://vicanek.de/articles/FreqFromTwoBins.pdf) if you choose a = b = 1.
>>I have uploaded an extended version of my writeup to facilitate
>>translation between my centered DFT notation and the standard DFT. Perhaps
>>you can get some inspiration from there for a formula based on Re[k] (and
>>some combination of formulas to optimize SNR)?
>>
>>Best
>>Martin
>>---------------------------------------
>>Posted through http://www.DSPRelated.com
>
>Hi Martin
>
>If a = b = 1, wouldn't my formula include Re[k]?
>
>I have tried to insert some of my test bin values into your formula.
>The result was not what I expected. I think I need to rotate the DFT
>bins before applying your formula. I will look into that.
>
>Thank you for the link to the updated document.
>I am sure I will find inspiration in that.
>
>This subject is so basic and pure that I don't understand why it
>hasn't been treated in DSP litterature.

There is actually a very large library of research on fine frequency
estimation, including a large library on interpolation of DFT
coefficients.   Three-term, two-term, and other techniques, like
dichotomous search, have really been beaten to death over the last
couple of decades, and much of the research goes back way before that.

That said, there are still a number of recent papers on refining
things for niche applications or corner cases, but the basic
optimizations don't seem to have changed much for a while.

---
This email has been checked for viruses by Avast antivirus software.
https://www.avast.com/antivirus

```
```On Mon, 06 Feb 2017 19:17:27 +0100, Michael Plet <me@home.com> wrote:

>Hi Group
>
>I have derived this estimator. The limited tests I have done shows
>high accuracy.
>
>
>Let k be the index of the DFT bin with the largest magnitude.
>Let j be the index of the DFT bin with the largest magnitude
>neigboring bin k.
>That is j=k-1 or j=k+1.
>
>Now let P=Tan(PI*k/N) and Q=Tan(PI*j/N)
>
>Then the normalized frequency is estimated by:
>
>Freq=(N/PI)*Arctan(Sqr(P*Q*(Im[k]*P-Im[j]*Q)/(Im[k]*Q-Im[j]*P)))
>
>where Sqr is Square root and Im[] is the imaginary part of X[].
>
>
>I was wondering if anyone has a test suite, so it would be possible to
>compare with other estimators?
>
>
>Regards,
>Michael Plet

Using trigonometrical identities I have managed to simplify my formula
(losing squares and square root).

The new formulation is:

Let k be the index of the DFT bin with the largest magnitude.
Let j be the index of the DFT bin with the largest magnitude
neigboring bin k.
That is j=k-1 or j=k+1.

Now let

P=Sin(2*PI*j/N), Q=Sin(2*PI*k/N), R=Cos(2*PI*j/N) and S=Cos(2*PI*k/N)

Then the normalized frequency is estimated by:

Freq=(N/(2*PI))*Arccos((Im[k]*P*S-Im[j]*Q*R)/(Im[k]*P-Im[j]*Q))

where Im[] is the imaginary part of X[].

Michael Plet
```
```On Thu, 09 Feb 2017 21:16:17 +0100, Michael Plet <me@home.com> wrote:

>On Mon, 06 Feb 2017 19:17:27 +0100, Michael Plet <me@home.com> wrote:
>
>>Hi Group
>>
>>I have derived this estimator. The limited tests I have done shows
>>high accuracy.
>>
>>
>>Let k be the index of the DFT bin with the largest magnitude.
>>Let j be the index of the DFT bin with the largest magnitude
>>neigboring bin k.
>>That is j=k-1 or j=k+1.
>>
>>Now let P=Tan(PI*k/N) and Q=Tan(PI*j/N)
>>
>>Then the normalized frequency is estimated by:
>>
>>Freq=(N/PI)*Arctan(Sqr(P*Q*(Im[k]*P-Im[j]*Q)/(Im[k]*Q-Im[j]*P)))
>>
>>where Sqr is Square root and Im[] is the imaginary part of X[].
>>
>>
>>I was wondering if anyone has a test suite, so it would be possible to
>>compare with other estimators?
>>
>>
>>Regards,
>>Michael Plet
>
>
>Using trigonometrical identities I have managed to simplify my formula
>(losing squares and square root).
>
>The new formulation is:
>
>
>Let k be the index of the DFT bin with the largest magnitude.
>Let j be the index of the DFT bin with the largest magnitude
>neigboring bin k.
>That is j=k-1 or j=k+1.
>
>Now let
>
>P=Sin(2*PI*j/N), Q=Sin(2*PI*k/N), R=Cos(2*PI*j/N) and S=Cos(2*PI*k/N)
>
>Then the normalized frequency is estimated by:
>
>Freq=(N/(2*PI))*Arccos((Im[k]*P*S-Im[j]*Q*R)/(Im[k]*P-Im[j]*Q))
>
>where Im[] is the imaginary part of X[].
>
>
>Michael Plet

Is there any advantage to this method?   It still seems very high in
computational complexity to me.

---
This email has been checked for viruses by Avast antivirus software.
https://www.avast.com/antivirus

```
```On Fri, 10 Feb 2017 15:08:40 GMT, eric.jacobsen@ieee.org wrote:

>On Thu, 09 Feb 2017 21:16:17 +0100, Michael Plet <me@home.com> wrote:
>
>>On Mon, 06 Feb 2017 19:17:27 +0100, Michael Plet <me@home.com> wrote:
>>
>>>Hi Group
>>>
>>>I have derived this estimator. The limited tests I have done shows
>>>high accuracy.
>>>
>>>
>>>Let k be the index of the DFT bin with the largest magnitude.
>>>Let j be the index of the DFT bin with the largest magnitude
>>>neigboring bin k.
>>>That is j=k-1 or j=k+1.
>>>
>>>Now let P=Tan(PI*k/N) and Q=Tan(PI*j/N)
>>>
>>>Then the normalized frequency is estimated by:
>>>
>>>Freq=(N/PI)*Arctan(Sqr(P*Q*(Im[k]*P-Im[j]*Q)/(Im[k]*Q-Im[j]*P)))
>>>
>>>where Sqr is Square root and Im[] is the imaginary part of X[].
>>>
>>>
>>>I was wondering if anyone has a test suite, so it would be possible to
>>>compare with other estimators?
>>>
>>>
>>>Regards,
>>>Michael Plet
>>
>>
>>Using trigonometrical identities I have managed to simplify my formula
>>(losing squares and square root).
>>
>>The new formulation is:
>>
>>
>>Let k be the index of the DFT bin with the largest magnitude.
>>Let j be the index of the DFT bin with the largest magnitude
>>neigboring bin k.
>>That is j=k-1 or j=k+1.
>>
>>Now let
>>
>>P=Sin(2*PI*j/N), Q=Sin(2*PI*k/N), R=Cos(2*PI*j/N) and S=Cos(2*PI*k/N)
>>
>>Then the normalized frequency is estimated by:
>>
>>Freq=(N/(2*PI))*Arccos((Im[k]*P*S-Im[j]*Q*R)/(Im[k]*P-Im[j]*Q))
>>
>>where Im[] is the imaginary part of X[].
>>
>>
>>Michael Plet
>
>Is there any advantage to this method?   It still seems very high in
>computational complexity to me.
>
>
>
>---
>This email has been checked for viruses by Avast antivirus software.
>https://www.avast.com/antivirus

The only advantage is that it gives the exact frequency in the
noiseless case and perform well in noisy cases.
```
```>>
>>Now let
>>
>>P=Sin(2*PI*j/N), Q=Sin(2*PI*k/N), R=Cos(2*PI*j/N) and S=Cos(2*PI*k/N)
>>
>>Then the normalized frequency is estimated by:
>>
>>Freq=(N/(2*PI))*Arccos((Im[k]*P*S-Im[j]*Q*R)/(Im[k]*P-Im[j]*Q))
>>
>>where Im[] is the imaginary part of X[].
>>
>>
>>Michael Plet
>
>Is there any advantage to this method?   It still seems very high in
>computational complexity to me.
>
>
>
>---
>This email has been checked for viruses by Avast antivirus software.
>https://www.avast.com/antivirus

Really?

Considering that P, S, Q, R, P*S, Q*R, and N/(2*PI) can all be
precomputed, the actual evaluation requires fewer computations, and is

In addition, since it is an exact solution, not an approximation which
does better with larger N, in any given application a DFT with a smaller N
can probably be used greatly reducing the number of overall calculations

Ced
---------------------------------------
Posted through http://www.DSPRelated.com
```
```On Fri, 10 Feb 2017 09:51:24 -0600, "Cedron" <103185@DSPRelated>
wrote:

>>>
>>>Now let
>>>
>>>P=Sin(2*PI*j/N), Q=Sin(2*PI*k/N), R=Cos(2*PI*j/N) and S=Cos(2*PI*k/N)
>>>
>>>Then the normalized frequency is estimated by:
>>>
>>>Freq=(N/(2*PI))*Arccos((Im[k]*P*S-Im[j]*Q*R)/(Im[k]*P-Im[j]*Q))
>>>
>>>where Im[] is the imaginary part of X[].
>>>
>>>
>>>Michael Plet
>>
>>Is there any advantage to this method?   It still seems very high in
>>computational complexity to me.
>>
>>
>>
>>---
>>This email has been checked for viruses by Avast antivirus software.
>>https://www.avast.com/antivirus
>
>Really?
>
>Considering that P, S, Q, R, P*S, Q*R, and N/(2*PI) can all be
>precomputed, the actual evaluation requires fewer computations, and is

My estimator isn't really the benchmark these days unless one is
really looking for low computational complexity or especially only
cares about low computational complexity for low SNR cases.   Candan's
estimator is still what most people use as a research reference, and
it's still pretty hard to beat in the broad sense.

The arccos is the computational hurlde, but clearly not
insurmountable.   Even in a LUT it may be undesirable in constrained
applications, though.  Requirements often get in the way.

>In addition, since it is an exact solution, not an approximation which
>does better with larger N, in any given application a DFT with a smaller N
>can probably be used greatly reducing the number of overall calculations

Remember that many people don't think "exact" is relevant to
interpolations made from estimates.  The real issues is just whether
the performance vs computational complexity tradeoff works for a
particular application space of interest.   No single estimator does
all things for all people.  I was just asking what Michael thinks the
advantage of this particular estimator might be.   Several existing
estimators these days are derived from closed-form derivations, so
that's not really a distinction.

---
This email has been checked for viruses by Avast antivirus software.
https://www.avast.com/antivirus

```
```>>
>>Really?
>>
>>Considering that P, S, Q, R, P*S, Q*R, and N/(2*PI) can all be
>>precomputed, the actual evaluation requires fewer computations, and is
>
>My estimator isn't really the benchmark these days unless one is
>really looking for low computational complexity or especially only
>cares about low computational complexity for low SNR cases.   Candan's
>estimator is still what most people use as a research reference, and
>it's still pretty hard to beat in the broad sense.
>

Your estimator isn't as low computational complexity as you imply.  Yours
(Candan's and mine as well) requires a complex division.  A complex
division takes 8 multiplies, 1 divide, and 3 sums.  Plet's (and my 2 bin
real) uses a real division and that's a big advantage.

If they want the best, people should be using my formula.  I refer you to
Figure 4 in Julien's comparison paper and the data below.  Mine beats
Candan's in every sense except number of computations where it is roughly
comparable.

For low SNR cases it is obvious that doing extra calculations for extra
precision is meaningless.

>The arccos is the computational hurlde, but clearly not
>insurmountable.   Even in a LUT it may be undesirable in constrained
>applications, though.  Requirements often get in the way.
>

If one doesn't mind a little approximation, which shouldn't matter in most
practical applications, a LUT solution for arccos is going to be faster
than a complex division.  I'm not sure if a acos call might not be faster
either, but I don't feel like testing it.

>>In addition, since it is an exact solution, not an approximation which
>>does better with larger N, in any given application a DFT with a smaller
N
>>can probably be used greatly reducing the number of overall
calculations
>
>Remember that many people don't think "exact" is relevant to
>interpolations made from estimates.  The real issues is just whether
>the performance vs computational complexity tradeoff works for a
>particular application space of interest.   No single estimator does
>all things for all people.  I was just asking what Michael thinks the
>advantage of this particular estimator might be.   Several existing
>estimators these days are derived from closed-form derivations, so
>that's not really a distinction.
>

I don't think of a frequency calculation as an interpolation.  That
concept comes from thinking of the discrete case as a dirac delta version
of the continuous case.  I know that how DSP is mostly taught, but it also
a source of many misconceptions, such as leakage being defined by the sinc
function.

Again, I want to make the point that if a little bit more calculation in
the frequency equation allows you to reduce the size of the DFT then a
relatively huge amount of calculations is saved.

Finally, closed form is not synonymous with exact.  These are the exact
solutions I am aware of:

Mine: 2 bin real, 3 bin real, 2 bin complex, 3 bin complex (actually I
have n bin version of both real and complex)

Martin Vicanek: 2 bin real

Michael Plet: two versions of 2 bin real

Candan 2013: 3 bin complex, but you would know it was exact from his
derivation.  It is mathematically equivalent to my 3 bin complex which
preceded his.

Do you know of any others that you can add to the list?

Ced

===========================================================

The sample count is 100
and the run size is 100
The phase value is  0.785398

Errors are shown at 100x actual value

Target Noise Level = 0.000

Freq Candan2013      CD3Bin          NewMP2Bin
---- -------------   -------------   -------------
3.0   0.000  0.000    0.000  0.000   -0.000  0.000
3.1  -0.165  0.000   -0.000  0.000   -0.000  0.000
3.2  -0.509  0.000    0.000  0.000   -0.000  0.000
3.3  -0.711  0.000    0.000  0.000   -0.000  0.000
3.4  -0.532  0.000   -0.000  0.000    0.000  0.000
3.5  -0.001  0.000    0.000  0.000   -0.000  0.000
3.6  -0.336  0.000   -0.000  0.000    0.000  0.000
3.7  -0.437  0.000    0.000  0.000   -0.000  0.000
3.8  -0.303  0.000    0.000  0.000   -0.000  0.000
3.9  -0.096  0.000    0.000  0.000    0.000  0.000

Target Noise Level = 0.001

Freq Candan2013      CD3Bin          NewMP2Bin
---- -------------   -------------   -------------
3.0   0.001  0.010    0.001  0.011    0.003  0.018
3.1  -0.166  0.010    0.000  0.010    0.000  0.012
3.2  -0.508  0.011    0.001  0.010   -0.000  0.008
3.3  -0.710  0.010    0.001  0.009    0.001  0.008
3.4  -0.534  0.013   -0.002  0.012   -0.001  0.009
3.5  -0.004  0.018   -0.002  0.015   -0.001  0.011
3.6  -0.337  0.013   -0.001  0.015    0.002  0.015
3.7  -0.438  0.011   -0.001  0.012   -0.004  0.066
3.8  -0.305  0.011   -0.001  0.012   -0.005  0.068
3.9  -0.095  0.011    0.000  0.011   -0.002  0.028

Target Noise Level = 0.010

Freq Candan2013      CD3Bin          NewMP2Bin
---- -------------   -------------   -------------
3.0   0.005  0.105    0.004  0.109   -0.002  0.183
3.1  -0.173  0.103   -0.009  0.098   -0.007  0.104
3.2  -0.493  0.117    0.015  0.109    0.008  0.098
3.3  -0.708  0.120    0.003  0.106    0.004  0.084
3.4  -0.543  0.126   -0.011  0.110   -0.005  0.092
3.5  -0.023  0.166   -0.019  0.144   -0.014  0.113
3.6  -0.323  0.134    0.014  0.152   -0.020  0.173
3.7  -0.421  0.124    0.019  0.135    0.051  0.601
3.8  -0.302  0.118    0.002  0.128    0.027  0.733
3.9  -0.093  0.103    0.005  0.106   -0.050  0.336

Target Noise Level = 0.100

Freq Candan2013      CD3Bin          NewMP2Bin
---- -------------   -------------   -------------
3.0   0.218  1.009    0.227  1.017    0.178  1.715
3.1  -0.047  1.161    0.093  1.123    0.027  1.263
3.2  -0.584  1.177   -0.089  1.074   -0.127  0.894
3.3  -0.649  1.288    0.035  1.145   -0.048  0.801
3.4  -0.554  1.441   -0.048  1.272   -0.067  0.976
3.5  -0.187  1.769   -0.137  1.506    0.021  1.160
3.6  -0.241  1.190    0.105  1.344   -0.233  1.755
3.7  -0.445  1.252   -0.031  1.377   -0.467  6.959
3.8  -0.314  1.063   -0.019  1.159    1.331  7.047
3.9  -0.104  0.943   -0.026  0.998    0.657  2.764

---------------------------------------
Posted through http://www.DSPRelated.com
```
```>
>  Candan 2013: 3 bin complex, but you would know it was exact from his
>derivation.

I should proof read more carefully. That should say "wouldn't know it".

Ced
---------------------------------------
Posted through http://www.DSPRelated.com
```