I have attached an Excel sheet that I made with the equations needed to design a barometric altitude sensor.
For lower altitudes and for small altitude deltas, one can assume that there is a linear relationship between pressure and altitude. Please see the attached Excel sheet for more details.
The SCP1000 pressure sensor seems like a very good candidate for a pressure sensor:
http://www.sparkfun.com/commerce/product_info.php?products_id=8161
The manufacturer gives these data:
Resolution: 3 Pascal (the fineness to which the sensor can be read - nothing to do do with the accuracy)
Relative accuracy (how close to the real pressure the sensor measurements are) over temp range +10C to +40C and 600 hPa to 1200 hPa: +/- 50 Pa. (It is not given in the datasheet if this is sensor to sensor or for the same sensor)
The manufacturer does not give noise data. However sources on the Internet may indicate that the sensor has an inherent noise of +/- 50 Pa.
Datasheet: http://www.sparkfun.com/datasheets/Components/SCP1000-D01.pdf
Sample code: http://www.sparkfun.com/datasheets/Components/SCP1000-D01.pdf
Another guy using the sensor for a variometer for paragliding: http://www.pixelproc.net/varios.html
The noise may be filtered out using an averaging filter however. Even a kalman filter and two sensors may be used.
baro_height_measurement_unit.xls
Comments
1-exp(-Dt/Tgps) has nothing to do with the altitude.
The first order filter has the continuous form: Tgps.d(h_filt)/dt = - h_filt + h
where h is the input of the filter and h_filt is the output
If you want to transform this continuous form into the discrete form you would want to write :
Tgps.(h_filt(n) - h_filt(n-1))/Dt = - h_filt(n-1) + h(n)
h_filt(n) = h_filt(n-1) . (1-Dt/Tgps) + h(n).Dt/Tgps
h_filt(n) = h_filt(n-1) + Dt/Tgps . (h(n) - h_filt(n-1))
Thus the expression :
h_filt += Dt/Tgps . (h(n) - h_filt(n-1))
But this expression is OK if Dt is very small compared to Tgsp.
If not you have to integrate the above expression over Dt with the input constant and you get :
h_filt += (1-exp(-Dt/Tgps)) . (h(n) - h_filt(n-1))
Note : 1-exp(x) ~ x if x is very small
In our case it is true that 25ms << 2s (EM406 time delay) so I could have directly written Dt/Tgps
Paul
I don't have my PC and I have a lot of difficulties trying to write a comment on DIY Drones with my smart phone. I will answer you tomorrow.
The gps alt is delayed so we must delay the predicted alt with the same time lag before we compare them. The exponential term is required just for the computation of the time filter. I will explain that tomorrow.
Regards,
Paul
In case you would like to take a look at it and explain. Perhaps I use wrong Dt and / or Tgps values?
UFO_MAN
I assume the idea is that the when the data from the pressure sensor is to be compared to the data from the GPS, the data from the GPS is already t time periods old (it was sampled t time periods ago).
I assume the 1-exp(-Dt/Tgps) has something to do with that the value of the height from the pressure sensor has increased or decreased as a continuous compounding (therefore the natural logaritm base) since the GPS data was taken and that the GPS data therefore needs to be transformed forward in time to match the pressure sensor data. However, when I try to use Dt 25ms and Tgps 1000ms , I get a GPSFILT value of 0,0487. This is a very small value. I cant imagine that the GPS data has to be corrected that much.Can you explain in more detail? Have I got it wrong?
I also assume the correction happens every 25ms, but that the data is actually used every GPS sample, so it is transffered from the 25ms loop to the 1sec loop.
UFO_MAN
Here is the plot of GPSFILT per Dt/Tgps
This filter has some sense if and only if the mean GPS altitude is reliable. If this mean has some kind of offset that is greater than the simple baro altimeter error then this filter can be discarded.
I am not a specialist, but I suppose that the error of the simple baro altimeter is due to the pressure sensor gain (effect of temperature and aging), to the atmosphere temperature (we could compensate that with a measurement of the temperature)
So the main benefit of the filter (if the mean GPS altitude is reliable) is to auto calibrate the baro altimeter.
GPSFILT = 1 - exp(-Dt/Tgps)
where Dt is the pressure sensor sampling rate (= 25ms for example),
and Tgps is the GPS time response (to be identified).
This first order filter is required because the GPS altitude has some delay. So the Kalman filter output, the estimated altitude, must be delayed by the same amount of time before computing the correction.
Of course the altitude is zeroed during initialization. This way h will represent the AGL.
The 1st order filter output (state) h_filt is also zeroed.
I have worked out all the matrix computations in order to minimize the computational load.
Paul
Can you specify the plots you want?
The code is programmed in C. Most of the Matlab S-function can be throwned away, it is the interface with the Simulink integration engine.
Here is the core:
/* Prediction phase (each 25 ms) */
/* ====================== */
{
/* Input = calibarated pressure sensor */
/* ---------------------------------------------- */
/* Pressure variation P(n) - P(n-1)*/
dP = pressure - pressure_old;
/* Predicted altitude */
/* ------------------ */
/* Note: sensitivity coefficient state[1] stays unchanged */
state[0] += state[1]*dP;
/* Predicted covariance matrix */
/* --------------------------- */
/* Save Vhc */
f_buff = variance[2];
/* Vhc = Vhc + Vc*dP */
variance[2] += variance[1]*dP;
/* Vc = Vc + Qc */
variance[1] += Q[1];
/* Vh = Vh + 2.Vhc.dP + Vc.(dP)^2 + Qh */
variance[0] += (f_buff + variance[2])*dP + Q[0];
/* Phase (align) estimated altitude with GPS altitude (GPS has some delay) */
h_filt += GPSFILT * (state[0] - h_filt);
/* Save pressure */
pressure_old = pressure;
}
/* Correction phase (at GPS sampling rate) */
/* ========================= */
{
/* Residual covariance: RC = Vhp + R */
/* --------------------------------- */
RC = variance[0] + GPS_NOISE;
/* Kalman gain vector */
/* ------------------ */
/* Kalman altitude gain: Kh = Vhp/RC */
Kh = variance[0]/RC;
/* Kalman coefficient gain: Kc = Vhc/RC */
Kc = variance[2]/RC;
/* Altitude error */
/* ----------------- */
dh = hgps - h_filt;
/* Corrected state */
/* --------------- */
state[0] += Kh*dh; // h = hp + Kh*dh
state[1] += Kc*dh; // c = cp + Kc*dh (note: cp = c)
/* Corrected state covariance matrix */
/* --------------------------------- */
variance[0] -= Kh*variance[0];
variance[1] -= Kc*variance[2];
variance[2] -= Kh*variance[2];
}
-
1
-
2
-
3
-
4
-
5
of 6 Next