Empirically Optimized FFT bin interpolator

Here is an FFT bin interpolator I have been using with great success.  It is
based on the standard quadratic interpolation formula:

float CalcFFTInterpolation(float yPrecedingPower, float yMaxPower, float
yFollowingPower)
   // yMaxPower is the power amplitude in a local peak
   // yPrecedingPower is the power amplitude of the preceding bin
   // yFollowingPower is the power amplitude of the following bin
{
  float c = 2 * yMaxPower - yFollowingPower - yPrecedingPower;
  if(c)
  {
    int index = (yMaxPower - yPrecedingPower) / c * 64;
    if(index > 63) index = 63;
      c = InterFFTLookup[index];  //  0...63 covers -.50 to +.49 bins
  }
  return c;
}

The returned value of the function is the bin index offset that you should add
to the index where yMaxPower was found to get the interpolated frequency of this
peak.  Everything up to the lookup table is standard quadratic interpolation.
But the lookup table takes the results of the quadratic interpolation and
improves it.  The lookup table was empirically generated by synthesizing some
pure sine waves of various in-between frequencies and just seeing what the
quadratic interpolation yielded.  Since I generated the original time series
signals, I know exactly what the right answer should be.  So here is the lookup
table:

float InterFFTLookup[] = {
-0.500,-0.496,-0.491,-0.487,-0.482,-0.478,-0.473,-0.469,
-0.464,-0.460,-0.455,-0.451,-0.446,-0.442,-0.437,-0.428,
-0.420,-0.411,-0.402,-0.393,-0.385,-0.376,-0.361,-0.345,
-0.330,-0.314,-0.288,-0.262,-0.228,-0.187,-0.119,-0.041,
0.000,0.041,0.119,0.190,0.230,0.264,0.290,0.315,
0.331,0.348,0.364,0.378,0.387,0.396,0.406,0.415,
0.424,0.433,0.440,0.444,0.448,0.453,0.457,0.461,
0.465,0.470,0.474,0.478,0.483,0.487,0.491,0.495};

Subsequent tests with real-world signals in a piano-tuning application have
shown this interpolator to be accurate in most cases to within 0.05 bins - more
accurate than the straight quadratic - but not so accurate that it would justify
going to a larger lookup table with greater resolution.  Perhaps more could be
done using amplitudes from other nearby bins.  But in my application the
real-world nature of those signals makes more distant bins much less dependable.
Also, my application focuses on frequencies that range from to 0.02 to 0.35 of
the Nyquist frequency.  Perhaps this interpolator is not so good outside of that
range.  But I would be interested in seeing what you people think of it.


Robert Scott
Ypsilanti, Michigan

On 28 Dec 2006 11:36:05 -0800, "Ron N." <rhnlogic@yahoo.com> wrote:

>..My guess is that the fractional bin offset table above will be
>very similar to a linear interpolation plus a lateral correction
>factor for the difference between the transform of your
>fft window (a sinc for a rectangular window) and a parabola.
>Thus the optimal table will vary depending on the window used,
>and should probably be recalculated for each different window
>function or shape.

FYI, the FFTs use in my development were 8192-point and 65536-point.  No
windowing was used.


Robert Scott
Ypsilanti, Michigan

Maybe you can help me. I've posted in this forum about an easy DFT or
FFT code. I need to know the amplitude and the phases of the sub-waves
of a waveform. Have you something of similar?

