Qucs-core
0.0.19
|
00001 /* 00002 * interpolator.cpp - interpolator class implementation 00003 * 00004 * Copyright (C) 2009 Stefan Jahn <stefan@lkcc.org> 00005 * 00006 * This is free software; you can redistribute it and/or modify 00007 * it under the terms of the GNU General Public License as published by 00008 * the Free Software Foundation; either version 2, or (at your option) 00009 * any later version. 00010 * 00011 * This software is distributed in the hope that it will be useful, 00012 * but WITHOUT ANY WARRANTY; without even the implied warranty of 00013 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the 00014 * GNU General Public License for more details. 00015 * 00016 * You should have received a copy of the GNU General Public License 00017 * along with this package; see the file COPYING. If not, write to 00018 * the Free Software Foundation, Inc., 51 Franklin Street - Fifth Floor, 00019 * Boston, MA 02110-1301, USA. 00020 * 00021 * $Id$ 00022 * 00023 */ 00024 00025 #if HAVE_CONFIG_H 00026 # include <config.h> 00027 #endif 00028 00029 #include <stdio.h> 00030 #include <stdlib.h> 00031 #include <string.h> 00032 #include <assert.h> 00033 00034 #include "poly.h" 00035 #include "spline.h" 00036 #include "object.h" 00037 #include "vector.h" 00038 #include "interpolator.h" 00039 00040 namespace qucs { 00041 00042 // Constructor creates an instance of the interpolator class. 00043 interpolator::interpolator () { 00044 rsp = isp = NULL; 00045 rx = ry = NULL; 00046 cy = NULL; 00047 repeat = dataType = interpolType = length = 0; 00048 duration = 0.0; 00049 } 00050 00051 00052 // Destructor deletes an instance of the interpolator class. 00053 interpolator::~interpolator () { 00054 delete rsp; 00055 delete isp; 00056 free (rx); 00057 free (ry); 00058 free (cy); 00059 } 00060 00061 // Cleans up vector storage. 00062 void interpolator::cleanup (void) { 00063 if (rx) { free (rx); rx = NULL; } 00064 if (ry) { free (ry); ry = NULL; } 00065 if (cy) { free (cy); cy = NULL; } 00066 } 00067 00068 // Pass real interpolation datapoints as pointers. 00069 void interpolator::vectors (nr_double_t * y, nr_double_t * x, int len) { 00070 int len1 = len; 00071 int len2 = 2 + len * sizeof (nr_double_t); 00072 if (len > 0) { 00073 ry = (nr_double_t *) malloc (len2 * sizeof (nr_double_t)); 00074 memcpy (ry, y, len1 * sizeof (nr_double_t)); 00075 } 00076 if (len > 0) { 00077 rx = (nr_double_t *) malloc (len2 * sizeof (nr_double_t)); 00078 memcpy (rx, x, len1 * sizeof (nr_double_t)); 00079 } 00080 00081 dataType = (DATA_REAL & DATA_MASK_TYPE); 00082 length = len; 00083 } 00084 00085 // Pass real interpolation datapoints as vectors. 00086 void interpolator::rvectors (vector * y, vector * x) { 00087 int len = y->getSize (); 00088 int len1 = len; 00089 int len2 = 2 + len * sizeof (nr_double_t); 00090 cleanup (); 00091 if (len > 0) { 00092 ry = (nr_double_t *) malloc (len2 * sizeof (nr_double_t)); 00093 for (int i = 0; i < len1; i++) ry[i] = real (y->get (i)); 00094 } 00095 if (len > 0) { 00096 rx = (nr_double_t *) malloc (len2 * sizeof (nr_double_t)); 00097 for (int i = 0; i < len1; i++) rx[i] = real (x->get (i)); 00098 } 00099 00100 dataType = (DATA_REAL & DATA_MASK_TYPE); 00101 length = len; 00102 } 00103 00104 // Pass complex interpolation datapoints as pointers. 00105 void interpolator::vectors (nr_complex_t * y, nr_double_t * x, int len) { 00106 int len1 = len; 00107 int len2 = 2 + len; 00108 cleanup (); 00109 if (len > 0) { 00110 cy = (nr_complex_t *) malloc (len2 * sizeof (nr_complex_t)); 00111 memcpy (cy, y, len1 * sizeof (nr_complex_t)); 00112 } 00113 if (len > 0) { 00114 rx = (nr_double_t *) malloc (len2 * sizeof (nr_double_t)); 00115 memcpy (rx, x, len1 * sizeof (nr_double_t)); 00116 } 00117 00118 dataType = (DATA_COMPLEX & DATA_MASK_TYPE); 00119 length = len; 00120 } 00121 00122 // Pass complex interpolation datapoints as vectors. 00123 void interpolator::cvectors (vector * y, vector * x) { 00124 int len = y->getSize (); 00125 int len1 = len; 00126 int len2 = 2 + len; 00127 cleanup (); 00128 if (len > 0) { 00129 cy = (nr_complex_t *) malloc (len2 * sizeof (nr_complex_t)); 00130 for (int i = 0; i < len1; i++) cy[i] = y->get (i); 00131 } 00132 if (len > 0) { 00133 rx = (nr_double_t *) malloc (len2 * sizeof (nr_double_t)); 00134 for (int i = 0; i < len1; i++) rx[i] = real (x->get (i)); 00135 } 00136 00137 dataType = (DATA_COMPLEX & DATA_MASK_TYPE); 00138 length = len; 00139 } 00140 00141 // Prepares interpolator instance, e.g. setups spline object. 00142 void interpolator::prepare (int interpol, int repitition, int domain) { 00143 interpolType = interpol; 00144 dataType |= (domain & DATA_MASK_DOMAIN); 00145 repeat = repitition; 00146 00147 // preparations for cyclic interpolations 00148 if (repeat & REPEAT_YES) { 00149 duration = rx[length - 1] - rx[0]; 00150 // set first value to the end of the value vector 00151 if (cy) cy[length - 1] = cy[0]; 00152 if (ry) ry[length - 1] = ry[0]; 00153 } 00154 00155 // preparations for polar complex data 00156 if (cy != NULL && (domain & DATA_POLAR) && length > 1) { 00157 // unwrap phase of complex data vector 00158 vector ang = vector (length); 00159 for (int i = 0; i < length; i++) ang (i) = arg (cy[i]); 00160 ang = unwrap (ang); 00161 // store complex data 00162 for (int i = 0; i < length; i++) { 00163 cy[i] = nr_complex_t (abs (cy[i]), real (ang (i))); 00164 } 00165 } 00166 00167 // preparations spline interpolations 00168 if (interpolType & INTERPOL_CUBIC) { 00169 00170 // prepare complex vector interpolation using splines 00171 if (cy != NULL) { 00172 // create splines if necessary 00173 delete rsp; 00174 delete isp; 00175 rsp = new spline (SPLINE_BC_NATURAL); 00176 isp = new spline (SPLINE_BC_NATURAL); 00177 if (repeat & REPEAT_YES) { 00178 rsp->setBoundary (SPLINE_BC_PERIODIC); 00179 isp->setBoundary (SPLINE_BC_PERIODIC); 00180 } 00181 // prepare data vectors 00182 vector rv = vector (length); 00183 vector iv = vector (length); 00184 vector rt = vector (length); 00185 for (int i = 0; i < length; i++) { 00186 rv (i) = real (cy[i]); 00187 iv (i) = imag (cy[i]); 00188 rt (i) = rx[i]; 00189 } 00190 // pass data vectors to splines and construct these 00191 rsp->vectors (rv, rt); 00192 isp->vectors (iv, rt); 00193 rsp->construct (); 00194 isp->construct (); 00195 } 00196 00197 // prepare real vector interpolation using spline 00198 else { 00199 delete rsp; 00200 rsp = new spline (SPLINE_BC_NATURAL); 00201 if (repeat & REPEAT_YES) rsp->setBoundary (SPLINE_BC_PERIODIC); 00202 rsp->vectors (ry, rx, length); 00203 rsp->construct (); 00204 } 00205 } 00206 } 00207 00208 /* The function performs a binary search on the ascending sorted 00209 x-vector in order to return the left-hand-side index pointer into 00210 the x-vector based on the given value. */ 00211 int interpolator::findIndex (nr_double_t x) { 00212 int lo = 0; 00213 int hi = length; 00214 int av; 00215 while (lo < hi) { 00216 av = lo + ((hi - lo) / 2); 00217 if (x >= rx[av]) 00218 lo = av + 1; 00219 else 00220 // can't be hi = av-1: here rx[av] >= x, 00221 // so hi can't be < av if rx[av] == x 00222 hi = av; 00223 } 00224 // hi == lo, using hi or lo depends on taste 00225 if (lo <= length && lo > 0 && x >= rx[lo - 1]) 00226 return lo - 1; // found 00227 else 00228 return 0; // not found 00229 } 00230 00231 /* Return the left-hand-side index pointer into the x-vector based on 00232 the given value. Returns 0 or 'length' if the value is beyond the 00233 x-vectors scope. */ 00234 int interpolator::findIndex_old (nr_double_t x) { 00235 int idx = 0; 00236 for (int i = 0; i < length; i++) { 00237 if (x >= rx[i]) idx = i; 00238 } 00239 return idx; 00240 } 00241 00242 /* Computes simple linear interpolation value for the given values. */ 00243 nr_double_t interpolator::linear (nr_double_t x, 00244 nr_double_t x1, nr_double_t x2, 00245 nr_double_t y1, nr_double_t y2) { 00246 if (x1 == x2) 00247 return (y1 + y2) / 2; 00248 else 00249 return ((x2 - x) * y1 + (x - x1) * y2) / (x2 - x1); 00250 } 00251 00252 /* Returns real valued linear interpolation. */ 00253 nr_double_t interpolator::rlinear (nr_double_t x, int idx) { 00254 return linear (x, rx[idx], rx[idx+1], ry[idx], ry[idx+1]); 00255 } 00256 00257 /* Returns complex valued linear interpolation. */ 00258 nr_complex_t interpolator::clinear (nr_double_t x, int idx) { 00259 nr_double_t x1, x2, r, i; 00260 nr_complex_t y1, y2; 00261 x1 = rx[idx]; x2 = rx[idx+1]; 00262 y1 = cy[idx]; y2 = cy[idx+1]; 00263 r = linear (x, x1, x2, real (y1), real (y2)); 00264 i = linear (x, x1, x2, imag (y1), imag (y2)); 00265 return nr_complex_t (r, i); 00266 } 00267 00268 /* This function interpolates for real values. Returns the linear 00269 interpolation of the real y-vector for the given value in the 00270 x-vector. */ 00271 nr_double_t interpolator::rinterpolate (nr_double_t x) { 00272 int idx = -1; 00273 nr_double_t res = 0.0; 00274 00275 // no chance to interpolate 00276 if (length <= 0) { 00277 return res; 00278 } 00279 // no interpolation necessary 00280 else if (length == 1) { 00281 res = ry[0]; 00282 return res; 00283 } 00284 else if (repeat & REPEAT_YES) 00285 x = x - std::floor (x / duration) * duration; 00286 00287 // linear interpolation 00288 if (interpolType & INTERPOL_LINEAR) { 00289 idx = findIndex (x); 00290 // dependency variable in scope or beyond 00291 if (x == rx[idx]) 00292 res = ry[idx]; 00293 // dependency variable is beyond scope; use last tangent 00294 else { 00295 if (idx == length - 1) idx--; 00296 res = rlinear (x, idx); 00297 } 00298 } 00299 // cubic spline interpolation 00300 else if (interpolType & INTERPOL_CUBIC) { 00301 // evaluate spline functions 00302 res = rsp->evaluate (x).f0; 00303 } 00304 else if (interpolType & INTERPOL_HOLD) { 00305 // find appropriate dependency index 00306 idx = findIndex (x); 00307 res = ry[idx]; 00308 } 00309 return res; 00310 } 00311 00312 /* This function interpolates for complex values. Returns the complex 00313 interpolation of the real y-vector for the given value in the 00314 x-vector. */ 00315 nr_complex_t interpolator::cinterpolate (nr_double_t x) { 00316 int idx = -1; 00317 nr_complex_t res = 0.0; 00318 00319 // no chance to interpolate 00320 if (length <= 0) { 00321 return res; 00322 } 00323 // no interpolation necessary 00324 else if (length == 1) { 00325 res = cy[0]; 00326 return res; 00327 } 00328 else if (repeat & REPEAT_YES) 00329 x = x - std::floor (x / duration) * duration; 00330 00331 // linear interpolation 00332 if (interpolType & INTERPOL_LINEAR) { 00333 idx = findIndex (x); 00334 // dependency variable in scope or beyond 00335 if (x == rx[idx]) 00336 res = cy[idx]; 00337 // dependency variable is beyond scope; use last tangent 00338 else { 00339 if (idx == length - 1) idx--; 00340 res = clinear (x, idx); 00341 } 00342 } 00343 // cubic spline interpolation 00344 else if (interpolType & INTERPOL_CUBIC) { 00345 // evaluate spline functions 00346 nr_double_t r = rsp->evaluate (x).f0; 00347 nr_double_t i = isp->evaluate (x).f0; 00348 res = nr_complex_t (r, i); 00349 } 00350 else if (interpolType & INTERPOL_HOLD) { 00351 // find appropriate dependency index 00352 idx = findIndex (x); 00353 res = cy[idx]; 00354 } 00355 00356 // depending on the data type return appropriate complex value 00357 if (dataType & DATA_POLAR) 00358 return std::polar (real (res), imag (res)); 00359 else 00360 return res; 00361 } 00362 00363 } // namespace qucs