Strange form of increasing when derivating sine waves

2 vues (au cours des 30 derniers jours)
Bruce Rogers
Bruce Rogers le 9 Juin 2021
Commenté : Bruce Rogers le 11 Juin 2021
Hello everyone, I'm trying to get coordinates for a robot movement from a sine wave. Now when I change the duration of a period, the coordinates seem to drift, as shown in the picture figure 5. Can anyone explain? Or is there a mistake in my code? Thanks for your help, I'm thankful for every idea.
I'm using Matlab R2020b academic use.
clear;
b = 4; %Factor for changing the duration of period
t = 0.05; % Frequenz
A = 50; %Amplitude
x = 0:t:10*pi;
y = A*sin(b*x);
T = (1:length(x))*t;
figure(1)
plot(T,y,'.')
xlabel('time')
ylabel('distance')
%%
vel = gradient(y)./gradient(x);
figure(2)
plot(T,vel,'.')
xlabel('time')
ylabel('velocity')
%%
acc = gradient(vel)./gradient(x);
figure(3)
plot(T,acc,'.')
xlabel('time')
ylabel('acceleration')
acc = acc';
%xlswrite('Acceleration.xlsx',acc);
acc = acc';
%%
%velocity Cosinus Figure
v = zeros(1,length(acc));
v(1,1) = b*A * cos(0);
l_v = 1:length(acc);
for i = 2:length(acc)
v(i) = ((acc(i)+acc(i-1))/2) * t + v(i-1);
end
figure(4)
scatter(T,v,'.')
xlabel('time')
ylabel('velocity')
v = v';
%xlswrite('Velocity.xlsx',v);
v = v';
%%
%distance Sinus Figure
s = zeros(1,length(v));
s(1,1) = b*sin(0);
for i = 2:length(v)
s(i) = ((v(i)+v(i-1))/2) * t + s(i-1);
end
figure(5)
scatter(T,s,'.')
xlabel('time')
ylabel('distance')
%%
s = s';
%xlswrite('idealwaves.xlsx',s)
  4 commentaires
Jan
Jan le 10 Juin 2021
The term (v(i)+v(i-1)) / 2 cannot be the maximum value of v. It is likely that this has more influence on the positive part than on the negative part (or vice versa). Therefore such drifts are typical for numerical integrations of accelerations or velocities. You are not the first person, who observes this.
There is no magic solution. To reduce the influence of drifting, you need absolute position values as a calibration.
J. Alex Lee
J. Alex Lee le 10 Juin 2021
by the way you should be able to accomplish your numerical integration with cumtrapz, instead of manually summing
s = b*sin(0) + cumtrapz(T,v)
If you want to understand the "drift" you can systematically reduce your value of t (sample faster) to see how well your integrated position matches the analytical expectation.

Connectez-vous pour commenter.

Réponse acceptée

Paul
Paul le 10 Juin 2021
I think you're seeing the effect of the discrete time step combined with numerical approximations to the differentiation and the integration. Try decreasing t. Also, I noticed few other things in the code that may be of interest.
clear;
b = 4; %Factor for changing the duration of period
t = 0.05; % Frequenz % time step
A = 50; %Amplitude
x = 0:t:10*pi;
y = A*sin(b*x);
%T = (1:length(x))*t;
T = (0:length(x)-1)*t; % corrected, why use T and not x, x is the independent variable?
figure(1)
plot(T,y)
xlabel('time')
ylabel('distance')
%%
%vel = gradient(y)./gradient(x);
vel = gradient(y,t);
% vel = b*A*cos(b*x); % use this for better accuracy
figure(2)
plot(T,vel)
xlabel('time')
ylabel('velocity')
%%
% acc = gradient(vel)./gradient(x);
acc = gradient(vel,t);
% acc = -b^2*A*sin(b*x); % use this for better accuracy
figure(3)
plot(T,acc)
xlabel('time')
ylabel('acceleration')
acc = acc';
%xlswrite('Acceleration.xlsx',acc);
acc = acc';
%%
%velocity Cosinus Figure
v = zeros(1,length(acc));
% v(1,1) = b*A * cos(0);
v(1,1) = b*A * cos(b*0);
l_v = 1:length(acc);
for i = 2:length(acc)
v(i) = ((acc(i)+acc(i-1))/2) * t + v(i-1);
end
figure(4)
plot(T,v,T,vel,'o')
xlabel('time')
ylabel('velocity')
v = v';
%xlswrite('Velocity.xlsx',v);
v = v';
%%
%distance Sinus Figure
s = zeros(1,length(v));
% s(1,1) = b*sin(0);
s(1,1) = A*sin(b*0); % amplitude is A
for i = 2:length(v)
s(i) = ((v(i)+v(i-1))/2) * t + s(i-1);
end
figure(5)
plot(T,s,T,y,'o')
xlabel('time')
ylabel('distance')
figure(6)
plot(T,s-y) % error
  1 commentaire
Bruce Rogers
Bruce Rogers le 11 Juin 2021
Hey Paul, that's great, thank you very much! I like the last figure, it shows the error when the drift appears

Connectez-vous pour commenter.

Plus de réponses (0)

Produits

Community Treasure Hunt

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

Start Hunting!

Translated by