"Interpolating" to estimate exact frequency of peak from FFT results 2

[electricpete]

Vibration programs such as E-monitor have some algorithm built in that allows them to estimate the frequency of a peak to a finer resolution than the bin-width.

They use only the magnitude of the spectrum in the vicinity of the peak (no phase information is saved by E-monitor).

What kind of algorithm might be used for this purpose?

I am pretty sure that it must be either a spline or a best-fit line through the magnitudes in vicinity of the peak. The question is which one (fit or spline) and what form to use.

I tried it using a simple assumed quadratic form:
Y = c0 + c1*X + c2*X^2;
(where Y is magnitude and X is frequency)
With three points we can solve the three constants.

The three points are
(XLeft,Yleft), (XCenter,YCenter), (XRight,Right)

Where XLeft,XCenter and XRight are equally spaced (as would be the case for FFT frequency bins) and YCenter is assumed larger than XLeft and XRight

The excercize of taking the derivative of the polynomial and setting to 0 to find the peak is pretty trivial, but is given here:

http://home.comcast.net/~electricpete/eng-tips/InterpolatePeakFreq.pdf

A spreadsheet implementing this approach (just on three points demonstration purposes) is shown here:

http://home.comcast.net/~electricpete/eng-tips/InterpolatePeakFreqSpreadsheet.xls

My main questions:
1 – Does anyone know any other algorithm for accomplishing this task
2 – If using a polynomial spline, would the above quadratic form be correct? Or some other form?
3 - Does the quadratic form linked above seem reasonable?

=====================================
Eng-tips forums: The best place on the web for engineering discussions.

 

[GregLocock]

OK, for pure sinusoids. Now try it with a noisy signal. I haven't got the faintest idea what will happen.

Cheers

Greg Locock

SIGlease see FAQ731-376 for tips on how to make the best use of Eng-Tips.

 

[electricpete]

I will try it.

I certainly didn't mean to imply that all those decimal places would necessarily be meaningful. I would summarize by saying that the DTFT reconstruction method gives us the ability to obtain an arbitrarily fine resolution for estimating the frequency of peaks. It does not guarantee accuracy. And of course when two peaks are very close together, the situation is more difficult.

=====================================
Eng-tips forums: The best place on the web for engineering discussions.

 

[GregLocock]

Yes, that's exactly what I was getting at. In my time domain example you may be able to estimate how many cycles of the sinusoid are present to any degree of resolution, but the presence of other signals would rapidly mask that extra resolution.

I suppose another way of looking at it is that Fourier describes a possible harmonic structure for the time history, but there is no proof that it is the only harmonic structure with a 1:1 correspondence. (is there?)

Cheers

Greg Locock

SIGlease see FAQ731-376 for tips on how to make the best use of Eng-Tips.

 

[electricpete]

A revised spreadsheet:

http://home.comcast.net/~electricpete/eng-tips/InterpolatePeakFreqSpreadsheetTrial2.xls

I created an iterative routine to estimate the frequency of the peak to any arbitrarily specified resolution. (Function findpeakfreq in the vba editor).

I applied it to a sinusoid at 77.07626094... hz (irrational number) which had added broadband noise as Greg suggested. (sin (w*t) + rand()-0.5)

I applied my iterative routine with specified resolution of 0.001hz and the resulting frequency had an error of 0.11 hz (reducing the specified resolution further did not improve the results). So, the accuracy of our estimate in this case is limited by the characteristics of the input signal, not by our ability to resolve the output of the FFT.

Simple quadratic interpolation in this case gave an error approx 0.1hz (slightly better than the iterative...I think that was just luck and they are approx the same for this signal).

Still, considering that the binwidth is 4hz and you would get an error up to 2hz by just picking the highest bin center, I think 0.11 hz is a pretty good improvement in error.

A few more comments about the quadratic interpolation.

1 - While testing my iterative routine on the noiseless signal (Trail77Reconstruct), I initially got extremely inaccurate results when I specified a very tight tolerance. Further investigation revealed that the output of the iterative routine was fine, but the final quadratic interpolation was causing huge errors. The quadratic interpolation was predicting frequencies outside of the upper/lower frequencies of the input data!. In fact the quadratic routine as originally set up was very susceptible to the effect of roundoff errors when deltaF <<< F (the frequency spacing between the three points was a tiny tiny fraction of the actual frequencies). This was corrected in the quadratic interpolation routine of my newest spreadsheet by first performing the calculation as if the center frequency is 0, and then adding the actual center frequency to the output of that calculation.

