How do I filter my accelerometer data so that when I integrate it I have no drift?

Hi all,

I'm attempting to integrate my accelerometer data x2 to get displacement. I need to design a filter in order to avoid drift and I've done an fft of my data to view the signal in the frequency domain and don't really know where to go from there. I implemented a high pass filter but ut it dropped my data to a mean of around zero when I should be at least be having around 1 G in my data.

I see a have an asymmetrical signal at 0 Hz but don't know how to interpret it.

My raw acceleration data is above and when I integrate my raw data I have discrepencies.

How should I design my filter?

Answer 1

[Star Strider]

The fft plot does not make sense.

Before you do any signal processing (including the fft), first be certain that the signal is sampled with consistent sampling intervals. One way to do that is to calculate:

sigstd = std(diff(t))

where ‘t’ is the time vector. It should be close to zero. If not, then use the resample function to resample it to a consistent sampling frequency.

I would re-calculate the fft, and then use a bandpass filter (such as this one) to filter the signal. That will eliminiate the gravity offset (that would otherwise be integrated as well). Determine the passband from the fft result.

Fs = 125;                                           % Supply Correct Sampling Frequency
Fn = Fs/2;
Wp = [0.03 0.2]/Fn;                                 % Supply Correct Passband
Ws = Wp.*[0.95 1.05];                               % Calculated Stopband
Rs = 50;
Rp =  1;
[n,Wn] = ellipord(Wp,Ws,Rp,Rs);
[z,p,k] = ellip(n,Rp,Rs,Wp);                        % Design Elliptic Filter
[sos,g] = zp2sos(z,p,k);

figure
freqz(sos, 2^16, Fs)
% set(subplot(2,1,1), 'XLim',[0 0.5])                 % Zoom Axis (Optional) Choose The Correct Limnits
% set(subplot(2,1,2), 'XLim',[0 0.5])                 % Zoom Axis (Optional) Choose The Correct Limnits

signal_filt = filtfilt(sos, g, signal);

It will be necessary to experiment with this to get the result you want.

.

[Ian Larson]

Thanks for the response. My accelerometer was used to calculate the deflection of a beam, so i'm not sure If I would need to eliminate the gravity offset. my ans from calculating sigstd was very close to zero, in the order of e-15. My code now is a little messy from forum searching and when I use a filter or detrend my acceleration plot oscilates from about 0 m/s^2 instead of around 10, which is desired.

% Acceleration data 
% acceleration is in m/s^2        


tStep = 1.620*10.0^-3;  % Length of each time step
N = length(acc)*tStep; %complete time 
t = 0:tStep:N;
t(end) = [];
N = length(t);
dt = mean(diff(t));       % Average dt
fs = 1/dt;                  % Frequency [Hz] or sampling rate
% some additionnal high pass filtering
N = 4;
fc = .1; % Hz
[B,A] = butter(N,2*fc/fs,'high');  
%acc=filtfilt(B,A,acc);
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
acc2 = filter(B,A,acc); 
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
velocity = cumtrapz(dt,acc2);
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
disp = cumtrapz(dt,velocity);    
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
figure(1); 
accnew=detrend(acc)
plot(t,acc)
title('Acceleration VS Time') 
xlabel('time (s)') 
ylabel('m/s^2') 
%%
transformacc=fft(acc); 
figure(2) 
f=linspace(0,fs,length(acc));
plot(f,abs(transformacc)) 
xlabel('Frequency (Hz)')  
xlim([0, 300000])
ylabel('magnitude')
%% 
filteredsig=filter(Hd,acc); 
figure(3); 
plot(t,filteredsig)
figure(4); 
plot(t,velocity);
title('Velocity VS time') 
xlabel('time (s)') 
ylabel('Velocity (m/s)')
figure(5);
plot(t,disp); 
title('Displacement VS Time') 
xlabel('time (s)') 
ylabel('Displacement (m)') 
grid on;
%%
newfft=fft(filteredsig); 
figure(2) 
fnew=linspace(0,fs,length(filteredsig));
plot(fnew,abs(filteredsig)) 
xlabel('Frequency (Hz)')  
xlim([0, 300000]);
ylabel('magnitude') 

[Star Strider]

My pleasure!

It would help to have the original data.

The highpass filter will eliminate the gravity offset.

Creating a highpass filter using my ellpitic filter would be preferable to a Butterworth design. Use filtfilt, not filter, since filtfilt is phase-neutral (no phase distortion or delay).

The second-order-section realisation creates a stable filter. The transfer function realisation does not.

There appears to be high-frequency noise in the signal (the ‘constant’ sections at the beginning of the signal appear to display high-frequency noise). so a bandpass filter might be more appropriate.

I do not see any specific trend in the signal, so detrending it may not be appropriate.

Answer 2

[Walter Roberson]

I need to design a filter in order to avoid drift

That is not going to work. Errors in accelerometers are not symmetric. Any system of position measurement that depends only on accelerometers will either exhibit drift or else will be inaccurate for actual measurements.

Suppose your accelerometer can measure down to 1/100 g. Suppose you have measurement noise that is mostly on the order of 1/1000 g. Suppose 1 time in 10000, the measurement noise adds up to +1/100 g . If the measurement system was continuous, then two -1/200g noises would balance out the +1/100g. But -1/200g is below the measurement limit so anything of less magnitude than -1/100g is not going to counter-balance the +1/100g noise.

This is the sort of reason that missile navigation systems for longer distances never use dead reconning based just on acceleration measurements.

[Ian Larson]

Thanks, I appreciate the answer. That was insightful.

[Walter Roberson]

https://en.wikipedia.org/wiki/Inertial_navigation_system#Drift_rate

"Even the best accelerometers, with a standard error of 10 micro-g, would accumulate a 50-meter error within 17 minutes"