How to take first and second order derivative of discrete data?
73 vues (au cours des 30 derniers jours)
Afficher commentaires plus anciens
Hi, I have two discrete signals with constant time intervel dt=0.01, one is y(t) (displacement), one is Cl(t) (force). Each signal is a vector with 200000 elements, and I don't know the function expression of them, so I use numerical difference technique in matlab to take first and second derivatives. I use the second-order accuracy central difference method, i.e. gradient() function in matlab.
The code is:
load t.mat y.mat Cl.mat
% take 1st and 2nd derivative of y(t)
dy = gradient(y)./gradient(t);
ddy = gradient(dy)./gradient(t);
% take 1st and 2nd derivative of Cl(t)
dCl = gradient(Cl)./gradient(t);
ddCl = gradient(dCl)./gradient(t);
% plot y, y', y''
subplot(3,2,1)
plot(t,y,'k-'); ylabel('y');
subplot(3,2,3)
plot(t,dy,'r-'); ylabel('y\prime');
subplot(3,2,5)
plot(t,ddy,'r-'); ylabel('y\prime\prime');
% plot Cl, Cl', Cl''
subplot(3,2,2)
plot(t,Cl,'k-'); ylabel('Cl');
subplot(3,2,4)
plot(t,dCl,'r-'); ylabel('Cl\prime'); xlim([200 2000])
subplot(3,2,6)
plot(t,ddCl,'r-'); ylabel('Cl\prime\prime'); xlim([200 2000])
The results is shown below.
![](https://www.mathworks.com/matlabcentral/answers/uploaded_files/1233022/image.png)
The direct numerical difference of signal 1 (y) seems to obtain good results, while the second derivative of signal 2 shows weired shapes, it is not like noise but with some repeated pattern.
How can I filter or smooth the second signal to obtain correct derivatives? Ang suggestion is welcome! Thanks!!!
Réponse acceptée
Askic V
le 16 Déc 2022
Hello Ying,
you wrote that second derivative of signal 2 looks weird. I assume that is ddCl (also based on how plot looks like).
I can offer applying Kalman filtering function I wrote in earlier topic:
load t.mat; load y.mat; load Cl.mat
% take 1st and 2nd derivative of y(t)
dy = gradient(y)./gradient(t);
ddy = gradient(dy)./gradient(t);
% take 1st and 2nd derivative of Cl(t)
dCl = gradient(Cl)./gradient(t);
ddCl = gradient(dCl)./gradient(t);
% use kalman filter
dt = mean(diff(t)); % sample time dt = t(2)-t(1) if constant
Cl_filt = kalm_estim(Cl, dt);
Cl_k = Cl_filt(1,:).'; % Kalman estimation of original signal
dCl_k = Cl_filt(2,:).'; % Kalman estimation of 1st deriv
ddCl_k = Cl_filt(3,:).'; % Kalman estimation of 2nd deriv
% plot y, y', y''
subplot(3,2,1)
plot(t,y,'k-'); ylabel('y');
subplot(3,2,3)
plot(t,dy,'r-'); ylabel('y\prime');
subplot(3,2,5)
plot(t,ddy,'r-'); ylabel('y\prime\prime');
% plot Cl, Cl', Cl''
subplot(3,2,2)
plot(t,Cl,'k-'); ylabel('Cl');
subplot(3,2,4)
plot(t,dCl,'r-'); ylabel('Cl\prime'); xlim([200 2000])
subplot(3,2,6)
plot(t,ddCl,'r-'); ylabel('Cl\prime\prime'); xlim([200 2000])
figure (2)
subplot(3,2,1)
plot(t,Cl,'k-'); ylabel('Cl'); xlim([200 400]); ylim([-2 2])
subplot(3,2,3)
plot(t,dCl,'r-'); ylabel('Cl\prime'); xlim([200 400]); ylim([-0.8 0.8])
subplot(3,2,5)
plot(t,ddCl,'r-'); ylabel('Cl\prime\prime'); xlim([200 400]); ylim([-15 15])
% plot Cl_k, Cl_k', Cl_k''
subplot(3,2,2)
plot(t,Cl_k,'k-'); ylabel('Cl_k'); xlim([200 400]); ylim([-2 2])
subplot(3,2,4)
plot(t,dCl_k,'r-'); ylabel('Cl_k\prime'); xlim([200 400]); ylim([-0.8 0.8])
subplot(3,2,6)
plot(t,ddCl_k,'r-'); ylabel('Cl_k\prime\prime'); xlim([200 400]); ylim([-0.5 0.5])
function y_filt = kalm_estim(y, dt)
y_noise = y;
% Now implement and apply Kalman filter
T = dt;
% State space model of a system
% dx/dt = A*x, y = C*x;
A = [1 T T^2/2; 0 1 T; 0 0 1];
C = [1 0 0];
% process variance
Q = 0.1*[T^4/4 T^3/2 T^2/2; T^3/2 T^2 T; T^2/2 T 1];
% sensor noise variance
std_dev = 0.1;
R = std_dev^2;
% initial state estimate variance
P0 = [ 1 0 0; 0 1 0;0 0 1];
% Initial value of state vector
% x(1) - position
% x(2) - velocity
% x(3)- accelerartion
x = [0;0;0];
% Design filter
P = P0;
N = length(y_noise);
y_filt = zeros(3, N);
for i = 1:N
xp = A*x;
Pp = A*P*A' + Q;
S = C*Pp*C' + R;
K = Pp*C' / S;
residual = y_noise(i) - C*xp;
x = xp+K*residual;
P = Pp - K*C*Pp;
y_filt(:,i) = x;
end
end
This code will produce results shown in the attached picture. This looks more like one would expect.
![](https://www.mathworks.com/matlabcentral/answers/uploaded_files/1233052/image.png)
3 commentaires
Askic V
le 16 Déc 2022
Kalman filter is an optimal state estimator. If the system is non linear, the filter is not optimal, but it doesn't mean it cannot provide good results and it can be near optimal. In any case, it can produce better result than trying to filter out high frequency noise by trial and error design. It is my first try when I have measurements which I need to use to estimate/calculate higher order quantities. It is especially good for estimating velocity and acceleration from position measurements.
Of course, there are Nonlinear Kalman filters such as Extended Kalman Filter, but I don't have much experience with it.
Plus de réponses (0)
Voir également
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!