sinusoidfitter

Postscript file available

Non-linear least-squares fit sinusoid to data

include "sinusoidfitter.g"

Constructors
Functions

eval	Evaluate a sinusoid (usually the result of a fit)
fit	Fit a sinusoid to data
setstate	Set the state of the fitter

Description

sinusoidfitter is a class that does a non-linear least-squares fit of a sinusoid of a specified order to data. An optional estimate of the errors in the data may be provided. The state of the fitter, which includes the parameters of the sinusoid, estimated errors, the covariance matrix, the number of iterations actually done, and whether the fit converged or not, is available at any time and is set whenever the fit function is invoked.. A function to evaluate the fitted sinusoid is also provided. The controlling values of the fitter (initial guesses for the parameters, the maximum number of iterations, and the convergence criteria) can be set by the user.

At present, only fits to real data without constraints are supported.

The fits are carried out with double precision.

The form of the sinusoid is: Acos(2 $\pi$ (x-x₀)/p)

It is best to remove at a mean value from the data before the fit since the sinusoid has zero mean.

The fit is not particularly robust against aliasing so care should be taken to inspect the result at all times. The fitter does best when the initial period guess is longer than the suspected best fit.

The computation is actually performed in a C++ executable, not in Glish. Details about the algorithms are available in the documentation of the underlying C++ fitting classes.

sinusoidfitterdemo() and sinusoidfittertest() functions are available.

Example

First we need to get access to the sinusoid fitting declarations and create a fitting tool:

  include "sinusoidfitter.g"
  fitter := sinusoidfitter()

Now lets manufacture some data from the sinusoid

y = 3cos2 $\displaystyle \pi$ (x-8)/10

  x := [0:99]/5.0  # 0,0.2,0.4,0.6,0.8,1.0,1.2,...,19.8
  y := 3*cos(2*pi*(x-8)/10)

The fitter needs a reasonable initial guess (this is obviously contrived, since we know what the sinusoid is).

  guess := [=];
  guess.amplitude := 1.0
  guess.period := 13.0
  guess.x0 := 0.0
  fitter.setstate(guess)

Now, lets perform the fit:

  ok := fitter.fit(state, x, y)
  # and examination of the state record shows that the fit
  # has not yet converged.  However, the current values of
  # the sinusoid parameters appear to be converging.
  # So we attempt the fit again.  Since we have not provided
  # any new initial guess via the setstate member, the fit
  # will use the most recent values as its starting point.
  ok := fitter.fit(state, x, y)
  # the fit has now converged according to the state record.

The state argument is a record which will contain the state of the fitter at the end of the fit (either maxiter will have been reached or the fit will have converged). The other two arguments, x and y, are both input arguments and are not changed by using this function. We did not provide an error estimate $\sigma$ , so the software assumes the errors are identical for each data point and $\sigma$ = 1. The state record will indicate if the fit converged or not. The arguments of the fit function are:

ok: The return parameter. It is false only when some severe error has occurred. If a fit has merely failed to converge, the state record will show that and the return value will still be true.
state: The state of the fitter after the fit. See the description of the fit function for a full listing of all of the fields in the state record.
x,y: The data - x and y must be vectors having the same number of elements.
sigma: An optional argument indicating the errors in the data values, y. When absent, it is assumed that all y values have the same error.

The other function in this class is one to evaluate a sinusoid. It is the sinusoid given by the current state of the fitter. Usually this will be the sinusoid you have just fit, but it could be any sinusoid you specify. For example, to find out the largest deviation of the fit values from the data values we could do the following:

  ok := fitter.eval(y2, x)
  print max(abs(y2 - y))

ok: The return value is True if the fit succeeds, False otherwise. (Output)
y2: The sinusoid evaluated at x. (Output)
x: Where to evaluate the sinusoid. (Input)

See Also
sinusoidfitterdemo sinusoidfittertest

Next: sinusoidfitter.fit - Function Up: mathematics - Module Previous: polyfittertest - Function Contents Index

	Package	utility
	Module	mathematics

sinusoidfitter - Tool