# Removing drift from integration of accelerometer data

I am trying to get a positional data from the accelerometer data using the following steps:

1. Re-zero the accelerometer value
2. Removing mean from accelerometer value
3. First integration of accelerometer value to get velocity
4. Removing mean from velocity
5. Integrate velocity to get position

, and the following is the result from every integration value:

Can anyone suggest on how to remove the drift appear on the velocity value? Some suggest to push the velocity value to zero when accelerometer value is near to zero but I don't think this is very robust.

• "Some suggest to push the velocity value to zero when accelerometer value is near to zero but I don't think this is very robust." That would be wrong. If you're cruising at constant speed in a car, acceleration is zero but velocity is not. – endolith Sep 27 '16 at 21:29

Double integration amplifies any offsets, non-linearities and noise. These can't be removed without the use of some type of external reference point measurements (e.g. not from the accelerometer) or other information. Calibration might be able to remove most of the ones that do not vary with time (temperature, etc.)

For instance, if the size of the room is known, you can robustly clamp position to be inside the walls and above the floor.

• I also have gyroscope value. would that in any way, help? – Sharah Sep 27 '16 at 21:42
• @Sharah : Doubtful. You need an absolute position measurement (or relative to the room or Earth, etc.). Maybe from the camera (photometrics with known size markers). Or differential GPS. Or human input if accurate. etc. – hotpaw2 Sep 27 '16 at 23:40

I can say a couple of things about this problem since I have been working on a similar one for a few weeks. First thing to note is that you should subtract the mean of a section of noise from your data, not the entire section.

As below, moving from figure 1 to figure 2:

You could also fit a 0 degree polynominal and subtract from the original data:

A = num - np.poly1d(np.polyfit(time[k1:k2], num[k1:k2], degree))(time)


Remember you are finding the mean/0 degree constant value from a section of noise, but applying it over the whole range of data.

Also, an important thing to note is that the gyroscope values you mentioned DO help you, if you want a more complicated method of removing gravity, gyroscope allows you to track the direction which is receiving gravity.

My accelerometer fortunately outputted "quaternion" units, which are a hybrid of the accelerometer and gyroscope values, and allows a dynamic and instantaneous calculation of the acceleration due to gravity fraction that each x, y, z direction is receiving at all times.

This is calculated according to the equations below:

Xgrav = 2*((Qw*Qx) + (Qy*Qz))
Ygrav = 2*((Qx*Qz) - (Qw*Qy))
Zgrav = ((Qw*Qw)-(Qx*Qx)-(Qy*Qy)+(Qz*Qz))


Note that the order of those equations applying to x, y, z respectively will depend on the setup orientation of your accelerometer.

Again, apologies if quaternion values are irrelevant to your datasets. The accelerometer I used is the GCDC HAM IMU which has Ax, Ay, Az, Gx, Gy, Gz, and Qw, Qx, Qy, Qz.

Then... there will still inevitably be noise and integration errors in the accelerometer data, as well as, what is basically referred to as the accelerometer sensor "losing it's reference frame" and not really being able to tell which direction in the data to add the magnitude of the forces that it's receiving onto.

There are a lot of options available regarding filters and regression corrections. I have looked at moving-average smoothing, bandpass filters, wiener filters etc. Overall, those are going to distort the actual magnitude of the peaks of your acceleration data, and give you incorrect max/min value. Also, the accelerometer signal is not periodic in nature, and, a bandpass filter will turn it into that.

So, I ended up going more for linear regression corrections, more than algorithmic filters. For example, if you can tell the point where your accelerometer has basically "lost it's reference frame" by looking at the intercept on the x-data where the drifting noise points to, you can fit a linear line to the drift, split the data, subtract that line from the original data up to the point where you think the drift started, and join the data and carry on...

As below:

Basically, this:

Becomes this:

In python using scipy:

l1, l2 = = 140000,141000

R = linregress(time[l1:l2], Vz[l1:l2])

inter = -R[1]/R[0]
start = np.where(dfVwave['time'] > inter)[0][0]

transformed = Vz[start:] - (R[0]*time[start:] + R[1])
Vznew = np.concatenate((Vz[:start], transformed))


Code is essentially:

• select slice of drifting linear noise only

• fit linear line to this

• solve for where this intersects with the x-axis and assume that is where the drift started

• subtract this line from your original data up unto that point

• combine this data with preceding original data

Now, the actual accuracy of doing this to your data is probably questionable. So, you basically need to already have in mind what the correct answer/realistic visualizations for your data is.

Another related noise correction strategy, different to using averaging filters or bandpass or kalman, is wavelets. The theory is very complex but pywavelets makes it very easy to implement. Best of all is it does not easily distort your data. You can reduce right down to the barebones of your signal and it will still keep a very similar max value.

The basic wavelet commands are wave deconstruct and wave reconstruct, and then you choose the number of levels wavelets that you want to omit to remove the right amount of noise while preserving the signal you want.

def wavelets(data, wavelet, uselevels, mode):

levels = (np.floor(np.log2(Axcorrected.shape[0]))).astype(int)

omit = levels - uselevels

coeffs = pywt.wavedec(data, wavelet, level=levels)

A = pywt.waverec(coeffs[:-omit] + [None] * omit, wavelet, mode=mode)

return A

Axwave = wavelets(data = Axcorrected, wavelet = 'haar', uselevels = 7, mode = 'zero')


Note: the levels command here yielded 12 for my dataset

Also wavelet can implement an absolute value threshold of which to remove values under.

def wave3(data, mult, mode):

coeffs = pywt.wavedec(data, 'haar', level=12)

threshold = np.std(data[k1:k2])*np.sqrt(2*np.log2(data.size))*mult

new_coeffs = map(lambda data: pywt.threshold(data, threshold, mode=mode, substitute = 0), coeffs)

return pywt.waverec(list(new_coeffs), 'haar')

Vxwave3 = wave3(data = Axcorrected, mult = 10, mode='hard')


So overall you can end up going from having your data look like: