Qucs-core
0.0.19
|
00001 /* 00002 * vector.cpp - vector class implementation 00003 * 00004 * Copyright (C) 2003, 2004, 2005, 2006, 2007 Stefan Jahn <stefan@lkcc.org> 00005 * Copyright (C) 2006, 2007 Gunther Kraut <gn.kraut@t-online.de> 00006 * 00007 * This is free software; you can redistribute it and/or modify 00008 * it under the terms of the GNU General Public License as published by 00009 * the Free Software Foundation; either version 2, or (at your option) 00010 * any later version. 00011 * 00012 * This software is distributed in the hope that it will be useful, 00013 * but WITHOUT ANY WARRANTY; without even the implied warranty of 00014 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the 00015 * GNU General Public License for more details. 00016 * 00017 * You should have received a copy of the GNU General Public License 00018 * along with this package; see the file COPYING. If not, write to 00019 * the Free Software Foundation, Inc., 51 Franklin Street - Fifth Floor, 00020 * Boston, MA 02110-1301, USA. 00021 * 00022 * $Id$ 00023 * 00024 */ 00025 00026 #if HAVE_CONFIG_H 00027 # include <config.h> 00028 #endif 00029 00030 #include <limits> 00031 00032 #include <stdio.h> 00033 #include <stdlib.h> 00034 #include <string.h> 00035 #include <cmath> 00036 #include <float.h> 00037 #include <assert.h> 00038 00039 #include "real.h" 00040 #include "complex.h" 00041 #include "object.h" 00042 #include "logging.h" 00043 #include "strlist.h" 00044 #include "vector.h" 00045 #include "consts.h" 00046 00047 namespace qucs { 00048 00049 // Constructor creates an unnamed instance of the vector class. 00050 vector::vector () : object () { 00051 capacity = size = 0; 00052 data = NULL; 00053 dependencies = NULL; 00054 origin = NULL; 00055 requested = 0; 00056 next = prev = nullptr; 00057 } 00058 00059 /* Constructor creates an unnamed instance of the vector class with a 00060 given initial size. */ 00061 vector::vector (int s) : object () { 00062 assert (s >= 0); 00063 capacity = size = s; 00064 data = s > 0 ? (nr_complex_t *) 00065 calloc (capacity, sizeof (nr_complex_t)) : NULL; 00066 dependencies = NULL; 00067 origin = NULL; 00068 requested = 0; 00069 next = prev = nullptr; 00070 } 00071 00072 /* Constructor creates an unnamed instance of the vector class with a 00073 given initial size and content. */ 00074 vector::vector (int s, nr_complex_t val) : object () { 00075 assert (s >= 0); 00076 capacity = size = s; 00077 data = s > 0 ? (nr_complex_t *) 00078 calloc (capacity, sizeof (nr_complex_t)) : NULL; 00079 for (int i = 0; i < s; i++) data[i] = val; 00080 dependencies = NULL; 00081 origin = NULL; 00082 requested = 0; 00083 next = prev = nullptr; 00084 } 00085 00086 // Constructor creates an named instance of the vector class. 00087 vector::vector (const std::string &n) : object (n) { 00088 capacity = size = 0; 00089 data = NULL; 00090 dependencies = NULL; 00091 origin = NULL; 00092 requested = 0; 00093 next = prev = nullptr; 00094 } 00095 00096 /* This constructor creates a named instance of the vector class with 00097 a given initial size. */ 00098 vector::vector (const std::string &n, int s) : object (n) { 00099 assert (s >= 0); 00100 capacity = size = s; 00101 data = s > 0 ? (nr_complex_t *) 00102 calloc (capacity, sizeof (nr_complex_t)) : NULL; 00103 dependencies = NULL; 00104 origin = NULL; 00105 requested = 0; 00106 next = prev = nullptr; 00107 } 00108 00109 /* The copy constructor creates a new instance based on the given 00110 vector object. */ 00111 vector::vector (const vector & v) : object (v) { 00112 size = v.size; 00113 capacity = v.capacity; 00114 data = (nr_complex_t *) malloc (sizeof (nr_complex_t) * capacity); 00115 memcpy (data, v.data, sizeof (nr_complex_t) * size); 00116 dependencies = v.dependencies ? new strlist (*v.dependencies) : NULL; 00117 origin = v.origin ? strdup (v.origin) : NULL; 00118 requested = v.requested; 00119 next = v.next; 00120 prev = v.prev; 00121 } 00122 00123 /* The assignment copy constructor creates a new instance based on the 00124 given vector object. It copies the data only and leaves any other 00125 properties untouched. */ 00126 const vector& vector::operator=(const vector & v) { 00127 if (&v != this) { 00128 size = v.size; 00129 capacity = v.capacity; 00130 if (data) { free (data); data = NULL; } 00131 if (capacity > 0) { 00132 data = (nr_complex_t *) malloc (sizeof (nr_complex_t) * capacity); 00133 if (size > 0) memcpy (data, v.data, sizeof (nr_complex_t) * size); 00134 } 00135 } 00136 return *this; 00137 } 00138 00139 // Destructor deletes a vector object. 00140 vector::~vector () { 00141 free (data); 00142 delete dependencies; 00143 free (origin); 00144 } 00145 00146 // Returns data dependencies. 00147 strlist * vector::getDependencies (void) { 00148 return dependencies; 00149 } 00150 00151 // Sets the data dependencies. 00152 void vector::setDependencies (strlist * s) { 00153 delete dependencies; 00154 dependencies = s; 00155 } 00156 00157 /* The function appends a new complex data item to the end of the 00158 vector and ensures that the vector can hold the increasing number 00159 of data items. */ 00160 void vector::add (nr_complex_t c) { 00161 if (data == NULL) { 00162 size = 0; capacity = 64; 00163 data = (nr_complex_t *) malloc (sizeof (nr_complex_t) * capacity); 00164 } 00165 else if (size >= capacity) { 00166 capacity *= 2; 00167 data = (nr_complex_t *) realloc (data, sizeof (nr_complex_t) * capacity); 00168 } 00169 data[size++] = c; 00170 } 00171 00172 /* This function appends the given vector to the vector. */ 00173 void vector::add (vector * v) { 00174 if (v != NULL) { 00175 if (data == NULL) { 00176 size = 0; capacity = v->getSize (); 00177 data = (nr_complex_t *) malloc (sizeof (nr_complex_t) * capacity); 00178 } 00179 else if (size + v->getSize () > capacity) { 00180 capacity += v->getSize (); 00181 data = (nr_complex_t *) realloc (data, sizeof (nr_complex_t) * capacity); 00182 } 00183 for (int i = 0; i < v->getSize (); i++) data[size++] = v->get (i); 00184 } 00185 } 00186 00187 // Returns the complex data item at the given position. 00188 nr_complex_t vector::get (int i) { 00189 return data[i]; 00190 } 00191 00192 void vector::set (nr_double_t d, int i) { 00193 data[i] = nr_complex_t (d); 00194 } 00195 00196 void vector::set (const nr_complex_t z, int i) { 00197 data[i] = nr_complex_t (z); 00198 } 00199 00200 // The function returns the current size of the vector. 00201 int vector::getSize (void) const { 00202 return size; 00203 } 00204 00205 int vector::checkSizes (vector v1, vector v2) { 00206 if (v1.getSize () != v2.getSize ()) { 00207 logprint (LOG_ERROR, "vector '%s' and '%s' have different sizes\n", 00208 v1.getName (), v2.getName ()); 00209 return 0; 00210 } 00211 return 1; 00212 } 00213 00214 // searches the maximum value of the vector elements. 00215 // complex numbers in the 1. and 4. quadrant are counted as "abs(c)". 00216 // complex numbers in the 2. and 3. quadrant are counted as "-abs(c)". 00217 nr_double_t vector::maximum (void) { 00218 nr_complex_t c; 00219 nr_double_t d, max_D = -std::numeric_limits<nr_double_t>::max(); 00220 for (int i = 0; i < getSize (); i++) { 00221 c = data[i]; 00222 d = fabs (arg (c)) < pi_over_2 ? abs (c) : -abs (c); 00223 if (d > max_D) max_D = d; 00224 } 00225 return max_D; 00226 } 00227 00228 // searches the minimum value of the vector elements. 00229 // complex numbers in the 1. and 4. quadrant are counted as "abs(c)". 00230 // complex numbers in the 2. and 3. quadrant are counted as "-abs(c)". 00231 nr_double_t vector::minimum (void) { 00232 nr_complex_t c; 00233 nr_double_t d, min_D = +std::numeric_limits<nr_double_t>::max(); 00234 for (int i = 0; i < getSize (); i++) { 00235 c = data[i]; 00236 d = fabs (arg (c)) < pi_over_2 ? abs (c) : -abs (c); 00237 if (d < min_D) min_D = d; 00238 } 00239 return min_D; 00240 } 00241 00242 /* Unwraps a phase vector in radians. Adds +/- 2*Pi if consecutive 00243 values jump about |Pi|. */ 00244 vector unwrap (vector v, nr_double_t tol, nr_double_t step) { 00245 vector result (v.getSize ()); 00246 nr_double_t add = 0; 00247 result (0) = v (0); 00248 for (int i = 1; i < v.getSize (); i++) { 00249 nr_double_t diff = real (v (i) - v (i-1)); 00250 if (diff > +tol) { 00251 add -= step; 00252 } else if (diff < -tol) { 00253 add += step; 00254 } 00255 result (i) = v (i) + add; 00256 } 00257 return result; 00258 } 00259 00260 nr_complex_t sum (vector v) { 00261 nr_complex_t result (0.0); 00262 for (int i = 0; i < v.getSize (); i++) result += v.get (i); 00263 return result; 00264 } 00265 00266 nr_complex_t prod (vector v) { 00267 nr_complex_t result (1.0); 00268 for (int i = 0; i < v.getSize (); i++) result *= v.get (i); 00269 return result; 00270 } 00271 00272 nr_complex_t avg (vector v) { 00273 nr_complex_t result (0.0); 00274 for (int i = 0; i < v.getSize (); i++) result += v.get (i); 00275 return result / (nr_double_t) v.getSize (); 00276 } 00277 00278 vector signum (vector v) { 00279 vector result (v); 00280 for (int i = 0; i < v.getSize (); i++) result.set (signum (v.get (i)), i); 00281 return result; 00282 } 00283 00284 vector sign (vector v) { 00285 vector result (v); 00286 for (int i = 0; i < v.getSize (); i++) result.set (sign (v.get (i)), i); 00287 return result; 00288 } 00289 00290 vector xhypot (vector v, const nr_complex_t z) { 00291 vector result (v); 00292 for (int i = 0; i < v.getSize (); i++) result.set (xhypot (v.get(i), z), i); 00293 return result; 00294 } 00295 00296 vector xhypot (vector v, const nr_double_t d) { 00297 vector result (v); 00298 for (int i = 0; i < v.getSize (); i++) result.set (xhypot (v.get(i), d), i); 00299 return result; 00300 } 00301 00302 vector xhypot (const nr_complex_t z, vector v) { 00303 vector result (v); 00304 for (int i = 0; i < v.getSize (); i++) result.set (xhypot (z, v.get (i)), i); 00305 return result; 00306 } 00307 00308 vector xhypot (const nr_double_t d, vector v) { 00309 vector result (v); 00310 for (int i = 0; i < v.getSize (); i++) result.set (xhypot (d, v.get (i)), i); 00311 return result; 00312 } 00313 00314 vector xhypot (vector v1, vector v2) { 00315 int j, i, n, len, len1 = v1.getSize (), len2 = v2.getSize (); 00316 if (len1 >= len2) { 00317 assert (len1 % len2 == 0); 00318 len = len1; 00319 } else { 00320 assert (len2 % len1 == 0); 00321 len = len2; 00322 } 00323 vector res (len); 00324 for (j = i = n = 0; n < len; n++) { 00325 res (n) = xhypot (v1 (i), v2 (j)); 00326 if (++i >= len1) i = 0; if (++j >= len2) j = 0; 00327 } 00328 return res; 00329 } 00330 00331 vector sinc (vector v) { 00332 vector result (v); 00333 for (int i = 0; i < v.getSize (); i++) result.set (sinc (v.get (i)), i); 00334 return result; 00335 } 00336 00337 vector abs (vector v) { 00338 vector result (v); 00339 for (int i = 0; i < v.getSize (); i++) result.set (abs (v.get (i)), i); 00340 return result; 00341 } 00342 00343 vector norm (vector v) { 00344 vector result (v); 00345 for (int i = 0; i < v.getSize (); i++) result.set (norm (v.get (i)), i); 00346 return result; 00347 } 00348 00349 vector arg (vector v) { 00350 vector result (v); 00351 for (int i = 0; i < v.getSize (); i++) result.set (arg (v.get (i)), i); 00352 return result; 00353 } 00354 00355 vector real (vector v) { 00356 vector result (v); 00357 for (int i = 0; i < v.getSize (); i++) result.set (real (v.get (i)), i); 00358 return result; 00359 } 00360 00361 vector imag (vector v) { 00362 vector result (v); 00363 for (int i = 0; i < v.getSize (); i++) result.set (imag (v.get (i)), i); 00364 return result; 00365 } 00366 00367 vector conj (vector v) { 00368 vector result (v); 00369 for (int i = 0; i < v.getSize (); i++) result.set (conj (v.get (i)), i); 00370 return result; 00371 } 00372 00373 vector dB (vector v) { 00374 vector result (v); 00375 for (int i = 0; i < v.getSize (); i++) 00376 result.set (10.0 * std::log10 (norm (v.get (i))), i); 00377 return result; 00378 } 00379 00380 vector sqrt (vector v) { 00381 vector result (v); 00382 for (int i = 0; i < v.getSize (); i++) result.set (sqrt (v.get (i)), i); 00383 return result; 00384 } 00385 00386 vector exp (vector v) { 00387 vector result (v); 00388 for (int i = 0; i < v.getSize (); i++) result.set (exp (v.get (i)), i); 00389 return result; 00390 } 00391 00392 vector limexp (vector v) { 00393 vector result (v); 00394 for (int i = 0; i < v.getSize (); i++) result.set (limexp (v.get (i)), i); 00395 return result; 00396 } 00397 00398 vector log (vector v) { 00399 vector result (v); 00400 for (int i = 0; i < v.getSize (); i++) result.set (log (v.get (i)), i); 00401 return result; 00402 } 00403 00404 vector log10 (vector v) { 00405 vector result (v); 00406 for (int i = 0; i < v.getSize (); i++) result.set (log10 (v.get (i)), i); 00407 return result; 00408 } 00409 00410 vector log2 (vector v) { 00411 vector result (v); 00412 for (int i = 0; i < v.getSize (); i++) result.set (log2 (v.get (i)), i); 00413 return result; 00414 } 00415 00416 vector pow (vector v, const nr_complex_t z) { 00417 vector result (v); 00418 for (int i = 0; i < v.getSize (); i++) result.set (pow (v.get(i), z), i); 00419 return result; 00420 } 00421 00422 vector pow (vector v, const nr_double_t d) { 00423 vector result (v); 00424 for (int i = 0; i < v.getSize (); i++) result.set (pow (v.get(i), d), i); 00425 return result; 00426 } 00427 00428 vector pow (const nr_complex_t z, vector v) { 00429 vector result (v); 00430 for (int i = 0; i < v.getSize (); i++) result.set (pow (z, v.get (i)), i); 00431 return result; 00432 } 00433 00434 vector pow (const nr_double_t d, vector v) { 00435 vector result (v); 00436 for (int i = 0; i < v.getSize (); i++) result.set (pow (d, v.get (i)), i); 00437 return result; 00438 } 00439 00440 vector pow (vector v1, vector v2) { 00441 int j, i, n, len, len1 = v1.getSize (), len2 = v2.getSize (); 00442 if (len1 >= len2) { 00443 assert (len1 % len2 == 0); 00444 len = len1; 00445 } else { 00446 assert (len2 % len1 == 0); 00447 len = len2; 00448 } 00449 vector res (len); 00450 for (j = i = n = 0; n < len; n++) { 00451 res (n) = pow (v1 (i), v2 (j)); 00452 if (++i >= len1) i = 0; if (++j >= len2) j = 0; 00453 } 00454 return res; 00455 } 00456 00457 vector sin (vector v) { 00458 vector result (v); 00459 for (int i = 0; i < v.getSize (); i++) result.set (sin (v.get (i)), i); 00460 return result; 00461 } 00462 00463 vector asin (vector v) { 00464 vector result (v); 00465 for (int i = 0; i < v.getSize (); i++) result.set (asin (v.get (i)), i); 00466 return result; 00467 } 00468 00469 vector acos (vector v) { 00470 vector result (v); 00471 for (int i = 0; i < v.getSize (); i++) result.set (acos (v.get (i)), i); 00472 return result; 00473 } 00474 00475 vector cos (vector v) { 00476 vector result (v); 00477 for (int i = 0; i < v.getSize (); i++) result.set (cos (v.get (i)), i); 00478 return result; 00479 } 00480 00481 vector tan (vector v) { 00482 vector result (v); 00483 for (int i = 0; i < v.getSize (); i++) result.set (tan (v.get (i)), i); 00484 return result; 00485 } 00486 00487 vector atan (vector v) { 00488 vector result (v); 00489 for (int i = 0; i < v.getSize (); i++) result.set (atan (v.get (i)), i); 00490 return result; 00491 } 00492 00493 vector cot (vector v) { 00494 vector result (v); 00495 for (int i = 0; i < v.getSize (); i++) result.set (cot (v.get (i)), i); 00496 return result; 00497 } 00498 00499 vector acot (vector v) { 00500 vector result (v); 00501 for (int i = 0; i < v.getSize (); i++) result.set (acot (v.get (i)), i); 00502 return result; 00503 } 00504 00505 vector sinh (vector v) { 00506 vector result (v); 00507 for (int i = 0; i < v.getSize (); i++) result.set (sinh (v.get (i)), i); 00508 return result; 00509 } 00510 00511 vector asinh (vector v) { 00512 vector result (v); 00513 for (int i = 0; i < v.getSize (); i++) result.set (asinh (v.get (i)), i); 00514 return result; 00515 } 00516 00517 vector cosh (vector v) { 00518 vector result (v); 00519 for (int i = 0; i < v.getSize (); i++) result.set (cosh (v.get (i)), i); 00520 return result; 00521 } 00522 00523 vector sech (vector v) { 00524 vector result (v); 00525 for (int i = 0; i < v.getSize (); i++) result.set (sech (v.get (i)), i); 00526 return result; 00527 } 00528 00529 vector cosech (vector v) { 00530 vector result (v); 00531 for (int i = 0; i < v.getSize (); i++) result.set (cosech (v.get (i)), i); 00532 return result; 00533 } 00534 00535 vector acosh (vector v) { 00536 vector result (v); 00537 for (int i = 0; i < v.getSize (); i++) result.set (acosh (v.get (i)), i); 00538 return result; 00539 } 00540 00541 vector asech (vector v) { 00542 vector result (v); 00543 for (int i = 0; i < v.getSize (); i++) result.set (asech (v.get (i)), i); 00544 return result; 00545 } 00546 00547 vector tanh (vector v) { 00548 vector result (v); 00549 for (int i = 0; i < v.getSize (); i++) result.set (tanh (v.get (i)), i); 00550 return result; 00551 } 00552 00553 vector atanh (vector v) { 00554 vector result (v); 00555 for (int i = 0; i < v.getSize (); i++) result.set (atanh (v.get (i)), i); 00556 return result; 00557 } 00558 00559 vector coth (vector v) { 00560 vector result (v); 00561 for (int i = 0; i < v.getSize (); i++) result.set (coth (v.get (i)), i); 00562 return result; 00563 } 00564 00565 vector acoth (vector v) { 00566 vector result (v); 00567 for (int i = 0; i < v.getSize (); i++) result.set (acoth (v.get (i)), i); 00568 return result; 00569 } 00570 00571 // converts impedance to reflexion coefficient 00572 vector ztor (vector v, nr_complex_t zref) { 00573 vector result (v); 00574 for (int i = 0; i < v.getSize (); i++) result (i) = ztor (v (i), zref); 00575 return result; 00576 } 00577 00578 // converts admittance to reflexion coefficient 00579 vector ytor (vector v, nr_complex_t zref) { 00580 vector result (v); 00581 for (int i = 0; i < v.getSize (); i++) result (i) = ytor (v (i), zref); 00582 return result; 00583 } 00584 00585 // converts reflexion coefficient to impedance 00586 vector rtoz (vector v, nr_complex_t zref) { 00587 vector result (v); 00588 for (int i = 0; i < v.getSize (); i++) result (i) = rtoz (v (i), zref); 00589 return result; 00590 } 00591 00592 // converts reflexion coefficient to admittance 00593 vector rtoy (vector v, nr_complex_t zref) { 00594 vector result (v); 00595 for (int i = 0; i < v.getSize (); i++) result (i) = rtoy (v (i), zref); 00596 return result; 00597 } 00598 00599 // differentiates 'var' with respect to 'dep' exactly 'n' times 00600 vector diff (vector var, vector dep, int n) { 00601 int k, xi, yi, exchange = 0; 00602 vector x, y; 00603 // exchange dependent and independent variable if necessary 00604 if (var.getSize () < dep.getSize ()) { 00605 x = vector (var); 00606 y = vector (dep); 00607 exchange++; 00608 } 00609 else { 00610 x = vector (dep); 00611 y = vector (var); 00612 } 00613 assert (y.getSize () % x.getSize () == 0 && x.getSize () >= 2); 00614 vector result (y); 00615 nr_complex_t c; 00616 00617 for (k = 0; k < n; k++) { // differentiate n times 00618 for (yi = xi = 0; yi < y.getSize (); yi++, xi++) { 00619 if (xi == x.getSize ()) xi = 0; // roll through independent vector 00620 if (xi == 0) { 00621 c = (y.get (yi + 1) - y.get (yi)) / (x.get (xi + 1) - x.get (xi)); 00622 } else if (xi == x.getSize () - 1) { 00623 c = (y.get (yi) - y.get (yi - 1)) / (x.get (xi) - x.get (xi - 1)); 00624 } 00625 else { 00626 c = 00627 ((y.get (yi) - y.get (yi - 1)) / (x.get (xi) - x.get (xi - 1)) + 00628 (y.get (yi + 1) - y.get (yi)) / (x.get (xi + 1) - x.get (xi))) / 00629 2.0; 00630 } 00631 result.set (exchange ? 1.0 / c : c, yi); 00632 } 00633 y = result; 00634 } 00635 return result; 00636 } 00637 00638 vector vector::operator=(const nr_complex_t c) { 00639 for (int i = 0; i < size; i++) data[i] = c; 00640 return *this; 00641 } 00642 00643 vector vector::operator=(const nr_double_t d) { 00644 for (int i = 0; i < size; i++) data[i] = d; 00645 return *this; 00646 } 00647 00648 vector vector::operator+=(vector v) { 00649 int i, n, len = v.getSize (); 00650 assert (size % len == 0); 00651 for (i = n = 0; i < size; i++) { data[i] += v (n); if (++n >= len) n = 0; } 00652 return *this; 00653 } 00654 00655 vector vector::operator+=(const nr_complex_t c) { 00656 for (int i = 0; i < size; i++) data[i] += c; 00657 return *this; 00658 } 00659 00660 vector vector::operator+=(const nr_double_t d) { 00661 for (int i = 0; i < size; i++) data[i] += d; 00662 return *this; 00663 } 00664 00665 vector operator+(vector v1, vector v2) { 00666 int len1 = v1.getSize (), len2 = v2.getSize (); 00667 vector result; 00668 if (len1 >= len2) { 00669 result = v1; 00670 result += v2; 00671 } else { 00672 result = v2; 00673 result += v1; 00674 } 00675 return result; 00676 } 00677 00678 vector operator+(vector v, const nr_complex_t c) { 00679 vector result (v); 00680 result += c; 00681 return result; 00682 } 00683 00684 vector operator+(const nr_complex_t c, vector v) { 00685 return v + c; 00686 } 00687 00688 vector operator+(vector v, const nr_double_t d) { 00689 vector result (v); 00690 result += d; 00691 return result; 00692 } 00693 00694 vector operator+(const nr_double_t d, vector v) { 00695 return v + d; 00696 } 00697 00698 vector vector::operator-() { 00699 vector result (size); 00700 for (int i = 0; i < size; i++) result (i) = -data[i]; 00701 return result; 00702 } 00703 00704 vector vector::operator-=(vector v) { 00705 int i, n, len = v.getSize (); 00706 assert (size % len == 0); 00707 for (i = n = 0; i < size; i++) { data[i] -= v (n); if (++n >= len) n = 0; } 00708 return *this; 00709 } 00710 00711 vector vector::operator-=(const nr_complex_t c) { 00712 for (int i = 0; i < size; i++) data[i] -= c; 00713 return *this; 00714 } 00715 00716 vector vector::operator-=(const nr_double_t d) { 00717 for (int i = 0; i < size; i++) data[i] -= d; 00718 return *this; 00719 } 00720 00721 vector operator-(vector v1, vector v2) { 00722 int len1 = v1.getSize (), len2 = v2.getSize (); 00723 vector result; 00724 if (len1 >= len2) { 00725 result = v1; 00726 result -= v2; 00727 } else { 00728 result = -v2; 00729 result += v1; 00730 } 00731 return result; 00732 } 00733 00734 vector operator-(vector v, const nr_complex_t c) { 00735 vector result (v); 00736 result -= c; 00737 return result; 00738 } 00739 00740 vector operator-(vector v, const nr_double_t d) { 00741 vector result (v); 00742 result -= d; 00743 return result; 00744 } 00745 00746 vector operator-(const nr_complex_t c, vector v) { 00747 vector result (-v); 00748 result += c; 00749 return result; 00750 } 00751 00752 vector operator-(const nr_double_t d, vector v) { 00753 vector result (-v); 00754 result += d; 00755 return result; 00756 } 00757 00758 vector vector::operator*=(vector v) { 00759 int i, n, len = v.getSize (); 00760 assert (size % len == 0); 00761 for (i = n = 0; i < size; i++) { data[i] *= v (n); if (++n >= len) n = 0; } 00762 return *this; 00763 } 00764 00765 vector vector::operator*=(const nr_complex_t c) { 00766 for (int i = 0; i < size; i++) data[i] *= c; 00767 return *this; 00768 } 00769 00770 vector vector::operator*=(const nr_double_t d) { 00771 for (int i = 0; i < size; i++) data[i] *= d; 00772 return *this; 00773 } 00774 00775 vector operator*(vector v1, vector v2) { 00776 int len1 = v1.getSize (), len2 = v2.getSize (); 00777 vector result; 00778 if (len1 >= len2) { 00779 result = v1; 00780 result *= v2; 00781 } else { 00782 result = v2; 00783 result *= v1; 00784 } 00785 return result; 00786 } 00787 00788 vector operator*(vector v, const nr_complex_t c) { 00789 vector result (v); 00790 result *= c; 00791 return result; 00792 } 00793 00794 vector operator*(vector v, const nr_double_t d) { 00795 vector result (v); 00796 result *= d; 00797 return result; 00798 } 00799 00800 vector operator*(const nr_complex_t c, vector v) { 00801 return v * c; 00802 } 00803 00804 vector operator*(const nr_double_t d, vector v) { 00805 return v * d; 00806 } 00807 00808 vector vector::operator/=(vector v) { 00809 int i, n, len = v.getSize (); 00810 assert (size % len == 0); 00811 for (i = n = 0; i < size; i++) { data[i] /= v (n); if (++n >= len) n = 0; } 00812 return *this; 00813 } 00814 00815 vector vector::operator/=(const nr_complex_t c) { 00816 for (int i = 0; i < size; i++) data[i] /= c; 00817 return *this; 00818 } 00819 00820 vector vector::operator/=(const nr_double_t d) { 00821 for (int i = 0; i < size; i++) data[i] /= d; 00822 return *this; 00823 } 00824 00825 vector operator/(vector v1, vector v2) { 00826 int len1 = v1.getSize (), len2 = v2.getSize (); 00827 vector result; 00828 if (len1 >= len2) { 00829 assert (len1 % len2 == 0); 00830 result = v1; 00831 result /= v2; 00832 } else { 00833 assert (len2 % len1 == 0); 00834 result = 1 / v2; 00835 result *= v1; 00836 } 00837 return result; 00838 } 00839 00840 vector operator/(vector v, const nr_complex_t c) { 00841 vector result (v); 00842 result /= c; 00843 return result; 00844 } 00845 00846 vector operator/(vector v, const nr_double_t d) { 00847 vector result (v); 00848 result /= d; 00849 return result; 00850 } 00851 00852 vector operator/(const nr_complex_t c, vector v) { 00853 vector result (v); 00854 result = c; 00855 result /= v; 00856 return result; 00857 } 00858 00859 vector operator/(const nr_double_t d, vector v) { 00860 vector result (v); 00861 result = d; 00862 result /= v; 00863 return result; 00864 } 00865 00866 vector operator%(vector v, const nr_complex_t z) { 00867 int len = v.getSize (); 00868 vector result (len); 00869 for (int i = 0; i < len; i++) result (i) = v (i) % z; 00870 return result; 00871 } 00872 00873 vector operator%(vector v, const nr_double_t d) { 00874 int len = v.getSize (); 00875 vector result (len); 00876 for (int i = 0; i < len; i++) result (i) = v (i) % d; 00877 return result; 00878 } 00879 00880 vector operator%(const nr_complex_t z, vector v) { 00881 int len = v.getSize (); 00882 vector result (len); 00883 for (int i = 0; i < len; i++) result (i) = z % v (i); 00884 return result; 00885 } 00886 00887 vector operator%(const nr_double_t d, vector v) { 00888 int len = v.getSize (); 00889 vector result (len); 00890 for (int i = 0; i < len; i++) result (i) = d % v (i); 00891 return result; 00892 } 00893 00894 vector operator%(vector v1, vector v2) { 00895 int j, i, n, len, len1 = v1.getSize (), len2 = v2.getSize (); 00896 if (len1 >= len2) { 00897 assert (len1 % len2 == 0); 00898 len = len1; 00899 } else { 00900 assert (len2 % len1 == 0); 00901 len = len2; 00902 } 00903 vector res (len); 00904 for (j = i = n = 0; n < len; n++) { 00905 res (n) = v1 (i) % v2 (j); 00906 if (++i >= len1) i = 0; if (++j >= len2) j = 0; 00907 } 00908 return res; 00909 } 00910 00911 /* This function reverses the order of the data list. */ 00912 void vector::reverse (void) { 00913 nr_complex_t * buffer = (nr_complex_t *) 00914 malloc (sizeof (nr_complex_t) * size); 00915 for (int i = 0; i < size; i++) buffer[i] = data[size - 1 - i]; 00916 free (data); 00917 data = buffer; 00918 capacity = size; 00919 } 00920 00921 // Sets the origin (the analysis) of the vector. 00922 void vector::setOrigin (const char * n) { 00923 free (origin); 00924 origin = n ? strdup (n) : NULL; 00925 } 00926 00927 // Returns the origin (the analysis) of the vector. 00928 char * vector::getOrigin (void) { 00929 return origin; 00930 } 00931 00932 /* The function returns the number of entries with the given value 00933 deviating no more than the given epsilon. */ 00934 int vector::contains (nr_complex_t val, nr_double_t eps) { 00935 int count = 0; 00936 for (int i = 0; i < size; i++) { 00937 if (abs (data[i] - val) <= eps) count++; 00938 } 00939 return count; 00940 } 00941 00942 // Sorts the vector either in ascending or descending order. 00943 void vector::sort (bool ascending) { 00944 nr_complex_t t; 00945 for (int i = 0; i < size; i++) { 00946 for (int n = 0; n < size - 1; n++) { 00947 if (ascending ? data[n] > data[n+1] : data[n] < data[n+1]) { 00948 t = data[n]; 00949 data[n] = data[n+1]; 00950 data[n+1] = t; 00951 } 00952 } 00953 } 00954 } 00955 00956 /* The function creates a linear stepped vector of values starting at 00957 the given start value, ending with the given stop value and 00958 containing points elements. */ 00959 vector linspace (nr_double_t start, nr_double_t stop, int points) { 00960 vector result (points); 00961 nr_double_t val, step = (stop - start) / (points - 1); 00962 for (int i = 0; i < points; i++) { 00963 val = start + (i * step); 00964 if (i != 0 && fabs (val) < fabs (step) / 4 && fabs (val) < std::numeric_limits<nr_double_t>::epsilon()) 00965 val = 0.0; 00966 result.set (val, i); 00967 } 00968 return result; 00969 } 00970 00971 /* The function creates a logarithmic stepped vector of values 00972 starting at the given start value, ending with the given stop value 00973 and containing points elements. */ 00974 vector logspace (nr_double_t start, nr_double_t stop, int points) { 00975 assert (start * stop > 0); 00976 vector result (points); 00977 nr_double_t step, first, last, d; 00978 00979 // ensure the last value being larger than the first 00980 if (fabs (start) > fabs (stop)) { 00981 first = fabs (stop); 00982 last = fabs (start); 00983 } 00984 else { 00985 first = fabs (start); 00986 last = fabs (stop); 00987 } 00988 // check direction and sign of values 00989 d = fabs (start) > fabs (stop) ? -1 : 1; 00990 // compute logarithmic step size 00991 step = (::log (last) - ::log (first)) / (points - 1); 00992 for (int i = 0, j = points - 1; i < points; i++, j--) { 00993 if (d > 0) 00994 result.set (start * ::exp (step * i), i); 00995 else 00996 result.set (stop * ::exp (step * i), j); 00997 } 00998 return result; 00999 } 01000 01001 vector cumsum (vector v) { 01002 vector result (v); 01003 nr_complex_t val (0.0); 01004 for (int i = 0; i < v.getSize (); i++) { 01005 val += v.get (i); 01006 result.set (val, i); 01007 } 01008 return result; 01009 } 01010 01011 vector cumavg (vector v) { 01012 vector result (v); 01013 nr_complex_t val (0.0); 01014 for (int i = 0; i < v.getSize (); i++) { 01015 val = (val * (nr_double_t) i + v.get (i)) / (i + 1.0); 01016 result.set (val, i); 01017 } 01018 return result; 01019 } 01020 01021 vector cumprod (vector v) { 01022 vector result (v); 01023 nr_complex_t val (1.0); 01024 for (int i = 0; i < v.getSize (); i++) { 01025 val *= v.get (i); 01026 result.set (val, i); 01027 } 01028 return result; 01029 } 01030 01031 vector ceil (vector v) { 01032 vector result (v); 01033 for (int i = 0; i < v.getSize (); i++) result.set (ceil (v.get (i)), i); 01034 return result; 01035 } 01036 01037 vector fix (vector v) { 01038 vector result (v); 01039 for (int i = 0; i < v.getSize (); i++) result.set (fix (v.get (i)), i); 01040 return result; 01041 } 01042 01043 vector floor (vector v) { 01044 vector result (v); 01045 for (int i = 0; i < v.getSize (); i++) result.set (floor (v.get (i)), i); 01046 return result; 01047 } 01048 01049 vector round (vector v) { 01050 vector result (v); 01051 for (int i = 0; i < v.getSize (); i++) result.set (round (v.get (i)), i); 01052 return result; 01053 } 01054 01055 vector sqr (vector v) { 01056 vector result (v); 01057 for (int i = 0; i < v.getSize (); i++) result.set (sqr (v.get (i)), i); 01058 return result; 01059 } 01060 01061 vector step (vector v) { 01062 vector result (v); 01063 for (int i = 0; i < v.getSize (); i++) result.set (step (v.get (i)), i); 01064 return result; 01065 } 01066 01067 static nr_double_t integrate_n (vector v) { /* using trapezoidal rule */ 01068 nr_double_t result = 0.0; 01069 for (int i = 1; i < v.getSize () - 1; i++) result += norm (v.get (i)); 01070 result += 0.5 * norm (v.get (0)); 01071 result += 0.5 * norm (v.get (v.getSize () - 1)); 01072 return result; 01073 } 01074 01075 nr_double_t vector::rms (void) { 01076 nr_double_t result = std::sqrt (integrate_n (*this) / getSize ()); 01077 return result; 01078 } 01079 01080 nr_double_t vector::variance (void) { 01081 nr_double_t result = 0.0; 01082 nr_complex_t average = avg (*this); 01083 for (int i = 0; i < getSize (); i++) result += norm (get (i) - average); 01084 if (getSize () > 1) 01085 return result / (getSize () - 1); 01086 return 0.0; 01087 } 01088 01089 nr_double_t vector::stddev (void) { 01090 return std::sqrt (variance ()); 01091 } 01092 01093 vector erf (vector v) { 01094 vector result (v); 01095 for (int i = 0; i < v.getSize (); i++) result.set (erf (v.get (i)), i); 01096 return result; 01097 } 01098 01099 vector erfc (vector v) { 01100 vector result (v); 01101 for (int i = 0; i < v.getSize (); i++) result.set (erfc (v.get (i)), i); 01102 return result; 01103 } 01104 01105 vector erfinv (vector v) { 01106 vector result (v); 01107 for (int i = 0; i < v.getSize (); i++) result.set (erfinv (v.get (i)), i); 01108 return result; 01109 } 01110 01111 vector erfcinv (vector v) { 01112 vector result (v); 01113 for (int i = 0; i < v.getSize (); i++) result.set (erfcinv (v.get (i)), i); 01114 return result; 01115 } 01116 01117 vector rad2deg (vector v) { 01118 vector result (v); 01119 for (int i = 0; i < v.getSize (); i++) result.set (rad2deg (v.get (i)), i); 01120 return result; 01121 } 01122 01123 vector deg2rad (vector v) { 01124 vector result (v); 01125 for (int i = 0; i < v.getSize (); i++) result.set (deg2rad (v.get (i)), i); 01126 return result; 01127 } 01128 01129 vector i0 (vector v) { 01130 vector result (v); 01131 for (int i = 0; i < v.getSize (); i++) result.set (i0 (v.get (i)), i); 01132 return result; 01133 } 01134 01135 vector jn (const int n, vector v) { 01136 vector result (v); 01137 for (int i = 0; i < v.getSize (); i++) result.set (jn (n, v.get (i)), i); 01138 return result; 01139 } 01140 01141 vector yn (const int n, vector v) { 01142 vector result (v); 01143 for (int i = 0; i < v.getSize (); i++) result.set (yn (n, v.get (i)), i); 01144 return result; 01145 } 01146 01147 vector polar (const nr_complex_t a, vector v) { 01148 vector result (v); 01149 for (int i = 0; i < v.getSize (); i++) result.set (qucs::polar (a, v.get (i)), i); 01150 return result; 01151 } 01152 01153 vector polar (vector v, const nr_complex_t p) { 01154 vector result (v); 01155 for (int i = 0; i < v.getSize (); i++) result.set (qucs::polar (v.get (i), p), i); 01156 return result; 01157 } 01158 01159 vector polar (vector a, vector p) { 01160 int j, i, n, len, len1 = a.getSize (), len2 = p.getSize (); 01161 if (len1 >= len2) { 01162 assert (len1 % len2 == 0); 01163 len = len1; 01164 } else { 01165 assert (len2 % len1 == 0); 01166 len = len2; 01167 } 01168 vector res (len); 01169 for (j = i = n = 0; n < len; n++) { 01170 res (n) = qucs::polar (a (i), p (j)); 01171 if (++i >= len1) i = 0; if (++j >= len2) j = 0; 01172 } 01173 return res; 01174 } 01175 01176 vector atan2 (const nr_double_t y, vector v) { 01177 vector result (v); 01178 for (int i = 0; i < v.getSize (); i++) 01179 result.set (atan2 (y, v.get (i)), i); 01180 return result; 01181 } 01182 01183 vector atan2 (vector v, const nr_double_t x) { 01184 vector result (v); 01185 for (int i = 0; i < v.getSize (); i++) 01186 result.set (atan2 (v.get (i), x) , i); 01187 return result; 01188 } 01189 01190 vector atan2 (vector y, vector x) { 01191 int j, i, n, len, len1 = y.getSize (), len2 = x.getSize (); 01192 if (len1 >= len2) { 01193 assert (len1 % len2 == 0); 01194 len = len1; 01195 } else { 01196 assert (len2 % len1 == 0); 01197 len = len2; 01198 } 01199 vector res (len); 01200 for (j = i = n = 0; n < len; n++) { 01201 res (n) = atan2 (y (i), x (j)); 01202 if (++i >= len1) i = 0; if (++j >= len2) j = 0; 01203 } 01204 return res; 01205 } 01206 01207 vector w2dbm (vector v) { 01208 vector result (v); 01209 for (int i = 0; i < v.getSize (); i++) 01210 result.set (10.0 * log10 (v.get (i) / 0.001), i); 01211 return result; 01212 } 01213 01214 vector dbm2w (vector v) { 01215 vector result (v); 01216 for (int i = 0; i < v.getSize (); i++) 01217 result.set (0.001 * pow (10.0 , v.get (i) / 10.0), i); 01218 return result; 01219 } 01220 01221 nr_double_t integrate (vector v, const nr_double_t h) { 01222 nr_double_t s = real (v.get (0) ) / 2; 01223 for (int i = 1; i < v.getSize () - 1; i++) 01224 s += real (v.get (i)); 01225 return (s + real (v.get (v.getSize () - 1) ) / 2) * h; 01226 } 01227 01228 nr_complex_t integrate (vector v, const nr_complex_t h) { 01229 nr_complex_t s; 01230 s = v.get (0) / 2.0; 01231 for (int i = 1; i < v.getSize () - 1; i++) 01232 s += v.get (i); 01233 return (s + v.get (v.getSize () - 1) / 2.0) * h; 01234 } 01235 01236 vector dbm (vector v, const nr_complex_t z) { 01237 vector result (v); 01238 for (int i = 0; i < v.getSize (); i++) 01239 result.set (10.0 * log10 (norm (v.get (i)) / conj (z) / 0.001), i); 01240 return result; 01241 } 01242 01243 vector runavg (const nr_complex_t x, const int n) { 01244 vector result (n); 01245 for (int i = 0; i < n; i++) result.set (x, i); 01246 return result; 01247 } 01248 01249 vector runavg (vector v, const int n) { 01250 nr_complex_t s (0.0), y; 01251 int len = v.getSize () - n + 1, i; 01252 vector result (len); 01253 for (i = 0; i < n; i++) s += v.get (i); 01254 y = s / (nr_double_t) n; // first running average value 01255 result.set (y, 0); 01256 for (i = 0; i < len - 1; i++) { 01257 y += (v.get (i + n) - v.get (i)) / (nr_double_t) n; 01258 result.set (y, i + 1); 01259 } 01260 return result; 01261 } 01262 01263 #ifdef DEBUG 01264 // Debug function: Prints the vector object. 01265 void vector::print (void) { 01266 for (int r = 0; r < size; r++) { 01267 fprintf (stderr, "%+.2e%+.2ei\n", (double) real (get (r)), 01268 (double) imag (get (r))); 01269 } 01270 } 01271 #endif /* DEBUG */ 01272 01273 } // namespace qucs