001 /* 002 * Licensed to the Apache Software Foundation (ASF) under one or more 003 * contributor license agreements. See the NOTICE file distributed with 004 * this work for additional information regarding copyright ownership. 005 * The ASF licenses this file to You under the Apache License, Version 2.0 006 * (the "License"); you may not use this file except in compliance with 007 * the License. You may obtain a copy of the License at 008 * 009 * http://www.apache.org/licenses/LICENSE-2.0 010 * 011 * Unless required by applicable law or agreed to in writing, software 012 * distributed under the License is distributed on an "AS IS" BASIS, 013 * WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied. 014 * See the License for the specific language governing permissions and 015 * limitations under the License. 016 */ 017 018 package org.apache.commons.math.optimization.fitting; 019 020 import org.apache.commons.math.FunctionEvaluationException; 021 import org.apache.commons.math.exception.util.LocalizedFormats; 022 import org.apache.commons.math.optimization.DifferentiableMultivariateVectorialOptimizer; 023 import org.apache.commons.math.optimization.OptimizationException; 024 import org.apache.commons.math.util.FastMath; 025 026 /** This class implements a curve fitting specialized for sinusoids. 027 * <p>Harmonic fitting is a very simple case of curve fitting. The 028 * estimated coefficients are the amplitude a, the pulsation ω and 029 * the phase φ: <code>f (t) = a cos (ω t + φ)</code>. They are 030 * searched by a least square estimator initialized with a rough guess 031 * based on integrals.</p> 032 * @version $Revision: 1073158 $ $Date: 2011-02-21 22:46:52 +0100 (lun. 21 f??vr. 2011) $ 033 * @since 2.0 034 */ 035 public class HarmonicFitter { 036 037 /** Fitter for the coefficients. */ 038 private final CurveFitter fitter; 039 040 /** Values for amplitude, pulsation ω and phase φ. */ 041 private double[] parameters; 042 043 /** Simple constructor. 044 * @param optimizer optimizer to use for the fitting 045 */ 046 public HarmonicFitter(final DifferentiableMultivariateVectorialOptimizer optimizer) { 047 this.fitter = new CurveFitter(optimizer); 048 parameters = null; 049 } 050 051 /** Simple constructor. 052 * <p>This constructor can be used when a first guess of the 053 * coefficients is already known.</p> 054 * @param optimizer optimizer to use for the fitting 055 * @param initialGuess guessed values for amplitude (index 0), 056 * pulsation ω (index 1) and phase φ (index 2) 057 */ 058 public HarmonicFitter(final DifferentiableMultivariateVectorialOptimizer optimizer, 059 final double[] initialGuess) { 060 this.fitter = new CurveFitter(optimizer); 061 this.parameters = initialGuess.clone(); 062 } 063 064 /** Add an observed weighted (x,y) point to the sample. 065 * @param weight weight of the observed point in the fit 066 * @param x abscissa of the point 067 * @param y observed value of the point at x, after fitting we should 068 * have P(x) as close as possible to this value 069 */ 070 public void addObservedPoint(double weight, double x, double y) { 071 fitter.addObservedPoint(weight, x, y); 072 } 073 074 /** Fit an harmonic function to the observed points. 075 * @return harmonic function best fitting the observed points 076 * @throws OptimizationException if the sample is too short or if 077 * the first guess cannot be computed 078 */ 079 public HarmonicFunction fit() throws OptimizationException { 080 081 // shall we compute the first guess of the parameters ourselves ? 082 if (parameters == null) { 083 final WeightedObservedPoint[] observations = fitter.getObservations(); 084 if (observations.length < 4) { 085 throw new OptimizationException(LocalizedFormats.INSUFFICIENT_OBSERVED_POINTS_IN_SAMPLE, 086 observations.length, 4); 087 } 088 089 HarmonicCoefficientsGuesser guesser = new HarmonicCoefficientsGuesser(observations); 090 guesser.guess(); 091 parameters = new double[] { 092 guesser.getGuessedAmplitude(), 093 guesser.getGuessedPulsation(), 094 guesser.getGuessedPhase() 095 }; 096 097 } 098 099 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 }