MATLAB Answers

How to apply a fit on this shape?

10 views (last 30 days)
Niklas Kurz
Niklas Kurz on 15 Jul 2020
Edited: Niklas Kurz on 19 Jul 2020
Since my data doens't really count as a function cftool is not an option. I wonder whether this is still faisable to fit:
(I tried several elliptical fits but they don't seem to work. May be because I dealt with it a little clumsily)

  0 Comments

Sign in to comment.

Accepted Answer

John D'Errico
John D'Errico on 15 Jul 2020
Edited: John D'Errico on 15 Jul 2020
What does it mean to "fit" it to you? What will you do with that fit? The data is not a function, just a set of points, that happen to represent a closed curve in the plane.
A simple solution might be to convert to polar coordinates.
help cart2pol
Now, fit the function r(theta). Be careful, since r is a periodic function. It should have the property that r(0) == r(2*pi). So a polynomial fit would be a bad idea. Instead, perhaps fit this as a sum.
r(theta) = a0 + a1*sin(theta) + b1*cos(theta) + a2*sin(2*theta) + b2*cos(2*theta) + ...
Essentially, I am describing what is essentially a Fourier series there (that is, not an fft, but a Fourier seies.)
As it is constructed, the function ALWAYS has the property of being periodic. And, while you could, in theory, trivially use the curve fitting toolbox for this task, there is no reason to do so. All of the parameters can be estimated using one call to backslash.
x = x(:);
y = y(:);
mux = mean(x);
muy = mean(y);
[Theta,R] = cart2pol(x- mux,y-muy);
[Theta,ind] = sort(Theta);
R = R(ind);
plot(Theta,R,'o')
n = numel(x);
M = 6;
A = [ones(n,1),sin(Theta*(1:M)),cos(Theta*(1:M))];
coeffs = A\R;
Rhat = A*coeffs;
hold on
plot(Theta,Rhat)
Larger values of M will yield a better approximation. A perfect circle would have m == 0. Because your curve has an oval shape, you need to have a few extra terms, but not many. M == 6 will require the estimation of 13 coefficients for the series.
I lack your data, so I generated some random crap data as an example, using ginput.
plot(x,y,'o')
Nothing exciting. Vaguely ovoid. I'm actually surprised it was that good.
R as a function of Theta seems reasonable.
And now, we can reconstruct the original curve.
[xhat,yhat] = pol2cart(Theta,Rhat);
plot(x,y,'bo',xhat + mux,yhat + muy,'-xr')
Surprisingly good for a simple, moderately low order Fourier series approximation, given the data came from my shaky eye/mouse coordinated hand.

  6 Comments

Show 3 older comments
Niklas Kurz
Niklas Kurz on 17 Jul 2020
If I may ask additionally: How to determine the area beneath the curve? Pretending the ellipse would be above the x axis like in the data I attached as well. Since the values are repeating trapz(X,Y) doesn't seem right.
John D'Errico
John D'Errico on 19 Jul 2020
Why does this question seem to be veering in random directions, with me trying to chase a moving target?
If the goal is to find the area betwen the lower lobe of that closed curve, and the x-axis, there is a problem. You need to find the right and left hand end points to integrate between. A problem is the data is relatively noisy, but also, the noise is not really true noise.
And given this data, the code I wrote needs quite a few terms to perform even reasonably well.
FSfit(x,y,100,1);
While the fit does not look too bad when viewed as an entire plot, we can see it shows serious oscillations in places. Another option is to use a smoothing spline, but that looks just as bad.
I'd suggest a big part of the problem is the data is heavily quantized in the x variable. For example, there are 9287 data points in your curve. However, x is clearly not any sort of continuous variable.
numel(unique(x))
ans =
125
numel(unique(y))
ans =
9287
One thing that I often do whenever someone asks a question about how to model some data, the first thing I often ask is what do you want to get out of it? I should have pushed for that here. If the goal is merely to find the area between the lower lobe of that curve and the x axis, this very much impacts how I would want to approach the problem.
The scheme I showed you how to implement will be pretty robust to approximate the area contained within the closed curve. That is because it is insensitive to high frequency oscillations in the estimated curve, because those oscillations have high frequency, but small amplitude.
[~,~,~,~,~,~,area] = FSfit(x,y,10,0)
area =
13979.5620217219
>> [~,~,~,~,~,~,area] = FSfit(x,y,20,0)
area =
13726.1735377685
>> [~,~,~,~,~,~,area] = FSfit(x,y,30,0)
area =
13708.1851809597
>> [~,~,~,~,~,~,area] = FSfit(x,y,100,0)
area =
13707.1667535602
Given the noise, the prediction seems to be fairly stable as long as you use 20 terms in the model, or so.
But to find the area under that curve will be more difficult if we are working with a model of the form I built. So IF my goal was to compute the area blow the curve and above the x axis, I would do it like this:
xmin = min(x);
xmax = max(x);
xadj = [x;xmin;xmax];
yadj = [y;0;0];
E1 = convhull(x,y);
E2 = convhull(xadj,yadj);
polyarea(xadj(E2),yadj(E2)) - polyarea(x(E1),y(E1))
ans =
4811.42632861809
So we approximate the area under the curve as the difference of two areas, thus the area of the region in purple here:
plot(polyshape(xadj(E2),yadj(E2)))
hold on
plot(polyshape(x(E1),y(E1)))
Niklas Kurz
Niklas Kurz on 19 Jul 2020
Oh.My.Gosh, you're on of the best! I'm like 100000% satisfied now. Excuse me for my randomness, sometimes even I don't know what I'm searching for.

Sign in to comment.

More Answers (0)

Community Treasure Hunt

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

Start Hunting!