How to take first and second order derivative of discrete data?

Question

Ying Wu 2022-12-16

1
链接

此问题的直接链接

https://ww2.mathworks.cn/matlabcentral/answers/1879472-how-to-take-first-and-second-order-derivative-of-discrete-data

评论： Ying Wu 2022-12-16

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.

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!!!

1 个评论
显示 -1更早的评论隐藏 -1更早的评论

Ying Wu 2022-12-16

@Askic V Hi Askic, could you please take a look at this more detailed case? Thanks for your help!

请先登录，再进行评论。

请先登录，再回答此问题。

Answer 1

Askic V 2022-12-16

1
链接

此回答的直接链接

https://ww2.mathworks.cn/matlabcentral/answers/1879472-how-to-take-first-and-second-order-derivative-of-discrete-data#answer_1129267

在 MATLAB Online 中打开

Kalman filter.png

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.

3 个评论
显示 1更早的评论隐藏 1更早的评论

Askic V 2022-12-16

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.

Ying Wu 2022-12-16

@Askic V Thanks for your detailed explanation and information!

请先登录，再进行评论。

How to take first and second order derivative of discrete data?

1 个评论
显示 -1更早的评论隐藏 -1更早的评论

采纳的回答

3 个评论
显示 1更早的评论隐藏 1更早的评论

更多回答（0 个）

另请参阅

类别

标签

产品

版本

Community Treasure Hunt

How to take first and second order derivative of discrete data?

1 个评论 显示 -1更早的评论隐藏 -1更早的评论

采纳的回答

3 个评论 显示 1更早的评论隐藏 1更早的评论

更多回答（0 个）

另请参阅

类别

标签

产品

版本

Community Treasure Hunt

1 个评论
显示 -1更早的评论隐藏 -1更早的评论

3 个评论
显示 1更早的评论隐藏 1更早的评论