1/* 2 * Licensed to the Apache Software Foundation (ASF) under one or more 3 * contributor license agreements. See the NOTICE file distributed with 4 * this work for additional information regarding copyright ownership. 5 * The ASF licenses this file to You under the Apache License, Version 2.0 6 * (the "License"); you may not use this file except in compliance with 7 * the License. You may obtain a copy of the License at 8 * 9 * http://www.apache.org/licenses/LICENSE-2.0 10 * 11 * Unless required by applicable law or agreed to in writing, software 12 * distributed under the License is distributed on an "AS IS" BASIS, 13 * WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied. 14 * See the License for the specific language governing permissions and 15 * limitations under the License. 16 */ 17 18package org.apache.commons.math.optimization.fitting; 19 20import org.apache.commons.math.FunctionEvaluationException; 21import org.apache.commons.math.exception.util.LocalizedFormats; 22import org.apache.commons.math.optimization.DifferentiableMultivariateVectorialOptimizer; 23import org.apache.commons.math.optimization.OptimizationException; 24import org.apache.commons.math.util.FastMath; 25 26/** This class implements a curve fitting specialized for sinusoids. 27 * <p>Harmonic fitting is a very simple case of curve fitting. The 28 * estimated coefficients are the amplitude a, the pulsation ω and 29 * the phase φ: <code>f (t) = a cos (ω t + φ)</code>. They are 30 * searched by a least square estimator initialized with a rough guess 31 * based on integrals.</p> 32 * @version $Revision: 1073158 $ $Date: 2011-02-21 22:46:52 +0100 (lun. 21 févr. 2011) $ 33 * @since 2.0 34 */ 35public class HarmonicFitter { 36 37 /** Fitter for the coefficients. */ 38 private final CurveFitter fitter; 39 40 /** Values for amplitude, pulsation ω and phase φ. */ 41 private double[] parameters; 42 43 /** Simple constructor. 44 * @param optimizer optimizer to use for the fitting 45 */ 46 public HarmonicFitter(final DifferentiableMultivariateVectorialOptimizer optimizer) { 47 this.fitter = new CurveFitter(optimizer); 48 parameters = null; 49 } 50 51 /** Simple constructor. 52 * <p>This constructor can be used when a first guess of the 53 * coefficients is already known.</p> 54 * @param optimizer optimizer to use for the fitting 55 * @param initialGuess guessed values for amplitude (index 0), 56 * pulsation ω (index 1) and phase φ (index 2) 57 */ 58 public HarmonicFitter(final DifferentiableMultivariateVectorialOptimizer optimizer, 59 final double[] initialGuess) { 60 this.fitter = new CurveFitter(optimizer); 61 this.parameters = initialGuess.clone(); 62 } 63 64 /** Add an observed weighted (x,y) point to the sample. 65 * @param weight weight of the observed point in the fit 66 * @param x abscissa of the point 67 * @param y observed value of the point at x, after fitting we should 68 * have P(x) as close as possible to this value 69 */ 70 public void addObservedPoint(double weight, double x, double y) { 71 fitter.addObservedPoint(weight, x, y); 72 } 73 74 /** Fit an harmonic function to the observed points. 75 * @return harmonic function best fitting the observed points 76 * @throws OptimizationException if the sample is too short or if 77 * the first guess cannot be computed 78 */ 79 public HarmonicFunction fit() throws OptimizationException { 80 81 // shall we compute the first guess of the parameters ourselves ? 82 if (parameters == null) { 83 final WeightedObservedPoint[] observations = fitter.getObservations(); 84 if (observations.length < 4) { 85 throw new OptimizationException(LocalizedFormats.INSUFFICIENT_OBSERVED_POINTS_IN_SAMPLE, 86 observations.length, 4); 87 } 88 89 HarmonicCoefficientsGuesser guesser = new HarmonicCoefficientsGuesser(observations); 90 guesser.guess(); 91 parameters = new double[] { 92 guesser.getGuessedAmplitude(), 93 guesser.getGuessedPulsation(), 94 guesser.getGuessedPhase() 95 }; 96 97 } 98 99 try { 100 double[] fitted = fitter.fit(new ParametricHarmonicFunction(), parameters); 101 return new HarmonicFunction(fitted[0], fitted[1], fitted[2]); 102 } catch (FunctionEvaluationException fee) { 103 // should never happen 104 throw new RuntimeException(fee); 105 } 106 107 } 108 109 /** Parametric harmonic function. */ 110 private static class ParametricHarmonicFunction implements ParametricRealFunction { 111 112 /** {@inheritDoc} */ 113 public double value(double x, double[] parameters) { 114 final double a = parameters[0]; 115 final double omega = parameters[1]; 116 final double phi = parameters[2]; 117 return a * FastMath.cos(omega * x + phi); 118 } 119 120 /** {@inheritDoc} */ 121 public double[] gradient(double x, double[] parameters) { 122 final double a = parameters[0]; 123 final double omega = parameters[1]; 124 final double phi = parameters[2]; 125 final double alpha = omega * x + phi; 126 final double cosAlpha = FastMath.cos(alpha); 127 final double sinAlpha = FastMath.sin(alpha); 128 return new double[] { cosAlpha, -a * x * sinAlpha, -a * sinAlpha }; 129 } 130 131 } 132 133} 134