MATLAB Newsgroup

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

"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

"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

"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

> 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

"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

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

>

> 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

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

> >

>

You can think of your watch list as threads that you have bookmarked.

You can add tags, authors, threads, and even search results to your watch list. This way you can easily keep track of topics that you're interested in. To view your watch list, click on the "My Newsreader" link.

To add items to your watch list, click the "add to watch list" link at the bottom of any page.

To add search criteria to your watch list, search for the desired term in the search box. Click on the "Add this search to my watch list" link on the search results page.

You can also add a tag to your watch list by searching for the tag with the directive "tag:tag_name" where tag_name is the name of the tag you would like to watch.

To add an author to your watch list, go to the author's profile page and click on the "Add this author to my watch list" link at the top of the page. You can also add an author to your watch list by going to a thread that the author has posted to and clicking on the "Add this author to my watch list" link. You will be notified whenever the author makes a post.

To add a thread to your watch list, go to the thread page and click the "Add this thread to my watch list" link at the top of the page.

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.

The newsgroups are a worldwide forum that is open to everyone. Newsgroups are used to discuss a huge range of topics, make announcements, and trade files.

Discussions are threaded, or grouped in a way that allows you to read a posted message and all of its replies in chronological order. This makes it easy to follow the thread of the conversation, and to see what’s already been said before you post your own reply or make a new posting.

Newsgroup content is distributed by servers hosted by various organizations on the Internet. Messages are exchanged and managed using open-standard protocols. No single entity “owns” the newsgroups.

There are thousands of newsgroups, each addressing a single topic or area of interest. The MATLAB Central Newsreader posts and displays messages in the comp.soft-sys.matlab newsgroup.

**MATLAB Central**

You can use the integrated newsreader at the MATLAB Central website to read and post messages in this newsgroup. MATLAB Central is hosted by MathWorks.

Messages posted through the MATLAB Central Newsreader are seen by everyone using the newsgroups, regardless of how they access the newsgroups. There are several advantages to using MATLAB Central.

**One Account**

Your MATLAB Central account is tied to your MathWorks Account for easy access.

**Use the Email Address of Your Choice**

The MATLAB Central Newsreader allows you to define an alternative email address as your posting address, avoiding clutter in your primary mailbox and reducing spam.

**Spam Control**

Most newsgroup spam is filtered out by the MATLAB Central Newsreader.

**Tagging**

Messages can be tagged with a relevant label by any signed-in user. Tags can be used as keywords to find particular files of interest, or as a way to categorize your bookmarked postings. You may choose to allow others to view your tags, and you can view or search others’ tags as well as those of the community at large. Tagging provides a way to see both the big trends and the smaller, more obscure ideas and applications.

**Watch lists**

Setting up watch lists allows you to be notified of updates made to postings selected by author, thread, or any search variable. Your watch list notifications can be sent by email (daily digest or immediate), displayed in My Newsreader, or sent via RSS feed.

- Use a newsreader through your school, employer, or internet service provider
- Pay for newsgroup access from a commercial provider
- Use Google Groups
- Mathforum.org provides a newsreader with access to the comp.soft sys.matlab newsgroup
- Run your own server. For typical instructions, see: http://www.slyck.com/ng.php?page=2