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.optimization.OptimizationException; 021 022 /** This class guesses harmonic coefficients from a sample. 023 024 * <p>The algorithm used to guess the coefficients is as follows:</p> 025 026 * <p>We know f (t) at some sampling points t<sub>i</sub> and want to find a, 027 * ω and φ such that f (t) = a cos (ω t + φ). 028 * </p> 029 * 030 * <p>From the analytical expression, we can compute two primitives : 031 * <pre> 032 * If2 (t) = ∫ f<sup>2</sup> = a<sup>2</sup> × [t + S (t)] / 2 033 * If'2 (t) = ∫ f'<sup>2</sup> = a<sup>2</sup> ω<sup>2</sup> × [t - S (t)] / 2 034 * where S (t) = sin (2 (ω t + φ)) / (2 ω) 035 * </pre> 036 * </p> 037 * 038 * <p>We can remove S between these expressions : 039 * <pre> 040 * If'2 (t) = a<sup>2</sup> ω<sup>2</sup> t - ω<sup>2</sup> If2 (t) 041 * </pre> 042 * </p> 043 * 044 * <p>The preceding expression shows that If'2 (t) is a linear 045 * combination of both t and If2 (t): If'2 (t) = A × t + B × If2 (t) 046 * </p> 047 * 048 * <p>From the primitive, we can deduce the same form for definite 049 * integrals between t<sub>1</sub> and t<sub>i</sub> for each t<sub>i</sub> : 050 * <pre> 051 * 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>)) 052 * </pre> 053 * </p> 054 * 055 * <p>We can find the coefficients A and B that best fit the sample 056 * to this linear expression by computing the definite integrals for 057 * each sample points. 058 * </p> 059 * 060 * <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 061 * coefficients A and B that minimize a least square criterion 062 * ∑ (z<sub>i</sub> - z (x<sub>i</sub>, y<sub>i</sub>))<sup>2</sup> are given by these expressions:</p> 063 * <pre> 064 * 065 * ∑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> 066 * A = ------------------------ 067 * ∑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> 068 * 069 * ∑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> 070 * B = ------------------------ 071 * ∑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> 072 * </pre> 073 * </p> 074 * 075 * 076 * <p>In fact, we can assume both a and ω are positive and 077 * compute them directly, knowing that A = a<sup>2</sup> ω<sup>2</sup> and that 078 * B = - ω<sup>2</sup>. The complete algorithm is therefore:</p> 079 * <pre> 080 * 081 * for each t<sub>i</sub> from t<sub>1</sub> to t<sub>n-1</sub>, compute: 082 * f (t<sub>i</sub>) 083 * 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>) 084 * x<sub>i</sub> = t<sub>i</sub> - t<sub>1</sub> 085 * y<sub>i</sub> = ∫ f<sup>2</sup> from t<sub>1</sub> to t<sub>i</sub> 086 * z<sub>i</sub> = ∫ f'<sup>2</sup> from t<sub>1</sub> to t<sub>i</sub> 087 * 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> 088 * end for 089 * 090 * |-------------------------- 091 * \ | ∑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> 092 * a = \ | ------------------------ 093 * \| ∑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> 094 * 095 * 096 * |-------------------------- 097 * \ | ∑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> 098 * ω = \ | ------------------------ 099 * \| ∑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> 100 * 101 * </pre> 102 * </p> 103 104 * <p>Once we know ω, we can compute: 105 * <pre> 106 * fc = ω f (t) cos (ω t) - f' (t) sin (ω t) 107 * fs = ω f (t) sin (ω t) + f' (t) cos (ω t) 108 * </pre> 109 * </p> 110 111 * <p>It appears that <code>fc = a ω cos (φ)</code> and 112 * <code>fs = -a ω sin (φ)</code>, so we can use these 113 * expressions to compute φ. The best estimate over the sample is 114 * given by averaging these expressions. 115 * </p> 116 117 * <p>Since integrals and means are involved in the preceding 118 * estimations, these operations run in O(n) time, where n is the 119 * number of measurements.</p> 120 121 * @version $Revision: 786479 $ $Date: 2009-06-19 08:36:16 -0400 (Fri, 19 Jun 2009) $ 122 * @since 2.0 123 124 */ 125 public class HarmonicCoefficientsGuesser { 126 127 /** Sampled observations. */ 128 private final WeightedObservedPoint[] observations; 129 130 /** Guessed amplitude. */ 131 private double a; 132 133 /** Guessed pulsation ω. */ 134 private double omega; 135 136 /** Guessed phase φ. */ 137 private double phi; 138 139 /** Simple constructor. 140 * @param observations sampled observations 141 */ 142 public HarmonicCoefficientsGuesser(WeightedObservedPoint[] observations) { 143 this.observations = observations.clone(); 144 a = Double.NaN; 145 omega = Double.NaN; 146 } 147 148 /** Estimate a first guess of the coefficients. 149 * @exception OptimizationException if the sample is too short or if 150 * the first guess cannot be computed (when the elements under the 151 * square roots are negative). 152 * */ 153 public void guess() throws OptimizationException { 154 sortObservations(); 155 guessAOmega(); 156 guessPhi(); 157 } 158 159 /** Sort the observations with respect to the abscissa. 160 */ 161 private void sortObservations() { 162 163 // Since the samples are almost always already sorted, this 164 // method is implemented as an insertion sort that reorders the 165 // elements in place. Insertion sort is very efficient in this case. 166 WeightedObservedPoint curr = observations[0]; 167 for (int j = 1; j < observations.length; ++j) { 168 WeightedObservedPoint prec = curr; 169 curr = observations[j]; 170 if (curr.getX() < prec.getX()) { 171 // the current element should be inserted closer to the beginning 172 int i = j - 1; 173 WeightedObservedPoint mI = observations[i]; 174 while ((i >= 0) && (curr.getX() < mI.getX())) { 175 observations[i + 1] = mI; 176 if (i-- != 0) { 177 mI = observations[i]; 178 } else { 179 mI = null; 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("unable to first guess the harmonic coefficients"); 242 } 243 a = Math.sqrt(c1 / c2); 244 omega = Math.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 = Math.cos(omegaX); 269 double sine = Math.sin(omegaX); 270 fcMean += omega * currentY * cosine - currentYPrime * sine; 271 fsMean += omega * currentY * sine + currentYPrime * cosine; 272 273 } 274 275 phi = Math.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 }