From: Thomas Kragh <tkragh@NOSPAMeecs.umich.edu>
Subject: Re: What's the best acceleration approximation for discrete data?
Date: Sat, 24 Jul 1999 01:51:07 GMT
Newsgroups: sci.math.num-analysis
Keywords: filtering noise from data

Rob Jones wrote:

> Hi,
>
> Has anyone got any acceleration approximation methods, given a discrete
> position output, which are less prone to stochastic noise than the
> difference method which I'm currently employing?

If you can characterize your noise, and your system is (reasonably) linear,
the Kalman Filter is your best bet.  Actually, it's the optimal solution.


> (I'm trying to damp down
>  resonance in a hyper-accurate postioning system using acceleration
> feedback, but my acceleration approximation has too much noise riding on
> it)

For real time control with feedback, any phase delay will kill you.  So,
high-order models / polynominal fits / etc. are out, since there is phase lag
inherent with waiting for all those samples to come into before your filter
spits out an estimate.

>
>  current method is:
>
>  acceleration(estimate) = x(n+2) -2x(n+1) + x(n) / Ts^2
>

Standard formula for the 2nd difference.  However, I think you mixed your
+/-.. For real-time stuff, it should be x(n-2) -2x(n-1) + x(n) / Ts^2, since
you're looking at past values.  The above finite-difference scheme is
actually a discrete-time filter in its own right.  Take the z-transform, you
get

H(z) = 1 - 2z^(-1) + z^(-2)

And you can get the frequency response in matlab by

>>B = [1 -2 1];
>>freqz(B, 1);

You'll notice a gain of +40db/dec (as expected for a 2x-derivative), and a
linear phase response (as expected for an FIR filter).  Note, an ideal
2x-derivative will have  a frequency response of

H(w) = -w^2

for a gain of +40db/dec, and phase = -180deg (constant, not linear phase).


Now, for a few comments on other people's suggestions:

     One of my favorites, although it may not be computationally
     reasonable or even needed for your problem, is to Fourier transform

     the data, zero out the small, high frequency terms, and do an
     inverse transform with the rest.

UGH! NOT a good thing to do, unless you need a quick & dirty fix.  Hard
Clipping in the frequency domain will lead to sinc(t) ringing in the time
domain, due to the Fourier x-form pair relationship between the unit pulse
and sinc function.

     Using a running-mean filter (i.e. a "square window": w(j) = 1/N, N:
     window
     length ) in the time domain would give you sidelobes in the
     frequency
     domain (the Gibbs phenomenon). Using a cosine window (w(j) =
     1/2(1-Cos[2
     \pi j/(N-1)]), j=0,..,N-1) reduces these sidelobes much better.


Better idea for a quick & dirty filtering.  There are a whole bunch of
windowing schemes for filter design, of which cosine filter is one.  Pick one
like Cosine or Hanning.  If you have an idea of your noise bandwidth, simply
set your filter cross-over to about that, and pre-filter the noisy input
signal before differentiation.

     Then you should be able to translate this into state
     space and design an observer.  If you want to get really fancy
     about it, you
     can even set it up as a Kalman filter or Wiener filter to optimize
     noise
     response....

Now we're getting closer.  Ideally, a discrete-time Kalman filter will be
your best best.  FYI, a Weiner Filter is simply a special case of a Kalman
filter for stationary noise statistics.

Last comments.  If you have stationary noise, whose statistics can be defined
(i.e. do you have a power spectrum measurement available?), do the kalman
filter.  If you only have a rough idea of the noise bandwidth, then form a
discrete filter that low-pass-filters and then differentiates.

p.s.  Also try the comp.dsp newsgroup.  They live for this kind of stuff.

Tom Kragh
remove obvious spam-block to reply.
[HTML portion deleted -- djr]