given a set of x,y points I need a procedure to give me the
a values for a polynomium fitting the points.

I think the procedure head would look like this

procedure polynomiumfit(const x : array of real;
                        const y : array of real;
                        const order : integer;
                        var   a : array of real);

y(x) = a0 + a1* x + a2 * x*x + a3*x*x*x .....

A least-square fit is ok :
 

Can anyone point to a function/procedure like this?
}

PROGRAM LineFit;
Uses    Crt, Graph;

CONST   MaxK   = 10;     { max. number of fitting functions }
        MaxN   = 400;    { max. number of xy-points }
        { MaxK*MaxN*SizeOf(Real) must be smaller then 64K }

TYPE    Real = Double;   { it's faster and more exact }
        Vector = ARRAY[1..MaxN] OF Real;
        Matrix = ARRAY[1..MaxK,1..MaxK] OF Real;
 

PROCEDURE Gauss(n : Integer;  M: matrix;  Y: Vector; VAR X : Vector;
                     VAR det : Real;  eps : Real);
{ Solution a Linaer-Algebraic-Equitation-Sytem (M*x=y)    }
{ with 'Gauss-Jordan-Elimination' and Pivoting            }
LABEL      10, 20 ;
VAR        i, j, k, maxindex        : Integer;
           change                   : Integer;
           maxel, temp, wert        : Real;
BEGIN
  change:=0;
  FOR i:=1 TO n-1 DO
    BEGIN
      maxindex:=i;    maxel:=Abs(M[i,i]);
      FOR j:=i+1 TO n DO
       BEGIN   IF Abs(M[j,i]) > maxel
                  THEN BEGIN    maxel:=Abs(M[j,i]);
                                maxindex:=j;
                       END;
       END;
      IF maxindex <> i THEN
        BEGIN   FOR j:=i TO n DO
                  BEGIN  temp:=M[i,j];
                         M[i,j]:=M[maxindex,j];
                         M[maxindex,j]:=temp;
                   END;
                 temp:=Y[i];   Y[i]:=Y[maxindex];
                 Y[maxindex]:=temp;
                 Inc(change);
        END;
      IF Abs(M[i,i]) < eps THEN GOTO 10;
      FOR j:=i+1 TO n DO
        BEGIN  wert:=M[j,i]/M[i,i];
               FOR k:=i+1 TO n DO
                   M[j,k]:=M[j,k] - wert*M[i,k];
               Y[j]:=Y[j] - wert*Y[i];
        END;
    END;
  IF Abs(M[n,n]) < eps THEN GOTO 10;
  det:=1;
  FOR i:=n DOWNTO 1 DO
    BEGIN  wert:=0;
           FOR j:=i+1 TO n DO
             wert:=wert + M[i,j]*X[j];
           X[i]:=(Y[i] - wert)/M[i,i];
           det:=det*M[i,i];
    END;
  IF Odd(change) THEN det:=-det;
  GOTO 20;
  10 :  det:=0;    { Determinant<eps ==> System maybe singular }
  20 :   ;
END;  {---- Gauss ----}
 

PROCEDURE FKT(k:Integer; xarg: Real; VAR PHI: Vector);
{ for example Polynomials }
VAR   i  : Integer;
BEGIN
  PHI[1]:=1;    { polynomials }
  FOR i:=2 TO k DO PHI[i]:=PHI[i-1]*xarg;

{ Another set of 'better' (for the fitting problem) functions can be used too }
{ These functions must be 'linearly independent' (polynomials are sure) }
{
  PHI[1]:=1;
  PHI[2]:=1/(xarg+1);
  PHI[3]:=xarg/(xarg*xarg+1);
  PHI[4]:=sin(xarg);
  PHI[5]:=exp(xarg);
  ....
  ....
 }
END;   {---- FKT ----}
 

FUNCTION CalcFitting(k : Integer; A: Vector; xwert : Real) : Real;
VAR   PHI  : Vector;
      sum  : Real;
      i    : Integer;
BEGIN
  FKT(k,xwert,PHI);
  sum:=0;
  FOR i:=1 TO k DO sum:=sum + A[i]*PHI[i];
  CalcFitting:=sum;
END;
 

PROCEDURE DoFitting(n,k : Integer; X,Y : Vector; VAR A : Vector);
TYPE FType       = ARRAY[1..MaxN,1..MaxK] OF Real;
VAR  i, j, l     : Integer;
     PHI, B      : Vector;
     F           : ^FType;
     M           : ^Matrix;
     det         : Real;
BEGIN
  New(M);   New(F);

  FOR i:=1 TO n DO         { Prepare the approximation }
   BEGIN
    FKT(k,X[i],PHI);
    FOR j:=1 TO k DO F^[i,j]:=PHI[j];
   END;
  FOR j:=1 TO k DO         { Build the matrix of the LinEqu. }
   FOR i:=1 TO k DO
    BEGIN  M^[i,j]:=0;
      FOR l:=1 TO N DO
        M^[i,j]:=M^[i,j] +  F^[l,i]*F^[l,j];
      M^[j,i]:=M^[i,j];
    END;
  FOR i:=1 TO k DO         { Build the right side of the LinEqu.}
   BEGIN  B[i]:=0;
     FOR l:=1 TO n DO B[i]:=B[i] + F^[l,i]*Y[l];
   END;

  Gauss(k, M^, B, A,  det, 1.0e-15);    { Solve the system }

  IF det=0 THEN { Missing Errorhandling };

  Dispose(F);  Dispose(M);
END;   {---- DoFitting ----}
 

VAR   X, Y, A                          : Vector;
      n, i, k, GraphDriver, GraphMode  : Integer;
      delta, Zero, xpos, ypos          : Real;
      ch                               : Char;

BEGIN     { ---- MAIN ---- }
  n:=200;      { 100 Points }
  GraphDriver:=Detect;
  InitGraph(GraphDriver,GraphMode,'');

  delta:=(GetMaxX)/(n-1);  Zero:=GetMaxY/2;
  FOR i:=1 TO n DO      { Build an Example }
   BEGIN
    x[i]:=(i-1)*delta;
    y[i]:=Zero + (Sin(X[i]/150)*Zero*0.5)*(X[i]/400) + (Random-0.5)*25 ;
   END;

  delta:=GetMaxX/500;    { Show 7 line-fittings with polynomials }
  FOR k:=1 TO 7 DO       { k = (degree of the polynomial)-1 }
   BEGIN
    ClearDevice;
    SetColor(Red);          { Show the Points }
    FOR i:=1 TO n DO Circle(Round(x[i]),Round(y[i]),2);
    DoFitting(n,k, X,Y, A); { Do the calculation }
    FOR i:=0 TO 500 do      { Show the curves with 501 points }
     BEGIN
       xpos:=i*delta;
       ypos:=CalcFitting(k,A,xpos);
       PutPixel(Round(xpos),Round(ypos),((9+k) MOD 15)+1);
     END;
    ch:=ReadKey;
   END;
  CloseGraph;
END.

This page hosted by  Get your own Free Homepage