Qucs-core
0.0.19
|
00001 /* 00002 * matrix.cpp - matrix class implementation 00003 * 00004 * Copyright (C) 2003-2009 Stefan Jahn <stefan@lkcc.org> 00005 * 00006 * This is free software; you can redistribute it and/or modify 00007 * it under the terms of the GNU General Public License as published by 00008 * the Free Software Foundation; either version 2, or (at your option) 00009 * any later version. 00010 * 00011 * This software is distributed in the hope that it will be useful, 00012 * but WITHOUT ANY WARRANTY; without even the implied warranty of 00013 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the 00014 * GNU General Public License for more details. 00015 * 00016 * You should have received a copy of the GNU General Public License 00017 * along with this package; see the file COPYING. If not, write to 00018 * the Free Software Foundation, Inc., 51 Franklin Street - Fifth Floor, 00019 * Boston, MA 02110-1301, USA. 00020 * 00021 * $Id$ 00022 * 00023 */ 00082 #if HAVE_CONFIG_H 00083 # include <config.h> 00084 #endif 00085 00086 #include <assert.h> 00087 #include <stdio.h> 00088 #include <cstdlib> 00089 #include <string.h> 00090 #include <cmath> 00091 00092 #include "logging.h" 00093 #include "object.h" 00094 #include "complex.h" 00095 #include "vector.h" 00096 #include "matrix.h" 00097 00098 namespace qucs { 00099 00104 matrix::matrix () { 00105 rows = 0; 00106 cols = 0; 00107 data = NULL; 00108 } 00109 00117 matrix::matrix (int s) { 00118 rows = cols = s; 00119 data = (s > 0) ? new nr_complex_t[s * s] : NULL; 00120 } 00121 00122 /* \brief Creates a matrix 00123 00124 Constructor creates an unnamed instance of the matrix class with a 00125 certain number of rows and columns. 00126 \param[in] r number of rows 00127 \param[in] c number of column 00128 \todo Why not r and c constant 00129 \todo Assert r >= 0 and c >= 0 00130 */ 00131 matrix::matrix (int r, int c) { 00132 rows = r; 00133 cols = c; 00134 data = (r > 0 && c > 0) ? new nr_complex_t[r * c] : NULL; 00135 } 00136 00137 /* \brief copy constructor 00138 00139 The copy constructor creates a new instance based on the given 00140 matrix object. 00141 \todo Add assert tests 00142 */ 00143 matrix::matrix (const matrix & m) { 00144 rows = m.rows; 00145 cols = m.cols; 00146 data = NULL; 00147 00148 // copy matrix elements 00149 if (rows > 0 && cols > 0) { 00150 data = new nr_complex_t[rows * cols]; 00151 memcpy (data, m.data, sizeof (nr_complex_t) * rows * cols); 00152 } 00153 } 00154 00164 const matrix& matrix::operator=(const matrix & m) { 00165 if (&m != this) { 00166 rows = m.rows; 00167 cols = m.cols; 00168 if (data) { 00169 delete[] data; 00170 data = NULL; 00171 } 00172 if (rows > 0 && cols > 0) { 00173 data = new nr_complex_t[rows * cols]; 00174 memcpy (data, m.data, sizeof (nr_complex_t) * rows * cols); 00175 } 00176 } 00177 return *this; 00178 } 00179 00184 matrix::~matrix () { 00185 delete[] data; 00186 } 00187 00194 nr_complex_t matrix::get (int r, int c) { 00195 return data[r * cols + c]; 00196 } 00197 00205 void matrix::set (int r, int c, nr_complex_t z) { 00206 data[r * cols + c] = z; 00207 } 00208 00209 #ifdef DEBUG 00210 00211 void matrix::print (void) { 00212 for (int r = 0; r < rows; r++) { 00213 for (int c = 0; c < cols; c++) { 00214 fprintf (stderr, "%+.2e,%+.2e ", (double) real (get (r, c)), 00215 (double) imag (get (r, c))); 00216 } 00217 fprintf (stderr, "\n"); 00218 } 00219 } 00220 #endif /* DEBUG */ 00221 00228 matrix operator + (matrix a, matrix b) { 00229 assert (a.getRows () == b.getRows () && a.getCols () == b.getCols ()); 00230 00231 matrix res (a.getRows (), a.getCols ()); 00232 for (int r = 0; r < a.getRows (); r++) 00233 for (int c = 0; c < a.getCols (); c++) 00234 res.set (r, c, a.get (r, c) + b.get (r, c)); 00235 return res; 00236 } 00237 00243 matrix matrix::operator += (matrix a) { 00244 assert (a.getRows () == rows && a.getCols () == cols); 00245 00246 int r, c, i; 00247 for (i = 0, r = 0; r < a.getRows (); r++) 00248 for (c = 0; c < a.getCols (); c++, i++) 00249 data[i] += a.get (r, c); 00250 return *this; 00251 } 00252 00259 matrix operator - (matrix a, matrix b) { 00260 assert (a.getRows () == b.getRows () && a.getCols () == b.getCols ()); 00261 00262 matrix res (a.getRows (), a.getCols ()); 00263 for (int r = 0; r < a.getRows (); r++) 00264 for (int c = 0; c < a.getCols (); c++) 00265 res.set (r, c, a.get (r, c) - b.get (r, c)); 00266 return res; 00267 } 00268 00270 matrix matrix::operator - () { 00271 matrix res (getRows (), getCols ()); 00272 int r, c, i; 00273 for (i = 0, r = 0; r < getRows (); r++) 00274 for (c = 0; c < getCols (); c++, i++) 00275 res.set (r, c, -data[i]); 00276 return res; 00277 } 00278 00283 matrix matrix::operator -= (matrix a) { 00284 assert (a.getRows () == rows && a.getCols () == cols); 00285 int r, c, i; 00286 for (i = 0, r = 0; r < a.getRows (); r++) 00287 for (c = 0; c < a.getCols (); c++, i++) 00288 data[i] -= a.get (r, c); 00289 return *this; 00290 } 00291 00298 matrix operator * (matrix a, nr_complex_t z) { 00299 matrix res (a.getRows (), a.getCols ()); 00300 for (int r = 0; r < a.getRows (); r++) 00301 for (int c = 0; c < a.getCols (); c++) 00302 res.set (r, c, a.get (r, c) * z); 00303 return res; 00304 } 00305 00313 matrix operator * (nr_complex_t z, matrix a) { 00314 return a * z; 00315 } 00316 00323 matrix operator * (matrix a, nr_double_t d) { 00324 matrix res (a.getRows (), a.getCols ()); 00325 for (int r = 0; r < a.getRows (); r++) 00326 for (int c = 0; c < a.getCols (); c++) 00327 res.set (r, c, a.get (r, c) * d); 00328 return res; 00329 } 00330 00338 matrix operator * (nr_double_t d, matrix a) { 00339 return a * d; 00340 } 00341 00348 matrix operator / (matrix a, nr_complex_t z) { 00349 matrix res (a.getRows (), a.getCols ()); 00350 for (int r = 0; r < a.getRows (); r++) 00351 for (int c = 0; c < a.getCols (); c++) 00352 res.set (r, c, a.get (r, c) / z); 00353 return res; 00354 } 00355 00362 matrix operator / (matrix a, nr_double_t d) { 00363 matrix res (a.getRows (), a.getCols ()); 00364 for (int r = 0; r < a.getRows (); r++) 00365 for (int c = 0; c < a.getCols (); c++) 00366 res.set (r, c, a.get (r, c) / d); 00367 return res; 00368 } 00369 00378 matrix operator * (matrix a, matrix b) { 00379 assert (a.getCols () == b.getRows ()); 00380 00381 int r, c, i, n = a.getCols (); 00382 nr_complex_t z; 00383 matrix res (a.getRows (), b.getCols ()); 00384 for (r = 0; r < a.getRows (); r++) { 00385 for (c = 0; c < b.getCols (); c++) { 00386 for (i = 0, z = 0; i < n; i++) z += a.get (r, i) * b.get (i, c); 00387 res.set (r, c, z); 00388 } 00389 } 00390 return res; 00391 } 00392 00399 matrix operator + (matrix a, nr_complex_t z) { 00400 matrix res (a.getRows (), a.getCols ()); 00401 for (int r = 0; r < a.getRows (); r++) 00402 for (int c = 0; c < a.getCols (); c++) 00403 res.set (r, c, a.get (r, c) + z); 00404 return res; 00405 } 00406 00414 matrix operator + (nr_complex_t z, matrix a) { 00415 return a + z; 00416 } 00417 00424 matrix operator + (matrix a, nr_double_t d) { 00425 matrix res (a.getRows (), a.getCols ()); 00426 for (int r = 0; r < a.getRows (); r++) 00427 for (int c = 0; c < a.getCols (); c++) 00428 res.set (r, c, a.get (r, c) + d); 00429 return res; 00430 } 00431 00439 matrix operator + (nr_double_t d, matrix a) { 00440 return a + d; 00441 } 00442 00450 matrix operator - (matrix a, nr_complex_t z) { 00451 return -z + a; 00452 } 00453 00461 matrix operator - (nr_complex_t z, matrix a) { 00462 return -a + z; 00463 } 00464 00472 matrix operator - (matrix a, nr_double_t d) { 00473 return -d + a; 00474 } 00475 00483 matrix operator - (nr_double_t d, matrix a) { 00484 return -a + d; 00485 } 00486 00492 matrix transpose (matrix a) { 00493 matrix res (a.getCols (), a.getRows ()); 00494 for (int r = 0; r < a.getRows (); r++) 00495 for (int c = 0; c < a.getCols (); c++) 00496 res.set (c, r, a.get (r, c)); 00497 return res; 00498 } 00499 00505 matrix conj (matrix a) { 00506 matrix res (a.getRows (), a.getCols ()); 00507 for (int r = 0; r < a.getRows (); r++) 00508 for (int c = 0; c < a.getCols (); c++) 00509 res.set (r, c, conj (a.get (r, c))); 00510 return res; 00511 } 00512 00522 matrix adjoint (matrix a) { 00523 return transpose (conj (a)); 00524 } 00525 00531 matrix abs (matrix a) { 00532 matrix res (a.getRows (), a.getCols ()); 00533 for (int r = 0; r < a.getRows (); r++) 00534 for (int c = 0; c < a.getCols (); c++) 00535 res.set (r, c, abs (a.get (r, c))); 00536 return res; 00537 } 00538 00542 matrix dB (matrix a) { 00543 matrix res (a.getRows (), a.getCols ()); 00544 for (int r = 0; r < a.getRows (); r++) 00545 for (int c = 0; c < a.getCols (); c++) 00546 res.set (r, c, dB (a.get (r, c))); 00547 return res; 00548 } 00549 00555 matrix arg (matrix a) { 00556 matrix res (a.getRows (), a.getCols ()); 00557 for (int r = 0; r < a.getRows (); r++) 00558 for (int c = 0; c < a.getCols (); c++) 00559 res.set (r, c, arg (a.get (r, c))); 00560 return res; 00561 } 00562 00568 matrix real (matrix a) { 00569 matrix res (a.getRows (), a.getCols ()); 00570 for (int r = 0; r < a.getRows (); r++) 00571 for (int c = 0; c < a.getCols (); c++) 00572 res.set (r, c, real (a.get (r, c))); 00573 return res; 00574 } 00575 00581 matrix imag (matrix a) { 00582 matrix res (a.getRows (), a.getCols ()); 00583 for (int r = 0; r < a.getRows (); r++) 00584 for (int c = 0; c < a.getCols (); c++) 00585 res.set (r, c, imag (a.get (r, c))); 00586 return res; 00587 } 00588 00592 matrix sqr (matrix a) { 00593 return a * a; 00594 } 00595 00603 matrix eye (int rs, int cs) { 00604 matrix res (rs, cs); 00605 for (int r = 0; r < res.getRows (); r++) 00606 for (int c = 0; c < res.getCols (); c++) 00607 if (r == c) res.set (r, c, 1); 00608 return res; 00609 } 00610 00616 matrix eye (int s) { 00617 return eye (s, s); 00618 } 00619 00624 matrix diagonal (qucs::vector diag) { 00625 int size = diag.getSize (); 00626 matrix res (size); 00627 for (int i = 0; i < size; i++) res (i, i) = diag (i); 00628 return res; 00629 } 00630 00631 // Compute n-th power of the given matrix. 00632 matrix pow (matrix a, int n) { 00633 matrix res; 00634 if (n == 0) { 00635 res = eye (a.getRows (), a.getCols ()); 00636 } 00637 else { 00638 res = a = n < 0 ? inverse (a) : a; 00639 for (int i = 1; i < std::abs (n); i++) 00640 res = res * a; 00641 } 00642 return res; 00643 } 00644 00657 nr_complex_t cofactor (matrix a, int u, int v) { 00658 matrix res (a.getRows () - 1, a.getCols () - 1); 00659 int r, c, ra, ca; 00660 for (ra = r = 0; r < res.getRows (); r++, ra++) { 00661 if (ra == u) ra++; 00662 for (ca = c = 0; c < res.getCols (); c++, ca++) { 00663 if (ca == v) ca++; 00664 res.set (r, c, a.get (ra, ca)); 00665 } 00666 } 00667 nr_complex_t z = detLaplace (res); 00668 return ((u + v) & 1) ? -z : z; 00669 } 00670 00686 nr_complex_t detLaplace (matrix a) { 00687 assert (a.getRows () == a.getCols ()); 00688 int s = a.getRows (); 00689 nr_complex_t res = 0; 00690 if (s > 1) { 00691 /* always use the first row for sub-determinant, but you should 00692 use the row or column with most zeros in it */ 00693 int r = 0; 00694 for (int i = 0; i < s; i++) { 00695 res += a.get (r, i) * cofactor (a, r, i); 00696 } 00697 return res; 00698 } 00699 /* 1 by 1 matrix */ 00700 else if (s == 1) { 00701 return a (0, 0); 00702 } 00703 /* 0 by 0 matrix */ 00704 return 1; 00705 } 00706 00717 nr_complex_t detGauss (matrix a) { 00718 assert (a.getRows () == a.getCols ()); 00719 nr_double_t MaxPivot; 00720 nr_complex_t f, res; 00721 matrix b; 00722 int i, c, r, pivot, n = a.getCols (); 00723 00724 // return special matrix cases 00725 if (n == 0) return 1; 00726 if (n == 1) return a (0, 0); 00727 00728 // make copy of original matrix 00729 b = matrix (a); 00730 00731 // triangulate the matrix 00732 for (res = 1, i = 0; i < n; i++) { 00733 // find maximum column value for pivoting 00734 for (MaxPivot = 0, pivot = r = i; r < n; r++) { 00735 if (abs (b.get (r, i)) > MaxPivot) { 00736 MaxPivot = abs (b.get (r, i)); 00737 pivot = r; 00738 } 00739 } 00740 // exchange rows if necessary 00741 assert (MaxPivot != 0); 00742 if (i != pivot) { b.exchangeRows (i, pivot); res = -res; } 00743 // compute new rows and columns 00744 for (r = i + 1; r < n; r++) { 00745 f = b.get (r, i) / b.get (i, i); 00746 for (c = i + 1; c < n; c++) { 00747 b.set (r, c, b.get (r, c) - f * b.get (i, c)); 00748 } 00749 } 00750 } 00751 00752 // now compute determinant by multiplying diagonal 00753 for (i = 0; i < n; i++) res *= b.get (i, i); 00754 return res; 00755 } 00756 00762 nr_complex_t det (matrix a) { 00763 #if 0 00764 return detLaplace (a); 00765 #else 00766 return detGauss (a); 00767 #endif 00768 } 00769 00779 matrix inverseLaplace (matrix a) { 00780 matrix res (a.getRows (), a.getCols ()); 00781 nr_complex_t d = detLaplace (a); 00782 assert (abs (d) != 0); // singular matrix 00783 for (int r = 0; r < a.getRows (); r++) 00784 for (int c = 0; c < a.getCols (); c++) 00785 res.set (r, c, cofactor (a, c, r) / d); 00786 return res; 00787 } 00788 00798 matrix inverseGaussJordan (matrix a) { 00799 nr_double_t MaxPivot; 00800 nr_complex_t f; 00801 matrix b, e; 00802 int i, c, r, pivot, n = a.getCols (); 00803 00804 // create temporary matrix and the result matrix 00805 b = matrix (a); 00806 e = eye (n); 00807 00808 // create the eye matrix in 'b' and the result in 'e' 00809 for (i = 0; i < n; i++) { 00810 // find maximum column value for pivoting 00811 for (MaxPivot = 0, pivot = r = i; r < n; r++) { 00812 if (abs (b (r, i)) > MaxPivot) { 00813 MaxPivot = abs (b (r, i)); 00814 pivot = r; 00815 } 00816 } 00817 // exchange rows if necessary 00818 assert (MaxPivot != 0); // singular matrix 00819 if (i != pivot) { 00820 b.exchangeRows (i, pivot); 00821 e.exchangeRows (i, pivot); 00822 } 00823 00824 // compute current row 00825 for (f = b (i, i), c = 0; c < n; c++) { 00826 b (i, c) /= f; 00827 e (i, c) /= f; 00828 } 00829 00830 // compute new rows and columns 00831 for (r = 0; r < n; r++) { 00832 if (r != i) { 00833 for (f = b (r, i), c = 0; c < n; c++) { 00834 b (r, c) -= f * b (i, c); 00835 e (r, c) -= f * e (i, c); 00836 } 00837 } 00838 } 00839 } 00840 return e; 00841 } 00842 00847 matrix inverse (matrix a) { 00848 #if 0 00849 return inverseLaplace (a); 00850 #else 00851 return inverseGaussJordan (a); 00852 #endif 00853 } 00854 00890 matrix stos (matrix s, qucs::vector zref, qucs::vector z0) { 00891 int d = s.getRows (); 00892 matrix e, r, a; 00893 00894 assert (d == s.getCols () && d == z0.getSize () && d == zref.getSize ()); 00895 00896 e = eye (d); 00897 r = diagonal ((z0 - zref) / (z0 + zref)); 00898 a = diagonal (sqrt (z0 / zref) / (z0 + zref)); 00899 return inverse (a) * (s - r) * inverse (e - r * s) * a; 00900 } 00901 00910 matrix stos (matrix s, nr_complex_t zref, nr_complex_t z0) { 00911 int d = s.getRows (); 00912 return stos (s, qucs::vector (d, zref), qucs::vector (d, z0)); 00913 } 00914 00923 matrix stos (matrix s, nr_double_t zref, nr_double_t z0) { 00924 return stos (s, nr_complex_t (zref, 0), nr_complex_t (z0, 0)); 00925 } 00926 00935 matrix stos (matrix s, qucs::vector zref, nr_complex_t z0) { 00936 return stos (s, zref, qucs::vector (zref.getSize (), z0)); 00937 } 00938 00947 matrix stos (matrix s, nr_complex_t zref, qucs::vector z0) { 00948 return stos (s, qucs::vector (z0.getSize (), zref), z0); 00949 } 00950 00973 matrix stoz (matrix s, qucs::vector z0) { 00974 int d = s.getRows (); 00975 matrix e, zref, gref; 00976 00977 assert (d == s.getCols () && d == z0.getSize ()); 00978 00979 e = eye (d); 00980 zref = diagonal (z0); 00981 gref = diagonal (sqrt (real (1 / z0))); 00982 return inverse (gref) * inverse (e - s) * (s * zref + zref) * gref; 00983 } 00984 00992 matrix stoz (matrix s, nr_complex_t z0) { 00993 return stoz (s, qucs::vector (s.getRows (), z0)); 00994 } 00995 01018 matrix ztos (matrix z, qucs::vector z0) { 01019 int d = z.getRows (); 01020 matrix e, zref, gref; 01021 01022 assert (d == z.getCols () && d == z0.getSize ()); 01023 01024 e = eye (d); 01025 zref = diagonal (z0); 01026 gref = diagonal (sqrt (real (1 / z0))); 01027 return gref * (z - zref) * inverse (z + zref) * inverse (gref); 01028 } 01029 01037 matrix ztos (matrix z, nr_complex_t z0) { 01038 return ztos (z, qucs::vector (z.getRows (), z0)); 01039 } 01040 01050 matrix ztoy (matrix z) { 01051 assert (z.getRows () == z.getCols ()); 01052 return inverse (z); 01053 } 01054 01082 matrix stoy (matrix s, qucs::vector z0) { 01083 int d = s.getRows (); 01084 matrix e, zref, gref; 01085 01086 assert (d == s.getCols () && d == z0.getSize ()); 01087 01088 e = eye (d); 01089 zref = diagonal (z0); 01090 gref = diagonal (sqrt (real (1 / z0))); 01091 return inverse (gref) * inverse (s * zref + zref) * (e - s) * gref; 01092 } 01093 01101 matrix stoy (matrix s, nr_complex_t z0) { 01102 return stoy (s, qucs::vector (s.getRows (), z0)); 01103 } 01104 01133 matrix ytos (matrix y, qucs::vector z0) { 01134 int d = y.getRows (); 01135 matrix e, zref, gref; 01136 01137 assert (d == y.getCols () && d == z0.getSize ()); 01138 01139 e = eye (d); 01140 zref = diagonal (z0); 01141 gref = diagonal (sqrt (real (1 / z0))); 01142 return gref * (e - zref * y) * inverse (e + zref * y) * inverse (gref); 01143 } 01151 matrix ytos (matrix y, nr_complex_t z0) { 01152 return ytos (y, qucs::vector (y.getRows (), z0)); 01153 } 01181 matrix stoa (matrix s, nr_complex_t z1, nr_complex_t z2) { 01182 nr_complex_t d = s (0, 0) * s (1, 1) - s (0, 1) * s (1, 0); 01183 nr_complex_t n = 2.0 * s (1, 0) * sqrt (fabs (real (z1) * real (z2))); 01184 matrix a (2); 01185 01186 assert (s.getRows () >= 2 && s.getCols () >= 2); 01187 01188 a.set (0, 0, (conj (z1) + z1 * s (0, 0) - 01189 conj (z1) * s (1, 1) - z1 * d) / n); 01190 a.set (0, 1, (conj (z1) * conj (z2) + z1 * conj (z2) * s (0, 0) + 01191 conj (z1) * z2 * s (1, 1) + z1 * z2 * d) / n); 01192 a.set (1, 0, (1.0 - s (0, 0) - s (1, 1) + d) / n); 01193 a.set (1, 1, (conj (z2) - conj (z2) * s (0, 0) + 01194 z2 * s (1, 1) - z2 * d) / n); 01195 return a; 01196 } 01197 01198 01222 matrix atos (matrix a, nr_complex_t z1, nr_complex_t z2) { 01223 nr_complex_t d = 2.0 * sqrt (fabs (real (z1) * real (z2))); 01224 nr_complex_t n = a (0, 0) * z2 + a (0, 1) + 01225 a (1, 0) * z1 * z2 + a (1, 1) * z1; 01226 matrix s (2); 01227 01228 assert (a.getRows () >= 2 && a.getCols () >= 2); 01229 01230 s.set (0, 0, (a (0, 0) * z2 + a (0, 1) 01231 - a (1, 0) * conj (z1) * z2 - a (1, 1) * conj (z1)) / n); 01232 s.set (0, 1, (a (0, 0) * a (1, 1) - 01233 a (0, 1) * a (1, 0)) * d / n); 01234 s.set (1, 0, d / n); 01235 s.set (1, 1, (a (1, 1) * z1 - a (0, 0) * conj (z2) + 01236 a (0, 1) - a (1, 0) * z1 * conj (z2)) / n); 01237 return s; 01238 } 01239 01269 matrix stoh (matrix s, nr_complex_t z1, nr_complex_t z2) { 01270 nr_complex_t n = s (0, 1) * s (1, 0); 01271 nr_complex_t d = (1.0 - s (0, 0)) * (1.0 + s (1, 1)) + n; 01272 matrix h (2); 01273 01274 assert (s.getRows () >= 2 && s.getCols () >= 2); 01275 01276 h.set (0, 0, ((1.0 + s (0, 0)) * (1.0 + s (1, 1)) - n) * z1 / d); 01277 h.set (0, 1, +2.0 * s (0, 1) / d); 01278 h.set (1, 0, -2.0 * s (1, 0) / d); 01279 h.set (1, 1, ((1.0 - s (0, 0)) * (1.0 - s (1, 1)) - n) / z2 / d); 01280 return h; 01281 } 01282 01307 matrix htos (matrix h, nr_complex_t z1, nr_complex_t z2) { 01308 nr_complex_t n = h (0, 1) * h (1, 0); 01309 nr_complex_t d = (1.0 + h (0, 0) / z1) * (1.0 + z2 * h (1, 1)) - n; 01310 matrix s (2); 01311 01312 assert (h.getRows () >= 2 && h.getCols () >= 2); 01313 01314 s.set (0, 0, ((h (0, 0) / z1 - 1.0) * (1.0 + z2 * h (1, 1)) - n) / d); 01315 s.set (0, 1, +2.0 * h (0, 1) / d); 01316 s.set (1, 0, -2.0 * h (1, 0) / d); 01317 s.set (1, 1, ((1.0 + h (0, 0) / z1) * (1.0 - z2 * h (1, 1)) + n) / d); 01318 return s; 01319 } 01320 01321 /*\brief Converts scattering parameters to second hybrid matrix. 01322 \bug{Programmed formulae are valid only for Z real} 01323 \bug{Not documented and references} 01324 \param[in] s Scattering matrix 01325 \param[in] z1 impedance at input 1 01326 \param[in] z2 impedance at input 2 01327 \return second hybrid matrix 01328 \note Assert 2 by 2 matrix 01329 \todo Why not s,z1,z2 const 01330 */ 01331 matrix stog (matrix s, nr_complex_t z1, nr_complex_t z2) { 01332 nr_complex_t n = s (0, 1) * s (1, 0); 01333 nr_complex_t d = (1.0 + s (0, 0)) * (1.0 - s (1, 1)) + n; 01334 matrix g (2); 01335 01336 assert (s.getRows () >= 2 && s.getCols () >= 2); 01337 01338 g.set (0, 0, ((1.0 - s (0, 0)) * (1.0 - s (1, 1)) - n) / z1 / d); 01339 g.set (0, 1, -2.0 * s (0, 1) / d); 01340 g.set (1, 0, +2.0 * s (1, 0) / d); 01341 g.set (1, 1, ((1.0 + s (0, 0)) * (1.0 + s (1, 1)) - n) * z2 / d); 01342 return g; 01343 } 01344 01345 /*\brief Converts second hybrid matrix to scattering parameters. 01346 \bug{Programmed formulae are valid only for Z real} 01347 \bug{Not documented and references} 01348 \param[in] g second hybrid matrix 01349 \param[in] z1 impedance at input 1 01350 \param[in] z2 impedance at input 2 01351 \return scattering matrix 01352 \note Assert 2 by 2 matrix 01353 \todo Why not g,z1,z2 const 01354 */ 01355 matrix gtos (matrix g, nr_complex_t z1, nr_complex_t z2) { 01356 nr_complex_t n = g (0, 1) * g (1, 0); 01357 nr_complex_t d = (1.0 + g (0, 0) * z1) * (1.0 + g (1, 1) / z2) - n; 01358 matrix s (2); 01359 01360 assert (g.getRows () >= 2 && g.getCols () >= 2); 01361 01362 s.set (0, 0, ((1.0 - g (0, 0) * z1) * (1.0 + g (1, 1) / z2) + n) / d); 01363 s.set (0, 1, -2.0 * g (0, 1) / d); 01364 s.set (1, 0, +2.0 * g (1, 0) / d); 01365 s.set (1, 1, ((g (0, 0) * z1 + 1.0) * (g (1, 1) / z2 - 1.0) - n) / d); 01366 return s; 01367 } 01368 01380 matrix ytoz (matrix y) { 01381 assert (y.getRows () == y.getCols ()); 01382 return inverse (y); 01383 } 01384 01404 matrix cytocs (matrix cy, matrix s) { 01405 matrix e = eye (s.getRows ()); 01406 01407 assert (cy.getRows () == cy.getCols () && s.getRows () == s.getCols () && 01408 cy.getRows () == s.getRows ()); 01409 01410 return (e + s) * cy * adjoint (e + s) / 4; 01411 } 01412 01430 matrix cstocy (matrix cs, matrix y) { 01431 matrix e = eye (y.getRows ()); 01432 01433 assert (cs.getRows () == cs.getCols () && y.getRows () == y.getCols () && 01434 cs.getRows () == y.getRows ()); 01435 01436 return (e + y) * cs * adjoint (e + y); 01437 } 01438 01457 matrix cztocs (matrix cz, matrix s) { 01458 matrix e = eye (s.getRows ()); 01459 01460 assert (cz.getRows () == cz.getCols () && s.getRows () == s.getCols () && 01461 cz.getRows () == s.getRows ()); 01462 01463 return (e - s) * cz * adjoint (e - s) / 4; 01464 } 01465 01483 matrix cstocz (matrix cs, matrix z) { 01484 assert (cs.getRows () == cs.getCols () && z.getRows () == z.getCols () && 01485 cs.getRows () == z.getRows ()); 01486 matrix e = eye (z.getRows ()); 01487 return (e + z) * cs * adjoint (e + z); 01488 } 01489 01507 matrix cztocy (matrix cz, matrix y) { 01508 assert (cz.getRows () == cz.getCols () && y.getRows () == y.getCols () && 01509 cz.getRows () == y.getRows ()); 01510 01511 return y * cz * adjoint (y); 01512 } 01513 01531 matrix cytocz (matrix cy, matrix z) { 01532 assert (cy.getRows () == cy.getCols () && z.getRows () == z.getCols () && 01533 cy.getRows () == z.getRows ()); 01534 return z * cy * adjoint (z); 01535 } 01536 01543 void matrix::exchangeRows (int r1, int r2) { 01544 nr_complex_t * s = new nr_complex_t[cols]; 01545 int len = sizeof (nr_complex_t) * cols; 01546 01547 assert (r1 >= 0 && r2 >= 0 && r1 < rows && r2 < rows); 01548 01549 memcpy (s, &data[r1 * cols], len); 01550 memcpy (&data[r1 * cols], &data[r2 * cols], len); 01551 memcpy (&data[r2 * cols], s, len); 01552 delete[] s; 01553 } 01554 01561 void matrix::exchangeCols (int c1, int c2) { 01562 nr_complex_t s; 01563 01564 assert (c1 >= 0 && c2 >= 0 && c1 < cols && c2 < cols); 01565 01566 for (int r = 0; r < rows * cols; r += cols) { 01567 s = data[r + c1]; 01568 data[r + c1] = data[r + c2]; 01569 data[r + c2] = s; 01570 } 01571 } 01572 01594 matrix twoport (matrix m, char in, char out) { 01595 assert (m.getRows () >= 2 && m.getCols () >= 2); 01596 nr_complex_t d; 01597 matrix res (2); 01598 01599 switch (in) { 01600 case 'Y': 01601 switch (out) { 01602 case 'Y': // Y to Y 01603 res = m; 01604 break; 01605 case 'Z': // Y to Z 01606 d = m (0, 0) * m (1, 1) - m (0, 1) * m (1, 0); 01607 res.set (0, 0, m (1, 1) / d); 01608 res.set (0, 1, -m (0, 1) / d); 01609 res.set (1, 0, -m (1, 0) / d); 01610 res.set (1, 1, m (0, 0) / d); 01611 break; 01612 case 'H': // Y to H 01613 d = m (0, 0); 01614 res.set (0, 0, 1.0 / d); 01615 res.set (0, 1, -m (0, 1) / d); 01616 res.set (1, 0, m (1, 0) / d); 01617 res.set (1, 1, m (1, 1) - m (0, 1) * m (1, 0) / d); 01618 break; 01619 case 'G': // Y to G 01620 d = m (1, 1); 01621 res.set (0, 0, m (0, 0) - m (0, 1) * m (1, 0) / d); 01622 res.set (0, 1, m (0, 1) / d); 01623 res.set (1, 0, -m (1, 0) / d); 01624 res.set (1, 1, 1.0 / d); 01625 break; 01626 case 'A': // Y to A 01627 d = m (1, 0); 01628 res.set (0, 0, -m (1, 1) / d); 01629 res.set (0, 1, -1.0 / d); 01630 res.set (1, 0, m (0, 1) - m (1, 1) * m (0, 0) / d); 01631 res.set (1, 1, -m (0, 0) / d); 01632 break; 01633 case 'S': // Y to S 01634 res = ytos (m); 01635 break; 01636 } 01637 break; 01638 case 'Z': 01639 switch (out) { 01640 case 'Y': // Z to Y 01641 d = m (0, 0) * m (1, 1) - m (0, 1) * m (1, 0); 01642 res.set (0, 0, m (1, 1) / d); 01643 res.set (0, 1, -m (0, 1) / d); 01644 res.set (1, 0, -m (1, 0) / d); 01645 res.set (1, 1, m (0, 0) / d); 01646 break; 01647 case 'Z': // Z to Z 01648 res = m; 01649 break; 01650 case 'H': // Z to H 01651 d = m (1, 1); 01652 res.set (0, 0, m (0, 0) - m (0, 1) * m (1, 0) / d); 01653 res.set (0, 1, m (0, 1) / d); 01654 res.set (1, 0, -m (1, 0) / d); 01655 res.set (1, 1, 1.0 / d); 01656 break; 01657 case 'G': // Z to G 01658 d = m (0, 0); 01659 res.set (0, 0, 1.0 / d); 01660 res.set (0, 1, -m (0, 1) / d); 01661 res.set (1, 0, m (1, 0) / d); 01662 res.set (1, 1, m (1, 1) - m (0, 1) * m (1, 0) / d); 01663 break; 01664 case 'A': // Z to A 01665 d = m (1, 0); 01666 res.set (0, 0, m (0, 0) / d); 01667 res.set (0, 1, m (0, 0) * m (1, 1) / d - m (0, 1)); 01668 res.set (1, 0, 1.0 / d); 01669 res.set (1, 1, m (1, 1) / d); 01670 break; 01671 case 'S': // Z to S 01672 res = ztos (m); 01673 break; 01674 } 01675 break; 01676 case 'H': 01677 switch (out) { 01678 case 'Y': // H to Y 01679 d = m (0, 0); 01680 res.set (0, 0, 1.0 / d); 01681 res.set (0, 1, -m (0, 1) / d); 01682 res.set (1, 0, m (1, 0) / d); 01683 res.set (1, 1, m (1, 1) - m (0, 1) * m.get(2, 1) / d); 01684 break; 01685 case 'Z': // H to Z 01686 d = m (1, 1); 01687 res.set (0, 0, m (0, 0) - m (0, 1) * m (1, 0) / d); 01688 res.set (0, 1, m (0, 1) / d); 01689 res.set (1, 0, -m (1, 0) / d); 01690 res.set (1, 1, 1.0 / d); 01691 break; 01692 case 'H': // H to H 01693 res = m; 01694 break; 01695 case 'G': // H to G 01696 d = m (0, 0) * m (1, 1) - m (0, 1) * m (1, 0); 01697 res.set (0, 0, m (1, 1) / d); 01698 res.set (0, 1, -m (0, 1) / d); 01699 res.set (1, 0, -m (1, 0) / d); 01700 res.set (1, 1, m (0, 0) / d); 01701 break; 01702 case 'A': // H to A 01703 d = m (1, 0); 01704 res.set (0, 0, m (0, 1) - m (0, 0) * m (1, 1) / d); 01705 res.set (0, 1, -m (0, 0) / d); 01706 res.set (1, 0, -m (1, 1) / d); 01707 res.set (1, 1, -1.0 / d); 01708 break; 01709 case 'S': // H to S 01710 res = htos (m); 01711 break; 01712 } 01713 break; 01714 case 'G': 01715 switch (out) { 01716 case 'Y': // G to Y 01717 d = m (1, 1); 01718 res.set (0, 0, m (0, 0) - m (0, 1) * m (1, 0) / d); 01719 res.set (0, 1, m (0, 1) / d); 01720 res.set (1, 0, -m (1, 0) / d); 01721 res.set (1, 1, 1.0 / d); 01722 break; 01723 case 'Z': // G to Z 01724 d = m (0, 0); 01725 res.set (0, 0, 1.0 / d); 01726 res.set (0, 1, -m (0, 1) / d); 01727 res.set (1, 0, m (1, 0) / d); 01728 res.set (1, 1, m (1, 1) - m (0, 1) * m (1, 0) / d); 01729 break; 01730 case 'H': // G to H 01731 d = m (0, 0) * m (1, 1) - m (0, 1) * m (1, 0); 01732 res.set (0, 0, m (1, 1) / d); 01733 res.set (0, 1, -m (0, 1) / d); 01734 res.set (1, 0, -m (1, 0) / d); 01735 res.set (1, 1, m (0, 0) / d); 01736 break; 01737 case 'G': // G to G 01738 res = m; 01739 break; 01740 case 'A': // G to A 01741 d = m (1, 0); 01742 res.set (0, 0, 1.0 / d); 01743 res.set (0, 1, m (1, 1) / d); 01744 res.set (1, 0, m (0, 0) / d); 01745 res.set (1, 1, m (0, 0) * m (1, 1) / d - m (0, 1)); 01746 break; 01747 case 'S': // G to S 01748 res = gtos (m); 01749 break; 01750 } 01751 break; 01752 case 'A': 01753 switch (out) { 01754 case 'Y': // A to Y 01755 d = m (0, 1); 01756 res.set (0, 0, m (1, 1) / d); 01757 res.set (0, 1, m (1, 0) - m (0, 0) * m (1, 1) / d); 01758 res.set (1, 0, -1.0 / d); 01759 res.set (1, 1, m (0, 0) / d); 01760 break; 01761 case 'Z': // A to Z 01762 d = m (1, 0); 01763 res.set (0, 0, m (0, 0) / d); 01764 res.set (0, 1, m (0, 0) * m (1, 1) / d - m (0, 1)); 01765 res.set (1, 0, 1.0 / d); 01766 res.set (1, 1, m (1, 1) / d); 01767 break; 01768 case 'H': // A to H 01769 d = m (1, 1); 01770 res.set (0, 0, m (0, 1) / d); 01771 res.set (0, 1, m (0, 0) - m (0, 1) * m (1, 0) / d); 01772 res.set (1, 0, -1.0 / d); 01773 res.set (1, 1, m (1, 0) / d); 01774 break; 01775 case 'G': // A to G 01776 d = m (0, 0); 01777 res.set (0, 0, m (1, 0) / d); 01778 res.set (0, 1, m (1, 0) * m (0, 1) / d - m (1, 1)); 01779 res.set (1, 0, 1.0 / d); 01780 res.set (1, 1, m (0, 1) / d); 01781 break; 01782 case 'A': // A to A 01783 res = m; 01784 break; 01785 case 'S': // A to S 01786 res = atos (m); 01787 break; 01788 } 01789 break; 01790 case 'S': 01791 switch (out) { 01792 case 'S': // S to S 01793 res = m; 01794 break; 01795 case 'T': // S to T 01796 d = m (1, 0); 01797 res.set (0, 0, m (0, 1) - m (0, 0) * m (1, 1) / d); 01798 res.set (0, 1, m (0, 0) / d); 01799 res.set (1, 0, -m (1, 1) / d); 01800 res.set (0, 1, 1.0 / d); 01801 break; 01802 case 'A': // S to A 01803 res = stoa (m); 01804 break; 01805 case 'H': // S to H 01806 res = stoh (m); 01807 break; 01808 case 'G': // S to G 01809 res = stog (m); 01810 break; 01811 case 'Y': // S to Y 01812 res = stoy (m); 01813 break; 01814 case 'Z': // S to Z 01815 res = stoz (m); 01816 break; 01817 } 01818 break; 01819 case 'T': 01820 switch (out) { 01821 case 'S': // T to S 01822 d = m (1, 1); 01823 res.set (0, 0, m (0, 1) / d); 01824 res.set (0, 1, m (0, 0) - m (0, 1) * m (1, 0) / d); 01825 res.set (1, 0, 1.0 / d); 01826 res.set (0, 1, -m (1, 0) / d); 01827 break; 01828 case 'T': // T to T 01829 res = m; 01830 break; 01831 } 01832 break; 01833 } 01834 return res; 01835 } 01836 01854 nr_double_t rollet (matrix m) { 01855 assert (m.getRows () >= 2 && m.getCols () >= 2); 01856 nr_double_t res; 01857 res = (1 - norm (m (0, 0)) - norm (m (1, 1)) + norm (det (m))) / 01858 2 / abs (m (0, 1) * m (1, 0)); 01859 return res; 01860 } 01861 01862 /* Computes stability measure B1 of the given S-parameter matrix. */ 01863 nr_double_t b1 (matrix m) { 01864 assert (m.getRows () >= 2 && m.getCols () >= 2); 01865 nr_double_t res; 01866 res = 1 + norm (m (0, 0)) - norm (m (1, 1)) - norm (det (m)); 01867 return res; 01868 } 01869 01870 01871 matrix rad2deg (matrix a) { 01872 matrix res (a.getRows (), a.getCols ()); 01873 for (int r = 0; r < a.getRows (); r++) 01874 for (int c = 0; c < a.getCols (); c++) 01875 res.set (r, c, rad2deg (a.get (r, c))); 01876 return res; 01877 } 01878 01879 matrix deg2rad (matrix a) { 01880 matrix res (a.getRows (), a.getCols ()); 01881 for (int r = 0; r < a.getRows (); r++) 01882 for (int c = 0; c < a.getCols (); c++) 01883 res.set (r, c, deg2rad (a.get (r, c))); 01884 return res; 01885 } 01886 01887 } // namespace qucs