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.MathRuntimeException; 022 import org.apache.commons.math.optimization.DifferentiableMultivariateVectorialOptimizer; 023 import org.apache.commons.math.optimization.OptimizationException; 024 025 /** This class implements a curve fitting specialized for sinusoids. 026 * <p>Harmonic fitting is a very simple case of curve fitting. The 027 * estimated coefficients are the amplitude a, the pulsation ω and 028 * the phase φ: <code>f (t) = a cos (ω t + φ)</code>. They are 029 * searched by a least square estimator initialized with a rough guess 030 * based on integrals.</p> 031 * @version $Revision: 786479 $ $Date: 2009-06-19 08:36:16 -0400 (Fri, 19 Jun 2009) $ 032 * @since 2.0 033 */ 034 public class HarmonicFitter { 035 036 /** Fitter for the coefficients. */ 037 private final CurveFitter fitter; 038 039 /** Values for amplitude, pulsation ω and phase φ. */ 040 private double[] parameters; 041 042 /** Simple constructor. 043 * @param optimizer optimizer to use for the fitting 044 */ 045 public HarmonicFitter(final DifferentiableMultivariateVectorialOptimizer optimizer) { 046 this.fitter = new CurveFitter(optimizer); 047 parameters = null; 048 } 049 050 /** Simple constructor. 051 * <p>This constructor can be used when a first guess of the 052 * coefficients is already known.</p> 053 * @param optimizer optimizer to use for the fitting 054 * @param initialGuess guessed values for amplitude (index 0), 055 * pulsation ω (index 1) and phase φ (index 2) 056 */ 057 public HarmonicFitter(final DifferentiableMultivariateVectorialOptimizer optimizer, 058 final double[] initialGuess) { 059 this.fitter = new CurveFitter(optimizer); 060 this.parameters = initialGuess.clone(); 061 } 062 063 /** Add an observed weighted (x,y) point to the sample. 064 * @param weight weight of the observed point in the fit 065 * @param x abscissa of the point 066 * @param y observed value of the point at x, after fitting we should 067 * have P(x) as close as possible to this value 068 */ 069 public void addObservedPoint(double weight, double x, double y) { 070 fitter.addObservedPoint(weight, x, y); 071 } 072 073 /** Fit an harmonic function to the observed points. 074 * @return harmonic function best fitting the observed points 075 * @throws OptimizationException if the sample is too short or if 076 * the first guess cannot be computed 077 */ 078 public HarmonicFunction fit() throws OptimizationException { 079 try { 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("sample contains {0} observed points, at least {1} are required", 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 double[] fitted = fitter.fit(new ParametricHarmonicFunction(), parameters); 100 return new HarmonicFunction(fitted[0], fitted[1], fitted[2]); 101 102 } catch (FunctionEvaluationException fee) { 103 // this should never happen 104 throw MathRuntimeException.createInternalError(fee); 105 } 106 } 107 108 /** Parametric harmonic function. */ 109 private static class ParametricHarmonicFunction implements ParametricRealFunction { 110 111 /** {@inheritDoc} */ 112 public double value(double x, double[] parameters) { 113 final double a = parameters[0]; 114 final double omega = parameters[1]; 115 final double phi = parameters[2]; 116 return a * Math.cos(omega * x + phi); 117 } 118 119 /** {@inheritDoc} */ 120 public double[] gradient(double x, double[] parameters) { 121 final double a = parameters[0]; 122 final double omega = parameters[1]; 123 final double phi = parameters[2]; 124 final double alpha = omega * x + phi; 125 final double cosAlpha = Math.cos(alpha); 126 final double sinAlpha = Math.sin(alpha); 127 return new double[] { cosAlpha, -a * x * sinAlpha, -a * sinAlpha }; 128 } 129 130 } 131 132 }