2 – When examining the summary tab (all cases except for noisey case), we can see that the quadratic interplation gets better and better AS A FRACTION OF THE INPUT FREQUENCY SPACING when the input spacing is closer and closer around the true peak. (see graph on the right of summary table). i.e. if we decrease our input spacing by 10, we don't just get a factor of 10 improvement in our output, we get maybe a factor of 100 absolute improvement in performance of the quad interpolation. I believe this is because the true curve looks much smoother when we zoom into a small area around the peak and is much better fit by a quadratic polynomial. If we zoom out, the curve is more complex and not fit as well by a quad polynomial.


=====================================
Eng-tips forums: The best place on the web for engineering discussions.

 

[electricpete]

Before we leave this topic, here is my attempt at a summary: (open to any comments or suggestions)

The ability to resolve the frequency from the output of FFT is seemingly almost limitless if we have only a single-frequency input. The time windowing does have the effect of shifting the energy to the left and right in frequency, but that shifting is symmetric and does not change the center of the energy from the sinusoid.

When we have multiple frequencies (which includes noise), the time-windowed output of those frequencies may overlap. In this case our accuracy in estimating peak frequencies suffer and we may be unable to distinguish closely spaced frequencies.

The frequency spreading which occurs is proportional to the binwidth (inversely proportional to number of samples at a given sample rate). It also depends to a lesser extent on the type of window used.

If we decrease the binwidth by increasing the number of time samples (at a given sample rate), then there is less spreading and less overlap and a better ability to resolve frequencies and accurately estimate frequencies without interference from other signals.

Once we have fixed the number of time samples, there are three things discussed above that can increase the resolution in estimating a peak.

1 – Quadratic interpolation. The least computationally intensive and biggest bang for the computational buck. Also well suited to improve the output of one of the other methods.

2 – Zero-padding. More computationally intensive. Gives a higher resolution for all peaks on the whole spectrum. (Note).

