2002 NARAM R&D Project

My 1st place winning C division project for the 2002 NARAM was "Barometric Apogee Detection Using the Kalman Filter". Kind of a high highfalutin name but the idea is very simple - improving apogee detection in altimeters used in model and high power rockets.

The full text of the report is in a Adobe Acrobat file here.

If you have an RDAS and would like to try the filter out on your data, download this program

This is a command line program ("DOS") that runs in a Windows "DOS" window. It is a UNIX style filter program that reads from standard input and writes to standard output. To use it type:

rdaskal3 < yourdatafile > outputfile

This will input "yourdatafile" which you must export from the RDAS program. Be sure to use the interpreted output if you want the correct units. The filtered output will be in "outputfile". To look at the data, use GNUPlot which is included with the RDASPlot program.

After starting up GNUPlot, change directories to where your data file is. If you are not sure how to do this just use the built in help to learn how.

To plot the filtered altitude data use:

plot "outputfile" using 1:3 with lines

To plot estimated velocity use:

plot "outputfile" using 1:4 with lines

To plot the estimated acceleration use:

plot "outputfile" using 1:5 with lines

To zoom in on specific times use:

plot [15:25] "outputfile" using 1:3 with lines

Which will plot the data from 15 to 25 seconds.

This will give one of my favorite plots when zoomed in around apogee.

set grid
set y2tics
set xlabel "Time from First Motion (seconds)"
set ylabel "Altitude (feet)"
set y2label "Velocity (feet/second)"
plot [15:25] "outputfile" using 1:3 title "Filtered Altitude" with lines
replot "outputfile" using 1:2 title "Unfiltered Altitude" with lines
replot "outputfile" using 1:4 axes x1y2 title "Velocity" with lines

The result is filtered and unfiltered altitude data along with the velocity estimate all in one picture. It should be really easy to see where apogee is on this one. That should be more than enough to get you started.

Here is a try at coding the filter for the PIC16 processors. While this code has been through the MPASM mill, it has not been run and likely contains numerous horrible awful bugs. For educational purposes only. Use at your own risk.

The week before NARAM when I should have been building competition rockets, I was still messing around with the Kalman filter. I updated the filter so that it would use both the pressure altitude and accelerometer measurements in estimating altitude, velocity, and acceleration. The improvement was amazing. I was able to use a ratio of acceleration measurement noise to model noise of 100 to 1 or less. I usually run the pressure only filter at 50,000 to 1. The result is really good acceleration and velocity estimates.

Because of the low tuning ratio, the filter still gets upset during Mach transitions but the velocity estimate doesn't decrease below zero on the data sets that I have. This still concerned me so I added some extra code.

The data sets I have for Mach 1+ flights exhibit a step increase in pressure altitude when going above Mach 1 and a step decrease when dropping below. The increase doesn't worry me as that tends to increase the estimated velocity. The decrease is a major worry since that decreases the velocity estimate. So I added some code to test for:

A large difference between measured and estimated altitude.

Velocity estimate between 1200 and 800 fps.

Negative acceleration.

When these conditions are present, the correction based on the difference between estimated and measured altitude is set to zero. This prevents the velocity and acceleration estimates from being corrupted. Then the estimated altitude is set to the measured.

The result is that during this period of time the algorithm is basically that of an integrating accelerometer. The accuracy isn't as good but as soon as the speed drops below 800 fps, the pressure measurement kicks in again and all is well.

Here is an example using a data set from a hybrid flight. (Not mine!) Note that the altitude exceeded the range of the pressure sensor. Which reminds me. It might seem obvious but perhaps it isn't. If the measurement range of the sensors is exceeded, the filter has problems. The problem with the pressure measurement is handled by checking to see if it is above its valid range and then ignore it. Effectively the filter reverts to integrating the accelerometer. It works pretty well since time of apogee as determined by this filter agrees to within one second with the GPS data from the flight.

The period between 25 and 30 seconds is where the Mach transition occurred and you can clearly see the altitude tracking the measured altitude through the transition. But there is no impact on the velocity estimate. After looking at the first plot I made, I went back and removed the test for negative acceleration. It produces a much better velocity estimate during motor burn this way.

