Robert Penoyer has an M.S. in electrical engineering from the University of Southern California. He has 17 years experience in electronic circuit and system design and simulation. He has been programming in C for five years and working with alpha-beta filters for four years. Bob can be reached on GEnie at RPENOYER. or on Compuserve at 71603,1335.## Child of the Kalman Filter

## The Classic a-b Filter

where:

s= position

v= velocitya = multiplier on the position measurement error (0a1)

b = multiplier on the velocity measurement error (0b1)

T= time between measurements

nparameter value forthistime increment

pnpredicted parameter value forthistime increment

p(n+1) predicted parameter value for thenexttime increment

mmeasured parameter^ estimated value (corrected prediction)

Though these equations use the terms

positionandvelocity, any parameters that are related by the derivative operator (e.g.,d/dx) can be substituted. The equations form a recursive relationship and are used formally as follows:Given a new position measurement, an estimate is made of the true position and velocity for the current time increment using Equation 1 and Equation 2. The estimates from Equation 1 and Equation 2 are used to predict the position and velocity for the next time increment using Equation 3 and Equation 4. When the next time increment occurs, i.e., when a new position measurement is made, the process is repeated. This method requires that

Tbe constant or nearly so.The a and b terms give the filter its name. These fixed coefficients replace the optimized coefficients of the Kalman filter.

## Simplifying the a-b Equations

While fine as given, Equation 1 through Equation 4 can be massaged to make the programming effort simpler. However, beware that any time a derivative is used on noisy measurements, the derivative tends to multiply the effects of the noise. Because the velocity term in the a-b filter is related to the position term by the derivative operator, the velocity estimate can be very poor. Therefore, I assumed that the velocity estimate of this a-b filter would not be used to provide a velocity estimate outside the equations. That is, I used the a-b filter only to estimate position from some set of position measurements. (If you should want to use the velocity term as an estimator of true velocity, consider smoothing it with a second a-b filter.)

Because I was not going to use the velocity term outside the equations, I could afford to modify it. First, I rearranged Equation 2 by multiplying through by

T:I further modified Equation 5 such that

Twas absorbed intov.This is a valid simplification since I do not usevdirectly. Consider these facts aboutT:

Toperates only onv

Every place voccurs,Toccurs (with the understanding that Equation 4 is also multiplied byT)Given these facts, the absorption of

Tintovhas no effect on the predicted and estimated positions. Each of thevvariables is simply scaled byT. Consider, also, these facts aboutT:

Tis the time between position measurements

Tdepends upon the operating speed of the computer system (assuming its speed determines the time between measurements)The filter's characteristics (i.e., its predictions and estimations of

s) are not affected byT. That is, though a faster computer yields a faster filter, the damping of the filter is not affected. This will be made clear in the geometric presentation that follows. WhenTis absorbed intov, Equation 5 becomes:As noted above, Equation 6 simply causes

vto be scaled byT.To begin the recursive estimation sequence, the initial estimate of position, s^o, is taken to be the first measured position,

sm0. Set the initial velocity estimate, v^o, to 0 unless you have a good reason for doing otherwise. Given these starting conditions, the equations can be rearranged, placing Equation 3 and Equation 4 first. The indices in these two equations then becomepnrather thanp(n+1) and n-1 rather thannso that these equations are used to predict the parameters for the current time increment based on the previous time increment. (Note that, as originally formulated, the equations were used to predict the parameter value for the next time increment based on the current time increment.) Also, since Equation 4 provides no new information, it can be ignored — the old estimate of velocity becomes the new predicted velocity.Following the elimination of Equation 4 and the replacement of

vpnwith v^n-1 in Equation 6, the reformulation of the original set of equations and their operating sequence now becomes:When the above sequence begins, the entry point is Equation 7 with s^o set equal to the first value of

smn(i.e.,sm0) and v^o set equal to 0. This set of equations and conditions are the basic smoothing algorithm. The actual programming implementation adds one more step.The parenthetical term in both Equation 1 and Equation 8 is the difference between the measured position and the predicted position. This difference is the measurement error. Instead of calculating it twice, once for each equation, a single measurement error term can be calculated following Equation 7:

