Branch data Line data Source code
1 : : /*
2 : : * vector.cpp - vector class implementation
3 : : *
4 : : * Copyright (C) 2003, 2004, 2005, 2006, 2007 Stefan Jahn <stefan@lkcc.org>
5 : : * Copyright (C) 2006, 2007 Gunther Kraut <gn.kraut@t-online.de>
6 : : *
7 : : * This is free software; you can redistribute it and/or modify
8 : : * it under the terms of the GNU General Public License as published by
9 : : * the Free Software Foundation; either version 2, or (at your option)
10 : : * any later version.
11 : : *
12 : : * This software is distributed in the hope that it will be useful,
13 : : * but WITHOUT ANY WARRANTY; without even the implied warranty of
14 : : * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
15 : : * GNU General Public License for more details.
16 : : *
17 : : * You should have received a copy of the GNU General Public License
18 : : * along with this package; see the file COPYING. If not, write to
19 : : * the Free Software Foundation, Inc., 51 Franklin Street - Fifth Floor,
20 : : * Boston, MA 02110-1301, USA.
21 : : *
22 : : * $Id$
23 : : *
24 : : */
25 : :
26 : : #if HAVE_CONFIG_H
27 : : # include <config.h>
28 : : #endif
29 : :
30 : : #include <limits>
31 : :
32 : : #include <stdio.h>
33 : : #include <stdlib.h>
34 : : #include <string.h>
35 : : #include <cmath>
36 : : #include <float.h>
37 : : #include <assert.h>
38 : :
39 : : #include "real.h"
40 : : #include "complex.h"
41 : : #include "object.h"
42 : : #include "logging.h"
43 : : #include "strlist.h"
44 : : #include "vector.h"
45 : : #include "consts.h"
46 : :
47 : : namespace qucs {
48 : :
49 : : // Constructor creates an unnamed instance of the vector class.
50 : 146 : vector::vector () : object () {
51 : 146 : capacity = size = 0;
52 : 146 : data = NULL;
53 : 146 : dependencies = NULL;
54 : 146 : origin = NULL;
55 : 146 : requested = 0;
56 : 146 : }
57 : :
58 : : /* Constructor creates an unnamed instance of the vector class with a
59 : : given initial size. */
60 : 174 : vector::vector (int s) : object () {
61 [ - + ][ # # ]: 174 : assert (s >= 0);
62 : 174 : capacity = size = s;
63 : : data = s > 0 ? (nr_complex_t *)
64 [ + - ][ # # ]: 174 : calloc (capacity, sizeof (nr_complex_t)) : NULL;
65 : 174 : dependencies = NULL;
66 : 174 : origin = NULL;
67 : 174 : requested = 0;
68 : 174 : }
69 : :
70 : : /* Constructor creates an unnamed instance of the vector class with a
71 : : given initial size and content. */
72 : 1070 : vector::vector (int s, nr_complex_t val) : object () {
73 [ - + ][ # # ]: 1070 : assert (s >= 0);
74 : 1070 : capacity = size = s;
75 : : data = s > 0 ? (nr_complex_t *)
76 [ + - ][ # # ]: 1070 : calloc (capacity, sizeof (nr_complex_t)) : NULL;
77 [ + + ][ # # ]: 4570 : for (int i = 0; i < s; i++) data[i] = val;
78 : 1070 : dependencies = NULL;
79 : 1070 : origin = NULL;
80 : 1070 : requested = 0;
81 : 1070 : }
82 : :
83 : : // Constructor creates an named instance of the vector class.
84 : 1126 : vector::vector (const char * n) : object (n) {
85 : 1126 : capacity = size = 0;
86 : 1126 : data = NULL;
87 : 1126 : dependencies = NULL;
88 : 1126 : origin = NULL;
89 : 1126 : requested = 0;
90 : 1126 : }
91 : :
92 : : /* This constructor creates a named instance of the vector class with
93 : : a given initial size. */
94 : 0 : vector::vector (const char * n, int s) : object (n) {
95 [ # # ][ # # ]: 0 : assert (s >= 0);
96 : 0 : capacity = size = s;
97 : : data = s > 0 ? (nr_complex_t *)
98 [ # # ][ # # ]: 0 : calloc (capacity, sizeof (nr_complex_t)) : NULL;
99 : 0 : dependencies = NULL;
100 : 0 : origin = NULL;
101 : 0 : requested = 0;
102 : 0 : }
103 : :
104 : : /* The copy constructor creates a new instance based on the given
105 : : vector object. */
106 : 9223 : vector::vector (const vector & v) : object (v) {
107 : 9223 : size = v.size;
108 : 9223 : capacity = v.capacity;
109 : 9223 : data = (nr_complex_t *) malloc (sizeof (nr_complex_t) * capacity);
110 : 9223 : memcpy (data, v.data, sizeof (nr_complex_t) * size);
111 [ + + ][ + - ]: 9223 : dependencies = v.dependencies ? new strlist (*v.dependencies) : NULL;
[ + - ][ # # ]
[ # # ][ # # ]
112 [ + + ][ + - ]: 9223 : origin = v.origin ? strdup (v.origin) : NULL;
[ # # ][ # # ]
113 : 9223 : requested = v.requested;
114 : 9223 : }
115 : :
116 : : /* The assignment copy constructor creates a new instance based on the
117 : : given vector object. It copies the data only and leaves any other
118 : : properties untouched. */
119 : 42 : const vector& vector::operator=(const vector & v) {
120 [ + - ]: 42 : if (&v != this) {
121 : 42 : size = v.size;
122 : 42 : capacity = v.capacity;
123 [ + + ]: 42 : if (data) { free (data); data = NULL; }
124 [ + - ]: 42 : if (capacity > 0) {
125 : 42 : data = (nr_complex_t *) malloc (sizeof (nr_complex_t) * capacity);
126 [ + - ]: 42 : if (size > 0) memcpy (data, v.data, sizeof (nr_complex_t) * size);
127 : : }
128 : : }
129 : 42 : return *this;
130 : : }
131 : :
132 : : // Destructor deletes a vector object.
133 : 11733 : vector::~vector () {
134 [ + - ][ # # ]: 11733 : if (data) free (data);
135 [ + + ][ + - ]: 11733 : if (dependencies) delete dependencies;
[ + - ][ # # ]
[ # # ][ # # ]
136 [ + + ][ # # ]: 11733 : if (origin) free (origin);
137 [ - + ][ # # ]: 13229 : }
138 : :
139 : : // Returns data dependencies.
140 : 90296 : strlist * vector::getDependencies (void) {
141 : 90296 : return dependencies;
142 : : }
143 : :
144 : : // Sets the data dependencies.
145 : 927 : void vector::setDependencies (strlist * s) {
146 [ + + ][ + - ]: 927 : if (dependencies) delete dependencies;
147 : 927 : dependencies = s;
148 : 927 : }
149 : :
150 : : /* The function appends a new complex data item to the end of the
151 : : vector and ensures that the vector can hold the increasing number
152 : : of data items. */
153 : 602265 : void vector::add (nr_complex_t c) {
154 [ + + ]: 602265 : if (data == NULL) {
155 : 1235 : size = 0; capacity = 64;
156 : 1235 : data = (nr_complex_t *) malloc (sizeof (nr_complex_t) * capacity);
157 : : }
158 [ + + ]: 601030 : else if (size >= capacity) {
159 : 2236 : capacity *= 2;
160 : 2236 : data = (nr_complex_t *) realloc (data, sizeof (nr_complex_t) * capacity);
161 : : }
162 : 602265 : data[size++] = c;
163 : 602265 : }
164 : :
165 : : /* This function appends the given vector to the vector. */
166 : 0 : void vector::add (vector * v) {
167 [ # # ]: 0 : if (v != NULL) {
168 [ # # ]: 0 : if (data == NULL) {
169 : 0 : size = 0; capacity = v->getSize ();
170 : 0 : data = (nr_complex_t *) malloc (sizeof (nr_complex_t) * capacity);
171 : : }
172 [ # # ]: 0 : else if (size + v->getSize () > capacity) {
173 : 0 : capacity += v->getSize ();
174 : 0 : data = (nr_complex_t *) realloc (data, sizeof (nr_complex_t) * capacity);
175 : : }
176 [ # # ]: 0 : for (int i = 0; i < v->getSize (); i++) data[size++] = v->get (i);
177 : : }
178 : 0 : }
179 : :
180 : : // Returns the complex data item at the given position.
181 : 830476 : nr_complex_t vector::get (int i) {
182 : 830476 : return data[i];
183 : : }
184 : :
185 : 124038 : void vector::set (nr_double_t d, int i) {
186 : 124038 : data[i] = nr_complex_t (d);
187 : 124038 : }
188 : :
189 : 10862 : void vector::set (const nr_complex_t z, int i) {
190 : 10862 : data[i] = nr_complex_t (z);
191 : 10862 : }
192 : :
193 : : // The function returns the current size of the vector.
194 : 751446 : int vector::getSize (void) const {
195 : 751446 : return size;
196 : : }
197 : :
198 : 0 : int vector::checkSizes (vector v1, vector v2) {
199 [ # # ]: 0 : if (v1.getSize () != v2.getSize ()) {
200 : : logprint (LOG_ERROR, "vector '%s' and '%s' have different sizes\n",
201 : 0 : v1.getName (), v2.getName ());
202 : 0 : return 0;
203 : : }
204 : 0 : return 1;
205 : : }
206 : :
207 : : // searches the maximum value of the vector elements.
208 : : // complex numbers in the 1. and 4. quadrant are counted as "abs(c)".
209 : : // complex numbers in the 2. and 3. quadrant are counted as "-abs(c)".
210 : 0 : nr_double_t vector::maximum (void) {
211 : 0 : nr_complex_t c;
212 : 0 : nr_double_t d, max_D = -std::numeric_limits<nr_double_t>::max();
213 [ # # ]: 0 : for (int i = 0; i < getSize (); i++) {
214 : 0 : c = data[i];
215 [ # # ]: 0 : d = fabs (arg (c)) < M_PI_2 ? abs (c) : -abs (c);
216 [ # # ]: 0 : if (d > max_D) max_D = d;
217 : : }
218 : 0 : return max_D;
219 : : }
220 : :
221 : : // searches the minimum value of the vector elements.
222 : : // complex numbers in the 1. and 4. quadrant are counted as "abs(c)".
223 : : // complex numbers in the 2. and 3. quadrant are counted as "-abs(c)".
224 : 2 : nr_double_t vector::minimum (void) {
225 : 2 : nr_complex_t c;
226 : 2 : nr_double_t d, min_D = +std::numeric_limits<nr_double_t>::max();
227 [ + + ]: 4 : for (int i = 0; i < getSize (); i++) {
228 : 2 : c = data[i];
229 [ + - ]: 2 : d = fabs (arg (c)) < M_PI_2 ? abs (c) : -abs (c);
230 [ + - ]: 2 : if (d < min_D) min_D = d;
231 : : }
232 : 2 : return min_D;
233 : : }
234 : :
235 : : /* Unwraps a phase vector in radians. Adds +/- 2*Pi if consecutive
236 : : values jump about |Pi|. */
237 : 2 : vector unwrap (vector v, nr_double_t tol, nr_double_t step) {
238 [ + - ]: 2 : vector result (v.getSize ());
239 : 2 : nr_double_t add = 0;
240 : 2 : result (0) = v (0);
241 [ + + ]: 916 : for (int i = 1; i < v.getSize (); i++) {
242 [ + - ]: 914 : nr_double_t diff = real (v (i) - v (i-1));
243 [ + + ]: 914 : if (diff > +tol) {
244 : 2 : add -= step;
245 [ - + ]: 912 : } else if (diff < -tol) {
246 : 0 : add += step;
247 : : }
248 : 914 : result (i) = v (i) + add;
249 : : }
250 : 2 : return result;
251 : : }
252 : :
253 : 2 : nr_complex_t sum (vector v) {
254 : 2 : nr_complex_t result (0.0);
255 [ + + ]: 4 : for (int i = 0; i < v.getSize (); i++) result += v.get (i);
256 : 2 : return result;
257 : : }
258 : :
259 : 0 : nr_complex_t prod (vector v) {
260 : 0 : nr_complex_t result (1.0);
261 [ # # ]: 0 : for (int i = 0; i < v.getSize (); i++) result *= v.get (i);
262 : 0 : return result;
263 : : }
264 : :
265 : 0 : nr_complex_t avg (vector v) {
266 : 0 : nr_complex_t result (0.0);
267 [ # # ]: 0 : for (int i = 0; i < v.getSize (); i++) result += v.get (i);
268 : 0 : return result / (nr_double_t) v.getSize ();
269 : : }
270 : :
271 : 0 : vector signum (vector v) {
272 : 0 : vector result (v);
273 [ # # ][ # # ]: 0 : for (int i = 0; i < v.getSize (); i++) result.set (signum (v.get (i)), i);
274 : 0 : return result;
275 : : }
276 : :
277 : 0 : vector sign (vector v) {
278 : 0 : vector result (v);
279 [ # # ][ # # ]: 0 : for (int i = 0; i < v.getSize (); i++) result.set (sign (v.get (i)), i);
280 : 0 : return result;
281 : : }
282 : :
283 : 0 : vector xhypot (vector v, const nr_complex_t z) {
284 : 0 : vector result (v);
285 [ # # ][ # # ]: 0 : for (int i = 0; i < v.getSize (); i++) result.set (xhypot (v.get(i), z), i);
286 : 0 : return result;
287 : : }
288 : :
289 : 0 : vector xhypot (vector v, const nr_double_t d) {
290 : 0 : vector result (v);
291 [ # # ][ # # ]: 0 : for (int i = 0; i < v.getSize (); i++) result.set (xhypot (v.get(i), d), i);
292 : 0 : return result;
293 : : }
294 : :
295 : 0 : vector xhypot (const nr_complex_t z, vector v) {
296 : 0 : vector result (v);
297 [ # # ][ # # ]: 0 : for (int i = 0; i < v.getSize (); i++) result.set (xhypot (z, v.get (i)), i);
298 : 0 : return result;
299 : : }
300 : :
301 : 0 : vector xhypot (const nr_double_t d, vector v) {
302 : 0 : vector result (v);
303 [ # # ][ # # ]: 0 : for (int i = 0; i < v.getSize (); i++) result.set (xhypot (d, v.get (i)), i);
304 : 0 : return result;
305 : : }
306 : :
307 : 0 : vector xhypot (vector v1, vector v2) {
308 : 0 : int j, i, n, len, len1 = v1.getSize (), len2 = v2.getSize ();
309 [ # # ]: 0 : if (len1 >= len2) {
310 [ # # ]: 0 : assert (len1 % len2 == 0);
311 : 0 : len = len1;
312 : : } else {
313 [ # # ]: 0 : assert (len2 % len1 == 0);
314 : 0 : len = len2;
315 : : }
316 : 0 : vector res (len);
317 [ # # ]: 0 : for (j = i = n = 0; n < len; n++) {
318 [ # # ]: 0 : res (n) = xhypot (v1 (i), v2 (j));
319 [ # # ][ # # ]: 0 : if (++i >= len1) i = 0; if (++j >= len2) j = 0;
320 : : }
321 : 0 : return res;
322 : : }
323 : :
324 : 0 : vector sinc (vector v) {
325 : 0 : vector result (v);
326 [ # # ][ # # ]: 0 : for (int i = 0; i < v.getSize (); i++) result.set (sinc (v.get (i)), i);
327 : 0 : return result;
328 : : }
329 : :
330 : 26 : vector abs (vector v) {
331 : 26 : vector result (v);
332 [ + + ]: 38068 : for (int i = 0; i < v.getSize (); i++) result.set (abs (v.get (i)), i);
333 : 26 : return result;
334 : : }
335 : :
336 : 0 : vector norm (vector v) {
337 : 0 : vector result (v);
338 [ # # ][ # # ]: 0 : for (int i = 0; i < v.getSize (); i++) result.set (norm (v.get (i)), i);
339 : 0 : return result;
340 : : }
341 : :
342 : 6 : vector arg (vector v) {
343 : 6 : vector result (v);
344 [ + + ]: 2423 : for (int i = 0; i < v.getSize (); i++) result.set (arg (v.get (i)), i);
345 : 6 : return result;
346 : : }
347 : :
348 : 1071 : vector real (vector v) {
349 : 1071 : vector result (v);
350 [ + + ]: 4721 : for (int i = 0; i < v.getSize (); i++) result.set (real (v.get (i)), i);
351 : 1071 : return result;
352 : : }
353 : :
354 : 1 : vector imag (vector v) {
355 : 1 : vector result (v);
356 [ + + ]: 151 : for (int i = 0; i < v.getSize (); i++) result.set (imag (v.get (i)), i);
357 : 1 : return result;
358 : : }
359 : :
360 : 0 : vector conj (vector v) {
361 : 0 : vector result (v);
362 [ # # ]: 0 : for (int i = 0; i < v.getSize (); i++) result.set (conj (v.get (i)), i);
363 : 0 : return result;
364 : : }
365 : :
366 : 32 : vector dB (vector v) {
367 : 32 : vector result (v);
368 [ + + ]: 12324 : for (int i = 0; i < v.getSize (); i++)
369 [ + - ]: 12292 : result.set (10.0 * std::log10 (norm (v.get (i))), i);
370 : 32 : return result;
371 : : }
372 : :
373 : 1070 : vector sqrt (vector v) {
374 : 1070 : vector result (v);
375 [ + - ][ + + ]: 4570 : for (int i = 0; i < v.getSize (); i++) result.set (sqrt (v.get (i)), i);
376 : 1070 : return result;
377 : : }
378 : :
379 : 0 : vector exp (vector v) {
380 : 0 : vector result (v);
381 [ # # ][ # # ]: 0 : for (int i = 0; i < v.getSize (); i++) result.set (exp (v.get (i)), i);
382 : 0 : return result;
383 : : }
384 : :
385 : 0 : vector limexp (vector v) {
386 : 0 : vector result (v);
387 [ # # ][ # # ]: 0 : for (int i = 0; i < v.getSize (); i++) result.set (limexp (v.get (i)), i);
388 : 0 : return result;
389 : : }
390 : :
391 : 0 : vector log (vector v) {
392 : 0 : vector result (v);
393 [ # # ][ # # ]: 0 : for (int i = 0; i < v.getSize (); i++) result.set (log (v.get (i)), i);
394 : 0 : return result;
395 : : }
396 : :
397 : 2 : vector log10 (vector v) {
398 : 2 : vector result (v);
399 [ + - ][ + + ]: 142 : for (int i = 0; i < v.getSize (); i++) result.set (log10 (v.get (i)), i);
400 : 2 : return result;
401 : : }
402 : :
403 : 0 : vector log2 (vector v) {
404 : 0 : vector result (v);
405 [ # # ][ # # ]: 0 : for (int i = 0; i < v.getSize (); i++) result.set (log2 (v.get (i)), i);
406 : 0 : return result;
407 : : }
408 : :
409 : 0 : vector pow (vector v, const nr_complex_t z) {
410 : 0 : vector result (v);
411 [ # # ][ # # ]: 0 : for (int i = 0; i < v.getSize (); i++) result.set (pow (v.get(i), z), i);
412 : 0 : return result;
413 : : }
414 : :
415 : 0 : vector pow (vector v, const nr_double_t d) {
416 : 0 : vector result (v);
417 [ # # ][ # # ]: 0 : for (int i = 0; i < v.getSize (); i++) result.set (pow (v.get(i), d), i);
418 : 0 : return result;
419 : : }
420 : :
421 : 0 : vector pow (const nr_complex_t z, vector v) {
422 : 0 : vector result (v);
423 [ # # ][ # # ]: 0 : for (int i = 0; i < v.getSize (); i++) result.set (pow (z, v.get (i)), i);
424 : 0 : return result;
425 : : }
426 : :
427 : 0 : vector pow (const nr_double_t d, vector v) {
428 : 0 : vector result (v);
429 [ # # ][ # # ]: 0 : for (int i = 0; i < v.getSize (); i++) result.set (pow (d, v.get (i)), i);
430 : 0 : return result;
431 : : }
432 : :
433 : 0 : vector pow (vector v1, vector v2) {
434 : 0 : int j, i, n, len, len1 = v1.getSize (), len2 = v2.getSize ();
435 [ # # ]: 0 : if (len1 >= len2) {
436 [ # # ]: 0 : assert (len1 % len2 == 0);
437 : 0 : len = len1;
438 : : } else {
439 [ # # ]: 0 : assert (len2 % len1 == 0);
440 : 0 : len = len2;
441 : : }
442 : 0 : vector res (len);
443 [ # # ]: 0 : for (j = i = n = 0; n < len; n++) {
444 [ # # ]: 0 : res (n) = pow (v1 (i), v2 (j));
445 [ # # ][ # # ]: 0 : if (++i >= len1) i = 0; if (++j >= len2) j = 0;
446 : : }
447 : 0 : return res;
448 : : }
449 : :
450 : 9 : vector sin (vector v) {
451 : 9 : vector result (v);
452 [ + - ][ + + ]: 6315 : for (int i = 0; i < v.getSize (); i++) result.set (sin (v.get (i)), i);
453 : 9 : return result;
454 : : }
455 : :
456 : 0 : vector asin (vector v) {
457 : 0 : vector result (v);
458 [ # # ][ # # ]: 0 : for (int i = 0; i < v.getSize (); i++) result.set (asin (v.get (i)), i);
459 : 0 : return result;
460 : : }
461 : :
462 : 0 : vector acos (vector v) {
463 : 0 : vector result (v);
464 [ # # ][ # # ]: 0 : for (int i = 0; i < v.getSize (); i++) result.set (acos (v.get (i)), i);
465 : 0 : return result;
466 : : }
467 : :
468 : 0 : vector cos (vector v) {
469 : 0 : vector result (v);
470 [ # # ][ # # ]: 0 : for (int i = 0; i < v.getSize (); i++) result.set (cos (v.get (i)), i);
471 : 0 : return result;
472 : : }
473 : :
474 : 0 : vector tan (vector v) {
475 : 0 : vector result (v);
476 [ # # ][ # # ]: 0 : for (int i = 0; i < v.getSize (); i++) result.set (tan (v.get (i)), i);
477 : 0 : return result;
478 : : }
479 : :
480 : 0 : vector atan (vector v) {
481 : 0 : vector result (v);
482 [ # # ][ # # ]: 0 : for (int i = 0; i < v.getSize (); i++) result.set (atan (v.get (i)), i);
483 : 0 : return result;
484 : : }
485 : :
486 : 0 : vector cot (vector v) {
487 : 0 : vector result (v);
488 [ # # ][ # # ]: 0 : for (int i = 0; i < v.getSize (); i++) result.set (cot (v.get (i)), i);
489 : 0 : return result;
490 : : }
491 : :
492 : 0 : vector acot (vector v) {
493 : 0 : vector result (v);
494 [ # # ][ # # ]: 0 : for (int i = 0; i < v.getSize (); i++) result.set (acot (v.get (i)), i);
495 : 0 : return result;
496 : : }
497 : :
498 : 0 : vector sinh (vector v) {
499 : 0 : vector result (v);
500 [ # # ][ # # ]: 0 : for (int i = 0; i < v.getSize (); i++) result.set (sinh (v.get (i)), i);
501 : 0 : return result;
502 : : }
503 : :
504 : 0 : vector asinh (vector v) {
505 : 0 : vector result (v);
506 [ # # ][ # # ]: 0 : for (int i = 0; i < v.getSize (); i++) result.set (asinh (v.get (i)), i);
507 : 0 : return result;
508 : : }
509 : :
510 : 0 : vector cosh (vector v) {
511 : 0 : vector result (v);
512 [ # # ][ # # ]: 0 : for (int i = 0; i < v.getSize (); i++) result.set (cosh (v.get (i)), i);
513 : 0 : return result;
514 : : }
515 : :
516 : 0 : vector sech (vector v) {
517 : 0 : vector result (v);
518 [ # # ][ # # ]: 0 : for (int i = 0; i < v.getSize (); i++) result.set (sech (v.get (i)), i);
519 : 0 : return result;
520 : : }
521 : :
522 : 0 : vector cosech (vector v) {
523 : 0 : vector result (v);
524 [ # # ][ # # ]: 0 : for (int i = 0; i < v.getSize (); i++) result.set (cosech (v.get (i)), i);
525 : 0 : return result;
526 : : }
527 : :
528 : 0 : vector acosh (vector v) {
529 : 0 : vector result (v);
530 [ # # ][ # # ]: 0 : for (int i = 0; i < v.getSize (); i++) result.set (acosh (v.get (i)), i);
531 : 0 : return result;
532 : : }
533 : :
534 : 0 : vector asech (vector v) {
535 : 0 : vector result (v);
536 [ # # ][ # # ]: 0 : for (int i = 0; i < v.getSize (); i++) result.set (asech (v.get (i)), i);
537 : 0 : return result;
538 : : }
539 : :
540 : 0 : vector tanh (vector v) {
541 : 0 : vector result (v);
542 [ # # ][ # # ]: 0 : for (int i = 0; i < v.getSize (); i++) result.set (tanh (v.get (i)), i);
543 : 0 : return result;
544 : : }
545 : :
546 : 0 : vector atanh (vector v) {
547 : 0 : vector result (v);
548 [ # # ][ # # ]: 0 : for (int i = 0; i < v.getSize (); i++) result.set (atanh (v.get (i)), i);
549 : 0 : return result;
550 : : }
551 : :
552 : 0 : vector coth (vector v) {
553 : 0 : vector result (v);
554 [ # # ][ # # ]: 0 : for (int i = 0; i < v.getSize (); i++) result.set (coth (v.get (i)), i);
555 : 0 : return result;
556 : : }
557 : :
558 : 0 : vector acoth (vector v) {
559 : 0 : vector result (v);
560 [ # # ][ # # ]: 0 : for (int i = 0; i < v.getSize (); i++) result.set (acoth (v.get (i)), i);
561 : 0 : return result;
562 : : }
563 : :
564 : : // converts impedance to reflexion coefficient
565 : 0 : vector ztor (vector v, nr_complex_t zref) {
566 : 0 : vector result (v);
567 [ # # ][ # # ]: 0 : for (int i = 0; i < v.getSize (); i++) result (i) = ztor (v (i), zref);
568 : 0 : return result;
569 : : }
570 : :
571 : : // converts admittance to reflexion coefficient
572 : 0 : vector ytor (vector v, nr_complex_t zref) {
573 : 0 : vector result (v);
574 [ # # ][ # # ]: 0 : for (int i = 0; i < v.getSize (); i++) result (i) = ytor (v (i), zref);
575 : 0 : return result;
576 : : }
577 : :
578 : : // converts reflexion coefficient to impedance
579 : 3 : vector rtoz (vector v, nr_complex_t zref) {
580 : 3 : vector result (v);
581 [ + - ][ + + ]: 764 : for (int i = 0; i < v.getSize (); i++) result (i) = rtoz (v (i), zref);
582 : 3 : return result;
583 : : }
584 : :
585 : : // converts reflexion coefficient to admittance
586 : 0 : vector rtoy (vector v, nr_complex_t zref) {
587 : 0 : vector result (v);
588 [ # # ][ # # ]: 0 : for (int i = 0; i < v.getSize (); i++) result (i) = rtoy (v (i), zref);
589 : 0 : return result;
590 : : }
591 : :
592 : : // differentiates 'var' with respect to 'dep' exactly 'n' times
593 : 2 : vector diff (vector var, vector dep, int n) {
594 : 2 : int k, xi, yi, exchange = 0;
595 [ + - ][ + - ]: 2 : vector x, y;
596 : : // exchange dependent and independent variable if necessary
597 [ - + ]: 2 : if (var.getSize () < dep.getSize ()) {
598 [ # # ][ # # ]: 0 : x = vector (var);
599 [ # # ][ # # ]: 0 : y = vector (dep);
600 : 0 : exchange++;
601 : : }
602 : : else {
603 [ + - ][ + - ]: 2 : x = vector (dep);
604 [ + - ][ + - ]: 2 : y = vector (var);
605 : : }
606 [ + - ][ - + ]: 2 : assert (y.getSize () % x.getSize () == 0 && x.getSize () >= 2);
[ - + ]
607 [ + - ]: 2 : vector result (y);
608 : 2 : nr_complex_t c;
609 : :
610 [ + + ]: 4 : for (k = 0; k < n; k++) { // differentiate n times
611 [ + + ]: 918 : for (yi = xi = 0; yi < y.getSize (); yi++, xi++) {
612 [ - + ]: 916 : if (xi == x.getSize ()) xi = 0; // roll through independent vector
613 [ + + ]: 916 : if (xi == 0) {
614 [ + - ][ + - ]: 2 : c = (y.get (yi + 1) - y.get (yi)) / (x.get (xi + 1) - x.get (xi));
[ + - ]
615 [ + + ]: 914 : } else if (xi == x.getSize () - 1) {
616 [ + - ][ + - ]: 2 : c = (y.get (yi) - y.get (yi - 1)) / (x.get (xi) - x.get (xi - 1));
[ + - ]
617 : : }
618 : : else {
619 : : c =
620 [ + - ][ + - ]: 912 : ((y.get (yi) - y.get (yi - 1)) / (x.get (xi) - x.get (xi - 1)) +
621 [ + - ][ + - ]: 1824 : (y.get (yi + 1) - y.get (yi)) / (x.get (xi + 1) - x.get (xi))) /
[ + - ][ + - ]
622 : 912 : 2.0;
623 : : }
624 [ - + ][ # # ]: 916 : result.set (exchange ? 1.0 / c : c, yi);
[ - + ][ # # ]
625 : : }
626 : 2 : y = result;
627 : : }
628 [ + - ][ + - ]: 2 : return result;
629 : : }
630 : :
631 : 0 : vector vector::operator=(const nr_complex_t c) {
632 [ # # ]: 0 : for (int i = 0; i < size; i++) data[i] = c;
633 : 0 : return *this;
634 : : }
635 : :
636 : 1070 : vector vector::operator=(const nr_double_t d) {
637 [ + + ]: 4570 : for (int i = 0; i < size; i++) data[i] = d;
638 : 1070 : return *this;
639 : : }
640 : :
641 : 2 : vector vector::operator+=(vector v) {
642 : 2 : int i, n, len = v.getSize ();
643 [ - + ]: 2 : assert (size % len == 0);
644 [ + + ][ + + ]: 661 : for (i = n = 0; i < size; i++) { data[i] += v (n); if (++n >= len) n = 0; }
645 : 2 : return *this;
646 : : }
647 : :
648 : 0 : vector vector::operator+=(const nr_complex_t c) {
649 [ # # ]: 0 : for (int i = 0; i < size; i++) data[i] += c;
650 : 0 : return *this;
651 : : }
652 : :
653 : 1 : vector vector::operator+=(const nr_double_t d) {
654 [ + + ]: 101 : for (int i = 0; i < size; i++) data[i] += d;
655 : 1 : return *this;
656 : : }
657 : :
658 : 2 : vector operator+(vector v1, vector v2) {
659 : 2 : int len1 = v1.getSize (), len2 = v2.getSize ();
660 : 2 : vector result;
661 [ + - ]: 2 : if (len1 >= len2) {
662 : 2 : result = v1;
663 [ + - ][ + - ]: 2 : result += v2;
[ + - ][ + - ]
664 : : } else {
665 : 0 : result = v2;
666 [ # # ][ # # ]: 0 : result += v1;
[ # # ][ # # ]
667 : : }
668 : 2 : return result;
669 : : }
670 : :
671 : 0 : vector operator+(vector v, const nr_complex_t c) {
672 : 0 : vector result (v);
673 [ # # ][ # # ]: 0 : result += c;
674 : 0 : return result;
675 : : }
676 : :
677 : 0 : vector operator+(const nr_complex_t c, vector v) {
678 [ # # ]: 0 : return v + c;
679 : : }
680 : :
681 : 1 : vector operator+(vector v, const nr_double_t d) {
682 : 1 : vector result (v);
683 [ + - ][ + - ]: 1 : result += d;
684 : 1 : return result;
685 : : }
686 : :
687 : 0 : vector operator+(const nr_double_t d, vector v) {
688 [ # # ]: 0 : return v + d;
689 : : }
690 : :
691 : 10 : vector vector::operator-() {
692 : 10 : vector result (size);
693 [ + + ]: 3305 : for (int i = 0; i < size; i++) result (i) = -data[i];
694 : 10 : return result;
695 : : }
696 : :
697 : 18 : vector vector::operator-=(vector v) {
698 : 18 : int i, n, len = v.getSize ();
699 [ - + ]: 18 : assert (size % len == 0);
700 [ + + ][ + + ]: 36758 : for (i = n = 0; i < size; i++) { data[i] -= v (n); if (++n >= len) n = 0; }
701 : 18 : return *this;
702 : : }
703 : :
704 : 0 : vector vector::operator-=(const nr_complex_t c) {
705 [ # # ]: 0 : for (int i = 0; i < size; i++) data[i] -= c;
706 : 0 : return *this;
707 : : }
708 : :
709 : 14 : vector vector::operator-=(const nr_double_t d) {
710 [ + + ]: 3300 : for (int i = 0; i < size; i++) data[i] -= d;
711 : 14 : return *this;
712 : : }
713 : :
714 : 18 : vector operator-(vector v1, vector v2) {
715 : 18 : int len1 = v1.getSize (), len2 = v2.getSize ();
716 : 18 : vector result;
717 [ + - ]: 18 : if (len1 >= len2) {
718 : 18 : result = v1;
719 [ + - ][ + - ]: 18 : result -= v2;
[ + - ][ + - ]
720 : : } else {
721 [ # # ][ # # ]: 0 : result = -v2;
722 [ # # ][ # # ]: 0 : result += v1;
[ # # ][ # # ]
723 : : }
724 : 18 : return result;
725 : : }
726 : :
727 : 0 : vector operator-(vector v, const nr_complex_t c) {
728 : 0 : vector result (v);
729 [ # # ][ # # ]: 0 : result -= c;
730 : 0 : return result;
731 : : }
732 : :
733 : 14 : vector operator-(vector v, const nr_double_t d) {
734 : 14 : vector result (v);
735 [ + - ][ + - ]: 14 : result -= d;
736 : 14 : return result;
737 : : }
738 : :
739 : 0 : vector operator-(const nr_complex_t c, vector v) {
740 : 0 : vector result (-v);
741 [ # # ][ # # ]: 0 : result += c;
742 : 0 : return result;
743 : : }
744 : :
745 : 0 : vector operator-(const nr_double_t d, vector v) {
746 : 0 : vector result (-v);
747 [ # # ][ # # ]: 0 : result += d;
748 : 0 : return result;
749 : : }
750 : :
751 : 0 : vector vector::operator*=(vector v) {
752 : 0 : int i, n, len = v.getSize ();
753 [ # # ]: 0 : assert (size % len == 0);
754 [ # # ][ # # ]: 0 : for (i = n = 0; i < size; i++) { data[i] *= v (n); if (++n >= len) n = 0; }
755 : 0 : return *this;
756 : : }
757 : :
758 : 0 : vector vector::operator*=(const nr_complex_t c) {
759 [ # # ]: 0 : for (int i = 0; i < size; i++) data[i] *= c;
760 : 0 : return *this;
761 : : }
762 : :
763 : 25 : vector vector::operator*=(const nr_double_t d) {
764 [ + + ]: 12012 : for (int i = 0; i < size; i++) data[i] *= d;
765 : 25 : return *this;
766 : : }
767 : :
768 : 0 : vector operator*(vector v1, vector v2) {
769 : 0 : int len1 = v1.getSize (), len2 = v2.getSize ();
770 : 0 : vector result;
771 [ # # ]: 0 : if (len1 >= len2) {
772 : 0 : result = v1;
773 [ # # ][ # # ]: 0 : result *= v2;
[ # # ][ # # ]
774 : : } else {
775 : 0 : result = v2;
776 [ # # ][ # # ]: 0 : result *= v1;
[ # # ][ # # ]
777 : : }
778 : 0 : return result;
779 : : }
780 : :
781 : 0 : vector operator*(vector v, const nr_complex_t c) {
782 : 0 : vector result (v);
783 [ # # ][ # # ]: 0 : result *= c;
784 : 0 : return result;
785 : : }
786 : :
787 : 25 : vector operator*(vector v, const nr_double_t d) {
788 : 25 : vector result (v);
789 [ + - ][ + - ]: 25 : result *= d;
790 : 25 : return result;
791 : : }
792 : :
793 : 0 : vector operator*(const nr_complex_t c, vector v) {
794 [ # # ]: 0 : return v * c;
795 : : }
796 : :
797 : 25 : vector operator*(const nr_double_t d, vector v) {
798 [ + - ]: 25 : return v * d;
799 : : }
800 : :
801 : 1083 : vector vector::operator/=(vector v) {
802 : 1083 : int i, n, len = v.getSize ();
803 [ - + ]: 1083 : assert (size % len == 0);
804 [ + + ][ + + ]: 7101 : for (i = n = 0; i < size; i++) { data[i] /= v (n); if (++n >= len) n = 0; }
805 : 1083 : return *this;
806 : : }
807 : :
808 : 3 : vector vector::operator/=(const nr_complex_t c) {
809 [ + + ]: 1377 : for (int i = 0; i < size; i++) data[i] /= c;
810 : 3 : return *this;
811 : : }
812 : :
813 : 7 : vector vector::operator/=(const nr_double_t d) {
814 [ + + ]: 2366 : for (int i = 0; i < size; i++) data[i] /= d;
815 : 7 : return *this;
816 : : }
817 : :
818 : 13 : vector operator/(vector v1, vector v2) {
819 : 13 : int len1 = v1.getSize (), len2 = v2.getSize ();
820 : 13 : vector result;
821 [ + - ]: 13 : if (len1 >= len2) {
822 [ - + ]: 13 : assert (len1 % len2 == 0);
823 : 13 : result = v1;
824 [ + - ][ + - ]: 13 : result /= v2;
[ + - ][ + - ]
825 : : } else {
826 [ # # ]: 0 : assert (len2 % len1 == 0);
827 [ # # ][ # # ]: 0 : result = 1 / v2;
[ # # ][ # # ]
828 [ # # ][ # # ]: 0 : result *= v1;
[ # # ][ # # ]
829 : : }
830 : 13 : return result;
831 : : }
832 : :
833 : 3 : vector operator/(vector v, const nr_complex_t c) {
834 : 3 : vector result (v);
835 [ + - ][ + - ]: 3 : result /= c;
836 : 3 : return result;
837 : : }
838 : :
839 : 7 : vector operator/(vector v, const nr_double_t d) {
840 : 7 : vector result (v);
841 [ + - ][ + - ]: 7 : result /= d;
842 : 7 : return result;
843 : : }
844 : :
845 : 0 : vector operator/(const nr_complex_t c, vector v) {
846 : 0 : vector result (v);
847 [ # # ][ # # ]: 0 : result = c;
848 [ # # ][ # # ]: 0 : result /= v;
[ # # ][ # # ]
849 : 0 : return result;
850 : : }
851 : :
852 : 1070 : vector operator/(const nr_double_t d, vector v) {
853 : 1070 : vector result (v);
854 [ + - ][ + - ]: 1070 : result = d;
855 [ + - ][ + - ]: 1070 : result /= v;
[ + - ][ + - ]
856 : 1070 : return result;
857 : : }
858 : :
859 : 0 : vector operator%(vector v, const nr_complex_t z) {
860 : 0 : int len = v.getSize ();
861 : 0 : vector result (len);
862 [ # # ][ # # ]: 0 : for (int i = 0; i < len; i++) result (i) = v (i) % z;
863 : 0 : return result;
864 : : }
865 : :
866 : 0 : vector operator%(vector v, const nr_double_t d) {
867 : 0 : int len = v.getSize ();
868 : 0 : vector result (len);
869 [ # # ][ # # ]: 0 : for (int i = 0; i < len; i++) result (i) = v (i) % d;
870 : 0 : return result;
871 : : }
872 : :
873 : 0 : vector operator%(const nr_complex_t z, vector v) {
874 : 0 : int len = v.getSize ();
875 : 0 : vector result (len);
876 [ # # ][ # # ]: 0 : for (int i = 0; i < len; i++) result (i) = z % v (i);
877 : 0 : return result;
878 : : }
879 : :
880 : 0 : vector operator%(const nr_double_t d, vector v) {
881 : 0 : int len = v.getSize ();
882 : 0 : vector result (len);
883 [ # # ][ # # ]: 0 : for (int i = 0; i < len; i++) result (i) = d % v (i);
884 : 0 : return result;
885 : : }
886 : :
887 : 0 : vector operator%(vector v1, vector v2) {
888 : 0 : int j, i, n, len, len1 = v1.getSize (), len2 = v2.getSize ();
889 [ # # ]: 0 : if (len1 >= len2) {
890 [ # # ]: 0 : assert (len1 % len2 == 0);
891 : 0 : len = len1;
892 : : } else {
893 [ # # ]: 0 : assert (len2 % len1 == 0);
894 : 0 : len = len2;
895 : : }
896 : 0 : vector res (len);
897 [ # # ]: 0 : for (j = i = n = 0; n < len; n++) {
898 [ # # ]: 0 : res (n) = v1 (i) % v2 (j);
899 [ # # ][ # # ]: 0 : if (++i >= len1) i = 0; if (++j >= len2) j = 0;
900 : : }
901 : 0 : return res;
902 : : }
903 : :
904 : : /* This function reverses the order of the data list. */
905 : 0 : void vector::reverse (void) {
906 : : nr_complex_t * buffer = (nr_complex_t *)
907 : 0 : malloc (sizeof (nr_complex_t) * size);
908 [ # # ]: 0 : for (int i = 0; i < size; i++) buffer[i] = data[size - 1 - i];
909 : 0 : free (data);
910 : 0 : data = buffer;
911 : 0 : capacity = size;
912 : 0 : }
913 : :
914 : : // Sets the origin (the analysis) of the vector.
915 : 1051 : void vector::setOrigin (char * n) {
916 [ - + ]: 1051 : if (origin) free (origin);
917 [ + - ]: 1051 : origin = n ? strdup (n) : NULL;
918 : 1051 : }
919 : :
920 : : // Returns the origin (the analysis) of the vector.
921 : 85863 : char * vector::getOrigin (void) {
922 : 85863 : return origin;
923 : : }
924 : :
925 : : /* The function returns the number of entries with the given value
926 : : deviating no more than the given epsilon. */
927 : 0 : int vector::contains (nr_complex_t val, nr_double_t eps) {
928 : 0 : int count = 0;
929 [ # # ]: 0 : for (int i = 0; i < size; i++) {
930 [ # # ]: 0 : if (abs (data[i] - val) <= eps) count++;
931 : : }
932 : 0 : return count;
933 : : }
934 : :
935 : : // Sorts the vector either in ascending or descending order.
936 : 0 : void vector::sort (bool ascending) {
937 : 0 : nr_complex_t t;
938 [ # # ]: 0 : for (int i = 0; i < size; i++) {
939 [ # # ]: 0 : for (int n = 0; n < size - 1; n++) {
940 [ # # ][ # # ]: 0 : if (ascending ? data[n] > data[n+1] : data[n] < data[n+1]) {
[ # # ][ # # ]
941 : 0 : t = data[n];
942 : 0 : data[n] = data[n+1];
943 : 0 : data[n+1] = t;
944 : : }
945 : : }
946 : : }
947 : 0 : }
948 : :
949 : : /* The function creates a linear stepped vector of values starting at
950 : : the given start value, ending with the given stop value and
951 : : containing points elements. */
952 : 134 : vector linspace (nr_double_t start, nr_double_t stop, int points) {
953 : 134 : vector result (points);
954 : 134 : nr_double_t val, step = (stop - start) / (points - 1);
955 [ + + ]: 60216 : for (int i = 0; i < points; i++) {
956 : 60082 : val = start + (i * step);
957 [ + + ][ + + ]: 60082 : if (i != 0 && fabs (val) < fabs (step) / 4 && fabs (val) < std::numeric_limits<nr_double_t>::epsilon())
[ + - ][ + + ]
958 : 4 : val = 0.0;
959 : 60082 : result.set (val, i);
960 : : }
961 : 134 : return result;
962 : : }
963 : :
964 : : /* The function creates a logarithmic stepped vector of values
965 : : starting at the given start value, ending with the given stop value
966 : : and containing points elements. */
967 : 25 : vector logspace (nr_double_t start, nr_double_t stop, int points) {
968 [ - + ]: 25 : assert (start * stop > 0);
969 : 25 : vector result (points);
970 : : nr_double_t step, first, last, d;
971 : :
972 : : // ensure the last value being larger than the first
973 [ - + ]: 25 : if (fabs (start) > fabs (stop)) {
974 : 0 : first = fabs (stop);
975 : 0 : last = fabs (start);
976 : : }
977 : : else {
978 : 25 : first = fabs (start);
979 : 25 : last = fabs (stop);
980 : : }
981 : : // check direction and sign of values
982 [ - + ]: 25 : d = fabs (start) > fabs (stop) ? -1 : 1;
983 : : // compute logarithmic step size
984 : 25 : step = (::log (last) - ::log (first)) / (points - 1);
985 [ + + ]: 7430 : for (int i = 0, j = points - 1; i < points; i++, j--) {
986 [ + - ]: 7405 : if (d > 0)
987 : 7405 : result.set (start * ::exp (step * i), i);
988 : : else
989 : 0 : result.set (stop * ::exp (step * i), j);
990 : : }
991 : 25 : return result;
992 : : }
993 : :
994 : 0 : vector cumsum (vector v) {
995 [ # # ]: 0 : vector result (v);
996 : 0 : nr_complex_t val (0.0);
997 [ # # ]: 0 : for (int i = 0; i < v.getSize (); i++) {
998 : 0 : val += v.get (i);
999 : 0 : result.set (val, i);
1000 : : }
1001 : 0 : return result;
1002 : : }
1003 : :
1004 : 0 : vector cumavg (vector v) {
1005 [ # # ]: 0 : vector result (v);
1006 : 0 : nr_complex_t val (0.0);
1007 [ # # ]: 0 : for (int i = 0; i < v.getSize (); i++) {
1008 : 0 : val = (val * (nr_double_t) i + v.get (i)) / (i + 1.0);
1009 : 0 : result.set (val, i);
1010 : : }
1011 : 0 : return result;
1012 : : }
1013 : :
1014 : 0 : vector cumprod (vector v) {
1015 [ # # ]: 0 : vector result (v);
1016 : 0 : nr_complex_t val (1.0);
1017 [ # # ]: 0 : for (int i = 0; i < v.getSize (); i++) {
1018 : 0 : val *= v.get (i);
1019 : 0 : result.set (val, i);
1020 : : }
1021 : 0 : return result;
1022 : : }
1023 : :
1024 : 0 : vector ceil (vector v) {
1025 : 0 : vector result (v);
1026 [ # # ][ # # ]: 0 : for (int i = 0; i < v.getSize (); i++) result.set (ceil (v.get (i)), i);
1027 : 0 : return result;
1028 : : }
1029 : :
1030 : 0 : vector fix (vector v) {
1031 : 0 : vector result (v);
1032 [ # # ][ # # ]: 0 : for (int i = 0; i < v.getSize (); i++) result.set (fix (v.get (i)), i);
1033 : 0 : return result;
1034 : : }
1035 : :
1036 : 0 : vector floor (vector v) {
1037 : 0 : vector result (v);
1038 [ # # ][ # # ]: 0 : for (int i = 0; i < v.getSize (); i++) result.set (floor (v.get (i)), i);
1039 : 0 : return result;
1040 : : }
1041 : :
1042 : 0 : vector round (vector v) {
1043 : 0 : vector result (v);
1044 [ # # ][ # # ]: 0 : for (int i = 0; i < v.getSize (); i++) result.set (round (v.get (i)), i);
1045 : 0 : return result;
1046 : : }
1047 : :
1048 : 0 : vector sqr (vector v) {
1049 : 0 : vector result (v);
1050 [ # # ][ # # ]: 0 : for (int i = 0; i < v.getSize (); i++) result.set (sqr (v.get (i)), i);
1051 : 0 : return result;
1052 : : }
1053 : :
1054 : 0 : vector step (vector v) {
1055 : 0 : vector result (v);
1056 [ # # ][ # # ]: 0 : for (int i = 0; i < v.getSize (); i++) result.set (step (v.get (i)), i);
1057 : 0 : return result;
1058 : : }
1059 : :
1060 : 0 : static nr_double_t integrate_n (vector v) { /* using trapezoidal rule */
1061 : 0 : nr_double_t result = 0.0;
1062 [ # # ]: 0 : for (int i = 1; i < v.getSize () - 1; i++) result += norm (v.get (i));
1063 : 0 : result += 0.5 * norm (v.get (0));
1064 : 0 : result += 0.5 * norm (v.get (v.getSize () - 1));
1065 : 0 : return result;
1066 : : }
1067 : :
1068 : 0 : nr_double_t vector::rms (void) {
1069 [ # # ]: 0 : nr_double_t result = std::sqrt (integrate_n (*this) / getSize ());
1070 : 0 : return result;
1071 : : }
1072 : :
1073 : 0 : nr_double_t vector::variance (void) {
1074 : 0 : nr_double_t result = 0.0;
1075 [ # # ][ # # ]: 0 : nr_complex_t average = avg (*this);
[ # # ]
1076 [ # # ][ # # ]: 0 : for (int i = 0; i < getSize (); i++) result += norm (get (i) - average);
[ # # ]
1077 [ # # ]: 0 : if (getSize () > 1)
1078 : 0 : return result / (getSize () - 1);
1079 : 0 : return 0.0;
1080 : : }
1081 : :
1082 : 0 : nr_double_t vector::stddev (void) {
1083 : 0 : return std::sqrt (variance ());
1084 : : }
1085 : :
1086 : 0 : vector erf (vector v) {
1087 : 0 : vector result (v);
1088 [ # # ][ # # ]: 0 : for (int i = 0; i < v.getSize (); i++) result.set (erf (v.get (i)), i);
1089 : 0 : return result;
1090 : : }
1091 : :
1092 : 0 : vector erfc (vector v) {
1093 : 0 : vector result (v);
1094 [ # # ][ # # ]: 0 : for (int i = 0; i < v.getSize (); i++) result.set (erfc (v.get (i)), i);
1095 : 0 : return result;
1096 : : }
1097 : :
1098 : 0 : vector erfinv (vector v) {
1099 : 0 : vector result (v);
1100 [ # # ][ # # ]: 0 : for (int i = 0; i < v.getSize (); i++) result.set (erfinv (v.get (i)), i);
1101 : 0 : return result;
1102 : : }
1103 : :
1104 : 0 : vector erfcinv (vector v) {
1105 : 0 : vector result (v);
1106 [ # # ][ # # ]: 0 : for (int i = 0; i < v.getSize (); i++) result.set (erfcinv (v.get (i)), i);
1107 : 0 : return result;
1108 : : }
1109 : :
1110 : 0 : vector i0 (vector v) {
1111 : 0 : vector result (v);
1112 [ # # ][ # # ]: 0 : for (int i = 0; i < v.getSize (); i++) result.set (i0 (v.get (i)), i);
1113 : 0 : return result;
1114 : : }
1115 : :
1116 : 0 : vector jn (const int n, vector v) {
1117 : 0 : vector result (v);
1118 [ # # ][ # # ]: 0 : for (int i = 0; i < v.getSize (); i++) result.set (jn (n, v.get (i)), i);
1119 : 0 : return result;
1120 : : }
1121 : :
1122 : 0 : vector yn (const int n, vector v) {
1123 : 0 : vector result (v);
1124 [ # # ][ # # ]: 0 : for (int i = 0; i < v.getSize (); i++) result.set (yn (n, v.get (i)), i);
1125 : 0 : return result;
1126 : : }
1127 : :
1128 : 0 : vector polar (const nr_complex_t a, vector v) {
1129 : 0 : vector result (v);
1130 [ # # ][ # # ]: 0 : for (int i = 0; i < v.getSize (); i++) result.set (qucs::polar (a, v.get (i)), i);
1131 : 0 : return result;
1132 : : }
1133 : :
1134 : 0 : vector polar (vector v, const nr_complex_t p) {
1135 : 0 : vector result (v);
1136 [ # # ][ # # ]: 0 : for (int i = 0; i < v.getSize (); i++) result.set (qucs::polar (v.get (i), p), i);
1137 : 0 : return result;
1138 : : }
1139 : :
1140 : 0 : vector polar (vector a, vector p) {
1141 : 0 : int j, i, n, len, len1 = a.getSize (), len2 = p.getSize ();
1142 [ # # ]: 0 : if (len1 >= len2) {
1143 [ # # ]: 0 : assert (len1 % len2 == 0);
1144 : 0 : len = len1;
1145 : : } else {
1146 [ # # ]: 0 : assert (len2 % len1 == 0);
1147 : 0 : len = len2;
1148 : : }
1149 : 0 : vector res (len);
1150 [ # # ]: 0 : for (j = i = n = 0; n < len; n++) {
1151 [ # # ]: 0 : res (n) = qucs::polar (a (i), p (j));
1152 [ # # ][ # # ]: 0 : if (++i >= len1) i = 0; if (++j >= len2) j = 0;
1153 : : }
1154 : 0 : return res;
1155 : : }
1156 : :
1157 : 0 : vector atan2 (const nr_double_t y, vector v) {
1158 : 0 : vector result (v);
1159 [ # # ]: 0 : for (int i = 0; i < v.getSize (); i++)
1160 [ # # ]: 0 : result.set (atan2 (y, v.get (i)), i);
1161 : 0 : return result;
1162 : : }
1163 : :
1164 : 0 : vector atan2 (vector v, const nr_double_t x) {
1165 : 0 : vector result (v);
1166 [ # # ]: 0 : for (int i = 0; i < v.getSize (); i++)
1167 [ # # ]: 0 : result.set (atan2 (v.get (i), x) , i);
1168 : 0 : return result;
1169 : : }
1170 : :
1171 : 0 : vector atan2 (vector y, vector x) {
1172 : 0 : int j, i, n, len, len1 = y.getSize (), len2 = x.getSize ();
1173 [ # # ]: 0 : if (len1 >= len2) {
1174 [ # # ]: 0 : assert (len1 % len2 == 0);
1175 : 0 : len = len1;
1176 : : } else {
1177 [ # # ]: 0 : assert (len2 % len1 == 0);
1178 : 0 : len = len2;
1179 : : }
1180 : 0 : vector res (len);
1181 [ # # ]: 0 : for (j = i = n = 0; n < len; n++) {
1182 [ # # ]: 0 : res (n) = atan2 (y (i), x (j));
1183 [ # # ][ # # ]: 0 : if (++i >= len1) i = 0; if (++j >= len2) j = 0;
1184 : : }
1185 : 0 : return res;
1186 : : }
1187 : :
1188 : 0 : vector w2dbm (vector v) {
1189 : 0 : vector result (v);
1190 [ # # ]: 0 : for (int i = 0; i < v.getSize (); i++)
1191 [ # # ]: 0 : result.set (10.0 * log10 (v.get (i) / 0.001), i);
1192 : 0 : return result;
1193 : : }
1194 : :
1195 : 0 : vector dbm2w (vector v) {
1196 : 0 : vector result (v);
1197 [ # # ]: 0 : for (int i = 0; i < v.getSize (); i++)
1198 [ # # ]: 0 : result.set (0.001 * pow (10.0 , v.get (i) / 10.0), i);
1199 : 0 : return result;
1200 : : }
1201 : :
1202 : 0 : nr_double_t integrate (vector v, const nr_double_t h) {
1203 : 0 : nr_double_t s = real (v.get (0) ) / 2;
1204 [ # # ]: 0 : for (int i = 1; i < v.getSize () - 1; i++)
1205 : 0 : s += real (v.get (i));
1206 : 0 : return (s + real (v.get (v.getSize () - 1) ) / 2) * h;
1207 : : }
1208 : :
1209 : 0 : nr_complex_t integrate (vector v, const nr_complex_t h) {
1210 : 0 : nr_complex_t s;
1211 : 0 : s = v.get (0) / 2.0;
1212 [ # # ]: 0 : for (int i = 1; i < v.getSize () - 1; i++)
1213 : 0 : s += v.get (i);
1214 : 0 : return (s + v.get (v.getSize () - 1) / 2.0) * h;
1215 : : }
1216 : :
1217 : 0 : vector dbm (vector v, const nr_complex_t z) {
1218 : 0 : vector result (v);
1219 [ # # ]: 0 : for (int i = 0; i < v.getSize (); i++)
1220 [ # # ][ # # ]: 0 : result.set (10.0 * log10 (norm (v.get (i)) / conj (z) / 0.001), i);
[ # # ]
1221 : 0 : return result;
1222 : : }
1223 : :
1224 : 0 : vector runavg (const nr_complex_t x, const int n) {
1225 : 0 : vector result (n);
1226 [ # # ]: 0 : for (int i = 0; i < n; i++) result.set (x, i);
1227 : 0 : return result;
1228 : : }
1229 : :
1230 : 0 : vector runavg (vector v, const int n) {
1231 : 0 : nr_complex_t s (0.0), y;
1232 : 0 : int len = v.getSize () - n + 1, i;
1233 [ # # ]: 0 : vector result (len);
1234 [ # # ]: 0 : for (i = 0; i < n; i++) s += v.get (i);
1235 : 0 : y = s / (nr_double_t) n; // first running average value
1236 : 0 : result.set (y, 0);
1237 [ # # ]: 0 : for (i = 0; i < len - 1; i++) {
1238 [ # # ]: 0 : y += (v.get (i + n) - v.get (i)) / (nr_double_t) n;
1239 : 0 : result.set (y, i + 1);
1240 : : }
1241 : 0 : return result;
1242 : : }
1243 : :
1244 : : #ifdef DEBUG
1245 : : // Debug function: Prints the vector object.
1246 : 0 : void vector::print (void) {
1247 [ # # ]: 0 : for (int r = 0; r < size; r++) {
1248 : : fprintf (stderr, "%+.2e%+.2ei\n", (double) real (get (r)),
1249 [ # # ]: 0 : (double) imag (get (r)));
1250 : : }
1251 : 0 : }
1252 : : #endif /* DEBUG */
1253 : :
1254 : : } // namespace qucs
|
