I faced the same problem.

Drift in your velocity and displacement or position comes from the numerical integration, which is an approximation, even when you use better techniques like Simpson's 1-4-1 integration rule instead of trapezoidal rule. 

Noise are high frequency signals where as drifts are low or very low frequency signals.

If drift is not desirable use high pass filter for your result (velocity and displacement) with appropriate cutoff frequencies to avoid it

If noise is also present in your result along with drift, use band-pass filter with appropriate pass and stop frequencies.

Instead of using band-pass filter, using low pass filter for measurement (acceleration) to remove noise and then high pass filter for integrated signals (velocity and position) yields better integrated signals.

You have to be aware of Nyquists law of sampling. Lets say the fastest change of values in your system has the frequency f. The sampling frequency needs to be at least double of that value(2f) to get valid data. This is the theoretical minimum, if your system can't meet this the data won't be valid.

The other thing that could be hampering your attempts is data noise. You can reduce the need for a very high sampling frequency by filtering out higher frequency noise. In that way you will be compliant with Nyquist, and won't need to get better data.

I faced the same problem.

Drift in your velocity and displacement or position comes from the numerical integration, which is an approximation, even when you use better techniques like Simpson's 1-4-1 integration rule instead of trapezoidal rule. 

Noise are high frequency signals where as drifts are low or very low frequency signals.

If drift is not desirable use high pass filter for your result (velocity and displacement) with appropriate cutoff frequencies to avoid it

If noise is also present in your result along with drift, use band-pass filter with appropriate pass and stop frequencies.

Instead of using band-pass filter, using low pass filter for measurement (acceleration) to remove noise and then high pass filter for integrated signals (velocity and position) yields better integrated signals.

If the drift problem is caused by some bias in your acceleration measurement, you could also try one of the following:

- Identify the bias by keeping your sensor non-accelerated and record the measurement signal. This gives you additional information on the sensor noise as well. If you know the bias, substract it from your acceleration data before integrating it.

- assuming constant bias over time you may free your results from bias influence by substracting a linear (or quadratic) regression of the integrated data from itself.

This is a part of my function in Matlab.

I hope, that it will be helpfull for someone.

[b,a]=butter(2,f/Fs,'low'); 

[d,c]=butter(2,f/Fs,'high');

a=filter(b,a,signal);

v=cumtrapz(t,a);

v=filter(d,c,v);

s=cumtrapz(t,v);

where:

f- signal frequency.

Fs- sampling frequency

signal - acceleration signal (input)

a - acceleration vector without noise (after filtering)

v,s,t - velocity, displacement and time vectors

Robert.

Robert what is the difference between the two frequencies 

As Robert stated, a Butterworth band-pass filter with the desired cut-off frequencies will solve the drift and the noise problems if designed properly. In his case, he applied two different filters, one high-pass and other low-pass to achieve the same result. The cut-off frequencies used where the same in both filters and the parameter f/Fs will be the normalized frequency of interest.

On the other hand, Gunter's method could also work in some cases, although a new undesired bias may appeard in other cases. Basically this will set the value of the signal to the real value it should have been measured.

Hi, i had similar problems. I think Butterworth band-pass filter works very well for my case. First, i used FFT to convert acceleration data from time domain to frequency domain. I figured out the cut-off frequency based on the FFT result. Afterwards, I applied a band-pass filter to take out the unwanted data. Then i did integration of filtered acceleration data, got velocity and displacement. Here is the matlab coding which i used. If you have any questions, you can email me by [email protected].

load('C:\Users\mm395\Documents\MATLAB\REC0819ch2.mat'); % load acceleration data

figure;plot(REC0819ch2(:,1),REC0819ch2(:,2));

title('Acceleration data');

fft_accel = fft(REC0819ch2(:,2)); % Compute the FFT of the acceleration data

L_accel = length(REC0819ch2(:,2)); % Calculate the number of the data

P_2 = abs(fft_accel/L_accel); % Compute the two-sided spectrum P2

P_1 = P_2(1:L_accel/2+1); % Compute the sing-sided spectrum P1

P_1(2:end-1) = 2*P_1(2:end-1);

Fs = 3200; %1/(137/L_accel);%

f_accel = Fs*(0:(L_accel/2))/L_accel;

figure;plot(f_accel,P_1);

title('Single-Sided Amplitude Spectrum of Acceleration Data');

xlabel('f (Hz)');

ylabel('|P1(f)|');

% %apply low pass filter

% fc = 100; %cut off frequency

% fs = Fs;

% [b,a] = butter(6,fc/(fs/2));

% freqz(b,a);

%

% %apply high pass filter

% [z,p,k] = butter(9,300/500,'high');

% sos = zp2sos(z,p,k);

% fvtool(sos,'Analysis','freq')

% %apply band pass filter

d = fdesign.bandpass('N,F3dB1,F3dB2',10,1,20,3200);

Hd = design(d,'butter');

% Frequency analysis of Filtered acceleration data

%filtered_accel = filter(b, a, REC0819ch2(:,2) );

%

filtered_accel = filter(Hd, REC0819ch2(:,2) );

figure;plot(REC0819ch2(:,1),filtered_accel);

title(' Filtered Acceleration data');

fft_accel_f = fft(filtered_accel); % Compute the FFT of the acceleration data

