Qucs-core  0.0.19
interpolator.cpp
Go to the documentation of this file.
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