Re: fitting points to a sin function

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@xxxxxxxxxxxxxxxxxx>
wrote in message <fpu2l6$5ts$1@xxxxxxxxxxxxxxxxxx>...
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


.

Relevant Pages

  • More on serendipity
    ... I agree that if one fits the linear ... in the presence of noise if one fit both models would b1 ... equations and nonlinear regression equations is not empty, ...
    (sci.stat.math)
  • Re: Jeff Hawkins Q&A
    ... >> features and aspects of complexity and non linearity but don't really ... Computer power isn't the issue for nonlinear complexity. ... But linear complexity of itself has nothing to ... to the mechanization of intelligence for that reason. ...
    (comp.ai.philosophy)
  • Re: Any ideas how to linearize (and evaluate) this nonlinear equation?
    ... >nutrient, X0 is the total amount of the nutrient contained in all ... >point for the nonlinear one. ... if I had done a linear ... These four kinds of confidence intervals which you ...
    (sci.math.num-analysis)
  • Re: Is the Universe really that complicated?
    ... analysis vs abstraction problem (aka the egg/chicken ... Let's reformulate it into "Is the universe a linear ... nonlinear system and that no problem is really linear ... One odd case would be the matrix algebra where matrix ...
    (sci.math)
  • Re: Digicams With MF Film Quality
    ... so I don't see how any mapping can be linear ... Assume we know how a Dirac delta function gets imaged (ie we know the ... The blurred image is then a sum of PSFs. ... The positivity constraint makes the method nonlinear. ...
    (rec.photo.digital)
Loading