diff options
Diffstat (limited to 'src/main/java/org/apache/commons/math3/fitting/HarmonicFitter.java')
| -rw-r--r-- | src/main/java/org/apache/commons/math3/fitting/HarmonicFitter.java | 386 |
1 files changed, 386 insertions, 0 deletions
diff --git a/src/main/java/org/apache/commons/math3/fitting/HarmonicFitter.java b/src/main/java/org/apache/commons/math3/fitting/HarmonicFitter.java new file mode 100644 index 0000000..1a41398 --- /dev/null +++ b/src/main/java/org/apache/commons/math3/fitting/HarmonicFitter.java @@ -0,0 +1,386 @@ +/* + * 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.math3.fitting; + +import org.apache.commons.math3.analysis.function.HarmonicOscillator; +import org.apache.commons.math3.exception.MathIllegalStateException; +import org.apache.commons.math3.exception.NumberIsTooSmallException; +import org.apache.commons.math3.exception.ZeroException; +import org.apache.commons.math3.exception.util.LocalizedFormats; +import org.apache.commons.math3.optim.nonlinear.vector.MultivariateVectorOptimizer; +import org.apache.commons.math3.util.FastMath; + +/** + * Class that implements a curve fitting specialized for sinusoids. + * + * <p>Harmonic fitting is a very simple case of curve fitting. The estimated coefficients are the + * amplitude a, the pulsation ω and the phase φ: <code>f (t) = a cos (ω t + φ) + * </code>. They are searched by a least square estimator initialized with a rough guess based on + * integrals. + * + * @since 2.0 + * @deprecated As of 3.3. Please use {@link HarmonicCurveFitter} and {@link WeightedObservedPoints} + * instead. + */ +@Deprecated +public class HarmonicFitter extends CurveFitter<HarmonicOscillator.Parametric> { + /** + * Simple constructor. + * + * @param optimizer Optimizer to use for the fitting. + */ + public HarmonicFitter(final MultivariateVectorOptimizer optimizer) { + super(optimizer); + } + + /** + * Fit an harmonic function to the observed points. + * + * @param initialGuess First guess values in the following order: + * <ul> + * <li>Amplitude + * <li>Angular frequency + * <li>Phase + * </ul> + * + * @return the parameters of the harmonic function that best fits the observed points (in the + * same order as above). + */ + public double[] fit(double[] initialGuess) { + return fit(new HarmonicOscillator.Parametric(), initialGuess); + } + + /** + * Fit an harmonic function to the observed points. An initial guess will be automatically + * computed. + * + * @return the parameters of the harmonic function that best fits the observed points (see the + * other {@link #fit(double[]) fit} method. + * @throws NumberIsTooSmallException if the sample is too short for the the first guess to be + * computed. + * @throws ZeroException if the first guess cannot be computed because the abscissa range is + * zero. + */ + public double[] fit() { + return fit((new ParameterGuesser(getObservations())).guess()); + } + + /** + * This class guesses harmonic coefficients from a sample. + * + * <p>The algorithm used to guess the coefficients is as follows: + * + * <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>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>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>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>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>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>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: + * + * <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>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: + * + * <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>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>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>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 { + /** Amplitude. */ + private final double a; + + /** Angular frequency. */ + private final double omega; + + /** Phase. */ + private final double phi; + + /** + * Simple constructor. + * + * @param observations Sampled observations. + * @throws NumberIsTooSmallException if the sample is too short. + * @throws ZeroException if the abscissa range is zero. + * @throws MathIllegalStateException when the guessing procedure cannot produce sensible + * results. + */ + public ParameterGuesser(WeightedObservedPoint[] observations) { + if (observations.length < 4) { + throw new NumberIsTooSmallException( + LocalizedFormats.INSUFFICIENT_OBSERVED_POINTS_IN_SAMPLE, + observations.length, + 4, + true); + } + + final WeightedObservedPoint[] sorted = sortObservations(observations); + + final double aOmega[] = guessAOmega(sorted); + a = aOmega[0]; + omega = aOmega[1]; + + phi = guessPhi(sorted); + } + + /** + * Gets an estimation of the parameters. + * + * @return the guessed parameters, in the following order: + * <ul> + * <li>Amplitude + * <li>Angular frequency + * <li>Phase + * </ul> + */ + public double[] guess() { + return new double[] {a, omega, phi}; + } + + /** + * Sort the observations with respect to the abscissa. + * + * @param unsorted Input observations. + * @return the input observations, sorted. + */ + private WeightedObservedPoint[] sortObservations(WeightedObservedPoint[] unsorted) { + final WeightedObservedPoint[] observations = unsorted.clone(); + + // 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]; + } + } + + return observations; + } + + /** + * Estimate a first guess of the amplitude and angular frequency. This method assumes that + * the {@link #sortObservations(WeightedObservedPoint[])} method has been called previously. + * + * @param observations Observations, sorted w.r.t. abscissa. + * @throws ZeroException if the abscissa range is zero. + * @throws MathIllegalStateException when the guessing procedure cannot produce sensible + * results. + * @return the guessed amplitude (at index 0) and circular frequency (at index 1). + */ + private double[] guessAOmega(WeightedObservedPoint[] observations) { + final double[] aOmega = new double[2]; + + // initialize the sums for the linear model between the two integrals + double sx2 = 0; + double sy2 = 0; + double sxy = 0; + double sxz = 0; + double syz = 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) || (c2 / c3 < 0)) { + final int last = observations.length - 1; + // Range of the observations, assuming that the + // observations are sorted. + final double xRange = observations[last].getX() - observations[0].getX(); + if (xRange == 0) { + throw new ZeroException(); + } + aOmega[1] = 2 * Math.PI / xRange; + + double yMin = Double.POSITIVE_INFINITY; + double yMax = Double.NEGATIVE_INFINITY; + for (int i = 1; i < observations.length; ++i) { + final double y = observations[i].getY(); + if (y < yMin) { + yMin = y; + } + if (y > yMax) { + yMax = y; + } + } + aOmega[0] = 0.5 * (yMax - yMin); + } else { + if (c2 == 0) { + // In some ill-conditioned cases (cf. MATH-844), the guesser + // procedure cannot produce sensible results. + throw new MathIllegalStateException(LocalizedFormats.ZERO_DENOMINATOR); + } + + aOmega[0] = FastMath.sqrt(c1 / c2); + aOmega[1] = FastMath.sqrt(c2 / c3); + } + + return aOmega; + } + + /** + * Estimate a first guess of the phase. + * + * @param observations Observations, sorted w.r.t. abscissa. + * @return the guessed phase. + */ + private double guessPhi(WeightedObservedPoint[] observations) { + // initialize the means + double fcMean = 0; + double fsMean = 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; + } + + return FastMath.atan2(-fsMean, fcMean); + } + } +} |