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 18 package org.apache.commons.math.optimization.fitting; 19 20 import org.apache.commons.math.exception.util.LocalizedFormats; 21 import org.apache.commons.math.optimization.OptimizationException; 22 import org.apache.commons.math.util.FastMath; 23 24 /** This class guesses harmonic coefficients from a sample. 25 26 * <p>The algorithm used to guess the coefficients is as follows:</p> 27 28 * <p>We know f (t) at some sampling points t<sub>i</sub> and want to find a, 29 * ω and φ such that f (t) = a cos (ω t + φ). 30 * </p> 31 * 32 * <p>From the analytical expression, we can compute two primitives : 33 * <pre> 34 * If2 (t) = ∫ f<sup>2</sup> = a<sup>2</sup> × [t + S (t)] / 2 35 * If'2 (t) = ∫ f'<sup>2</sup> = a<sup>2</sup> ω<sup>2</sup> × [t - S (t)] / 2 36 * where S (t) = sin (2 (ω t + φ)) / (2 ω) 37 * </pre> 38 * </p> 39 * 40 * <p>We can remove S between these expressions : 41 * <pre> 42 * If'2 (t) = a<sup>2</sup> ω<sup>2</sup> t - ω<sup>2</sup> If2 (t) 43 * </pre> 44 * </p> 45 * 46 * <p>The preceding expression shows that If'2 (t) is a linear 47 * combination of both t and If2 (t): If'2 (t) = A × t + B × If2 (t) 48 * </p> 49 * 50 * <p>From the primitive, we can deduce the same form for definite 51 * integrals between t<sub>1</sub> and t<sub>i</sub> for each t<sub>i</sub> : 52 * <pre> 53 * 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>)) 54 * </pre> 55 * </p> 56 * 57 * <p>We can find the coefficients A and B that best fit the sample 58 * to this linear expression by computing the definite integrals for 59 * each sample points. 60 * </p> 61 * 62 * <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 63 * coefficients A and B that minimize a least square criterion 64 * ∑ (z<sub>i</sub> - z (x<sub>i</sub>, y<sub>i</sub>))<sup>2</sup> are given by these expressions:</p> 65 * <pre> 66 * 67 * ∑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> 68 * A = ------------------------ 69 * ∑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> 70 * 71 * ∑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> 72 * B = ------------------------ 73 * ∑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> 74 * </pre> 75 * </p> 76 * 77 * 78 * <p>In fact, we can assume both a and ω are positive and 79 * compute them directly, knowing that A = a<sup>2</sup> ω<sup>2</sup> and that 80 * B = - ω<sup>2</sup>. The complete algorithm is therefore:</p> 81 * <pre> 82 * 83 * for each t<sub>i</sub> from t<sub>1</sub> to t<sub>n-1</sub>, compute: 84 * f (t<sub>i</sub>) 85 * 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>) 86 * x<sub>i</sub> = t<sub>i</sub> - t<sub>1</sub> 87 * y<sub>i</sub> = ∫ f<sup>2</sup> from t<sub>1</sub> to t<sub>i</sub> 88 * z<sub>i</sub> = ∫ f'<sup>2</sup> from t<sub>1</sub> to t<sub>i</sub> 89 * 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> 90 * end for 91 * 92 * |-------------------------- 93 * \ | ∑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> 94 * a = \ | ------------------------ 95 * \| ∑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> 96 * 97 * 98 * |-------------------------- 99 * \ | ∑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 } 301
