FFT / Inverse FFT producing small compounding error

[ChrisLudwig]

Hello All!

I am using an FFT and inverse FFT as a sort of band pass filter to inspect the vibratory content between a low cutoff frequency and a high cutoff frequency.

My reason for the use of FFT and not other bandpass filtering methods is performance.  I am running filtering and limit checking accels in a control room and my program consumes 25% of the CPU on our higher end workstations.  I am afraid that if I gave our users the ability to specify high ordered digital filters that I'd consume more than 50% of the CPU.  I've found that can't push a windows machine faster than about 50% or latency becomes a problem and I start filling buffers and dropping frames of data.

I have a test signal which is 5g at 30Hz (we think in half peak to peak amplitude, not RMS at my company).  The signal is acquired at ~3145Hz and came from a lab calibrated accel mounted to a calibrated portable shaker.  When I run the test signal through my algorithm, I find it has a compounding error of some sort, like the FFT is not initializing correctly or something.  Often times, I see a definite repeating pattern in the deviation from the correct amplitude.

In the attached vi, I have a 13 second data sample that should have a amplitude of 5-5.03g, but instead, the inverse FFT produced signal has an amplitude stepping up to 5.6g, well over the envelope of the source waveform.  Looking at the inverse FFT reconstructed signal, there appears to be definite artifacts at the beginning and end of each analyzed segment of the signal.  That is why I am thinking about a problem with initialization of state information or reentrancy, but I can't see how I am doing anything wrong in that regard.

My target accuracy for stable signals such as this is less than 5% error with a goal of getting even better. Once a possible exceedance is identified, the user has the option to pass the data through a better algorithm to confirm or reject the exceedance, but I can't keep flagging false positives as each must be cleared by the operator. 

Can anyone advise what I might be doing wrong? 

[altenbach]

Your segments don't contain integer numbers of cycles of your fundamental frequency (most easily seen by looking at the mean of the segments, which is never zero in your case). This means you get spectral leakage and all kinds of artifacts.

Why don't you filter the entire data at once?

[ChrisLudwig]

I am running tests all day long with a need to keep updating the operator's display at around 1Hz.  I have to chunk out the data and process in pieces.  Trying to get integer cycles of the fundamental frequency would also be a problem as the shaft speeds sweep up from zero to the design speed during start up and shutdown.  Once operating, some shafts will still change speed by 10% or more.  Other shafts have no design speed at all and vary by 50% in operation.  Still other accels are in places where they see vibrations sourced from 2 or 3 different sources and therefore have 2-3 fundamental frequencies, none of which ever coincide.

In a few cases I have shaft azimuth signals, but that is a minority of cases.  I can usually detect the strongest fundamental frequency, but not always as the fundamental frequency is not always the highest magnitude on an FFT.  I don't think that I can identify integer cycles on the complex waveform when I get really complicated waveforms. 

My FFTs for spectrum visualization use Flat Top windows and do not suffer from any identifiable leakage errors.  Do I have similar options for the inverse FFT?

[Kevin_Price]

I'm not an expert in these things but have always thought that direct application of most any "reasonable" time domain filter is more computationally efficient than taking the scenic route into the frequency domain and back.

What kind of data bandwidth are you dealing with (sample rates, channel counts, etc.)?   How did you do the time-domain filtering in LabVIEW?  The simple filtering functions might be unnecessarily recomputing the coefficients every call.

-Kevin P

[ChrisLudwig]

We are typically looking at about 35 accelerometers at a max acquisition rate of 3145Hz.  Most run at that rate, but a few run slower.  Our acquisition system does not give us very much choice on rates.  For each accel, we typically track magnitudes with tracking filters at 2-4 different frequencies.  These normally track the primary fundamental frequencies of the system in part that the accel is installed on.  We also normally have an average of 1 band pass filter per accel.

I am running an acquisition rate of 3145 on ~35 accels with a total demand of ~105 tracking filters giving amplitudes at a discrete frequency and ~35 band pass filters over some defined static frequency band.  Behind the scenes, I do two FFTs on every accel, one with a flattop window.  That FFT is used to extract magnitudes at frequencies of interest.  I do another FFT with a user selectable window since they almost always insist on looking at a Hanning windowed FFT.  I then do another unwindowed FFT plus inverse FFT for each desired band pass filter.  So, summing that up, I'd expect to be doing 105 FFTs per second and 35 Inverse FFTs per second with the current usage of the program.

