functions::QuadraticFit Class Reference

Functor(-ish) class to calculate fit of a quadratic to N data-points. More...

#include <QuadraticFit.h>

Public Member Functions
	QuadraticFit (int length)
	constructor Taking the number of samples that will be fit to when using this instance
template<typename InputIterator >
void	Fit (InputIterator i_in, double &a, double &b, double &c)
	Fit a given set of samples.
int	GetSize () const
Private Attributes
int	fN
double	fDet
double	fX2
double	fX4

Detailed Description

Functor(-ish) class to calculate fit of a quadratic to N data-points.

Author:: Ben Krikler

Assumes that the x values are spaced evenly and centred at 0 so that the sum over odd powers of x is 0. Uses a formula derived from least squares regression

Definition at line 14 of file QuadraticFit.h.

Constructor & Destructor Documentation

functions::QuadraticFit::QuadraticFit ( int length ) [inline]

constructor Taking the number of samples that will be fit to when using this instance

Definition at line 18 of file QuadraticFit.h.

References fDet, fN, fX2, and fX4.

00018                           :fN(length+(length+1)%2),fX2(0),fX4(0){
00019    int range=(fN-1)/2;
00020    for(double x= -range; x <= range; ++x){
00021      fX2+=x*x;
00022      fX4+=x*x*x*x;
00023    }
00024    fDet=fX4*fN*fX2 - fX2*fX2*fX2;
00025   }

Member Function Documentation

template<typename InputIterator >

void functions::QuadraticFit::Fit	(	InputIterator	i_in,
		double &	a,
		double &	b,
		double &	c
	)			`[inline]`

Fit a given set of samples.

The input iterator should be the first sample to fit. Since we assume the central sample is the origin, (which simplifies the maths considerably), the initial iterator should be set back N/2 samples from the centre of the fit, where N is the number of samples to fit. For example, to fit a peak at X using 6 samples, the input iterator should point to X-3.

Definition at line 34 of file QuadraticFit.h.

References fDet, fN, fX2, and fX4.

Referenced by TemplateConvolver::FitPeak().

00034                                                                {
00035    double Sy=0, Sxy=0, Sx2y=0, y=0;
00036    int range=(fN-1)/2;
00037    for(double x= -range; x <= range; ++x){
00038      y=*i_in;
00039      Sy+=y;
00040      Sxy+=x*y;
00041      Sx2y+=x*x*y;
00042      //DEBUG_VALUE(x, y,Sy, Sxy, Sx2y);
00043      ++i_in;
00044    }
00045    double Da=Sx2y*fX2*fN  - Sy  *fX2*fX2;
00046    double Db=Sxy *fX4*fN  - Sxy *fX2*fX2;
00047    double Dc=Sy  *fX4*fX2 - Sx2y*fX2*fX2;
00048    a=Da/fDet;
00049    b=Db/fDet;
00050    c=Dc/fDet;
00051   }