from acceleration to velocity: problem of integration

Hello everybody,

I explain the situation: I did a simple movement with my inertial measurement unit (IMU). I slide it manually along the x direction 2 times (once along -x and then along + x). No other movement has been realised along the 2 other axes (y and z). NB: z is vertical and x and y are in the horizontal plane.

Then I get the acceleration: (the gravity g along z is removed)

Then I integrate the acceleration to get the velocity:

for i = 2:length(acc)
    linVel(i,:) = linVel(i-1,:) + linAcc(i,:) * samplePeriod;
end

This is what I get:

Then, I filter with a high pass filter(butterworth, cut off freq =0,1 Hz) to remove the drift. I get:

order = 1;
filtCutOff = 0.1;
[b, a] = butter(order, (2*filtCutOff)/(1/samplePeriod), 'high');
linVelHP = filtfilt(b, a, linVel);

But there is still a problem... Why do I have 3 "trends" on this x curve? How do you explain it? Can we correct it? It may be a problem if I want the position...

And why do we use a high pass filter? Intuitively, I would have used a low pass filter but I saw somebody using a high pass filter (here: https://github.com/xioTechnologies/Oscillatory-Motion-Tracking-With-x-IMU ) to get the velocity from the acceleration...

NB: when I plot the pwelch (power spectral density) of the accel, there is nothing above 0,8 Hz... here:

Thanks to help me !!

[Minialoe]

any idea?

Answer 1

[Are Mjaavatten]

I have followed this question, hoping that somebody with better signal processing skills that myself would come up with an elegant solution. Since this has not happened, I decided to give it a try. My solution is not very general, but seems to solve your specific problem.

I make use of a run length encoding routine, commonly used for compacting graphics. There are many such routines on the File Exchange. I use the function rle by Said Bourezg.

The philosophy is to set the value of the acceleration (x) to zero whever it stays below a given threshold for a specified number of samples.

function x = removebias(x,biaslevel, period)
% REMOVEBIAS: Set x to zero whenever values stay below biaslevel
%   for a a minimum period
% Uses run length encoder rle by Said BOUREZG, from Mathworks File Exchange
    b = rle(abs(x) < biaslevel);
    nperiods = length(b)/2;
    b = reshape(b,2,nperiods);
    pos = [1,cumsum(b(2,:))];
    for i = 1:nperiods
        if b(1,i) && pos(i+1)-pos(i) >= period
            x(pos(i):pos(i+1)) = 0;
        end
    end
end

Testing with synthetic data:

% Create sample data:
amplitude = 8;          % Spike amplitude
bias = [0.5,0.1,0.3]; 
noise = 0.1;  % Noise amplitude
spike = zeros(20,1);
spike(1:10)= linspace(0,amplitude,10);
spike(11:20) = linspace(amplitude,0,10);
x =zeros(1500,1);
x(1:500) = bias(1)+noise*randn(500,1);
x(501:1000) = bias(2)+noise*randn(500,1);
x(1001:1500) = bias(3)+noise*randn(500,1);
x(491:510) = spike;
x(991:1010) = -spike;
y = removebias(x,1,5);
plot([x,y])

[Minialoe]

Hi, your algorithm seems to work on my velocity (instead of 1.5 for the biaslevel, I put 0.2 and it deletes the slopes... but it erodes also a little bit the base of my spikes). When I integrate one more time, what I get, the position, seems to be not too bad.

But it is of course not the best "signal processing" solution... Thanks anyway

[Are Mjaavatten]

I agree it is an ad-hoc solution, but it may be good enough in some cases.

I noticed that I did not submit my final version of the removebias function, so I have made an edit. The function now contains a "period" argument that ensures that shorter periods of low values are not automatically set to zero. Setting the period to 1 will give the same answer as before,

Answer 2

[Minialoe]

hello, I want to know if you want other elements or a clarification to answer my question ? thanks