emn= smn- spnThus, the actual implementation of the a-b filter is:

This formulation is very simple. It lends itself to simple coding and to a simple geometric explanation of the a-b filter's operation.

## A Geometric Explanation of the a-b Filter

Figure 1 shows that, given the starting point, s^n-1, and the velocity, v^n-1, a straight line predicts the new position,

spn. This is the operation performed by Equation 7. Since v^ includes theTmultiplier, the units of v^ (i.e., Tv^) are distance, not velocity.In general, the measured position will not coincide with

spn. Equation 9 quantifies the measurement error. This error is due primarily to noise. However, some of the measurement error legitimately may be due to a change in velocity. Therefore, since some of the difference may be valid, all the difference cannot be ignored. Thus, a, a positive multiplier that generally has a magnitude smaller than 1, is used to add a fraction of the measurement error to the predicted position to improve the position estimate. Equation 1 yields the final estimate of position. The function of Equation 1 is illustrated in Figure 2.From the geometry of Figure 1 and Figure 2, it can be seen that if

vis varied inversely asTis varied, then none of the details in the figures (i.e., positions or slopes) will change. That is, if the data and measurements are speeded up equally, the details of the figures will not be affected. Hence, variations inTaffect only the filter's speed; damping is not affected.

## Relating b to a

Just as a determines the fraction of the measurement error applied to the position prediction, b determines the fraction of the measurement error applied to the velocity prediction. Equation 8 uses b for this purpose. Benedict and Bordner have derived a relationship between a and b that optimizes the filter's performance:

This relationship is based on the condition that the best estimate be provided for a system that is tracking noisy information that is changing at a constant velocity. This may or may not be satisfactory for your particular application.

Equation 10 yields a slightly underdamped system. Being underdamped implies that the system will "ring" somewhat. That is, if the data makes a sudden step in position, the estimate of position will overshoot the true position at least once before settling down to an accurate estimate. A more heavily damped system might perform better in your application. Benedict and Bordner offer an a-b relationship for a critically damped system. If their relationship is rearranged to solve for b, it will be found that damping is not critical. To get closer to critical damping, Equation 11 multiplies the rearranged relationship by 0.8.

Equation 11 yields a system that is nearly critically damped. (In a critically-damped system, the output agrees with the input as quickly as possible with no ringing.) By using Equation 11 instead of Equation 10, you can reduce — but not eliminate — the ringing response. Because a is never permitted to exceed the range 0a1, the quantity under the radical is never negative.

Consider Equation 10 and Equation 11 merely as reference points. As long as the limits 0a1 and 0b1 are not violated, you are free to select the values of a and b to suit your needs.

## A Programming Example in C

I begin by estimating measured data uncorrupted by noise. Figure 3 contains a profile of the data that begins at a steady state, ramps to a second steady state, then steps to a third steady state. (Normally, you would not filter clean, uncorrupted data like this for an estimated value, but I am using it to demonstrate the performance of the a-b filter against a ramp and a step.) The profile might represent a sequence of measured voltages, currents, distances, altitudes, volumes of storage tank contents, disk-drive head positions, highway traffic volumes, percentage of

xin a vat full ofy, etc.I wrote and compiled

ALPHABET.C(Listing 1) using the Borland C++, Version 3.1 compiler configured for MS-DOS. I verified it on a system running MS-DOS, version 5.0. The plotting functions assume a standard VGA monitor (640 x 480 pixels). If your monitor is different, modify the#definestatements forPLOTxandSTEADY_STATExby simple scaling.profile(Listing 1) generates the data profile of Figure 3.

