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:
To your velocity looking like:
Appears like some improvement...
Also, another thing is that if you have gyroscope values outputted AND the motion that you want to visualize is a swinging motion at constant radius, then the gyroscope values themselves alone can give you an excellent representation of velocity, and you can ignore the accelerometer values all together.
See below a comparison of using accelerometer vs gyroscope values to calculate velocity and distance travelled:
Versus for the gyroscope:
Much less little drift, very little noise. It tracks the x-axis right up until the movement, tracks the target motion and shows the small deviations from constant angular velocity, and snaps right back to the x-axis after the movement is complete. This is without any filters or corrections etc.
Wow, what a guy...roscope.