The program I made was designed to replace a hardware in the loop offering that primarily provided FFTs on screen and a waterfall or 2D color map time history.  To that end, this application started life as a FFT program with 2D colormap which then automated limit tracking at discrete frequencies by observing the FFT output and then used FFT + inverse FFT for band pass calculations.  Never take a customer spec too seriously though.  After I made it, they decided they didn't really need to be looking at the FFT displays all the time and said they'd prefer a time history plot of all tracked limits and would rarely need to look at the FFT display.

So to the big questions:

I am a noob at this.  My background was in onboard software for pilot displays that do not perform a dynamics function.  I was asked to take on this project because I could complete it in about 1/3rd of the time that C++ programmers were quoting.  I'll run some test cases comparing time domain filters against the current FFT filters to see how the performance compares, but what exactly is a "reasonable" time domain filter?  I'll be looking at 3, 5, and 7 pole butterworth filters and probably a tone extraction function as well to get vibes at discreet frequencies, but what is the optimum number of samples for those filter types?  For an FFT, it is easy.  If I want to see 0-1500Hz, I need to sample at 3000Hz (Nyquist).  It therefore works out that about 1 second of data should be analyzed for an FFT, unless you are intentionally oversampling to get more solutions per second.  A better way to say it is to analyze a sample count that is twice your highest anticipated frequency, then round that up to the next highest power of 2 for computational efficiency.

[Kevin_Price]

The modest-order Butterworth filter you mentioned falls *easily* into the realm that I'd consider reasonable.  A high order filter (dunno an exact threshold, but am suspicious when getting into very far into double digits) often gives diminishing returns for its extra processing demands.

These are just crude personal rules of thumb, no extensive theory or broad-ranging experience.  You're pretty deep in your project and can just pay attention to the results of your experiments.

As to the specific question of optimizing the # samples to pass to these algorithms: the time domain filters aren't so sensitive the way an FFT is.  They maintain internal state so that feeding data in chunks still causes them to produce output as though you fed it the entire continuous data set.  Just don't set "init=True" every time.

Also, as mentioned before, look for the filter forms that let you calculate coefficients once and then simply apply them repeatedly rather than re-calculating every iteration.  That'd be a potential source of inefficiency.

-Kevin P

[ChrisLudwig]

I understand what you are saying about butterworth filters with the coefficients and preserving state from one execution to the next.

I did find this though regarding the same for FFTs:

Calling Mathematics and Signal Processing VIs from a Reentrant VI

What state information is an FFT or inverse FFT trying to preserve?  I could only imagine phasing info, not magnitude info.  I see that both can be re-initialized, but I don't know why.  I have my FFT-Inverse FFT pair inside a loop running through multiple parameters so I am definitely not preserving state on those two functions.

[Kevin_Price]

I honestly don't know what state information the FFT function would preserve.  Maybe internally it has some dependency on the point-by-point FFT function which *does* need to preserve state?   Maybe it stores the windowing function values for re-use if subsequently called with same # samples & window type?   Dunno.

-Kevin P

[jamiva]

There is no compounding error with the FFT. It is doing exactly what it is supposed to do. An FFT assumes that the waveform is PERIODIC. If you do not have an integral number of cycles, you will get discontinuities between the end of one period and the beginning of the next. As altenbach said "you get spectral leakage and all kinds of artifacts"

 "Often times, I see a definite repeating pattern in the deviation from the correct amplitude."

"Looking at the inverse FFT reconstructed signal, there appears to be definite artifacts at the beginning and end of each analyzed segment of the signal."

Here is an example of doing an FFT=>bandpass Filter=>inverse FFT of a simple sine wave pattern with discontinuities before and after filtering. (p.s., FFT was done on a 2048 pt sign wave segment. The plots below show 3 concatenated segments )

I would go with the time domain filtering. Chunking up your data in groups of 2048 points will almost certainly give you discontinuities.