ALPHABET .Cis a simple demonstration of the performance of the a-b filter for two sets of a-b values. In each case, I give a an arbitrary value of 0.25; thus, 25% of each measurement error is used as a correction to the predicted position. Over the first part of the program, b is determined by Equation 10, the so-called optimal value. (It may not be optimal for your application.) The second part of the program uses the b value determined by Equation 11, the so-called near-critically damped value.First, the program plots the data profile. Once it plots the profile, the program pauses until you press a key. After you press a key, the program plots the response of the a-b filter to the data profile over the profile. You will see slight overshoots as the profile begins and ends the ramp. You will notice a significant overshoot following the step. Pressing another key plots the filter's output by itself.

Press another key to produce the signal corrupted by noise.

gauss_noisegenerates a mean-zero, gaussian (normal) noise function. To do so, it uses the Box-Muller method (Press, et al. 1988, MathSoft, Inc. 1988). The noise is given an arbitrary standard deviation of s = 10. The mean is set to m = 0 to comply with the requirement that the mean of the corrupting noise be zero. If the mean were non-zero, the a-b filter's response would be offset by the mean value of the noise. Ingauss_noisethe variablemeanis retained, despite being set to zero, to clarify its part in the computation of the noise function.Continue pressing a key to run through all the combinations of clean signal, noisy signal, and smoothing provided by Listing 1.

## Filter Performance

Several figures illustrate the filter's performance. Figure 4 compares the original data profile with the resultant output due to the a-b relationship of Equation 10. The filter tracks the ramp well. It tracks the step well, too, except for the overshoot following the step. The ringing of the underdamped filter is apparent here. Figure 5 illustrates the data profile corrupted by noise. The noisy data remains the same throughout a single run of the program so that each step of the program operates on the same noise. The noise is different each time the program runs so that the noise you will see if you run it will be different from that shown in Figure 5. Figure 6 compares the filtered noisy data with the original, uncorrupted data. Keep in mind that the intent is to recover the original, uncorrupted data as accurately as possible, thus the comparison. Figure 7 and Figure 8 use the more heavily damped filter due to the a-b relationship of Equation 11.

You might find it useful to modify Listing 1 by supplying your own values for a and b. Just be sure that neither is permitted to be negative or to exceed a value of 1.0.

alpha_betaperforms the actual a-b filter algorithm. Three things should be said about this function. First, if it were used in real time, all variables could be scalar. No arrays would be needed. As each estimated position (or whatever) value was calculated it could be returned to the calling function. Only the most recent values of variabless_estandv_estneed to be preserved. Therefore, if the algorithm were to be modified to filter measurements one call at a time,s_estandv_estshould be made static variables.Second, keep in mind that the method used here has absorbed the term

Tvinto the single variable,v. If you want to make use ofvas well ass, be sure to dividev(actually,Tv) byT, the fixed period over which each measurement is made. As implemented, the b/Tterm of Equation 2 has been change to simply b, thus eliminating any need for a division. If you would rather preserve the true scaling ofv, you can absorbTinto the calculation of b, avoiding repeated divisions byT.Finally, you will have to make your choice of a based on the nature of the data you are attempting to smooth. A large a (near 1.0) results in very little smoothing and very little noise reduction. A small a (near 0.0) results in a lot of smoothing, i.e., a lot of noise reduction and very rounded estimates of data steps. Your choice must be based as much on your good judgment as it is on the data to be smoothed. Equation 10 or Equation 11 will usually provide a suitable b value, though you might want to experiment with this too. Let experience and testing be your guide.

## Tracking

Tracking, another important application of the a-b filter, differs from smoothing by centering a tracking window on the predicted position of the measured data as determined by Equation 7. When the measured data falls inside the window, that data is taken as a valid measurement and used to update the filter's position estimate. Data that falls outside the window is ignored.

When the data is ignored for failing to fall inside the track window, the calculation of Equation 9 cannot be completed. Instead, ignore Equation 9 and assign 0 as the error value. Then carry out the other calculations. Because the a-b filter incorporates a velocity term, you use it to predict where the data should have been. When the error value is zero, the velocity estimate of Equation 8 is not changed so that the velocity estimate remains constant. That is, the velocity estimate, based on the history of calculations up to the data that failed to fall within the track window, is used to update the position prediction of Equation 7 and, hence, the position of the track window. If the data falls inside the track window once again, the normal calculation of measurement error, Equation 9, is resumed and tracking updates continue. If the system goes through too many iterations without the data falling inside the track window, the program must declare loss of track. You must determine for your application what constitutes too many misses.

