MATLAB Fitting Function

Question

I am trying to fit a line to some data without using polyfit and polyval. I got some good help already on how to implement this and I have gotten it to work with a simple sin function. However, when applied to the function I am trying to fit, it does not work. Here is my code:

clear all 
clc

lb=0.001; %lowerbound of data 
ub=10; %upperbound of data 
step=.1; %step-size through data 
a=.03; 
la=1482/120000; %1482 is speed of sound in water and 120kHz
ep1=.02;
ep2=.1;
x=lb:step:ub;
r_sq_des=0.90; %desired value of r^2 for the fit of data without noise present 


i=1; 

for x=lb:step:ub 
     G(i,1)= abs(sin((a/la)*pi*x*(sqrt(1+(1/x)^2)-1)));
     N(i,1)=2*rand()-1; 
     Ghat(i,1)=(1+ep1*N(i,1))*G(i,1)+ep2*N(i,1); 
     r(i,1)=x; 
     i=i+1; 
end 

x=r;
y=G;
V=[x.^0];
Vfit=[x.^0];


for i=1:1:1000
    V = [x.^i V];
    c = V \ y;

    Vfit = [x.^i Vfit];
    yFit=Vfit*c;

    plot(x,y,'o',x,yFit,'--')

    drawnow
    pause
end

The first two sections are just defining variables and the function. The second for loop is where I am making the fit. As you can see, I have it pause after every nth order in order to see the fit.

Answer 1

As dennis mentioned, a different set of basis functions might do better. However you can improve the polynomial fit with QR factorisation, rather than just \ to solve the matrix equation. It is a badly conditioned problem no matter what you do however, and using smooth basis functions wont allow you to accurately reproduce the sharp corners in the actual function.

clear all
close all
clc

lb=0.001; %lowerbound of data 
ub=10; %upperbound of data 
step=.1; %step-size through data 
a=.03; 
la=1482/120000; %1482 is speed of sound in water and 120kHz
ep1=.02;
ep2=.1;
x=logspace(log10(lb),log10(ub),100)';
r_sq_des=0.90; %desired value of r^2 for the fit of data without noise present 

y=abs(sin(a/la*pi*x.*(sqrt(1+(1./x).^2)-1)));
N=2*rand(size(x))-1;
Ghat=(1+ep1*N).*y+ep2*N;

V=[x.^0];

xs=(lb:.01:ub)';
Vfit=[xs.^0];

for i=1:1:20%length(x)-1
    V = [x.^i V];
    Vfit = [xs.^i Vfit];
    [Q,R]=qr(V,0);
    c = R\(Q'*y);
    yFit=Vfit*c;

    plot(x,y,'o',xs,yFit)
    axis([0 10 0 1])
    drawnow
    pause
end

Answer 2

I changed your fit formula a bit, I got the same answers but quickly got a warning that the matrix was singular. No sense in continuing past the point that the inversion is singular.

Depending on what you are doing you can usually change out variables or change domains. This doesn't do a lot better, but it seemed to help a little bit.

I increased the number of samples by a factor of 10 since the initial part of the curve didn't look sampled highly enough.

I added a weighting variable but it is set to equal weight in the code below. Attempts to deweight the tail didn't help as much as I hoped.

Probably not really a solution, but perhaps will help with a few more knobs/variables.

    ...
    step=.01; %step-size through data 
    ...
    x=r;
    y=G;
    t=x.*sqrt(1+x.^(-2));
    t=log(t);
    V=[ t.^0];
    w=ones(size(t));


    for i=1:1:1000
         %  Trying to solve for value of c
         %  c that 
         %   yhat = V*c  approximates y
         %  or     y = V*c
         %      V'*y = V'*V * c
         %         c = (V'*V) \  V'*y
         V = [t.^i V];
         c =  (V'*diag(w.^2)*V ) \ (V'*diag(w.^2)*y) ;

         yFit=V*c;
         subplot(211)
         plot(t,y,'o',t,yFit,'--')
         subplot(212)
         plot(x,y,'o',x,yFit,'--')

         drawnow
         pause
    end

It looks like more of a frequency estimation problem, and trying to fit a unknown frequency with polynomial tends to be touch and go. Replacing the polynomial basis with a quick sin/cos basis didn't seem to do to bad.

    V = [sin(t*i) cos(t*i) V];

Unless you specifically need a polynomial basis, you can apply your knowledge of the problem domain to find other potential basis functions for your fit, or to attempt to make the domain in which you are performing the fit more linear.