Robert Scott wrote:
> Here is an FFT bin interpolator I have been using with great success.  It is
> based on the standard quadratic interpolation formula:
>
> float CalcFFTInterpolation(float yPrecedingPower, float yMaxPower, float
> yFollowingPower)
>    // yMaxPower is the power amplitude in a local peak
>    // yPrecedingPower is the power amplitude of the preceding bin
>    // yFollowingPower is the power amplitude of the following bin
> {
>   float c = 2 * yMaxPower - yFollowingPower - yPrecedingPower;
>   if(c)
>   {
>     int index = (yMaxPower - yPrecedingPower) / c * 64;
>     if(index > 63) index = 63;
>       c = InterFFTLookup[index];  //  0...63 covers -.50 to +.49 bins
>   }
>   return c;
> }
>
> The returned value of the function is the bin index offset that you should add
> to the index where yMaxPower was found to get the interpolated frequency of this
> peak.  Everything up to the lookup table is standard quadratic interpolation.
> But the lookup table takes the results of the quadratic interpolation and
> improves it.  The lookup table was empirically generated by synthesizing some
> pure sine waves of various in-between frequencies and just seeing what the
> quadratic interpolation yielded.  Since I generated the original time series
> signals, I know exactly what the right answer should be.  So here is the lookup
> table:
>
> float InterFFTLookup[] = {
> -0.500,-0.496,-0.491,-0.487,-0.482,-0.478,-0.473,-0.469,
> -0.464,-0.460,-0.455,-0.451,-0.446,-0.442,-0.437,-0.428,
> -0.420,-0.411,-0.402,-0.393,-0.385,-0.376,-0.361,-0.345,
> -0.330,-0.314,-0.288,-0.262,-0.228,-0.187,-0.119,-0.041,
> 0.000,0.041,0.119,0.190,0.230,0.264,0.290,0.315,
> 0.331,0.348,0.364,0.378,0.387,0.396,0.406,0.415,
> 0.424,0.433,0.440,0.444,0.448,0.453,0.457,0.461,
> 0.465,0.470,0.474,0.478,0.483,0.487,0.491,0.495};
>
> Subsequent tests with real-world signals in a piano-tuning application have
> shown this interpolator to be accurate in most cases to within 0.05 bins - more
> accurate than the straight quadratic - but not so accurate that it would justify
> going to a larger lookup table with greater resolution.  Perhaps more could be
> done using amplitudes from other nearby bins.  But in my application the
> real-world nature of those signals makes more distant bins much less dependable.
> Also, my application focuses on frequencies that range from to 0.02 to 0.35 of
> the Nyquist frequency.  Perhaps this interpolator is not so good outside of that
> range.  But I would be interested in seeing what you people think of it.
>
>
> Robert Scott
> Ypsilanti, Michigan

My guess is that the fractional bin offset table above will be
very similar to a linear interpolation plus a lateral correction
factor for the difference between the transform of your
fft window (a sinc for a rectangular window) and a parabola.
Thus the optimal table will vary depending on the window used,
and should probably be recalculated for each different window
function or shape.

Another method I've read about to pre-correct for parabolic
interpolation is to use the transform of a parabola as the
fft window.


IMHO. YMMV.
-- 
rhn A.T nicholson d.0.t C-o-M

Robert Scott wrote:
> FYI, the FFTs use in my development were 8192-point and 65536-point.  No
> windowing was used.

By no windowing I assume you mean a rectangular window
(whose transform is a sinc function), since 65536 is a finite
number.

The main lobe of a sinc is a different shape from that of an
unside-down parabola; thus the need for your correction
factors.  If you do a least squares fit of a rotated sinc
directly to the complex data, you can use more points, and
perhaps get a more accurate or robust estimate.


IMHO. YMMV.
-- 
rhn A.T nicholson d.0.t C-o-M

"Ron N." <rhnlogic@yahoo.com> wrote in message 
news:1167337278.960522.51310@73g2000cwn.googlegroups.com...
> Robert Scott wrote:
>> FYI, the FFTs use in my development were 8192-point and 65536-point.  No
>> windowing was used.
>
> By no windowing I assume you mean a rectangular window
> (whose transform is a sinc function), since 65536 is a finite
> number.
As I understand it, FFT( Data, n) acts as if the n Data points are repeated. 
Windows are frequently used to minimize discontinuity between the end of 
Data and start of the next "copy" of Data.  If you are using a 65536 point 
rectangular window on a 65536 point data set, does it do anything?  What is 
the width of the sinc functions main lobe?

>
> The main lobe of a sinc is a different shape from that of an
> unside-down parabola; thus the need for your correction
> factors.  If you do a least squares fit of a rotated sinc
> directly to the complex data, you can use more points, and
> perhaps get a more accurate or robust estimate.
>
>
> IMHO. YMMV.
> -- 
> rhn A.T nicholson d.0.t C-o-M
> 

Philip Martel wrote:
> "Ron N." <rhnlogic@yahoo.com> wrote in message
> news:1167337278.960522.51310@73g2000cwn.googlegroups.com...
> > Robert Scott wrote:
> >> FYI, the FFTs use in my development were 8192-point and 65536-point.  No
> >> windowing was used.
> >
> > By no windowing I assume you mean a rectangular window
> > (whose transform is a sinc function), since 65536 is a finite
> > number.
> As I understand it, FFT( Data, n) acts as if the n Data points are repeated.
> Windows are frequently used to minimize discontinuity between the end of
> Data and start of the next "copy" of Data.  If you are using a 65536 point
> rectangular window on a 65536 point data set, does it do anything?  What is
> the width of the sinc functions main lobe?