Once you determine that track has been lost, you need a search-and-acquire algorithm to reacquire the data. During acquisition, you might find it useful to scale a (and b) from a larger value (near 1.0) when the target is first acquired, down to a more suitable value (e.g., 0.25 or so) to maintain track. (Note that initially setting s^o = smo is equivalent to starting with an alpha symbol value of 1.0.) The down scaling can be accomplished over a few iterations. This strategy can speed the convergence of the position estimate to the true data value.

You will have to determine how wide the track window should be for your application. You will also have to figure out how to determine that track has been lost and how the data is to be reacquired.

## References

Benedict, T.R., and G.W. Bordner. July 1962. "Synthesis of an Optimal Set of Radar Track-While-Scan Smoothing Equations,"

IRE Transactions On Automatic Control,pp. 27-32. New York, NY: The Institute of Radio Engineers.Gelb, Arthur. 1974.

Applied Optimal Estimation,Chapter 4. Cambridge, Massachusetts: The M.I.T. Press.Lewis, Frank L. 1986.

Optimal Estimation: With an Introduction to Stochastic Control Theory,pp. 88-89. New York, NY: John Wiley & Sons.Meditch, J.S. 1969.

Stochastic Optimal Linear Estimation and Control, Chapter 5. New York, NY: McGraw-Hill, Inc.Mendel, Jerry M. 1973.

Discrete Techniques of Parameter Estimation, pp. 159-162. New York, NY: Marcel Dekker, Inc.Morris, Guy V. 1988.

Airborne Pulsed Doppler Radar,Section 12.7.1. Norwood, MA: Artech House, Inc.Press, William H., et al. 1988.

Numerical Recipes in C, pp. 216-217. New York, NY: Cambridge University Press.MathSoft, Inc. Summer 1988. "Industrial Engineering Prof. Uses MathCAD for Simulations, Statistics, and Teaching,"

MathCAD User's Journal, Vol. 2, No. 3, p. 4. Cambridge, MA: MathSoft, Inc.

## Figure 1

The last position estimate is used to predict the current position. The prediction is based on the last estimate of rate.

## Figure 2

A fraction of the measurement error is added to the current position prediction to improve the position estimate.

## Figure 3

Clean data profile

## Figure 4

Clean data profile with underdamped alpha-beta smoothing

## Figure 5

Noisy profile

## Figure 6

Underdamped alpha-beta smoothing of noisy profile with clean data profile

## Figure 7

Clean data profile with near-critically damped alpha-beta smoothing

## Figure 8

Near-critically damped alpha-beta smoothing of noisy profile with clean data profile

## Listing 1