3- Reconstruction of the DTFT. The most computationally intensive (in my implementation, anyway... I don't think this can be broken down into a simple convolution which could be processed by FFT multiplication). Gives unlimited ability to improve the resolution (Note).

Note - The accuracy of frequency estimation under 1, 2, 3 are all limited by the amount of overlap of adjacent frequencies... that aspect can only be improved by increasing number of samples or altering the window type.


=====================================
Eng-tips forums: The best place on the web for engineering discussions.

 

[sreid]

Now a little late but I happened on this (at a familar website) while I was looking for something else.

http://www.dspguru.com/howto/tech/peakfft.htm

 

[electricpete]

Interesting. They mentioned the quadratic method. And then several other methods that use a similar interpolation algorithm... just a little more sophisticated than the quadratic method. No mention of zero-padding (2 above) or reconstruction (3 above), which I believe are among the most precise methods (also computationally intensive).

=====================================
Eng-tips forums: The best place on the web for engineering discussions.

 

[VisiGoth]

All very interesting. But the "best" answer is really quite simple. The FFT is providing maximal compaction of data already. Do not zero pad. Do not window. Those techniques are for their own purpose. If you do neither you can use a zero bias answer that is also lowest variance.
Simply linearly interpolate between the largest picket and its largest adjacent neighbor. If the frequency source is truly a tone this is the "best". If you window then other ideas are better, but only approach this one in quality. The most commonly used one is the parabolic fit. It is robust, works with many window and padding methods and works with spectral densities (not perfect tones).

Try all of the advice you have been given in this thread for bias. To do this make up some test data. Use a pure tone of the entire duration of the sample. Add no noise. Run a loop over maybe 1% of your bin width. You know the correct answer and you will find what each of your techniques give as answers. The zero bias method I gave will give an exact answer to within numerical limits. If you use the power instead of the voltage magnitude of the pickets or many other methods you will see a bias curve that has a cyclic nature. You might be surprised just how bad some of those home brew methods are with no noise!

Now try the same thing by adding only Gaussian noise. If you want to be clever you can subtract out the bias curve you developed and you are left with only the noise component. You can then judge the various methods against these two metrics. You can make a single metric say by taking the max abs of the bias and add two sigma and see what is best.

Now, if the signal is not a tone and the noise is not Gausssian and there are nonlinear issues then other methods might be best for you. These other things are why you might use the simple parabolic fit method.


jsolar

 

[electricpete]

If I understand you correctly , you're saying that "linear interplation" can beat any of the methods discussed above for estimating the frequency based on a set of samples of a pure sinusoid with no noise.

I do not agree. Go to the following link:

http://home.comcast.net/~electricpete/eng-tips/InterpolatePeakFreqSpreadsheetTrial3.xls

Look at tab labeled Trial77NWReconstruct

It is a 77hz sinusoid input sampled at 2000hz for 512 samples.
The samples are given in column D. No window is applied. (or rectangular window depending on your terminology preference).

The Fourier results are in columns F / G /H. Specifically Column F is frequency, Column G is the complex Fourier coefficients, and column H is the magnitudes.

The values of column F / G / H in the vicinity of the peak are as follows:
70.3125 31.2239842348268+22.0919489360351i 0.149410432
74.21875 72.5199163858084+54.5972239034363i 0.354587574
78.125 -173.19941295231-138.434693741231i 0.866115267
82.03125 -37.4162769033782-31.6873355665261i 0.191528382
85.9375 -20.3504449522914-18.2295011520029i 0.106723963

Method 3 above (reconstruction of DTFT from DFT using the algorithm linked in my 6 Feb 08 0:26 message) gives the following estimate: 76.99679064 hz. (actual value is 77 hz).

(Note that method 2 will give the same results as method 3 if we pad the FFT long enough.)

Please provide an analysis of the above data using linear interpolation.

Or else provide your own example analysed by linear interpolation with either TWF or FFT data available so I can analyse the same data using method 3.

Thanks

=====================================
Eng-tips forums: The best place on the web for engineering discussions.

 

[hacksaw]

pete,

are the spreadsheet functions accurate enough for this?

I have experienced some surprizing levels of roundoff/truncation in excel for less complex computations.

 

[electricpete]

The function indpeakfreq that fills in cell L4 accepts a tolerance argument (it stops searching for the peak when the interval gets down to that tolerance). For large tolerances, the error is close to the tolerance. For small tolerances, I no longer come close to the tolerance. For example in the link above I specified 0.000001 tolerance, and the result 76.99679064 is 0.003209364 off the actual. I attribute this to roundoff error.

I don't think that computing to such tolerances has much practical purplse. I did think it would help settle my disagreement with Visigoth. Personally, I think it should be obvious without need for any example calculation that reconstruction of the DTFT is the most accurate way to estimate the frequency. But obvious is in the eye of the beholder.

I will look forward to a response.

=====================================
Eng-tips forums: The best place on the web for engineering discussions.

 

[electricpete]

A few followup comments (I find it handy to keep my notes on this subject together in one place here in this thread)

1 - I spent awhile looking at sreid's link. It had a lot of results but not much discussion, details, proof. Most of the links were dead. Here is a page that discusses very much the same issues and same algorithms and provides links to detailed analysis / comparison:

http://www.ericjacobsen.org/fe.htm

2 - If the time samples are available, there is another way to accomlish "Reconstruction of the DTFT" (where the DTFT is a continuous function of frequency whose samples at discrete points are given by the DFT/FFT). We simply use
X(w) = Sum(x[n]*exp(-j*w*n]
where w is the trial value of the frequency variable of the DTFT.

Not only is it much computationally less intensive than what I linked above from Dr. Mitra, but for me it is quite a bit more intuitively satisfying because we are simply constructing the DTFT from it's definition.

Another intuitive benefit is that we can easily see the relationship that the DFT/FFT are samples of the DTFT by comparing the above to the definition of DFT:
X(k) = Sum(x[n]*exp(-j*n*(2*pi*k/N)]

It reminds us that the DFT X(k) is simply samples of the DTFT X(w) at discrete points w = (2*pi*k/N) If you look only at those samples (DFT), then you are looking through a picket-fence and seeing only portions of the true DTFT / spectrum.

=====================================
Eng-tips forums: The best place on the web for engineering discussions.