summaryrefslogtreecommitdiffstats
path: root/scilab/modules/cacsd/macros/tf2ss.sci
blob: d33e554d0c3f713d2cf562bbee2d16e6b38a008e (plain)
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
// Scilab ( http://www.scilab.org/ ) - This file is part of Scilab
// Copyright (C) 2010 - INRIA - Serge Steer
//
// This file must be used under the terms of the CeCILL.
// This source file is licensed as described in the file COPYING, which
// you should have received as part of this distribution.  The terms
// are also available at
// http://www.cecill.info/licences/Licence_CeCILL_V2-en.txt

function [sl]=tf2ss(h,tol)
// Transfer function to state-space.
//Syntax : sl=tf2ss(h [,tol])
// h = transfer matrix
// sl = linear system in state-space representation (syslin list)
//!

  [lhs,rhs]=argn(0)

  if type(h)<=2 then
    sl=syslin([],[],[],[],h);
    return;
  end
  [num,den]=h(['num','den']);
  //
  [nd,md]=size(den)
  a=[];b=[];c=[];d=[];n1=0; // s=[]
  for k=1:md
    for l=1:nd
      [r,q]=pdiv(num(l,k),den(l,k));
      dk(l)=q;
      num(l,k)=r;
    end
    if nd<>1 then
      [nk,pp]=cmndred(num(:,k),den(:,k));
    else
      pp=den(k);nk=num(k);
    end;

    slk=cont_frm(nk,pp);
    [ak,bk,ck,dk1]=slk(2:5);
    // [s sk]
    [n2,m2]=size(bk);
    if n2<>0 then
      a(n1+1:n1+n2,n1+1:n1+n2)=ak;
      b(n1+1:n1+n2,k)=bk;
      c=[c ck];
      n1=n1+n2;
    else
      if n1<>0 then b(n1,k)=0;end
    end;
    d=[d dk];
  end;

  if degree(d)==0 then d=coeff(d),end

  if n1<>0 then
    nrmb=norm(b,1);
    nrmc=norm(c,1);
    fact=sqrt(nrmc*nrmb);
    b=b*fact/nrmb;
    c=c*fact/nrmc;
    atmp=a;
    [a,u]=balanc(a);
    //next lines commented out to fix bug 3796
    //   if rcond(u)< %eps*n1*n1*100
    //     nn=size(a,1);u=eye(nn,nn);a=atmp;
    //   end
    
    //apply transformation u without matrix inversion
    [k,l]=find(u<>0) //get the permutation
    u=u(k,l);c=c(:,k)*u; b=diag(1 ./diag(u))*b(k,:);

    if rhs<2 then
      [no,u]=contr(a',c',%eps);
    else
      [no,u]=contr(a',c',tol);
    end
    u=u(:,1:no);
    a=u'*a*u;b=u'*b;c=c*u;
    sl=syslin(h('dt'),a,b,c,d);
  else
    sl=syslin(h('dt'),[],[],[],d)
  end
endfunction