Removing drift from noisy accelerometer data

Hi all,

I am using the sensors within my phone to generate a CSV file of accelerations in 3-axis (x,y,z). I have now imported the data to matlab using the CSVread funtion and have began processing the data.

I have applied a filter to reduce some of the noise from the signal however upon integration the signal still drifts. The code I am using is shown below, for simplicitys sake I am only showing data from one axis. Any help is appreciated

clear; close; clc;
D=csvread('test20m3.csv');
t=D(:,1);                                                                           %Define time
XAccRaw=D(:,5);                                                                     %Define X acceleration
XAcc=XAccRaw*9.81;                                                                  %Convert to m/s^2
d=designfilt('lowpassfir','filterorder',10,'CutOffFrequency',10,'SampleRate',100);  %Lowpass FIR filter
AX=filtfilt(d,XAcc);                                                                %Apply filter to data
VX=cumtrapz(t,AX);                                                                  %Integrate acceleration to get velocity
SX=cumtrapz(t,VX);                                                                  %Integrate velocity to get displacement
figure(1);
plot(t,SX);
xlabel(Time (s));
ylabel(Displacement (m))

Answer 1

[Image Analyst]

I think this works great. It's well commented, practically self documenting but if you have any questions, ask.

% Initialization steps.
clc;    % Clear the command window.
fprintf('Beginning to run %s.m ...\n', mfilename);
close all;  % Close all figures (except those of imtool.)
clear;  % Erase all existing variables. Or clearvars if you want.
workspace;  % Make sure the workspace panel is showing.
format long g;
format compact;

% Load data and plot it.
hFig1 = figure;
sa = load('acceleration.mat')
accel = sa.acceleration1;
st = load('time.mat')
t = st.time1;
plot(t, accel, 'b.-', 'MarkerSize', 9);
grid on;
hold on;
fontSize = 20;
xlabel('Time', 'FontSize', fontSize);
ylabel('Acceleration', 'FontSize', fontSize);
title('Original Signal', 'FontSize', fontSize);
hFig1.WindowState = 'maximized'; % Maximize the figure window.
% Draw a line at y=0
yline(0, 'LineWidth', 2);

% A moving trend is influenced by the huge outliers, so get rid of those first.
% Find outliers
outlierIndexes = isoutlier(accel);
plot(t(outlierIndexes), accel(outlierIndexes), 'ro', 'MarkerSize', 15);
% Extract the good data.
tGood = t(~outlierIndexes);
accelGood = accel(~outlierIndexes);
% plot(t(~outlierIndexes), accel(~outlierIndexes), 'mo', 'MarkerSize', 10); % Plot circles around the good data.
% Do a Savitzky-Golay filter (moving quadratic).
windowWidth = 51; % Smaller for tighter following of original data, bigger for smoother curve.
smoothedy = sgolayfilt(accelGood, 2, windowWidth);
hold on;
plot(tGood, smoothedy, 'r-', 'LineWidth', 2);
legend('Original Signal', 'X axis', 'Outliers', 'Smoothed Signal');

% Now it looks pretty reasonable since we didn't include the outliers.
% But smoothedy has fewer points so if we're to subtract it from the original
% we have to fill in the missing points.
smoothedy = interp1(tGood, smoothedy, t);
% Now subtract the smoothed signal to get the variation
signal = accel - smoothedy;
% Plot it.
hFig2 = figure;
plot(t, signal, 'b.-', 'MarkerSize', 9);
grid on;
hold on;
title('Corrected Signal', 'FontSize', fontSize);
xlabel('Time', 'FontSize', fontSize);
ylabel('Acceleration', 'FontSize', fontSize);
hFig2.Units = 'normalized';
hFig2.Position = [.2, .2, .5, .5]; % Size the figure window.
% Draw a line at y=0
yline(0, 'LineWidth', 2);

[Jack Upton]

Hi this is excellent, could you please explain what the smoothedy parts are for? I tried subtractung this from my data and the signal seems to become noisy again. Where is the source for this code, as i would like to reference some parts.

Edit: I think i understand now, I can now apply a filter to smooth this data, as the drift is now removed, correct?

[Image Analyst]

The smoothedy is the red curve, which is the trend. If you subtract that from the original signal, you get the signal without the trend. If you consider that noise, rather than signal, then just go with smoothedy as the signal. You didn't really say what was noise and what was signal that you want to extract from the original signal.

[Jack Upton]

The real data which should be represented would be a straight line at y=0 as the device was stationary when recording data, What I am trying to do is elimate as much noise as possible so when I am integrating to find position, the position is as close to 0 as possible, but im having a lot of trouble doing so

