Branch data Line data Source code
1 : : /*
2 : : * matrix.cpp - matrix class implementation
3 : : *
4 : : * Copyright (C) 2003-2009 Stefan Jahn <stefan@lkcc.org>
5 : : *
6 : : * This is free software; you can redistribute it and/or modify
7 : : * it under the terms of the GNU General Public License as published by
8 : : * the Free Software Foundation; either version 2, or (at your option)
9 : : * any later version.
10 : : *
11 : : * This software is distributed in the hope that it will be useful,
12 : : * but WITHOUT ANY WARRANTY; without even the implied warranty of
13 : : * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
14 : : * GNU General Public License for more details.
15 : : *
16 : : * You should have received a copy of the GNU General Public License
17 : : * along with this package; see the file COPYING. If not, write to
18 : : * the Free Software Foundation, Inc., 51 Franklin Street - Fifth Floor,
19 : : * Boston, MA 02110-1301, USA.
20 : : *
21 : : * $Id$
22 : : *
23 : : */
24 : : /*!\file matrix.cpp
25 : : \brief Dense matrix class implementation
26 : :
27 : : References:
28 : :
29 : : [1] Power Waves and the Scattering Matrix
30 : : Kurokawa, K.
31 : : Microwave Theory and Techniques, IEEE Transactions on,
32 : : Vol.13, Iss.2, Mar 1965
33 : : Pages: 194- 202
34 : :
35 : : [2] A Rigorous Technique for Measuring the Scattering Matrix of
36 : : a Multiport Device with a 2-Port Network Analyzer
37 : : John C. TIPPET, Ross A. SPECIALE
38 : : Microwave Theory and Techniques, IEEE Transactions on,
39 : : Vol.82, Iss.5, May 1982
40 : : Pages: 661- 666
41 : :
42 : : [3] Comments on "A Rigorous Techique for Measuring the Scattering
43 : : Matrix of a Multiport Device with a Two-Port Network Analyzer"
44 : : Dropkin, H.
45 : : Microwave Theory and Techniques, IEEE Transactions on,
46 : : Vol. 83, Iss.1, Jan 1983
47 : : Pages: 79 - 81
48 : :
49 : : [4] Arbitrary Impedance
50 : : "Accurate Measurements In Almost Any
51 : : Impedance Environment"
52 : : in Scropion Application note
53 : : Anritsu
54 : : online(2007/07/30) http://www.eu.anritsu.com/files/11410-00284B.pdf
55 : :
56 : : [5] Conversions between S, Z, Y, H, ABCD, and T parameters
57 : : which are valid for complex source and load impedances
58 : : Frickey, D.A.
59 : : Microwave Theory and Techniques, IEEE Transactions on
60 : : Vol. 42, Iss. 2, Feb 1994
61 : : pages: 205 - 211
62 : : doi: 10.1109/22.275248
63 : :
64 : : [6] Comments on "Conversions between S, Z, Y, h, ABCD,
65 : : and T parameters which are valid for complex source and load impedances" [and reply]
66 : : Marks, R.B.; Williams, D.F.; Frickey, D.A.
67 : : Microwave Theory and Techniques, IEEE Transactions on,
68 : : Vol.43, Iss.4, Apr 1995
69 : : Pages: 914- 915
70 : : doi: 10.1109/22.375247
71 : :
72 : : [7] Wave Techniques for Noise Modeling and Measurement
73 : : S. W. Wedge and D. B. Rutledge,
74 : : IEEE Transactions on Microwave Theory and Techniques,
75 : : vol. 40, no. 11, Nov. 1992.
76 : : pages 2004-2012,
77 : : doi: 10.1109/22.168757
78 : : Author copy online (2007/07/31)
79 : : http://authors.library.caltech.edu/6226/01/WEDieeetmtt92.pdf
80 : :
81 : : */
82 : : #if HAVE_CONFIG_H
83 : : # include <config.h>
84 : : #endif
85 : :
86 : : #include <assert.h>
87 : : #include <stdio.h>
88 : : #include <cstdlib>
89 : : #include <string.h>
90 : : #include <cmath>
91 : :
92 : : #include "logging.h"
93 : : #include "object.h"
94 : : #include "complex.h"
95 : : #include "vector.h"
96 : : #include "matrix.h"
97 : :
98 : : namespace qucs {
99 : :
100 : : /*!\brief Create an empty matrix
101 : :
102 : : Constructor creates an unnamed instance of the matrix class.
103 : : */
104 : 13273 : matrix::matrix () {
105 : 13273 : rows = 0;
106 : 13273 : cols = 0;
107 : 13273 : data = NULL;
108 : 13273 : }
109 : :
110 : : /*!\brief Creates a square matrix
111 : :
112 : : Constructor creates an unnamed instance of the matrix class with a
113 : : certain number of rows and columns. Particularly creates a square matrix.
114 : : \param[in] s number of rows or colums of square matrix
115 : : \todo Why not s const?
116 : : */
117 : 5950 : matrix::matrix (int s) {
118 : 5950 : rows = cols = s;
119 [ + - ][ + + ]: 108950 : data = (s > 0) ? new nr_complex_t[s * s] : NULL;
[ # # ][ # # ]
120 : 5950 : }
121 : :
122 : : /* \brief Creates a matrix
123 : :
124 : : Constructor creates an unnamed instance of the matrix class with a
125 : : certain number of rows and columns.
126 : : \param[in] r number of rows
127 : : \param[in] c number of column
128 : : \todo Why not r and c constant
129 : : \todo Assert r >= 0 and c >= 0
130 : : */
131 : 15913 : matrix::matrix (int r, int c) {
132 : 15913 : rows = r;
133 : 15913 : cols = c;
134 [ + - ][ + - ]: 257152 : data = (r > 0 && c > 0) ? new nr_complex_t[r * c] : NULL;
[ + + ][ # # ]
[ # # ][ # # ]
135 : 15913 : }
136 : :
137 : : /* \brief copy constructor
138 : :
139 : : The copy constructor creates a new instance based on the given
140 : : matrix object.
141 : : \todo Add assert tests
142 : : */
143 : 13130 : matrix::matrix (const matrix & m) {
144 : 13130 : rows = m.rows;
145 : 13130 : cols = m.cols;
146 : 13130 : data = NULL;
147 : :
148 : : // copy matrix elements
149 [ + - ][ + - ]: 13130 : if (rows > 0 && cols > 0) {
[ # # ][ # # ]
150 [ + + ][ # # ]: 255730 : data = new nr_complex_t[rows * cols];
151 : 13130 : memcpy (data, m.data, sizeof (nr_complex_t) * rows * cols);
152 : : }
153 : 13130 : }
154 : :
155 : : /*!\brief Assignment operator
156 : :
157 : : The assignment copy constructor creates a new instance based on the
158 : : given matrix object.
159 : :
160 : : \param[in] m object to copy
161 : : \return assigned object
162 : : \note m = m is safe
163 : : */
164 : 13273 : const matrix& matrix::operator=(const matrix & m) {
165 [ + - ]: 13273 : if (&m != this) {
166 : 13273 : rows = m.rows;
167 : 13273 : cols = m.cols;
168 [ - + ]: 13273 : if (data) {
169 [ # # ]: 0 : delete[] data;
170 : 0 : data = NULL;
171 : : }
172 [ + - ][ + - ]: 13273 : if (rows > 0 && cols > 0) {
173 [ + + ]: 148912 : data = new nr_complex_t[rows * cols];
174 : 13273 : memcpy (data, m.data, sizeof (nr_complex_t) * rows * cols);
175 : : }
176 : : }
177 : 13273 : return *this;
178 : : }
179 : :
180 : : /*!\bried Destructor
181 : :
182 : : Destructor deletes a matrix object.
183 : : */
184 : 48266 : matrix::~matrix () {
185 [ + - ][ + - ]: 48266 : if (data) delete[] data;
[ # # ][ # # ]
186 : 48266 : }
187 : :
188 : : /*!\brief Returns the matrix element at the given row and column.
189 : : \param[in] r row number
190 : : \param[in] c column number
191 : : \todo Why not inline and synonymous of ()
192 : : \todo c and r const
193 : : */
194 : 1394550 : nr_complex_t matrix::get (int r, int c) {
195 : 1394550 : return data[r * cols + c];
196 : : }
197 : :
198 : : /*!\brief Sets the matrix element at the given row and column.
199 : : \param[in] r row number
200 : : \param[in] c column number
201 : : \param[in] z complex number to assign
202 : : \todo Why not inline and synonymous of ()
203 : : \todo r c and z const
204 : : */
205 : 248237 : void matrix::set (int r, int c, nr_complex_t z) {
206 : 248237 : data[r * cols + c] = z;
207 : 248237 : }
208 : :
209 : : #ifdef DEBUG
210 : : /*!\brief Debug function: Prints the matrix object */
211 : 0 : void matrix::print (void) {
212 [ # # ]: 0 : for (int r = 0; r < rows; r++) {
213 [ # # ]: 0 : for (int c = 0; c < cols; c++) {
214 : : fprintf (stderr, "%+.2e,%+.2e ", (double) real (get (r, c)),
215 [ # # ]: 0 : (double) imag (get (r, c)));
216 : : }
217 : 0 : fprintf (stderr, "\n");
218 : : }
219 : 0 : }
220 : : #endif /* DEBUG */
221 : :
222 : : /*!\brief Matrix addition.
223 : : \param[a] first matrix
224 : : \param[b] second matrix
225 : : \note assert same size
226 : : \todo a and b are const
227 : : */
228 : 1210 : matrix operator + (matrix a, matrix b) {
229 [ + - ][ - + ]: 1210 : assert (a.getRows () == b.getRows () && a.getCols () == b.getCols ());
[ - + ]
230 : :
231 : 1210 : matrix res (a.getRows (), a.getCols ());
232 [ + + ]: 6110 : for (int r = 0; r < a.getRows (); r++)
233 [ + + ]: 34700 : for (int c = 0; c < a.getCols (); c++)
234 [ + - ]: 29800 : res.set (r, c, a.get (r, c) + b.get (r, c));
235 : 1210 : return res;
236 : : }
237 : :
238 : : /*!\brief Intrinsic matrix addition.
239 : : \param[in] a matrix to add
240 : : \note assert same size
241 : : \todo a is const
242 : : */
243 : 0 : matrix matrix::operator += (matrix a) {
244 [ # # ][ # # ]: 0 : assert (a.getRows () == rows && a.getCols () == cols);
[ # # ]
245 : :
246 : : int r, c, i;
247 [ # # ]: 0 : for (i = 0, r = 0; r < a.getRows (); r++)
248 [ # # ]: 0 : for (c = 0; c < a.getCols (); c++, i++)
249 : 0 : data[i] += a.get (r, c);
250 : 0 : return *this;
251 : : }
252 : :
253 : : /*!\brief Matrix subtraction.
254 : : \param[a] first matrix
255 : : \param[b] second matrix
256 : : \note assert same size
257 : : \todo a and b are const
258 : : */
259 : 1070 : matrix operator - (matrix a, matrix b) {
260 [ + - ][ - + ]: 1070 : assert (a.getRows () == b.getRows () && a.getCols () == b.getCols ());
[ - + ]
261 : :
262 : 1070 : matrix res (a.getRows (), a.getCols ());
263 [ + + ]: 4570 : for (int r = 0; r < a.getRows (); r++)
264 [ + + ]: 19300 : for (int c = 0; c < a.getCols (); c++)
265 [ + - ]: 15800 : res.set (r, c, a.get (r, c) - b.get (r, c));
266 : 1070 : return res;
267 : : }
268 : :
269 : : /*!\brief Unary minus. */
270 : 0 : matrix matrix::operator - () {
271 : 0 : matrix res (getRows (), getCols ());
272 : : int r, c, i;
273 [ # # ]: 0 : for (i = 0, r = 0; r < getRows (); r++)
274 [ # # ]: 0 : for (c = 0; c < getCols (); c++, i++)
275 : 0 : res.set (r, c, -data[i]);
276 : 0 : return res;
277 : : }
278 : :
279 : : /*!\brief Intrinsic matrix subtraction.
280 : : \param[in] a matrix to substract
281 : : \note assert same size
282 : : */
283 : 0 : matrix matrix::operator -= (matrix a) {
284 [ # # ][ # # ]: 0 : assert (a.getRows () == rows && a.getCols () == cols);
[ # # ]
285 : : int r, c, i;
286 [ # # ]: 0 : for (i = 0, r = 0; r < a.getRows (); r++)
287 [ # # ]: 0 : for (c = 0; c < a.getCols (); c++, i++)
288 : 0 : data[i] -= a.get (r, c);
289 : 0 : return *this;
290 : : }
291 : :
292 : : /*!\brief Matrix scaling complex version
293 : : \param[in] a matrix to scale
294 : : \param[in] z scaling complex
295 : : \return Scaled matrix
296 : : \todo Why not a and z const
297 : : */
298 : 0 : matrix operator * (matrix a, nr_complex_t z) {
299 : 0 : matrix res (a.getRows (), a.getCols ());
300 [ # # ]: 0 : for (int r = 0; r < a.getRows (); r++)
301 [ # # ]: 0 : for (int c = 0; c < a.getCols (); c++)
302 [ # # ]: 0 : res.set (r, c, a.get (r, c) * z);
303 : 0 : return res;
304 : : }
305 : :
306 : : /*!\brief Matrix scaling complex version (different order)
307 : : \param[in] a matrix to scale
308 : : \param[in] z scaling complex
309 : : \return Scaled matrix
310 : : \todo Why not a and z const
311 : : \todo Why not inline
312 : : */
313 : 0 : matrix operator * (nr_complex_t z, matrix a) {
314 [ # # ]: 0 : return a * z;
315 : : }
316 : :
317 : : /*!\brief Matrix scaling complex version
318 : : \param[in] a matrix to scale
319 : : \param[in] d scaling real
320 : : \return Scaled matrix
321 : : \todo Why not d and a const
322 : : */
323 : 70 : matrix operator * (matrix a, nr_double_t d) {
324 : 70 : matrix res (a.getRows (), a.getCols ());
325 [ + + ]: 770 : for (int r = 0; r < a.getRows (); r++)
326 [ + + ]: 7700 : for (int c = 0; c < a.getCols (); c++)
327 : 7000 : res.set (r, c, a.get (r, c) * d);
328 : 70 : return res;
329 : : }
330 : :
331 : : /*!\brief Matrix scaling real version (different order)
332 : : \param[in] a matrix to scale
333 : : \param[in] d scaling real
334 : : \return Scaled matrix
335 : : \todo Why not inline?
336 : : \todo Why not d and a const
337 : : */
338 : 0 : matrix operator * (nr_double_t d, matrix a) {
339 [ # # ]: 0 : return a * d;
340 : : }
341 : :
342 : : /*!\brief Matrix scaling division by complex version
343 : : \param[in] a matrix to scale
344 : : \param[in] z scaling complex
345 : : \return Scaled matrix
346 : : \todo Why not a and z const
347 : : */
348 : 0 : matrix operator / (matrix a, nr_complex_t z) {
349 : 0 : matrix res (a.getRows (), a.getCols ());
350 [ # # ]: 0 : for (int r = 0; r < a.getRows (); r++)
351 [ # # ]: 0 : for (int c = 0; c < a.getCols (); c++)
352 [ # # ]: 0 : res.set (r, c, a.get (r, c) / z);
353 : 0 : return res;
354 : : }
355 : :
356 : : /*!\brief Matrix scaling division by real version
357 : : \param[in] a matrix to scale
358 : : \param[in] d scaling real
359 : : \return Scaled matrix
360 : : \todo Why not a and d const
361 : : */
362 : 70 : matrix operator / (matrix a, nr_double_t d) {
363 : 70 : matrix res (a.getRows (), a.getCols ());
364 [ + + ]: 770 : for (int r = 0; r < a.getRows (); r++)
365 [ + + ]: 7700 : for (int c = 0; c < a.getCols (); c++)
366 : 7000 : res.set (r, c, a.get (r, c) / d);
367 : 70 : return res;
368 : : }
369 : :
370 : : /*! Matrix multiplication.
371 : :
372 : : Dumb and not optimized matrix multiplication
373 : : \param[a] first matrix
374 : : \param[b] second matrix
375 : : \note assert compatibility
376 : : \todo a and b are const
377 : : */
378 : 4290 : matrix operator * (matrix a, matrix b) {
379 [ - + ]: 4290 : assert (a.getCols () == b.getRows ());
380 : :
381 : 4290 : int r, c, i, n = a.getCols ();
382 : 4290 : nr_complex_t z;
383 [ + - ]: 4290 : matrix res (a.getRows (), b.getCols ());
384 [ + + ]: 20790 : for (r = 0; r < a.getRows (); r++) {
385 [ + + ]: 104700 : for (c = 0; c < b.getCols (); c++) {
386 [ + - ][ + + ]: 720600 : for (i = 0, z = 0; i < n; i++) z += a.get (r, i) * b.get (i, c);
387 : 88200 : res.set (r, c, z);
388 : : }
389 : : }
390 : 4290 : return res;
391 : : }
392 : :
393 : : /*!\brief Complex scalar addition.
394 : : \param[in] a matrix
395 : : \param[in] z complex to add
396 : : \todo Move near other +
397 : : \todo a and z are const
398 : : */
399 : 0 : matrix operator + (matrix a, nr_complex_t z) {
400 : 0 : matrix res (a.getRows (), a.getCols ());
401 [ # # ]: 0 : for (int r = 0; r < a.getRows (); r++)
402 [ # # ]: 0 : for (int c = 0; c < a.getCols (); c++)
403 [ # # ]: 0 : res.set (r, c, a.get (r, c) + z);
404 : 0 : return res;
405 : : }
406 : :
407 : : /*!\brief Complex scalar addition different order.
408 : : \param[in] a matrix
409 : : \param[in] z complex to add
410 : : \todo Move near other +
411 : : \todo a and z are const
412 : : \todo Why not inline
413 : : */
414 : 0 : matrix operator + (nr_complex_t z, matrix a) {
415 [ # # ]: 0 : return a + z;
416 : : }
417 : :
418 : : /*!\brief Real scalar addition.
419 : : \param[in] a matrix
420 : : \param[in] d real to add
421 : : \todo Move near other +
422 : : \todo a and d are const
423 : : */
424 : 0 : matrix operator + (matrix a, nr_double_t d) {
425 : 0 : matrix res (a.getRows (), a.getCols ());
426 [ # # ]: 0 : for (int r = 0; r < a.getRows (); r++)
427 [ # # ]: 0 : for (int c = 0; c < a.getCols (); c++)
428 : 0 : res.set (r, c, a.get (r, c) + d);
429 : 0 : return res;
430 : : }
431 : :
432 : : /*!\brief Real scalar addition different order.
433 : : \param[in] a matrix
434 : : \param[in] d real to add
435 : : \todo Move near other +
436 : : \todo a and d are const
437 : : \todo Why not inline
438 : : */
439 : 0 : matrix operator + (nr_double_t d, matrix a) {
440 [ # # ]: 0 : return a + d;
441 : : }
442 : :
443 : : /*!\brief Complex scalar substraction
444 : : \param[in] a matrix
445 : : \param[in] z complex to add
446 : : \todo Move near other +
447 : : \todo a and z are const
448 : : \todo Why not inline
449 : : */
450 : 0 : matrix operator - (matrix a, nr_complex_t z) {
451 [ # # ]: 0 : return -z + a;
452 : : }
453 : :
454 : : /*!\brief Complex scalar substraction different order
455 : : \param[in] a matrix
456 : : \param[in] z complex to add
457 : : \todo Move near other +
458 : : \todo a and z are const
459 : : \todo Why not inline
460 : : */
461 : 0 : matrix operator - (nr_complex_t z, matrix a) {
462 [ # # ]: 0 : return -a + z;
463 : : }
464 : :
465 : : /*!\brief Real scalar substraction
466 : : \param[in] a matrix
467 : : \param[in] z real to add
468 : : \todo Move near other +
469 : : \todo a and z are const
470 : : \todo Why not inline
471 : : */
472 : 0 : matrix operator - (matrix a, nr_double_t d) {
473 [ # # ]: 0 : return -d + a;
474 : : }
475 : :
476 : : /*!\brief Real scalar substraction different order
477 : : \param[in] a matrix
478 : : \param[in] z real to add
479 : : \todo Move near other +
480 : : \todo a and z are const
481 : : \todo Why not inline
482 : : */
483 : 0 : matrix operator - (nr_double_t d, matrix a) {
484 [ # # ]: 0 : return -a + d;
485 : : }
486 : :
487 : : /*!\brief Matrix transposition
488 : : \param[in] a Matrix to transpose
489 : : \todo add transpose in place
490 : : \todo a is const
491 : : */
492 : 70 : matrix transpose (matrix a) {
493 : 70 : matrix res (a.getCols (), a.getRows ());
494 [ + + ]: 770 : for (int r = 0; r < a.getRows (); r++)
495 [ + + ]: 7700 : for (int c = 0; c < a.getCols (); c++)
496 : 7000 : res.set (c, r, a.get (r, c));
497 : 70 : return res;
498 : : }
499 : :
500 : : /*!\brief Conjugate complex matrix.
501 : : \param[in] a Matrix to conjugate
502 : : \todo add conj in place
503 : : \todo a is const
504 : : */
505 : 70 : matrix conj (matrix a) {
506 : 70 : matrix res (a.getRows (), a.getCols ());
507 [ + + ]: 770 : for (int r = 0; r < a.getRows (); r++)
508 [ + + ]: 7700 : for (int c = 0; c < a.getCols (); c++)
509 : 7000 : res.set (r, c, conj (a.get (r, c)));
510 : 70 : return res;
511 : : }
512 : :
513 : : /*!\brief adjoint matrix
514 : :
515 : : The function returns the adjoint complex matrix. This is also
516 : : called the adjugate or transpose conjugate.
517 : : \param[in] a Matrix to transpose
518 : : \todo add adjoint in place
519 : : \todo Do not lazy and avoid conj and transpose copy
520 : : \todo a is const
521 : : */
522 : 70 : matrix adjoint (matrix a) {
523 [ + - ][ + - ]: 70 : return transpose (conj (a));
524 : : }
525 : :
526 : : /*!\brief Computes magnitude of each matrix element.
527 : : \param[in] a matrix
528 : : \todo add abs in place
529 : : \todo a is const
530 : : */
531 : 0 : matrix abs (matrix a) {
532 : 0 : matrix res (a.getRows (), a.getCols ());
533 [ # # ]: 0 : for (int r = 0; r < a.getRows (); r++)
534 [ # # ]: 0 : for (int c = 0; c < a.getCols (); c++)
535 : 0 : res.set (r, c, abs (a.get (r, c)));
536 : 0 : return res;
537 : : }
538 : :
539 : : /*!\brief Computes magnitude in dB of each matrix element.
540 : : \param[in] a matrix
541 : : */
542 : 0 : matrix dB (matrix a) {
543 : 0 : matrix res (a.getRows (), a.getCols ());
544 [ # # ]: 0 : for (int r = 0; r < a.getRows (); r++)
545 [ # # ]: 0 : for (int c = 0; c < a.getCols (); c++)
546 [ # # ]: 0 : res.set (r, c, dB (a.get (r, c)));
547 : 0 : return res;
548 : : }
549 : :
550 : : /*!\brief Computes the argument of each matrix element.
551 : : \param[in] a matrix
552 : : \todo add arg in place
553 : : \todo a is const
554 : : */
555 : 0 : matrix arg (matrix a) {
556 : 0 : matrix res (a.getRows (), a.getCols ());
557 [ # # ]: 0 : for (int r = 0; r < a.getRows (); r++)
558 [ # # ]: 0 : for (int c = 0; c < a.getCols (); c++)
559 : 0 : res.set (r, c, arg (a.get (r, c)));
560 : 0 : return res;
561 : : }
562 : :
563 : : /*!\brief Real part matrix.
564 : : \param[in] a matrix
565 : : \todo add real in place
566 : : \todo a is const
567 : : */
568 : 0 : matrix real (matrix a) {
569 : 0 : matrix res (a.getRows (), a.getCols ());
570 [ # # ]: 0 : for (int r = 0; r < a.getRows (); r++)
571 [ # # ]: 0 : for (int c = 0; c < a.getCols (); c++)
572 : 0 : res.set (r, c, real (a.get (r, c)));
573 : 0 : return res;
574 : : }
575 : :
576 : : /*!\brief Imaginary part matrix.
577 : : \param[in] a matrix
578 : : \todo add imag in place
579 : : \todo a is const
580 : : */
581 : 0 : matrix imag (matrix a) {
582 : 0 : matrix res (a.getRows (), a.getCols ());
583 [ # # ]: 0 : for (int r = 0; r < a.getRows (); r++)
584 [ # # ]: 0 : for (int c = 0; c < a.getCols (); c++)
585 : 0 : res.set (r, c, imag (a.get (r, c)));
586 : 0 : return res;
587 : : }
588 : :
589 : : /*!\brief Multiply a matrix by itself
590 : : \param[in] a matrix
591 : : */
592 : 0 : matrix sqr (matrix a) {
593 [ # # ][ # # ]: 0 : return a * a;
594 : : }
595 : :
596 : : /*!\brief Create identity matrix with specified number of rows and columns.
597 : : \param[in] rs row number
598 : : \param[in] cs column number
599 : : \todo Avoid res.get*
600 : : \todo Use memset
601 : : \todo rs, cs are const
602 : : */
603 : 3280 : matrix eye (int rs, int cs) {
604 : 3280 : matrix res (rs, cs);
605 [ + + ]: 14480 : for (int r = 0; r < res.getRows (); r++)
606 [ + + ]: 65600 : for (int c = 0; c < res.getCols (); c++)
607 [ + + ]: 54400 : if (r == c) res.set (r, c, 1);
608 : 3280 : return res;
609 : : }
610 : :
611 : : /*!\brief Create a square identity matrix
612 : : \param[in] s row or column number of square matrix
613 : : \todo Do not by lazy and implement it
614 : : \todo s is const
615 : : */
616 : 3280 : matrix eye (int s) {
617 : 3280 : return eye (s, s);
618 : : }
619 : :
620 : : /*!\brief Create a diagonal matrix from a vector
621 : : \param[in] diag vector to write on the diagonal
622 : : \todo diag is const
623 : : */
624 : 2140 : matrix diagonal (qucs::vector diag) {
625 : 2140 : int size = diag.getSize ();
626 : 2140 : matrix res (size);
627 [ + + ]: 9140 : for (int i = 0; i < size; i++) res (i, i) = diag (i);
628 : 2140 : return res;
629 : : }
630 : :
631 : : // Compute n-th power of the given matrix.
632 : 0 : matrix pow (matrix a, int n) {
633 : 0 : matrix res;
634 [ # # ]: 0 : if (n == 0) {
635 [ # # ][ # # ]: 0 : res = eye (a.getRows (), a.getCols ());
636 : : }
637 : : else {
638 [ # # ][ # # ]: 0 : res = a = n < 0 ? inverse (a) : a;
[ # # ][ # # ]
[ # # ][ # # ]
[ # # ][ # # ]
639 [ # # ]: 0 : for (int i = 1; i < std::abs (n); i++)
640 [ # # ][ # # ]: 0 : res = res * a;
[ # # ][ # # ]
641 : : }
642 : 0 : return res;
643 : : }
644 : :
645 : : /*!\brief Computes the complex cofactor of the given determinant
646 : :
647 : : The cofactor is the determinant obtained by deleting the row and column
648 : : of a given element of a matrix or determinant.
649 : : The cofactor is preceded by a + or - sign depending of the sign
650 : : of \f$(-1)^(u+v)\f$
651 : :
652 : : \bug This algortihm is recursive! Stack overfull!
653 : : \todo ((u + v) & 1) is cryptic use (u + v)% 2
654 : : \todo #ifdef 0
655 : : \todo static?
656 : : */
657 : 0 : nr_complex_t cofactor (matrix a, int u, int v) {
658 [ # # ]: 0 : matrix res (a.getRows () - 1, a.getCols () - 1);
659 : : int r, c, ra, ca;
660 [ # # ]: 0 : for (ra = r = 0; r < res.getRows (); r++, ra++) {
661 [ # # ]: 0 : if (ra == u) ra++;
662 [ # # ]: 0 : for (ca = c = 0; c < res.getCols (); c++, ca++) {
663 [ # # ]: 0 : if (ca == v) ca++;
664 : 0 : res.set (r, c, a.get (ra, ca));
665 : : }
666 : : }
667 [ # # ][ # # ]: 0 : nr_complex_t z = detLaplace (res);
668 [ # # ]: 0 : return ((u + v) & 1) ? -z : z;
669 : : }
670 : :
671 : : /*!\brief Compute determinant of the given matrix using Laplace expansion.
672 : :
673 : : The Laplace expansion of the determinant of
674 : : an n by n square matrix a expresses
675 : : the determinant of a as a sum of n determinants of (n-1) by (n-1)
676 : : sub-matrices of a. There are 2n such expressions, one for each row
677 : : and column of a.
678 : :
679 : : See Wikipedia http://en.wikipedia.org/wiki/Laplace_expansion
680 : : \param[in] a matrix
681 : : \bug This algortihm is recursive! Stack overfull!
682 : : \note assert square matrix
683 : : \todo #ifdef 0
684 : : \todo static ?
685 : : */
686 : 0 : nr_complex_t detLaplace (matrix a) {
687 [ # # ]: 0 : assert (a.getRows () == a.getCols ());
688 : 0 : int s = a.getRows ();
689 : 0 : nr_complex_t res = 0;
690 [ # # ]: 0 : if (s > 1) {
691 : : /* always use the first row for sub-determinant, but you should
692 : : use the row or column with most zeros in it */
693 : 0 : int r = 0;
694 [ # # ]: 0 : for (int i = 0; i < s; i++) {
695 [ # # ][ # # ]: 0 : res += a.get (r, i) * cofactor (a, r, i);
[ # # ]
696 : : }
697 : 0 : return res;
698 : : }
699 : : /* 1 by 1 matrix */
700 [ # # ]: 0 : else if (s == 1) {
701 : 0 : return a (0, 0);
702 : : }
703 : : /* 0 by 0 matrix */
704 : 0 : return 1;
705 : : }
706 : :
707 : : /*!\brief Compute determinant Gaussian algorithm
708 : :
709 : : Compute determinant of the given matrix using the Gaussian
710 : : algorithm. This means to triangulate the matrix and multiply all
711 : : the diagonal elements.
712 : : \param[in] a matrix
713 : : \note assert square matrix
714 : : \todo static ?
715 : : \todo a const?
716 : : */
717 : 0 : nr_complex_t detGauss (matrix a) {
718 [ # # ]: 0 : assert (a.getRows () == a.getCols ());
719 : : nr_double_t MaxPivot;
720 : 0 : nr_complex_t f, res;
721 : 0 : matrix b;
722 : 0 : int i, c, r, pivot, n = a.getCols ();
723 : :
724 : : // return special matrix cases
725 [ # # ]: 0 : if (n == 0) return 1;
726 [ # # ]: 0 : if (n == 1) return a (0, 0);
727 : :
728 : : // make copy of original matrix
729 [ # # ][ # # ]: 0 : b = matrix (a);
730 : :
731 : : // triangulate the matrix
732 [ # # ]: 0 : for (res = 1, i = 0; i < n; i++) {
733 : : // find maximum column value for pivoting
734 [ # # ]: 0 : for (MaxPivot = 0, pivot = r = i; r < n; r++) {
735 [ # # ]: 0 : if (abs (b.get (r, i)) > MaxPivot) {
736 : 0 : MaxPivot = abs (b.get (r, i));
737 : 0 : pivot = r;
738 : : }
739 : : }
740 : : // exchange rows if necessary
741 [ # # ]: 0 : assert (MaxPivot != 0);
742 [ # # ][ # # ]: 0 : if (i != pivot) { b.exchangeRows (i, pivot); res = -res; }
743 : : // compute new rows and columns
744 [ # # ]: 0 : for (r = i + 1; r < n; r++) {
745 [ # # ]: 0 : f = b.get (r, i) / b.get (i, i);
746 [ # # ]: 0 : for (c = i + 1; c < n; c++) {
747 [ # # ][ # # ]: 0 : b.set (r, c, b.get (r, c) - f * b.get (i, c));
748 : : }
749 : : }
750 : : }
751 : :
752 : : // now compute determinant by multiplying diagonal
753 [ # # ]: 0 : for (i = 0; i < n; i++) res *= b.get (i, i);
754 : 0 : return res;
755 : : }
756 : :
757 : : /*!\brief Compute determinant of the given matrix.
758 : : \param[in] a matrix
759 : : \return Complex determinant
760 : : \todo a const?
761 : : */
762 : 0 : nr_complex_t det (matrix a) {
763 : : #if 0
764 : : return detLaplace (a);
765 : : #else
766 [ # # ]: 0 : return detGauss (a);
767 : : #endif
768 : : }
769 : :
770 : : /*!\brief Compute inverse matrix using Laplace expansion
771 : :
772 : : Compute inverse matrix of the given matrix using Laplace expansion.
773 : : \param[in] a matrix to invert
774 : : \todo Static?
775 : : \bug recursive! Stack overflow
776 : : \todo a const?
777 : : \todo #ifdef 0
778 : : */
779 : 0 : matrix inverseLaplace (matrix a) {
780 [ # # ]: 0 : matrix res (a.getRows (), a.getCols ());
781 [ # # ][ # # ]: 0 : nr_complex_t d = detLaplace (a);
782 [ # # ]: 0 : assert (abs (d) != 0); // singular matrix
783 [ # # ]: 0 : for (int r = 0; r < a.getRows (); r++)
784 [ # # ]: 0 : for (int c = 0; c < a.getCols (); c++)
785 [ # # ][ # # ]: 0 : res.set (r, c, cofactor (a, c, r) / d);
[ # # ]
786 : 0 : return res;
787 : : }
788 : :
789 : : /*!\brief Compute inverse matrix using Gauss-Jordan elimination
790 : :
791 : : Compute inverse matrix of the given matrix by Gauss-Jordan
792 : : elimination.
793 : : \todo a const?
794 : : \todo static?
795 : : \note assert non singulat matix
796 : : \param[in] a matrix to invert
797 : : */
798 : 2140 : matrix inverseGaussJordan (matrix a) {
799 : : nr_double_t MaxPivot;
800 : 2140 : nr_complex_t f;
801 : 2140 : matrix b, e;
802 : 2140 : int i, c, r, pivot, n = a.getCols ();
803 : :
804 : : // create temporary matrix and the result matrix
805 [ + - ][ + - ]: 2140 : b = matrix (a);
806 [ + - ][ + - ]: 2140 : e = eye (n);
807 : :
808 : : // create the eye matrix in 'b' and the result in 'e'
809 [ + + ]: 9140 : for (i = 0; i < n; i++) {
810 : : // find maximum column value for pivoting
811 [ + + ]: 26300 : for (MaxPivot = 0, pivot = r = i; r < n; r++) {
812 [ + + ]: 19300 : if (abs (b (r, i)) > MaxPivot) {
813 : 7700 : MaxPivot = abs (b (r, i));
814 : 7700 : pivot = r;
815 : : }
816 : : }
817 : : // exchange rows if necessary
818 [ - + ]: 7000 : assert (MaxPivot != 0); // singular matrix
819 [ + + ]: 7000 : if (i != pivot) {
820 [ + - ]: 500 : b.exchangeRows (i, pivot);
821 [ + - ]: 500 : e.exchangeRows (i, pivot);
822 : : }
823 : :
824 : : // compute current row
825 [ + + ]: 38600 : for (f = b (i, i), c = 0; c < n; c++) {
826 : 31600 : b (i, c) /= f;
827 : 31600 : e (i, c) /= f;
828 : : }
829 : :
830 : : // compute new rows and columns
831 [ + + ]: 38600 : for (r = 0; r < n; r++) {
832 [ + + ]: 31600 : if (r != i) {
833 [ + + ]: 193800 : for (f = b (r, i), c = 0; c < n; c++) {
834 [ + - ]: 169200 : b (r, c) -= f * b (i, c);
835 [ + - ]: 169200 : e (r, c) -= f * e (i, c);
836 : : }
837 : : }
838 : : }
839 : : }
840 : 2140 : return e;
841 : : }
842 : :
843 : : /*!\brief Compute inverse matrix
844 : : \param[in] a matrix to invert
845 : : \todo a is const
846 : : */
847 : 2140 : matrix inverse (matrix a) {
848 : : #if 0
849 : : return inverseLaplace (a);
850 : : #else
851 [ + - ]: 2140 : return inverseGaussJordan (a);
852 : : #endif
853 : : }
854 : :
855 : : /*!\brief S params to S params
856 : :
857 : : Convert scattering parameters with the reference impedance 'zref'
858 : : to scattering parameters with the reference impedance 'z0'.
859 : :
860 : : Detail are given in [1], under equation (32)
861 : :
862 : : New scatering matrix \f$S'\f$ is:
863 : : \f[
864 : : S'=A^{-1}(S-\Gamma^+)(I-\Gamma S)^{-1}A^+
865 : : \f]
866 : : Where x^+ is the adjoint (or complex tranposate) of x,
867 : : I the identity matrix and \f$A\f$ is diagonal the matrix such as:
868 : : \f$ \Gamma_i= r_i \f$ and \f$A\f$ the diagonal matrix such
869 : : as:
870 : : \f[
871 : : A_i = \frac{(1-r_i^*)\sqrt{|1-r_ir_i^*|}}{|1-r_i|}
872 : : \f]
873 : : Where \f$x*\f$ is the complex conjugate of \f$x\f$
874 : : and \f$r_i\f$ is wave reflexion coefficient of \f$Z_i'\f$ with respect
875 : : to \f$Z_i^*\f$ (where \f$Z_i'\f$ is the new impedance and
876 : : \f$Z_i\f$ is the old impedance), ie:
877 : : \f[
878 : : r_i = \frac{Z_i'-Z_i}{Z_i'-Z_i^*}
879 : : \f]
880 : :
881 : : \param[in] s original S matrix
882 : : \param[in] zref original reference impedance
883 : : \param[in] z0 new reference impedance
884 : : \bug This formula is valid only for real z!
885 : : \todo Correct documentation about standing waves [1-4]
886 : : \todo Implement Speciale implementation [2-3] if applicable
887 : : \return Renormalized scattering matrix
888 : : \todo s, zref and z0 const
889 : : */
890 : 0 : matrix stos (matrix s, qucs::vector zref, qucs::vector z0) {
891 : 0 : int d = s.getRows ();
892 : 0 : matrix e, r, a;
893 : :
894 [ # # ][ # # ]: 0 : assert (d == s.getCols () && d == z0.getSize () && d == zref.getSize ());
[ # # ][ # # ]
[ # # ][ # # ]
895 : :
896 [ # # ][ # # ]: 0 : e = eye (d);
897 [ # # ][ # # ]: 0 : r = diagonal ((z0 - zref) / (z0 + zref));
[ # # ][ # # ]
[ # # ][ # # ]
[ # # ][ # # ]
[ # # ][ # # ]
[ # # ][ # # ]
[ # # ][ # # ]
[ # # ][ # # ]
898 [ # # ][ # # ]: 0 : a = diagonal (sqrt (z0 / zref) / (z0 + zref));
[ # # ][ # # ]
[ # # ][ # # ]
[ # # ][ # # ]
[ # # ][ # # ]
[ # # ][ # # ]
[ # # ][ # # ]
[ # # ][ # # ]
[ # # ][ # # ]
899 [ # # ][ # # ]: 0 : return inverse (a) * (s - r) * inverse (e - r * s) * a;
[ # # ][ # # ]
[ # # ][ # # ]
[ # # ][ # # ]
[ # # ][ # # ]
[ # # ][ # # ]
[ # # ][ # # ]
[ # # ]
900 : : }
901 : :
902 : : /*!\brief S renormalization with all part identic
903 : : \param[in] s original S matrix
904 : : \param[in] zref original reference impedance
905 : : \param[in] z0 new reference impedance
906 : : \todo Why not inline
907 : : \return Renormalized scattering matrix
908 : : \todo s, zref and z0 const
909 : : */
910 : 0 : matrix stos (matrix s, nr_complex_t zref, nr_complex_t z0) {
911 : 0 : int d = s.getRows ();
912 [ # # ][ # # ]: 0 : return stos (s, qucs::vector (d, zref), qucs::vector (d, z0));
[ # # ][ # # ]
913 : : }
914 : :
915 : : /*!\brief S renormalization with all part identic and real
916 : : \param[in] s original S matrix
917 : : \param[in] zref original reference impedance
918 : : \param[in] z0 new reference impedance
919 : : \todo Why not inline
920 : : \return Renormalized scattering matrix
921 : : \todo s, zref and z0 const
922 : : */
923 : 0 : matrix stos (matrix s, nr_double_t zref, nr_double_t z0) {
924 [ # # ][ # # ]: 0 : return stos (s, nr_complex_t (zref, 0), nr_complex_t (z0, 0));
925 : : }
926 : :
927 : : /*!\brief S renormalization (variation)
928 : : \param[in] s original S matrix
929 : : \param[in] zref original reference impedance
930 : : \param[in] z0 new reference impedance
931 : : \todo Why not inline
932 : : \return Renormalized scattering matrix
933 : : \todo s, zref and z0 const
934 : : */
935 : 0 : matrix stos (matrix s, qucs::vector zref, nr_complex_t z0) {
936 [ # # ][ # # ]: 0 : return stos (s, zref, qucs::vector (zref.getSize (), z0));
[ # # ][ # # ]
937 : : }
938 : :
939 : : /*!\brief S renormalization (variation)
940 : : \param[in] s original S matrix
941 : : \param[in] zref original reference impedance
942 : : \param[in] z0 new reference impedance
943 : : \todo Why not inline
944 : : \todo s, zref and z0 const
945 : : \return Renormalized scattering matrix
946 : : */
947 : 0 : matrix stos (matrix s, nr_complex_t zref, qucs::vector z0) {
948 [ # # ][ # # ]: 0 : return stos (s, qucs::vector (z0.getSize (), zref), z0);
[ # # ][ # # ]
[ # # ]
949 : : }
950 : :
951 : : /*!\brief Scattering parameters to impedance matrix
952 : :
953 : : Convert scattering parameters to impedance matrix.
954 : : According to [1] eq (19):
955 : : \f[
956 : : Z=F^{-1} (I- S)^{-1} (SG + G^+) F
957 : : \f]
958 : : Where \f$S\f$ is the scattering matrix, \f$x^+\f$
959 : : is the adjoint of x, I the identity matrix. The matrix
960 : : F and G are diagonal matrix defined by:
961 : : \f{align*}
962 : : F_i&=\frac{1}{2\sqrt{\Re\text{e}\; Z_i}} \\
963 : : G_i&=Z_i
964 : : \f}
965 : : \param[in] s Scattering matrix
966 : : \param[in] z0 Normalisation impedance
967 : : \note We could safely drop the \f$1/2\f$ in \f$F\f$ because we compute
968 : : \f$FXF^{-1}\f$ and therefore \f$1/2\f$ will simplify.
969 : : \bug not correct if zref is complex
970 : : \todo s, z0 const
971 : : \return Impedance matrix
972 : : */
973 : 0 : matrix stoz (matrix s, qucs::vector z0) {
974 : 0 : int d = s.getRows ();
975 : 0 : matrix e, zref, gref;
976 : :
977 [ # # ][ # # ]: 0 : assert (d == s.getCols () && d == z0.getSize ());
[ # # ][ # # ]
978 : :
979 [ # # ][ # # ]: 0 : e = eye (d);
980 [ # # ][ # # ]: 0 : zref = diagonal (z0);
[ # # ][ # # ]
981 [ # # ][ # # ]: 0 : gref = diagonal (sqrt (real (1 / z0)));
[ # # ][ # # ]
[ # # ][ # # ]
[ # # ][ # # ]
[ # # ][ # # ]
982 [ # # ][ # # ]: 0 : return inverse (gref) * inverse (e - s) * (s * zref + zref) * gref;
[ # # ][ # # ]
[ # # ][ # # ]
[ # # ][ # # ]
[ # # ][ # # ]
[ # # ][ # # ]
[ # # ][ # # ]
[ # # ]
983 : : }
984 : :
985 : : /*!\brief Scattering parameters to impedance matrix identic case
986 : : \param[in] s Scattering matrix
987 : : \param[in] z0 Normalisation impedance
988 : : \return Impedance matrix
989 : : \todo Why not inline?
990 : : \todo s and z0 const?
991 : : */
992 : 0 : matrix stoz (matrix s, nr_complex_t z0) {
993 [ # # ][ # # ]: 0 : return stoz (s, qucs::vector (s.getRows (), z0));
994 : : }
995 : :
996 : : /*!\brief Convert impedance matrix scattering parameters.
997 : :
998 : : Convert scattering parameters to impedance matrix.
999 : : According to [1] eq (18):
1000 : : \f[
1001 : : S=F(Z-G^+)(Z+G)^{-1} F^{-1}
1002 : : \f]
1003 : : Where \f$Z\f$ is the scattering matrix, \f$x^+\f$
1004 : : is the adjoint of x, I the identity matrix. The matrix
1005 : : F and G are diagonal matrix defined by:
1006 : : \f{align*}
1007 : : F_i&=\frac{1}{2\sqrt{\Re\text{e}\; Z_i}} \\
1008 : : G_i&=Z_i
1009 : : \f}
1010 : : \param[in] Z Impedance matrix
1011 : : \param[in] z0 Normalisation impedance
1012 : : \return Scattering matrix
1013 : : \note We could safely drop the \f$1/2\f$ in \f$F\f$ because we compute
1014 : : \f$FXF^{-1}\f$ and therefore \f$1/2\f$ will simplify.
1015 : : \bug not correct if zref is complex
1016 : : \todo z and z0 const?
1017 : : */
1018 : 600 : matrix ztos (matrix z, qucs::vector z0) {
1019 : 600 : int d = z.getRows ();
1020 : 600 : matrix e, zref, gref;
1021 : :
1022 [ + - ][ + - ]: 600 : assert (d == z.getCols () && d == z0.getSize ());
[ - + ][ - + ]
1023 : :
1024 [ + - ][ + - ]: 600 : e = eye (d);
1025 [ + - ][ + - ]: 600 : zref = diagonal (z0);
[ + - ][ + - ]
1026 [ + - ][ + - ]: 600 : gref = diagonal (sqrt (real (1 / z0)));
[ + - ][ + - ]
[ + - ][ + - ]
[ + - ][ + - ]
[ + - ][ + - ]
1027 [ + - ][ + - ]: 600 : return gref * (z - zref) * inverse (z + zref) * inverse (gref);
[ + - ][ + - ]
[ + - ][ + - ]
[ + - ][ + - ]
[ + - ][ + - ]
[ + - ][ + - ]
[ + - ]
1028 : : }
1029 : :
1030 : : /*!\brief Convert impedance matrix to scattering parameters identic case
1031 : : \param[in] Z Impedance matrix
1032 : : \param[in] z0 Normalisation impedance
1033 : : \return Scattering matrix
1034 : : \todo Why not inline
1035 : : \todo z and z0 const
1036 : : */
1037 : 600 : matrix ztos (matrix z, nr_complex_t z0) {
1038 [ + - ][ + - ]: 600 : return ztos (z, qucs::vector (z.getRows (), z0));
1039 : : }
1040 : :
1041 : : /*!\brief impedance matrix to admittance matrix.
1042 : :
1043 : : Convert impedance matrix to admittance matrix. By definition
1044 : : \f$Y=Z^{-1}\f$
1045 : : \param[in] z impedance matrix
1046 : : \return Admittance matrix
1047 : : \todo Why not inline
1048 : : \todo z const
1049 : : */
1050 : 0 : matrix ztoy (matrix z) {
1051 [ # # ]: 0 : assert (z.getRows () == z.getCols ());
1052 [ # # ]: 0 : return inverse (z);
1053 : : }
1054 : :
1055 : : /*!\brief Scattering parameters to admittance matrix.
1056 : :
1057 : : Convert scattering parameters to admittance matrix.
1058 : : According to [1] eq (19):
1059 : : \f[
1060 : : Z=F^{-1} (I- S)^{-1} (SG + G^+) F
1061 : : \f]
1062 : : Where \f$S\f$ is the scattering matrix, \f$x^+\f$
1063 : : is the adjoint of x, I the identity matrix. The matrix
1064 : : F and G are diagonal matrix defined by:
1065 : : \f{align*}
1066 : : F_i&=\frac{1}{2\sqrt{\Re\text{e}\; Z_i}} \\
1067 : : G_i&=Z_i
1068 : : \f}
1069 : : Using the well know formula \f$(AB)^{-1}=B^{1}A^{1}\f$,
1070 : : we derivate:
1071 : : \f[
1072 : : Y=F^{-1} (SG+G^+)^{-1} (I-S) F
1073 : : \f]
1074 : : \param[in] s Scattering matrix
1075 : : \param[in] z0 Normalisation impedance
1076 : : \note We could safely drop the \f$1/2\f$ in \f$F\f$ because we compute
1077 : : \f$FXF^{-1}\f$ and therefore \f$1/2\f$ will simplify.
1078 : : \bug not correct if zref is complex
1079 : : \todo s and z0 const
1080 : : \return Admittance matrix
1081 : : */
1082 : 0 : matrix stoy (matrix s, qucs::vector z0) {
1083 : 0 : int d = s.getRows ();
1084 : 0 : matrix e, zref, gref;
1085 : :
1086 [ # # ][ # # ]: 0 : assert (d == s.getCols () && d == z0.getSize ());
[ # # ][ # # ]
1087 : :
1088 [ # # ][ # # ]: 0 : e = eye (d);
1089 [ # # ][ # # ]: 0 : zref = diagonal (z0);
[ # # ][ # # ]
1090 [ # # ][ # # ]: 0 : gref = diagonal (sqrt (real (1 / z0)));
[ # # ][ # # ]
[ # # ][ # # ]
[ # # ][ # # ]
[ # # ][ # # ]
1091 [ # # ][ # # ]: 0 : return inverse (gref) * inverse (s * zref + zref) * (e - s) * gref;
[ # # ][ # # ]
[ # # ][ # # ]
[ # # ][ # # ]
[ # # ][ # # ]
[ # # ][ # # ]
[ # # ][ # # ]
[ # # ]
1092 : : }
1093 : :
1094 : : /*!\brief Convert scattering pto adminttance parameters identic case
1095 : : \param[in] S Scattering matrix
1096 : : \param[in] z0 Normalisation impedance
1097 : : \return Admittance matrix
1098 : : \todo Why not inline
1099 : : \todo s and z0 const
1100 : : */
1101 : 0 : matrix stoy (matrix s, nr_complex_t z0) {
1102 [ # # ][ # # ]: 0 : return stoy (s, qucs::vector (s.getRows (), z0));
1103 : : }
1104 : :
1105 : : /*!\brief Admittance matrix to scattering parameters
1106 : :
1107 : : Convert admittance matrix to scattering parameters.
1108 : : Using the same methodology as [1] eq (16-19), but writing
1109 : : (16) as \f$i=Yv\f$, ie
1110 : : \f[
1111 : : S=F(I-G^+Y)(I-GY)^{-1}F^{-1}
1112 : : \f]
1113 : : Where \f$S\f$ is the scattering matrix, \f$x^+\f$
1114 : : is the adjoint of x, I the identity matrix. The matrix
1115 : : F and G are diagonal matrix defined by:
1116 : : \f{align*}
1117 : : F_i&=\frac{1}{2\sqrt{\Re\text{e}\; Z_i}} \\
1118 : : G_i&=Z_i
1119 : : \f}
1120 : : Using the well know formula \f$(AB)^{-1}=B^{1}A^{1}\f$,
1121 : : we derivate:
1122 : : \f[
1123 : : Y=F^{-1} (SG+G^+)^{-1} (I-S) F
1124 : : \f]
1125 : : \param[in] y admittance matrix
1126 : : \param[in] z0 Normalisation impedance
1127 : : \note We could safely drop the \f$1/2\f$ in \f$F\f$ because we compute
1128 : : \f$FXF^{-1}\f$ and therefore \f$1/2\f$ will simplify.
1129 : : \bug not correct if zref is complex
1130 : : \todo why not y and z0 const
1131 : : \return Scattering matrix
1132 : : */
1133 : 470 : matrix ytos (matrix y, qucs::vector z0) {
1134 : 470 : int d = y.getRows ();
1135 : 470 : matrix e, zref, gref;
1136 : :
1137 [ + - ][ + - ]: 470 : assert (d == y.getCols () && d == z0.getSize ());
[ - + ][ - + ]
1138 : :
1139 [ + - ][ + - ]: 470 : e = eye (d);
1140 [ + - ][ + - ]: 470 : zref = diagonal (z0);
[ + - ][ + - ]
1141 [ + - ][ + - ]: 470 : gref = diagonal (sqrt (real (1 / z0)));
[ + - ][ + - ]
[ + - ][ + - ]
[ + - ][ + - ]
[ + - ][ + - ]
1142 [ + - ][ + - ]: 470 : return gref * (e - zref * y) * inverse (e + zref * y) * inverse (gref);
[ + - ][ + - ]
[ + - ][ + - ]
[ + - ][ + - ]
[ + - ][ + - ]
[ + - ][ + - ]
[ + - ][ + - ]
[ + - ][ + - ]
[ + - ]
1143 : : }
1144 : : /*!\brief Convert Admittance matrix to scattering parameters identic case
1145 : : \param[in] y Admittance matrix
1146 : : \param[in] z0 Normalisation impedance
1147 : : \return Scattering matrix
1148 : : \todo Why not inline
1149 : : \todo y and z0 const
1150 : : */
1151 : 470 : matrix ytos (matrix y, nr_complex_t z0) {
1152 [ + - ][ + - ]: 470 : return ytos (y, qucs::vector (y.getRows (), z0));
1153 : : }
1154 : : /*!\brief Converts chain matrix to scattering parameters.
1155 : :
1156 : : Converts scattering parameters to chain matrix.
1157 : : Formulae are given by [5] tab 1. and are remembered here:
1158 : :
1159 : : \f{align*}
1160 : : A&=\frac{(Z_{01}^*+S_{11}Z_{01})(1-S_{22})
1161 : : +S_{12}S_{21}Z_{01}}{\Delta} \\
1162 : : B&=\frac{(Z_{01}^*+S_{11}Z_{01})(Z_{02}^*+S_{22}Z_{02})
1163 : : -S_{12}S_{21}Z_{01}Z_{02}}{\Delta} \\
1164 : : C&=\frac{(1-S_{11})(1-S_{22})
1165 : : -S_{12}S_{21}}{\Delta} \\
1166 : : D&=\frac{(1-S_{11})(Z_{02}^*+S_{22}Z_{02})
1167 : : +S_{12}S_{21}Z_{02}}{\Delta}
1168 : : \f}
1169 : : Where:
1170 : : \f[
1171 : : \Delta = 2 S_{21}\sqrt{\Re\text{e}\;Z_{01}\Re\text{e}\;Z_{02}}
1172 : : \f]
1173 : : \bug Do not need fabs
1174 : : \param[in] s Scattering matrix
1175 : : \param[in] z1 impedance at input 1
1176 : : \param[in] z2 impedance at input 2
1177 : : \return Chain matrix
1178 : : \note Assert 2 by 2 matrix
1179 : : \todo Why not s,z1,z2 const
1180 : : */
1181 : 0 : matrix stoa (matrix s, nr_complex_t z1, nr_complex_t z2) {
1182 [ # # ][ # # ]: 0 : nr_complex_t d = s (0, 0) * s (1, 1) - s (0, 1) * s (1, 0);
[ # # ]
1183 [ # # ]: 0 : nr_complex_t n = 2.0 * s (1, 0) * sqrt (fabs (real (z1) * real (z2)));
1184 [ # # ]: 0 : matrix a (2);
1185 : :
1186 [ # # ][ # # ]: 0 : assert (s.getRows () >= 2 && s.getCols () >= 2);
[ # # ]
1187 : :
1188 [ # # ]: 0 : a.set (0, 0, (conj (z1) + z1 * s (0, 0) -
1189 [ # # ]: 0 : conj (z1) * s (1, 1) - z1 * d) / n);
[ # # # # ]
[ # # ][ # # ]
[ # # ]
1190 [ # # ][ # # ]: 0 : a.set (0, 1, (conj (z1) * conj (z2) + z1 * conj (z2) * s (0, 0) +
[ # # ]
1191 [ # # ][ # # ]: 0 : conj (z1) * z2 * s (1, 1) + z1 * z2 * d) / n);
[ # # ][ # # ]
[ # # ][ # # ]
[ # # ][ # # ]
1192 [ # # ][ # # ]: 0 : a.set (1, 0, (1.0 - s (0, 0) - s (1, 1) + d) / n);
[ # # ]
1193 [ # # ]: 0 : a.set (1, 1, (conj (z2) - conj (z2) * s (0, 0) +
1194 [ # # ]: 0 : z2 * s (1, 1) - z2 * d) / n);
[ # # # # ]
[ # # ][ # # ]
[ # # ]
1195 : 0 : return a;
1196 : : }
1197 : :
1198 : :
1199 : : /*!\brief Converts chain matrix to scattering parameters.
1200 : :
1201 : : Converts chain matrix to scattering parameters
1202 : : Formulae are given by [5] and are remembered here:
1203 : : \f{align*}
1204 : : S_{11}&=\frac{AZ_{02}+B-CZ_{01}^*Z_{02}-DZ_{01}^*}{\Delta} \\
1205 : : S_{12}&=\frac{2(AD-BC)
1206 : : (\Re\text{e}\;Z_{01}\Re\text{e}\;Z_{02})^{1/2}}
1207 : : {\Delta}\\
1208 : : S_{21}&=\frac{2(\Re\text{e}\;Z_{01}\Re\text{e}\;Z_{02})^{1/2}}{\Delta}\\
1209 : : S_{22}&=\frac{-AZ_{02}^*+B-CZ_{01}^*Z_{02}+DZ_{01}}{\Delta}
1210 : : \f}
1211 : : Where:
1212 : : \f[
1213 : : \Delta =AZ_{02}+B+CZ_{01}Z_{02}-DZ_{01}
1214 : : \f]
1215 : : \param[in] a Chain matrix
1216 : : \param[in] z1 impedance at input 1
1217 : : \param[in] z2 impedance at input 2
1218 : : \return Scattering matrix
1219 : : \bug Do not use fabs
1220 : : \todo a, z1, z2 const
1221 : : */
1222 : 0 : matrix atos (matrix a, nr_complex_t z1, nr_complex_t z2) {
1223 [ # # ]: 0 : nr_complex_t d = 2.0 * sqrt (fabs (real (z1) * real (z2)));
1224 [ # # ]: 0 : nr_complex_t n = a (0, 0) * z2 + a (0, 1) +
1225 [ # # ][ # # ]: 0 : a (1, 0) * z1 * z2 + a (1, 1) * z1;
[ # # ][ # # ]
[ # # ][ # # ]
1226 [ # # ]: 0 : matrix s (2);
1227 : :
1228 [ # # ][ # # ]: 0 : assert (a.getRows () >= 2 && a.getCols () >= 2);
[ # # ]
1229 : :
1230 [ # # ]: 0 : s.set (0, 0, (a (0, 0) * z2 + a (0, 1)
1231 [ # # ][ # # ]: 0 : - a (1, 0) * conj (z1) * z2 - a (1, 1) * conj (z1)) / n);
[ # # ][ # # ]
[ # # ][ # # ]
[ # # ]
1232 : 0 : s.set (0, 1, (a (0, 0) * a (1, 1) -
1233 [ # # # # ]: 0 : a (0, 1) * a (1, 0)) * d / n);
[ # # ][ # # ]
[ # # ]
1234 [ # # ]: 0 : s.set (1, 0, d / n);
1235 [ # # ][ # # ]: 0 : s.set (1, 1, (a (1, 1) * z1 - a (0, 0) * conj (z2) +
1236 [ # # ][ # # ]: 0 : a (0, 1) - a (1, 0) * z1 * conj (z2)) / n);
[ # # ][ # # ]
[ # # ][ # # ]
1237 : 0 : return s;
1238 : : }
1239 : :
1240 : : /*!\brief Converts scattering parameters to hybrid matrix.
1241 : :
1242 : : Converts chain matrix to scattering parameters
1243 : : Formulae are given by [5] and are remembered here:
1244 : : \f{align*}
1245 : : h_{11}&=\frac{(Z_{01}^*+S_{11}Z_{01})
1246 : : (Z_{02}^*+S_{22}Z_{02})
1247 : : -S_{12}S_{21}Z_{01}Z_{02}}{\Delta}\\
1248 : : h_{12}&=\frac{2S_{12}
1249 : : (\Re\text{e}\;Z_{01} \Re\text{e}\;Z_{02})^\frac{1}{2}}
1250 : : {\Delta} \\
1251 : : h_{21}&=\frac{-2S_{21}
1252 : : (\Re\text{e}\;Z_{01} \Re\text{e}\;Z_{02})^\frac{1}{2}}
1253 : : {\Delta} \\
1254 : :
1255 : : h_{22}&=\frac{(1-S_{11})(1-S_{22})-S_{12}S_{21}}{\Delta}
1256 : : \f}
1257 : : Where \f$\Delta\f$ is:
1258 : : \f[
1259 : : \Delta=(1-S_{11})(Z_{02}^*+S_{22}Z_{02})+S_{12}S_{21}Z_{02}
1260 : : \f]
1261 : : \bug{Programmed formulae are valid only for Z real}
1262 : : \param[in] s Scattering matrix
1263 : : \param[in] z1 impedance at input 1
1264 : : \param[in] z2 impedance at input 2
1265 : : \return hybrid matrix
1266 : : \note Assert 2 by 2 matrix
1267 : : \todo Why not s,z1,z2 const
1268 : : */
1269 : 0 : matrix stoh (matrix s, nr_complex_t z1, nr_complex_t z2) {
1270 [ # # ]: 0 : nr_complex_t n = s (0, 1) * s (1, 0);
1271 [ # # ][ # # ]: 0 : nr_complex_t d = (1.0 - s (0, 0)) * (1.0 + s (1, 1)) + n;
1272 [ # # ]: 0 : matrix h (2);
1273 : :
1274 [ # # ][ # # ]: 0 : assert (s.getRows () >= 2 && s.getCols () >= 2);
[ # # ]
1275 : :
1276 [ # # ][ # # ]: 0 : h.set (0, 0, ((1.0 + s (0, 0)) * (1.0 + s (1, 1)) - n) * z1 / d);
[ # # ][ # # ]
1277 [ # # ]: 0 : h.set (0, 1, +2.0 * s (0, 1) / d);
1278 [ # # ]: 0 : h.set (1, 0, -2.0 * s (1, 0) / d);
1279 [ # # ][ # # ]: 0 : h.set (1, 1, ((1.0 - s (0, 0)) * (1.0 - s (1, 1)) - n) / z2 / d);
[ # # ][ # # ]
1280 : 0 : return h;
1281 : : }
1282 : :
1283 : : /*!\brief Converts hybrid matrix to scattering parameters.
1284 : :
1285 : : Formulae are given by [5] and are remembered here:
1286 : : \f{align*}
1287 : : S_{11}&=\frac{(h_{11}-Z_{01}^*)(1+h_{22}Z_{02})-h_{12}h_{21}Z_{02}}
1288 : : {\Delta}\\
1289 : : S_{12}&=\frac{-2h_{12}(\Re\text{e}\;Z_{01}\Re\text{e}\;Z_{02})^\frac{1}{2}}
1290 : : {\Delta}\\
1291 : : S_{21}&=\frac{-2h_{21}(\Re\text{e}\;Z_{01}\Re\text{e}\;Z_{02})^\frac{1}{2}}
1292 : : {\Delta}\\
1293 : : S_{22}&=\frac{(h_{11}+Z_{01})(1-h_{22}Z_{02}^*)-h_{12}h_{21}Z_{02}^*}
1294 : : {\Delta}
1295 : : \f}
1296 : : Where \f$\Delta\f$ is:
1297 : : \f[
1298 : : \Delta=(Z_{01}+h_{11})(1+h_{22}Z_{02})-h_{12}h_{21}Z_{02}
1299 : : \f]
1300 : : \param[in] h hybrid matrix
1301 : : \param[in] z1 impedance at input 1
1302 : : \param[in] z2 impedance at input 2
1303 : : \return scattering matrix
1304 : : \note Assert 2 by 2 matrix
1305 : : \todo Why not h,z1,z2 const
1306 : : */
1307 : 0 : matrix htos (matrix h, nr_complex_t z1, nr_complex_t z2) {
1308 [ # # ]: 0 : nr_complex_t n = h (0, 1) * h (1, 0);
1309 [ # # ][ # # ]: 0 : nr_complex_t d = (1.0 + h (0, 0) / z1) * (1.0 + z2 * h (1, 1)) - n;
[ # # ][ # # ]
1310 [ # # ]: 0 : matrix s (2);
1311 : :
1312 [ # # ][ # # ]: 0 : assert (h.getRows () >= 2 && h.getCols () >= 2);
[ # # ]
1313 : :
1314 [ # # ][ # # ]: 0 : s.set (0, 0, ((h (0, 0) / z1 - 1.0) * (1.0 + z2 * h (1, 1)) - n) / d);
[ # # ][ # # ]
[ # # ]
1315 [ # # ]: 0 : s.set (0, 1, +2.0 * h (0, 1) / d);
1316 [ # # ]: 0 : s.set (1, 0, -2.0 * h (1, 0) / d);
1317 [ # # ][ # # ]: 0 : s.set (1, 1, ((1.0 + h (0, 0) / z1) * (1.0 - z2 * h (1, 1)) + n) / d);
[ # # ][ # # ]
[ # # ]
1318 : 0 : return s;
1319 : : }
1320 : :
1321 : : /*\brief Converts scattering parameters to second hybrid matrix.
1322 : : \bug{Programmed formulae are valid only for Z real}
1323 : : \bug{Not documented and references}
1324 : : \param[in] s Scattering matrix
1325 : : \param[in] z1 impedance at input 1
1326 : : \param[in] z2 impedance at input 2
1327 : : \return second hybrid matrix
1328 : : \note Assert 2 by 2 matrix
1329 : : \todo Why not s,z1,z2 const
1330 : : */
1331 : 0 : matrix stog (matrix s, nr_complex_t z1, nr_complex_t z2) {
1332 [ # # ]: 0 : nr_complex_t n = s (0, 1) * s (1, 0);
1333 [ # # ][ # # ]: 0 : nr_complex_t d = (1.0 + s (0, 0)) * (1.0 - s (1, 1)) + n;
1334 [ # # ]: 0 : matrix g (2);
1335 : :
1336 [ # # ][ # # ]: 0 : assert (s.getRows () >= 2 && s.getCols () >= 2);
[ # # ]
1337 : :
1338 [ # # ][ # # ]: 0 : g.set (0, 0, ((1.0 - s (0, 0)) * (1.0 - s (1, 1)) - n) / z1 / d);
[ # # ][ # # ]
1339 [ # # ]: 0 : g.set (0, 1, -2.0 * s (0, 1) / d);
1340 [ # # ]: 0 : g.set (1, 0, +2.0 * s (1, 0) / d);
1341 [ # # ][ # # ]: 0 : g.set (1, 1, ((1.0 + s (0, 0)) * (1.0 + s (1, 1)) - n) * z2 / d);
[ # # ][ # # ]
1342 : 0 : return g;
1343 : : }
1344 : :
1345 : : /*\brief Converts second hybrid matrix to scattering parameters.
1346 : : \bug{Programmed formulae are valid only for Z real}
1347 : : \bug{Not documented and references}
1348 : : \param[in] g second hybrid matrix
1349 : : \param[in] z1 impedance at input 1
1350 : : \param[in] z2 impedance at input 2
1351 : : \return scattering matrix
1352 : : \note Assert 2 by 2 matrix
1353 : : \todo Why not g,z1,z2 const
1354 : : */
1355 : 0 : matrix gtos (matrix g, nr_complex_t z1, nr_complex_t z2) {
1356 [ # # ]: 0 : nr_complex_t n = g (0, 1) * g (1, 0);
1357 [ # # ][ # # ]: 0 : nr_complex_t d = (1.0 + g (0, 0) * z1) * (1.0 + g (1, 1) / z2) - n;
[ # # ][ # # ]
1358 [ # # ]: 0 : matrix s (2);
1359 : :
1360 [ # # ][ # # ]: 0 : assert (g.getRows () >= 2 && g.getCols () >= 2);
[ # # ]
1361 : :
1362 [ # # ][ # # ]: 0 : s.set (0, 0, ((1.0 - g (0, 0) * z1) * (1.0 + g (1, 1) / z2) + n) / d);
[ # # ][ # # ]
[ # # ]
1363 [ # # ]: 0 : s.set (0, 1, -2.0 * g (0, 1) / d);
1364 [ # # ]: 0 : s.set (1, 0, +2.0 * g (1, 0) / d);
1365 [ # # ][ # # ]: 0 : s.set (1, 1, ((g (0, 0) * z1 + 1.0) * (g (1, 1) / z2 - 1.0) - n) / d);
[ # # ][ # # ]
[ # # ]
1366 : 0 : return s;
1367 : : }
1368 : :
1369 : : /*!\brief Convert admittance matrix to impedance matrix.
1370 : :
1371 : : Convert \f$Y\f$ matrix to \f$Z\f$ matrix using well known relation
1372 : : \f$Z=Y^{-1}\f$
1373 : : \param[in] y admittance matrix
1374 : : \return Impedance matrix
1375 : : \note Check if y matrix is a square matrix
1376 : : \todo Why not inline
1377 : : \todo y const
1378 : : \todo move near ztoy()
1379 : : */
1380 : 0 : matrix ytoz (matrix y) {
1381 [ # # ]: 0 : assert (y.getRows () == y.getCols ());
1382 [ # # ]: 0 : return inverse (y);
1383 : : }
1384 : :
1385 : : /*!\brief Admittance noise correlation matrix to S-parameter noise
1386 : : correlation matrix
1387 : :
1388 : : Converts admittance noise correlation matrix to S-parameter noise
1389 : : correlation matrix. According to [7] fig 2:
1390 : : \f[
1391 : : C_s=\frac{1}{4}(I+S)C_y(I+S)^+
1392 : : \f]
1393 : : Where \f$C_s\f$ is the scattering noise correlation matrix,
1394 : : \f$C_y\f$ the admittance noise correlation matrix, \f$I\f$
1395 : : the identity matrix and \f$S\f$ the scattering matrix
1396 : : of device. \f$x^+\f$ is the adjoint of \f$x\f$
1397 : : \warning cy matrix and s matrix are assumed to be normalized
1398 : : \param[in] cy Admittance noise correlation
1399 : : \param[in] s S parameter matrix of device
1400 : : \return S-parameter noise correlation matrix
1401 : : \note Assert compatiblity of matrix
1402 : : \todo cy s const
1403 : : */
1404 : 70 : matrix cytocs (matrix cy, matrix s) {
1405 [ + - ]: 70 : matrix e = eye (s.getRows ());
1406 : :
1407 [ + - ][ + - ]: 70 : assert (cy.getRows () == cy.getCols () && s.getRows () == s.getCols () &&
[ - + ]
1408 [ - + ]: 70 : cy.getRows () == s.getRows ());
1409 : :
1410 [ + - ][ + - ]: 70 : return (e + s) * cy * adjoint (e + s) / 4;
[ + - ][ + - ]
[ + - ][ + - ]
[ + - ][ + - ]
[ + - ][ + - ]
[ + - ]
1411 : : }
1412 : :
1413 : : /*!\brief Converts S-parameter noise correlation matrix to admittance noise
1414 : : correlation matrix.
1415 : :
1416 : : According to [7] fig 2:
1417 : : \f[
1418 : : C_y=(I+Y)C_s(I+Y)^+
1419 : : \f]
1420 : : Where \f$C_s\f$ is the scattering noise correlation matrix,
1421 : : \f$C_y\f$ the admittance noise correlation matrix, \f$I\f$
1422 : : the identity matrix and \f$S\f$ the scattering matrix
1423 : : of device. \f$x^+\f$ is the adjoint of \f$x\f$
1424 : : \warning cs matrix and y matrix are assumed to be normalized
1425 : : \param[in] cs S parameter noise correlation
1426 : : \param[in] y Admittance matrix of device
1427 : : \return admittance noise correlation matrix
1428 : : \todo cs, y const
1429 : : */
1430 : 0 : matrix cstocy (matrix cs, matrix y) {
1431 [ # # ]: 0 : matrix e = eye (y.getRows ());
1432 : :
1433 [ # # ][ # # ]: 0 : assert (cs.getRows () == cs.getCols () && y.getRows () == y.getCols () &&
[ # # ]
1434 [ # # ]: 0 : cs.getRows () == y.getRows ());
1435 : :
1436 [ # # ][ # # ]: 0 : return (e + y) * cs * adjoint (e + y);
[ # # ][ # # ]
[ # # ][ # # ]
[ # # ][ # # ]
[ # # ][ # # ]
1437 : : }
1438 : :
1439 : : /*!\brief Converts impedance noise correlation matrix to S-parameter noise
1440 : : correlation matrix.
1441 : :
1442 : : According to [7] fig 2:
1443 : : \f[
1444 : : C_s=\frac{1}{4}(I-S)C_z(I-S)
1445 : : \f]
1446 : : Where \f$C_s\f$ is the scattering noise correlation matrix,
1447 : : \f$C_z\f$ the impedance noise correlation matrix, \f$I\f$
1448 : : the identity matrix and \f$S\f$ the scattering matrix
1449 : : of device. \f$x^+\f$ is the adjoint of \f$x\f$
1450 : : \warning Cz matrix and s matrix are assumed to be normalized
1451 : : \param[in] cz Impedance noise correlation
1452 : : \param[in] s S parameter matrix of device
1453 : : \return S-parameter noise correlation matrix
1454 : : \note Assert compatiblity of matrix
1455 : : \todo cz, s const
1456 : : */
1457 : 0 : matrix cztocs (matrix cz, matrix s) {
1458 [ # # ]: 0 : matrix e = eye (s.getRows ());
1459 : :
1460 [ # # ][ # # ]: 0 : assert (cz.getRows () == cz.getCols () && s.getRows () == s.getCols () &&
[ # # ]
1461 [ # # ]: 0 : cz.getRows () == s.getRows ());
1462 : :
1463 [ # # ][ # # ]: 0 : return (e - s) * cz * adjoint (e - s) / 4;
[ # # ][ # # ]
[ # # ][ # # ]
[ # # ][ # # ]
[ # # ][ # # ]
[ # # ]
1464 : : }
1465 : :
1466 : : /*!\brief Converts S-parameter noise correlation matrix to impedance noise
1467 : : correlation matrix.
1468 : :
1469 : : According to [7] fig 2:
1470 : : \f[
1471 : : C_z=(I+Z)C_s(I+Z)^+
1472 : : \f]
1473 : : Where \f$C_s\f$ is the scattering noise correlation matrix,
1474 : : \f$C_z\f$ the impedance noise correlation matrix, \f$I\f$
1475 : : the identity matrix and \f$S\f$ the scattering matrix
1476 : : of device. \f$x^+\f$ is the adjoint of \f$x\f$
1477 : : \warning cs matrix and y matrix are assumed to be normalized
1478 : : \param[in] cs S parameter noise correlation
1479 : : \param[in] z Impedance matrix of device
1480 : : \return Impedance noise correlation matrix
1481 : : \todo cs, z const
1482 : : */
1483 : 0 : matrix cstocz (matrix cs, matrix z) {
1484 [ # # ][ # # ]: 0 : assert (cs.getRows () == cs.getCols () && z.getRows () == z.getCols () &&
[ # # ]
1485 [ # # ]: 0 : cs.getRows () == z.getRows ());
1486 [ # # ]: 0 : matrix e = eye (z.getRows ());
1487 [ # # ][ # # ]: 0 : return (e + z) * cs * adjoint (e + z);
[ # # ][ # # ]
[ # # ][ # # ]
[ # # ][ # # ]
[ # # ][ # # ]
1488 : : }
1489 : :
1490 : : /*!\brief Converts impedance noise correlation matrix to admittance noise
1491 : : correlation matrix.
1492 : :
1493 : : According to [7] fig 2:
1494 : : \f[
1495 : : C_y=YC_zY^+
1496 : : \f]
1497 : : Where \f$C_z\f$ is the impedance correlation matrix,
1498 : : \f$I\f$ the identity matrix and \f$C_y\f$ the admittance noise
1499 : : correlation matrix.
1500 : : \f$x^+\f$ is the adjoint of \f$x\f$
1501 : : \warning cz matrix and y matrix are assumed to be normalized
1502 : : \param[in] cz impedance noise correlation
1503 : : \param[in] y Admittance matrix of device
1504 : : \return admittance noise correlation matrix
1505 : : \todo cs, y const
1506 : : */
1507 : 0 : matrix cztocy (matrix cz, matrix y) {
1508 [ # # ][ # # ]: 0 : assert (cz.getRows () == cz.getCols () && y.getRows () == y.getCols () &&
[ # # ]
1509 [ # # ]: 0 : cz.getRows () == y.getRows ());
1510 : :
1511 [ # # ][ # # ]: 0 : return y * cz * adjoint (y);
[ # # ][ # # ]
[ # # ]
1512 : : }
1513 : :
1514 : : /*!\brief Converts admittance noise correlation matrix to impedance noise
1515 : : correlation matrix.
1516 : :
1517 : : According to [7] fig 2:
1518 : : \f[
1519 : : C_z=ZC_yZ^+
1520 : : \f]
1521 : : Where \f$C_z\f$ is the impedance correlation matrix,
1522 : : \f$I\f$ the identity matrix and \f$C_y\f$ the admittance noise
1523 : : correlation matrix.
1524 : : \f$x^+\f$ is the adjoint of \f$x\f$
1525 : : \warning cy matrix and z matrix are assumed to be normalized
1526 : : \param[in] cy Admittance noise correlation
1527 : : \param[in] z Impedance matrix of device
1528 : : \return Impedance noise correlation matrix
1529 : : \todo cs, z const
1530 : : */
1531 : 0 : matrix cytocz (matrix cy, matrix z) {
1532 [ # # ][ # # ]: 0 : assert (cy.getRows () == cy.getCols () && z.getRows () == z.getCols () &&
[ # # ]
1533 [ # # ]: 0 : cy.getRows () == z.getRows ());
1534 [ # # ][ # # ]: 0 : return z * cy * adjoint (z);
[ # # ][ # # ]
[ # # ]
1535 : : }
1536 : :
1537 : : /*!\brief The function swaps the given rows with each other.
1538 : : \param[in] r1 source row
1539 : : \param[in] r2 destination row
1540 : : \note assert not out of bound r1 and r2
1541 : : \todo r1 and r2 const
1542 : : */
1543 : 1000 : void matrix::exchangeRows (int r1, int r2) {
1544 [ + + ]: 3800 : nr_complex_t * s = new nr_complex_t[cols];
1545 : 1000 : int len = sizeof (nr_complex_t) * cols;
1546 : :
1547 [ + - ][ + - ]: 1000 : assert (r1 >= 0 && r2 >= 0 && r1 < rows && r2 < rows);
[ + - ][ - + ]
[ - + ]
1548 : :
1549 : 1000 : memcpy (s, &data[r1 * cols], len);
1550 : 1000 : memcpy (&data[r1 * cols], &data[r2 * cols], len);
1551 : 1000 : memcpy (&data[r2 * cols], s, len);
1552 [ + - ]: 1000 : delete[] s;
1553 : 1000 : }
1554 : :
1555 : : /*!\brief The function swaps the given column with each other.
1556 : : \param[in] c1 source column
1557 : : \param[in] c2 destination column
1558 : : \note assert not out of bound r1 and r2
1559 : : \todo c1 and c2 const
1560 : : */
1561 : 0 : void matrix::exchangeCols (int c1, int c2) {
1562 : 0 : nr_complex_t s;
1563 : :
1564 [ # # ][ # # ]: 0 : assert (c1 >= 0 && c2 >= 0 && c1 < cols && c2 < cols);
[ # # ][ # # ]
[ # # ]
1565 : :
1566 [ # # ]: 0 : for (int r = 0; r < rows * cols; r += cols) {
1567 : 0 : s = data[r + c1];
1568 : 0 : data[r + c1] = data[r + c2];
1569 : 0 : data[r + c2] = s;
1570 : : }
1571 : 0 : }
1572 : :
1573 : : /*!\brief Generic conversion matrix
1574 : :
1575 : : This function converts 2x2 matrices from any of the matrix forms Y,
1576 : : Z, H, G and A to any other. Also converts S<->(A, T, H, Y and Z)
1577 : : matrices.
1578 : : Convertion assumed:
1579 : :
1580 : : Y->Y, Y->Z, Y->H, Y->G, Y->A, Y->S,
1581 : : Z->Y, Z->Z, Z->H, Z->G, Z->A, Z->S,
1582 : : H->Y, H->Z, H->H, H->G, H->A, H->S,
1583 : : G->Y, G->Z, G->H, G->G, G->A, G->S,
1584 : : A->Y, A->Z, A->H, A->G, A->A, A->S,
1585 : : S->Y, S->Z, S->H, S->G, S->A, S->S,
1586 : : S->T,T->T,T->S
1587 : : \note assert 2x2 matrix
1588 : : \param[in] m base matrix
1589 : : \param[in] in matrix
1590 : : \param[in] out matrix
1591 : : \return matrix given by format out
1592 : : \todo m, in, out const
1593 : : */
1594 : 0 : matrix twoport (matrix m, char in, char out) {
1595 [ # # ][ # # ]: 0 : assert (m.getRows () >= 2 && m.getCols () >= 2);
[ # # ]
1596 : 0 : nr_complex_t d;
1597 [ # # ]: 0 : matrix res (2);
1598 : :
1599 [ # # # # : 0 : switch (in) {
# # # # ]
1600 : : case 'Y':
1601 [ # # # # : 0 : switch (out) {
# # # ]
1602 : : case 'Y': // Y to Y
1603 [ # # ]: 0 : res = m;
1604 : 0 : break;
1605 : : case 'Z': // Y to Z
1606 [ # # ][ # # ]: 0 : d = m (0, 0) * m (1, 1) - m (0, 1) * m (1, 0);
[ # # ]
1607 [ # # ]: 0 : res.set (0, 0, m (1, 1) / d);
1608 [ # # ]: 0 : res.set (0, 1, -m (0, 1) / d);
1609 [ # # ]: 0 : res.set (1, 0, -m (1, 0) / d);
1610 [ # # ]: 0 : res.set (1, 1, m (0, 0) / d);
1611 : 0 : break;
1612 : : case 'H': // Y to H
1613 : 0 : d = m (0, 0);
1614 : 0 : res.set (0, 0, 1.0 / d);
1615 [ # # ]: 0 : res.set (0, 1, -m (0, 1) / d);
1616 [ # # ]: 0 : res.set (1, 0, m (1, 0) / d);
1617 [ # # ][ # # ]: 0 : res.set (1, 1, m (1, 1) - m (0, 1) * m (1, 0) / d);
[ # # ]
1618 : 0 : break;
1619 : : case 'G': // Y to G
1620 : 0 : d = m (1, 1);
1621 [ # # ][ # # ]: 0 : res.set (0, 0, m (0, 0) - m (0, 1) * m (1, 0) / d);
[ # # ]
1622 [ # # ]: 0 : res.set (0, 1, m (0, 1) / d);
1623 [ # # ]: 0 : res.set (1, 0, -m (1, 0) / d);
1624 : 0 : res.set (1, 1, 1.0 / d);
1625 : 0 : break;
1626 : : case 'A': // Y to A
1627 : 0 : d = m (1, 0);
1628 [ # # ]: 0 : res.set (0, 0, -m (1, 1) / d);
1629 : 0 : res.set (0, 1, -1.0 / d);
1630 [ # # ][ # # ]: 0 : res.set (1, 0, m (0, 1) - m (1, 1) * m (0, 0) / d);
[ # # ]
1631 [ # # ]: 0 : res.set (1, 1, -m (0, 0) / d);
1632 : 0 : break;
1633 : : case 'S': // Y to S
1634 [ # # ][ # # ]: 0 : res = ytos (m);
[ # # ]
1635 : 0 : break;
1636 : : }
1637 : 0 : break;
1638 : : case 'Z':
1639 [ # # # # : 0 : switch (out) {
# # # ]
1640 : : case 'Y': // Z to Y
1641 [ # # ][ # # ]: 0 : d = m (0, 0) * m (1, 1) - m (0, 1) * m (1, 0);
[ # # ]
1642 [ # # ]: 0 : res.set (0, 0, m (1, 1) / d);
1643 [ # # ]: 0 : res.set (0, 1, -m (0, 1) / d);
1644 [ # # ]: 0 : res.set (1, 0, -m (1, 0) / d);
1645 [ # # ]: 0 : res.set (1, 1, m (0, 0) / d);
1646 : 0 : break;
1647 : : case 'Z': // Z to Z
1648 [ # # ]: 0 : res = m;
1649 : 0 : break;
1650 : : case 'H': // Z to H
1651 : 0 : d = m (1, 1);
1652 [ # # ][ # # ]: 0 : res.set (0, 0, m (0, 0) - m (0, 1) * m (1, 0) / d);
[ # # ]
1653 [ # # ]: 0 : res.set (0, 1, m (0, 1) / d);
1654 [ # # ]: 0 : res.set (1, 0, -m (1, 0) / d);
1655 : 0 : res.set (1, 1, 1.0 / d);
1656 : 0 : break;
1657 : : case 'G': // Z to G
1658 : 0 : d = m (0, 0);
1659 : 0 : res.set (0, 0, 1.0 / d);
1660 [ # # ]: 0 : res.set (0, 1, -m (0, 1) / d);
1661 [ # # ]: 0 : res.set (1, 0, m (1, 0) / d);
1662 [ # # ][ # # ]: 0 : res.set (1, 1, m (1, 1) - m (0, 1) * m (1, 0) / d);
[ # # ]
1663 : 0 : break;
1664 : : case 'A': // Z to A
1665 : 0 : d = m (1, 0);
1666 [ # # ]: 0 : res.set (0, 0, m (0, 0) / d);
1667 [ # # ][ # # ]: 0 : res.set (0, 1, m (0, 0) * m (1, 1) / d - m (0, 1));
[ # # ]
1668 : 0 : res.set (1, 0, 1.0 / d);
1669 [ # # ]: 0 : res.set (1, 1, m (1, 1) / d);
1670 : 0 : break;
1671 : : case 'S': // Z to S
1672 [ # # ][ # # ]: 0 : res = ztos (m);
[ # # ]
1673 : 0 : break;
1674 : : }
1675 : 0 : break;
1676 : : case 'H':
1677 [ # # # # : 0 : switch (out) {
# # # ]
1678 : : case 'Y': // H to Y
1679 : 0 : d = m (0, 0);
1680 : 0 : res.set (0, 0, 1.0 / d);
1681 [ # # ]: 0 : res.set (0, 1, -m (0, 1) / d);
1682 [ # # ]: 0 : res.set (1, 0, m (1, 0) / d);
1683 [ # # ][ # # ]: 0 : res.set (1, 1, m (1, 1) - m (0, 1) * m.get(2, 1) / d);
[ # # ]
1684 : 0 : break;
1685 : : case 'Z': // H to Z
1686 : 0 : d = m (1, 1);
1687 [ # # ][ # # ]: 0 : res.set (0, 0, m (0, 0) - m (0, 1) * m (1, 0) / d);
[ # # ]
1688 [ # # ]: 0 : res.set (0, 1, m (0, 1) / d);
1689 [ # # ]: 0 : res.set (1, 0, -m (1, 0) / d);
1690 : 0 : res.set (1, 1, 1.0 / d);
1691 : 0 : break;
1692 : : case 'H': // H to H
1693 [ # # ]: 0 : res = m;
1694 : 0 : break;
1695 : : case 'G': // H to G
1696 [ # # ][ # # ]: 0 : d = m (0, 0) * m (1, 1) - m (0, 1) * m (1, 0);
[ # # ]
1697 [ # # ]: 0 : res.set (0, 0, m (1, 1) / d);
1698 [ # # ]: 0 : res.set (0, 1, -m (0, 1) / d);
1699 [ # # ]: 0 : res.set (1, 0, -m (1, 0) / d);
1700 [ # # ]: 0 : res.set (1, 1, m (0, 0) / d);
1701 : 0 : break;
1702 : : case 'A': // H to A
1703 : 0 : d = m (1, 0);
1704 [ # # ][ # # ]: 0 : res.set (0, 0, m (0, 1) - m (0, 0) * m (1, 1) / d);
[ # # ]
1705 [ # # ]: 0 : res.set (0, 1, -m (0, 0) / d);
1706 [ # # ]: 0 : res.set (1, 0, -m (1, 1) / d);
1707 : 0 : res.set (1, 1, -1.0 / d);
1708 : 0 : break;
1709 : : case 'S': // H to S
1710 [ # # ][ # # ]: 0 : res = htos (m);
[ # # ]
1711 : 0 : break;
1712 : : }
1713 : 0 : break;
1714 : : case 'G':
1715 [ # # # # : 0 : switch (out) {
# # # ]
1716 : : case 'Y': // G to Y
1717 : 0 : d = m (1, 1);
1718 [ # # ][ # # ]: 0 : res.set (0, 0, m (0, 0) - m (0, 1) * m (1, 0) / d);
[ # # ]
1719 [ # # ]: 0 : res.set (0, 1, m (0, 1) / d);
1720 [ # # ]: 0 : res.set (1, 0, -m (1, 0) / d);
1721 : 0 : res.set (1, 1, 1.0 / d);
1722 : 0 : break;
1723 : : case 'Z': // G to Z
1724 : 0 : d = m (0, 0);
1725 : 0 : res.set (0, 0, 1.0 / d);
1726 [ # # ]: 0 : res.set (0, 1, -m (0, 1) / d);
1727 [ # # ]: 0 : res.set (1, 0, m (1, 0) / d);
1728 [ # # ][ # # ]: 0 : res.set (1, 1, m (1, 1) - m (0, 1) * m (1, 0) / d);
[ # # ]
1729 : 0 : break;
1730 : : case 'H': // G to H
1731 [ # # ][ # # ]: 0 : d = m (0, 0) * m (1, 1) - m (0, 1) * m (1, 0);
[ # # ]
1732 [ # # ]: 0 : res.set (0, 0, m (1, 1) / d);
1733 [ # # ]: 0 : res.set (0, 1, -m (0, 1) / d);
1734 [ # # ]: 0 : res.set (1, 0, -m (1, 0) / d);
1735 [ # # ]: 0 : res.set (1, 1, m (0, 0) / d);
1736 : 0 : break;
1737 : : case 'G': // G to G
1738 [ # # ]: 0 : res = m;
1739 : 0 : break;
1740 : : case 'A': // G to A
1741 : 0 : d = m (1, 0);
1742 : 0 : res.set (0, 0, 1.0 / d);
1743 [ # # ]: 0 : res.set (0, 1, m (1, 1) / d);
1744 [ # # ]: 0 : res.set (1, 0, m (0, 0) / d);
1745 [ # # ][ # # ]: 0 : res.set (1, 1, m (0, 0) * m (1, 1) / d - m (0, 1));
[ # # ]
1746 : 0 : break;
1747 : : case 'S': // G to S
1748 [ # # ][ # # ]: 0 : res = gtos (m);
[ # # ]
1749 : 0 : break;
1750 : : }
1751 : 0 : break;
1752 : : case 'A':
1753 [ # # # # : 0 : switch (out) {
# # # ]
1754 : : case 'Y': // A to Y
1755 : 0 : d = m (0, 1);
1756 [ # # ]: 0 : res.set (0, 0, m (1, 1) / d);
1757 [ # # ][ # # ]: 0 : res.set (0, 1, m (1, 0) - m (0, 0) * m (1, 1) / d);
[ # # ]
1758 : 0 : res.set (1, 0, -1.0 / d);
1759 [ # # ]: 0 : res.set (1, 1, m (0, 0) / d);
1760 : 0 : break;
1761 : : case 'Z': // A to Z
1762 : 0 : d = m (1, 0);
1763 [ # # ]: 0 : res.set (0, 0, m (0, 0) / d);
1764 [ # # ][ # # ]: 0 : res.set (0, 1, m (0, 0) * m (1, 1) / d - m (0, 1));
[ # # ]
1765 : 0 : res.set (1, 0, 1.0 / d);
1766 [ # # ]: 0 : res.set (1, 1, m (1, 1) / d);
1767 : 0 : break;
1768 : : case 'H': // A to H
1769 : 0 : d = m (1, 1);
1770 [ # # ]: 0 : res.set (0, 0, m (0, 1) / d);
1771 [ # # ][ # # ]: 0 : res.set (0, 1, m (0, 0) - m (0, 1) * m (1, 0) / d);
[ # # ]
1772 : 0 : res.set (1, 0, -1.0 / d);
1773 [ # # ]: 0 : res.set (1, 1, m (1, 0) / d);
1774 : 0 : break;
1775 : : case 'G': // A to G
1776 : 0 : d = m (0, 0);
1777 [ # # ]: 0 : res.set (0, 0, m (1, 0) / d);
1778 [ # # ][ # # ]: 0 : res.set (0, 1, m (1, 0) * m (0, 1) / d - m (1, 1));
[ # # ]
1779 : 0 : res.set (1, 0, 1.0 / d);
1780 [ # # ]: 0 : res.set (1, 1, m (0, 1) / d);
1781 : 0 : break;
1782 : : case 'A': // A to A
1783 [ # # ]: 0 : res = m;
1784 : 0 : break;
1785 : : case 'S': // A to S
1786 [ # # ][ # # ]: 0 : res = atos (m);
[ # # ]
1787 : 0 : break;
1788 : : }
1789 : 0 : break;
1790 : : case 'S':
1791 [ # # # # : 0 : switch (out) {
# # # # ]
1792 : : case 'S': // S to S
1793 [ # # ]: 0 : res = m;
1794 : 0 : break;
1795 : : case 'T': // S to T
1796 : 0 : d = m (1, 0);
1797 [ # # ][ # # ]: 0 : res.set (0, 0, m (0, 1) - m (0, 0) * m (1, 1) / d);
[ # # ]
1798 [ # # ]: 0 : res.set (0, 1, m (0, 0) / d);
1799 [ # # ]: 0 : res.set (1, 0, -m (1, 1) / d);
1800 : 0 : res.set (0, 1, 1.0 / d);
1801 : 0 : break;
1802 : : case 'A': // S to A
1803 [ # # ][ # # ]: 0 : res = stoa (m);
[ # # ]
1804 : 0 : break;
1805 : : case 'H': // S to H
1806 [ # # ][ # # ]: 0 : res = stoh (m);
[ # # ]
1807 : 0 : break;
1808 : : case 'G': // S to G
1809 [ # # ][ # # ]: 0 : res = stog (m);
[ # # ]
1810 : 0 : break;
1811 : : case 'Y': // S to Y
1812 [ # # ][ # # ]: 0 : res = stoy (m);
[ # # ]
1813 : 0 : break;
1814 : : case 'Z': // S to Z
1815 [ # # ][ # # ]: 0 : res = stoz (m);
[ # # ]
1816 : 0 : break;
1817 : : }
1818 : 0 : break;
1819 : : case 'T':
1820 [ # # # ]: 0 : switch (out) {
1821 : : case 'S': // T to S
1822 : 0 : d = m (1, 1);
1823 [ # # ]: 0 : res.set (0, 0, m (0, 1) / d);
1824 [ # # ][ # # ]: 0 : res.set (0, 1, m (0, 0) - m (0, 1) * m (1, 0) / d);
[ # # ]
1825 : 0 : res.set (1, 0, 1.0 / d);
1826 [ # # ]: 0 : res.set (0, 1, -m (1, 0) / d);
1827 : 0 : break;
1828 : : case 'T': // T to T
1829 [ # # ]: 0 : res = m;
1830 : 0 : break;
1831 : : }
1832 : 0 : break;
1833 : : }
1834 : 0 : return res;
1835 : : }
1836 : :
1837 : : /*!\brief Compute the Rollet stabilty factor
1838 : :
1839 : : The function returns the Rollet stability factor (\f$K\f) of the given
1840 : : S-parameter matrix:
1841 : : \[
1842 : : K=\frac{1-|S_{11}|^2-|S_{22}|^2+|\delta|^2}{2|S_{12}S_{21}|}
1843 : : \]
1844 : : Where:
1845 : : \[
1846 : : \Delta=S_{11}S_{22}-S_{12}S_{21}
1847 : : \]
1848 : : \param[in] m S parameter matrix
1849 : : \return Rollet factor
1850 : : \note Assert 2x2 matrix
1851 : : \todo m const?
1852 : : \todo Rewrite with abs and expand det. It is cleaner.
1853 : : */
1854 : 0 : nr_double_t rollet (matrix m) {
1855 [ # # ][ # # ]: 0 : assert (m.getRows () >= 2 && m.getCols () >= 2);
[ # # ]
1856 : : nr_double_t res;
1857 [ # # ][ # # ]: 0 : res = (1 - norm (m (0, 0)) - norm (m (1, 1)) + norm (det (m))) /
1858 [ # # ]: 0 : 2 / abs (m (0, 1) * m (1, 0));
1859 : 0 : return res;
1860 : : }
1861 : :
1862 : : /* Computes stability measure B1 of the given S-parameter matrix. */
1863 : 0 : nr_double_t b1 (matrix m) {
1864 [ # # ][ # # ]: 0 : assert (m.getRows () >= 2 && m.getCols () >= 2);
[ # # ]
1865 : : nr_double_t res;
1866 [ # # ][ # # ]: 0 : res = 1 + norm (m (0, 0)) - norm (m (1, 1)) - norm (det (m));
1867 : 0 : return res;
1868 : : }
1869 : :
1870 : : } // namespace qucs
|
