Discover MakerZone

MATLAB and Simulink resources for Arduino, LEGO, and Raspberry Pi

Learn more

Discover what MATLAB® can do for your career.

Opportunities for recent engineering grads.

Apply Today

Thread Subject:
fitting points to a sin function

Subject: fitting points to a sin function

From: Rafael Herrejon

Date: 25 Feb, 2008 09:44:06

Message: 1 of 9

Hello, I dont know what I am doing wrong but I cant fit a
cloud of points to a sin function.

Im using around 100 points but for simplicity ill write
this question with just 6 points to fix

I have used fminsearch as follow

xdata=(1.3485:0.0145:1.4065)';
ydata=[ 1.0472; 0.5411; 0; 0.7752; 1.0472];
[estimates2, model2] = fitcurvedemo2(xdata,ydata)

where

function [estimates2, model2] = fitcurvedemo(xdata, ydata)
start_point2 = [0.6 1 0.6];
model2 = @expfun2;
estimates2 = fminsearch(model2, start_point2);

function [sse, Fitted] = expfun2(params)
        A = params(1);
        B = params(2);
        C = params(2);
        Fitted = A .* sin(B * xdata)+C;
        ErrorVector = Fitted - ydata;
        sse = sum(ErrorVector .^ 2);
    end
end

Even though it is visible that a sin function could fit, im
not able to get the parameters correctly

I have tried also with the fit function of curve fitting
toolbox

fopts=fitoptions('Method','NonlinearLeastSquares',...
                'Lower',[0,0],...
                'Upper',[Inf,max(xdata)],...
                'Startpoint',[.6 1 .6],...
                'Algorithm', 'Levenberg-Marquardt');
ftype=fittype('a*sin(b*x)+c','options',fopts);
sinufit=fit(xdata,ydata,ftype);

as well as the cftool obtaining same wrong results. I
REALLY need your help. I have been trying this for weeks
and I dont know what to do anymore.

btw, if u are interested to try with the whole data, here
it is.

xdata=(0:0.0145:1.4065)';
ydata=[ 0.5411; 0; 0; 0.9626;
0.5411; 0.5054; 0; 1.1279; 0.5054;
            0.4759; 0.7227; 0.8957; 0.4759;
0.4759; 0.8957; 0.6797; 0.4510; 0.4759;
            1.1102; 0; 0.4510; 0.6797;
0.9273; 0.4297; 0.4297; 0.9273; 0.6435;
            0.4111; 0.4297; 0.9273; 0.5857;
0.5621; 0.6126; 0.8810; 0.4111; 0.5621;
            0.8810; 0.7227; 0.5411; 0;
0.9480; 0.4111; 0.5411; 0.5857; 0.9480;
            0.5411; 0.5411; 0.7227; 0.7227;
0.5411; 0.3948; 0.9480; 0.5621; 0.6435;
                 0; 0.8810; 0; 0.5411;
0.7227; 0.8810; 0.5621; 0.5621; 0.8810;
            0.7565; 0.4111; 0; 1.0989;
0.6126; 0.4111; 0.4297; 1.1593; 0.4510;
            0.4297; 0.6435; 1.0472;
0; 0; 0.6797; 0.9818; 0; 0;
            0.8957; 0.7227; 0; 0.4759;
1.0472; 0.7227; 0; 0.5054; 1.1279;
            0.7227; 0; 0.7227; 1.0472;
0.5411; 0; 0.7752; 1.0472];

really appreciate your help.

Best Regards. Rafael

Subject: fitting points to a sin function

From: John D'Errico

Date: 25 Feb, 2008 11:10:22

Message: 2 of 9

"Rafael Herrejon" <rafael.erasethis@ic.is.tohoku.ac.jp> wrote in message
<fpu2l6$5ts$1@fred.mathworks.com>...
> Hello, I dont know what I am doing wrong but I cant fit a
> cloud of points to a sin function.

(snip)

> function [estimates2, model2] = fitcurvedemo(xdata, ydata)
> start_point2 = [0.6 1 0.6];
> model2 = @expfun2;
> estimates2 = fminsearch(model2, start_point2);
>
> function [sse, Fitted] = expfun2(params)
> A = params(1);
> B = params(2);
> C = params(2);
> Fitted = A .* sin(B * xdata)+C;
> ErrorVector = Fitted - ydata;
> sse = sum(ErrorVector .^ 2);
> end
> end
>
> Even though it is visible that a sin function could fit, im
> not able to get the parameters correctly

