Hi Group
I have deduced/discovered some formulas
for interpolating between FFT bins when
calculating frequency and magnitude.
I have written a "formal" pdf document.
For posting I have made this text version.
I hope it still makes sense.
I would appreciate you comments and thoughts
very much. Do you agree with my reasoning?
TIA, Michael
Here is the document: (Sorry, it's a bit long)
---------------------
T h e S i n c E s t i m a t o r
___________________________________
by Michael Plet
1. FFT bins
____________
Consider a sine wave with a frequency of 860 Hz and an amplitude of 1.
Suppose we don't know the frequency and the amplitude.
One way of calculating these is using the Fast Fourier Transform (FFT).
Let's do that - We collect 256 samples and perform the transformation.
This leaves us with 256 pairs of complex numbers in two arrays:
Re(i) and Im(i), i=0,...,255.
For each of the array indices/bins we calculate the magnitude of the signal:
Mag(i) = Square Root(Re(i)^2 + Im(i)^2)
Now the bin with the largest magnitude represents the frequency of the
input signal like this:
Frequency = (SampleRate * BinNo) / FFTSize = (44100 * BinNo) / 256
Graphically it looks like this:
(Picture not available in text version)
The bin representing the input frequency has a calculated magnitude of 1.00.
All other bins have almost zero magnitude.
The above is a lucky case because the input frequency almost coincides with
the frequency represented by bin 5:
Bin 5 Freq. = (44100 * 5) / 256 = 861.3
Let's see what happens when our input frequency is 920 Hz and keep all
other conditions as before. Here is the graph:
(Picture not available in text version)
The peak is still at the bin corresponding to the frequency 861.3 Hz.
This time the calculated magnitude of the peak is 0.83.
In other words both the frequency and the magnitude are wrong.
The reason is that the input frequency lies between bins.
2. Quadratic Interpolation
___________________________
The neighboring bins seem to indicate that the true frequency is higher than
861.3 Hz because the next magnitude is larger than the previous.
The obvious (and standard) way around this is to interpolate the spectrum
by fitting a parabola through the peak and the neighboring bins.
Graphically it looks like this:
(Picture not available in text version)
If we call the three magnitudes m0, m1 and m2 (m1 being the peak),
the calculation gives
Estimated true Bin = PeakBinNo - (m2 - m0) / (2 * (m2 - 2 * m1 + m0))
In this case the values are m0 = 0.2227, m1 = 0.8317 and m2 = 0.4136.
As in the first example PeakBinNo is 5. Inserting the values gives
Estimted true Bin = 5.0929, corresponding to a frequency of 877.3 Hz.
The interpolated magnitude is
Estimated true Magnitude = m1 - (m2 - m0)^2 / (8 * (m2 - 2 * m1 + m0))
= 0.8361
These results could be improved by increasing the number of samples used
for the FFT.
Another way would be to improve the interpolation itself.
3. Sinc Interpolation
______________________
The FFT transforms a rectangular window to a Sinc function.
This function has the definition
Sinc(x) = Sin(PI * x) / (PI * x), for x =/= 0
Sinc(x) = 1 , for x = 0
Drawing this function (actually the absolute value) along with our first
example looks like this:
(Picture not available in text version)
Notice that the value is zero for all bins except the peak.
Now, if we draw the Sinc function with our second example we get this
(Picture not available in text version)
It is very clear that it fits nicely to the values of the bins.
The question is how to calculate the frequency and magnitude of the true
peak.
From the graph it is intuitively clear (and it can be proven), that only the
peak
and it's largest neighbor, will be inside the area of the main peak of the
Sinc function.
Let's call these m0 and m1 with the corresponding bins b0 and b1.
Notice that the peak can be either m0 or m1. They are just the peak and it's
largest neighbor. Also note that b1 is the bin following b0.
So the points on the sinc curve are (b0, m0) and (b1, m1).
The true peak is the peak of the Sinc curve. Let's call this point (b', m').
Because the Sinc function has it's maximum at x = 0 , we must do a parallel
shift (in our minds) of the graph so that the maximum is at (0, m').
Subtracting b' from the other points makes these (b0 - b', m0) and (b1 - b',
m1).
As they are both points on the graph of the Sinc function they fulfill
m0 = Sin(PI * (b0 - b')) / (PI * (b0 - b')) I)
and
m1 = Sin(PI * (b1 - b')) / (PI * (b1 - b'))
as b1 = b0 + 1 and Sin(x + PI) = -Sin(x) we get
m1 = Sin(PI * (b0 - b') + PI) / (PI * (b0 + 1 - b'))
= Sin(PI * (b0 - b')) / (PI * (b' - b0 - 1)) II)
Noticing that the numerator of both I) and II) are the same, we combine the
equations
m0 * PI * (b0 - b') = m1 * PI * (b' - b0 - 1)
<=>
m0 * b0 - m0 * b' = m1 * b' - m1 * b0 - m1
<=>
b' = (m0 * b0 + m1 * b0 + m1) / (m0 + m1)
<=>
b' = b0 + (m1 / (m0 + m1)) "The Sinc Frequency Estimator"
Fortunately this has reduced to a very easy-to-calculate estimator.
Inserting our values from the second example (920 Hz) gives
b' = 5 + (0.4136 / (0.8317 + 0.4136)) = 5.3321
This corresponds to a frequency of (44100 * 5.3321) / 256 = 918.54 Hz.
A much better estimate than the Quadratic Estimator's 877.3 Hz.
Now for the magnitude - this is easy. Since the interpolating curve is just
a
scaled Sinc function, all points must fulfill
m(x) = ScaleFactor * Sinc(x - b') III)
The maximum peak is at x = b', so the ScaleFactor is m'
m(b') = m' = ScaleFactor
Inserting the point (b0, m0) into III) gives
m(b0) = m0 = m' * Sinc(b0 - b')
Now the estimate of the magnitude can be expressed as follows
m' = m0 / Sinc(b0 - b')
This can be "reduced" to
"The Sinc Magnitude Estimator"
m' = (PI * m0 * m1) / ((m0 + m1) * Sin((PI * m0) / (m0 + m1)))
Since Sin((PI * m0) / (m0 + m1)) = Sin((PI * m1) / (m0 + m1)) this can also
be written
m' = (PI * m0 * m1) / ((m0 + m1) * Sin((PI * m1) / (m0 + m1)))
Inserting the values gives us the estimate
m' = (3.1416 * 0.8317 * 0.4136)
/ ((0.8317 + 0.4136) * Sin((3.1416 * 0.8317) / (0.8317 + 0.4136)))
= 1.0807 / (1.2453 * Sin(2.0982))
= 1.0042
which is much closer to 1 (the true magnitude) than the quadratic estimate
of 0.8361.
4. Results
___________
A wide range of test has been run to test the Sinc Estimator varying both
the
frequency, amplitude and the FFT size.
All tests have shown that the Sinc Estimator is extremely stable over the
full
spectrum of frequencies from 0 Hz to 22050 Hz.
Increasing the FFT size will make the Quadratic Estimator more precise.
The Sinc Estimator is still better, but to a lesser degree.
Windows.
--------
When "windowing" the samples prior to doing the FFT, we are no longer
transforming a "rectangular window" (which is no window).
This affects the Sinc Estimator and it is not better than the Quadratic
Estimator
(In some cases it is worse).
Zero-padding.
-------------
The effect of zero-padding on the FFT result is very interesting.
The diagram shows the previous example but zero-padded:
(Picture not available in text version)
What we see here is "The Sinc Estimator" failing completely.
The reason is just as obvious: The sinc curve doesn't fit the bins.
The frequency estimate will be ok, but the magnitude estimate is far to big.
Even more interesting is the pattern of the bin magnitudes. They still
sketch
a sinc curve. This means that we could modify the sinc curve and thereby
the estimator. We could make it wider and shorter taking into account the
degree of zero-padding.
Not much would be gained by this though, because with this much
zero-padding,
the Quadratic Estimator would be just as precise.
5. Conclusion
______________
We have deduced "The Sinc Estimator" and shown it to be more precise than
the estimator normally used.
It has also been established that "The Sinc Estimator" should be used only
on
non-windowed signals with no zero-padding.
All calculations and graphs have been made using a computer program.
My Sinc Estmator
Started by ●May 17, 2011
Reply by ●May 17, 20112011-05-17
On 5/17/2011 2:43 PM, Michael P. wrote:> Hi Group > > I have deduced/discovered some formulas > for interpolating between FFT bins when > calculating frequency and magnitude. > > I have written a "formal" pdf document. > For posting I have made this text version. > I hope it still makes sense. > > I would appreciate you comments and thoughts > very much. Do you agree with my reasoning? > > TIA, Michael > > Here is the document: (Sorry, it's a bit long) > > --------------------- > > > > T h e S i n c E s t i m a t o r > ___________________________________ > > by Michael Plet > > > 1. FFT bins > ____________ > > Consider a sine wave with a frequency of 860 Hz and an amplitude of 1. > Suppose we don't know the frequency and the amplitude. > One way of calculating these is using the Fast Fourier Transform (FFT). > Let's do that - We collect 256 samples and perform the transformation. > This leaves us with 256 pairs of complex numbers in two arrays: > > Re(i) and Im(i), i=0,...,255. > > For each of the array indices/bins we calculate the magnitude of the > signal: > > Mag(i) = Square Root(Re(i)^2 + Im(i)^2) > > Now the bin with the largest magnitude represents the frequency of the > input signal like this: > > Frequency = (SampleRate * BinNo) / FFTSize = (44100 * BinNo) / 256 > > Graphically it looks like this: > (Picture not available in text version) > > The bin representing the input frequency has a calculated magnitude of > 1.00. > All other bins have almost zero magnitude. > > The above is a lucky case because the input frequency almost coincides with > the frequency represented by bin 5: > > Bin 5 Freq. = (44100 * 5) / 256 = 861.3 > > Let's see what happens when our input frequency is 920 Hz and keep all > other conditions as before. Here is the graph: > (Picture not available in text version) > > The peak is still at the bin corresponding to the frequency 861.3 Hz. > This time the calculated magnitude of the peak is 0.83. > In other words both the frequency and the magnitude are wrong. > The reason is that the input frequency lies between bins. > > > 2. Quadratic Interpolation > ___________________________ > > The neighboring bins seem to indicate that the true frequency is higher > than > 861.3 Hz because the next magnitude is larger than the previous. > The obvious (and standard) way around this is to interpolate the spectrum > by fitting a parabola through the peak and the neighboring bins. > Graphically it looks like this: > (Picture not available in text version) > > If we call the three magnitudes m0, m1 and m2 (m1 being the peak), > the calculation gives > > Estimated true Bin = PeakBinNo - (m2 - m0) / (2 * (m2 - 2 * m1 + m0)) > > In this case the values are m0 = 0.2227, m1 = 0.8317 and m2 = 0.4136. > As in the first example PeakBinNo is 5. Inserting the values gives > > Estimted true Bin = 5.0929, corresponding to a frequency of 877.3 Hz. > > The interpolated magnitude is > > Estimated true Magnitude = m1 - (m2 - m0)^2 / (8 * (m2 - 2 * m1 + m0)) > = 0.8361 > > These results could be improved by increasing the number of samples used > for the FFT. > Another way would be to improve the interpolation itself. > > > 3. Sinc Interpolation > ______________________ > > The FFT transforms a rectangular window to a Sinc function. > This function has the definition > > Sinc(x) = Sin(PI * x) / (PI * x), for x =/= 0 > > Sinc(x) = 1 , for x = 0 > > Drawing this function (actually the absolute value) along with our first > example looks like this: > (Picture not available in text version) > > Notice that the value is zero for all bins except the peak. > > Now, if we draw the Sinc function with our second example we get this > (Picture not available in text version) > > It is very clear that it fits nicely to the values of the bins. > The question is how to calculate the frequency and magnitude of the true > peak. > > From the graph it is intuitively clear (and it can be proven), that > only the peak > and it's largest neighbor, will be inside the area of the main peak of > the Sinc function. > Let's call these m0 and m1 with the corresponding bins b0 and b1. > Notice that the peak can be either m0 or m1. They are just the peak and > it's > largest neighbor. Also note that b1 is the bin following b0. > So the points on the sinc curve are (b0, m0) and (b1, m1). > > The true peak is the peak of the Sinc curve. Let's call this point (b', > m'). > Because the Sinc function has it's maximum at x = 0 , we must do a parallel > shift (in our minds) of the graph so that the maximum is at (0, m'). > Subtracting b' from the other points makes these (b0 - b', m0) and (b1 - > b', m1). > > As they are both points on the graph of the Sinc function they fulfill > > m0 = Sin(PI * (b0 - b')) / (PI * (b0 - b')) I) > > and > > m1 = Sin(PI * (b1 - b')) / (PI * (b1 - b')) > > as b1 = b0 + 1 and Sin(x + PI) = -Sin(x) we get > > m1 = Sin(PI * (b0 - b') + PI) / (PI * (b0 + 1 - b')) > > = Sin(PI * (b0 - b')) / (PI * (b' - b0 - 1)) II) > > Noticing that the numerator of both I) and II) are the same, we combine > the equations > > m0 * PI * (b0 - b') = m1 * PI * (b' - b0 - 1) > > <=> > > m0 * b0 - m0 * b' = m1 * b' - m1 * b0 - m1 > > <=> > > b' = (m0 * b0 + m1 * b0 + m1) / (m0 + m1) > > <=> > > > b' = b0 + (m1 / (m0 + m1)) "The Sinc Frequency Estimator" > > > Fortunately this has reduced to a very easy-to-calculate estimator. > Inserting our values from the second example (920 Hz) gives > > b' = 5 + (0.4136 / (0.8317 + 0.4136)) = 5.3321 > > This corresponds to a frequency of (44100 * 5.3321) / 256 = 918.54 Hz. > A much better estimate than the Quadratic Estimator's 877.3 Hz. > > Now for the magnitude - this is easy. Since the interpolating curve is > just a > scaled Sinc function, all points must fulfill > > m(x) = ScaleFactor * Sinc(x - b') III) > > The maximum peak is at x = b', so the ScaleFactor is m' > > m(b') = m' = ScaleFactor > > Inserting the point (b0, m0) into III) gives > > m(b0) = m0 = m' * Sinc(b0 - b') > > Now the estimate of the magnitude can be expressed as follows > > m' = m0 / Sinc(b0 - b') > > This can be "reduced" to > > > "The Sinc Magnitude Estimator" > > m' = (PI * m0 * m1) / ((m0 + m1) * Sin((PI * m0) / (m0 + m1))) > > > Since Sin((PI * m0) / (m0 + m1)) = Sin((PI * m1) / (m0 + m1)) this can > also be written > > m' = (PI * m0 * m1) / ((m0 + m1) * Sin((PI * m1) / (m0 + m1))) > > Inserting the values gives us the estimate > > m' = (3.1416 * 0.8317 * 0.4136) > / ((0.8317 + 0.4136) * Sin((3.1416 * 0.8317) / (0.8317 + 0.4136))) > > = 1.0807 / (1.2453 * Sin(2.0982)) > > = 1.0042 > > which is much closer to 1 (the true magnitude) than the quadratic > estimate of 0.8361. > > > 4. Results > ___________ > > A wide range of test has been run to test the Sinc Estimator varying > both the > frequency, amplitude and the FFT size. > All tests have shown that the Sinc Estimator is extremely stable over > the full > spectrum of frequencies from 0 Hz to 22050 Hz. > Increasing the FFT size will make the Quadratic Estimator more precise. > The Sinc Estimator is still better, but to a lesser degree. > > Windows. > -------- > When "windowing" the samples prior to doing the FFT, we are no longer > transforming a "rectangular window" (which is no window). > This affects the Sinc Estimator and it is not better than the Quadratic > Estimator > (In some cases it is worse). > > Zero-padding. > ------------- > The effect of zero-padding on the FFT result is very interesting. > The diagram shows the previous example but zero-padded: > (Picture not available in text version) > > What we see here is "The Sinc Estimator" failing completely. > The reason is just as obvious: The sinc curve doesn't fit the bins. > > The frequency estimate will be ok, but the magnitude estimate is far to > big. > > Even more interesting is the pattern of the bin magnitudes. They still > sketch > a sinc curve. This means that we could modify the sinc curve and thereby > the estimator. We could make it wider and shorter taking into account the > degree of zero-padding. > Not much would be gained by this though, because with this much > zero-padding, > the Quadratic Estimator would be just as precise. > > > 5. Conclusion > ______________ > > We have deduced "The Sinc Estimator" and shown it to be more precise than > the estimator normally used. > It has also been established that "The Sinc Estimator" should be used > only on > non-windowed signals with no zero-padding. > All calculations and graphs have been made using a computer program. > >Kinda tedious to read without the figures / .pdf What I got from a "quick read" is that you're matching an FFT of a pure sine wave to suitably shifted sinc values. Is that right? If so, the it's only good for pure sinusoids without noise? Is that right? Fred
Reply by ●May 17, 20112011-05-17
Michael P. wrote:> Hi Group > > I have deduced/discovered some formulas > for interpolating between FFT bins when > calculating frequency and magnitude. > > I have written a "formal" pdf document. > For posting I have made this text version. > I hope it still makes sense. > > I would appreciate you comments and thoughts > very much. Do you agree with my reasoning? > > TIA, Michael//------ "Fast, Accurate frequency estimators" E. Jacobsen, P. Kootsookos. http://www.ingelec.uns.edu.ar/pds2803/Materiales/Articulos/AnalisisFrecuencial/04205098.pdf //-------- Vladimir Vassilevsky DSP and Mixed Signal Design Consultant http://www.abvolt.com
Reply by ●May 17, 20112011-05-17
Michael P. <me@home.com> wrote: (snip)> T h e S i n c E s t i m a t o r > ___________________________________ > > by Michael Plet> 1. FFT bins > ____________> Consider a sine wave with a frequency of 860 Hz and an amplitude of 1. > Suppose we don't know the frequency and the amplitude. > One way of calculating these is using the Fast Fourier Transform (FFT). > Let's do that - We collect 256 samples and perform the transformation. > This leaves us with 256 pairs of complex numbers in two arrays:If you have one frequency component, then the FFT isn't usually the best way to find that frequency. If you need a more accurate determination of the frequency components of a complicated signal, then you need more samples. -- glen
Reply by ●May 18, 20112011-05-18
On May 17, 6:29�pm, Fred Marshall <fmarshallxremove_th...@acm.org> wrote:> On 5/17/2011 2:43 PM, Michael P. wrote: > > > > > Hi Group > > > I have deduced/discovered some formulas > > for interpolating between FFT bins when > > calculating frequency and magnitude. > > > I have written a "formal" pdf document. > > For posting I have made this text version. > > I hope it still makes sense. > > > I would appreciate you comments and thoughts > > very much. Do you agree with my reasoning? > > > TIA, Michael > > > Here is the document: (Sorry, it's a bit long) > > > --------------------- > > > T h e S i n c E s t i m a t o r > > ___________________________________ > > > by Michael Plet > > > 1. FFT bins > > ____________ > > > Consider a sine wave with a frequency of 860 Hz and an amplitude of 1. > > Suppose we don't know the frequency and the amplitude. > > One way of calculating these is using the Fast Fourier Transform (FFT). > > Let's do that - We collect 256 samples and perform the transformation. > > This leaves us with 256 pairs of complex numbers in two arrays: > > > Re(i) and Im(i), i=0,...,255. > > > For each of the array indices/bins we calculate the magnitude of the > > signal: > > > Mag(i) = Square Root(Re(i)^2 + Im(i)^2) > > > Now the bin with the largest magnitude represents the frequency of the > > input signal like this: > > > Frequency = (SampleRate * BinNo) / FFTSize = (44100 * BinNo) / 256 > > > Graphically it looks like this: > > (Picture not available in text version) > > > The bin representing the input frequency has a calculated magnitude of > > 1.00. > > All other bins have almost zero magnitude. > > > The above is a lucky case because the input frequency almost coincides with > > the frequency represented by bin 5: > > > Bin 5 Freq. = (44100 * 5) / 256 = 861.3 > > > Let's see what happens when our input frequency is 920 Hz and keep all > > other conditions as before. Here is the graph: > > (Picture not available in text version) > > > The peak is still at the bin corresponding to the frequency 861.3 Hz. > > This time the calculated magnitude of the peak is 0.83. > > In other words both the frequency and the magnitude are wrong. > > The reason is that the input frequency lies between bins. > > > 2. Quadratic Interpolation > > ___________________________ > > > The neighboring bins seem to indicate that the true frequency is higher > > than > > 861.3 Hz because the next magnitude is larger than the previous. > > The obvious (and standard) way around this is to interpolate the spectrum > > by fitting a parabola through the peak and the neighboring bins. > > Graphically it looks like this: > > (Picture not available in text version) > > > If we call the three magnitudes m0, m1 and m2 (m1 being the peak), > > the calculation gives > > > Estimated true Bin = PeakBinNo - (m2 - m0) / (2 * (m2 - 2 * m1 + m0)) > > > In this case the values are m0 = 0.2227, m1 = 0.8317 and m2 = 0.4136. > > As in the first example PeakBinNo is 5. Inserting the values gives > > > Estimted true Bin = 5.0929, corresponding to a frequency of 877.3 Hz. > > > The interpolated magnitude is > > > Estimated true Magnitude = m1 - (m2 - m0)^2 / (8 * (m2 - 2 * m1 + m0)) > > = 0.8361 > > > These results could be improved by increasing the number of samples used > > for the FFT. > > Another way would be to improve the interpolation itself. > > > 3. Sinc Interpolation > > ______________________ > > > The FFT transforms a rectangular window to a Sinc function. > > This function has the definition > > > Sinc(x) = Sin(PI * x) / (PI * x), for x =/= 0 > > > Sinc(x) = 1 , for x = 0 > > > Drawing this function (actually the absolute value) along with our first > > example looks like this: > > (Picture not available in text version) > > > Notice that the value is zero for all bins except the peak. > > > Now, if we draw the Sinc function with our second example we get this > > (Picture not available in text version) > > > It is very clear that it fits nicely to the values of the bins. > > The question is how to calculate the frequency and magnitude of the true > > peak. > > > �From the graph it is intuitively clear (and it can be proven), that > > only the peak > > and it's largest neighbor, will be inside the area of the main peak of > > the Sinc function. > > Let's call these m0 and m1 with the corresponding bins b0 and b1. > > Notice that the peak can be either m0 or m1. They are just the peak and > > it's > > largest neighbor. Also note that b1 is the bin following b0. > > So the points on the sinc curve are (b0, m0) and (b1, m1). > > > The true peak is the peak of the Sinc curve. Let's call this point (b', > > m'). > > Because the Sinc function has it's maximum at x = 0 , we must do a parallel > > shift (in our minds) of the graph so that the maximum is at (0, m'). > > Subtracting b' from the other points makes these (b0 - b', m0) and (b1 - > > b', m1). > > > As they are both points on the graph of the Sinc function they fulfill > > > m0 = Sin(PI * (b0 - b')) / (PI * (b0 - b')) I) > > > and > > > m1 = Sin(PI * (b1 - b')) / (PI * (b1 - b')) > > > as b1 = b0 + 1 and Sin(x + PI) = -Sin(x) we get > > > m1 = Sin(PI * (b0 - b') + PI) / (PI * (b0 + 1 - b')) > > > = Sin(PI * (b0 - b')) / (PI * (b' - b0 - 1)) II) > > > Noticing that the numerator of both I) and II) are the same, we combine > > the equations > > > m0 * PI * (b0 - b') = m1 * PI * (b' - b0 - 1) > > > <=> > > > m0 * b0 - m0 * b' = m1 * b' - m1 * b0 - m1 > > > <=> > > > b' = (m0 * b0 + m1 * b0 + m1) / (m0 + m1) > > > <=> > > > b' = b0 + (m1 / (m0 + m1)) "The Sinc Frequency Estimator" > > > Fortunately this has reduced to a very easy-to-calculate estimator. > > Inserting our values from the second example (920 Hz) gives > > > b' = 5 + (0.4136 / (0.8317 + 0.4136)) = 5.3321 > > > This corresponds to a frequency of (44100 * 5.3321) / 256 = 918.54 Hz. > > A much better estimate than the Quadratic Estimator's 877.3 Hz. > > > Now for the magnitude - this is easy. Since the interpolating curve is > > just a > > scaled Sinc function, all points must fulfill > > > m(x) = ScaleFactor * Sinc(x - b') III) > > > The maximum peak is at x = b', so the ScaleFactor is m' > > > m(b') = m' = ScaleFactor > > > Inserting the point (b0, m0) into III) gives > > > m(b0) = m0 = m' * Sinc(b0 - b') > > > Now the estimate of the magnitude can be expressed as follows > > > m' = m0 / Sinc(b0 - b') > > > This can be "reduced" to > > > "The Sinc Magnitude Estimator" > > > m' = (PI * m0 * m1) / ((m0 + m1) * Sin((PI * m0) / (m0 + m1))) > > > Since Sin((PI * m0) / (m0 + m1)) = Sin((PI * m1) / (m0 + m1)) this can > > also be written > > > m' = (PI * m0 * m1) / ((m0 + m1) * Sin((PI * m1) / (m0 + m1))) > > > Inserting the values gives us the estimate > > > m' = (3.1416 * 0.8317 * 0.4136) > > / ((0.8317 + 0.4136) * Sin((3.1416 * 0.8317) / (0.8317 + 0.4136))) > > > = 1.0807 / (1.2453 * Sin(2.0982)) > > > = 1.0042 > > > which is much closer to 1 (the true magnitude) than the quadratic > > estimate of 0.8361. > > > 4. Results > > ___________ > > > A wide range of test has been run to test the Sinc Estimator varying > > both the > > frequency, amplitude and the FFT size. > > All tests have shown that the Sinc Estimator is extremely stable over > > the full > > spectrum of frequencies from 0 Hz to 22050 Hz. > > Increasing the FFT size will make the Quadratic Estimator more precise. > > The Sinc Estimator is still better, but to a lesser degree. > > > Windows. > > -------- > > When "windowing" the samples prior to doing the FFT, we are no longer > > transforming a "rectangular window" (which is no window). > > This affects the Sinc Estimator and it is not better than the Quadratic > > Estimator > > (In some cases it is worse). > > > Zero-padding. > > ------------- > > The effect of zero-padding on the FFT result is very interesting. > > The diagram shows the previous example but zero-padded: > > (Picture not available in text version) > > > What we see here is "The Sinc Estimator" failing completely. > > The reason is just as obvious: The sinc curve doesn't fit the bins. > > > The frequency estimate will be ok, but the magnitude estimate is far to > > big. > > > Even more interesting is the pattern of the bin magnitudes. They still > > sketch > > a sinc curve. This means that we could modify the sinc curve and thereby > > the estimator. We could make it wider and shorter taking into account the > > degree of zero-padding. > > Not much would be gained by this though, because with this much > > zero-padding, > > the Quadratic Estimator would be just as precise. > > > 5. Conclusion > > ______________ > > > We have deduced "The Sinc Estimator" and shown it to be more precise than > > the estimator normally used. > > It has also been established that "The Sinc Estimator" should be used > > only on > > non-windowed signals with no zero-padding. > > All calculations and graphs have been made using a computer program. > > Kinda tedious to read without the figures / .pdf > > What I got from a "quick read" is that you're matching an FFT of a pure > sine wave to suitably shifted sinc values. �Is that right? > If so, the it's only good for pure sinusoids without noise? �Is that right? > > FredFred, given a pure sinusoid that is between 2 FFT bins the structure of the magnitude of the near by FFT bins is deterministic. He's trying to use that structure and interpolate to find where the peak would be and also it's associated frequency. The method is very sensitive to the window - since in the Frequency domain the window is convolved with the data. So I guess you could see the filter as shift The method works well for pure sinusoids where they are separated by enough FFT bins. How many bins depends on the length of the Filter and the height of the sidelobes of the window being used - to prevent sprectal leakage from nearby sinusoids. The method works well with wideband noise, since in the frequency domain the energy is spread out across frequency bins - compared to the sinusoid. Cheers, Dave
Reply by ●May 18, 20112011-05-18
Have you looked at http://www.dsprelated.com/showcode/164.php and http://www.dsprelated.com/showcode/159.php ? Jesus.
Reply by ●May 18, 20112011-05-18
Dave <dspguy2@netscape.net> wrote: (big snip on FFT interpolation)> Fred, given a pure sinusoid that is between 2 FFT bins the structure > of the magnitude of the near by FFT bins is deterministic. He's trying > to use that structure and interpolate to find where the peak would be > and also it's associated frequency. The method is very sensitive to > the window - since in the Frequency domain the window is convolved > with the data. So I guess you could see the filter as shiftBut if it is just one sinusoid, I would expect something like least squares to be faster, and get a more accurate result. Even for a small number of sinusoids, it should work.> The method works well for pure sinusoids where they are separated by > enough FFT bins. How many bins depends on the length of the Filter and > the height of the sidelobes of the window being used - to prevent > sprectal leakage from nearby sinusoids. The method works well with > wideband noise, since in the frequency domain the energy is spread out > across frequency bins - compared to the sinusoid.Least squares should work even if they are close, though it is harder to get it right. Some years ago I was doing least squares fits on sums of Gaussians. Sometimes it was sensitive to the initial values (guess), other times it would converge even with bad guesses. More sample points will make both least squares and FFT give better results for a sum of sinusoids. -- glen
Reply by ●May 18, 20112011-05-18
On Tue, 17 May 2011 23:43:46 +0200, "Michael P." <me@home.com> wrote:>Hi Group > >I have deduced/discovered some formulas >for interpolating between FFT bins when >calculating frequency and magnitude. > >I have written a "formal" pdf document. >For posting I have made this text version. >I hope it still makes sense. > >I would appreciate you comments and thoughts >very much. Do you agree with my reasoning? > >TIA, Michael > >Here is the document: (Sorry, it's a bit long) > >--------------------- > > > > T h e S i n c E s t i m a t o r > ___________________________________ > > by Michael Plet > > >1. FFT bins >____________ > >Consider a sine wave with a frequency of 860 Hz and an amplitude of 1. >Suppose we don't know the frequency and the amplitude. >One way of calculating these is using the Fast Fourier Transform (FFT). >Let's do that - We collect 256 samples and perform the transformation. >This leaves us with 256 pairs of complex numbers in two arrays: > > Re(i) and Im(i), i=0,...,255. > >For each of the array indices/bins we calculate the magnitude of the signal: > > Mag(i) = Square Root(Re(i)^2 + Im(i)^2) > >Now the bin with the largest magnitude represents the frequency of the >input signal like this: > >Frequency = (SampleRate * BinNo) / FFTSize = (44100 * BinNo) / 256 > >Graphically it looks like this: >(Picture not available in text version) > >The bin representing the input frequency has a calculated magnitude of 1.00. >All other bins have almost zero magnitude. > >The above is a lucky case because the input frequency almost coincides with >the frequency represented by bin 5: > > Bin 5 Freq. = (44100 * 5) / 256 = 861.3 > >Let's see what happens when our input frequency is 920 Hz and keep all >other conditions as before. Here is the graph: >(Picture not available in text version) > >The peak is still at the bin corresponding to the frequency 861.3 Hz. >This time the calculated magnitude of the peak is 0.83. >In other words both the frequency and the magnitude are wrong. >The reason is that the input frequency lies between bins. > > >2. Quadratic Interpolation >___________________________ > >The neighboring bins seem to indicate that the true frequency is higher than >861.3 Hz because the next magnitude is larger than the previous. >The obvious (and standard) way around this is to interpolate the spectrum >by fitting a parabola through the peak and the neighboring bins. >Graphically it looks like this: >(Picture not available in text version) > >If we call the three magnitudes m0, m1 and m2 (m1 being the peak), >the calculation gives > > Estimated true Bin = PeakBinNo - (m2 - m0) / (2 * (m2 - 2 * m1 + m0)) > >In this case the values are m0 = 0.2227, m1 = 0.8317 and m2 = 0.4136. >As in the first example PeakBinNo is 5. Inserting the values gives > > Estimted true Bin = 5.0929, corresponding to a frequency of 877.3 Hz. > >The interpolated magnitude is > > Estimated true Magnitude = m1 - (m2 - m0)^2 / (8 * (m2 - 2 * m1 + m0)) > = 0.8361 > >These results could be improved by increasing the number of samples used >for the FFT. >Another way would be to improve the interpolation itself. > > >3. Sinc Interpolation >______________________ > >The FFT transforms a rectangular window to a Sinc function. >This function has the definition > > Sinc(x) = Sin(PI * x) / (PI * x), for x =/= 0 > > Sinc(x) = 1 , for x = 0 > >Drawing this function (actually the absolute value) along with our first >example looks like this: >(Picture not available in text version) > >Notice that the value is zero for all bins except the peak. > >Now, if we draw the Sinc function with our second example we get this >(Picture not available in text version) > >It is very clear that it fits nicely to the values of the bins. >The question is how to calculate the frequency and magnitude of the true >peak. > >From the graph it is intuitively clear (and it can be proven), that only the >peak >and it's largest neighbor, will be inside the area of the main peak of the >Sinc function. >Let's call these m0 and m1 with the corresponding bins b0 and b1. >Notice that the peak can be either m0 or m1. They are just the peak and it's >largest neighbor. Also note that b1 is the bin following b0. >So the points on the sinc curve are (b0, m0) and (b1, m1). > >The true peak is the peak of the Sinc curve. Let's call this point (b', m'). >Because the Sinc function has it's maximum at x = 0 , we must do a parallel >shift (in our minds) of the graph so that the maximum is at (0, m'). >Subtracting b' from the other points makes these (b0 - b', m0) and (b1 - b', >m1). > >As they are both points on the graph of the Sinc function they fulfill > > m0 = Sin(PI * (b0 - b')) / (PI * (b0 - b')) I) > >and > > m1 = Sin(PI * (b1 - b')) / (PI * (b1 - b')) > >as b1 = b0 + 1 and Sin(x + PI) = -Sin(x) we get > > m1 = Sin(PI * (b0 - b') + PI) / (PI * (b0 + 1 - b')) > > = Sin(PI * (b0 - b')) / (PI * (b' - b0 - 1)) II) > >Noticing that the numerator of both I) and II) are the same, we combine the >equations > > m0 * PI * (b0 - b') = m1 * PI * (b' - b0 - 1) > > <=> > > m0 * b0 - m0 * b' = m1 * b' - m1 * b0 - m1 > > <=> > > b' = (m0 * b0 + m1 * b0 + m1) / (m0 + m1) > > <=> > > > b' = b0 + (m1 / (m0 + m1)) "The Sinc Frequency Estimator" > > >Fortunately this has reduced to a very easy-to-calculate estimator. >Inserting our values from the second example (920 Hz) gives > > b' = 5 + (0.4136 / (0.8317 + 0.4136)) = 5.3321 > >This corresponds to a frequency of (44100 * 5.3321) / 256 = 918.54 Hz. >A much better estimate than the Quadratic Estimator's 877.3 Hz. > >Now for the magnitude - this is easy. Since the interpolating curve is just >a >scaled Sinc function, all points must fulfill > > m(x) = ScaleFactor * Sinc(x - b') III) > >The maximum peak is at x = b', so the ScaleFactor is m' > > m(b') = m' = ScaleFactor > >Inserting the point (b0, m0) into III) gives > > m(b0) = m0 = m' * Sinc(b0 - b') > >Now the estimate of the magnitude can be expressed as follows > > m' = m0 / Sinc(b0 - b') > >This can be "reduced" to > > > "The Sinc Magnitude Estimator" > > m' = (PI * m0 * m1) / ((m0 + m1) * Sin((PI * m0) / (m0 + m1))) > > >Since Sin((PI * m0) / (m0 + m1)) = Sin((PI * m1) / (m0 + m1)) this can also >be written > > m' = (PI * m0 * m1) / ((m0 + m1) * Sin((PI * m1) / (m0 + m1))) > >Inserting the values gives us the estimate > > m' = (3.1416 * 0.8317 * 0.4136) > / ((0.8317 + 0.4136) * Sin((3.1416 * 0.8317) / (0.8317 + 0.4136))) > > = 1.0807 / (1.2453 * Sin(2.0982)) > > = 1.0042 > >which is much closer to 1 (the true magnitude) than the quadratic estimate >of 0.8361. > > >4. Results >___________ > >A wide range of test has been run to test the Sinc Estimator varying both >the >frequency, amplitude and the FFT size. >All tests have shown that the Sinc Estimator is extremely stable over the >full >spectrum of frequencies from 0 Hz to 22050 Hz. >Increasing the FFT size will make the Quadratic Estimator more precise. >The Sinc Estimator is still better, but to a lesser degree. > >Windows. >-------- >When "windowing" the samples prior to doing the FFT, we are no longer >transforming a "rectangular window" (which is no window). >This affects the Sinc Estimator and it is not better than the Quadratic >Estimator >(In some cases it is worse). > >Zero-padding. >------------- >The effect of zero-padding on the FFT result is very interesting. >The diagram shows the previous example but zero-padded: >(Picture not available in text version) > >What we see here is "The Sinc Estimator" failing completely. >The reason is just as obvious: The sinc curve doesn't fit the bins. > >The frequency estimate will be ok, but the magnitude estimate is far to big. > >Even more interesting is the pattern of the bin magnitudes. They still >sketch >a sinc curve. This means that we could modify the sinc curve and thereby >the estimator. We could make it wider and shorter taking into account the >degree of zero-padding. >Not much would be gained by this though, because with this much >zero-padding, >the Quadratic Estimator would be just as precise. > > >5. Conclusion >______________ > >We have deduced "The Sinc Estimator" and shown it to be more precise than >the estimator normally used. >It has also been established that "The Sinc Estimator" should be used only >on >non-windowed signals with no zero-padding. >All calculations and graphs have been made using a computer program.Without any comparative performance plots, etc., it's kind of tough to comment on how this compares to similar techniques. There's a quick survey of some related techniques here, and matlab code is included. http://www.ericjacobsen.org/fe2/fe2.htm These sorts of things are fairly well covered in the literature, but there is still a tradeoff between some of the fine details in performance and things like estimator complexity, etc. Eric Jacobsen http://www.ericjacobsen.org http://www.dsprelated.com/blogs-1//Eric_Jacobsen.php
Reply by ●May 18, 20112011-05-18
"glen herrmannsfeldt" <gah@ugcs.caltech.edu> wrote in message news:iqv42u$e4b$1@dont-email.me...> Michael P. <me@home.com> wrote: > > (snip) > >> T h e S i n c E s t i m a t o r >> ___________________________________ >> >> by Michael Plet > >> 1. FFT bins >> ____________ > >> Consider a sine wave with a frequency of 860 Hz and an amplitude of 1. >> Suppose we don't know the frequency and the amplitude. >> One way of calculating these is using the Fast Fourier Transform (FFT). >> Let's do that - We collect 256 samples and perform the transformation. >> This leaves us with 256 pairs of complex numbers in two arrays: > > If you have one frequency component, then the FFT isn't usually > the best way to find that frequency. If you need a more accurate > determination of the frequency components of a complicated signal, > then you need more samples. > > -- glenYou are probably right. But I think of this as a theoretical exercise. So I am just interested in the problem of interpolating between bins. Speed doesn't matter either. Michael
Reply by ●May 18, 20112011-05-18
"Michael P." <me@home.com> wrote in message news:TcednVNkjp6OcU_QnZ2dnUVZ8iadnZ2d@giganews.com...> Hi Group > > I have deduced/discovered some formulas > for interpolating between FFT bins when > calculating frequency and magnitude. > > I have written a "formal" pdf document. > For posting I have made this text version. > I hope it still makes sense. > > I would appreciate you comments and thoughts > very much. Do you agree with my reasoning? > > TIA, Michael > > Here is the document: (Sorry, it's a bit long) >Hi Group. Thank you for all your replies so far! I know that there are other methods and formulas. My formulas were deduced using pure mathematics from the fact that FFT transforms a rectangular window to a Sinc curve. I just wanted your opinion about my formulas. I understand that it is difficult to read without the graphs. I will find a way of uploading the pdf file somewhere. Michael
