"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.

 

[electricpete]

by the way "h" is my notation for the bin width (spacing between X values)

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

 

[IRstuff]

The algorithm should be based on a sinc function, similar to what's done in image processing with the Lanczos interpolation.


TTFN

FAQ731-376

 

[GregLocock]

Look up the picket fence shape for your particular window, then you should be able to see how to build an interpolator. You could even compare the frequency estimate for two different windows, or no window at all (uniform).

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

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

 

[electricpete]

IRStuff - Great point. I remember that if we sample a bandlimited continuous signal below it's Nyquist frequency, then we can reconstruct the original continuous signal from the samples by a sum of weighted/shift sync functions. It makes sense that there would be a dual property where we can reconstruct the continuous DTFT from it's discrete samples (the DFT). I'll give that a try.

Greg - I'll have to think about that.

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

 

[electricpete]

With IRStuff's comment to focus my search, I was able to find the formula using google... slide 39 here:

http://web.engr.oregonstate.edu/~luca/ece464-564/LecturesPPT/Ch3_2.ppt

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.

 

[electricpete]

To Greg's comment about the uncertainy principle, I think that once the continuous frequency function is reconstructed, the peaks will still be somewhat spread out by the time windowing that is assumed somewhere along the line. The longer the sample, the narrower the peaks (this reflects the uncertainy principle). But still I think we can develop an estimate of the frequency corresponding to the center/top of the fattened peak with much more precision/accuracy than the bin width.

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

 

[GregLocock]

If we go back to the time domain, FFT says I can't use the information that there are (for example) 76 complete sine waves in my 1 second long frame of data, compared with 76.25

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

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

 

[IRstuff]

Do you have an example of the interpolation done by the program, i.e., the raw data and the answer it came up with?

TTFN

FAQ731-376

 

[MikeyP]

I'm no expert on anything other than DFT signal processing, but the MUSIC-type algorithms may be able to help you out here.

M

--
Dr Michael F Platten

 

[electricpete]

Thanks Mike - I'll look at that if I get a chance.

Here is the raw data including frequencies and magnitudes for vibration velocity on an electric motor.

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

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://home.comcast.net/~electricpete/eng-tips/InterpExample.jpg

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://home.comcast.net/~electricpete/eng-tips/InterpolatePeakFreqSpreadsheet7200.xls

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

 

[electricpete]

Actually, looking at the format of the output raw , I notice the very first bin is labeled 67.5 which would be the right boundary of the bin, not the center of the bin. If the numbers listed were the center of the bins, the output should have centered bin frequency labels: 0.5*(67.5), 1.5*67.5, 2.5*67.5 etc instead of 1*(67.5), 2*67.5, 3*67.5.

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.

 

[electricpete]

I must be missing something about how the frequency numbers relate to the bins. It doesn't make sense that they would label the bins with a frequency corresponding to the right boundary of the bin. I'll look at it some more tonight.

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

 

[electricpete]

ok, I think I have it figured out.
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.

 

[electricpete]

The zero frequency bin is -0.5*BW -> 0.5 * BW. It perhaps is discarded.

Correction, that bin information is captured in the bin labeled 54000 = 800 * 67.5 which is the center frequency (N/2)*BW for the bin spanning (N/2-0.5)*BW -> (N/2+0.5)*BW (note N is number of samples)


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

 

[IRstuff]

This is one of those problems that have multitude of "answers," all of them partially correct ;-)

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

 

[sreid]

I don't know if this is any help (or applicable)but "Zero Padding" collected data and doing the FFT on this larger data set is sometimes cslled "Sharpening." An example would be, collect 64 data points, pad with 64 zeros, do a 128 point FFT. This doubles the number of frequency bins.

 

[WCFoiles]

Do you think they incorporate the phase values in the calculation? These are available to those who calculate the fft.

Regards,

Bill

 

[electricpete]

I have done some review and fiddling around with numerical examples and I think I now know more than I ever wanted to know about estimating frequency of a peak from the FFT.

The following spreadsheet compares various methods for estimating the "exact" frequency of a peak from FFT spectral results.

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

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.engr.oregonstate.edu/~luca/ece464-564/LecturesPPT/Ch3_2.ppt

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.

 

[electricpete]

By the way Bill - the last method (reconstructing DTFT) uses phase. As I mentioned my earlier post (6 Feb 08 0:26), I'm sure Emonitor does note save the phase data for the FFT. I suspect they use something similar to the quadratic interpolation.

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