"Interpolating" to estimate exact frequency of peak from FFT results
"Interpolating" to estimate exact frequency of peak from FFT results
(OP)
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://hom e.comcast. net/~elect ricpete/en g-tips/Int erpolatePe akFreq.pdf
A spreadsheet implementing this approach (just on three points demonstration purposes) is shown here:
h ttp://home .comcast.n et/~electr icpete/eng -tips/Inte rpolatePea kFreqSprea dsheet.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?
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://hom
A spreadsheet implementing this approach (just on three points demonstration purposes) is shown here:
h
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.





RE: "Interpolating" to estimate exact frequency of peak from FFT results
=====================================
Eng-tips forums: The best place on the web for engineering discussions.
RE: "Interpolating" to estimate exact frequency of peak from FFT results
TTFN
FAQ731-376: Eng-Tips.com Forum Policies
RE: "Interpolating" to estimate exact frequency of peak from FFT results
I'd add that this apparently breaks the Heisenberg uncertainity rule 1=BT, it may be OK because we know something that the FFT doesn't - it really is a stationary signal, with good S/N.
I've also wondered if some sort of heterodyning approach might be possible, in the time domain.
Cheers
Greg Locock
SIG:Please see FAQ731-376: Eng-Tips.com Forum Policies for tips on how to make the best use of Eng-Tips.
RE: "Interpolating" to estimate exact frequency of peak from FFT results
Greg - I'll have to think about that.
=====================================
Eng-tips forums: The best place on the web for engineering discussions.
RE: "Interpolating" to estimate exact frequency of peak from FFT results
http://web
It converts from the DFT / FFT (function of discrete frequency variable k) to the DTFT (function of continuous frequency variable w.... w goes from 0 to pi as k goes from 0 to N/2).
I guess if I know the max lies in a certain interval I can perform the calculation iteratively at frequencies within a smaller and smaller interval to pinpoint more exactltly the frequency of the peak.
Still I wonder if there is an easier way. I am positive the Entek's Emonitor database (version 3.1) does not store the phase information for an FFT. The hardware / software was bult at a time when memory was apparently at a premium and I think they felt a condition montoring vib analyst only needed the magnitude, not the phase. (time waveform is not saved unless you set up to save it... separate from saved fft).
But yet Entek seems to do a very good job at their peak estimation (somehow without phase). Any easier ways to do it?
=====================================
Eng-tips forums: The best place on the web for engineering discussions.
RE: "Interpolating" to estimate exact frequency of peak from FFT results
=====================================
Eng-tips forums: The best place on the web for engineering discussions.
RE: "Interpolating" to estimate exact frequency of peak from FFT results
But, if I know that the signal is sinusoidal then I can be confident that the two frequences that the two frequencies are 76 and 76.25
All Fourier says is, well, maybe.
Cheers
Greg Locock
SIG:Please see FAQ731-376: Eng-Tips.com Forum Policies for tips on how to make the best use of Eng-Tips.
RE: "Interpolating" to estimate exact frequency of peak from FFT results
TTFN
FAQ731-376: Eng-Tips.com Forum Policies
RE: "Interpolating" to estimate exact frequency of peak from FFT results
M
--
Dr Michael F Platten
RE: "Interpolating" to estimate exact frequency of peak from FFT results
Here is the raw data including frequencies and magnitudes for vibration velocity on an electric motor.
h
The motor is expected to have a vibration peak very near 7200 cpm (twice line frequency).
The bin width is 67.5 cpm. The data points in the vicinity of 7200 cpm are:
f (cpm) v (ips)
7155 0.05894788
7222.5 0.06896067
7290 0.02033701
The Entek E-monitor program calculates / labels the peak as 7195.7, which is pretty good (certainly better than 7222.5) as shown here:
http
I did try cutting/pasting that data into my quadratic interpolation spreadsheet (linked below). The frequency comes out to be 7200.276 cpm which I believe is probably a better estimate than 7195.7. I'm not sure if that is dumb luck or not. I do have the feeling that the true continuous curve will resemble a "smooth" curve passing through the discrete data. Here is the quadratic spreadsheet with this data filled in:
http:
=====================================
Eng-tips forums: The best place on the web for engineering discussions.
RE: "Interpolating" to estimate exact frequency of peak from FFT results
So if the frequencies listed in the spreadsheet really represent a frequency label identifying the right side of the bin (instead of center of the bin), then the centered frequencies should be 67.5/2 lower and my quadratic estimate would be 67.5/2 lower... and way off. Maybe it was just dumb luck that I got 7200 the first time. Hmmm.
=====================================
Eng-tips forums: The best place on the web for engineering discussions.
RE: "Interpolating" to estimate exact frequency of peak from FFT results
=====================================
Eng-tips forums: The best place on the web for engineering discussions.
RE: "Interpolating" to estimate exact frequency of peak from FFT results
The zero frequency bin is -0.5*BW -> 0.5 * BW. It perhaps is discarded.
The next bin is 0.5*BW -> 1.5 * BW. It is labeled based on its center 1*BW
The next bin is 1.5*BW -> 2.5 * BW. It is labeled based on its center 2*BW
That is kind of obvious, but somehow I got myself confused.
So, disregard my posts 6 Feb 08 9:23 and 6 Feb 08 9:34
So back to the quadratic method seems to have worked pretty well... again still may be dumb luck.
I may try some other examples. Next time I will pick a 4-pole motor so there is less chance that the 2*LF might be polluted by a harmonic of running speed.
=====================================
Eng-tips forums: The best place on the web for engineering discussions.
RE: "Interpolating" to estimate exact frequency of peak from FFT results
=====================================
Eng-tips forums: The best place on the web for engineering discussions.
RE: "Interpolating" to estimate exact frequency of peak from FFT results
There are two other approaches that can be taken:
> star tracker slope intercept (not sure what it's really called). You take the slopes of both sides and find the intersection = 7199.8533
> centroid. find the center of mass of the 6 points around the peak = 7191.0810
The first approach is supposedly better, at least, that's what the star tracker people claim. They purposely defocus their telescopes to blur out the stars to cover multiple pixels and use the slope-intercepts as the actual location of the stars. It allows them to calculate star positions to better than 1/50th of a pixel when SNR>10
TTFN
FAQ731-376: Eng-Tips.com Forum Policies
RE: "Interpolating" to estimate exact frequency of peak from FFT results
RE: "Interpolating" to estimate exact frequency of peak from FFT results
Regards,
Bill
RE: "Interpolating" to estimate exact frequency of peak from FFT results
The following spreadsheet compares various methods for estimating the "exact" frequency of a peak from FFT spectral results.
http:/
In all cases the sampling rate was 2000 hz (0.0005 seconds between samples).
512 samples were used for 512 point FFT (except for zero-padded case).
The FFT output (more specifically, the graphed half of the output) is 256 points representing amplitude at frequencies 0..1000 hz
The bin width is 1000hz / 256 = 3.90625
Through experimentation, I found that the first bin corresponds to frequency centered on 0*binwidth (different than Emonitor output whose first bin is centered on 1*binwidth).
(By the way, E-monitor does in fact discard that first zero-frequency bin.... it is not reflected in at N/2 as I previously said.... it would be at N.... but the spectrum stops at N/2).
Several sinusoids were input at frequencies between 74 and 79.1 hz
The relevant frequency bin centers in that neighborhood are: 70.3125, 74.21875, 78.125, and 82.03125
Results of estimating the frequencies from the FFT are shown in the summary tab (the tab all the way to the right)
The first five rows of the summary (tabs Trial74 through Trial 79.1) use a Hanning-windowed input at various frequencies. The quadratic interpolation to estimate frequency gives a ~ 1.5– 5% error as fraction of binwidth.
The next two rows/tabs/trials use no windowing (or "rectangle windowing" if you prefer). The quadratic interpolation to estimate frequency gives a ~ 15 – 20% error as fraction of binwidth.
The next two trials (tabs Trial77ZeroPad and Trial77.3ZeroPad) use zero-padding as suggested by sreid. (I applied the Hanning window first and the zero-padding afterwards). 1536 0's are added to the end of the 512 samples to give a total 2048 point FFT. This gives an effective binwidth which is a factor of 4 lower (although I have left the bin width the same in the table so we can compare errors as fraction of bin-width on an apples-to-apples basis). Then the quadratic interpolation was applied and the resulting frequency estimation errors were in the neighborhood 0.1% of original (nonzeropadded) binwidth.
I would have thought that for zero-padding plus quadratic interpolation, we would expect a reduction of 25% (for reduced binwidth) times a reduction of 1.5-5% (for quadratic interpolation) which would give a reduction to 0.375% - 1.25% error. We did much better than that (only 0.1% error). Apparently there is a synergistic effect from combining these two methods. I think maybe (?) the quadratic method works better when the amplitudes on either side of the highest amplitude are not too far down... gets less accurate the further down they are compared to the center.
The last method (tab Trial 77Reconstruct) was similar to IRStuff's suggestion. As discussed above, the DTFT (continuous function of frequency) can be reconstructed from the DFT/FFT (discrete function of frequency) using the method described on slide 39 here:
http://web
It is not too hard to implement in vba, but two things gave my problems initially:
1 – If you want to store a complex variable in vba (without putting it into a spreadsheet), you have to use a string variable. The complex functions like imsum take strings as input arguments and return strings as an output (nothing in the help says anything about that... but when I pulled a complex number from spreadsheet into a vba variant variable, the type became variant/string).
2 – The ratio of sin(x)/sin(x/N) needs special attention when x is 0, and this x (w*N-2*pi*k) can be 0 at certain k for certain w (w is tied to the frequency you want to interpolate at).
The function was finally programmed successfully in the spreadsheet as a function named InterpAmplitudeFromFFT which you can see using Tools/Macro/Visual Basic editor (note the quadratic interpolation is there also, under function interppeakfreq).
I reconstructed the DTFT at 0.02 hz increments (from the FFT at 4hz binwidth) and got a very smooth curve. Then I used quadratic interpolation to finish off the job based on the peak of that 0.02-hz resolution curve. This gave a predicted frequency of 77.00001015 against an actual frequency of 77. (0.00026% of original bin width error). I'm sure I could have done better if I wanted by using finer intervals or a recursive algorithm to narrow down the interval.
On a completely unrelated side note – It is interesting to compare the spectrum of Trial77 (Hanning window) and Trial77NW (no window). The top portion of the Hanning windowed peak is actually FATTER than the top portion of the unwindowed peak as you can see from the fact that the adjacent bin amplitudes are higher in the windowed than the unwindowed spectrum. (I think this is a reflection of the fact that the Hanning window has a wider main lobe than the rectangular window). Of course the base on the unwindowed peak is much fatter than the base of the windowed peak (this is a reflection of the fact that the rectangular window has much higher sidelobes than the Hanning window).
=====================================
Eng-tips forums: The best place on the web for engineering discussions.
RE: "Interpolating" to estimate exact frequency of peak from FFT results
=====================================
Eng-tips forums: The best place on the web for engineering discussions.
RE: "Interpolating" to estimate exact frequency of peak from FFT results
Cheers
Greg Locock
SIG:Please see FAQ731-376: Eng-Tips.com Forum Policies for tips on how to make the best use of Eng-Tips.
RE: "Interpolating" to estimate exact frequency of peak from FFT results
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.
RE: "Interpolating" to estimate exact frequency of peak from FFT results
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
SIG:Please see FAQ731-376: Eng-Tips.com Forum Policies for tips on how to make the best use of Eng-Tips.
RE: "Interpolating" to estimate exact frequency of peak from FFT results
http://
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.
RE: "Interpolating" to estimate exact frequency of peak from FFT results
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.
RE: "Interpolating" to estimate exact frequency of peak from FFT results
http://www.dspguru.com/howto/tech/peakfft.htm
RE: "Interpolating" to estimate exact frequency of peak from FFT results
=====================================
Eng-tips forums: The best place on the web for engineering discussions.
RE: "Interpolating" to estimate exact frequency of peak from FFT results
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
RE: "Interpolating" to estimate exact frequency of peak from FFT results
I do not agree. Go to the following link:
http://
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.
RE: "Interpolating" to estimate exact frequency of peak from FFT results
are the spreadsheet functions accurate enough for this?
I have experienced some surprizing levels of roundoff/truncation in excel for less complex computations.
RE: "Interpolating" to estimate exact frequency of peak from FFT results
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.
RE: "Interpolating" to estimate exact frequency of peak from FFT results
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.