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

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

[Ying Wu]

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

Answer 1

[Askic V]

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.

[Ying Wu]

@Askic V Hi Askic, thank you so much! The code seems to give a better filtering of singal 2. I have a quick question after checking the code and the principles of Kalman filter. It says this Kalman filter can only be used for linear system, since current state depends on state at last time and the current input. My question is if my signal is obtained from numerically solving Navier-Stokes equations (a nonlinear ODE system for fluids), then can the Kalman filter be used to smooth it?

[Askic V]

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]

@Askic V Thanks for your detailed explanation and information!