Oh. Here is a plot of the unfiltered data.

Another measurement error that is possible is exceeding the range of the accelerometer. While this is fairly difficult to do with a +/-50G range it is possible. A local flier flew their Shadow Composites Raven on an I435 a while back. He forgot to switch the AltAcc from the +/-25G range and the acceleration data was clipped badly. It also caused a very early deployment since the integrated velocity was much too low due to the clipping. Deployment occurred while the Raven was still traveling at about 250fps. Amazingly, it survived to fly again.

I was curious as to how the filter would handle this situation. So I altered the input code of the RDAS filter so that it could read data from the AltAcc. At the time the AltAcc deployed, the filter had recovered nicely from the acceleration over-range and was tracking altitude nicely and showing a velocity of 250fps. So it looks like it will recover from this situation. Unless you can exceed the range of the accelerometer until the pressure sensor exceeds its range in which case all bets are off. :-)

I have coded up the two measurement filter for the PIC processor and it was really quite easy. Total code space for the filter is 330 words which fits easily into just about any PIC. Execution time is well within the capability of all PICs. I hope to have this up and running on real hardware before very long.

I have changed the sampling rate on some of my RDAS data sets to as low as 2 samples/second. The two measurement filter still performed great. The only limit on accuracy of time of apogee was the low sample rate and even that could be worked around. Since velocity and acceleration are known, it is pretty easy to predict the time of apogee. Just do a quick check (velocity + acceleration/sample rate < 0) and then compute the actual time of apogee and set a timer interrupt.

I have moved the flight hardware/software stuff to its own page.

27 Oct. 2003

One of the big questions that I still have about the Kalman filter is what to do about the model noise covariance matrix. I wasn't really sure to do so I only used one term for the acceleration variance and left all the other 8 terms at zero. This has worked pretty well so far.

But then I noticed a pdf file on my hard drive with the tempting file name of "kalman.pdf". This turned out to be an article from Embedded Systems Programming magazine (June 2001, page 72) and it actually dealt with the covariance matrix issue and it even had an example! So I revisited my code and added the 8 missing terms.

I compared the results and the changes were underwhelming. For the most part, there was no change in the results. The Kalman gains changed slightly but had little effect on the performance. Then I ran it on the N hybrid data file. Now I could see a change. This data set is unique among the data sets I have in that the pressure sensor range is exceeded. I altered the code to check for this and the filter ignores the pressure measurement when this happens. The results were pretty good as the apogee time (which is what I am after) agreed with the recorded GPS data within one second.

But the modified code detects apogee (zero velocity) about 1.5 seconds earlier than the previous version. And this is in the direction that increases the difference with the GPS time. But, I could (and I will :-) argue that the GPS time is actually late.

If the GPS engine were simply reporting the fix it computes from the SV time deltas, there would be no problem. But I have seen some indications from my Garmin GPS units that they are also including a dynamic motion model to help smooth the positions.

I had noticed a couple of times that the GPS position would keep moving for a couple of seconds after I stopped moving. I didn't pay too much attention to that until now. I also noticed an interesting anomaly in the data I recorded at LDRS. If you look at the lower left corner of the graph, you will see that there is a data point along the road south of the turn to g back to the flying field. I did NOT miss that turn and then back up. So something else was going one that fooled the GPS receiver into thinking that I continued south instead of making a left turn. That could only be a dynamic model which could be as simple as a constant velocity model.

Such a model would result in the GPS position being slightly in error if the motion did not match the dynamic model precisely. Something like continuous acceleration when you assume constant velocity.

I have no way of knowing what is going on inside the GPS unit used with the RDAS but I suspect that it is not simply reporting a position fix. Instead I suspect that it is smoothing the data with a filter (Kalman I bet!) using a constant velocity model.

I have been playing around with the Kalman filter program for a while now and I finally realized that I have several versions for various altimeters. Some of which I don't even own! So I think I need to present all of them in some vaguely organized fashion. So here it is:

The "DOS" versions were compiled under Win98 using the djgpp compiler. They may or may not work on other versions of Windows.

In order to filter the code the data must be in a plain text file. Each program reads from standard input and writes to standard output so you will have to redirect appropriately.

Home