diff options
Diffstat (limited to 'src/main/java/org/apache/commons/math/optimization/fitting/HarmonicCoefficientsGuesser.java')
| -rw-r--r-- | src/main/java/org/apache/commons/math/optimization/fitting/HarmonicCoefficientsGuesser.java | 300 |
1 files changed, 300 insertions, 0 deletions
diff --git a/src/main/java/org/apache/commons/math/optimization/fitting/HarmonicCoefficientsGuesser.java b/src/main/java/org/apache/commons/math/optimization/fitting/HarmonicCoefficientsGuesser.java new file mode 100644 index 0000000..0796f67 --- /dev/null +++ b/src/main/java/org/apache/commons/math/optimization/fitting/HarmonicCoefficientsGuesser.java @@ -0,0 +1,300 @@ +/* + * Licensed to the Apache Software Foundation (ASF) under one or more + * contributor license agreements. See the NOTICE file distributed with + * this work for additional information regarding copyright ownership. + * The ASF licenses this file to You under the Apache License, Version 2.0 + * (the "License"); you may not use this file except in compliance with + * the License. You may obtain a copy of the License at + * + * http://www.apache.org/licenses/LICENSE-2.0 + * + * Unless required by applicable law or agreed to in writing, software + * distributed under the License is distributed on an "AS IS" BASIS, + * WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied. + * See the License for the specific language governing permissions and + * limitations under the License. + */ + +package org.apache.commons.math.optimization.fitting; + +import org.apache.commons.math.exception.util.LocalizedFormats; +import org.apache.commons.math.optimization.OptimizationException; +import org.apache.commons.math.util.FastMath; + +/** This class guesses harmonic coefficients from a sample. + + * <p>The algorithm used to guess the coefficients is as follows:</p> + + * <p>We know f (t) at some sampling points t<sub>i</sub> and want to find a, + * ω and φ such that f (t) = a cos (ω t + φ). + * </p> + * + * <p>From the analytical expression, we can compute two primitives : + * <pre> + * If2 (t) = ∫ f<sup>2</sup> = a<sup>2</sup> × [t + S (t)] / 2 + * If'2 (t) = ∫ f'<sup>2</sup> = a<sup>2</sup> ω<sup>2</sup> × [t - S (t)] / 2 + * where S (t) = sin (2 (ω t + φ)) / (2 ω) + * </pre> + * </p> + * + * <p>We can remove S between these expressions : + * <pre> + * If'2 (t) = a<sup>2</sup> ω<sup>2</sup> t - ω<sup>2</sup> If2 (t) + * </pre> + * </p> + * + * <p>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) + * </p> + * + * <p>From the primitive, we can deduce the same form for definite + * integrals between t<sub>1</sub> and t<sub>i</sub> for each t<sub>i</sub> : + * <pre> + * If2 (t<sub>i</sub>) - If2 (t<sub>1</sub>) = A × (t<sub>i</sub> - t<sub>1</sub>) + B × (If2 (t<sub>i</sub>) - If2 (t<sub>1</sub>)) + * </pre> + * </p> + * + * <p>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. + * </p> + * + * <p>For a bilinear expression z (x<sub>i</sub>, y<sub>i</sub>) = A × x<sub>i</sub> + B × y<sub>i</sub>, the + * coefficients A and B that minimize a least square criterion + * ∑ (z<sub>i</sub> - z (x<sub>i</sub>, y<sub>i</sub>))<sup>2</sup> are given by these expressions:</p> + * <pre> + * + * ∑y<sub>i</sub>y<sub>i</sub> ∑x<sub>i</sub>z<sub>i</sub> - ∑x<sub>i</sub>y<sub>i</sub> ∑y<sub>i</sub>z<sub>i</sub> + * A = ------------------------ + * ∑x<sub>i</sub>x<sub>i</sub> ∑y<sub>i</sub>y<sub>i</sub> - ∑x<sub>i</sub>y<sub>i</sub> ∑x<sub>i</sub>y<sub>i</sub> + * + * ∑x<sub>i</sub>x<sub>i</sub> ∑y<sub>i</sub>z<sub>i</sub> - ∑x<sub>i</sub>y<sub>i</sub> ∑x<sub>i</sub>z<sub>i</sub> + * B = ------------------------ + * ∑x<sub>i</sub>x<sub>i</sub> ∑y<sub>i</sub>y<sub>i</sub> - ∑x<sub>i</sub>y<sub>i</sub> ∑x<sub>i</sub>y<sub>i</sub> + * </pre> + * </p> + * + * + * <p>In fact, we can assume both a and ω are positive and + * compute them directly, knowing that A = a<sup>2</sup> ω<sup>2</sup> and that + * B = - ω<sup>2</sup>. The complete algorithm is therefore:</p> + * <pre> + * + * for each t<sub>i</sub> from t<sub>1</sub> to t<sub>n-1</sub>, compute: + * f (t<sub>i</sub>) + * f' (t<sub>i</sub>) = (f (t<sub>i+1</sub>) - f(t<sub>i-1</sub>)) / (t<sub>i+1</sub> - t<sub>i-1</sub>) + * x<sub>i</sub> = t<sub>i</sub> - t<sub>1</sub> + * y<sub>i</sub> = ∫ f<sup>2</sup> from t<sub>1</sub> to t<sub>i</sub> + * z<sub>i</sub> = ∫ f'<sup>2</sup> from t<sub>1</sub> to t<sub>i</sub> + * update the sums ∑x<sub>i</sub>x<sub>i</sub>, ∑y<sub>i</sub>y<sub>i</sub>, ∑x<sub>i</sub>y<sub>i</sub>, ∑x<sub>i</sub>z<sub>i</sub> and ∑y<sub>i</sub>z<sub>i</sub> + * end for + * + * |-------------------------- + * \ | ∑y<sub>i</sub>y<sub>i</sub> ∑x<sub>i</sub>z<sub>i</sub> - ∑x<sub>i</sub>y<sub>i</sub> ∑y<sub>i</sub>z<sub>i</sub> + * a = \ | ------------------------ + * \| ∑x<sub>i</sub>y<sub>i</sub> ∑x<sub>i</sub>z<sub>i</sub> - ∑x<sub>i</sub>x<sub>i</sub> ∑y<sub>i</sub>z<sub>i</sub> + * + * + * |-------------------------- + * \ | ∑x<sub>i</sub>y<sub>i</sub> ∑x<sub>i</sub>z<sub>i</sub> - ∑x<sub>i</sub>x<sub>i</sub> ∑y<sub>i</sub>z<sub>i</sub> + * ω = \ | ------------------------ + * \| ∑x<sub>i</sub>x<sub>i</sub> ∑y<sub>i</sub>y<sub>i</sub> - ∑x<sub>i</sub>y<sub>i</sub> ∑x<sub>i</sub>y<sub>i</sub> + * + * </pre> + * </p> + + * <p>Once we know ω, we can compute: + * <pre> + * fc = ω f (t) cos (ω t) - f' (t) sin (ω t) + * fs = ω f (t) sin (ω t) + f' (t) cos (ω t) + * </pre> + * </p> + + * <p>It appears that <code>fc = a ω cos (φ)</code> and + * <code>fs = -a ω sin (φ)</code>, so we can use these + * expressions to compute φ. The best estimate over the sample is + * given by averaging these expressions. + * </p> + + * <p>Since integrals and means are involved in the preceding + * estimations, these operations run in O(n) time, where n is the + * number of measurements.</p> + + * @version $Revision: 1056034 $ $Date: 2011-01-06 20:41:43 +0100 (jeu. 06 janv. 2011) $ + * @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 f<sup>2</sup> and f'<sup>2</sup> + // 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; + } + +} |