On Thu, 28 Dec 2006 17:21:53 GMT, ---@--- (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; >}Hi Robert, darn, I don't "speak C language". Is there any chance I can talk you into describing your algorithm with plain ol' algebra? Thanks, [-Rick-]
Empirically Optimized FFT bin interpolator
Started by ●December 28, 2006
Reply by ●December 29, 20062006-12-29
Reply by ●December 29, 20062006-12-29
On Fri, 29 Dec 2006 15:05:44 GMT, R.Lyons@_BOGUS_ieee.org (Rick Lyons) wrote:>Hi Robert, > > darn, I don't "speak C language". >Is there any chance I can talk you into >describing your algorithm with >plain ol' algebra?OK, here it is. After doing an FFT, calculate the power spectrum for each frequency as the square root of the sum of the squares of the real and complex part of each element of the FFT array. Call this the power array. Search through the power array for a local maximum. Suppose you find this maximum at index "i" in the power array. "i" relates to a frequency as usual for FFTs. But we want to estimate the frequency of this peak to better than just the closest whole number index, so we calculate an offset, "c" that ranges from -.5 to +.5 to add to "i". That is the interpolation, which is based on the power amplitude in the neighboring elements of the PowerArray: yMaxPower = PowerArray[i] yPrecedingPower = PowerArray[i-1] yFollowingPower = PowerArray[i+1] c = 2 * yMaxPower - yFollowingPower - yPrecedingPower index = (yMaxPower - yPrecedingPower) * 64 / c index is rounded down to the next lowest whole number, ranging from 0 to 63 c = InterFFTLookup [ index ] So all you need now is the table of 64 elements called InterFFTLookup. Here they are for index 0 to index 63: -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 Robert Scott Ypsilanti, Michigan
Reply by ●December 29, 20062006-12-29
"Ron N." <rhnlogic@yahoo.com> wrote in message news:1167350167.822574.57430@n51g2000cwc.googlegroups.com...> 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.)I'm trying to see this. Could you show me an example where the magnitude of the transform of a non-bin frequency sinusoid decreases and then increases as the magnitude of a sinc would do? I've used the following code in Octave T = (0:8191)/8192; D = sin(2*pi*f*T); FT = fft(D,8192); plot(0:8191,abs(FT),'-') with f set to numbers like 1000.05, 1000.1, 1000.2. As far as I can see, the magnitude after the peak is monotonicly decreasing towards 4096. Thanks, --Phil Martel> > > IMHO. YMMV. > -- > rhn A.T nicholson d.0.t C-o-M >
Reply by ●December 30, 20062006-12-30
Philip Martel wrote:> "Ron N." <rhnlogic@yahoo.com> wrote in message > news:1167350167.822574.57430@n51g2000cwc.googlegroups.com...>> 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.) > > I'm trying to see this. Could you show me an example where the > magnitude of the transform of a non-bin frequency sinusoid > decreases and then increases as the magnitude of a sinc would > do? I've used the following code in Octave T = (0:8191)/8192; > D = sin(2*pi*f*T); > FT = fft(D,8192); > plot(0:8191,abs(FT),'-') > > with f set to numbers like 1000.05, 1000.1, 1000.2. > As far as I can see, the magnitude after the peak is monotonicly > decreasing towards 4096.Taking the absolute value of one value between every two zero crossings is sort of extracting the envelope, which is of course monotonic. I don't know Octave but presumably Re(FT) = 0 and if so, plotting Im(FT) should show you the sinc. Martin -- Values of beeta will give rise to dom!
Reply by ●December 30, 20062006-12-30
On Fri, 29 Dec 2006 15:05:10 GMT, ---@--- (Robert Scott) wrote:>On Fri, 29 Dec 2006 15:05:44 GMT, R.Lyons@_BOGUS_ieee.org (Rick Lyons) wrote: > >>Hi Robert, >> >> darn, I don't "speak C language". >>Is there any chance I can talk you into >>describing your algorithm with >>plain ol' algebra? > >OK, here it is. After doing an FFT, calculate the power spectrum for each >frequency as the square root of the sum of the squares of the real and complex >part of each element of the FFT array. Call this the power array. > >Search through the power array for a local maximum. Suppose you find this >maximum at index "i" in the power array. "i" relates to a frequency as usual >for FFTs. But we want to estimate the frequency of this peak to better than >just the closest whole number index, so we calculate an offset, "c" that ranges >from -.5 to +.5 to add to "i". That is the interpolation, which is based on the >power amplitude in the neighboring elements of the PowerArray: > > yMaxPower = PowerArray[i] > yPrecedingPower = PowerArray[i-1] > yFollowingPower = PowerArray[i+1] > > c = 2 * yMaxPower - yFollowingPower - yPrecedingPower > index = (yMaxPower - yPrecedingPower) * 64 / c > index is rounded down to the next lowest whole number, ranging from 0 to 63 > c = InterFFTLookup [ index ] > >So all you need now is the table of 64 elements called InterFFTLookup. Here >they are for index 0 to index 63: > >-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 > > >Robert Scott >Ypsilanti, MichiganHi Robert, Thanks much for going to the trouble of "spoon feeding" a description of your algorithm to me. I appreciate it. Your algorithm looks very similar, but not equal, to "spectral peak location" algorithms mentioned here by Robert Bristow-Johnson and Eric Jacobsen. However, you're using a lookup table and they didn't. My plans are to fire up my MATLAB software in the next few days and model your algorithm, Robert, to learn more about its behavior/performance. By the way, isn't it true that if you compute the square root of the sum of the squares of the real and complex part of each element of the FFT results, that will give you "magnitude" samples rather than "power" samples? So I wonder, ... are you computing FFT magnitude samples or power samples? Thanks again Robert, [-Rick-] Happy New Year!
Reply by ●December 30, 20062006-12-30
On Sat, 30 Dec 2006 14:53:27 GMT, R.Lyons@_BOGUS_ieee.org (Rick Lyons) wrote:>By the way, isn't it true that if you compute the >square root of the sum of the squares of the real >and complex part of each element of the FFT results, >that will give you "magnitude" samples rather than >"power" samples? So I wonder, ... are you computing >FFT magnitude samples or power samples?I'm not sure of the terminology. But I have always used the magnitude as I described it. Whether that is the same thing as power or not I don't know. Robert Scott Ypsilanti, Michigan
Reply by ●December 30, 20062006-12-30
"Martin Eisenberg" <martin.eisenberg@udo.edu> wrote in message news:1167489538.532703@localhost...> Philip Martel wrote: >> "Ron N." <rhnlogic@yahoo.com> wrote in message >> news:1167350167.822574.57430@n51g2000cwc.googlegroups.com... > >>> 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.) >> >> I'm trying to see this. Could you show me an example where the >> magnitude of the transform of a non-bin frequency sinusoid >> decreases and then increases as the magnitude of a sinc would >> do? I've used the following code in Octave >>T = (0:8191)/8192; >> D = sin(2*pi*f*T); >> FT = fft(D,8192); >> plot(0:8191,abs(FT),'-') >> >> with f set to numbers like 1000.05, 1000.1, 1000.2. >> As far as I can see, the magnitude after the peak is monotonicly >> decreasing towards 4096. > > Taking the absolute value of one value between every two zero > crossings is sort of extracting the envelope, which is of course > monotonic. I don't know Octave but presumably Re(FT) = 0 and if so, > plotting Im(FT) should show you the sinc. >Octave is similar to Matlab. Re(FT) will not be 0 in general, since there are not an integer number of cycles. Can you clarify for me what the zero crossings of the FT of the rectangular window with N points (8192 in this case) should be? Thanks, --Phil Martel> > Martin > > -- > Values of beeta will give rise to dom!
Reply by ●December 30, 20062006-12-30
Robert Scott wrote:> On Sat, 30 Dec 2006 14:53:27 GMT, R.Lyons@_BOGUS_ieee.org (Rick Lyons) wrote: >> By the way, isn't it true that if you compute the >> square root of the sum of the squares of the real >> and complex part of each element of the FFT results, >> that will give you "magnitude" samples rather than >> "power" samples? So I wonder, ... are you computing >> FFT magnitude samples or power samples? > > I'm not sure of the terminology. But I have always used the magnitude as I > described it. Whether that is the same thing as power or not I don't know.Isn't "power" magnitude squared? (The presumed digital impedance being 1.) Jerry -- Engineering is the art of making what you want from things you can get. �����������������������������������������������������������������������
Reply by ●December 30, 20062006-12-30
"Jerry Avins" <jya@ieee.org> wrote in message news:jJSdnSH-TLqZCAvYnZ2dnUVZ_q3inZ2d@rcn.net...> Robert Scott wrote: >> On Sat, 30 Dec 2006 14:53:27 GMT, R.Lyons@_BOGUS_ieee.org >> (Rick Lyons) wrote: >>> By the way, isn't it true that if you compute the square >>> root of the sum of the squares of the real and complex >>> part of each element of the FFT results, that will give >>> you "magnitude" samples rather than "power" samples? >>> So I wonder, ... are you computing FFT magnitude samples >>> or power samples? >> >> I'm not sure of the terminology. But I have always used >> the magnitude as I >> described it. Whether that is the same thing as power or >> not I don't know. > > Isn't "power" magnitude squared? (The presumed digital > impedance being 1.)How about if the thing you are analyzing is measured in units of "power" or "energy"?
Reply by ●December 30, 20062006-12-30
John E. Hadstate wrote:> "Jerry Avins" <jya@ieee.org> wrote in message > news:jJSdnSH-TLqZCAvYnZ2dnUVZ_q3inZ2d@rcn.net... >> Robert Scott wrote: >>> On Sat, 30 Dec 2006 14:53:27 GMT, R.Lyons@_BOGUS_ieee.org >>> (Rick Lyons) wrote: >>>> By the way, isn't it true that if you compute the square >>>> root of the sum of the squares of the real and complex >>>> part of each element of the FFT results, that will give >>>> you "magnitude" samples rather than "power" samples? >>>> So I wonder, ... are you computing FFT magnitude samples >>>> or power samples? >>> I'm not sure of the terminology. But I have always used >>> the magnitude as I >>> described it. Whether that is the same thing as power or >>> not I don't know. >> Isn't "power" magnitude squared? (The presumed digital >> impedance being 1.) > > How about if the thing you are analyzing is measured in > units of "power" or "energy"?I don't think that energy has real and complex parts. As to power, I have no idea what the square root of (watts^2 + vars^2) might represent. Jerry -- Engineering is the art of making what you want from things you can get. �����������������������������������������������������������������������
