Mercurial > hg > de.mpg.mpiwg.itgroup.digilib.plugin

comparison libs/commons-math-2.1/docs/apidocs/src-html/org/apache/commons/math/optimization/fitting/HarmonicCoefficientsGuesser.html @ 13:cbf34dd4d7e6

Find changesets by keywords (author, files, the commit message), revision number or hash, or revset expression.
commons-math-2.1 added
author dwinter
date Tue, 04 Jan 2011 10:02:07 +0100
parents
children
comparison
equal deleted inserted replaced
12:970d26a94fb7 13:cbf34dd4d7e6
1 <HTML>
2 <BODY BGCOLOR="white">
3 <PRE>
4 <FONT color="green">001</FONT> /*<a name="line.1"></a>
5 <FONT color="green">002</FONT> * Licensed to the Apache Software Foundation (ASF) under one or more<a name="line.2"></a>
6 <FONT color="green">003</FONT> * contributor license agreements. See the NOTICE file distributed with<a name="line.3"></a>
7 <FONT color="green">004</FONT> * this work for additional information regarding copyright ownership.<a name="line.4"></a>
8 <FONT color="green">005</FONT> * The ASF licenses this file to You under the Apache License, Version 2.0<a name="line.5"></a>
9 <FONT color="green">006</FONT> * (the "License"); you may not use this file except in compliance with<a name="line.6"></a>
10 <FONT color="green">007</FONT> * the License. You may obtain a copy of the License at<a name="line.7"></a>
11 <FONT color="green">008</FONT> *<a name="line.8"></a>
12 <FONT color="green">009</FONT> * http://www.apache.org/licenses/LICENSE-2.0<a name="line.9"></a>
13 <FONT color="green">010</FONT> *<a name="line.10"></a>
14 <FONT color="green">011</FONT> * Unless required by applicable law or agreed to in writing, software<a name="line.11"></a>
15 <FONT color="green">012</FONT> * distributed under the License is distributed on an "AS IS" BASIS,<a name="line.12"></a>
16 <FONT color="green">013</FONT> * WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.<a name="line.13"></a>
17 <FONT color="green">014</FONT> * See the License for the specific language governing permissions and<a name="line.14"></a>
18 <FONT color="green">015</FONT> * limitations under the License.<a name="line.15"></a>
19 <FONT color="green">016</FONT> */<a name="line.16"></a>
20 <FONT color="green">017</FONT> <a name="line.17"></a>
21 <FONT color="green">018</FONT> package org.apache.commons.math.optimization.fitting;<a name="line.18"></a>
22 <FONT color="green">019</FONT> <a name="line.19"></a>
23 <FONT color="green">020</FONT> import org.apache.commons.math.optimization.OptimizationException;<a name="line.20"></a>
24 <FONT color="green">021</FONT> <a name="line.21"></a>
25 <FONT color="green">022</FONT> /** This class guesses harmonic coefficients from a sample.<a name="line.22"></a>
26 <FONT color="green">023</FONT> <a name="line.23"></a>
27 <FONT color="green">024</FONT> * &lt;p&gt;The algorithm used to guess the coefficients is as follows:&lt;/p&gt;<a name="line.24"></a>
28 <FONT color="green">025</FONT> <a name="line.25"></a>
29 <FONT color="green">026</FONT> * &lt;p&gt;We know f (t) at some sampling points t&lt;sub&gt;i&lt;/sub&gt; and want to find a,<a name="line.26"></a>
30 <FONT color="green">027</FONT> * &amp;omega; and &amp;phi; such that f (t) = a cos (&amp;omega; t + &amp;phi;).<a name="line.27"></a>
31 <FONT color="green">028</FONT> * &lt;/p&gt;<a name="line.28"></a>
32 <FONT color="green">029</FONT> *<a name="line.29"></a>
33 <FONT color="green">030</FONT> * &lt;p&gt;From the analytical expression, we can compute two primitives :<a name="line.30"></a>
34 <FONT color="green">031</FONT> * &lt;pre&gt;<a name="line.31"></a>
35 <FONT color="green">032</FONT> * If2 (t) = &amp;int; f&lt;sup&gt;2&lt;/sup&gt; = a&lt;sup&gt;2&lt;/sup&gt; &amp;times; [t + S (t)] / 2<a name="line.32"></a>
36 <FONT color="green">033</FONT> * If'2 (t) = &amp;int; f'&lt;sup&gt;2&lt;/sup&gt; = a&lt;sup&gt;2&lt;/sup&gt; &amp;omega;&lt;sup&gt;2&lt;/sup&gt; &amp;times; [t - S (t)] / 2<a name="line.33"></a>
37 <FONT color="green">034</FONT> * where S (t) = sin (2 (&amp;omega; t + &amp;phi;)) / (2 &amp;omega;)<a name="line.34"></a>
38 <FONT color="green">035</FONT> * &lt;/pre&gt;<a name="line.35"></a>
39 <FONT color="green">036</FONT> * &lt;/p&gt;<a name="line.36"></a>
40 <FONT color="green">037</FONT> *<a name="line.37"></a>
41 <FONT color="green">038</FONT> * &lt;p&gt;We can remove S between these expressions :<a name="line.38"></a>
42 <FONT color="green">039</FONT> * &lt;pre&gt;<a name="line.39"></a>
43 <FONT color="green">040</FONT> * If'2 (t) = a&lt;sup&gt;2&lt;/sup&gt; &amp;omega;&lt;sup&gt;2&lt;/sup&gt; t - &amp;omega;&lt;sup&gt;2&lt;/sup&gt; If2 (t)<a name="line.40"></a>
44 <FONT color="green">041</FONT> * &lt;/pre&gt;<a name="line.41"></a>
45 <FONT color="green">042</FONT> * &lt;/p&gt;<a name="line.42"></a>
46 <FONT color="green">043</FONT> *<a name="line.43"></a>
47 <FONT color="green">044</FONT> * &lt;p&gt;The preceding expression shows that If'2 (t) is a linear<a name="line.44"></a>
48 <FONT color="green">045</FONT> * combination of both t and If2 (t): If'2 (t) = A &amp;times; t + B &amp;times; If2 (t)<a name="line.45"></a>
49 <FONT color="green">046</FONT> * &lt;/p&gt;<a name="line.46"></a>
50 <FONT color="green">047</FONT> *<a name="line.47"></a>
51 <FONT color="green">048</FONT> * &lt;p&gt;From the primitive, we can deduce the same form for definite<a name="line.48"></a>
52 <FONT color="green">049</FONT> * integrals between t&lt;sub&gt;1&lt;/sub&gt; and t&lt;sub&gt;i&lt;/sub&gt; for each t&lt;sub&gt;i&lt;/sub&gt; :<a name="line.49"></a>
53 <FONT color="green">050</FONT> * &lt;pre&gt;<a name="line.50"></a>
54 <FONT color="green">051</FONT> * If2 (t&lt;sub&gt;i&lt;/sub&gt;) - If2 (t&lt;sub&gt;1&lt;/sub&gt;) = A &amp;times; (t&lt;sub&gt;i&lt;/sub&gt; - t&lt;sub&gt;1&lt;/sub&gt;) + B &amp;times; (If2 (t&lt;sub&gt;i&lt;/sub&gt;) - If2 (t&lt;sub&gt;1&lt;/sub&gt;))<a name="line.51"></a>
55 <FONT color="green">052</FONT> * &lt;/pre&gt;<a name="line.52"></a>
56 <FONT color="green">053</FONT> * &lt;/p&gt;<a name="line.53"></a>
57 <FONT color="green">054</FONT> *<a name="line.54"></a>
58 <FONT color="green">055</FONT> * &lt;p&gt;We can find the coefficients A and B that best fit the sample<a name="line.55"></a>
59 <FONT color="green">056</FONT> * to this linear expression by computing the definite integrals for<a name="line.56"></a>
60 <FONT color="green">057</FONT> * each sample points.<a name="line.57"></a>
61 <FONT color="green">058</FONT> * &lt;/p&gt;<a name="line.58"></a>
62 <FONT color="green">059</FONT> *<a name="line.59"></a>
63 <FONT color="green">060</FONT> * &lt;p&gt;For a bilinear expression z (x&lt;sub&gt;i&lt;/sub&gt;, y&lt;sub&gt;i&lt;/sub&gt;) = A &amp;times; x&lt;sub&gt;i&lt;/sub&gt; + B &amp;times; y&lt;sub&gt;i&lt;/sub&gt;, the<a name="line.60"></a>
64 <FONT color="green">061</FONT> * coefficients A and B that minimize a least square criterion<a name="line.61"></a>
65 <FONT color="green">062</FONT> * &amp;sum; (z&lt;sub&gt;i&lt;/sub&gt; - z (x&lt;sub&gt;i&lt;/sub&gt;, y&lt;sub&gt;i&lt;/sub&gt;))&lt;sup&gt;2&lt;/sup&gt; are given by these expressions:&lt;/p&gt;<a name="line.62"></a>
66 <FONT color="green">063</FONT> * &lt;pre&gt;<a name="line.63"></a>
67 <FONT color="green">064</FONT> *<a name="line.64"></a>
68 <FONT color="green">065</FONT> * &amp;sum;y&lt;sub&gt;i&lt;/sub&gt;y&lt;sub&gt;i&lt;/sub&gt; &amp;sum;x&lt;sub&gt;i&lt;/sub&gt;z&lt;sub&gt;i&lt;/sub&gt; - &amp;sum;x&lt;sub&gt;i&lt;/sub&gt;y&lt;sub&gt;i&lt;/sub&gt; &amp;sum;y&lt;sub&gt;i&lt;/sub&gt;z&lt;sub&gt;i&lt;/sub&gt;<a name="line.65"></a>
69 <FONT color="green">066</FONT> * A = ------------------------<a name="line.66"></a>
70 <FONT color="green">067</FONT> * &amp;sum;x&lt;sub&gt;i&lt;/sub&gt;x&lt;sub&gt;i&lt;/sub&gt; &amp;sum;y&lt;sub&gt;i&lt;/sub&gt;y&lt;sub&gt;i&lt;/sub&gt; - &amp;sum;x&lt;sub&gt;i&lt;/sub&gt;y&lt;sub&gt;i&lt;/sub&gt; &amp;sum;x&lt;sub&gt;i&lt;/sub&gt;y&lt;sub&gt;i&lt;/sub&gt;<a name="line.67"></a>
71 <FONT color="green">068</FONT> *<a name="line.68"></a>
72 <FONT color="green">069</FONT> * &amp;sum;x&lt;sub&gt;i&lt;/sub&gt;x&lt;sub&gt;i&lt;/sub&gt; &amp;sum;y&lt;sub&gt;i&lt;/sub&gt;z&lt;sub&gt;i&lt;/sub&gt; - &amp;sum;x&lt;sub&gt;i&lt;/sub&gt;y&lt;sub&gt;i&lt;/sub&gt; &amp;sum;x&lt;sub&gt;i&lt;/sub&gt;z&lt;sub&gt;i&lt;/sub&gt;<a name="line.69"></a>
73 <FONT color="green">070</FONT> * B = ------------------------<a name="line.70"></a>
74 <FONT color="green">071</FONT> * &amp;sum;x&lt;sub&gt;i&lt;/sub&gt;x&lt;sub&gt;i&lt;/sub&gt; &amp;sum;y&lt;sub&gt;i&lt;/sub&gt;y&lt;sub&gt;i&lt;/sub&gt; - &amp;sum;x&lt;sub&gt;i&lt;/sub&gt;y&lt;sub&gt;i&lt;/sub&gt; &amp;sum;x&lt;sub&gt;i&lt;/sub&gt;y&lt;sub&gt;i&lt;/sub&gt;<a name="line.71"></a>
75 <FONT color="green">072</FONT> * &lt;/pre&gt;<a name="line.72"></a>
76 <FONT color="green">073</FONT> * &lt;/p&gt;<a name="line.73"></a>
77 <FONT color="green">074</FONT> *<a name="line.74"></a>
78 <FONT color="green">075</FONT> *<a name="line.75"></a>
79 <FONT color="green">076</FONT> * &lt;p&gt;In fact, we can assume both a and &amp;omega; are positive and<a name="line.76"></a>
80 <FONT color="green">077</FONT> * compute them directly, knowing that A = a&lt;sup&gt;2&lt;/sup&gt; &amp;omega;&lt;sup&gt;2&lt;/sup&gt; and that<a name="line.77"></a>
81 <FONT color="green">078</FONT> * B = - &amp;omega;&lt;sup&gt;2&lt;/sup&gt;. The complete algorithm is therefore:&lt;/p&gt;<a name="line.78"></a>
82 <FONT color="green">079</FONT> * &lt;pre&gt;<a name="line.79"></a>
83 <FONT color="green">080</FONT> *<a name="line.80"></a>
84 <FONT color="green">081</FONT> * for each t&lt;sub&gt;i&lt;/sub&gt; from t&lt;sub&gt;1&lt;/sub&gt; to t&lt;sub&gt;n-1&lt;/sub&gt;, compute:<a name="line.81"></a>
85 <FONT color="green">082</FONT> * f (t&lt;sub&gt;i&lt;/sub&gt;)<a name="line.82"></a>
86 <FONT color="green">083</FONT> * f' (t&lt;sub&gt;i&lt;/sub&gt;) = (f (t&lt;sub&gt;i+1&lt;/sub&gt;) - f(t&lt;sub&gt;i-1&lt;/sub&gt;)) / (t&lt;sub&gt;i+1&lt;/sub&gt; - t&lt;sub&gt;i-1&lt;/sub&gt;)<a name="line.83"></a>
87 <FONT color="green">084</FONT> * x&lt;sub&gt;i&lt;/sub&gt; = t&lt;sub&gt;i&lt;/sub&gt; - t&lt;sub&gt;1&lt;/sub&gt;<a name="line.84"></a>
88 <FONT color="green">085</FONT> * y&lt;sub&gt;i&lt;/sub&gt; = &amp;int; f&lt;sup&gt;2&lt;/sup&gt; from t&lt;sub&gt;1&lt;/sub&gt; to t&lt;sub&gt;i&lt;/sub&gt;<a name="line.85"></a>
89 <FONT color="green">086</FONT> * z&lt;sub&gt;i&lt;/sub&gt; = &amp;int; f'&lt;sup&gt;2&lt;/sup&gt; from t&lt;sub&gt;1&lt;/sub&gt; to t&lt;sub&gt;i&lt;/sub&gt;<a name="line.86"></a>
90 <FONT color="green">087</FONT> * update the sums &amp;sum;x&lt;sub&gt;i&lt;/sub&gt;x&lt;sub&gt;i&lt;/sub&gt;, &amp;sum;y&lt;sub&gt;i&lt;/sub&gt;y&lt;sub&gt;i&lt;/sub&gt;, &amp;sum;x&lt;sub&gt;i&lt;/sub&gt;y&lt;sub&gt;i&lt;/sub&gt;, &amp;sum;x&lt;sub&gt;i&lt;/sub&gt;z&lt;sub&gt;i&lt;/sub&gt; and &amp;sum;y&lt;sub&gt;i&lt;/sub&gt;z&lt;sub&gt;i&lt;/sub&gt;<a name="line.87"></a>
91 <FONT color="green">088</FONT> * end for<a name="line.88"></a>
92 <FONT color="green">089</FONT> *<a name="line.89"></a>
93 <FONT color="green">090</FONT> * |--------------------------<a name="line.90"></a>
94 <FONT color="green">091</FONT> * \ | &amp;sum;y&lt;sub&gt;i&lt;/sub&gt;y&lt;sub&gt;i&lt;/sub&gt; &amp;sum;x&lt;sub&gt;i&lt;/sub&gt;z&lt;sub&gt;i&lt;/sub&gt; - &amp;sum;x&lt;sub&gt;i&lt;/sub&gt;y&lt;sub&gt;i&lt;/sub&gt; &amp;sum;y&lt;sub&gt;i&lt;/sub&gt;z&lt;sub&gt;i&lt;/sub&gt;<a name="line.91"></a>
95 <FONT color="green">092</FONT> * a = \ | ------------------------<a name="line.92"></a>
96 <FONT color="green">093</FONT> * \| &amp;sum;x&lt;sub&gt;i&lt;/sub&gt;y&lt;sub&gt;i&lt;/sub&gt; &amp;sum;x&lt;sub&gt;i&lt;/sub&gt;z&lt;sub&gt;i&lt;/sub&gt; - &amp;sum;x&lt;sub&gt;i&lt;/sub&gt;x&lt;sub&gt;i&lt;/sub&gt; &amp;sum;y&lt;sub&gt;i&lt;/sub&gt;z&lt;sub&gt;i&lt;/sub&gt;<a name="line.93"></a>
97 <FONT color="green">094</FONT> *<a name="line.94"></a>
98 <FONT color="green">095</FONT> *<a name="line.95"></a>
99 <FONT color="green">096</FONT> * |--------------------------<a name="line.96"></a>
100 <FONT color="green">097</FONT> * \ | &amp;sum;x&lt;sub&gt;i&lt;/sub&gt;y&lt;sub&gt;i&lt;/sub&gt; &amp;sum;x&lt;sub&gt;i&lt;/sub&gt;z&lt;sub&gt;i&lt;/sub&gt; - &amp;sum;x&lt;sub&gt;i&lt;/sub&gt;x&lt;sub&gt;i&lt;/sub&gt; &amp;sum;y&lt;sub&gt;i&lt;/sub&gt;z&lt;sub&gt;i&lt;/sub&gt;<a name="line.97"></a>
101 <FONT color="green">098</FONT> * &amp;omega; = \ | ------------------------<a name="line.98"></a>
102 <FONT color="green">099</FONT> * \| &amp;sum;x&lt;sub&gt;i&lt;/sub&gt;x&lt;sub&gt;i&lt;/sub&gt; &amp;sum;y&lt;sub&gt;i&lt;/sub&gt;y&lt;sub&gt;i&lt;/sub&gt; - &amp;sum;x&lt;sub&gt;i&lt;/sub&gt;y&lt;sub&gt;i&lt;/sub&gt; &amp;sum;x&lt;sub&gt;i&lt;/sub&gt;y&lt;sub&gt;i&lt;/sub&gt;<a name="line.99"></a>
103 <FONT color="green">100</FONT> *<a name="line.100"></a>
104 <FONT color="green">101</FONT> * &lt;/pre&gt;<a name="line.101"></a>
105 <FONT color="green">102</FONT> * &lt;/p&gt;<a name="line.102"></a>
106 <FONT color="green">103</FONT> <a name="line.103"></a>
107 <FONT color="green">104</FONT> * &lt;p&gt;Once we know &amp;omega;, we can compute:<a name="line.104"></a>
108 <FONT color="green">105</FONT> * &lt;pre&gt;<a name="line.105"></a>
109 <FONT color="green">106</FONT> * fc = &amp;omega; f (t) cos (&amp;omega; t) - f' (t) sin (&amp;omega; t)<a name="line.106"></a>
110 <FONT color="green">107</FONT> * fs = &amp;omega; f (t) sin (&amp;omega; t) + f' (t) cos (&amp;omega; t)<a name="line.107"></a>
111 <FONT color="green">108</FONT> * &lt;/pre&gt;<a name="line.108"></a>
112 <FONT color="green">109</FONT> * &lt;/p&gt;<a name="line.109"></a>
113 <FONT color="green">110</FONT> <a name="line.110"></a>
114 <FONT color="green">111</FONT> * &lt;p&gt;It appears that &lt;code&gt;fc = a &amp;omega; cos (&amp;phi;)&lt;/code&gt; and<a name="line.111"></a>
115 <FONT color="green">112</FONT> * &lt;code&gt;fs = -a &amp;omega; sin (&amp;phi;)&lt;/code&gt;, so we can use these<a name="line.112"></a>
116 <FONT color="green">113</FONT> * expressions to compute &amp;phi;. The best estimate over the sample is<a name="line.113"></a>
117 <FONT color="green">114</FONT> * given by averaging these expressions.<a name="line.114"></a>
118 <FONT color="green">115</FONT> * &lt;/p&gt;<a name="line.115"></a>
119 <FONT color="green">116</FONT> <a name="line.116"></a>
120 <FONT color="green">117</FONT> * &lt;p&gt;Since integrals and means are involved in the preceding<a name="line.117"></a>
121 <FONT color="green">118</FONT> * estimations, these operations run in O(n) time, where n is the<a name="line.118"></a>
122 <FONT color="green">119</FONT> * number of measurements.&lt;/p&gt;<a name="line.119"></a>
123 <FONT color="green">120</FONT> <a name="line.120"></a>
124 <FONT color="green">121</FONT> * @version $Revision: 786479 $ $Date: 2009-06-19 08:36:16 -0400 (Fri, 19 Jun 2009) $<a name="line.121"></a>
125 <FONT color="green">122</FONT> * @since 2.0<a name="line.122"></a>
126 <FONT color="green">123</FONT> <a name="line.123"></a>
127 <FONT color="green">124</FONT> */<a name="line.124"></a>
128 <FONT color="green">125</FONT> public class HarmonicCoefficientsGuesser {<a name="line.125"></a>
129 <FONT color="green">126</FONT> <a name="line.126"></a>
130 <FONT color="green">127</FONT> /** Sampled observations. */<a name="line.127"></a>
131 <FONT color="green">128</FONT> private final WeightedObservedPoint[] observations;<a name="line.128"></a>
132 <FONT color="green">129</FONT> <a name="line.129"></a>
133 <FONT color="green">130</FONT> /** Guessed amplitude. */<a name="line.130"></a>
134 <FONT color="green">131</FONT> private double a;<a name="line.131"></a>
135 <FONT color="green">132</FONT> <a name="line.132"></a>
136 <FONT color="green">133</FONT> /** Guessed pulsation &amp;omega;. */<a name="line.133"></a>
137 <FONT color="green">134</FONT> private double omega;<a name="line.134"></a>
138 <FONT color="green">135</FONT> <a name="line.135"></a>
139 <FONT color="green">136</FONT> /** Guessed phase &amp;phi;. */<a name="line.136"></a>
140 <FONT color="green">137</FONT> private double phi;<a name="line.137"></a>
141 <FONT color="green">138</FONT> <a name="line.138"></a>
142 <FONT color="green">139</FONT> /** Simple constructor.<a name="line.139"></a>
143 <FONT color="green">140</FONT> * @param observations sampled observations<a name="line.140"></a>
144 <FONT color="green">141</FONT> */<a name="line.141"></a>
145 <FONT color="green">142</FONT> public HarmonicCoefficientsGuesser(WeightedObservedPoint[] observations) {<a name="line.142"></a>
146 <FONT color="green">143</FONT> this.observations = observations.clone();<a name="line.143"></a>
147 <FONT color="green">144</FONT> a = Double.NaN;<a name="line.144"></a>
148 <FONT color="green">145</FONT> omega = Double.NaN;<a name="line.145"></a>
149 <FONT color="green">146</FONT> }<a name="line.146"></a>
150 <FONT color="green">147</FONT> <a name="line.147"></a>
151 <FONT color="green">148</FONT> /** Estimate a first guess of the coefficients.<a name="line.148"></a>
152 <FONT color="green">149</FONT> * @exception OptimizationException if the sample is too short or if<a name="line.149"></a>
153 <FONT color="green">150</FONT> * the first guess cannot be computed (when the elements under the<a name="line.150"></a>
154 <FONT color="green">151</FONT> * square roots are negative).<a name="line.151"></a>
155 <FONT color="green">152</FONT> * */<a name="line.152"></a>
156 <FONT color="green">153</FONT> public void guess() throws OptimizationException {<a name="line.153"></a>
157 <FONT color="green">154</FONT> sortObservations();<a name="line.154"></a>
158 <FONT color="green">155</FONT> guessAOmega();<a name="line.155"></a>
159 <FONT color="green">156</FONT> guessPhi();<a name="line.156"></a>
160 <FONT color="green">157</FONT> }<a name="line.157"></a>
161 <FONT color="green">158</FONT> <a name="line.158"></a>
162 <FONT color="green">159</FONT> /** Sort the observations with respect to the abscissa.<a name="line.159"></a>
163 <FONT color="green">160</FONT> */<a name="line.160"></a>
164 <FONT color="green">161</FONT> private void sortObservations() {<a name="line.161"></a>
165 <FONT color="green">162</FONT> <a name="line.162"></a>
166 <FONT color="green">163</FONT> // Since the samples are almost always already sorted, this<a name="line.163"></a>
167 <FONT color="green">164</FONT> // method is implemented as an insertion sort that reorders the<a name="line.164"></a>
168 <FONT color="green">165</FONT> // elements in place. Insertion sort is very efficient in this case.<a name="line.165"></a>
169 <FONT color="green">166</FONT> WeightedObservedPoint curr = observations[0];<a name="line.166"></a>
170 <FONT color="green">167</FONT> for (int j = 1; j &lt; observations.length; ++j) {<a name="line.167"></a>
171 <FONT color="green">168</FONT> WeightedObservedPoint prec = curr;<a name="line.168"></a>
172 <FONT color="green">169</FONT> curr = observations[j];<a name="line.169"></a>
173 <FONT color="green">170</FONT> if (curr.getX() &lt; prec.getX()) {<a name="line.170"></a>
174 <FONT color="green">171</FONT> // the current element should be inserted closer to the beginning<a name="line.171"></a>
175 <FONT color="green">172</FONT> int i = j - 1;<a name="line.172"></a>
176 <FONT color="green">173</FONT> WeightedObservedPoint mI = observations[i];<a name="line.173"></a>
177 <FONT color="green">174</FONT> while ((i &gt;= 0) &amp;&amp; (curr.getX() &lt; mI.getX())) {<a name="line.174"></a>
178 <FONT color="green">175</FONT> observations[i + 1] = mI;<a name="line.175"></a>
179 <FONT color="green">176</FONT> if (i-- != 0) {<a name="line.176"></a>
180 <FONT color="green">177</FONT> mI = observations[i];<a name="line.177"></a>
181 <FONT color="green">178</FONT> } else {<a name="line.178"></a>
182 <FONT color="green">179</FONT> mI = null;<a name="line.179"></a>
183 <FONT color="green">180</FONT> }<a name="line.180"></a>
184 <FONT color="green">181</FONT> }<a name="line.181"></a>
185 <FONT color="green">182</FONT> observations[i + 1] = curr;<a name="line.182"></a>
186 <FONT color="green">183</FONT> curr = observations[j];<a name="line.183"></a>
187 <FONT color="green">184</FONT> }<a name="line.184"></a>
188 <FONT color="green">185</FONT> }<a name="line.185"></a>
189 <FONT color="green">186</FONT> <a name="line.186"></a>
190 <FONT color="green">187</FONT> }<a name="line.187"></a>
191 <FONT color="green">188</FONT> <a name="line.188"></a>
192 <FONT color="green">189</FONT> /** Estimate a first guess of the a and &amp;omega; coefficients.<a name="line.189"></a>
193 <FONT color="green">190</FONT> * @exception OptimizationException if the sample is too short or if<a name="line.190"></a>
194 <FONT color="green">191</FONT> * the first guess cannot be computed (when the elements under the<a name="line.191"></a>
195 <FONT color="green">192</FONT> * square roots are negative).<a name="line.192"></a>
196 <FONT color="green">193</FONT> */<a name="line.193"></a>
197 <FONT color="green">194</FONT> private void guessAOmega() throws OptimizationException {<a name="line.194"></a>
198 <FONT color="green">195</FONT> <a name="line.195"></a>
199 <FONT color="green">196</FONT> // initialize the sums for the linear model between the two integrals<a name="line.196"></a>
200 <FONT color="green">197</FONT> double sx2 = 0.0;<a name="line.197"></a>
201 <FONT color="green">198</FONT> double sy2 = 0.0;<a name="line.198"></a>
202 <FONT color="green">199</FONT> double sxy = 0.0;<a name="line.199"></a>
203 <FONT color="green">200</FONT> double sxz = 0.0;<a name="line.200"></a>
204 <FONT color="green">201</FONT> double syz = 0.0;<a name="line.201"></a>
205 <FONT color="green">202</FONT> <a name="line.202"></a>
206 <FONT color="green">203</FONT> double currentX = observations[0].getX();<a name="line.203"></a>
207 <FONT color="green">204</FONT> double currentY = observations[0].getY();<a name="line.204"></a>
208 <FONT color="green">205</FONT> double f2Integral = 0;<a name="line.205"></a>
209 <FONT color="green">206</FONT> double fPrime2Integral = 0;<a name="line.206"></a>
210 <FONT color="green">207</FONT> final double startX = currentX;<a name="line.207"></a>
211 <FONT color="green">208</FONT> for (int i = 1; i &lt; observations.length; ++i) {<a name="line.208"></a>
212 <FONT color="green">209</FONT> <a name="line.209"></a>
213 <FONT color="green">210</FONT> // one step forward<a name="line.210"></a>
214 <FONT color="green">211</FONT> final double previousX = currentX;<a name="line.211"></a>
215 <FONT color="green">212</FONT> final double previousY = currentY;<a name="line.212"></a>
216 <FONT color="green">213</FONT> currentX = observations[i].getX();<a name="line.213"></a>
217 <FONT color="green">214</FONT> currentY = observations[i].getY();<a name="line.214"></a>
218 <FONT color="green">215</FONT> <a name="line.215"></a>
219 <FONT color="green">216</FONT> // update the integrals of f&lt;sup&gt;2&lt;/sup&gt; and f'&lt;sup&gt;2&lt;/sup&gt;<a name="line.216"></a>
220 <FONT color="green">217</FONT> // considering a linear model for f (and therefore constant f')<a name="line.217"></a>
221 <FONT color="green">218</FONT> final double dx = currentX - previousX;<a name="line.218"></a>
222 <FONT color="green">219</FONT> final double dy = currentY - previousY;<a name="line.219"></a>
223 <FONT color="green">220</FONT> final double f2StepIntegral =<a name="line.220"></a>
224 <FONT color="green">221</FONT> dx * (previousY * previousY + previousY * currentY + currentY * currentY) / 3;<a name="line.221"></a>
225 <FONT color="green">222</FONT> final double fPrime2StepIntegral = dy * dy / dx;<a name="line.222"></a>
226 <FONT color="green">223</FONT> <a name="line.223"></a>
227 <FONT color="green">224</FONT> final double x = currentX - startX;<a name="line.224"></a>
228 <FONT color="green">225</FONT> f2Integral += f2StepIntegral;<a name="line.225"></a>
229 <FONT color="green">226</FONT> fPrime2Integral += fPrime2StepIntegral;<a name="line.226"></a>
230 <FONT color="green">227</FONT> <a name="line.227"></a>
231 <FONT color="green">228</FONT> sx2 += x * x;<a name="line.228"></a>
232 <FONT color="green">229</FONT> sy2 += f2Integral * f2Integral;<a name="line.229"></a>
233 <FONT color="green">230</FONT> sxy += x * f2Integral;<a name="line.230"></a>
234 <FONT color="green">231</FONT> sxz += x * fPrime2Integral;<a name="line.231"></a>
235 <FONT color="green">232</FONT> syz += f2Integral * fPrime2Integral;<a name="line.232"></a>
236 <FONT color="green">233</FONT> <a name="line.233"></a>
237 <FONT color="green">234</FONT> }<a name="line.234"></a>
238 <FONT color="green">235</FONT> <a name="line.235"></a>
239 <FONT color="green">236</FONT> // compute the amplitude and pulsation coefficients<a name="line.236"></a>
240 <FONT color="green">237</FONT> double c1 = sy2 * sxz - sxy * syz;<a name="line.237"></a>
241 <FONT color="green">238</FONT> double c2 = sxy * sxz - sx2 * syz;<a name="line.238"></a>
242 <FONT color="green">239</FONT> double c3 = sx2 * sy2 - sxy * sxy;<a name="line.239"></a>
243 <FONT color="green">240</FONT> if ((c1 / c2 &lt; 0.0) || (c2 / c3 &lt; 0.0)) {<a name="line.240"></a>
244 <FONT color="green">241</FONT> throw new OptimizationException("unable to first guess the harmonic coefficients");<a name="line.241"></a>
245 <FONT color="green">242</FONT> }<a name="line.242"></a>
246 <FONT color="green">243</FONT> a = Math.sqrt(c1 / c2);<a name="line.243"></a>
247 <FONT color="green">244</FONT> omega = Math.sqrt(c2 / c3);<a name="line.244"></a>
248 <FONT color="green">245</FONT> <a name="line.245"></a>
249 <FONT color="green">246</FONT> }<a name="line.246"></a>
250 <FONT color="green">247</FONT> <a name="line.247"></a>
251 <FONT color="green">248</FONT> /** Estimate a first guess of the &amp;phi; coefficient.<a name="line.248"></a>
252 <FONT color="green">249</FONT> */<a name="line.249"></a>
253 <FONT color="green">250</FONT> private void guessPhi() {<a name="line.250"></a>
254 <FONT color="green">251</FONT> <a name="line.251"></a>
255 <FONT color="green">252</FONT> // initialize the means<a name="line.252"></a>
256 <FONT color="green">253</FONT> double fcMean = 0.0;<a name="line.253"></a>
257 <FONT color="green">254</FONT> double fsMean = 0.0;<a name="line.254"></a>
258 <FONT color="green">255</FONT> <a name="line.255"></a>
259 <FONT color="green">256</FONT> double currentX = observations[0].getX();<a name="line.256"></a>
260 <FONT color="green">257</FONT> double currentY = observations[0].getY();<a name="line.257"></a>
261 <FONT color="green">258</FONT> for (int i = 1; i &lt; observations.length; ++i) {<a name="line.258"></a>
262 <FONT color="green">259</FONT> <a name="line.259"></a>
263 <FONT color="green">260</FONT> // one step forward<a name="line.260"></a>
264 <FONT color="green">261</FONT> final double previousX = currentX;<a name="line.261"></a>
265 <FONT color="green">262</FONT> final double previousY = currentY;<a name="line.262"></a>
266 <FONT color="green">263</FONT> currentX = observations[i].getX();<a name="line.263"></a>
267 <FONT color="green">264</FONT> currentY = observations[i].getY();<a name="line.264"></a>
268 <FONT color="green">265</FONT> final double currentYPrime = (currentY - previousY) / (currentX - previousX);<a name="line.265"></a>
269 <FONT color="green">266</FONT> <a name="line.266"></a>
270 <FONT color="green">267</FONT> double omegaX = omega * currentX;<a name="line.267"></a>
271 <FONT color="green">268</FONT> double cosine = Math.cos(omegaX);<a name="line.268"></a>
272 <FONT color="green">269</FONT> double sine = Math.sin(omegaX);<a name="line.269"></a>
273 <FONT color="green">270</FONT> fcMean += omega * currentY * cosine - currentYPrime * sine;<a name="line.270"></a>
274 <FONT color="green">271</FONT> fsMean += omega * currentY * sine + currentYPrime * cosine;<a name="line.271"></a>
275 <FONT color="green">272</FONT> <a name="line.272"></a>
276 <FONT color="green">273</FONT> }<a name="line.273"></a>
277 <FONT color="green">274</FONT> <a name="line.274"></a>
278 <FONT color="green">275</FONT> phi = Math.atan2(-fsMean, fcMean);<a name="line.275"></a>
279 <FONT color="green">276</FONT> <a name="line.276"></a>
280 <FONT color="green">277</FONT> }<a name="line.277"></a>
281 <FONT color="green">278</FONT> <a name="line.278"></a>
282 <FONT color="green">279</FONT> /** Get the guessed amplitude a.<a name="line.279"></a>
283 <FONT color="green">280</FONT> * @return guessed amplitude a;<a name="line.280"></a>
284 <FONT color="green">281</FONT> */<a name="line.281"></a>
285 <FONT color="green">282</FONT> public double getGuessedAmplitude() {<a name="line.282"></a>
286 <FONT color="green">283</FONT> return a;<a name="line.283"></a>
287 <FONT color="green">284</FONT> }<a name="line.284"></a>
288 <FONT color="green">285</FONT> <a name="line.285"></a>
289 <FONT color="green">286</FONT> /** Get the guessed pulsation &amp;omega;.<a name="line.286"></a>
290 <FONT color="green">287</FONT> * @return guessed pulsation &amp;omega;<a name="line.287"></a>
291 <FONT color="green">288</FONT> */<a name="line.288"></a>
292 <FONT color="green">289</FONT> public double getGuessedPulsation() {<a name="line.289"></a>
293 <FONT color="green">290</FONT> return omega;<a name="line.290"></a>
294 <FONT color="green">291</FONT> }<a name="line.291"></a>
295 <FONT color="green">292</FONT> <a name="line.292"></a>
296 <FONT color="green">293</FONT> /** Get the guessed phase &amp;phi;.<a name="line.293"></a>
297 <FONT color="green">294</FONT> * @return guessed phase &amp;phi;<a name="line.294"></a>
298 <FONT color="green">295</FONT> */<a name="line.295"></a>
299 <FONT color="green">296</FONT> public double getGuessedPhase() {<a name="line.296"></a>
300 <FONT color="green">297</FONT> return phi;<a name="line.297"></a>
301 <FONT color="green">298</FONT> }<a name="line.298"></a>
302 <FONT color="green">299</FONT> <a name="line.299"></a>
303 <FONT color="green">300</FONT> }<a name="line.300"></a>
304
305
306
307
308
309
310
311
312
313
314
315
316
317
318
319
320
321
322
323
324
325
326
327
328
329
330
331
332
333
334
335
336
337
338
339
340
341
342
343
344
345
346
347
348
349
350
351
352
353
354
355
356
357
358
359
360
361
362
363
364 </PRE>
365 </BODY>
366 </HTML>