org.apache.commons.math.optimization.fitting
Class HarmonicCoefficientsGuesser

java.lang.Object
  extended by org.apache.commons.math.optimization.fitting.HarmonicCoefficientsGuesser

public class HarmonicCoefficientsGuesser
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.

Since:
2.0
Version:
$Revision: 786479 $ $Date: 2009-06-19 08:36:16 -0400 (Fri, 19 Jun 2009) $

Field Summary
private  double a
          Guessed amplitude.
private  WeightedObservedPoint[] observations
          Sampled observations.
private  double omega
          Guessed pulsation ω.
private  double phi
          Guessed phase φ.
 
Constructor Summary
HarmonicCoefficientsGuesser(WeightedObservedPoint[] observations)
          Simple constructor.
 
Method Summary
 double getGuessedAmplitude()
          Get the guessed amplitude a.
 double getGuessedPhase()
          Get the guessed phase φ.
 double getGuessedPulsation()
          Get the guessed pulsation ω.
 void guess()
          Estimate a first guess of the coefficients.
private  void guessAOmega()
          Estimate a first guess of the a and ω coefficients.
private  void guessPhi()
          Estimate a first guess of the φ coefficient.
private  void sortObservations()
          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

observations

private final WeightedObservedPoint[] observations
Sampled observations.


a

private double a
Guessed amplitude.


omega

private double omega
Guessed pulsation ω.


phi

private double phi
Guessed phase φ.

Constructor Detail

HarmonicCoefficientsGuesser

public HarmonicCoefficientsGuesser(WeightedObservedPoint[] observations)
Simple constructor.

Parameters:
observations - sampled observations
Method Detail

guess

public void guess()
           throws OptimizationException
Estimate a first guess of the coefficients.

Throws:
OptimizationException - if the sample is too short or if the first guess cannot be computed (when the elements under the square roots are negative).

sortObservations

private void sortObservations()
Sort the observations with respect to the abscissa.


guessAOmega

private void guessAOmega()
                  throws OptimizationException
Estimate a first guess of the a and ω coefficients.

Throws:
OptimizationException - if the sample is too short or if the first guess cannot be computed (when the elements under the square roots are negative).

guessPhi

private void guessPhi()
Estimate a first guess of the φ coefficient.


getGuessedAmplitude

public double getGuessedAmplitude()
Get the guessed amplitude a.

Returns:
guessed amplitude a;

getGuessedPulsation

public double getGuessedPulsation()
Get the guessed pulsation ω.

Returns:
guessed pulsation ω

getGuessedPhase

public double getGuessedPhase()
Get the guessed phase φ.

Returns:
guessed phase φ


Copyright (c) 2003-2010 Apache Software Foundation