← Back to branch summary

~ubuntu-branches/ubuntu/karmic/scilab/karmic

« back to all changes in this revision

Viewing changes to macros/percent/%ar_p.sci

  • Committer: Bazaar Package Importer
  • Author(s): Torsten Werner
  • Date: 2002-03-21 16:57:43 UTC
  • Revision ID: james.westby@ubuntu.com-20020321165743-e9mv12c1tb1plztg
Tags: upstream-2.6
Import upstream version 2.6

Show diffs side-by-side

added added

removed removed

Lines of Context:
macros/percent/%ar_p.sci
 
1
function %ar_p(l)
 
2
// to provide automatic print of ar objects with armap
 
3
// Copyright ENPC 
 
4
// 
 
5
// We use armp_p as the default function for printing
 
6
armap_p(l)
 
7
 
 
8
function armap(ar,out)
 
9
//used for display of armax objects 
 
10
//!
 
11
[lhs,rhs]=argn(0)
 
12
if rhs==1 then out=%io(2),end
 
13
deff('[ta]=%cv(x)',['[m,n]=size(x);';
 
14
                    'frmt=format();frmt=10**frmt(2)/maxi([1,norm(x)]);';
 
15
                    'x=round(frmt*x)/frmt;';
 
16
                    't=[];for k=1:m,t=[t;''|''],end;';
 
17
                    'ta=t;for k=1:n,';
 
18
                    '        aa=string(x(:,k)),';
 
19
                    '        for l=1:m,';
 
20
                    '           if part(aa(l),1)<>''-'' then ';
 
21
                    '               aa(l)='' ''+aa(l),';
 
22
                    '           end,';
 
23
                    '        end,';
 
24
                    '        n=maxi(length(aa)),';
 
25
                    '        aa=part(aa+blank,1:n),';
 
26
                    '        ta=ta+aa+part(blank,1),';
 
27
                    'end;ta=ta+t;'])
 
28
 
 
29
// D(x)=Ax + Bu
 
30
//-------------
 
31
write(out,' ')
 
32
if ar(1)<>'ar' then write(%io(2)," This is not an ARMAX");return;end;
 
33
//
 
34
[ny,vid]=size(ar(2));
 
35
if ar(3) == [] then 
 
36
        write(out,"  A(z^-1)y= D(z^-1) e(t)");
 
37
else    
 
38
        write(out,"  A(z^-1)y=B(z^-1)u + D(z^-1) e(t)");
 
39
end 
 
40
 
 
41
write(out,"  ");
 
42
M=['A','B','D'];
 
43
lll=lines();
 
44
lll=lll(1);
 
45
if ar(3)==[] then I=[1,3];else I=1:3 ;end 
 
46
for i=I;
 
47
  [a]=ar(i+1)
 
48
  [na,lna]=size(a)
 
49
  nli=int((lna/na));
 
50
  if i==2,nli=int((lna/ar(6)));end;
 
51
  blank=[];for k=1:na,blank=[blank;'           '],end
 
52
  blank1=part(blank,1:7)
 
53
//  blank1([na/2,na/2+1])=['  '+M(i)+'(s)= ']
 
54
  blank1([int(na/2)+1])=['  '+M(i)+'(x)=']
 
55
  t=blank1;
 
56
  nna=na;
 
57
  blank1=part(blank,1:6)
 
58
  if i==2;nna=ar(6);end
 
59
  for j=1:nli;
 
60
    ta=a(:,1+(j-1)*nna:j*nna);
 
61
    if t==[],t=blank+%cv(ta);else t=t+%cv(ta);end;
 
62
    str='x^'+string(j-1)
 
63
    if length(str)==3,str=str+' ';end
 
64
    if j<>nli;str=str+'+ ';else str=str+'  ';end
 
65
//    blank1([na/2,na/2+1])=['        ';str];
 
66
    blank1([int(na/2)+1])=[str];
 
67
    t=t+blank1
 
68
    if length(t(1))>=lll,write(out,t),write(out," ");t=[];end
 
69
  end
 
