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 $Revision$ $Date$ - * @since 2.0 - - */ -public class HarmonicCoefficientsGuesser { - - /** Sampled observations. */ - private final WeightedObservedPoint[] observations; - - /** Guessed amplitude. */ - private double a; - - /** Guessed pulsation ω. */ - private double omega; - - /** Guessed phase φ. */ - private double phi; - - /** Simple constructor. - * @param observations sampled observations - */ - public HarmonicCoefficientsGuesser(WeightedObservedPoint[] observations) { - this.observations = observations.clone(); - a = Double.NaN; - omega = Double.NaN; - } - - /** Estimate a first guess of the coefficients. - * @exception OptimizationException 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 OptimizationException { - sortObservations(); - guessAOmega(); - guessPhi(); - } - - /** Sort the observations with respect to the abscissa. - */ - private void sortObservations() { - - // Since the samples are almost always already sorted, this - // method is implemented as an insertion sort that reorders the - // elements in place. Insertion sort is very efficient in this case. - WeightedObservedPoint curr = observations[0]; - for (int j = 1; j < observations.length; ++j) { - WeightedObservedPoint prec = curr; - curr = observations[j]; - if (curr.getX() < prec.getX()) { - // the current element should be inserted closer to the beginning - int i = j - 1; - WeightedObservedPoint mI = observations[i]; - while ((i >= 0) && (curr.getX() < mI.getX())) { - observations[i + 1] = mI; - if (i-- != 0) { - mI = observations[i]; - } - } - observations[i + 1] = curr; - curr = observations[j]; - } - } - - } - - /** Estimate a first guess of the a and ω coefficients. - * @exception 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 guessAOmega() throws OptimizationException { - - // 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; - - double currentX = observations[0].getX(); - double currentY = observations[0].getY(); - double f2Integral = 0; - double fPrime2Integral = 0; - final double startX = currentX; - for (int i = 1; i < observations.length; ++i) { - - // one step forward - final double previousX = currentX; - final double previousY = currentY; - currentX = observations[i].getX(); - currentY = observations[i].getY(); - - // update the integrals of f2 and f'2 - // considering a linear model for f (and therefore constant f') - final double dx = currentX - previousX; - final double dy = currentY - previousY; - final double f2StepIntegral = - dx * (previousY * previousY + previousY * currentY + currentY * currentY) / 3; - final double fPrime2StepIntegral = dy * dy / dx; - - final double x = currentX - startX; - f2Integral += f2StepIntegral; - fPrime2Integral += fPrime2StepIntegral; - - sx2 += x * x; - sy2 += f2Integral * f2Integral; - sxy += x * f2Integral; - sxz += x * fPrime2Integral; - syz += f2Integral * fPrime2Integral; - - } - - // 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 OptimizationException(LocalizedFormats.UNABLE_TO_FIRST_GUESS_HARMONIC_COEFFICIENTS); - } - a = FastMath.sqrt(c1 / c2); - omega = FastMath.sqrt(c2 / c3); - - } - - /** Estimate a first guess of the φ coefficient. - */ - private void guessPhi() { - - // initialize the means - double fcMean = 0.0; - double fsMean = 0.0; - - double currentX = observations[0].getX(); - double currentY = observations[0].getY(); - for (int i = 1; i < observations.length; ++i) { - - // one step forward - final double previousX = currentX; - final double previousY = currentY; - currentX = observations[i].getX(); - currentY = observations[i].getY(); - final double currentYPrime = (currentY - previousY) / (currentX - previousX); - - double omegaX = omega * currentX; - double cosine = FastMath.cos(omegaX); - double sine = FastMath.sin(omegaX); - fcMean += omega * currentY * cosine - currentYPrime * sine; - fsMean += omega * currentY * sine + currentYPrime * cosine; - - } - - phi = FastMath.atan2(-fsMean, fcMean); - - } - - /** Get the guessed amplitude a. - * @return guessed amplitude a; - */ - public double getGuessedAmplitude() { - return a; - } - - /** Get the guessed pulsation ω. - * @return guessed pulsation ω - */ - public double getGuessedPulsation() { - return omega; - } - - /** Get the guessed phase φ. - * @return guessed phase φ - */ - public double getGuessedPhase() { - return phi; - } - -} diff --git a/src/main/java/org/apache/commons/math/optimization/fitting/HarmonicFitter.java b/src/main/java/org/apache/commons/math/optimization/fitting/HarmonicFitter.java index 02dd08efc..41261ac21 100644 --- a/src/main/java/org/apache/commons/math/optimization/fitting/HarmonicFitter.java +++ b/src/main/java/org/apache/commons/math/optimization/fitting/HarmonicFitter.java @@ -17,113 +17,327 @@ package org.apache.commons.math.optimization.fitting; +import org.apache.commons.math.optimization.DifferentiableMultivariateVectorialOptimizer; +import org.apache.commons.math.analysis.function.HarmonicOscillator; +import org.apache.commons.math.exception.ZeroException; import org.apache.commons.math.exception.NumberIsTooSmallException; import org.apache.commons.math.exception.util.LocalizedFormats; -import org.apache.commons.math.optimization.DifferentiableMultivariateVectorialOptimizer; -import org.apache.commons.math.optimization.OptimizationException; -import org.apache.commons.math.analysis.ParametricUnivariateRealFunction; import org.apache.commons.math.util.FastMath; -/** This class implements a curve fitting specialized for sinusoids. - *Harmonic fitting is a very simple case of curve fitting. The
+/**
+ * Class that implements a curve fitting specialized for sinusoids.
+ *
+ * Harmonic fitting is a very simple case of curve fitting. The
* estimated coefficients are the amplitude a, the pulsation ω and
* the phase φ: f (t) = a cos (ω t + φ). They are
* searched by a least square estimator initialized with a rough guess
- * based on integrals.
This constructor can be used when a first guess of the - * coefficients is already known.
- * @param optimizer optimizer to use for the fitting - * @param initialGuess guessed values for amplitude (index 0), - * pulsation ω (index 1) and phase φ (index 2) - */ - public HarmonicFitter(final DifferentiableMultivariateVectorialOptimizer optimizer, - final double[] initialGuess) { - this.fitter = new CurveFitter(optimizer); - this.parameters = initialGuess.clone(); - } - - /** Add an observed weighted (x,y) point to the sample. - * @param weight weight of the observed point in the fit - * @param x abscissa of the point - * @param y observed value of the point at x, after fitting we should - * have P(x) as close as possible to this value - */ - public void addObservedPoint(double weight, double x, double y) { - fitter.addObservedPoint(weight, x, y); + super(optimizer); } /** * Fit an harmonic function to the observed points. * - * @return harmonic Function that best fits the observed points. - * @throws NumberIsTooSmallException if the sample is too short or if - * the first guess cannot be computed. - * @throws OptimizationException + * @param initialGuess First guess values in the following order: + *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.
+ */ + public static class ParameterGuesser { + /** Sampled observations. */ + private final WeightedObservedPoint[] observations; + /** Amplitude. */ + private double a; + /** Angular frequency. */ + private double omega; + /** Phase. */ + private double phi; + + /** + * Simple constructor. + * @param observations sampled observations + * @throws NumberIsTooSmallException if the sample is too short or if + * the first guess cannot be computed. + */ + public ParameterGuesser(WeightedObservedPoint[] observations) { if (observations.length < 4) { throw new NumberIsTooSmallException(LocalizedFormats.INSUFFICIENT_OBSERVED_POINTS_IN_SAMPLE, observations.length, 4, true); } - HarmonicCoefficientsGuesser guesser = new HarmonicCoefficientsGuesser(observations); - guesser.guess(); - parameters = new double[] { - guesser.getGuessedAmplitude(), - guesser.getGuessedPulsation(), - guesser.getGuessedPhase() - }; + this.observations = observations.clone(); } - double[] fitted = fitter.fit(new ParametricHarmonicFunction(), parameters); - return new HarmonicFunction(fitted[0], fitted[1], fitted[2]); + /** + * Estimate a first guess of the coefficients. + * + * @return the guessed coefficients, in the following order: + *