Branch data Line data Source code
1 : : /*
2 : : * interpolator.cpp - interpolator class implementation
3 : : *
4 : : * Copyright (C) 2009 Stefan Jahn <stefan@lkcc.org>
5 : : *
6 : : * This is free software; you can redistribute it and/or modify
7 : : * it under the terms of the GNU General Public License as published by
8 : : * the Free Software Foundation; either version 2, or (at your option)
9 : : * any later version.
10 : : *
11 : : * This software is distributed in the hope that it will be useful,
12 : : * but WITHOUT ANY WARRANTY; without even the implied warranty of
13 : : * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
14 : : * GNU General Public License for more details.
15 : : *
16 : : * You should have received a copy of the GNU General Public License
17 : : * along with this package; see the file COPYING. If not, write to
18 : : * the Free Software Foundation, Inc., 51 Franklin Street - Fifth Floor,
19 : : * Boston, MA 02110-1301, USA.
20 : : *
21 : : * $Id$
22 : : *
23 : : */
24 : :
25 : : #if HAVE_CONFIG_H
26 : : # include <config.h>
27 : : #endif
28 : :
29 : : #include <stdio.h>
30 : : #include <stdlib.h>
31 : : #include <string.h>
32 : : #include <assert.h>
33 : :
34 : : #include "poly.h"
35 : : #include "spline.h"
36 : : #include "object.h"
37 : : #include "vector.h"
38 : : #include "interpolator.h"
39 : :
40 : : namespace qucs {
41 : :
42 : : // Constructor creates an instance of the interpolator class.
43 : 0 : interpolator::interpolator () {
44 : 0 : rsp = isp = NULL;
45 : 0 : rx = ry = NULL;
46 : 0 : cy = NULL;
47 : 0 : repeat = dataType = interpolType = length = 0;
48 : 0 : duration = 0.0;
49 : 0 : }
50 : :
51 : :
52 : : // Destructor deletes an instance of the interpolator class.
53 : 0 : interpolator::~interpolator () {
54 [ # # ][ # # ]: 0 : if (rsp) delete rsp;
[ # # ][ # # ]
55 [ # # ][ # # ]: 0 : if (isp) delete isp;
[ # # ][ # # ]
56 [ # # ][ # # ]: 0 : if (rx) free (rx);
57 [ # # ][ # # ]: 0 : if (ry) free (ry);
58 [ # # ][ # # ]: 0 : if (cy) free (cy);
59 : 0 : }
60 : :
61 : : // Cleans up vector storage.
62 : 0 : void interpolator::cleanup (void) {
63 [ # # ]: 0 : if (rx) { free (rx); rx = NULL; }
64 [ # # ]: 0 : if (ry) { free (ry); ry = NULL; }
65 [ # # ]: 0 : if (cy) { free (cy); cy = NULL; }
66 : 0 : }
67 : :
68 : : // Pass real interpolation datapoints as pointers.
69 : 0 : void interpolator::vectors (nr_double_t * y, nr_double_t * x, int len) {
70 : 0 : int len1 = len;
71 : 0 : int len2 = 2 + len * sizeof (nr_double_t);
72 [ # # ]: 0 : if (len > 0) {
73 : 0 : ry = (nr_double_t *) malloc (len2 * sizeof (nr_double_t));
74 : 0 : memcpy (ry, y, len1 * sizeof (nr_double_t));
75 : : }
76 [ # # ]: 0 : if (len > 0) {
77 : 0 : rx = (nr_double_t *) malloc (len2 * sizeof (nr_double_t));
78 : 0 : memcpy (rx, x, len1 * sizeof (nr_double_t));
79 : : }
80 : :
81 : 0 : dataType = (DATA_REAL & DATA_MASK_TYPE);
82 : 0 : length = len;
83 : 0 : }
84 : :
85 : : // Pass real interpolation datapoints as vectors.
86 : 0 : void interpolator::rvectors (vector * y, vector * x) {
87 : 0 : int len = y->getSize ();
88 : 0 : int len1 = len;
89 : 0 : int len2 = 2 + len * sizeof (nr_double_t);
90 : 0 : cleanup ();
91 [ # # ]: 0 : if (len > 0) {
92 : 0 : ry = (nr_double_t *) malloc (len2 * sizeof (nr_double_t));
93 [ # # ]: 0 : for (int i = 0; i < len1; i++) ry[i] = real (y->get (i));
94 : : }
95 [ # # ]: 0 : if (len > 0) {
96 : 0 : rx = (nr_double_t *) malloc (len2 * sizeof (nr_double_t));
97 [ # # ]: 0 : for (int i = 0; i < len1; i++) rx[i] = real (x->get (i));
98 : : }
99 : :
100 : 0 : dataType = (DATA_REAL & DATA_MASK_TYPE);
101 : 0 : length = len;
102 : 0 : }
103 : :
104 : : // Pass complex interpolation datapoints as pointers.
105 : 0 : void interpolator::vectors (nr_complex_t * y, nr_double_t * x, int len) {
106 : 0 : int len1 = len;
107 : 0 : int len2 = 2 + len;
108 : 0 : cleanup ();
109 [ # # ]: 0 : if (len > 0) {
110 : 0 : cy = (nr_complex_t *) malloc (len2 * sizeof (nr_complex_t));
111 : 0 : memcpy (cy, y, len1 * sizeof (nr_complex_t));
112 : : }
113 [ # # ]: 0 : if (len > 0) {
114 : 0 : rx = (nr_double_t *) malloc (len2 * sizeof (nr_double_t));
115 : 0 : memcpy (rx, x, len1 * sizeof (nr_double_t));
116 : : }
117 : :
118 : 0 : dataType = (DATA_COMPLEX & DATA_MASK_TYPE);
119 : 0 : length = len;
120 : 0 : }
121 : :
122 : : // Pass complex interpolation datapoints as vectors.
123 : 0 : void interpolator::cvectors (vector * y, vector * x) {
124 : 0 : int len = y->getSize ();
125 : 0 : int len1 = len;
126 : 0 : int len2 = 2 + len;
127 : 0 : cleanup ();
128 [ # # ]: 0 : if (len > 0) {
129 : 0 : cy = (nr_complex_t *) malloc (len2 * sizeof (nr_complex_t));
130 [ # # ]: 0 : for (int i = 0; i < len1; i++) cy[i] = y->get (i);
131 : : }
132 [ # # ]: 0 : if (len > 0) {
133 : 0 : rx = (nr_double_t *) malloc (len2 * sizeof (nr_double_t));
134 [ # # ]: 0 : for (int i = 0; i < len1; i++) rx[i] = real (x->get (i));
135 : : }
136 : :
137 : 0 : dataType = (DATA_COMPLEX & DATA_MASK_TYPE);
138 : 0 : length = len;
139 : 0 : }
140 : :
141 : : // Prepares interpolator instance, e.g. setups spline object.
142 : 0 : void interpolator::prepare (int interpol, int repitition, int domain) {
143 : 0 : interpolType = interpol;
144 : 0 : dataType |= (domain & DATA_MASK_DOMAIN);
145 : 0 : repeat = repitition;
146 : :
147 : : // preparations for cyclic interpolations
148 [ # # ]: 0 : if (repeat & REPEAT_YES) {
149 : 0 : duration = rx[length - 1] - rx[0];
150 : : // set first value to the end of the value vector
151 [ # # ]: 0 : if (cy) cy[length - 1] = cy[0];
152 [ # # ]: 0 : if (ry) ry[length - 1] = ry[0];
153 : : }
154 : :
155 : : // preparations for polar complex data
156 [ # # ][ # # ]: 0 : if (cy != NULL && (domain & DATA_POLAR) && length > 1) {
[ # # ]
157 : : // unwrap phase of complex data vector
158 [ # # ]: 0 : vector ang = vector (length);
159 [ # # ]: 0 : for (int i = 0; i < length; i++) ang (i) = arg (cy[i]);
160 [ # # ][ # # ]: 0 : ang = unwrap (ang);
[ # # ][ # # ]
[ # # ]
161 : : // store complex data
162 [ # # ]: 0 : for (int i = 0; i < length; i++) {
163 : 0 : cy[i] = nr_complex_t (abs (cy[i]), real (ang (i)));
164 [ # # ]: 0 : }
165 : : }
166 : :
167 : : // preparations spline interpolations
168 [ # # ]: 0 : if (interpolType & INTERPOL_CUBIC) {
169 : :
170 : : // prepare complex vector interpolation using splines
171 [ # # ]: 0 : if (cy != NULL) {
172 : : // create splines if necessary
173 [ # # ][ # # ]: 0 : if (rsp) delete rsp;
[ # # ]
174 [ # # ][ # # ]: 0 : if (isp) delete isp;
[ # # ]
175 [ # # ][ # # ]: 0 : rsp = new spline (SPLINE_BC_NATURAL);
176 [ # # ][ # # ]: 0 : isp = new spline (SPLINE_BC_NATURAL);
177 [ # # ]: 0 : if (repeat & REPEAT_YES) {
178 : 0 : rsp->setBoundary (SPLINE_BC_PERIODIC);
179 : 0 : isp->setBoundary (SPLINE_BC_PERIODIC);
180 : : }
181 : : // prepare data vectors
182 [ # # ]: 0 : vector rv = vector (length);
183 [ # # ]: 0 : vector iv = vector (length);
184 [ # # ]: 0 : vector rt = vector (length);
185 [ # # ]: 0 : for (int i = 0; i < length; i++) {
186 : 0 : rv (i) = real (cy[i]);
187 : 0 : iv (i) = imag (cy[i]);
188 : 0 : rt (i) = rx[i];
189 : : }
190 : : // pass data vectors to splines and construct these
191 [ # # ][ # # ]: 0 : rsp->vectors (rv, rt);
[ # # ][ # # ]
[ # # ]
192 [ # # ][ # # ]: 0 : isp->vectors (iv, rt);
[ # # ][ # # ]
[ # # ]
193 [ # # ]: 0 : rsp->construct ();
194 [ # # ][ # # ]: 0 : isp->construct ();
[ # # ][ # # ]
195 : : }
196 : :
197 : : // prepare real vector interpolation using spline
198 : : else {
199 [ # # ][ # # ]: 0 : if (rsp) delete rsp;
200 [ # # ]: 0 : rsp = new spline (SPLINE_BC_NATURAL);
201 [ # # ]: 0 : if (repeat & REPEAT_YES) rsp->setBoundary (SPLINE_BC_PERIODIC);
202 : 0 : rsp->vectors (ry, rx, length);
203 : 0 : rsp->construct ();
204 : : }
205 : : }
206 : 0 : }
207 : :
208 : : /* The function performs a binary search on the ascending sorted
209 : : x-vector in order to return the left-hand-side index pointer into
210 : : the x-vector based on the given value. */
211 : 0 : int interpolator::findIndex (nr_double_t x) {
212 : 0 : int lo = 0;
213 : 0 : int hi = length;
214 : : int av;
215 [ # # ]: 0 : while (lo < hi) {
216 : 0 : av = lo + ((hi - lo) / 2);
217 [ # # ]: 0 : if (x >= rx[av])
218 : 0 : lo = av + 1;
219 : : else
220 : : // can't be hi = av-1: here rx[av] >= x,
221 : : // so hi can't be < av if rx[av] == x
222 : 0 : hi = av;
223 : : }
224 : : // hi == lo, using hi or lo depends on taste
225 [ # # ][ # # ]: 0 : if (lo <= length && lo > 0 && x >= rx[lo - 1])
[ # # ]
226 : 0 : return lo - 1; // found
227 : : else
228 : 0 : return 0; // not found
229 : : }
230 : :
231 : : /* Return the left-hand-side index pointer into the x-vector based on
232 : : the given value. Returns 0 or 'length' if the value is beyond the
233 : : x-vectors scope. */
234 : 0 : int interpolator::findIndex_old (nr_double_t x) {
235 : 0 : int idx = 0;
236 [ # # ]: 0 : for (int i = 0; i < length; i++) {
237 [ # # ]: 0 : if (x >= rx[i]) idx = i;
238 : : }
239 : 0 : return idx;
240 : : }
241 : :
242 : : /* Computes simple linear interpolation value for the given values. */
243 : 0 : nr_double_t interpolator::linear (nr_double_t x,
244 : : nr_double_t x1, nr_double_t x2,
245 : : nr_double_t y1, nr_double_t y2) {
246 [ # # ]: 0 : if (x1 == x2)
247 : 0 : return (y1 + y2) / 2;
248 : : else
249 : 0 : return ((x2 - x) * y1 + (x - x1) * y2) / (x2 - x1);
250 : : }
251 : :
252 : : /* Returns real valued linear interpolation. */
253 : 0 : nr_double_t interpolator::rlinear (nr_double_t x, int idx) {
254 : 0 : return linear (x, rx[idx], rx[idx+1], ry[idx], ry[idx+1]);
255 : : }
256 : :
257 : : /* Returns complex valued linear interpolation. */
258 : 0 : nr_complex_t interpolator::clinear (nr_double_t x, int idx) {
259 : : nr_double_t x1, x2, r, i;
260 : 0 : nr_complex_t y1, y2;
261 : 0 : x1 = rx[idx]; x2 = rx[idx+1];
262 : 0 : y1 = cy[idx]; y2 = cy[idx+1];
263 : 0 : r = linear (x, x1, x2, real (y1), real (y2));
264 : 0 : i = linear (x, x1, x2, imag (y1), imag (y2));
265 : 0 : return nr_complex_t (r, i);
266 : : }
267 : :
268 : : /* This function interpolates for real values. Returns the linear
269 : : interpolation of the real y-vector for the given value in the
270 : : x-vector. */
271 : 0 : nr_double_t interpolator::rinterpolate (nr_double_t x) {
272 : 0 : int idx = -1;
273 : 0 : nr_double_t res = 0.0;
274 : :
275 : : // no chance to interpolate
276 [ # # ]: 0 : if (length <= 0) {
277 : 0 : return res;
278 : : }
279 : : // no interpolation necessary
280 [ # # ]: 0 : else if (length == 1) {
281 : 0 : res = ry[0];
282 : 0 : return res;
283 : : }
284 [ # # ]: 0 : else if (repeat & REPEAT_YES)
285 : 0 : x = x - std::floor (x / duration) * duration;
286 : :
287 : : // linear interpolation
288 [ # # ]: 0 : if (interpolType & INTERPOL_LINEAR) {
289 : 0 : idx = findIndex (x);
290 : : // dependency variable in scope or beyond
291 [ # # ]: 0 : if (x == rx[idx])
292 : 0 : res = ry[idx];
293 : : // dependency variable is beyond scope; use last tangent
294 : : else {
295 [ # # ]: 0 : if (idx == length - 1) idx--;
296 : 0 : res = rlinear (x, idx);
297 : : }
298 : : }
299 : : // cubic spline interpolation
300 [ # # ]: 0 : else if (interpolType & INTERPOL_CUBIC) {
301 : : // evaluate spline functions
302 : 0 : res = rsp->evaluate (x).f0;
303 : : }
304 [ # # ]: 0 : else if (interpolType & INTERPOL_HOLD) {
305 : : // find appropriate dependency index
306 : 0 : idx = findIndex (x);
307 : 0 : res = ry[idx];
308 : : }
309 : 0 : return res;
310 : : }
311 : :
312 : : /* This function interpolates for complex values. Returns the complex
313 : : interpolation of the real y-vector for the given value in the
314 : : x-vector. */
315 : 0 : nr_complex_t interpolator::cinterpolate (nr_double_t x) {
316 : 0 : int idx = -1;
317 : 0 : nr_complex_t res = 0.0;
318 : :
319 : : // no chance to interpolate
320 [ # # ]: 0 : if (length <= 0) {
321 : 0 : return res;
322 : : }
323 : : // no interpolation necessary
324 [ # # ]: 0 : else if (length == 1) {
325 : 0 : res = cy[0];
326 : 0 : return res;
327 : : }
328 [ # # ]: 0 : else if (repeat & REPEAT_YES)
329 : 0 : x = x - std::floor (x / duration) * duration;
330 : :
331 : : // linear interpolation
332 [ # # ]: 0 : if (interpolType & INTERPOL_LINEAR) {
333 : 0 : idx = findIndex (x);
334 : : // dependency variable in scope or beyond
335 [ # # ]: 0 : if (x == rx[idx])
336 : 0 : res = cy[idx];
337 : : // dependency variable is beyond scope; use last tangent
338 : : else {
339 [ # # ]: 0 : if (idx == length - 1) idx--;
340 : 0 : res = clinear (x, idx);
341 : : }
342 : : }
343 : : // cubic spline interpolation
344 [ # # ]: 0 : else if (interpolType & INTERPOL_CUBIC) {
345 : : // evaluate spline functions
346 [ # # ]: 0 : nr_double_t r = rsp->evaluate (x).f0;
347 [ # # ]: 0 : nr_double_t i = isp->evaluate (x).f0;
348 : 0 : res = nr_complex_t (r, i);
349 : : }
350 [ # # ]: 0 : else if (interpolType & INTERPOL_HOLD) {
351 : : // find appropriate dependency index
352 : 0 : idx = findIndex (x);
353 : 0 : res = cy[idx];
354 : : }
355 : :
356 : : // depending on the data type return appropriate complex value
357 [ # # ]: 0 : if (dataType & DATA_POLAR)
358 : 0 : return std::polar (real (res), imag (res));
359 : : else
360 : 0 : return res;
361 : : }
362 : :
363 : : } // namespace qucs
|
