octave_packages/control-2.3.52/__dss_bilin__.m

   1 ## Copyright (C) 2011   Lukas F. Reichlin
   2 ##
   3 ## This file is part of LTI Syncope.
   4 ##
   5 ## LTI Syncope is free software: you can redistribute it and/or modify
   6 ## it under the terms of the GNU General Public License as published by
   7 ## the Free Software Foundation, either version 3 of the License, or
   8 ## (at your option) any later version.
   9 ##
  10 ## LTI Syncope is distributed in the hope that it will be useful,
  11 ## but WITHOUT ANY WARRANTY; without even the implied warranty of
  12 ## MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
  13 ## GNU General Public License for more details.
  14 ##
  15 ## You should have received a copy of the GNU General Public License
  16 ## along with LTI Syncope.  If not, see <http://www.gnu.org/licenses/>.
  17
  18 ## -*- texinfo -*-
  19 ## 1.  Discrete -> continuous
  20 ## _
  21 ## E = alpha*E + A
  22 ## _
  23 ## A = beta (A - alpha*E)
  24 ## _
  25 ## B = sqrt(2*alpha*beta) * B
  26 ## _                                         -1
  27 ## C = sqrt(2*alpha*beta) * C * (alpha*E + A)  * E
  28 ## _                        -1
  29 ## D = D - C * (alpha*E + A)  * B
  30 ##
  31 ##
  32 ## 2.  Continuous -> discrete
  33 ## _
  34 ## E = beta*E - A
  35 ## _
  36 ## A = alpha (beta*E + A)
  37 ## _
  38 ## B = sqrt(2*alpha*beta) * B
  39 ## _                                        -1
  40 ## C = sqrt(2*alpha*beta) * C * (beta*E - A)  * E
  41 ## _                       -1
  42 ## D = D + C * (beta*E - A)  * B
  43
  44 ## Special thanks to Andras Varga for the formulae.
  45 ## Author: Lukas Reichlin <lukas.reichlin@gmail.com>
  46 ## Created: October 2011
  47 ## Version: 0.1
  48
  49 function [Ar, Br, Cr, Dr, Er] = __dss_bilin__ (A, B, C, D, E, beta, discrete)
  50
  51   if (discrete)
  52     EpA = E + A;
  53     s2b = sqrt (2*beta);
  54     if (rcond (EpA) < eps)
  55       error ("d2c: E+A singular");
  56     endif
  57     CiEpA = C / EpA;
  58
  59     Er = EpA;
  60     Ar = beta * (A - E);
  61     Br = s2b * B;
  62     Cr = s2b * CiEpA * E;
  63     Dr = D - CiEpA * B;
  64
  65     ## Er = E + A;
  66     ## Ar = beta * (A - E);
  67     ## Br = sqrt (2*beta) * B;
  68     ## Cr = sqrt (2*beta) * C / (E + A) * E;
  69     ## Dr = D - C / (E + A) * B;
  70   else
  71     bEmA = beta*E - A;
  72     s2b = sqrt (2*beta);
  73     if (rcond (bEmA) < eps)
  74       error ("c2d: beta*E-A singular");
  75     endif
  76     CibEmA = C / bEmA;
  77
  78     Er = bEmA;
  79     Ar = beta*E + A;
  80     Br = s2b * B;
  81     Cr = s2b * CibEmA * E;
  82     Dr = D + CibEmA * B;
  83
  84     ## Er = beta*E - A;
  85     ## Ar = beta*E + A;
  86     ## Br = sqrt (2*beta) * B;
  87     ## Cr = sqrt (2*beta) * C / (beta*E - A) * E;
  88     ## Dr = D + C / (beta*E - A) * B;
  89   endif
  90
  91 endfunction