View Javadoc

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.optimization.OptimizationException;
21  
22  /** This class guesses harmonic coefficients from a sample.
23  
24   * <p>The algorithm used to guess the coefficients is as follows:</p>
25  
26   * <p>We know f (t) at some sampling points t<sub>i</sub> and want to find a,
27   * &omega; and &phi; such that f (t) = a cos (&omega; t + &phi;).
28   * </p>
29   *
30   * <p>From the analytical expression, we can compute two primitives :
31   * <pre>
32   *     If2  (t) = &int; f<sup>2</sup>  = a<sup>2</sup> &times; [t + S (t)] / 2
33   *     If'2 (t) = &int; f'<sup>2</sup> = a<sup>2</sup> &omega;<sup>2</sup> &times; [t - S (t)] / 2
34   *     where S (t) = sin (2 (&omega; t + &phi;)) / (2 &omega;)
35   * </pre>
36   * </p>
37   *
38   * <p>We can remove S between these expressions :
39   * <pre>
40   *     If'2 (t) = a<sup>2</sup> &omega;<sup>2</sup> t - &omega;<sup>2</sup> If2 (t)
41   * </pre>
42   * </p>
43   *
44   * <p>The preceding expression shows that If'2 (t) is a linear
45   * combination of both t and If2 (t): If'2 (t) = A &times; t + B &times; If2 (t)
46   * </p>
47   *
48   * <p>From the primitive, we can deduce the same form for definite
49   * integrals between t<sub>1</sub> and t<sub>i</sub> for each t<sub>i</sub> :
50   * <pre>
51   *   If2 (t<sub>i</sub>) - If2 (t<sub>1</sub>) = A &times; (t<sub>i</sub> - t<sub>1</sub>) + B &times; (If2 (t<sub>i</sub>) - If2 (t<sub>1</sub>))
52   * </pre>
53   * </p>
54   *
55   * <p>We can find the coefficients A and B that best fit the sample
56   * to this linear expression by computing the definite integrals for
57   * each sample points.
58   * </p>
59   *
60   * <p>For a bilinear expression z (x<sub>i</sub>, y<sub>i</sub>) = A &times; x<sub>i</sub> + B &times; y<sub>i</sub>, the
61   * coefficients A and B that minimize a least square criterion
62   * &sum; (z<sub>i</sub> - z (x<sub>i</sub>, y<sub>i</sub>))<sup>2</sup> are given by these expressions:</p>
63   * <pre>
64   *
65   *         &sum;y<sub>i</sub>y<sub>i</sub> &sum;x<sub>i</sub>z<sub>i</sub> - &sum;x<sub>i</sub>y<sub>i</sub> &sum;y<sub>i</sub>z<sub>i</sub>
66   *     A = ------------------------
67   *         &sum;x<sub>i</sub>x<sub>i</sub> &sum;y<sub>i</sub>y<sub>i</sub> - &sum;x<sub>i</sub>y<sub>i</sub> &sum;x<sub>i</sub>y<sub>i</sub>
68   *
69   *         &sum;x<sub>i</sub>x<sub>i</sub> &sum;y<sub>i</sub>z<sub>i</sub> - &sum;x<sub>i</sub>y<sub>i</sub> &sum;x<sub>i</sub>z<sub>i</sub>
70   *     B = ------------------------
71   *         &sum;x<sub>i</sub>x<sub>i</sub> &sum;y<sub>i</sub>y<sub>i</sub> - &sum;x<sub>i</sub>y<sub>i</sub> &sum;x<sub>i</sub>y<sub>i</sub>
72   * </pre>
73   * </p>
74   *
75   *
76   * <p>In fact, we can assume both a and &omega; are positive and
77   * compute them directly, knowing that A = a<sup>2</sup> &omega;<sup>2</sup> and that
78   * B = - &omega;<sup>2</sup>. The complete algorithm is therefore:</p>
79   * <pre>
80   *
81   * for each t<sub>i</sub> from t<sub>1</sub> to t<sub>n-1</sub>, compute:
82   *   f  (t<sub>i</sub>)
83   *   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>)
84   *   x<sub>i</sub> = t<sub>i</sub> - t<sub>1</sub>
85   *   y<sub>i</sub> = &int; f<sup>2</sup> from t<sub>1</sub> to t<sub>i</sub>
86   *   z<sub>i</sub> = &int; f'<sup>2</sup> from t<sub>1</sub> to t<sub>i</sub>
87   *   update the sums &sum;x<sub>i</sub>x<sub>i</sub>, &sum;y<sub>i</sub>y<sub>i</sub>, &sum;x<sub>i</sub>y<sub>i</sub>, &sum;x<sub>i</sub>z<sub>i</sub> and &sum;y<sub>i</sub>z<sub>i</sub>
88   * end for
89   *
90   *            |--------------------------
91   *         \  | &sum;y<sub>i</sub>y<sub>i</sub> &sum;x<sub>i</sub>z<sub>i</sub> - &sum;x<sub>i</sub>y<sub>i</sub> &sum;y<sub>i</sub>z<sub>i</sub>
92   * a     =  \ | ------------------------
93   *           \| &sum;x<sub>i</sub>y<sub>i</sub> &sum;x<sub>i</sub>z<sub>i</sub> - &sum;x<sub>i</sub>x<sub>i</sub> &sum;y<sub>i</sub>z<sub>i</sub>
94   *
95   *
96   *            |--------------------------
97   *         \  | &sum;x<sub>i</sub>y<sub>i</sub> &sum;x<sub>i</sub>z<sub>i</sub> - &sum;x<sub>i</sub>x<sub>i</sub> &sum;y<sub>i</sub>z<sub>i</sub>
98   * &omega;     =  \ | ------------------------
99   *           \| &sum;x<sub>i</sub>x<sub>i</sub> &sum;y<sub>i</sub>y<sub>i</sub> - &sum;x<sub>i</sub>y<sub>i</sub> &sum;x<sub>i</sub>y<sub>i</sub>
100  *
101  * </pre>
102  * </p>
103 
104  * <p>Once we know &omega;, we can compute:
105  * <pre>
106  *    fc = &omega; f (t) cos (&omega; t) - f' (t) sin (&omega; t)
107  *    fs = &omega; f (t) sin (&omega; t) + f' (t) cos (&omega; t)
108  * </pre>
109  * </p>
110 
111  * <p>It appears that <code>fc = a &omega; cos (&phi;)</code> and
112  * <code>fs = -a &omega; sin (&phi;)</code>, so we can use these
113  * expressions to compute &phi;. 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 &omega;. */
134     private double omega;
135 
136     /** Guessed phase &phi;. */
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 &omega; 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 &phi; 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 &omega;.
287      * @return guessed pulsation &omega;
288      */
289     public double getGuessedPulsation() {
290         return omega;
291     }
292 
293     /** Get the guessed phase &phi;.
294      * @return guessed phase &phi;
295      */
296     public double getGuessedPhase() {
297         return phi;
298     }
299 
300 }