How to integrate acceleration into displacement data correctly?

12 次查看(过去 30 天)
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

回答(1 个)

Paul
Paul 2024-4-2
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 个评论
Alireza
Alireza 2024-4-2
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 2024-4-2
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.

请先登录,再进行评论。

类别

Help CenterFile Exchange 中查找有关 Particle & Nuclear Physics 的更多信息

Community Treasure Hunt

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

Start Hunting!

Translated by