RSPLINES

Robust B-Splines

Applying B-Spline calculation on noisy data may result in an unwanted oszillation or number of free params.

This class here sees persues a least squares approach on the input data. It merly tries to optimize the B-Spline on only a few "breaks" instead of fitting it on the complete data and in addition tries to minimize the minimum least squares error when fitting the spline through the input data cloud.

Reference

The implementation here follows the reference on the matlab file found here.

In addition to the plain ppline parameter calcultion there is support for:

• Integral and Differential
• Periodic extension
• Robust estimation with outlier detection (deduction parameter beta)
• Nonuniform distribution of "breaks".
• Different polynom orders (1 - x)

Interface description

Creation

   TRobustBSpline
   public
     constructor Create( Order : integer = 4; beta : double = 0; periodicExpansion : boolean = False); 
   end;

• Order: Spline Order, default is cubic spline with order = 4;
• Beta: Beta = 0 -> non robust least squares solving. If greater than 0 three iterations of weighted least squares are performed. Weights are computed from previous residuals. BETA close to 0 gives all data equal weighting. Increase BETA to reduce the influence from outlying data. BETA close to 1 may cause instability or rank deficiency.
• PeriodicExpansion: If set to true periodic data extension is performed.

Applying data

   procedure InitSpline(const x, y, breaks : TDoubleDynArray); overload;
   procedure InitSpline(x, y, breaks : IMatrix); overload; 

The initialization is performed in 2D. X and Y need to be of same length in case the dynamic array init funtion is used or two vectors (width=1, height=n) of same length. The x values need to sorted (low values first).

The breaks array/vector contains x values. These values need within the space that is covered by X. The maximum number of breaks is n-2 where n is the x vector length. The minimum length of the breaks is 1 + order that was used in the constructor.

Spline Evaluation

   function EvalSpline(const x : TDoubleDynArray) : TDoubleDynArray; overload;
   function EvalSpline(x : IMatrix) : IMatrix; overload;

   // calculates the jth derrivative of the spline at points x
   function DiffSpline(const x : TDoubleDynArray; j : integer) : TDoubleDynArray; overload;
   function DiffSpline(x : IMatrix; j : integer) : IMatrix; overload;
   function IntSpline( const intFrom, intTo : double) : double; overload;      // definitive integral
   function IntSpline(const x : TDoubleDynArray) : TDoubleDynArray; overload;  // indefinite integral @ breaks[0]
   function IntSpline(a : double; const x : TDoubleDynArray) : TDoubleDynArray; overload;  // indefinite integral
   function IntSpline(x : IMatrix) : IMatrix; overload;  // indefinite integral @ breaks[0]
   function IntSpline(a : double; x : IMatrix) : IMatrix; overload;  // indefinite integral

• EvalSpline: Evaluates the spline at points defined by x. Result is the y vector.
• DiffSpline: Returns the jth derrivative of the spline at points x.
• IntSpline: Returns the splines integral between two x values or returns the integral function values at points x. A is an arbitrary integral constant.

Exmaples

The examples shown here are from the test program. The example show how to generate simple input with random noise data. As input a linear combination of two sines is used. The input is sampled for 200 input values (with noise) - support points (aka break) are only defined on 10 x values. Note that the x values are not constantly spaced.

   // ##################################################
   // #### Class definition
   TTestRBSpline = class(TBaseMatrixTestCase)
   private
     fRnd : TRandomGenerator;
     fNoiseAmpl : double;
    
     procedure GenSplineInput(var x, y, breaks : IMatrix; noiseAmpl : double; cl : TDoubleMatrixClass);

     procedure SinFunc(var value : double);
   published
     procedure SingleThreadSpline;
  

   // ###################################################
   // ##### Here the implemenation part
   implementation 
   
   procedure TTestRBSpline.SinFunc(var value: double);
   begin
        value := Sin(Value) + sin(2*value) + fNoiseAmpl*fRnd.RandGauss;
   end;

   procedure TTestRBSpline.GenSplineInput(var x, y, breaks: IMatrix; noiseAmpl : double; cl : TDoubleMatrixClass);
   begin
        fNoiseAmpl := noiseAmpl;
        fRnd := TRandomGenerator.Create(raMersenneTwister);
        try
           x := cl.CreateRand(1, 200);
           x.ScaleInPlace(2*pi);
           x.SortInPlace(False);

           y := x.ElementwiseFunc( {$IFDEF FPC}@{$ENDIF}sinFunc );
        finally
               fRnd.Free;
        end;

        breaks := cl.CreateLinSpace(10, 0, 2*pi);
   end;

   
   // #########################################
   // #### main routine
   procedure TTestRBSpline.SingleThreadSpline;
   var x, y, breaks : IMatrix;
       rb : TRobustBSpline;
       xx, yy : IMatrix;
       counter: Integer;
   begin
        TRobustBSpline.DefMatrixClass := TDoubleMatrix;

        GenSplineInput(x, y, breaks, 0.2, TDoubleMatrix);
     
        rb := TRobustBSpline.Create;
        try
           rb.InitSpline(x, y, breaks);

           xx := TDoubleMatrix.CreateLinSpace(400, 0, 2*pi);
           yy := rb.EvalSpline(xx);
    
           for counter := 0 to yy.VecLen - 1 do
               Check( Abs( yy.Vec[counter] - ( sin(xx.Vec[counter]) + sin(2*xx.Vec[counter]) )) < 0.25, 'Spline evaluation failed @' + IntToStr(counter) );
        finally
               rb.Free;
        end;
   end;