MATLAB Answers

Spline interpolation for angular data

조회 수: 4(최근 30일)
mukesh bisht
mukesh bisht 2021년 8월 23일
댓글: William Rose 2021년 8월 24일
Hi,, I have 4 set of data points for which I am using spline interpolation to obtain closed curves. Can somebody explain why the two curves are not closed smoothly (i.e. the magnetta and red curve) even though I apply similar procedure for all curves. I have attached the data file, the figure and my Matlab code.
X1 = xlsread('data_points.xlsx', 'Sheet1','A2:A26'); Y1 = xlsread('data_points.xlsx', 'Sheet1','B2:B26');
X2 = xlsread('data_points.xlsx', 'Sheet1','D2:D26'); Y2 = xlsread('data_points.xlsx', 'Sheet1','E2:E26');
X3 = xlsread('data_points.xlsx', 'Sheet1','G2:G26'); Y3 = xlsread('data_points.xlsx', 'Sheet1','H2:H26');
X4 = xlsread('data_points.xlsx', 'Sheet1','J32:J26'); Y4 = xlsread('data_points.xlsx', 'Sheet1','K2:K26');
YY1 = [X1,Y1]; YY2 = [X2,Y2]; YY3 = [X3,Y3]; YY4 = [X4,Y4];
YY1 = transpose(YY1); YY2 = transpose(YY2); YY3 = transpose(YY3); YY4 = transpose(YY4);
XX = pi*[0:1/12:2];
% Periodic interpolation of cubic spline curve
pp1 = spline(XX,YY1); pp2 = spline(XX,YY2); pp3 = spline(XX,YY3); pp4 = spline(XX,YY4);
yy = ppval(pp1, linspace(0,2*pi,101));
plot(yy(1,:),yy(2,:),'-g'); axis equal; hold on
yy = ppval(pp2, linspace(0,2*pi,101));
plot(yy(1,:),yy(2,:),'-m'); hold on
yy = ppval(pp3, linspace(0,2*pi,101));
plot(yy(1,:),yy(2,:),'-r'); hold on
yy = ppval(pp4, linspace(0,2*pi,101));
plot(yy(1,:),yy(2,:),'-k'); hold on


William Rose
William Rose 2021년 8월 23일
The shape of the fitted spline at each point depends on points before and after that point. The curves do mot match smoothly at the ends because there are no points before point 1 and no points after point N (where N=last point in the list). To get a smooth wraparound, make a copy of hte first 5 points and tack them on after points 1 through N. You could do this in the Excel workbook or in Matlab. Suppose you do it in the Excel workbook. Then, if N1=25 in the original workbook, N2=30 in the extended version, and points 26-30 will be the same as points 1-5 respectively. The you do the spline fit on all 30 points.. Then you plot the fitted spline form points 3 to N1+3=28. This will close the curve, since point 3=point 28, and it should be smooth, since the point before and after point 3 (point 28) ar the same at the beginning as they are at the end.
Good luck!

William Rose
William Rose 2021년 8월 24일
I suspect my earlier proposal won't work as I suggested, because the x-values passed to spline() must be unique. That means you cannot pass two point wiht the same x value. But you coud add 2Pi to some of them. Still , my earlier suggesiton may not be the best way. I'll get back to that later. First:
When I ran your code, I got an error:
>> closedSpline1
Error using horzcat
Dimensions of arrays being concatenated are not consistent.
Error in closedSpline1 (line 8)
YY1 = [X1,Y1]; YY2 = [X2,Y2]; YY3 = [X3,Y3]; YY4 = [X4,Y4];
This happens because your line 4 of your posted code is
X4 = xlsread('data_points.xlsx', 'Sheet1','J32:J26'); Y4 = xlsread('data_points.xlsx', 'Sheet1','K2:K26');
so I corrected it to
X4 = xlsread('data_points.xlsx', 'Sheet1','J2:J26'); Y4 = xlsread('data_points.xlsx', 'Sheet1','K2:K26');
and it ran without error. I made some other changes: do the transposing at the same time as the concatenating, and plot the original oints as well as the splines, and add a legend. See attached code. This does not address your original question though. But I'm getting there.
  댓글 수: 1
William Rose
William Rose 2021년 8월 24일
Now I see that the last point in your spreadsheet is a duplicate of the first point.
I also see that there is a potentially confusing subtlety when fitting angular data as it is presented in the manual example and in your problem. The potentially confusing thing is that "x" has two different meanings. "x" is vector of phase angles, which we can think of as "theta", the vector Pi*(0:1/12:2) for example. X is also the vector of x-coordinates of the points associated with those phase angles: X1, X2, X3, X4. When we do a spline fit in this case, we are actually fitting two sets of splines: theta versus X1 and theta versus Y1. And like wise for (X2,Y2), (X3,Y3), etc. A simple-minded person like me can get confused about when x means theta and when x means the x-coordinate of each data point.
Another thing that is potentially confusing is that the angle of each data point on the Cartesian plane has nothing to do with the value of the associated phase angle. For example, the first elements of (X1,Y1), (X2,Y2), (X3,Y3), and (X4,Y4) are all associated with phase angle zero in your model. But the angles to those four points, on the Cartesian plane, are -0.47, -1.64, -2.28, and -2.49 radians, repectvely. You can verify that by computing atan2(X1,Y1), etc., in the Excel spreadsheet.
Back to your orignal question: why does the fitted spline look un-smooth , especially for the X3, Y3 data? I recommended that you add several "duplicate" points before the first and after the last, so the splines wrap around consistently. It turns out that does not totally fix th problem in this case. The reason is that spline fits sometimes overshoot and undershoot between points, and that is hapening in this case, for the phase-versus-X3 spline fit. See plots that are roduced by the attached version of your code, which I have modified a bit. The code runs on the attached spreadsheet, in which the first seven and last seven points are the same. Thus it has 31 rows of numbers, but only 24 of the rows are unique.
Plots 1 and 2 are made by the code. Plot 1 is like your plot. Plot 2 shows the spline fit of phase versus X and of phase versus Y.
Plots 3 and 4 are zoomed in versions of 1 and 2 showing the section of the X3,Y3 data where the spline fit looks un-smooth. Plot shows that the spline fit is smooth through the points, but it has undershoot and overshoot. Plot 4 show that this is due to the spline fit of phase versus X3 data.
I think a better approach would be to interpolate with interp1(), and use the 'makima' or 'pchip' option, because they don't undershoot or overshoot. See the help for interp1().

댓글을 달려면 로그인하십시오.

Community Treasure Hunt

Find the treasures in MATLAB Central and discover how the community can help you!

Start Hunting!

Translated by