Candan's Tweaks of Jacobsen's Frequency Approximation
Introduction
This is an article to hopefully give a better understanding of the Discrete Fourier Transform (DFT) by explaining how a tweak to a well known frequency approximation formula makes it better, and another tweak makes it exact. The first tweak is shown to be the first of a pattern and a novel approximation formula is made from the second. It only requires a few extra calculations beyond the original approximation to come up with an approximation suitable for most uses.
Angles from Bin Numbers
The DFT bins can be thought of as being placed evenly around the unit circle on the complex plane. The number of bins is designated as $N$, and the bins are index by the variable $k$. This creates a scale around the circumference where each unit multiple is a bin location.
The corresponding angle to the bin is called $\beta_k$. On the angle scale, the arc distance between two bins is $\beta_1$.
In a single pure complex tone, the frequency is determined by the angle step size $\alpha$ for each sample. When a DFT signal frame is defined, the tone will have a frequency of $f$ on the bin index scale.
In this way, the complex unit circle acts like a frequency dial with two scales on it. The bins are the marks and the frequency is the needle. The bin index units of $k$ and $f$ are cycles per frame, and the angle units are either radians per bin or radians per sample. $$ \begin{aligned} \alpha &= f \cdot \frac{ 2\pi }{ N } \\ \beta_k &= k \cdot \frac{ 2\pi }{ N } \\ \beta_1 &= \frac{ 2\pi }{ N } \\ \end{aligned} \tag {1} $$
Bin Values at a Bin
A complex pure tone signal having a frequency of $f$, an amplitude of $M$, and a phase of $\phi$ can be defined by: $$ S_n = M e^{ i ( \alpha n + \phi ) } = M e^{ i \phi } \left[e^{ i \alpha} \right]^n = M e^{ i \phi } \left[e^{ i \frac{ 2\pi }{ N } f} \right]^n = M e^{ i \phi } e^{ i 2\pi \frac{fn}{N} } \tag {2} $$ This makes it clear that the units of $\alpha$ are radians per sample and $\phi$ is a radian offset.
The formula for determining the bin value for a pure complex tone in a $1/N$ normalized DFT having a whole number of cycles per frame is: $$ Z_f = M e^{ i \phi } \tag {3} $$ All the other bin values are zero.
The formula for determining the bin value for a pure complex tone that isn't a whole number of cycles in a frame is: $$ Z_k = \frac{ M }{ N } e^{ i \phi } \cdot \frac{ 1 - e^{ i\alpha N } }{ 1 - e^{ i(\alpha - \beta_k ) } } \tag {4} $$ Here $k$ represents an integer spacing on the bin index scale, but is can also be considered as a fractional value for interpolation. However, these in between values do not match the values one would get from the DFT itself due to a simplifying assumption ($e^{i \beta_k N }= e^{i 2 \pi k} = 1 $ ) in the derivation [1] which is obviously not true for fractional values.
None of the bin values are zero.
Twist Factors
Any value on the complex unit circle acts as a twist factor when used in multiplication. The angle of the twisting is the angle of the value. Using the twist factors as proxies for the angles makes a lot of the following math simpler to follow. $$ \begin{aligned} a &= e^{ i(\alpha - \beta_w ) } = e^{ i \frac{ 2\pi }{ N } ( f - w ) } \\ b &= e^{ -i \beta_1 } = e^{ -i \frac{ 2\pi }{ N } } \\ \end{aligned} \tag {5} $$ The index $w$ is the closest whole number to the frequency $f$, so $a$ as a twist factor is expected to be close to $1+i0$, aka unity. The $b$ value is the back twist multiplier representing a backward rotation of one bin index arc.
Using these variables makes the expression of other bin values relative to a base bin simpler in the bin value formula (4). $$ Z_{w+d} = \frac{ M }{ N } e^{ i \phi } \cdot \frac{ 1 - e^{ i\alpha N } }{ 1 - a b^d } \tag {6} $$
Jacobsen's Frequency Approximation Formula
A good approximation formula for the frequency calculates the fraction offset from the peak bin index on the bin index scale [2]. A quotient is calculated from bin values. The real part of the quotient is then added to the peak bin index to get the approximate frequency value. $$ \begin{aligned} Q_J &= \frac{ Z_{w-1} - Z_{w+1} }{ - Z_{w-1} + 2 Z_{w} - Z_{w+1} } \\ f &\approx w + Re[Q_J] \\ \end{aligned} \tag {7} $$ It is a straightforward calculation that yields very close results with either larger N or frequencies near bin values. In many applications, the approximation error will be swamped by noise error or other interference and the extra calculations needed for an apparent exact value are likely not worth doing.
Why it Works
Since the DFT is a linear transform, multiplying the input signal by any scalar will simply result in the DFT spectrum also being multiplied by the same scalar. Both the phase value ($e^{i\phi}$), and the amplitude ($M$) parameters are factors on the signal (2), so they are also factors of the DFT spectrum (4). Any linear combination of bin values will also have them as a common factor. Therefore a quotient of linear combination sums will cancel those factors and they will be washed out leaving only frequency related values.
Plug (6) into (7) and crunch the proxies. $$ \begin{aligned} Q_J &= \frac{ \frac{\frac{ M }{ N } e^{ i \phi } \left( 1 - e^{ i\alpha N } \right)}{1 - ab^{-1} } - \frac{\frac{ M }{ N } e^{ i \phi } \left( 1 - e^{ i\alpha N } \right)}{1 - ab } } { - \frac{\frac{ M }{ N } e^{ i \phi } \left( 1 - e^{ i\alpha N } \right)}{1 - ab^{-1} } + \frac{2\frac{ M }{ N } e^{ i \phi } \left( 1 - e^{ i\alpha N } \right)}{1 - a } - \frac{\frac{ M }{ N } e^{ i \phi } \left( 1 - e^{ i\alpha N } \right)}{1 - ab } } \\ &= \frac{ \frac{1}{1 - ab^{-1} } - \frac{1}{1 - ab } } { - \frac{1}{1 - ab^{-1} } + \frac{2}{1 - a } - \frac{1}{1 - ab } } \\ &= \frac{ (1-a) \left[ (1 - ab) - (1 - ab^{-1}) \right] } { - \left[ (1 - a ) (1 - ab) \right] + 2 \left[ (1 - ab^{-1}) (1 - ab) \right] - \left[ (1 - ab^{-1}) (1 - a ) \right] } \\ &= \frac{ (1-a) a \left[ -b + b^{-1} \right] } { - \left[ 1 - a - ab + a^2b \right] + 2 \left[ 1 - ab^{-1} - ab + a^2 \right] - \left[ 1 - ab^{-1} - a + a^2b^{-1} \right] } \\ &= \frac{ (1-a) \left[ -b + b^{-1} \right] } { \left[ 1 + b - ab \right] + \left[ - 2b^{-1} - 2b + 2a \right] + \left[ b^{-1} + 1 - ab^{-1} \right] } \\ &= \frac{ (1-a) \left[ -b + b^{-1} \right] } { \left[ 2 - b - b^{-1} \right] + a \left[ -b + 2 - b^{-1} \right] } \\ &= \frac{ (1-a) }{ (1+a) } \cdot \frac{ \left[ -b + b^{-1} \right] }{ \left[ -b + 2 - b^{-1} \right] } \\ \end{aligned} \tag {8} $$ The final line shows the unknown frequency proxy ($a$) separated from the known bin arc size proxy ($b$). Since $Q_J$ is known, an equation for $a$ can be derived directly which is done in a following section. A slicker approach uses a return to the angle scale from the twist factors and a little algebra to derive a trigonometric expression. $$ \begin{aligned} Q_J &= \frac{ \left(1-e^{ i(\alpha - \beta_w ) } \right) } { \left(1+e^{ i(\alpha - \beta_w ) } \right) } \cdot \frac{ \left[ -e^{ -i \beta_1 } + e^{ i \beta_1 } \right] } { \left[ -e^{ -i \beta_1 } + 2 - e^{ i \beta_1 } \right] } \\ &= \frac { e^{ -i(\alpha - \beta_w )/2 } - e^{ i(\alpha - \beta_w )/2 } } { e^{ -i(\alpha - \beta_w )/2 } + e^{ i(\alpha - \beta_w )/2 } } \cdot \frac { \left( e^{ i \beta_1 / 2 } - e^{ -i \beta_1 / 2 } \right) \left( e^{ i \beta_1 / 2 } + e^{ -i \beta_1 / 2 } \right) } { -\left( e^{ i \beta_1 / 2 } - e^{ -i \beta_1 / 2 } \right)^2 } \\ &= \frac{ -2i\sin\left(\frac{\alpha - \beta_w}{2}\right) } {2\cos\left(\frac{\alpha - \beta_w}{2}\right)} \cdot \frac{2\cos\left(\frac{\beta_1}{2}\right)} {-2i\sin\left(\frac{\beta_1}{2}\right)} \\ &= \tan\left(\frac{\alpha - \beta_w}{2}\right) \cdot \frac{1}{\tan\left(\frac{\beta_1}{2}\right)} \\ \end{aligned} \tag {9} $$ From this form it is clear that the value of $Q_J$ is expected to be real. When calculated from actual data there may be an imaginary component. This should be discarded.
Cross multiplying and taking the arctan of both sides brings the equation to the angle scale. $$ \frac{ \alpha - \beta_w }{2} = \arctan \left( Re[Q_J] \cdot \tan\left(\frac{\beta_1}{2}\right) \right) \tag {10} $$ The value of $\alpha$ can now be calculated. $$ \alpha = \beta_w + 2 \arctan \left( Re[Q_J] \cdot \tan\left(\frac{\beta_1}{2}\right) \right) \tag {11} $$ Multiplying by the number of bins per radian ($N/(2\pi)$) converts the equation to the bin index scale. $$ f = w + \arctan \left( Re[Q_J] \cdot \tan\left(\frac{\pi}{N}\right) \right) \cdot \frac{N}{\pi} \tag {12} $$ This equation is exact since no approximations were taken.
For small values of $\theta$, $\tan(\theta) \approx \theta $ and for small values of $m$, $\arctan(m) \approx m $. The smaller, the closer. Therefore, as $N$ gets larger, the $\arctan$ and $\tan$ functions approach identity and act almost like factors of one. The $\pi$ and $N$ values then cancel out and Jacobsen's approximation (7) is the result. $$ \lim_{N \to \infty} \arctan \left( Re[Q_J] \cdot \tan\left(\frac{\pi}{N}\right) \right) \cdot \frac{N}{\pi} = Re[Q_J] \tag {13} $$ Jacobsen's formula is the limit case as $N$ goes to infinity. The values in Appendix II show that N doesn't have to get very large in order for the approximation error to get very small.
Comparison to Candan's Formulas
Candan's 2011 formula [3] is very similar. It is the first tweak of Jacobsen's quotient (7). There are approximations taken in its derivation. $$ \begin{aligned} Q_{C1} &= \frac{ \tan(\pi/N) }{ \pi/N } \cdot Re[Q_J] \\ f &\approx w + Q_{C1} \\ \end{aligned} \tag {14} $$ Since $Q_J$ is expected to be real, the 'Re' is inert in the pure tone case.
The final result has a constant multiplier times the Jacobsen offset. Comparing it to (12) shows the employment of the $\arctan{m} \approx m$ approximation for small $m$.
Candan's 2013 updated formula [4] is based on the 2011 one. A tweak of a tweak, so to speak. $$ \begin{aligned} Q_{C2} &= \frac{\arctan \left( Q_{C1} \cdot \frac{\pi}{N} \right) }{ \frac{\pi}{N} } \\ f &= w + Q_{C2} \\ \end{aligned} \tag {15} $$ The tweaks can be combined into a single equation. $$ \begin{aligned} f &= w + \arctan \left( \frac{ \tan\left(\frac{\pi}{N}\right) } { \frac{\pi}{N} } \cdot Re[Q_J] \cdot \frac{\pi}{N} \right) \cdot \frac{N}{\pi} \\ &= w + \arctan \left( Re[Q_J] \cdot \tan \left(\frac{\pi}{N} \right) \right) \cdot \frac{N}{\pi} \\ \end{aligned} \tag {16} $$ This one matches (12) exactly, although it is derived quite differently.
ArcTangent of a Small Angle
The power series solution for the ArcTangent function is simple to express as an infinite polynomial. $$ \begin{aligned} \arctan \left( m \right) &=\sum_{t=0} (-1)^t \frac{m^{(2t+1)}}{2t+1}\\ &= m - \frac{m^3}{3} + \frac{m^5}{5} - \frac{m^7}{7} + \frac{m^9}{9} - ...\\ &\approx m \left( 1 - m^2 \left( \frac{1}{3} - m^2 \frac{1}{5} \right) \right) \end{aligned} \tag {17} $$ The last line is a truncated easy to calculate version. The series expansion clearly shows why $\arctan(m) \approx m$ for small values of $m$.
A Closer Approximation
By incorporating the first two terms of the arctan approximation (17) into the exact frequency formula (12), a new approximation formula can be derived. It is one step up from Candan's 2011 formula. $$ \begin{aligned} f &\approx w + \left[ \left( Q_J \cdot \tan\left(\frac{\pi}{N}\right) \right) - \frac{ \left( Q_J \cdot \tan\left(\frac{\pi}{N}\right) \right)^3 }{3} \right] \cdot \frac{ N } {\pi} \\ &= w + Q_J \cdot \tan\left(\frac{\pi}{N}\right) \cdot \frac{ N } {\pi} \left[ 1 - \frac{ \left( Q_J \cdot \tan\left(\frac{\pi}{N}\right) \right)^2 }{3} \right] \\ &= w + Q_J \left[ K_1 - Q_J^2 K_2 \right] \end{aligned} \tag {18} $$ Where $K_1$ and $K_1$ get defined as: $$ \begin{aligned} K_1 &= \tan\left(\frac{\pi}{N}\right) \cdot \frac{ N } {\pi} \approx 1 + \frac{1}{3} \frac{\pi^2}{N^2}\\ K_2 &= \tan^2\left(\frac{\pi}{N}\right) \cdot \frac{ K_1 } {3} \approx \frac{1}{3} \frac{\pi^2}{N^2}\\ \end{aligned} \tag {19} $$ This adds just a few calculations but the results should be just as good as using (12) for almost all applications.
The values of $K_1$ and $K_2$ are calculated ahead of time and only needs to be done once, so there is no material advantage to using a series over a $\tan$ call if it is available.
Of course, another term or more could be used to make a closer approximation. This isn't practically necessary though.
Alternative Exact Equation
Jacobsen's quotient (7) can be expressed in terms of vector products. $$ Q_J =\frac{ \begin{bmatrix} 1 & 0 & -1 \\ \end{bmatrix} \vec Z }{ \begin{bmatrix} -1 & 2 & -1 \\ \end{bmatrix} \vec Z } \tag {20} $$ It can also be expressed using the value of $a$ and $b$ using (8). The frequency information can also be represented by another proxy $R$. $$ Q_J = \frac{ (1-a) } { (1+a) } \cdot \frac{ \left( b - b^{-1} \right) } { \left( b - 2 + b^{-1} \right) } = R \cdot \frac{ \left( b - b^{-1} \right) } { \left( b - 2 + b^{-1} \right) } \tag {21} $$ Put these together and solve for R. $$ R = \frac{ (1-a) }{(1+a) } = \frac { \left( b - 2 + b^{-1} \right) } { \left( b - b^{-1} \right) } \cdot \frac{ \begin{bmatrix} 1 & 0 & -1 \\ \end{bmatrix} \vec Z }{ \begin{bmatrix} -1 & 2 & -1 \\ \end{bmatrix} \vec Z } \tag {22} $$ The value of $R$ can be calculated from the right hand part of the equation. Then it can be used to find the values of $a$ from the middle part. First, cross multiply. $$ R+Ra = 1-a \tag {23} $$ Second, isolate $a$ on the left side. Then replace $R$ using the right hand known value from (22). $$ \begin{aligned} a &= \frac{ 1-R }{1+R } \\ &= \frac{ \left( \left( b - b^{-1} \right) \begin{bmatrix} -1 & 2 & -1 \\ \end{bmatrix} -\left( b - 2 + b^{-1} \right) \begin{bmatrix} 1 & 0 & -1 \\ \end{bmatrix} \right) \vec Z }{ \left( \left( b - b^{-1} \right) \begin{bmatrix} -1 & 2 & -1 \\ \end{bmatrix} +\left( b - 2 + b^{-1} \right) \begin{bmatrix} 1 & 0 & -1 \\ \end{bmatrix} \right) \vec Z } \\ &= \frac{ \begin{bmatrix} -2b+2 & 2b - 2b^{-1} & -2 + 2b^{-1} \\ \end{bmatrix} \cdot \vec Z }{ \begin{bmatrix} -2 + 2b^{-1} & 2b - 2b^{-1} & -2b + 2 \\ \end{bmatrix} \cdot \vec Z } \\ &= \frac{ \begin{bmatrix} -(b-1) & (b^2 - 1 ) b^{-1} & -(b-1)b^{-1} \\ \end{bmatrix} \cdot \vec Z }{ \begin{bmatrix} -(b-1)b^{-1} & (b^2 - 1 ) b^{-1} & -(b-1) \\ \end{bmatrix} \cdot \vec Z } \\ &= \frac{ \begin{bmatrix} -b & b + 1 & -1 \\ \end{bmatrix} \cdot \vec Z }{ \begin{bmatrix} -1 & b + 1 & -b \\ \end{bmatrix} \cdot \vec Z } \\ \end{aligned} \tag {24} $$ This is the same equation that is called "The [1,1,1] Squasher" in my previous three bin complex tone frequency formula article [5]. Both the numerator and denominator vectors will zero out any $\begin{bmatrix} 1 & 1 & 1 \end{bmatrix}$ component in $\vec Z$. This is also true for $Q_J$.
The dot products can be multiplied out. $$ \begin{aligned} a &= \frac{ \begin{bmatrix} -b & b + 1 & -1 \\ \end{bmatrix} \cdot \begin{bmatrix} Z_{w-1} \\ Z_{w} \\ Z_{w+1} \\ \end{bmatrix} }{ \begin{bmatrix} -1 & b + 1 & -b \\ \end{bmatrix} \cdot \begin{bmatrix} Z_{w-1} \\ Z_{w} \\ Z_{w+1} \\ \end{bmatrix} } \\ &= \frac{ - b Z_{w-1} + (b+1) Z_{w} - Z_{w+1} }{ - Z_{w-1} + (b+1) Z_{w} - b Z_{w+1} } \\ \end{aligned} \tag {25} $$ This version is my original complex formula submitted to Julien Arzi for his comparison report [6]. This report also includes Jacobsen's and Candan's formulas.
Recovering the value of $f$ from $a$ means needing to know the angle value and converting it back to the bin index scale. The angle is found by using the $\arg$ function. $$ \begin{aligned} \arg( a ) &= \alpha - \beta_k \\ \alpha &= \beta_k + \arg( a ) \\ \end{aligned} \tag {26} $$ In programming, the angle is found by taking the inverse tangent. In order to cover all four quadrants, the $\operatorname{atan2}$ function is usually used rather the regular $\operatorname{atan}$ function. However, since $a$ is near unity, this is not a concern. $$ \begin{aligned} f &= k + \arg( a ) \frac{ N }{ 2\pi }\\ &= k + \operatorname{atan2} \left( Im[a], Re[a] \right) \cdot \frac{ N }{ 2\pi } \\ &= k + \arctan \left( \frac{Im[a]}{Re[a]} \right) \cdot \frac{ N }{ 2\pi } \\ \end{aligned} \tag {27} $$ This formula is not included in the comparison in the appendices since it is not directly based on Jacobsen's quotient (7) and its results match the exact value.
This approach takes a lot more calculations than (12).
The returned arc length is expected to be less than half a bin arc width. For most applications, it won't take many terms to make the series answer close enough to the exact answer. This can save having to make a trig function call.
The Squasher vs. Candan 2013
This equation and Candan's tweak of a tweak are both derived from (8) directly without approximation, and are therefore mathematically equivalent for pure tones. In practice though, because of how extraneous misfitting information is discarded, they will differ very slightly in behavior. The difference occurs in how skew from error and perhaps other tones is discarded. In (12), the imaginary part of Jacobsen's quotient (7) is discarded. In (27) skewness is discarded by ignoring the magnitude of $a$.
In each case, the size of the discard is an indicator of fitness, but can usually be ignored. A zero discard value cannot be interpreted as an exact value having been found, it might just be good luck.
Conclusion
Here are the frequency equations based on Jacobsen's Quotient (7). The first one is exact, and the rest are ordered by strictly decreasing accuracy. $$ \begin{aligned} f &= w + \arctan \left( Re[Q_J] \cdot \tan(\frac{\pi}{N}) \right) \cdot \frac{N}{\pi} \\ &\approx w + Q_J \left[ K_1 - Q_J^2 K_2 \right] \\ &\approx w + Q_J K_1 \\ &\approx w + Q_J \\ \end{aligned} \tag {28} $$ The distinctions between these formulas become negligible as the number of samples in the signal frame, also being the number of bins in the DFT, reaches less than a typical usage size. The distinctions become important with small frames and high SNR signals and the need for highly accurate results.
The other exact formula, my original, got a fresh (but very similar) derivation, and is still more complicated to calculate.
The approximation formula derived in this article is novel to me.
Appendix I - Comparison Source Code
This code is written to demonstrate the formulas, not for speed.
#include <math.h>
#include <stdio.h>
//============================================================================
void BuildPatch( int N, double f, double phi, double M,
double S[128][2], double Z[3][2] )
{
//--- Find parameters
int w = (int) f;
double RadiansPerBin = 2.0 * M_PI / (double) N;
double alpha = f * RadiansPerBin;
//--- Build the signal
for( int n = 0; n < N; n++ )
{
double angle = (double) n * alpha + phi;
S[n][0] = M * cos( angle );
S[n][1] = M * sin( angle );
}
//--- Take the partial unnormalized DFT (index scale comments)
for( int d = -1; d <=1; d++ )
{
int k = w + d;
double angle = 0.0; // int kn = 0;
double arc = k * RadiansPerBin;
double sum_real = 0.0;
double sum_imag = 0.0;
for( int n = 0; n < N; n++ )
{
double cos_beta_kn = cos( angle ); // = beta_real[kn];
double sin_beta_kn = sin( angle ); // = beta_imag[kn];
sum_real += S[n][0] * cos_beta_kn + S[n][1] * sin_beta_kn;
sum_imag += S[n][1] * cos_beta_kn - S[n][0] * sin_beta_kn;
angle += arc; // kn += k; if( kn >= N ) kn -= N;
}
Z[d+1][0] = sum_real;
Z[d+1][1] = sum_imag;
}
}
//============================================================================
double CalculateJacobsenQuotient( double Z[3][2] )
{
double num_real = Z[0][0] - Z[2][0];
double num_imag = Z[0][1] - Z[2][1];
double den_real = -Z[0][0] + 2.0 * Z[1][0]- Z[2][0];
double den_imag = -Z[0][1] + 2.0 * Z[1][1]- Z[2][1];
double den2 = den_real * den_real + den_imag * den_imag;
double quotient_real = ( num_real * den_real + num_imag * den_imag )
/ den2;
// double quotient_imag = ( num_imag * den_real - num_real * den_imag )
// / den2;
return quotient_real;
}
//============================================================================
void PrintLine( double p, double QJ,
double K1, double K2,
double PiOverN, double Tan_PiOverN, double Z[3][2] )
{
double Candan = QJ * K1;
double Closer = QJ * ( K1 - QJ * QJ * K2 );
double Exact = atan( QJ * Tan_PiOverN ) / PiOverN;
printf( "%3.1f %8.5f%8.5f%8.5f%8.5f ",
p, QJ, Candan, Closer, Exact );
double normalize = PiOverN / M_PI; // aka 1/N
double normalize2 = normalize * normalize;
double sum2 = 0.0;
for( int i = 0; i < 3; i++ )
{
double Z2 = Z[i][0] * Z[i][0] + Z[i][1] * Z[i][1];
double v2 = Z2 * normalize2;
printf( " %7.5f", sqrt( v2 ) );
sum2 += v2;
}
double portion = sqrt( sum2 );
printf( " %6.3f\n", portion );
}
//============================================================================
int main( int ArgCount, char* ArgValues[] )
{
double S[128][2];
double Z[3][2];
int N = 8;
for( int l = 0; l < 3; l++ )
{
double PiOverN = M_PI / (double) N;
double Tan_PiOverN = tan( PiOverN );
double K1 = Tan_PiOverN / PiOverN;
double K2 = Tan_PiOverN * Tan_PiOverN * K1 / 3.0;
printf( "\n N = %d\n\n"
" p QJ Can2011 Closer Exact |z[-1]| |z[0]| |z[+1]| cover\n"
"--- ------- ------- ------- ------- ------- ------- ------- -----\n", N );
for( int h = 0; h <= 50; h += 10 )
{
double p = (double) h * 0.01;
double f = 3.0 + p;
BuildPatch( N, f, 0.123, 1.0, S, Z );
double QJ = CalculateJacobsenQuotient( Z );
PrintLine( p, QJ, K1, K2, PiOverN, Tan_PiOverN, Z );
}
printf( "\n" );
N *= 4;
}
//--- Exit
return 0;
}
//============================================================================
Appendix II - Comparison Results
Here is the output.
N = 8 p QJ Can2011 Closer Exact |z[-1]| |z[0]| |z[+1]| cover --- ------- ------- ------- ------- ------- ------- ------- ----- 0.0 -0.00000-0.00000-0.00000-0.00000 0.00000 1.00000 0.00000 1.000 0.1 0.09485 0.10005 0.10000 0.10000 0.09226 0.98388 0.11160 0.994 0.2 0.19000 0.20041 0.20000 0.20000 0.16184 0.93645 0.23776 0.980 0.3 0.28574 0.30140 0.29999 0.30000 0.20696 0.86038 0.37256 0.960 0.4 0.38237 0.40332 0.39995 0.40000 0.22753 0.75995 0.50925 0.943 0.5 0.48022 0.50653 0.49985 0.50000 0.22499 0.64073 0.64073 0.934 N = 32 p QJ Can2011 Closer Exact |z[-1]| |z[0]| |z[+1]| cover --- ------- ------- ------- ------- ------- ------- ------- ----- 0.0 -0.00000-0.00000-0.00000-0.00000 0.00000 1.00000 0.00000 1.000 0.1 0.09968 0.10000 0.10000 0.10000 0.08960 0.98365 0.10943 0.994 0.2 0.19938 0.20003 0.20000 0.20000 0.15628 0.93555 0.23411 0.977 0.3 0.29912 0.30009 0.30000 0.30000 0.19863 0.85852 0.36817 0.955 0.4 0.39892 0.40021 0.40000 0.40000 0.21692 0.75702 0.50484 0.935 0.5 0.49879 0.50040 0.50000 0.50000 0.21298 0.63688 0.63688 0.926 N = 128 p QJ Can2011 Closer Exact |z[-1]| |z[0]| |z[+1]| cover --- ------- ------- ------- ------- ------- ------- ------- ----- 0.0 0.00000 0.00000 0.00000 0.00000 0.00000 1.00000 0.00000 1.000 0.1 0.09998 0.10000 0.10000 0.10000 0.08943 0.98363 0.10930 0.994 0.2 0.19996 0.20000 0.20000 0.20000 0.15594 0.93549 0.23389 0.977 0.3 0.29995 0.30001 0.30000 0.30000 0.19812 0.85840 0.36790 0.955 0.4 0.39993 0.40001 0.40000 0.40000 0.21628 0.75684 0.50457 0.935 0.5 0.49992 0.50003 0.50000 0.50000 0.21225 0.63664 0.63664 0.925
References
[1] Dawg, Cedron, DFT Bin Value Formulas for Pure Complex Tones, 2017
[2] Jacobsen, Eric and Kootsookos, Peter, "Fast, Accurate Frequency Estimators", IEEE SIGNAL PROCESSING MAGAZINE, pp. 123-125, MAY 2007
[3] Candan, Cagatay, "A Method For Fine Resolution Frequency Estimation From Three DFT Samples", IEEE SIGNAL PROCESSING LETTERS, VOL. 18, NO. 6, pp. 351-354, JUNE 2011
[4] C. CANDAN, "Analysis and further improvement of fine resolution frequency estimation method from three DFT samples", 2013
[5] Cedron Dawg, "Three Bin Exact Frequency Formulas for a Pure Complex Tone in a DFT", 2017
[6] Arzi, Julien, "Comparison of different frequency estimation algorithms", 2015
- Comments
- Write a Comment Select to add a comment
To post reply to a comment, click on the 'reply' button attached to each comment. To post a new comment (not a reply to a comment) check out the 'Write a Comment' tab at the top of the comments.
Please login (on the right) if you already have an account on this platform.
Otherwise, please use this form to register (free) an join one of the largest online community for Electrical/Embedded/DSP/FPGA/ML engineers:
