org.apache.commons.math3.fitting
Class HarmonicFitter.ParameterGuesser

java.lang.Object
  extended by org.apache.commons.math3.fitting.HarmonicFitter.ParameterGuesser
Enclosing class:
HarmonicFitter

public static class HarmonicFitter.ParameterGuesser
extends java.lang.Object

This class guesses harmonic coefficients from a sample.

The algorithm used to guess the coefficients is as follows:

We know f (t) at some sampling points ti and want to find a, ω and φ such that f (t) = a cos (ω t + φ).

From the analytical expression, we can compute two primitives :

     If2  (t) = ∫ f2  = a2 × [t + S (t)] / 2
     If'2 (t) = ∫ f'2 = a2 ω2 × [t - S (t)] / 2
     where S (t) = sin (2 (ω t + φ)) / (2 ω)
 

We can remove S between these expressions :

     If'2 (t) = a2 ω2 t - ω2 If2 (t)
 

The preceding expression shows that If'2 (t) is a linear combination of both t and If2 (t): If'2 (t) = A × t + B × If2 (t)

From the primitive, we can deduce the same form for definite integrals between t1 and ti for each ti :

   If2 (ti) - If2 (t1) = A × (ti - t1) + B × (If2 (ti) - If2 (t1))
 

We can find the coefficients A and B that best fit the sample to this linear expression by computing the definite integrals for each sample points.

For a bilinear expression z (xi, yi) = A × xi + B × yi, the coefficients A and B that minimize a least square criterion ∑ (zi - z (xi, yi))2 are given by these expressions:


         ∑yiyi ∑xizi - ∑xiyi ∑yizi
     A = ------------------------
         ∑xixi ∑yiyi - ∑xiyi ∑xiyi

         ∑xixi ∑yizi - ∑xiyi ∑xizi
     B = ------------------------
         ∑xixi ∑yiyi - ∑xiyi ∑xiyi
 

In fact, we can assume both a and ω are positive and compute them directly, knowing that A = a2 ω2 and that B = - ω2. The complete algorithm is therefore:


 for each ti from t1 to tn-1, compute:
   f  (ti)
   f' (ti) = (f (ti+1) - f(ti-1)) / (ti+1 - ti-1)
   xi = ti - t1
   yi = ∫ f2 from t1 to ti
   zi = ∫ f'2 from t1 to ti
   update the sums ∑xixi, ∑yiyi, ∑xiyi, ∑xizi and ∑yizi
 end for

            |--------------------------
         \  | ∑yiyi ∑xizi - ∑xiyi ∑yizi
 a     =  \ | ------------------------
           \| ∑xiyi ∑xizi - ∑xixi ∑yizi


            |--------------------------
         \  | ∑xiyi ∑xizi - ∑xixi ∑yizi
 ω     =  \ | ------------------------
           \| ∑xixi ∑yiyi - ∑xiyi ∑xiyi

 

Once we know ω, we can compute:

    fc = ω f (t) cos (ω t) - f' (t) sin (ω t)
    fs = ω f (t) sin (ω t) + f' (t) cos (ω t)
 

It appears that fc = a ω cos (φ) and fs = -a ω sin (φ), so we can use these expressions to compute φ. The best estimate over the sample is given by averaging these expressions.

Since integrals and means are involved in the preceding estimations, these operations run in O(n) time, where n is the number of measurements.


Field Summary
private  double a
          Amplitude.
private  double omega
          Angular frequency.
private  double phi
          Phase.
 
Constructor Summary
HarmonicFitter.ParameterGuesser(WeightedObservedPoint[] observations)
          Simple constructor.
 
Method Summary
 double[] guess()
          Gets an estimation of the parameters.
private  double[] guessAOmega(WeightedObservedPoint[] observations)
          Estimate a first guess of the amplitude and angular frequency.
private  double guessPhi(WeightedObservedPoint[] observations)
          Estimate a first guess of the phase.
private  WeightedObservedPoint[] sortObservations(WeightedObservedPoint[] unsorted)
          Sort the observations with respect to the abscissa.
 
Methods inherited from class java.lang.Object
clone, equals, finalize, getClass, hashCode, notify, notifyAll, toString, wait, wait, wait
 

Field Detail

a

private final double a
Amplitude.


omega

private final double omega
Angular frequency.


phi

private final double phi
Phase.

Constructor Detail

HarmonicFitter.ParameterGuesser

public HarmonicFitter.ParameterGuesser(WeightedObservedPoint[] observations)
Simple constructor.

Parameters:
observations - Sampled observations.
Throws:
NumberIsTooSmallException - if the sample is too short.
ZeroException - if the abscissa range is zero.
MathIllegalStateException - when the guessing procedure cannot produce sensible results.
Method Detail

guess

public double[] guess()
Gets an estimation of the parameters.

Returns:
the guessed parameters, in the following order:
  • Amplitude
  • Angular frequency
  • Phase

sortObservations

private WeightedObservedPoint[] sortObservations(WeightedObservedPoint[] unsorted)
Sort the observations with respect to the abscissa.

Parameters:
unsorted - Input observations.
Returns:
the input observations, sorted.

guessAOmega

private double[] guessAOmega(WeightedObservedPoint[] observations)
Estimate a first guess of the amplitude and angular frequency. This method assumes that the #sortObservations() method has been called previously.

Parameters:
observations - Observations, sorted w.r.t. abscissa.
Returns:
the guessed amplitude (at index 0) and circular frequency (at index 1).
Throws:
ZeroException - if the abscissa range is zero.
MathIllegalStateException - when the guessing procedure cannot produce sensible results.

guessPhi

private double guessPhi(WeightedObservedPoint[] observations)
Estimate a first guess of the phase.

Parameters:
observations - Observations, sorted w.r.t. abscissa.
Returns:
the guessed phase.


Copyright (c) 2003-2013 Apache Software Foundation