How to integrate acceleration into displacement data correctly?

Hi all,

I am trying to extract displacement data from acceleration dataset. Even after detrending and high-pass filtering the result does not look good.

Any hints?

Cod eattached.

Fs = 1000; % sampling freq.
Ts = 1 / Fs; % sampling period or time step 
td = 1; % signal duration in seconds
L = Fs * td; % signal length = number of total sampes taken 
t = Ts * (0:L-1); % time vector 

x = 2.* sin(2*pi*25*t); % cosine function with freq. of: 50 Hz 
figure(1);
plot(t, x);

Y = fft(x, L); 
P2 = abs(Y/L); % re-scaling  
P1 = P2(1 : L/2+1); 
P1(2 : end-1) = 2*P1(2 : end-1); 
f1 = (Fs/L)*(0:L/2);  
figure(2);
plot(f1, P1)

%% Integration
x_detr = detrend(x);
velnf = cumtrapz(t, x); 
dispnf = cumtrapz(t, velnf);
% De-trending 
x_vel = cumtrapz(t, x_detr);
x_vel_detr = detrend(x_vel);
x_disp = cumtrapz(t, x_vel_detr);
% after de-trending now it is the time for filtering high pass filter
%
figure(3);
plot(t, x, 'k');
hold on;
plot(t, x_disp, 'b');
hold on;
plot(t, dispnf, 'r');

% HP (High Pass) filtering 
N = 2; % filter order 
cutoff_freq = 1; % cutoff frequency in Hz
[butB, butA] = butter(N, cutoff_freq / (Fs/2), 'high');
accf2 = filter(butB, butA, x_detr); % butter w/ filter  
accf3 = filtfilt(butB, butA, x_detr); % butter w/ filtfilt 
velf2 = cumtrapz(t, accf2); 
velf2 = detrend(velf2); 
disp2 = cumtrapz(velf2); 
velf3 = cumtrapz(t, accf3); 
velf3 = detrend(velf3);
disp3 = cumtrapz(velf3);

figure(4); 
plot(t, x, 'k');
hold on;
plot(t, disp2, 'b');
hold on;
plot(t, disp3); 
hold off

Answer 1

[Paul]

The integration seems correct based on the implied initial conditions.

Fs = 1000; % sampling freq.
Ts = 1 / Fs; % sampling period or time step 
td = 1; % signal duration in seconds
L = Fs * td; % signal length = number of total sampes taken 
t = Ts * (0:L-1); % time vector 

The acceleration has frequency of 25 Hz, not 50 Hz.

x = 2.* sin(2*pi*25*t); % cosine function with freq. of: 50 Hz
velnf = cumtrapz(t, x); 
figure
plot(t,velnf)

dispnf = cumtrapz(t, velnf);
figure
plot(t,dispnf)

We see that cumtrapz starts the integration at 0 for both velnf and dispnf

Let's get the exact solutions for velocity and displacement, don't forget the constants of integration

vel = int(2*sin(2*sym(pi)*25*sym('t')),sym('t')) + sym('C1')

vel = 

disp = simplify(int(vel,sym('t'))) + sym('C2')

disp = 

In accordance with the cumtrapz, we have to select C1 and C2 such that vel(0) = disp(0) = 0;

C1 = solve(subs(vel,sym('t'),0))

C1 = 

C2 = solve(subs(disp,[sym('t') sym('C1')],[0 C1]))

C2 = 

0

Now plot

figure
fplot(subs(vel,sym('C1'),C1),[0 1])

figure
fplot(subs(disp,[sym('C1') sym('C2')],[C1 C2]),[0 1])

which are essentially the same plots as using cumtrapz.

[Alireza]

thanks Paul,

then I think cumtrapz is not the best/correct one for integrating at least in time domain.

maybe i need to use trapz and loop through all the time steps.

[Paul]

cumtrapz may be fine if properly accounting for initial conditions. Hard to say without more information on what the differerence is between the computed and expected results.

[Alireza]

problem is I need to double integrate accelearion data to get displacement data. but looks like there is a drift after integration, besides the fft of displacement is basically empty. not sure if the problem is for the cumtrapz or what ... i have de-trended my dataset and applied high pass filter. imagine instead of x as a sinusoid function i have such data attached:

[Paul]

As shown, the "drift" in the displacement is because the velocity has a constant offset, called C1 above, and that constant offset in velocity turns into a C1*t in the displacement, which is the drift.

If you think the displacement shouldn't have that drift you can always try to detrend.

Fs = 1000; % sampling freq.
Ts = 1 / Fs; % sampling period or time step 
td = 1; % signal duration in seconds
L = Fs * td; % signal length = number of total sampes taken 
t = Ts * (0:L-1); % time vector 
x = 2.* sin(2*pi*25*t); % cosine function with freq. of: 50 Hz
velnf = cumtrapz(t, x);
dispnf = cumtrapz(t, velnf);
dispnf = detrend(dispnf,1);
figure
plot(t,dispnf)

There might be other options, but, again, I have no idea what the expected result for displacement actually is.

I looked at ACC.mat, but had not idea what all the variables meant or how they are to be processed.