How to calculate differential of scattered data?

Hi, I have time sertis of two signals, one is y, one is Cl, and I'd like to calculate the first and second derivative of these two signals.

Now I use the following code (the dataset has been attached):

dy   = gradient(y)./gradient(time);
dCl  = gradient(Cl)./gradient(time);
ddy  = gradient(dy)./gradient(time);
ddCl = gradient(dCl)./gradient(time);

But when I use the data to plot phase-plane, the results look weird.

X1 = [y, dy, ddy];
X2 = [Cl, dCl, ddCl];
figure(1)
patch([X1(:,1);nan], ...
      [X1(:,2);nan], ...
      [X1(:,3);nan], ...
      [X1(:,1)+X1(:,2);nan],...
       'EdgeColor','interp','Marker','none','MarkerFaceColor','flat','LineWidth',0.8,'FaceAlpha',1);
view(3);

figure(2)
patch([X2(:,1);nan], ...
      [X2(:,2);nan], ...
      [X2(:,3);nan], ...
      [X2(:,1)+X2(:,2);nan],...
       'EdgeColor','interp','Marker','none','MarkerFaceColor','flat','LineWidth',0.8,'FaceAlpha',1);
view(3);

I'd like to determine the intersection between curve of y-dy (down figure) or Cl-dCl (up figure) with z=0 plane, but the ddCl varies above and below frequently around z=0, so the intersection cannot be identified, like figures attached.

Is there any problem during the data processing procedure?

Thanks!

[Askic V]

Well, you need to use some kind of filter. If you just do first order approximation of derivation (also called the Euler method), it will just increase the noise in the data.

Please have a look at this small example:

dt = 0.1; % step size (sec)
t_dur = 10; % total t_dur
t = 0:dt:t_dur;
y_pure = sin(t);
std_dev = 0.1;
y_noise = y_pure + std_dev*randn(size(t));
subplot(4,1,1)
plot(t,y_pure)
title('Pure signal')

dy= gradient(y_pure)./gradient(t);
subplot(4,1,2)
plot(t,dy)
title('Derivation of pure signal');
subplot(4,1,3)
plot(t,y_noise)
title('Signal with noise');
dy_noise = gradient(y_noise)./gradient(t);
subplot(4,1,4)
plot(t,dy_noise)
title('Derivation of noisy signal');

As you can see, you are amost not able to recognize that the last signal is cos(t).

So if the data contains noise (and most likely it is always the case), some kind of filtering must be applied.

I think most propriate for these applications would be kalman filtering, but I hope some more experience will jump to help.

Answer 1

[Askic V]

Okay, here's what I do in order to avoid calculating derivation on the discrete data, especially when I know that they represent position, velocity and acceleration.

In fact, I used very similar code (Kalman filter implementation), when I had real measured position values (using MDT sensors), but I needed to calculate velocity and acceleration in order to use feed forward control.

dt = 0.1; % step size (sec)
t_dur = 15; % total t_dur
t = 0:dt:t_dur;
y_pure = sin(t);
dy_pure = cos(t);
ddy_pure = -sin(t);

std_dev = 0.05; % Std. deviation of Gauss noise
y_noise = y_pure + std_dev*randn(size(t)); % add noise to the signal

% Calculate velocity and acceleration
% using gradient or diff on the discrete data
dy_noise = gradient(y_noise)./gradient(t);
ddy_noise = gradient(dy_noise)./gradient(t);

dt = 0.1; % step size (sec)
t_dur = 15; % total t_dur
t = 0:dt:t_dur;
y_pure = sin(t);
dy_pure = cos(t);
ddy_pure = -sin(t);

std_dev = 0.05; % Std. deviation of Gauss noise
y_noise = y_pure + std_dev*randn(size(t)); % add noise to the signal

% Calculate velocity and acceleration
% using gradient or diff on the discrete data
dy_noise = gradient(y_noise)./gradient(t);
ddy_noise = gradient(dy_noise)./gradient(t);

% plot data
subplot(3,1,1)
plot(t,y_pure);
hold on
plot(t,y_noise,'r')
legend('pure signal', 'noisy signal');

subplot(3,1,2)
plot(t,dy_pure)
hold on
plot(t, dy_noise)
legend('pure velocity', 'calculated (noisy) velocity');

subplot(3,1,3)
plot(t,ddy_pure)
hold on
plot(t,ddy_noise)
legend('pure acc', 'calculated (noisy) acc');

% 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
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;1;1];

% Design filter
P = P0;
N = length(y_noise);
xhat_array = 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;
    xhat_array(:,i) = x;
end
figure
subplot(2,1,1)
% Plot noisy velocity
plot(t, dy_noise)
hold on
% Plot filtered velocity
plot(t, xhat_array(2,:),'r');
% Plot ideal velocity
plot(t, dy_pure,'g');
legend('Noisy velocity', 'Filered velocity', 'Ideal velocity');
% Plot noisy acceleration

subplot(2,1,2);
plot(t, ddy_noise)
hold on
% Plot filtered acceleration
plot(t, xhat_array(3,:),'r');
hold on
% Plot ideal acceleration (pure)
plot(t, ddy_pure,'g');
legend('Noisy acc', 'Est./Filt. acc', 'Ideal acc.');

Basically, instead of calculating derivatives using the Euler method (or by definition), you should use estimation techniques such as Kalman filtering. Kalman filter is optimal linear state estimator. You can see from the example, how even small amount of noise on the original data can practically totally corrupt the calculated derivations (especially second derivation i.e. acceleration).

Derivations always tend to amplify the noise.

If you're interested I suggest searching for prof. Dan Simon and his excellent book "Optimal State Estimation".

I hope this will be useful to you (and probably others).

[Ying Wu]

@Askic V Many thanks for providing such an detailed explaination and nice example. My signal y also stands for position, velocity and acceleration, although my position data is not measured but simulated. I would check the Kalman filtering method. I notice some literature calculate the derivative of discrete data using the difference method with second order accurate and then smooth it. Is this equivalent to using filter? Thanks!

[Askic V]

@Ying Wu i didn't check your data, but if they're simulated, then, most probably you will not be able to see a difference between calculating derivatives directly and using Kalman filter (unless you added noise, like I had).

In general, you could try to generate velocity and acceleration by calculating derivatives using the Euler method, and then smoothing by using moving average or different filters, but unless you are very proficient in filter design the result would be most likely worse than using Kalman filter.

Remember, Kalman filter is an optimal state estimator for linear systems. It indeed peforms smoothing (filtering) as can be seen from the plots, but usually applying separete filter after more noisy data is generated (noise amplified through numerical derivation) will not achieve better result.

[Ying Wu]

@Askic V Thanks for reply! I am sorry about the unclear description of my data. In fact, the "simulated" here not refers to numerically solving ode system (with explicit expressions) like using ode45() function in matlab. My simulations results are obtained from the computational fluid dynamics software, so the direct output is the postiion of moving body and other related quantities in flow field, which is similar to the measured position in experiments. So I think my data should also be noisy and needs to be filtered.

[Ying Wu]

@Askic V Hi Askic, I have double check my two signals, and use gradient() to take first and second-order derivatives without doing any filtering. The results are shown below.

signal 1: y, dy, ddy

signal 2: Cl, dCl, ddCl

My direct thinking is that the singal 2 itself seems not very noisy, but its second derivative of looks very noisy, so how can I use filter to solve this problem? Thanks!