(snip)

> as well as the cftool obtaining same wrong results. I
> REALLY need your help. I have been trying this for weeks
> and I dont know what to do anymore.
>
> btw, if u are interested to try with the whole data, here
> it is.
>
> xdata=(0:0.0145:1.4065)';
> ydata=[ 0.5411; 0; 0; 0.9626;

...

Why do you think that this arbitrary cloud
of points has a sine wave as a model?

Why do you think that this particular sine
model is correct?

What do the exact zeros in your data result
from? If they are missing values, then don't
you think that the modeling will have
difficulties dealing with missing data that
are replaced with exact zeros?

John

Subject: fitting points to a sin function

From: Rafael Herrejon

Date: 25 Feb, 2008 11:54:06

Message: 3 of 9

"John D'Errico" <woodchips@rochester.rr.com> wrote in
message <fpu7mu$j82$1@fred.mathworks.com>...
> "Rafael Herrejon" <rafael.erasethis@ic.is.tohoku.ac.jp>
wrote in message
> <fpu2l6$5ts$1@fred.mathworks.com>...
> > Hello, I dont know what I am doing wrong but I cant fit
a
> > cloud of points to a sin function.
>
> (snip)
>
> > function [estimates2, model2] = fitcurvedemo(xdata,
ydata)
> > start_point2 = [0.6 1 0.6];
> > model2 = @expfun2;
> > estimates2 = fminsearch(model2, start_point2);
> >
> > function [sse, Fitted] = expfun2(params)
> > A = params(1);
> > B = params(2);
> > C = params(2);
> > Fitted = A .* sin(B * xdata)+C;
> > ErrorVector = Fitted - ydata;
> > sse = sum(ErrorVector .^ 2);
> > end
> > end
> >
> > Even though it is visible that a sin function could
fit, im
> > not able to get the parameters correctly
>
> (snip)
>
> > as well as the cftool obtaining same wrong results. I
> > REALLY need your help. I have been trying this for
weeks
> > and I dont know what to do anymore.
> >
> > btw, if u are interested to try with the whole data,
here
> > it is.
> >
> > xdata=(0:0.0145:1.4065)';
> > ydata=[ 0.5411; 0; 0; 0.9626;
>
> ...
>
> Why do you think that this arbitrary cloud
> of points has a sine wave as a model?

This "arbitrary" cloud of points is the calculated angle of
a coin rotating, observed by a camera.

>
> Why do you think that this particular sine
> model is correct?
the coin is rotating on its z-axis, the angle would go from
0 to 90 degrees..

>
> What do the exact zeros in your data result
> from?
they are not missing values. The values are calculated with
theta=cos(b/a) where b is the minor axis of the ellipse and
a is the major axis of the ellipse, when the minor axis is
equal to the major axis, the angle of the coin is 0.

>If they are missing values, then don't
> you think that the modeling will have
> difficulties dealing with missing data that
> are replaced with exact zeros?
>
if you plot the last 20 points u can see that an abs(sin)
could fit. Again is least squares fitting, even the fitting
is not perfect, should be able to get some results closer
to the ones im getting

Subject: fitting points to a sin function

From: Roger Stafford

Date: 26 Feb, 2008 01:30:21

Message: 4 of 9

"Rafael Herrejon" <rafael.erasethis@ic.is.tohoku.ac.jp> wrote in message
<fpua8u$aqi$1@fred.mathworks.com>...
> "John D'Errico" <woodchips@rochester.rr.com> wrote in
> message <fpu7mu$j82$1@fred.mathworks.com>...
> > Why do you think that this arbitrary cloud
> > of points has a sine wave as a model?
>
> This "arbitrary" cloud of points is the calculated angle of
> a coin rotating, observed by a camera.
>
> > Why do you think that this particular sine
> > model is correct?

> the coin is rotating on its z-axis, the angle would go from
> 0 to 90 degrees..
>
> >
> > What do the exact zeros in your data result
> > from?

> they are not missing values. The values are calculated with
> theta=cos(b/a) where b is the minor axis of the ellipse and
> a is the major axis of the ellipse, when the minor axis is
> equal to the major axis, the angle of the coin is 0.
>
> >If they are missing values, then don't
> > you think that the modeling will have
> > difficulties dealing with missing data that
> > are replaced with exact zeros?
> >
> if you plot the last 20 points u can see that an abs(sin)
> could fit. Again is least squares fitting, even the fitting
> is not perfect, should be able to get some results closer
> to the ones im getting
---------
  If this data represents the angle of a rotating coin, do you actually mean
theta = arccos(b/a) instead of cos(b/a)?

  If xdata represents elapsed time and ydata is this theta, one would expect
that the total angle of revolution of the coin would be linear in time because
of conservation of angular momentum, provided you allowed the angle to go
beyond the artificial limits of 0 to pi/2 (and provided the rotation is about a
line in the coin instead of a coin wobbling about a skew line.) Such a straight
line is not a good candidate for fitting with a sine curve. By constraining it to
the range 0 to pi/2 it would behave like a ball bouncing between two parallel
y-limit walls following a straight line between limit. This is still not a very apt
candidate for fitting with a sine curve.

  If the major axis, a, remains fixed, I would think you would want to do a sine
fit as a function of time on the value, b/a, itself except that you would want
to reverse its sign each time it passed through zero.

  Finally, as a general comment, it seems to me that you have not allowed for
the very important possible phase shift in your three parameters, A, B, and C.
You ought to be using four parameters,

 A*sin(B*x+C) + D,

in my opinion, to allow for the full generality of a sine fit.

Roger Stafford

Subject: fitting points to a sin function

From: Rafael Herrejon

Date: 26 Feb, 2008 05:25:04

Message: 5 of 9

> If this data represents the angle of a rotating coin,
do you actually mean
> theta = arccos(b/a) instead of cos(b/a)?
>
Roger. thank you for your comments, you are right, i meant
theta = arccos(b/a). Ydata is the "angle" observed by a
camera in that instant t, and this angle is calculated with
the formula above, (Arccos(minor/major)), the coin is
rotating on its major axis, so the 0 < minor axis < major
axis. this result in theta being always positive and
0<theta<pi/2. So i have discrete values of theta(t) that
should be able to fit an |A*sin(B*x+C)+D|, or at least
thats what i thought, a i wrong?
If i make it b/a as a function of time, the values still
have a sinusuidal behaviour.

What should i do?
Best Regards, Rafael

Subject: fitting points to a sin function

From: Roger Stafford

Date: 26 Feb, 2008 08:05:03

Message: 6 of 9

"Rafael Herrejon" <rafael.erasethis@ic.is.tohoku.ac.jp> wrote in message
<fq07rg$jf3$1@fred.mathworks.com>...
> Roger. thank you for your comments, you are right, i meant
> theta = arccos(b/a). Ydata is the "angle" observed by a
> camera in that instant t, and this angle is calculated with
> the formula above, (Arccos(minor/major)), the coin is
> rotating on its major axis, so the 0 < minor axis < major
> axis. this result in theta being always positive and
> 0<theta<pi/2. So i have discrete values of theta(t) that
> should be able to fit an |A*sin(B*x+C)+D|, or at least
> thats what i thought, a i wrong?
> If i make it b/a as a function of time, the values still
> have a sinusuidal behaviour.
>
> What should i do?
> Best Regards, Rafael
--------
  As I said before, Rafael, I claim that if you plot the coin's angle as you have
defined it with the principal value of theta = arccos(b/a), then you will get a
saw-toothed curve consisting of straight-line segments that "bounce" back
and forth between 0 and pi/2. This is because the true angle of rotation in
terms of total rotation is not limited to the range from 0 to pi/2 nor 0 to 2*pi,
and as such would be a linear function of time which would give a straight
line graph. You would then be squeezing this angle into the narrow range, 0
to pi/2, to give the linear saw-toothed segments.

  Try the following to see what I mean. The yellow dotted line would give the
total rotation angle as a function of time through a couple of complete
rotations, and the red line would be the saw-toothed part. The blue curve is
b/a before changing it to an always positive quantity.

 time = linspace(0,1,100); % One second duration
 omega = 13; % Constant angular velocity in radians per sec. (~2 rps)
 totangle = omega*time; % Total angle of rotation in radians over 1 sec.
 a = 2.3; % Diameter of quarter in cm.
 b1 = a*cos(totangle); % "Signed" minor axis
 b2 = abs(b1); % Minor axis of ellipse, made always positive
 theta = acos(b2/a); % Take the principal value
 plot(time,theta,'r-',time,totangle,'y.',time,b1/a,'b-')

As you can see, neither the yellow nor red curves could be fitted easily into a
good sine curve fit. Only the blue curve has such a capability, and it is a
perfect sinusoidal-type curve.

  You haven't yet made it clear what it is you expect to accomplish with your
analysis. What kind of relationship do you expect to arise out of the
sinusoidal fitting that you couldn't get by simply measuring the rotation
period accurately?

Roger Stafford

Subject: fitting points to a sin function

From: Izhak Bucher

Date: 26 Feb, 2008 08:56:01

Message: 7 of 9

Rafael.

Here's some code to help you (below).
A short explanation. The problem can be broken into linear
least squares and nonlinear. Fitting a sine with a phase
can be accomplished via modeling the measurement as:
   y = A*sin(w*t)+B*cos(w*t);
this function is linear in A and B and nonlinear in w.
 Given w, you can solve for A and B.
As for your data, it comtains several harmonics (I have
checked via fft).
  here's the code
==============================
%fit_sin_example.m
% I. Bucher
%% measureements
xdata=(0:0.0145:1.4065)';
ydata=[ 0.5411; 0; 0; 0.9626;
0.5411; 0.5054; 0; 1.1279; 0.5054;
            0.4759; 0.7227; 0.8957; 0.4759;
0.4759; 0.8957; 0.6797; 0.4510; 0.4759;
            1.1102; 0; 0.4510; 0.6797;
0.9273; 0.4297; 0.4297; 0.9273; 0.6435;
            0.4111; 0.4297; 0.9273; 0.5857;
0.5621; 0.6126; 0.8810; 0.4111; 0.5621;
            0.8810; 0.7227; 0.5411; 0;
0.9480; 0.4111; 0.5411; 0.5857; 0.9480;
            0.5411; 0.5411; 0.7227; 0.7227;
0.5411; 0.3948; 0.9480; 0.5621; 0.6435;
                 0; 0.8810; 0; 0.5411;
0.7227; 0.8810; 0.5621; 0.5621; 0.8810;
            0.7565; 0.4111; 0; 1.0989;
0.6126; 0.4111; 0.4297; 1.1593; 0.4510;
            0.4297; 0.6435; 1.0472;
0; 0; 0.6797; 0.9818; 0; 0;
            0.8957; 0.7227; 0; 0.4759;
1.0472; 0.7227; 0; 0.5054; 1.1279;
            0.7227; 0; 0.7227; 1.0472;
0.5411; 0; 0.7752; 1.0472];
%
 
N=length(xdata); % points
 
%------> rough estimate of dominant frequency for initial
guess
 Y=fft(ydata-mean(ydata));
  imax = max(abs(Y(1:fix(N/2))));
  dx=1/(xdata(2)-xdata(1))/N;
  f0 = dx*(imax-1); % highest peak
   
   
 %--------> Now prepare for fit
 
 figure
op=optimset('display','none');
op.TolX=1e-18;
w0=2*pi*f0;
Nharq=2; % no. of harmonic to fit

%>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>
% >>> Call the optimisation process
% <<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<
  [w,op,fcost,J]=lsqnonlin('fitw',w0,[],
[],op,xdata,ydata,Nharq);

 [dy,csq,yfit,a]=fitw(w,xdata,ydata,Nharq);
 title( sprintf('estimated frequency %f \n',w/2/pi))
 fprintf('\results\n')
 for q=1:Nharq
     fprintf(' Harmonic # %d -> Estimated Amplitude %f
\n',q,norm(a(2*q-1:2*q)))
 end
 
 
 legend('measured','fit','error: y-yfit')
=================================================
%the function fitw.m
function [dy,csq,yfit,a]=fitw(w,x,y,Nharq)
 
nt=length(y);
argq=x(:)*w; % generate the function to fit with main true
excitation frequency

 
csq=[];
  for qi=1:Nharq,
    csq=[csq cos(qi*argq) sin(qi*argq)];
 end
 csq=[csq ones(nt,1)];
 
a=csq\y;
yfit=csq*a;
 dy=y-yfit;
h=plot([y yfit y-yfit]);, drawnow
%set(h(2),'marker','-');
set(h(1),'marker','.','markersize',5);



"Rafael Herrejon" <rafael.erasethis@ic.is.tohoku.ac.jp>
wrote in message <fpu2l6$5ts$1@fred.mathworks.com>...
> Hello, I dont know what I am doing wrong but I cant fit a
> cloud of points to a sin function.
>
> Im using around 100 points but for simplicity ill write
> this question with just 6 points to fix
>
> I have used fminsearch as follow
>
> xdata=(1.3485:0.0145:1.4065)';
> ydata=[ 1.0472; 0.5411; 0; 0.7752; 1.0472];
> [estimates2, model2] = fitcurvedemo2(xdata,ydata)
>
> where
>
> function [estimates2, model2] = fitcurvedemo(xdata,
ydata)
> start_point2 = [0.6 1 0.6];
> model2 = @expfun2;
> estimates2 = fminsearch(model2, start_point2);
>
> function [sse, Fitted] = expfun2(params)
> A = params(1);
> B = params(2);
> C = params(2);
> Fitted = A .* sin(B * xdata)+C;
> ErrorVector = Fitted - ydata;
> sse = sum(ErrorVector .^ 2);
> end
> end
>
> Even though it is visible that a sin function could fit,
im
> not able to get the parameters correctly
>
> I have tried also with the fit function of curve fitting
> toolbox
>
> fopts=fitoptions('Method','NonlinearLeastSquares',...
> 'Lower',[0,0],...
> 'Upper',[Inf,max(xdata)],...
> 'Startpoint',[.6 1 .6],...
> 'Algorithm', 'Levenberg-Marquardt');
> ftype=fittype('a*sin(b*x)+c','options',fopts);
> sinufit=fit(xdata,ydata,ftype);
>
> as well as the cftool obtaining same wrong results. I
> REALLY need your help. I have been trying this for weeks
> and I dont know what to do anymore.
>
> btw, if u are interested to try with the whole data, here
> it is.
>
> xdata=(0:0.0145:1.4065)';
> ydata=[ 0.5411; 0; 0; 0.9626;
> 0.5411; 0.5054; 0; 1.1279; 0.5054;
> 0.4759; 0.7227; 0.8957; 0.4759;
> 0.4759; 0.8957; 0.6797; 0.4510; 0.4759;
> 1.1102; 0; 0.4510; 0.6797;
> 0.9273; 0.4297; 0.4297; 0.9273; 0.6435;
> 0.4111; 0.4297; 0.9273; 0.5857;
> 0.5621; 0.6126; 0.8810; 0.4111; 0.5621;
> 0.8810; 0.7227; 0.5411; 0;
> 0.9480; 0.4111; 0.5411; 0.5857; 0.9480;
> 0.5411; 0.5411; 0.7227; 0.7227;
> 0.5411; 0.3948; 0.9480; 0.5621; 0.6435;
> 0; 0.8810; 0; 0.5411;
> 0.7227; 0.8810; 0.5621; 0.5621; 0.8810;
> 0.7565; 0.4111; 0; 1.0989;
> 0.6126; 0.4111; 0.4297; 1.1593; 0.4510;
> 0.4297; 0.6435; 1.0472;
> 0; 0; 0.6797; 0.9818; 0; 0;
> 0.8957; 0.7227; 0; 0.4759;
> 1.0472; 0.7227; 0; 0.5054; 1.1279;
> 0.7227; 0; 0.7227; 1.0472;
> 0.5411; 0; 0.7752; 1.0472];
>
> really appreciate your help.
>
> Best Regards. Rafael
>

Subject: fitting points to a sin function

From: Rafael Herrejon

Date: 27 Feb, 2008 22:15:04

Message: 8 of 9

> As I said before, Rafael, I claim that if you plot the
coin's angle as you have
> defined it with the principal value of theta = arccos
(b/a), then you will get a
> saw-toothed curve consisting of straight-line segments
that "bounce" back
> and forth between 0 and pi/2. This is because the true
angle of rotation in
> terms of total rotation is not limited to the range from
0 to pi/2 nor 0 to 2*pi,
> and as such would be a linear function of time which
would give a straight
> line graph. You would then be squeezing this angle into
the narrow range, 0
> to pi/2, to give the linear saw-toothed segments.
>
> Try the following to see what I mean. The yellow
dotted line would give the
> total rotation angle as a function of time through a
couple of complete
> rotations, and the red line would be the saw-toothed
part. The blue curve is
> b/a before changing it to an always positive quantity.
>
> time = linspace(0,1,100); % One second duration
> omega = 13; % Constant angular velocity in radians per
sec. (~2 rps)
> totangle = omega*time; % Total angle of rotation in
radians over 1 sec.
> a = 2.3; % Diameter of quarter in cm.
> b1 = a*cos(totangle); % "Signed" minor axis
> b2 = abs(b1); % Minor axis of ellipse, made always
positive
> theta = acos(b2/a); % Take the principal value
> plot(time,theta,'r-',time,totangle,'y.',time,b1/a,'b-')
>
> As you can see, neither the yellow nor red curves could
be fitted easily into a
> good sine curve fit. Only the blue curve has such a
capability, and it is a
> perfect sinusoidal-type curve.
>
> You haven't yet made it clear what it is you expect to
accomplish with your
> analysis. What kind of relationship do you expect to
arise out of the
> sinusoidal fitting that you couldn't get by simply
measuring the rotation
> period accurately?

Hello Roger, thanks for your answer, i could understand
your point with the example but let me try to explain what
i want to do a little more detailed.

Using a camera, i can see a rotating coin in the distance.
Then coin gets closer to the camera and then gets farther
again. The image being seen is a serie of ellipses and
circles (depending on the rotating angle). fitting an
ellipse to those images, i can get its major axis( or axis
of rotation A) and the minor axis(B). Seen by the
camera, '0'<B<A, (in reality this is impossible because the
coin width<>0 ).
When the camera sees the coin as a circle, the coin's
relative rotation angle is 0deg(or 180deg) and when it is
seen as a 'line' the coins angle would be 90deg.

This is the axis im able to get from the image processing

a=
[7;7;7;7;7;8;7;7;8;9;8;8;9;9;8;9;10;9;9;9;10;9;10;11;11;10;1
0;12;11;10;12;13;11;11;12;13;11;12;14;12;12;12;14;12;12;14;1
4;12;12;14;13;12;13;15;12;11;12;14;12;11;13;13;11;11;12;11;1
1;11;12;11;10;10;11;10;10;10;10;9;9;9;9;8;8;9;9;8;8;8;8;7;8;
8;8;6;7;7;7;6];

b=
[6;7;7;4;6;7;7;3;7;8;6;5;8;8;5;7;9;8;4;9;9;7;6;10;10;6;8;11;
10;6;10;11;9;7;11;11;7;9;12;12;7;11;12;10;7;12;12;9;9;12;12;
7;11;12;12;7;12;12;9;7;11;11;7;8;11;11;5;9;11;10;4;9;10;8;5;
10;10;7;5;9;9;5;6;9;8;4;6;8;7;3;6;8;6;3;6;7;5;3]


But as you can see the data is not so clean, due to
ilumination or other factors(distance to the camera, blur
in the image, etc).
anyways, if i calculate the relative angle as

theta = acos(b2/a);

i get

xdata=0:0.0145:1.4065;
ydata=[ 0.5411; 0; 0; 0.9626; 0.5411; 0.5054; 0; 1.1279;
0.5054; 0.4759; 0.7227; 0.8957; 0.4759; 0.4759; 0.8957;
0.6797; 0.4510; 0.4759; 1.1102; 0; 0.4510; 0.6797; 0.9273;
0.4297; 0.4297; 0.9273; 0.6435; 0.4111; 0.4297; 0.9273;
0.5857; 0.5621; 0.6126; 0.8810; 0.4111; 0.5621; 0.8810;
0.7227; 0.5411; 0; 0.9480; 0.4111; 0.5411; 0.5857; 0.9480;
0.5411; 0.5411; 0.7227; 0.7227; 0.5411; 0.3948; 0.9480;
0.5621; 0.6435; 0; 0.8810; 0; 0.5411; 0.7227; 0.8810;
0.5621; 0.5621; 0.8810; 0.7565; 0.4111; 0; 1.0989; 0.6126;
0.4111; 0.4297; 1.1593; 0.4510; 0.4297; 0.6435; 1.0472; 0;
0; 0.6797; 0.9818; 0; 0; 0.8957; 0.7227; 0; 0.4759; 1.0472;
0.7227; 0; 0.5054; 1.1279; 0.7227; 0; 0.7227; 1.0472;
0.5411; 0; 0.7752; 1.0472];


and then i was trying to fix a sinusoidal wave to this
values so i can predict the relative rotation angle further
in time.

The data is full of noise, but a simple sinusoidal wave
(A*sin(B*t+C)) can be seen using only the last 18 values of
data

xdata=1.16:0.0145:1.4065;
ydata=[ 0; 0.8957; 0.7227; 0; 0.4759; 1.0472; 0.7227; 0;
0.5054; 1.1279; 0.7227; 0; 0.7227; 1.0472; 0.5411; 0;
0.7752; 1.0472];
conttime=xdata(1):0.0005:xdata(size(xdata));
figure
hold on
plot(xdata,ydata,'*');
plot(conttime,abs(1.16*sin((17.5239*pi*conttime)+2.8)),'r-
');
hold off

The parameters A,B and C where found by trial and error and
i havent calculated the error, but i would like to know how
to get the best parameters with a fitting algorithm that
minimizes the error. Measuring the rotation the rotation
period acurately its impossible with this noisy data.

What do you recommend me to do?

thank you in advance

Subject: fitting points to a sin function

From: Rafael Herrejon

Date: 27 Feb, 2008 22:24:01

Message: 9 of 9

Izhak, thank you so much for your hard work, it fits
perfectly to my data but i guess i couldn express what i
need to get correctly, you can read it a little more
detailed here

http://www.mathworks.com/matlabcentral/newsreader/view_threa
d/164459#417830

maybe you could give me an idea of what should i do.

Best Regards, Rafael

> Rafael.
>
> Here's some code to help you (below).
> A short explanation. The problem can be broken into
linear
> least squares and nonlinear. Fitting a sine with a phase
> can be accomplished via modeling the measurement as:
> y = A*sin(w*t)+B*cos(w*t);
> this function is linear in A and B and nonlinear in w.
> Given w, you can solve for A and B.
> As for your data, it comtains several harmonics (I have
> checked via fft).
> here's the code
> ==============================
> %fit_sin_example.m
> % I. Bucher
> %% measureements
> xdata=(0:0.0145:1.4065)';
> ydata=[ 0.5411; 0; 0; 0.9626;
> 0.5411; 0.5054; 0; 1.1279; 0.5054;
> 0.4759; 0.7227; 0.8957; 0.4759;
> 0.4759; 0.8957; 0.6797; 0.4510; 0.4759;
> 1.1102; 0; 0.4510; 0.6797;
> 0.9273; 0.4297; 0.4297; 0.9273; 0.6435;
> 0.4111; 0.4297; 0.9273; 0.5857;
> 0.5621; 0.6126; 0.8810; 0.4111; 0.5621;
> 0.8810; 0.7227; 0.5411; 0;
> 0.9480; 0.4111; 0.5411; 0.5857; 0.9480;
> 0.5411; 0.5411; 0.7227; 0.7227;
> 0.5411; 0.3948; 0.9480; 0.5621; 0.6435;
> 0; 0.8810; 0; 0.5411;
> 0.7227; 0.8810; 0.5621; 0.5621; 0.8810;
> 0.7565; 0.4111; 0; 1.0989;
> 0.6126; 0.4111; 0.4297; 1.1593; 0.4510;
> 0.4297; 0.6435; 1.0472;
> 0; 0; 0.6797; 0.9818; 0; 0;
> 0.8957; 0.7227; 0; 0.4759;
> 1.0472; 0.7227; 0; 0.5054; 1.1279;
> 0.7227; 0; 0.7227; 1.0472;
> 0.5411; 0; 0.7752; 1.0472];
> %
>
> N=length(xdata); % points
>
> %------> rough estimate of dominant frequency for initial
> guess
> Y=fft(ydata-mean(ydata));
> imax = max(abs(Y(1:fix(N/2))));
> dx=1/(xdata(2)-xdata(1))/N;
> f0 = dx*(imax-1); % highest peak
>
>
> %--------> Now prepare for fit
>
> figure
> op=optimset('display','none');
> op.TolX=1e-18;
> w0=2*pi*f0;
> Nharq=2; % no. of harmonic to fit
>
> %>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>
> % >>> Call the optimisation process
> % <<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<
> [w,op,fcost,J]=lsqnonlin('fitw',w0,[],
> [],op,xdata,ydata,Nharq);
>
> [dy,csq,yfit,a]=fitw(w,xdata,ydata,Nharq);
> title( sprintf('estimated frequency %f \n',w/2/pi))
> fprintf('\results\n')
> for q=1:Nharq
> fprintf(' Harmonic # %d -> Estimated Amplitude %f
> \n',q,norm(a(2*q-1:2*q)))
> end
>
>
> legend('measured','fit','error: y-yfit')
> =================================================
> %the function fitw.m
> function [dy,csq,yfit,a]=fitw(w,x,y,Nharq)
>
> nt=length(y);
> argq=x(:)*w; % generate the function to fit with main
true
> excitation frequency
>
>
> csq=[];
> for qi=1:Nharq,
> csq=[csq cos(qi*argq) sin(qi*argq)];
> end
> csq=[csq ones(nt,1)];
>
> a=csq\y;
> yfit=csq*a;
> dy=y-yfit;
> h=plot([y yfit y-yfit]);, drawnow
> %set(h(2),'marker','-');
> set(h(1),'marker','.','markersize',5);
>
>
>
> "Rafael Herrejon" <rafael.erasethis@ic.is.tohoku.ac.jp>
> wrote in message <fpu2l6$5ts$1@fred.mathworks.com>...
> > Hello, I dont know what I am doing wrong but I cant fit
a
> > cloud of points to a sin function.
> >
> > Im using around 100 points but for simplicity ill write
> > this question with just 6 points to fix
> >
> > I have used fminsearch as follow
> >
> > xdata=(1.3485:0.0145:1.4065)';
> > ydata=[ 1.0472; 0.5411; 0; 0.7752; 1.0472];
> > [estimates2, model2] = fitcurvedemo2(xdata,ydata)
> >
> > where
> >
> > function [estimates2, model2] = fitcurvedemo(xdata,
> ydata)
> > start_point2 = [0.6 1 0.6];
> > model2 = @expfun2;
> > estimates2 = fminsearch(model2, start_point2);
> >
> > function [sse, Fitted] = expfun2(params)
> > A = params(1);
> > B = params(2);
> > C = params(2);
> > Fitted = A .* sin(B * xdata)+C;
> > ErrorVector = Fitted - ydata;
> > sse = sum(ErrorVector .^ 2);
> > end
> > end
> >
> > Even though it is visible that a sin function could
fit,
> im
> > not able to get the parameters correctly
> >
> > I have tried also with the fit function of curve
fitting
> > toolbox
> >
> > fopts=fitoptions('Method','NonlinearLeastSquares',...
> > 'Lower',[0,0],...
> > 'Upper',[Inf,max(xdata)],...
> > 'Startpoint',[.6 1 .6],...
> > 'Algorithm', 'Levenberg-Marquardt');
> > ftype=fittype('a*sin(b*x)+c','options',fopts);
> > sinufit=fit(xdata,ydata,ftype);
> >
> > as well as the cftool obtaining same wrong results. I
> > REALLY need your help. I have been trying this for
weeks
> > and I dont know what to do anymore.
> >
> > btw, if u are interested to try with the whole data,
here
> > it is.
> >
> > xdata=(0:0.0145:1.4065)';
> > ydata=[ 0.5411; 0; 0; 0.9626;
> > 0.5411; 0.5054; 0; 1.1279; 0.5054;
> > 0.4759; 0.7227; 0.8957; 0.4759;
> > 0.4759; 0.8957; 0.6797; 0.4510; 0.4759;
> > 1.1102; 0; 0.4510; 0.6797;
> > 0.9273; 0.4297; 0.4297; 0.9273; 0.6435;
> > 0.4111; 0.4297; 0.9273; 0.5857;
> > 0.5621; 0.6126; 0.8810; 0.4111; 0.5621;
> > 0.8810; 0.7227; 0.5411; 0;
> > 0.9480; 0.4111; 0.5411; 0.5857; 0.9480;
> > 0.5411; 0.5411; 0.7227; 0.7227;
> > 0.5411; 0.3948; 0.9480; 0.5621; 0.6435;
> > 0; 0.8810; 0; 0.5411;
> > 0.7227; 0.8810; 0.5621; 0.5621; 0.8810;
> > 0.7565; 0.4111; 0; 1.0989;
> > 0.6126; 0.4111; 0.4297; 1.1593; 0.4510;
> > 0.4297; 0.6435; 1.0472;
> > 0; 0; 0.6797; 0.9818; 0;
0;
> > 0.8957; 0.7227; 0; 0.4759;
> > 1.0472; 0.7227; 0; 0.5054; 1.1279;
> > 0.7227; 0; 0.7227; 1.0472;
> > 0.5411; 0; 0.7752; 1.0472];
> >
> > really appreciate your help.
> >
> > Best Regards. Rafael
> >
>

Tags for this Thread

What are tags?

A tag is like a keyword or category label associated with each thread. Tags make it easier for you to find threads of interest.

Anyone can tag a thread. Tags are public and visible to everyone.

Contact us