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

Question

Ian Larson 2022-4-10

0
链接

此问题的直接链接

https://ww2.mathworks.cn/matlabcentral/answers/1693210-how-do-i-filter-my-accelerometer-data-so-that-when-i-integrate-it-i-have-no-drift

评论： Walter Roberson 2022-4-10

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?

0 个评论
显示 -2更早的评论隐藏 -2更早的评论

请先登录，再进行评论。

请先登录，再回答此问题。

Answer 1

Star Strider 2022-4-10

1
链接

此回答的直接链接

https://ww2.mathworks.cn/matlabcentral/answers/1693210-how-do-i-filter-my-accelerometer-data-so-that-when-i-integrate-it-i-have-no-drift#answer_939535

在 MATLAB Online 中打开

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.

.

2 个评论
显示无隐藏无

Ian Larson 2022-4-10

在 MATLAB Online 中打开

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 2022-4-10

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 2022-4-10

1
链接

此回答的直接链接

https://ww2.mathworks.cn/matlabcentral/answers/1693210-how-do-i-filter-my-accelerometer-data-so-that-when-i-integrate-it-i-have-no-drift#answer_939605

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.