Generates the data profiles of Figure 3 through Figure 8/* ALPHABET.C Demonstrate the alpha-beta filter with a particular data profile, with and without gaussian noise, and with different alpha-beta values. The graphic features here assume the presence of a VGA monitor. This code was written for Borland C++, Version 3.1, coded for MS-DOS. */ #include <graphics.h> #include <stdlib.h> #include <stdio.h> #include <conio.h> #include <alloc.h> #include <time.h> #include <math.h> #define PLOT1 120 /* range of x for first part of plot */ #define PLOT2 310 /* range of x for second part of plot */ #define PLOT3 430 /* range of x for third part of plot */ #define PLOTEND 640 /* end of x range */ #define STEADY_STATE1 120.0 /* first steady state y value */ #define STEADY_STATE2 400.0 /* second steady state y value */ #define STEADY_STATE3 70.0 /* third steady state y value */ #define STEADY_STATEBOT 480.0 /* bottom limit for steady state values */ #define ALPHA 0.25 /* arbitrary alpha value */ void alpha_beta(double *y,double *y_smoothed, double alpha,double beta); void profile(double *y); void gauss_noise(double *noise); void plot_it(double *y); void message(char *string1,char *string2, char *string3); void main(void) { int x,graphdriver=DETECT,graphmode,errorcode; double *y_pure,*y_noisy,*y_smoothed,*noise; double alpha,beta; /**********************************/ /* allocate memory for the arrays */ /**********************************/ if((y_pure = (double *)malloc(PLOTEND*sizeof(double)))== NULL) { printf("\a"); cprintf("Unable to allocate space for array "); cprintf("y_pure[]\n\r"); exit(1); } if((y_noisy = (double *)malloc(PLOTEND*sizeof(double)))== NULL) { printf("\a"); cprintf("Unable to allocate space for array "); cprintf("y_noisy[]\n\r"); exit(1); } if((y_smoothed = (double *)malloc(PLOTEND*sizeof(double)))== NULL) { printf("\a"); cprintf("Unable to allocate space for array "); cprintf("y_smoothed[]\n\r"); exit(1); } if((noise = (double *)malloc(PLOTEND*sizeof(double)))==NULL) { printf("\a"); cprintf("Unable to allocate space for array "); cprintf("noise[]\n\r"); exit(1); } /****************************/ /* initialize graphics mode */ /****************************/ registerbgidriver(EGAVGA_driver); initgraph(&graphdriver,&graphmode,""); errorcode=graphresult(); if(errorcode != grOk) { printf("\a"); cprintf("Graphics error: %s\n\r", grapherrormsg(errorcode)); exit(1); } /*************************************************/ /* generate the pure and corrupted data profiles */ /*************************************************/ /* generate the clean profile */ profile(y_pure); /* generate the noise */ gauss_noise(noise); /* corrupt the clean data with noise */ for(x=0;x<PLOTEND;++x) y_noisy[x] = y_pure[x] + noise[x]; /****************************************/ /* plot a pure profile with underdamped */ /* alpha-beta smoothing */ /****************************************/ /* plot the clean profile */ message("Clean data profile","",""); plot_it(y_pure); /* compute and plot the smoothed result over the clean profile */ alpha = ALPHA; /* arbitrary alpha value */ beta = alpha*alpha / (2.0 - alpha); /* standard, underdamped beta value */ alpha_beta(y_pure,y_smoothed,alpha,beta); message("","with", "Underdamped alpha-beta smoothing"); plot_it(y_smoothed); /* clear the screen and replot the smoothed result */ cleardevice(); message("Underdamped alpha-beta smoothing of clean \ profile","",""); plot_it(y_smoothed); /*****************************************/ /* plot a noisy profile with underdamped */ /* alpha-beta smoothing */ /*****************************************/ /* clear the screen and plot the profile */ cleardevice(); message("Noisy profile","",""); plot_it(y_noisy); /* compute and plot the smoothed result over the noisy profile */ alpha_beta(y_noisy,y_smoothed,alpha,beta); message("","with", "Underdamped alpha-beta smoothing"); plot_it(y_smoothed); /* clear the screen and replot the smoothed result */ cleardevice(); message("Underdamped alpha-beta smoothing of noisy \ profile","",""); plot_it(y_smoothed); /* plot the pure profile over the estimate based on the noisy version */ message("","with","Clean data profile"); plot_it(y_pure); /********************************************/ /* plot a pure profile with near-critically */ /* damped alpha-beta smoothing */ /********************************************/ /* clear the screen and plot the clean profile */ cleardevice(); message("Clean data profile","",""); plot_it(y_pure); /* compute and plot the smoothed result over the clean profile */ alpha = ALPHA; /* arbitrary alpha value */ beta = 0.8*(2.0 - alpha*alpha - 2.0*sqrt(1 - alpha*alpha))/(alpha*alpha); alpha_beta(y_pure,y_smoothed,alpha,beta); message("","with","Near-critically damped \ alpha-beta smoothing"); plot_it(y_smoothed); /*********************************************/ /* plot a noisy profile with near-critically */ /* damped alpha-beta smoothing */ /*********************************************/ /* clear the screen and plot the noisy profile */ cleardevice(); message("Noisy profile","",""); plot_it(y_noisy); /* compute and plot the smoothed response over the noisy profile */ beta = 0.8*(2.0 - alpha*alpha - 2.0*sqrt(1 - alpha*alpha))/(alpha*alpha); alpha_beta(y_noisy,y_smoothed,alpha,beta); message("Noisy profile","with", "Near-critically damped alpha-beta \ smoothing"); plot_it(y_smoothed); /* clear the screen and replot the smoothed result */ cleardevice(); message("Near-critically damped alpha-beta \ smoothing of noisy profile","",""); plot_it(y_smoothed); /* plot the pure profile over the estimate based on the noisy version */ message(""," with","Clean data profile"); plot_it(y_pure); /* prepare to exit */ closegraph(); } void alpha_beta(double *y,double *y_smoothed, double alpha,double beta) /* demonstrate the alpha-beta algorithm */ { int i; double s_est,v_est,s_pred,error; /* initialize the filter */ s_est = y[0]; v_est = 0.0; /* loop through all of the position measurements */ for(i=0;i<PLOTEND;++i) { /***********************************/ /* the alpha-beta filter algorithm */ /* (with T absorbed into v) */ /***********************************/ /* predict the new position */ s_pred = s_est + v_est; /* compute the measurement error */ error = y[i] - s_pred; /* estimate the true position and velocity */ s_est = s_pred + alpha*error; v_est = v_est + beta*error; /*****************************/ /* update the smoothed array */ /*****************************/ y_smoothed[i] = s_est; } } void profile(double *y) /* generate the profile to be studied */ { int i; /* first steady-state portion */ for(i=0;i<PLOT1;++i) y[i] = STEADY_STATE1; /* ramp from first steady state level to second steady state level */ for(i=PLOT1;i<PLOT2;++i) y[i] = (i-PLOT1) * (STEADY_STATE2-STEADY_STATE1) / (PLOT2-PLOT1) + STEADY_STATE1; /* second steady-state portion */ for(i-PLOT2;i<PLOT3;++i) y[i] = STEADY_STATE2; /* step to third steady-state portion */ for(i-PLOT3;i<PLOTEND;++i) y[i] = STEADY_STATE3; } void gauss_noise(double *noise) /* generate mean-zero gaussian noise using Box-Muller method */ { int i; double drand1,drand2,pi,sigma=10.0,mean=0.0; /* define pi */ pi = 4.0 * atan(1.0); /* seed the random number generator */ randomize(); /* generate gaussian noise using the Box-Muller method */ for(i=0;i<PLOTEND;++i) { /* compute two double random numbers */ drand1 = (double)rand)()/RAND_MAX; drand2 = (double)rand()/RAND_MAX; /* compute the noise term */ if(drand1 < 1e-10) drand1 = 1e-10; /* avoid division by zero */ noise[i] = sqrt(2.0*log(1.0/drand1)) * cos(2.0 * pi * drand2) * sigma + mean; } } void plot_it(double *y) /* plot one member of array y for each horizontal pixel position */ { int x,y_old; /* plot the profile */ y_old = (int)(y[0]+0.5); for(x=0;x<PLOTEND;++x) { line(x,y_old,x,(int)(y[x]+0.5)); y_old = (int)(y[x]+0.5); } /* wait for key press */ getch(); } void message(char *string1,char *string2, char *string3) /* print a message of three lines centered at bottom of screen */ { int midx,maxy; /* get screen parameters */ midx=getmaxx()/2; maxy=getmaxy(); /* print pieces of message, centered */ settextjustify(CENTER_TEXT,TOP_TEXT); outtextxy(midx,maxy-30,string1); outtextxy(midx,maxy-20,string2); outtextxy(midx,maxy-10,string3); } /* End of File */