package org.spaceroots.mantissa.fitting;
import java.io.Serializable;
import org.spaceroots.mantissa.functions.FunctionException;
import org.spaceroots.mantissa.functions.ExhaustedSampleException;
import org.spaceroots.mantissa.functions.vectorial.SampledFunctionIterator;
import org.spaceroots.mantissa.functions.vectorial.VectorialValuedPair;
import org.spaceroots.mantissa.quadrature.vectorial.EnhancedSimpsonIntegratorSampler;
import org.spaceroots.mantissa.estimation.EstimationException;
/** 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.
* @version $Id: HarmonicCoefficientsGuesser.java 1709 2006-12-03 21:16:50Z luc $
* @author L. Maisonobe
*/
public class HarmonicCoefficientsGuesser
implements Serializable{
public HarmonicCoefficientsGuesser(AbstractCurveFitter.FitMeasurement[] measurements) {
this.measurements =
(AbstractCurveFitter.FitMeasurement[]) measurements.clone();
a = Double.NaN;
omega = Double.NaN;
}
/** Estimate a first guess of the coefficients.
* @exception ExhaustedSampleException if the sample is exhausted.
* @exception FunctionException if the integrator throws one.
* @exception EstimationException if the sample is too short or if
* the first guess cannot be computed (when the elements under the
* square roots are negative).
* */
public void guess()
throws ExhaustedSampleException, FunctionException, EstimationException {
guessAOmega();
guessPhi();
}
/** Estimate a first guess of the a and ω coefficients.
* @exception ExhaustedSampleException if the sample is exhausted.
* @exception FunctionException if the integrator throws one.
* @exception EstimationException 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 guessAOmega()
throws ExhaustedSampleException, FunctionException, EstimationException {
// initialize the sums for the linear model between the two integrals
double sx2 = 0.0;
double sy2 = 0.0;
double sxy = 0.0;
double sxz = 0.0;
double syz = 0.0;
// build the integrals sampler
F2FP2Iterator iter = new F2FP2Iterator(measurements);
SampledFunctionIterator sampler =
new EnhancedSimpsonIntegratorSampler(iter);
VectorialValuedPair p0 = sampler.nextSamplePoint();
double p0X = p0.x;
double[] p0Y = p0.y;
// get the points for the linear model
while (sampler.hasNext()) {
VectorialValuedPair point = sampler.nextSamplePoint();
double pX = point.x;
double[] pY = point.y;
double x = pX - p0X;
double y = pY[0] - p0Y[0];
double z = pY[1] - p0Y[1];
sx2 += x * x;
sy2 += y * y;
sxy += x * y;
sxz += x * z;
syz += y * z;
}
// compute the amplitude and pulsation coefficients
double c1 = sy2 * sxz - sxy * syz;
double c2 = sxy * sxz - sx2 * syz;
double c3 = sx2 * sy2 - sxy * sxy;
if ((c1 / c2 < 0.0) || (c2 / c3 < 0.0)) {
throw new EstimationException("unable to guess a first estimate");
}
a = Math.sqrt(c1 / c2);
omega = Math.sqrt(c2 / c3);
}
/** Estimate a first guess of the φ coefficient.
* @exception ExhaustedSampleException if the sample is exhausted.
* @exception FunctionException if the sampler throws one.
*/
private void guessPhi()
throws ExhaustedSampleException, FunctionException {
SampledFunctionIterator iter = new FFPIterator(measurements);
// initialize the means
double fcMean = 0.0;
double fsMean = 0.0;
while (iter.hasNext()) {
VectorialValuedPair point = iter.nextSamplePoint();
double omegaX = omega * point.x;
double cosine = Math.cos(omegaX);
double sine = Math.sin(omegaX);
fcMean += omega * point.y[0] * cosine - point.y[1] * sine;
fsMean += omega * point.y[0] * sine + point.y[1] * cosine;
}
phi = Math.atan2(-fsMean, fcMean);
}
public double getOmega() {
return omega;
}
public double getA() {
return a;
}
public double getPhi() {
return phi;
}
private AbstractCurveFitter.FitMeasurement[] measurements;
private double a;
private double omega;
private double phi;
private static final long serialVersionUID = 2400399048702758814L;
}
Ads help maintain this website.