Finding formulation of curve using coordinates in MATLAB

Question

Suppose we are given some coordinates of points in plane, say 20 points with their coordinates as:

X = [0.0809627469389454;
0.0812319917473947;
0.0814760690888197;
0.0816493550701358;
0.0817167028412858;
0.0816873091729635;
0.0813456103002472;
0.0805960479373602;
0.0793506286325874;
0.0768532682616457;
0.0730568739618012;
0.0688299169022524;
0.0629775052745549;
0.0555093488004414;
0.0469410412695979;
0.0385912069617207;
0.0296789727506768;
0.0200899022102430;
0.00973217977333483;
-0.00138286576734164]

and

Y=[-0.737109973623757;
-0.737164532535452;
-0.737292603740519;
-0.737517126138711;
-0.737930464518794;
-0.738537206858314;
-0.740008741381755;
-0.742712757707760;
-0.745966011330839;
-0.751634798125669;
-0.759200611372096;
-0.767277401090950;
-0.777095158716249;
-0.787200839430418;
-0.797371590291212;
-0.806336086919729;
-0.816124045611269;
-0.825660302601719;
-0.835939409428708;
-0.846471284598907]

These points are a movement coordinates and it is on a curve which the formulation is unknown. How one could find the formulation of the movement using given coordinates/data? I need its formulation to calculate its speed and acceleration which are first and second derivatives of the movement, respectively. Thanks in advance!!

Answer 1

I guess you are more interested in the speed and acceleration values instead of a curve fit of the original data, but you assume that it is easier to estimate speed and acceleration by differentiation of a mathematical function instead of working with real-world data?

If my guess is correct let me tell you that there is a better way than a rough approximation of the measured data because you throw away too much valuable information at the beginning of your estimations. The result would be strong deviations of the final speed and acceleration values.

Just simply estimate the discrete derivatives of the measured X- and Y-displacements and finally try a smoothing to reduce the sample errors. I assume that the X/Y-pairs are sampled with a constant time interval so I use this interval as the new X of the following six plots. You should replace it by the real time values (e.g. milliseconds). You can use this as a skeleton for your own analysis:

subplot(4,1,1)
plot(X,Y,'.-')
xlabel('Measured displacement data of the hand movement')

speedX = diff(X)
speedY = diff(Y)
speed = sqrt( speedX.^2 + speedY.^2 );

subplot(4,2,3)
plot(speedX)
ylabel('X')
subplot(4,2,5)
plot(speedY)
ylabel('Y')
subplot(4,2,7)
plot(speed)
ylabel('Combined X and Y')
xlabel('Speed')

accelerationX = diff(speedX)
accelerationY = diff(speedY)
acceleration = sqrt( accelerationX.^2 + accelerationY.^2 );

subplot(4,2,4)
plot(accelerationX)
subplot(4,2,6)
plot(accelerationY)
subplot(4,2,8)
plot(acceleration)
xlabel('Acceleration')

% Simple solution to remove the outliers
x1 = 1:length(speed);
p1=polyfit(x1',speed,5)
estSpeed=polyval(p1,x1)
subplot(4,2,7)
hold on
plot(x1,estSpeed,'r')
hold off

x2 = 1:length(acceleration);
p2=polyfit(x2',acceleration,4)
estAcceleration=polyval(p2,x2)
subplot(4,2,8)
hold on
plot(x2,estAcceleration,'r')
hold off

The following picture shows the resulting plots:

You will notice that even with these high polynomial orders I chose, the curve fit is still not really good. Lower orders really disturb the interesting trends in the speed and acceleration curves. So, it is up to you to find a better approximation method ;-)

Answer 2

By plotting the data, this looks polynomial. As such, I suggest you use the polyfit function. What this does is given a set of co-ordinates x and y, you specify what order of polynomial you believe the data best matches, and it finds the coefficients of the polynomial equation that best fits this data. This is performed by least-squares error minimization, but I'll leave that to you to read up on it. The mathematics are actually quite elegant.

From your data, I'm going to guess a third-order polynomial because of the inflection point at around x = 0.04. Therefore, all you have to do is this:

coeff = polyfit(X, Y, 3); %// Specify third order polynomial

The output we get is:

coeff =

  174.5793  -13.4016    1.2463   -0.8457

The coefficients for the equation are structured in descending order, starting from the third power, down to the intercept value. As such, this equation is actually:

Y = 174.5793*X^3 - 13.4016*X^2 + 1.2463*X - 0.8457

To check to see if this equation is right, we can do a visual inspection. We specify a bunch of points between the minimum and maximum X values, then use these points to evaluate the above function and see where this lines up on our graph. As such, we can do this:

xPoints = linspace(min(X), max(X)); %// Specify 100 points between the min and max
yPredict = polyval(coeff, xPoints); %// Find the corresponding y-values
plot(X,Y,'b.',xPoints,yPredict,'r'); %// Plot the original points with the predicted
grid;                                %// curve in red.

Let's go through the above code slowly. linspace creates a linear spacing of vectors between two points. The default amount of points is 100, but you can change this to whatever you want. As such, the first line will generate 100 points between the minimum and maximum points in X. After, polyval takes the coefficients of the polynomial equation that you created from polyfit, and a set of X points to generate a set of Y points based on the predicted equation we found. Once I do this, we then generate a plot that shows both the original data points and what the predicted curve looks like. I also throw in a grid for good measure.

This is the plot I get:

Obviously, you need to tune the order of the polynomial depending on your data. Bear in mind that we were lucky this time. In practice, data does not naturally flow like this. As such, you'll need to examine your data first and figure out what class of curve the data fits. Once you decide this, you proceed with the curve fitting as each class of curve follows a different strategy to find the best estimate of parameters.