> >This subject is so basic and pure that I don't understand why it >hasn't been treated in DSP litterature. > >MichaelIsn't that about the same thing that Galileo asked the Catholic Church? Ced --------------------------------------- Posted through http://www.DSPRelated.com
New frequency estimator from two DFT bins
Started by ●February 6, 2017
Reply by ●February 8, 20172017-02-08
Reply by ●February 8, 20172017-02-08
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. > >MichaelHi 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
Reply by ●February 8, 20172017-02-08
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
Reply by ●February 9, 20172017-02-09
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 PletUsing 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
Reply by ●February 10, 20172017-02-10
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 PletIs 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
Reply by ●February 10, 20172017-02-10
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/antivirusThe only advantage is that it gives the exact frequency in the noiseless case and perform well in noisy cases.
Reply by ●February 10, 20172017-02-10
>> >>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/antivirusReally? 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 less complicated, than your estimator. 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 that have to be made. Ced --------------------------------------- Posted through http://www.DSPRelated.com
Reply by ●February 10, 20172017-02-10
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 >less complicated, than your estimator.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 >that have to be made.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
Reply by ●February 10, 20172017-02-10
>> >>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 >>less complicated, than your estimator. > >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 smallerN>>can probably be used greatly reducing the number of overallcalculations>>that have to be made. > >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
Reply by ●February 10, 20172017-02-10
> > 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
