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.exception.util.LocalizedFormats; 021 import org.apache.commons.math.optimization.OptimizationException; 022 import org.apache.commons.math.util.FastMath; 023 024 /** This class guesses harmonic coefficients from a sample. 025 026 * <p>The algorithm used to guess the coefficients is as follows:</p> 027 028 * <p>We know f (t) at some sampling points t<sub>i</sub> and want to find a, 029 * ω and φ such that f (t) = a cos (ω t + φ). 030 * </p> 031 * 032 * <p>From the analytical expression, we can compute two primitives : 033 * <pre> 034 * If2 (t) = ∫ f<sup>2</sup> = a<sup>2</sup> × [t + S (t)] / 2 035 * If'2 (t) = ∫ f'<sup>2</sup> = a<sup>2</sup> ω<sup>2</sup> × [t - S (t)] / 2 036 * where S (t) = sin (2 (ω t + φ)) / (2 ω) 037 * </pre> 038 * </p> 039 * 040 * <p>We can remove S between these expressions : 041 * <pre> 042 * If'2 (t) = a<sup>2</sup> ω<sup>2</sup> t - ω<sup>2</sup> If2 (t) 043 * </pre> 044 * </p> 045 * 046 * <p>The preceding expression shows that If'2 (t) is a linear 047 * combination of both t and If2 (t): If'2 (t) = A × t + B × If2 (t) 048 * </p> 049 * 050 * <p>From the primitive, we can deduce the same form for definite 051 * integrals between t<sub>1</sub> and t<sub>i</sub> for each t<sub>i</sub> : 052 * <pre> 053 * 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>)) 054 * </pre> 055 * </p> 056 * 057 * <p>We can find the coefficients A and B that best fit the sample 058 * to this linear expression by computing the definite integrals for 059 * each sample points. 060 * </p> 061 * 062 * <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 063 * coefficients A and B that minimize a least square criterion 064 * ∑ (z<sub>i</sub> - z (x<sub>i</sub>, y<sub>i</sub>))<sup>2</sup> are given by these expressions:</p> 065 * <pre> 066 * 067 * ∑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> 068 * A = ------------------------ 069 * ∑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> 070 * 071 * ∑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> 072 * B = ------------------------ 073 * ∑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> 074 * </pre> 075 * </p> 076 * 077 * 078 * <p>In fact, we can assume both a and ω are positive and 079 * compute them directly, knowing that A = a<sup>2</sup> ω<sup>2</sup> and that 080 * B = - ω<sup>2</sup>. The complete algorithm is therefore:</p> 081 * <pre> 082 * 083 * for each t<sub>i</sub> from t<sub>1</sub> to t<sub>n-1</sub>, compute: 084 * f (t<sub>i</sub>) 085 * 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>) 086 * x<sub>i</sub> = t<sub>i</sub> - t<sub>1</sub> 087 * y<sub>i</sub> = ∫ f<sup>2</sup> from t<sub>1</sub> to t<sub>i</sub> 088 * z<sub>i</sub> = ∫ f'<sup>2</sup> from t<sub>1</sub> to t<sub>i</sub> 089 * 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> 090 * end for 091 * 092 * |-------------------------- 093 * \ | ∑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> 094 * a = \ | ------------------------ 095 * \| ∑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> 096 * 097 * 098 * |-------------------------- 099 * \ | ∑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> 100 * ω = \ | ------------------------ 101 * \| ∑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> 102 * 103 * </pre> 104 * </p> 105 106 * <p>Once we know ω, we can compute: 107 * <pre> 108 * fc = ω f (t) cos (ω t) - f' (t) sin (ω t) 109 * fs = ω f (t) sin (ω t) + f' (t) cos (ω t) 110 * </pre> 111 * </p> 112 113 * <p>It appears that <code>fc = a ω cos (φ)</code> and 114 * <code>fs = -a ω sin (φ)</code>, so we can use these 115 * expressions to compute φ. The best estimate over the sample is 116 * given by averaging these expressions. 117 * </p> 118 119 * <p>Since integrals and means are involved in the preceding 120 * estimations, these operations run in O(n) time, where n is the 121 * number of measurements.</p> 122 123 * @version $Revision: 1056034 $ $Date: 2011-01-06 20:41:43 +0100 (jeu. 06 janv. 2011) $ 124 * @since 2.0 125 126 */ 127 public class HarmonicCoefficientsGuesser { 128 129 /** Sampled observations. */ 130 private final WeightedObservedPoint[] observations; 131 132 /** Guessed amplitude. */ 133 private double a; 134 135 /** Guessed pulsation ω. */ 136 private double omega; 137 138 /** Guessed phase φ. */ 139 private double phi; 140 141 /** Simple constructor. 142 * @param observations sampled observations 143 */ 144 public HarmonicCoefficientsGuesser(WeightedObservedPoint[] observations) { 145 this.observations = observations.clone(); 146 a = Double.NaN; 147 omega = Double.NaN; 148 } 149 150 /** Estimate a first guess of the coefficients. 151 * @exception OptimizationException if the sample is too short or if 152 * the first guess cannot be computed (when the elements under the 153 * square roots are negative). 154 * */ 155 public void guess() throws OptimizationException { 156 sortObservations(); 157 guessAOmega(); 158 guessPhi(); 159 } 160 161 /** Sort the observations with respect to the abscissa. 162 */ 163 private void sortObservations() { 164 165 // Since the samples are almost always already sorted, this 166 // method is implemented as an insertion sort that reorders the 167 // elements in place. Insertion sort is very efficient in this case. 168 WeightedObservedPoint curr = observations[0]; 169 for (int j = 1; j < observations.length; ++j) { 170 WeightedObservedPoint prec = curr; 171 curr = observations[j]; 172 if (curr.getX() < prec.getX()) { 173 // the current element should be inserted closer to the beginning 174 int i = j - 1; 175 WeightedObservedPoint mI = observations[i]; 176 while ((i >= 0) && (curr.getX() < mI.getX())) { 177 observations[i + 1] = mI; 178 if (i-- != 0) { 179 mI = observations[i]; 180 } 181 } 182 observations[i + 1] = curr; 183 curr = observations[j]; 184 } 185 } 186 187 } 188 189 /** Estimate a first guess of the a and ω coefficients. 190 * @exception OptimizationException if the sample is too short or if 191 * the first guess cannot be computed (when the elements under the 192 * square roots are negative). 193 */ 194 private void guessAOmega() throws OptimizationException { 195 196 // initialize the sums for the linear model between the two integrals 197 double sx2 = 0.0; 198 double sy2 = 0.0; 199 double sxy = 0.0; 200 double sxz = 0.0; 201 double syz = 0.0; 202 203 double currentX = observations[0].getX(); 204 double currentY = observations[0].getY(); 205 double f2Integral = 0; 206 double fPrime2Integral = 0; 207 final double startX = currentX; 208 for (int i = 1; i < observations.length; ++i) { 209 210 // one step forward 211 final double previousX = currentX; 212 final double previousY = currentY; 213 currentX = observations[i].getX(); 214 currentY = observations[i].getY(); 215 216 // update the integrals of f<sup>2</sup> and f'<sup>2</sup> 217 // considering a linear model for f (and therefore constant f') 218 final double dx = currentX - previousX; 219 final double dy = currentY - previousY; 220 final double f2StepIntegral = 221 dx * (previousY * previousY + previousY * currentY + currentY * currentY) / 3; 222 final double fPrime2StepIntegral = dy * dy / dx; 223 224 final double x = currentX - startX; 225 f2Integral += f2StepIntegral; 226 fPrime2Integral += fPrime2StepIntegral; 227 228 sx2 += x * x; 229 sy2 += f2Integral * f2Integral; 230 sxy += x * f2Integral; 231 sxz += x * fPrime2Integral; 232 syz += f2Integral * fPrime2Integral; 233 234 } 235 236 // compute the amplitude and pulsation coefficients 237 double c1 = sy2 * sxz - sxy * syz; 238 double c2 = sxy * sxz - sx2 * syz; 239 double c3 = sx2 * sy2 - sxy * sxy; 240 if ((c1 / c2 < 0.0) || (c2 / c3 < 0.0)) { 241 throw new OptimizationException(LocalizedFormats.UNABLE_TO_FIRST_GUESS_HARMONIC_COEFFICIENTS); 242 } 243 a = FastMath.sqrt(c1 / c2); 244 omega = FastMath.sqrt(c2 / c3); 245 246 } 247 248 /** Estimate a first guess of the φ coefficient. 249 */ 250 private void guessPhi() { 251 252 // initialize the means 253 double fcMean = 0.0; 254 double fsMean = 0.0; 255 256 double currentX = observations[0].getX(); 257 double currentY = observations[0].getY(); 258 for (int i = 1; i < observations.length; ++i) { 259 260 // one step forward 261 final double previousX = currentX; 262 final double previousY = currentY; 263 currentX = observations[i].getX(); 264 currentY = observations[i].getY(); 265 final double currentYPrime = (currentY - previousY) / (currentX - previousX); 266 267 double omegaX = omega * currentX; 268 double cosine = FastMath.cos(omegaX); 269 double sine = FastMath.sin(omegaX); 270 fcMean += omega * currentY * cosine - currentYPrime * sine; 271 fsMean += omega * currentY * sine + currentYPrime * cosine; 272 273 } 274 275 phi = FastMath.atan2(-fsMean, fcMean); 276 277 } 278 279 /** Get the guessed amplitude a. 280 * @return guessed amplitude a; 281 */ 282 public double getGuessedAmplitude() { 283 return a; 284 } 285 286 /** Get the guessed pulsation ω. 287 * @return guessed pulsation ω 288 */ 289 public double getGuessedPulsation() { 290 return omega; 291 } 292 293 /** Get the guessed phase φ. 294 * @return guessed phase φ 295 */ 296 public double getGuessedPhase() { 297 return phi; 298 } 299 300 }