His data samples are from a piano.  I'm not sure if the
most skilled pianist in the world could play a note exactly
65536 samples in length without requiring any windowing
of the sound waveform.

A dft does not "act" on its own.  It's only a complex vector
transform whose results are up to interpretation.


IMHO. YMMV.
-- 
rhn A.T nicholson d.0.t C-o-M

"Ron N." <rhnlogic@yahoo.com> wrote in message 
news:1167342623.830346.125920@a3g2000cwd.googlegroups.com...
>
> Philip Martel wrote:
>> "Ron N." <rhnlogic@yahoo.com> wrote in message
>> news:1167337278.960522.51310@73g2000cwn.googlegroups.com...
>> > Robert Scott wrote:
>> >> FYI, the FFTs use in my development were 8192-point and 65536-point. 
>> >> No
>> >> windowing was used.
>> >
>> > By no windowing I assume you mean a rectangular window
>> > (whose transform is a sinc function), since 65536 is a finite
>> > number.
>> As I understand it, FFT( Data, n) acts as if the n Data points are 
>> repeated.
>> Windows are frequently used to minimize discontinuity between the end of
>> Data and start of the next "copy" of Data.  If you are using a 65536 
>> point
>> rectangular window on a 65536 point data set, does it do anything?  What 
>> is
>> the width of the sinc functions main lobe?
>
> His data samples are from a piano.  I'm not sure if the
> most skilled pianist in the world could play a note exactly
> 65536 samples in length without requiring any windowing
> of the sound waveform.
>
> A dft does not "act" on its own.  It's only a complex vector
> transform whose results are up to interpretation.
True, but I think almost everyone here understands what I meant.

I guess my point was that I have trouble with your equating "no window" with 
a "rectangular window".  There is no 0->1 or 1->0 transition in "no window", 
it is a constant 1.  The ft is an impulse; the dft is a pulse.

>
>
> IMHO. YMMV.
> -- 
> rhn A.T nicholson d.0.t C-o-M
> 

Philip Martel wrote:
> "Ron N." <rhnlogic@yahoo.com> wrote in message
> news:1167342623.830346.125920@a3g2000cwd.googlegroups.com...
> >
> > Philip Martel wrote:
> >> "Ron N." <rhnlogic@yahoo.com> wrote in message
> >> news:1167337278.960522.51310@73g2000cwn.googlegroups.com...
> >> > Robert Scott wrote:
> >> >> FYI, the FFTs use in my development were 8192-point and 65536-point.
> >> >> No
> >> >> windowing was used.
> >> >
> >> > By no windowing I assume you mean a rectangular window
> >> > (whose transform is a sinc function), since 65536 is a finite
> >> > number.
> >> As I understand it, FFT( Data, n) acts as if the n Data points are
> >> repeated.
> >> Windows are frequently used to minimize discontinuity between the end of
> >> Data and start of the next "copy" of Data.  If you are using a 65536
> >> point
> >> rectangular window on a 65536 point data set, does it do anything?  What
> >> is
> >> the width of the sinc functions main lobe?
> >
> > His data samples are from a piano.  I'm not sure if the
> > most skilled pianist in the world could play a note exactly
> > 65536 samples in length without requiring any windowing
> > of the sound waveform.
> >
> > A dft does not "act" on its own.  It's only a complex vector
> > transform whose results are up to interpretation.
> True, but I think almost everyone here understands what I meant.
>
> I guess my point was that I have trouble with your equating "no window" with
> a "rectangular window".

Understand that "no window" on a finite length fft of real
world data (non-synchronized) is the same as a rectangular
window whose length is the same as that of the fft aperture.
Knowing this is a key to understanding some of what an fft
will do to an arbitrary frequency periodic waveform.

This window will transform any non-bin frequency sinusoid
into a sampled sinc.  (Actually it will do the same to a bin
frequency sinusoid.  But that sinc will be positioned so as
to have zeros at every bin location except one.  So the
sampling will make that sinc look the same as an impulse.)


IMHO. YMMV.
-- 
rhn A.T nicholson d.0.t C-o-M

john john wrote:
> Maybe you can help me. I've posted in this forum about an easy DFT or
> FFT code. I need to know the amplitude and the phases of the sub-waves
> of a waveform. Have you something of similar?

If I recall correctly, you asked for code and got it. I didn't check to 
see that the code was useful, but considering the source and unless it 
was intended to be a joke, it probably was. "Sub-wave" isn't a term I'm 
familiar with. I can guess the meaning, but it's better if you explain 
what you mean.

Jerry
-- 
Engineering is the art of making what you want from things you can get.
�����������������������������������������������������������������������