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 }