/* * 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. * *

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. * * @since 2.0 * @deprecated As of 3.3. Please use {@link HarmonicCurveFitter} and {@link WeightedObservedPoints} * instead. */ @Deprecated public class HarmonicFitter extends CurveFitter { /** * 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: *

* * @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. * *

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 { /** 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: *

*/ 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 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) || (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); } } }