org.apache.commons.math.optimization.fitting
public class HarmonicCoefficientsGuesser extends java.lang.Object
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.
Modifier and Type | Field and Description |
---|---|
private double |
a
Guessed amplitude.
|
private WeightedObservedPoint[] |
observations
Sampled observations.
|
private double |
omega
Guessed pulsation ω.
|
private double |
phi
Guessed phase φ.
|
Constructor and Description |
---|
HarmonicCoefficientsGuesser(WeightedObservedPoint[] observations)
Simple constructor.
|
Modifier and Type | Method and Description |
---|---|
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.
|
private final WeightedObservedPoint[] observations
private double a
private double omega
private double phi
public HarmonicCoefficientsGuesser(WeightedObservedPoint[] observations)
observations
- sampled observationspublic void guess() throws OptimizationException
OptimizationException
- if the sample is too short or if
the first guess cannot be computed (when the elements under the
square roots are negative).private void sortObservations()
private void guessAOmega() throws OptimizationException
OptimizationException
- if the sample is too short or if
the first guess cannot be computed (when the elements under the
square roots are negative).private void guessPhi()
public double getGuessedAmplitude()
public double getGuessedPulsation()
public double getGuessedPhase()
Copyright (c) 2003-2016 Apache Software Foundation