70
write(out,t);
 
71
write(out," ");
 
72
end
 
73
write(%io(2),'  e(t)=Sig*w(t); w(t) '+string(ny)+'-dim white noise');
 
74
write(out," ");
 
75
  [a]=ar(7)
 
76
  [na,lna]=size(a)
 
77
  blank=[];for k=1:na,blank=[blank;'           '],end
 
78
  blank1=part(blank,1:8)
 
79
//  blank1([na/2,na/2+1])=['        ';'  Sig=  ']
 
80
  blank1([int(na/2)+1])=['  Sig=  ']
 
81
  t=blank1;
 
82
  ta=a;
 
83
  t=t+%cv(ta);
 
84
write(out,t);
 
85
write(out," ");
 
86
 
 
87
function armap_p(ar,out)
 
88
//used for display of armax objects 
 
89
//use Scilan Polynomial matrix display
 
90
//!
 
91
[lhs,rhs]=argn(0)
 
92
if rhs==1 then out=%io(2),end
 
93
deff('[ta]=%cv(x)',['[m,n]=size(x);';
 
94
                    'frmt=format();frmt=10**frmt(2)/maxi([1,norm(x)]);';
 
95
                    'x=round(frmt*x)/frmt;';
 
96
                    't=[];for k=1:m,t=[t;''|''],end;';
 
97
                    'ta=t;for k=1:n,';
 
98
                    '        aa=string(x(:,k)),';
 
99
                    '        for l=1:m,';
 
100
                    '           if part(aa(l),1)<>''-'' then ';
 
101
                    '               aa(l)='' ''+aa(l),';
 
102
                    '           end,';
 
103
                    '        end,';
 
104
                    '        n=maxi(length(aa)),';
 
105
                    '        aa=part(aa+blank,1:n),';
 
106
                    '        ta=ta+aa+part(blank,1),';
 
107
                    'end;ta=ta+t;'])
 
108
 
 
109
// D(x)=Ax + Bu
 
110
//-------------
 
111
write(out,' ')
 
112
if ar(1)<>'ar' then write(%io(2)," This is not an ARMAX");return;end;
 
113
//
 
114
[ny,vid]=size(ar(2));
 
115
if ar(3) == [] then 
 
116
        write(out,"  A(z^-1)y= D(z^-1) e(t)");
 
117
else    
 
118
        write(out,"  A(z^-1)y=B(z^-1)u + D(z^-1) e(t)");
 
119
end 
 
120
 
 
121
write(out,"  ");
 
122
M=['A','B','D'];
 
123
lll=lines();
 
124
lll=lll(1);
 
125
if ar(3)==[] then I=[1,3];else I=1:3 ;end 
 
126
A=inv_coeff(ar(2));
 
127
//print(out,A);
 
128
disp(A,'A(x) =');
 
129
write(out,"  ");
 
130
if ar(3)<>[] then 
 
131
        [mb,nb]=size(ar(3));
 
132
        B=inv_coeff(ar(3),(nb/ar('nu'))-1);
 
133
        //print(out,B);
 
134
        disp(B, 'B(x) =');
 
135
        write(out,"  ");
 
136
end
 
137
D=inv_coeff(ar(4));
 
138
//print(out,D);
 
139
disp(D, 'D(x)');
 
140
write(out," ");
 
141
write(%io(2),'  e(t)=Sig*w(t); w(t) '+string(ny)+'-dim white noise');
 
142
write(out," ");
 
143
  [a]=ar(7)
 
144
  [na,lna]=size(a)
 
145
  blank=[];for k=1:na,blank=[blank;'           '],end
 
146
  blank1=part(blank,1:8)
 
147
  //  blank1([na/2,na/2+1])=['        ';'  Sig=  ']
 
148
  blank1([int(na/2)+1])=['  Sig=  ']
 
149
  t=blank1;
 
150
  ta=a;
 
151
  t=t+%cv(ta);
 
152
write(out,t);
 
153
write(out," ");
 
154
 
 
155
 
 
156
 
 
157
 
 
158
 
 
159