How to calculate differential of scattered data?

Question

Ying Wu 2022-12-15

0
链接

此问题的直接链接

https://ww2.mathworks.cn/matlabcentral/answers/1878657-how-to-calculate-differential-of-scattered-data

编辑： Ying Wu 2022-12-16

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!

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

Askic V 2022-12-15

在 MATLAB Online 中打开

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 2022-12-15

0
链接

此回答的直接链接

https://ww2.mathworks.cn/matlabcentral/answers/1878657-how-to-calculate-differential-of-scattered-data#answer_1128457

编辑：Askic V 2022-12-15

在 MATLAB Online 中打开

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).

4 个评论
显示 2更早的评论隐藏 2更早的评论

Ying Wu 2022-12-15

@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 2022-12-16

编辑：Ying Wu 2022-12-16

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

请先登录，再进行评论。

How to calculate differential of scattered data?

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

采纳的回答

4 个评论
显示 2更早的评论隐藏 2更早的评论

更多回答（0 个）

另请参阅

类别

标签

产品

版本

Community Treasure Hunt

How to calculate differential of scattered data?

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

采纳的回答

4 个评论 显示 2更早的评论隐藏 2更早的评论

更多回答（0 个）

另请参阅

类别

标签

产品

版本

Community Treasure Hunt

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

4 个评论
显示 2更早的评论隐藏 2更早的评论