Azzera filtri
Azzera filtri

How to integrate acceleration into displacement data correctly?

11 visualizzazioni (ultimi 30 giorni)
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

Risposte (1)

Paul
Paul il 2 Apr 2024
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.
  4 Commenti
Alireza
Alireza il 2 Apr 2024
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
Paul il 2 Apr 2024
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.

Accedi per commentare.

Community Treasure Hunt

Find the treasures in MATLAB Central and discover how the community can help you!

Start Hunting!

Translated by