P_2_f = abs(fft_accel_f/L_accel); % Compute the two-sided spectrum P2

P_1_f = P_2_f(1:L_accel/2+1); % Compute the sing-sided spectrum P1

P_1_f(2:end-1) = 2*P_1_f(2:end-1);

Fs = 3200; %1/(137/L_accel);%

f_accel_f = Fs*(0:(L_accel/2))/L_accel;

figure;plot(f_accel_f,P_1_f);

title('Single-Sided Amplitude Spectrum of Filtered Acceleration Data');

%caculate the accumulative velocity

V_acce = trapz(filtered_accel);

AC_V_acce = cumtrapz(filtered_accel)*(1/Fs);

%T = table(time',AC_V_acce','VariableNames',{'Time','CumulativeDistance'});

figure; plot(REC0819ch2(:,1),AC_V_acce);

title('Speed');

xlabel('Time (s)');

ylabel('Velocity (m/s)');

% %apply low pass filter

% filtered_V_acce = filter(b, a, V_acce);

%

%

%caculate accumulative displacement

%D_acce = trapz(filtered_V_acce);

%AD_V_acce = cumtrapz(filtered_V_acce)*(1/Fs);

D_acce = trapz(AC_V_acce);

AD_V_acce = cumtrapz(AC_V_acce)*(1/Fs)*1000;

%T = table(time',AC_V_acce','VariableNames',{'Time','CumulativeDistance'});

figure; plot(REC0819ch2(:,1),AD_V_acce);

title('Displacement');

xlabel('Time (s)');

ylabel('Displacement (mm)');

dA = deltaT * x(n-1)

y(n) = (1 - g) * y(n-1) + (1 + g) *dA

x(n) input sequence

y(n) output sequence (integrated signal)

Where g is small and can be taken as R * deltaT

Where 1/R = integrator response time in sec

When you use 'cumtrapz' matlab function you will always have drift effect. In order to remove the drift, based in my experiments experience, I sugesst you to use a detrend operation (implemented in matlab 'detrend') and then a high pass filter with a specific and tested cuttoff frequency. I recomend to use Butterworth or Chebychev high pass filters.

Actually it s pretty hard to obtain velocity/position data from acceleration data. And the efficiency varies between different physical applications.

Let me summarize some might-be-useful methods to do that.

1 - When you integrate an acceleration data, you also integrate the noise existing in the signal. Therefore, you ll definitely have errors. Even the parameter might go to infinity due to the noise. To sum up, you should use an optimal Kalman Filter. And you have to optimize & calibrate it due to the application.

2 - You have to determine the boundaries & limitations you work with. For example, you may want to obtain -+2 mm position data.

3- You have to determine an optimal high pass filter.

I've tried a lot of methods, and I will try some more.

I always have some drifting in the beginning and end of the time signal, after high pass filtering, integration and high pass filtering again of at least 10Hz.

Right now, for the problem of estimating the displacements orbits from two orthogonal accelerometers, i´m discarding 10% of the samples in the beginning and end.

Im looking for a solution to do not have to discard samples.

The errors can have different sources.

The first and usually the most important as effect is the offset the transducer has. If you use an accelerometer you can easily find the offset value and deduct it from measurement by "measuring" the g-values with the accelerometer in + and - directions. The difference between the 2 readings is 2x the offset.

For a velocity transducer it is more difficult to calibrate it.

The different procedures are more or less interesting depending on the error level you can accept and tolerate.

As a general approach you can take successive "slices" of data and define a plynom of a 2nd or 3rd degree through 3 or 4 successive points. You can compute the integral for each slice and sum them. It is similar to the 1-4-1 Simpson rule but with variable distances between points. As an example with a trapezoïdal approach (linear through 2 points) the integral of sin(phy) from 0 to pi/2 is 0.9994 if you use a second order approach as described with NOT equal steps the obtained value is 1.0001087.

In general a second order is satisfactory if the steps are adjusted according to the curvature.

Another method for reducing noise and drift which hasn't been mentioned, is simply using more accelerometers and averaging their values before calculating velocity and position. The averaged noise should reduce due to cancelling of noise between each accelerometer. Obviously, it's not always possible to use more accelerometers, but if you can, it will improve the signal. You can still add filters etc to the averaged signal afterwards.

Manual detrend with Excel?

That's a question that's worth a couple of bucks! I think people don't explain this in books because they make money with it!? The answer is very simple and usually is not found in any software under the name or algorithm for Integral+some-method! The problem is related to integration of low frequency White Noise that yields a random walk. So the key idea is to use a low pass filter with a cutoff frequency chosen to remove the low frequency part with a fixed gain,G(s)=Tau/(Tau s +1). The implementation is even simpler, use a discrete first order filter using Jordan canonical form, i.e. y(k)= y(k-1) + alpha*(u(k-1) - y(k-1)), where alpha = h/Tau, where h=sampling interval, Tau is the filter time constant, f_c = 1/(2 pi Tau). You need 1 multiplication and two additions and get rid of the random walk if you need to implement in real time using a microcontroller. See attached a simulink file illustrating the problem and the solution. Remember: A low pass filter is a very good integrator with (90 degree phase shift) for frequencies well above the cutoff frequency. By the way, use a series capacitor in the analog front data acquisition as a high pass filter before applying the low pass filter as an integrator.