I did a 1 minute hover flight today with the down-facing LiDAR + video reinstated to track the performance of the dynamic gravity values from the Butterworth as temperature shifts during the flight:

60s hover stats

To be honest, it’s really hard for me to work out if the Butterworth is working well extracting gravity from acceleration as the IMU temperature falls. There are so many different interdependent sensors here, it’s hard to see the wood for the trees!

As an example, here’s my pythonised Butterworth class and usage:

####################################################################################################
#
# Butterwork IIR Filter calculator and actor - this is carried out in the earth frame as we are track
# gravity drift over time from 0, 0, 1 (the primer values for egx, egy and egz)
#
# Code is derived from http://www.exstrom.com/journal/sigproc/bwlpf.c
#
####################################################################################################
class BUTTERWORTH:
def __init__(self, sampling, cutoff, order, primer):
self.n = int(round(order / 2))
self.A = []
self.d1 = []
self.d2 = []
self.w0 = []
self.w1 = []
self.w2 = []
a = math.tan(math.pi * cutoff / sampling)
a2 = math.pow(a, 2.0)
for ii in range(0, self.n):
r = math.sin(math.pi * (2.0 * ii + 1.0) / (4.0 * self.n))
s = a2 + 2.0 * a * r + 1.0
self.A.append(a2 / s)
self.d1.append(2.0 * (1 - a2) / s)
self.d2.append(-(a2 - 2.0 * a * r + 1.0) / s)
self.w0.append(primer / (self.A[ii] * 4))
self.w1.append(primer / (self.A[ii] * 4))
self.w2.append(primer / (self.A[ii] * 4))
def filter(self, input):
for ii in range(0, self.n):
self.w0[ii] = self.d1[ii] * self.w1[ii] + self.d2[ii] * self.w2[ii] + input
output = self.A[ii] * (self.w0[ii] + 2.0 * self.w1[ii] + self.w2[ii])
self.w2[ii] = self.w1[ii]
self.w1[ii] = self.w0[ii]
return output
#-------------------------------------------------------------------------------------------
# Setup and prime the butterworth - 0.01Hz 8th order, primed with the stable measured above.
#-------------------------------------------------------------------------------------------
bfx = BUTTERWORTH(motion_rate, 0.01, 8, egx)
bfy = BUTTERWORTH(motion_rate, 0.01, 8, egy)
bfz = BUTTERWORTH(motion_rate, 0.01, 8, egz)
#---------------------------------------------------------------------------------------
# Low pass butterworth filter to account for long term drift to the IMU due to temperature
# change - this happens significantly in a cold environment.
#---------------------------------------------------------------------------------------
eax, eay, eaz = RotateVector(qax, qay, qaz, -pa, -ra, -ya)
egx = bfx.filter(eax)
egy = bfy.filter(eay)
egz = bfz.filter(eaz)
qgx, qgy, qgz = RotateVector(egx, egy, egz, pa, ra, ya)

Dynamic gravity is produced by rotating the quad frame accelerometer readings (qa(x|y|z)) back to earth frame values (ea(x|y|z)), passing it through the butterworth filter (eg(x|y|z)), and then rotating this back to the quad frame (qd(x|y|z)). This is then used to find velocity and distance by integrating (accelerometer – gravity) against time.

Sounds great, doesn’t it?

Trouble is, the angles used for the rotation above should be calculated from the quad frame gravity values qg(x|y|z). See the chicken / egg problem?

The code gets around this because angles for a long time have been set up initially on takeoff, and then updated using the gyro rotation tweaked as an approximation to the rotation angle increments. During flight, the qa(x|y|z) angles are fed in over a complimentary filter.

Thursday is forecast to be cold, sunny, and wind-free; I’ll be testing the above with a long GPS waypoint flights which so far lost stability at about the 20s point. Fingers crossed I’m right that the drift of the accelerometer, and hence increasing errors on distance and velocity resolves this. We shall see.

P.S. I’ve updated the code on GitHub as the Butterworth code is not having a negative effect, and may be commented out easily if not wanted.