1dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond/* 2dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond * Licensed to the Apache Software Foundation (ASF) under one or more 3dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond * contributor license agreements. See the NOTICE file distributed with 4dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond * this work for additional information regarding copyright ownership. 5dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond * The ASF licenses this file to You under the Apache License, Version 2.0 6dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond * (the "License"); you may not use this file except in compliance with 7dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond * the License. You may obtain a copy of the License at 8dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond * 9dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond * http://www.apache.org/licenses/LICENSE-2.0 10dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond * 11dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond * Unless required by applicable law or agreed to in writing, software 12dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond * distributed under the License is distributed on an "AS IS" BASIS, 13dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond * WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied. 14dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond * See the License for the specific language governing permissions and 15dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond * limitations under the License. 16dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond */ 17dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond 18dee0849a9704d532af0b550146cbafbaa6ee1d19Raymondpackage org.apache.commons.math.optimization.fitting; 19dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond 20dee0849a9704d532af0b550146cbafbaa6ee1d19Raymondimport org.apache.commons.math.exception.util.LocalizedFormats; 21dee0849a9704d532af0b550146cbafbaa6ee1d19Raymondimport org.apache.commons.math.optimization.OptimizationException; 22dee0849a9704d532af0b550146cbafbaa6ee1d19Raymondimport org.apache.commons.math.util.FastMath; 23dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond 24dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond/** This class guesses harmonic coefficients from a sample. 25dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond 26dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond * <p>The algorithm used to guess the coefficients is as follows:</p> 27dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond 28dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond * <p>We know f (t) at some sampling points t<sub>i</sub> and want to find a, 29dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond * ω and φ such that f (t) = a cos (ω t + φ). 30dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond * </p> 31dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond * 32dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond * <p>From the analytical expression, we can compute two primitives : 33dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond * <pre> 34dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond * If2 (t) = ∫ f<sup>2</sup> = a<sup>2</sup> × [t + S (t)] / 2 35dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond * If'2 (t) = ∫ f'<sup>2</sup> = a<sup>2</sup> ω<sup>2</sup> × [t - S (t)] / 2 36dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond * where S (t) = sin (2 (ω t + φ)) / (2 ω) 37dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond * </pre> 38dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond * </p> 39dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond * 40dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond * <p>We can remove S between these expressions : 41dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond * <pre> 42dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond * If'2 (t) = a<sup>2</sup> ω<sup>2</sup> t - ω<sup>2</sup> If2 (t) 43dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond * </pre> 44dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond * </p> 45dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond * 46dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond * <p>The preceding expression shows that If'2 (t) is a linear 47dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond * combination of both t and If2 (t): If'2 (t) = A × t + B × If2 (t) 48dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond * </p> 49dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond * 50dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond * <p>From the primitive, we can deduce the same form for definite 51dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond * integrals between t<sub>1</sub> and t<sub>i</sub> for each t<sub>i</sub> : 52dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond * <pre> 53dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond * 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>)) 54dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond * </pre> 55dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond * </p> 56dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond * 57dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond * <p>We can find the coefficients A and B that best fit the sample 58dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond * to this linear expression by computing the definite integrals for 59dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond * each sample points. 60dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond * </p> 61dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond * 62dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond * <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 63dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond * coefficients A and B that minimize a least square criterion 64dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond * ∑ (z<sub>i</sub> - z (x<sub>i</sub>, y<sub>i</sub>))<sup>2</sup> are given by these expressions:</p> 65dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond * <pre> 66dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond * 67dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond * ∑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> 68dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond * A = ------------------------ 69dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond * ∑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> 70dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond * 71dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond * ∑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> 72dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond * B = ------------------------ 73dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond * ∑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> 74dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond * </pre> 75dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond * </p> 76dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond * 77dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond * 78dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond * <p>In fact, we can assume both a and ω are positive and 79dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond * compute them directly, knowing that A = a<sup>2</sup> ω<sup>2</sup> and that 80dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond * B = - ω<sup>2</sup>. The complete algorithm is therefore:</p> 81dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond * <pre> 82dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond * 83dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond * for each t<sub>i</sub> from t<sub>1</sub> to t<sub>n-1</sub>, compute: 84dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond * f (t<sub>i</sub>) 85dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond * 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>) 86dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond * x<sub>i</sub> = t<sub>i</sub> - t<sub>1</sub> 87dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond * y<sub>i</sub> = ∫ f<sup>2</sup> from t<sub>1</sub> to t<sub>i</sub> 88dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond * z<sub>i</sub> = ∫ f'<sup>2</sup> from t<sub>1</sub> to t<sub>i</sub> 89dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond * 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> 90dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond * end for 91dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond * 92dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond * |-------------------------- 93dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond * \ | ∑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> 94dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond * a = \ | ------------------------ 95dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond * \| ∑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> 96dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond * 97dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond * 98dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond * |-------------------------- 99dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond * \ | ∑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> 100dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond * ω = \ | ------------------------ 101dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond * \| ∑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> 102dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond * 103dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond * </pre> 104dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond * </p> 105dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond 106dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond * <p>Once we know ω, we can compute: 107dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond * <pre> 108dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond * fc = ω f (t) cos (ω t) - f' (t) sin (ω t) 109dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond * fs = ω f (t) sin (ω t) + f' (t) cos (ω t) 110dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond * </pre> 111dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond * </p> 112dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond 113dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond * <p>It appears that <code>fc = a ω cos (φ)</code> and 114dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond * <code>fs = -a ω sin (φ)</code>, so we can use these 115dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond * expressions to compute φ. The best estimate over the sample is 116dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond * given by averaging these expressions. 117dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond * </p> 118dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond 119dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond * <p>Since integrals and means are involved in the preceding 120dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond * estimations, these operations run in O(n) time, where n is the 121dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond * number of measurements.</p> 122dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond 123dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond * @version $Revision: 1056034 $ $Date: 2011-01-06 20:41:43 +0100 (jeu. 06 janv. 2011) $ 124dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond * @since 2.0 125dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond 126dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond */ 127dee0849a9704d532af0b550146cbafbaa6ee1d19Raymondpublic class HarmonicCoefficientsGuesser { 128dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond 129dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond /** Sampled observations. */ 130dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond private final WeightedObservedPoint[] observations; 131dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond 132dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond /** Guessed amplitude. */ 133dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond private double a; 134dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond 135dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond /** Guessed pulsation ω. */ 136dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond private double omega; 137dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond 138dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond /** Guessed phase φ. */ 139dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond private double phi; 140dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond 141dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond /** Simple constructor. 142dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond * @param observations sampled observations 143dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond */ 144dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond public HarmonicCoefficientsGuesser(WeightedObservedPoint[] observations) { 145dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond this.observations = observations.clone(); 146dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond a = Double.NaN; 147dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond omega = Double.NaN; 148dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond } 149dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond 150dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond /** Estimate a first guess of the coefficients. 151dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond * @exception OptimizationException if the sample is too short or if 152dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond * the first guess cannot be computed (when the elements under the 153dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond * square roots are negative). 154dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond * */ 155dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond public void guess() throws OptimizationException { 156dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond sortObservations(); 157dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond guessAOmega(); 158dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond guessPhi(); 159dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond } 160dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond 161dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond /** Sort the observations with respect to the abscissa. 162dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond */ 163dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond private void sortObservations() { 164dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond 165dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond // Since the samples are almost always already sorted, this 166dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond // method is implemented as an insertion sort that reorders the 167dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond // elements in place. Insertion sort is very efficient in this case. 168dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond WeightedObservedPoint curr = observations[0]; 169dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond for (int j = 1; j < observations.length; ++j) { 170dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond WeightedObservedPoint prec = curr; 171dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond curr = observations[j]; 172dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond if (curr.getX() < prec.getX()) { 173dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond // the current element should be inserted closer to the beginning 174dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond int i = j - 1; 175dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond WeightedObservedPoint mI = observations[i]; 176dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond while ((i >= 0) && (curr.getX() < mI.getX())) { 177dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond observations[i + 1] = mI; 178dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond if (i-- != 0) { 179dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond mI = observations[i]; 180dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond } 181dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond } 182dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond observations[i + 1] = curr; 183dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond curr = observations[j]; 184dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond } 185dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond } 186dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond 187dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond } 188dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond 189dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond /** Estimate a first guess of the a and ω coefficients. 190dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond * @exception OptimizationException if the sample is too short or if 191dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond * the first guess cannot be computed (when the elements under the 192dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond * square roots are negative). 193dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond */ 194dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond private void guessAOmega() throws OptimizationException { 195dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond 196dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond // initialize the sums for the linear model between the two integrals 197dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond double sx2 = 0.0; 198dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond double sy2 = 0.0; 199dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond double sxy = 0.0; 200dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond double sxz = 0.0; 201dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond double syz = 0.0; 202dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond 203dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond double currentX = observations[0].getX(); 204dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond double currentY = observations[0].getY(); 205dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond double f2Integral = 0; 206dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond double fPrime2Integral = 0; 207dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond final double startX = currentX; 208dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond for (int i = 1; i < observations.length; ++i) { 209dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond 210dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond // one step forward 211dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond final double previousX = currentX; 212dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond final double previousY = currentY; 213dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond currentX = observations[i].getX(); 214dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond currentY = observations[i].getY(); 215dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond 216dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond // update the integrals of f<sup>2</sup> and f'<sup>2</sup> 217dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond // considering a linear model for f (and therefore constant f') 218dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond final double dx = currentX - previousX; 219dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond final double dy = currentY - previousY; 220dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond final double f2StepIntegral = 221dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond dx * (previousY * previousY + previousY * currentY + currentY * currentY) / 3; 222dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond final double fPrime2StepIntegral = dy * dy / dx; 223dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond 224dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond final double x = currentX - startX; 225dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond f2Integral += f2StepIntegral; 226dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond fPrime2Integral += fPrime2StepIntegral; 227dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond 228dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond sx2 += x * x; 229dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond sy2 += f2Integral * f2Integral; 230dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond sxy += x * f2Integral; 231dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond sxz += x * fPrime2Integral; 232dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond syz += f2Integral * fPrime2Integral; 233dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond 234dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond } 235dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond 236dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond // compute the amplitude and pulsation coefficients 237dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond double c1 = sy2 * sxz - sxy * syz; 238dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond double c2 = sxy * sxz - sx2 * syz; 239dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond double c3 = sx2 * sy2 - sxy * sxy; 240dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond if ((c1 / c2 < 0.0) || (c2 / c3 < 0.0)) { 241dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond throw new OptimizationException(LocalizedFormats.UNABLE_TO_FIRST_GUESS_HARMONIC_COEFFICIENTS); 242dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond } 243dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond a = FastMath.sqrt(c1 / c2); 244dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond omega = FastMath.sqrt(c2 / c3); 245dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond 246dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond } 247dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond 248dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond /** Estimate a first guess of the φ coefficient. 249dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond */ 250dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond private void guessPhi() { 251dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond 252dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond // initialize the means 253dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond double fcMean = 0.0; 254dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond double fsMean = 0.0; 255dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond 256dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond double currentX = observations[0].getX(); 257dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond double currentY = observations[0].getY(); 258dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond for (int i = 1; i < observations.length; ++i) { 259dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond 260dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond // one step forward 261dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond final double previousX = currentX; 262dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond final double previousY = currentY; 263dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond currentX = observations[i].getX(); 264dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond currentY = observations[i].getY(); 265dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond final double currentYPrime = (currentY - previousY) / (currentX - previousX); 266dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond 267dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond double omegaX = omega * currentX; 268dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond double cosine = FastMath.cos(omegaX); 269dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond double sine = FastMath.sin(omegaX); 270dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond fcMean += omega * currentY * cosine - currentYPrime * sine; 271dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond fsMean += omega * currentY * sine + currentYPrime * cosine; 272dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond 273dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond } 274dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond 275dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond phi = FastMath.atan2(-fsMean, fcMean); 276dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond 277dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond } 278dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond 279dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond /** Get the guessed amplitude a. 280dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond * @return guessed amplitude a; 281dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond */ 282dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond public double getGuessedAmplitude() { 283dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond return a; 284dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond } 285dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond 286dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond /** Get the guessed pulsation ω. 287dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond * @return guessed pulsation ω 288dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond */ 289dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond public double getGuessedPulsation() { 290dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond return omega; 291dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond } 292dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond 293dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond /** Get the guessed phase φ. 294dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond * @return guessed phase φ 295dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond */ 296dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond public double getGuessedPhase() { 297dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond return phi; 298dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond } 299dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond 300dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond} 301