15 views (last 30 days)

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)

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.

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

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

Start Hunting!
## 0 Comments

Sign in to comment.