[Image Analyst]

So the whole signal you showed us was 100% noise? But how different is your signal going to be when you have actual, real accelerations? Is it as big as the noise? A hundred times bigger? If the real signal is way bigger than the noise then a simple smoothing would probably be fine. If your signal is about as big as the noise then how can you know for a given amplitude if it's all noise, all signal, or some of each?

[Vijay HP]

Thanks a lot @Image Analyst, the code really helped

[Emily Keys]

Hi, I used this code for some acceleration data recodring a weight going across a board to get a quasi-static movement i have. It worked great except for when i integrated it twice to get displement and the post load had a massive amount of drift. Any ideas how to help fix that? thanks

[Image Analyst]

@Emily Keys not really sure without seeing your data.

If you have any more questions, then attach your data and code to read it in with the paperclip icon, in a new discussion thread, after you read this:

TUTORIAL: How to ask a question (on Answers) and get a fast answer

[Emily Keys]

Apologies for being too vague. This is the data as acceleration and then as the displacement in the second figure, with the orange line being the actual displacement measured and the blue being the method you provided applied to the acceleration and then integrated twice. As you can see in the figure the data drifts far off to negative 100 when it should go back to 0. Any idea why this happens would be greatly appreciated.

[Image Analyst]

Same reply (not sure why you ignored it). Start a new discussion thread and attach your data and code, especially the code to turn acceleration into displacement.

[Emily Keys]

Hi , I have done so except the data that I used is too large to upload to even in zip file format. Apologies again. https://uk.mathworks.com/matlabcentral/answers/1943334-fixing-post-load-for-acceleration-signal-with-savitzky-golay-filter

Answer 2

[Ameer Hamza]

If by drift, you mean the signal continually move away from the actual value then you can try to use detrend() function: https://www.mathworks.com/help/matlab/ref/detrend.html

[Jack Upton]

unfortunately this still doesn't portray what is actually going on and generates unexpcted data.

[Ameer Hamza]

Can you attach sample data in a mat file?

[Jack Upton]

Ive attatched the first 200 data points in a .mat files

[Ameer Hamza]

I am not sure if this will work for a general dataset, but for this dataset, this seems to reduce the noise. I used simple filtering functions, you can get better performance by carefully designing the filters.

plot(time1, acceleration1, 'DisplayName', 'raw');
hold on

acc1 = detrend(acceleration1);
acc1 = movmean(acc1,10);
plot(time1, acc1, 'LineWidth', 2, 'DisplayName', 'mov-average');

sys = tf(1, [1 2]);
acc2 = detrend(acceleration1);
acc2 = lsim(sys,acc2,time1,acc2(1));
plot(time1, acc2, 'LineWidth', 2, 'DisplayName', 'low-pass');

legend

[Jack Upton]

I’ll check this through properly when I arrive back home shortly. What results show if you integrate the filtered data you have generated once to get velocity, and again to get displacement

[Ameer Hamza]

Following images shows the velocity and displacement profiles using the raw signal and filtered signals.

Velocity:

Displacement:

[Jack Upton]

Thank you this is amazing! Im now going to apply it to the whole of my data and compare with other solutions i have found.

[Ameer Hamza]

The advantage of using the filters I suggested in the answer is that

1. If you use a movemean filter, it is quite easy to implement, so even if you want to use it with the accelerometer on an embedded system, it can be quickly implemented.
2. If you use the continuous-time filter, you can discretize it using c2d() command in MATLAB, and directly implement on an embedded system.

[Jack Upton]

How would i implement the c2d function? Im having trouble with my data still. For the whole length of the data it shows 200m of range when in reality the accelerometers were stationary

[Ameer Hamza]

If there the bias still persists in your sensor, then simple filtering techniques will not work. You might have to use the Kalman filter, which takes into consideration the sensor values and dynamical model of your system to predict the next state. It is a bit complicated to understand, but it is the way to go if there are biases in your sensor data, and you also want to measure displacement using the acceleration sensor data.

[Jack Upton]

I thought that would be the case as thats what I have read elsewhere,Ive tried to implement it previously using GPS as well as my accelerometer results but I have had no luck

[Ameer Hamza]

Yes, it can be a bit complicated because you need to estimate the model of your system. You can try to follow some example on Kalman filter in MATLAB: https://www.mathworks.com/discovery/kalman-filter.html

https://www.mathworks.com/help/control/examples/state-estimation-using-time-varying-kalman-filter.html