diff options
| author | priyanka | 2015-06-24 15:03:17 +0530 |
|---|---|---|
| committer | priyanka | 2015-06-24 15:03:17 +0530 |
| commit | b1f5c3f8d6671b4331cef1dcebdf63b7a43a3a2b (patch) | |
| tree | ab291cffc65280e58ac82470ba63fbcca7805165 /2048/DEPENDENCIES | |
| download | Scilab-TBC-Uploads-b1f5c3f8d6671b4331cef1dcebdf63b7a43a3a2b.tar.gz Scilab-TBC-Uploads-b1f5c3f8d6671b4331cef1dcebdf63b7a43a3a2b.tar.bz2 Scilab-TBC-Uploads-b1f5c3f8d6671b4331cef1dcebdf63b7a43a3a2b.zip | |
initial commit / add all books
Diffstat (limited to '2048/DEPENDENCIES')
37 files changed, 1655 insertions, 0 deletions
diff --git a/2048/DEPENDENCIES/armac1.sci b/2048/DEPENDENCIES/armac1.sci new file mode 100755 index 000000000..75b8ee7c3 --- /dev/null +++ b/2048/DEPENDENCIES/armac1.sci @@ -0,0 +1,17 @@ +// Scilab description of an ARMAX process
+// Form:
+// A(q) y(t) = [B(q)/F(q)] u(t-nk) + [C(q)/D(q)] e(t)
+
+// [A(q)*F(q)*D(q)] y(t) = [B(q)*D(q)] u(t-nk) + [C(q)*F(q)]e(t)
+// A1(q) = [A(q)*F(q)*D(q)]
+// B1(q) = [B(q)*D(q)]
+// D1(q) = [C(q)*F(q)]
+
+function process_ar = armac1(a,b,c,d,f,sig)
+ny = 1; nu =1;
+a1 = convol(convol(a,f),d);
+b1 = convol(b,d);
+d1 = convol(c,f);
+process_ar = armac(a1,b1,d1,ny,nu,sig);
+
+endfunction;
diff --git a/2048/DEPENDENCIES/arx.sci b/2048/DEPENDENCIES/arx.sci new file mode 100755 index 000000000..038217da2 --- /dev/null +++ b/2048/DEPENDENCIES/arx.sci @@ -0,0 +1,110 @@ +// ARX model parameter estimation
+// Computes Covariance matrix
+// Computes Noisevariance of a process
+
+function [thetaN,covt,nvar,res] = arx(data,na,nb,nk)
+az = max(na,nb+nk-1);
+zer = zeros(az,1);
+zd = data;
+// Zeros appended
+zd1(:,1) = [zer; zd(:,1)];
+zd1(:,2) = [zer; zd(:,2)];
+[r,c] = size(zd1);
+t = az+1:r;
+yt = zd1(:,1); ut = zd1(:,2);
+yt1 = yt'; ut1 = ut'; // row vector
+len1 = length(yt1);
+yt2 = zeros(1,len1-az); ut2 = zeros(1,len1-az);
+
+// arx(Data,[na nb nk])
+ for i=1:na
+ yt2 = [yt2; -yt1(t-i)];
+ end;
+ for i=nk:nb+nk-1
+ ut2 = [ut2; ut1(t-i)];
+ end;
+[r1_a,c1_a] = size(yt2); [r2_a,c2_a] = size(ut2);
+phit = [yt2(2:r1_a,:); ut2(2:r2_a,:)];
+m1 = phit*phit';
+[qm,rm] = qr(m1);
+m2 = phit*zd(:,1);
+thetaN = inv(rm)*qm'*m2;
+// thetaN = inv(m1)*m2;
+// thetaN = m1\m2;
+
+[r1,c1] = size(thetaN);
+a = thetaN(1:na); b = thetaN(na+1:r1);
+
+// Sum of squared residuals
+
+yhat = phit'*thetaN;
+res = zd(:,1) - yhat;
+N = length(res);
+q = rank(phit);
+ssr = res'*res;
+sig2 = ssr/(N-q);
+nvar = sqrt(sig2);
+cov = inv(m1);
+covt = diag(cov);
+
+a = thetaN(1:na); b = thetaN(na+1:r1);
+cova = covt(1:na); covb = covt(na+1:r1);
+x = poly(0,'x');
+disp('Discrete time model: A(x)y(t) = B(x)u(t) + e(t)');
+A = poly([1 a'],'x','coeff');
+cov_a = [0 cova'];
+b1 = zeros(1,nk);
+B = poly([b1 b'],'x','coeff');
+cov_b = covb';
+
+
+
+//////////////////////////////////////
+//////////////////////////////////////
+
+////////// Model Display /////////////
+
+function disp_mod(N,covN)
+len = length(covN);
+B1 = pol2str(N);
+ind = strindex(B1,['+','-']);
+ind = ind - 1;
+B2 = strsplit(B1,ind);
+covB = string(covN);
+
+ if ascii(B2(1)) == 32
+ B2 = B2(2:len+1);
+ end;
+
+ B3(1) = ' ';
+ for i=1:len
+ B3(i) = strsubst(B2(i),'*x','(+-' + covB(i) + ')*x');
+ end;
+
+ B4 = B3(1);
+
+ for i=2:len
+ B4 = B4 + ' ' + B3(i);
+ end;
+
+disp(B4);
+endfunction;
+//////////////////////////////////////
+//////////////////////////////////////
+
+ if na==0
+ disp('A(x) = 1');
+ else
+ disp('A(x) = ');
+ disp_mod(A,cov_a);
+ end;
+
+ if nb~=0
+ disp('B(x) = ');
+ disp_mod(B,cov_b);
+ end;
+
+
+
+endfunction;
+
diff --git a/2048/DEPENDENCIES/ch_pol.sci b/2048/DEPENDENCIES/ch_pol.sci new file mode 100755 index 000000000..fa5fa8af2 --- /dev/null +++ b/2048/DEPENDENCIES/ch_pol.sci @@ -0,0 +1,10 @@ +// function [phi,psi] = ch_pol(N,epsilon)
+// Returns desired characteristic polynomial and numerator
+// N = rise time in number of sample times
+// epsilon = overshoot as a fraction of ss gain
+
+function [phi,psi] = ch_pol(N,epsilon)
+omega = %pi/2/N; r = epsilon^(omega/%pi);
+phi = [1 -2*r*cos(omega) r^2];
+psi = [1-r*cos(omega) (r^2-r*cos(omega))];
+endfunction;
diff --git a/2048/DEPENDENCIES/cindep.sci b/2048/DEPENDENCIES/cindep.sci new file mode 100755 index 000000000..3ffd243d4 --- /dev/null +++ b/2048/DEPENDENCIES/cindep.sci @@ -0,0 +1,35 @@ +// function b = cindep( S,gap)
+// Used in XD + YN = C. All rows except the last of are assumed to
+// be independent. The aim is to check if the last row is dependent on the
+// rest and if so how. The coefficients of dependence are sent in b.
+function b = cindep( S,gap)
+
+if argn(2) == 1
+ gap = 1.0e8;
+end
+eps = 2.2204e-016;
+[rows,cols] = size(S);
+if rows > cols
+ ind = 0;
+else
+ sigma = svd(S);
+ len = length(sigma);
+ if (sigma(len)/sigma(1) <= (eps*max(i,cols)))
+ ind = 0; //not independent
+ else
+ if or(sigma(1:len-1) ./sigma(2:len)>=gap)
+ ind = 0; // not dependent
+ else
+ ind = 1; //independent
+ end
+ end
+end
+if ind
+ b = [];
+else
+ b = S(rows,:)/S(1:rows-1,:);
+ b = makezero(b,gap);
+end
+endfunction
+
+
diff --git a/2048/DEPENDENCIES/clcoef.sci b/2048/DEPENDENCIES/clcoef.sci new file mode 100755 index 000000000..249a5370e --- /dev/null +++ b/2048/DEPENDENCIES/clcoef.sci @@ -0,0 +1,25 @@ +// H. Kwakernaak, July, 1990
+// Modified by Kannan Moudgalya in Nov. 1992
+
+function [P,degP] = clcoef(Q,degQ)
+
+[rQ,cQ] = polsize(Q,degQ);
+
+if and(and(Q==0))
+ P = zeros(rQ,cQ);
+ degP = 0;
+else
+ P = Q; degP = degQ; rP = rQ; cP = cQ;
+ j = degP+1;
+ while j >= 0
+ X = P(:,(j-1)*cP+1:j*cP)
+ if max(sum(abs(X'))) < (1e-8)*max(sum(abs(P)))
+ P = P(:,1:(j-1)*cP);
+ degP = degP-1;
+ else
+ j = 0;
+ end
+ j = j-1;
+ end
+end
+endfunction
diff --git a/2048/DEPENDENCIES/colsplit.sci b/2048/DEPENDENCIES/colsplit.sci new file mode 100755 index 000000000..ec4cafecc --- /dev/null +++ b/2048/DEPENDENCIES/colsplit.sci @@ -0,0 +1,36 @@ +// colsplit
+// The command
+// [P1,degP1,P2,degP2] = colsplit(P,degP,p1,p2)
+// produces two polynomial matrix P1 and P2. P1 consists of the first
+// p1 columns of P and P2 consists of the remaining p2 columns of P.
+
+// H. Kwakernaak, July, 1990
+
+
+function [P1,degP1,P2,degP2] = colsplit(P,degP,p1,p2)
+
+if isempty(P)
+ P1 = []; P2 = [];
+ degP1 = 0; degP2 = 0;
+ return;
+end
+
+[rP,cP] = polsize(P,degP);
+if p1 < 0 | p1 > cP | p2 < 0 | p2 > cP | p1+p2 ~= cP
+ error('colsplit: Inconsistent numbers of columns');
+end
+rP1 = rP; rP2 = rP; cP1 = p1; cP2 = p2;
+degP1= degP; degP2 = degP;
+
+if p1 == 0
+ P1 == []; P2 = P;
+elseif p2 == 0
+ P1 = P; P2 = [];
+else
+ P1 = zeros(rP1,(degP1+1)*cP1); P2 = zeros(rP2,(degP2+1)*cP2);
+ for i = 1:degP+1
+ P1(:,(i-1)*cP1+1:i*cP1) = P(:,(i-1)*cP+1:(i-1)*cP+cP1);
+ P2(:,(i-1)*cP2+1:i*cP2) = P(:,(i-1)*cP+cP1+1:i*cP);
+ end
+end
+endfunction;
diff --git a/2048/DEPENDENCIES/cosfil_ip.sci b/2048/DEPENDENCIES/cosfil_ip.sci new file mode 100755 index 000000000..b16ba4b27 --- /dev/null +++ b/2048/DEPENDENCIES/cosfil_ip.sci @@ -0,0 +1,15 @@ +// Input arguments are co efficients of numerator and denominator
+// polynomials in ascending powers of z^-1
+
+// Scicos/Xcos blocks need input polynomials
+// with positive powers of z
+
+function [nume,deno] = cosfil_ip(num,den)
+exec('polyno.sci',-1);
+[Nn,Nd] = polyno(num,'z');
+[Dn,Dd] = polyno(den,'z');
+nume = Nn*Dd;
+deno = Nd*Dn;
+
+endfunction;
+
diff --git a/2048/DEPENDENCIES/covar_m.sci b/2048/DEPENDENCIES/covar_m.sci new file mode 100755 index 000000000..bd84831b2 --- /dev/null +++ b/2048/DEPENDENCIES/covar_m.sci @@ -0,0 +1,46 @@ +//User defined equivalent function to Matlab covar function
+//For discrete time domain only
+//Uses Lyapunov's equation for computation
+//W: noise intensity (scalar)
+//Matlab result for unstable systems:
+//Warning: Unstable models have infinite covariance
+
+function P = covar_m(sys,W)
+a = roots(sys('den'));
+b = length(a);
+c = abs(a) > 1;
+if b ~= 0 then
+ for i = 1:b
+ if c(i) == %t then
+ disp('Warning: System being unstable has infinite covariance');
+ P = %inf;
+ else
+ s = tf2ss(sys);
+ [A,B,C,D] = s(2:5);
+ //Sylvester and Lyapunov solver
+ task = 2; flag = [1 0]; tran = 1;
+ Q1 = -B*W*B';
+ Q = linmeq(task,A,Q1,flag,tran)
+ P = C*Q*C' + D*W*D';
+ end;
+ end;
+else
+ s = tf2ss(sys);
+ [A,B,C,D] = s(2:5);
+ //Sylvester and Lyapunov solver
+ task = 2; flag = [1 0]; tran = 1;
+ Q1 = -B*W*B';
+ Q = linmeq(task,A,Q1,flag,tran)
+ P = C*Q*C' + D*W*D';
+end;
+endfunction;
+
+//if d==0 | c==%f
+// disp('Calc');
+//else
+// disp('Unstable');
+//end;
+// Above logic can also solve our purpose
+// But it gives incorrect answer if roots are [1 -1 1] or
+// [1 -1 -1] ....
+
diff --git a/2048/DEPENDENCIES/covf.sci b/2048/DEPENDENCIES/covf.sci new file mode 100755 index 000000000..2cf7e74d6 --- /dev/null +++ b/2048/DEPENDENCIES/covf.sci @@ -0,0 +1,30 @@ +// z: [y u] (Column matrix)
+// ---- ----
+// | |
+// | ryy(0) ryy(1) ... ryy(len) |
+// | |
+// v = | ruy(0) ruy(1) ... ruy(len) |
+// | |
+// | ryu(0) ryu(1) ... ryu(len) |
+// | |
+// | ruu(0) ruu(1) ... ruu(len) |
+// | |
+// ---- ----
+// No. of lags
+
+function v = covf(z,n)
+[r c] = size(z);
+y1 = z(:,1);
+u1 = z(:,2);
+ if n>r
+ // disp('Incorrect number of lags');
+ n = r;
+ else
+ r = n;
+ end;
+ryu = corr(y1,u1,r);
+ruy = corr(u1,y1,r);;
+ryy = corr(y1,r);
+ruu = corr(u1,r);
+v = [ryy; ruy; ryu; ruu]
+endfunction;
diff --git a/2048/DEPENDENCIES/cra.sci b/2048/DEPENDENCIES/cra.sci new file mode 100755 index 000000000..be0d6e2f2 --- /dev/null +++ b/2048/DEPENDENCIES/cra.sci @@ -0,0 +1,50 @@ +
+// [ir,r,cl] = cra(z,M,n)
+// M: No. of lags
+// n: Order of pre whitening filter
+// z is of the form:
+// z = [y u]
+// y and u are column matrices
+
+function [ir,r,cl] = cra(z,varargin)
+len = length(varargin);
+
+ if len==0,M = 20;n = 10;
+ elseif len==1,M = varargin(1);n=10;
+ else M = varargin(1);n = varargin(2);
+ end;
+
+[ro,co] = size(z);
+
+a1 = armax1(n,-1,0,z(:,2)',zeros(1,ro));
+[A,B,D] = arma2p(a1);
+a11 = coeff(A);
+z22(1,:) = filt(a11,1,z(:,1)');
+z22(2,:) = filt(a11,1,z(:,2)');
+Covr = covf(z22',M+1);
+scir=Covr(4,1);
+sccf=sqrt(Covr(1,1)*Covr(4,1));
+
+ir = Covr(2,:)'/scir;
+
+r(:,1) = -M:M;
+r(:,2) = [Covr(1,M+1:-1:1) Covr(1,2:M+1)]'
+r(:,3) = [Covr(4,M+1:-1:1) Covr(4,2:M+1)]'
+rhoyu = Covr(3,:)'/sccf;
+rhouy = Covr(2,:)'/sccf;
+r(:,4) = [rhoyu(M+1:-1:1); rhouy(2:M+1)];
+
+sdreu=2.58*sqrt(r(:,3)'*r(:,2))/scir/sqrt(ro)*ones(2*M+1,1);
+cl=sdreu(1);
+
+timeax=[0:length(ir)-1];
+plot(timeax,ir,'bo');
+plot2d3(timeax,ir,style = 2);
+plot(timeax,cl*ones(1,length(ir)),'b-.');
+plot(timeax,-cl*ones(1,length(ir)),'b-.');
+xtitle('Impulse response estimate','Lags');
+xgrid();
+
+endfunction;
+
+
diff --git a/2048/DEPENDENCIES/deconvol.sci b/2048/DEPENDENCIES/deconvol.sci new file mode 100755 index 000000000..b3bc8b632 --- /dev/null +++ b/2048/DEPENDENCIES/deconvol.sci @@ -0,0 +1,21 @@ +function q = deconvol(b,a)
+
+ if a(1) == 0
+ error('1st coefficient of A must be nonzero');
+ end;
+
+exec('filt.sci',-1);
+
+nb = length(b); na = length(a);
+
+ if na>nb
+ q = 0;
+ else
+ q = filt(b,a,[1 zeros(1,nb-na)]);
+ end;
+q = clean(q);
+endfunction;
+
+
+
+
diff --git a/2048/DEPENDENCIES/delc2d.sci b/2048/DEPENDENCIES/delc2d.sci new file mode 100755 index 000000000..c038b64b2 --- /dev/null +++ b/2048/DEPENDENCIES/delc2d.sci @@ -0,0 +1,48 @@ +// Discretizing a tf with delay
+// Exact solution
+// Applicable for a 1st order system
+
+// Ref.: pg.287,Digital Control,Prof.Kannan Moudgalya
+
+// D: Delay
+// TF: e^(-Ds) OR e^(-Ds)
+// ------------ ------------ (gen.)
+// tau*s + 1 tau*s + a
+
+//D = kTs + D' (gen.)
+// G: TF with delay component
+// G1: TF with zero delay
+// Required because G cannot be directly used in Scilab
+// Coefficients are returned for ascending powers of z^-1
+
+function [B,A,k1] = delc2d(G,G1,Ts)
+D = G.iodelay;
+d = coeff(G1('den'));
+ if d(1) == 1
+ tau = d(2);
+ mu = 1;
+ else
+ tau = d(2)/d(1);
+ mu = 1/d(1);
+ end;
+
+k = floor(D/Ts);
+Dpri = D - k*Ts;
+Dis = ((%z*(1 - (exp(-(Ts - Dpri)/tau)) ) )+ (exp(-(Ts - Dpri)/tau) - exp(-Ts/tau) ))/ ((%z^(k+1))*(%z - exp(-Ts/tau)))
+Dis1 = Dis*mu;
+disp('Warning: Exact discretization of first order system only');
+k1 = degree(Dis1('den')) - degree(Dis1('num'));
+B = coeff(Dis1('num'));
+A = coeff(Dis1('den'));
+B = flip(B);
+A = flip(A);
+endfunction;
+
+
+
+
+
+
+
+
+
diff --git a/2048/DEPENDENCIES/ext.sci b/2048/DEPENDENCIES/ext.sci new file mode 100755 index 000000000..469e68f2a --- /dev/null +++ b/2048/DEPENDENCIES/ext.sci @@ -0,0 +1,14 @@ +// function [B,degB] = ext(A,degA,k,l)
+// EXTRACTS THE (k,l) ELEMENT OF A polynomial matrix A into B
+
+function [B,degB] = ext(A,degA,k,l)
+
+[rA,cA] = polsize(A,degA);
+degB = degA;
+B = zeros(1,degB+1);
+for m = 0:degB
+ B(1,m+1) = A(k,(m*cA)+l);
+end
+
+[B,degB] = clcoef(B,degB);
+endfunction;
diff --git a/2048/DEPENDENCIES/filt.sci b/2048/DEPENDENCIES/filt.sci new file mode 100755 index 000000000..335a18626 --- /dev/null +++ b/2048/DEPENDENCIES/filt.sci @@ -0,0 +1,28 @@ +// Filters a data sequence using a digital filter
+// a(1)*y(n) = b(1)*x(n) + b(2)*x(n-1) + ... + b(nb+1)*x(n-nb)
+// - a(2)*y(n-1) - ... - a(na+1)*y(n-na)
+
+// x: input signal
+// x has maximum length
+function y = filt(b,a,x)
+l = length(a);
+m = length(b);
+n = length(x);
+b = [b zeros(1,n-m)];
+a = [a zeros(1,n-l)];
+
+y1(1) = b(1)*x(1);
+y2(1) = 0;
+y = y1 + y2;
+y = y/a(1);
+ for i=2:n
+ y1 = b(1:i) .* x(i:-1:1);
+ z1 = sum(y1);
+ y2 = a(2:i).*y(i-1:-1:1);
+ z2 = sum(y2);
+ y3 = z1 - z2;
+ y3 = y3/a(1);
+ y = [y y3];
+ end;
+
+endfunction;
diff --git a/2048/DEPENDENCIES/indep.sci b/2048/DEPENDENCIES/indep.sci new file mode 100755 index 000000000..860abae31 --- /dev/null +++ b/2048/DEPENDENCIES/indep.sci @@ -0,0 +1,38 @@ +// function b = indep(S,gap)
+// determines the first row that is dependent on the previous rows of S.
+// The coefficients of dependence is returned in b
+function b = indep( S,gap)
+
+if argn(2) == 1
+ gap = 1.0e8;
+ end
+[rows,cols] = size(S);
+ind = 1;
+i = 2;
+eps = 2.2204e-016;
+while ind & i <= rows
+ sigma = svd(S(1:i,:));
+ len = length(sigma);
+ if(sigma(len)/sigma(1) < (eps*max(i,cols)))
+ ind =0;
+ else
+ shsig = [sigma(2:len);sigma(len)];
+ if or( (sigma ./shsig) > gap)
+ ind = 0;
+ else
+ ind = 1;
+ i = i+1;
+ end
+ end
+
+end
+if ind
+ b =[];
+
+else
+ c = S(i,:)/S(1:i-1,:);
+ c = makezero(c,gap);
+ b = [-c 1];
+end
+endfunction
+
diff --git a/2048/DEPENDENCIES/iodelay.sci b/2048/DEPENDENCIES/iodelay.sci new file mode 100755 index 000000000..1b70d6265 --- /dev/null +++ b/2048/DEPENDENCIES/iodelay.sci @@ -0,0 +1,360 @@ +function Hd=iodelay(H,d)
+// Author: Serge Steer, Copyright INRIA
+ if or(size(H.num)<>size(d)) then error('Dimension mismatch'),end
+ Hd=mlist(['rd','H','iodelay'],H,d)
+endfunction
+
+//element wise product overloading functions
+function H=%rd_x_rd(H1,H2)
+// Author: Serge Steer, Copyright INRIA
+ H=mlist(['rd','H','iodelay'],H1.H.*H2.H,H1.iodelay+ H2.iodelay)
+endfunction
+function H=%rd_x_r(H1,H2)
+// Author: Serge Steer, Copyright INRIA
+ H=mlist(['rd','H','iodelay'],H1.H.*H2,H1.iodelay)
+endfunction
+function H=%r_x_rd(H1,H2)
+// Author: Serge Steer, Copyright INRIA
+ H=mlist(['rd','H','iodelay'],H1.*H2.H,H2.iodelay)
+endfunction
+function H=%rd_x_s(H1,H2)
+// Author: Serge Steer, Copyright INRIA
+ H=mlist(['rd','H','iodelay'],H1.H.*H2,H1.iodelay)
+endfunction
+function H=%s_x_rd(H1,H2)
+// Author: Serge Steer, Copyright INRIA
+ H=mlist(['rd','H','iodelay'],H1.*H2.H,H2.iodelay)
+endfunction
+
+//column concatenation overloading functions
+function H=%rd_c_rd(H1,H2)
+// Author: Serge Steer, Copyright INRIA
+ H=mlist(['rd','H','iodelay'],[H1.H H2.H],[H1.iodelay,H2.iodelay])
+endfunction
+function H=%rd_c_r(H1,H2)
+// Author: Serge Steer, Copyright INRIA
+ H=mlist(['rd','H','iodelay'],[H1.H H2],[H1.iodelay,zeros(H2.num)])
+endfunction
+function H=%r_c_rd(H1,H2)
+// Author: Serge Steer, Copyright INRIA
+ H=mlist(['rd','H','iodelay'],[H1 H2.H],[zeros(H1.num),H2.iodelay])
+endfunction
+function H=%rd_c_s(H1,H2)
+// Author: Serge Steer, Copyright INRIA
+ H=mlist(['rd','H','iodelay'],[H1.H H2],[H1.iodelay,zeros(H2.num)])
+endfunction
+function H=%s_c_rd(H1,H2)
+// Author: Serge Steer, Copyright INRIA
+ H=mlist(['rd','H','iodelay'],[H1 H2.H],[zeros(H1.num),H2.iodelay])
+endfunction
+
+//row concatenation overloading functions
+function H=%rd_f_rd(H1,H2)
+// Author: Serge Steer, Copyright INRIA
+ H=mlist(['rd','H','iodelay'],[H1.H;H2.H],[H1.iodelay;H2.iodelay])
+endfunction
+function H=%rd_f_r(H1,H2)
+// Author: Serge Steer, Copyright INRIA
+ H=mlist(['rd','H','iodelay'],[H1.H;H2],[H1.iodelay;zeros(H2.num)])
+endfunction
+function H=%r_f_rd(H1,H2)
+// Author: Serge Steer, Copyright INRIA
+ H=mlist(['rd','H','iodelay'],[H1;H2.H],[zeros(H1.num);H2.iodelay])
+endfunction
+function H=%rd_f_s(H1,H2)
+// Author: Serge Steer, Copyright INRIA
+ H=mlist(['rd','H','iodelay'],[H1.H;H2],[H1.iodelay;zeros(H2.num)])
+endfunction
+function H=%s_f_rd(H1,H2)
+// Author: Serge Steer, Copyright INRIA
+ H=mlist(['rd','H','iodelay'],[H1;H2.H],[zeros(H1.num);H2.iodelay])
+endfunction
+
+//matrix product overloading functions
+function H=%rd_m_rd(H1,H2)
+// Author: Serge Steer, Copyright INRIA
+ D1=H1.iodelay
+ D2=H2.iodelay
+ D=D1*D2;
+ for k=1;size(D2,2)
+ for l=1:size(D1,1)
+ d=D1(l,:)+(D2(:,k)')
+ if or(d(1)<>d) then error('Delays mismatched'),end
+ D(l,k)=d
+ end
+ end
+ H=mlist(['rd','H','iodelay'],H1.H*H2.H,D)
+endfunction
+function H=%rd_m_r(H1,H2)
+// Author: Serge Steer, Copyright INRIA
+ D1=H1.iodelay
+ for l=1:size(D1,1)
+ if or(D1(l,1)<>D1(l,:)) then error('Delays mismatched'),end
+ end
+ H=mlist(['rd','','iodelay'],H1.H*H2,D1)
+endfunction
+function H=%rd_m_s(H1,H2)
+// Author: Serge Steer, Copyright INRIA
+ D1=H1.iodelay
+ for l=1:size(D1,1)
+ if or(D1(l,1)<>D1(l,:)) then error('Delays mismatched'),end
+ end
+ H=mlist(['rd','H','iodelay'],H1.H*H2,D1)
+endfunction
+function H=%r_m_rd(H1,H2)
+// Author: Serge Steer, Copyright INRIA
+ D2=H2.iodelay
+ for l=1:size(D2,2)
+ if or(D2(1,l)<>D2(:,l)) then error('Delays mismatched'),end
+ end
+ H=mlist(['rd','H','iodelay'],H1*H2.H,D2)
+endfunction
+function H=%s_m_rd(H1,H2)
+// Author: Serge Steer, Copyright INRIA
+ D2=H2.iodelay
+ for l=1:size(D2,2)
+ if or(D2(1,l)<>D2(:,l)) then error('Delays mismatched'),end
+ end
+ H=mlist(['rd','H','iodelay'],H1*H2.H,D2)
+endfunction
+
+//Addition overloading functions
+function H=%rd_a_rd(H1,H2)
+// Author: Serge Steer, Copyright INRIA
+ D1=H1.iodelay
+ D2=H2.iodelay
+ if or(D1<>D2) then error('Delay mismatch'),end
+ H=mlist(['rd','H','iodelay'],H1.H+H2.H,D1)
+endfunction
+function H=%rd_a_r(H1,H2)
+// Author: Serge Steer, Copyright INRIA
+ D1=H1.iodelay
+ if or(D1<>0) then error('Delay mismatch'),end
+ H=mlist(['rd','H','iodelay'],H1.H+H2,D1)
+endfunction
+function H=%rd_a_s(H1,H2)
+// Author: Serge Steer, Copyright INRIA
+ D1=H1.iodelay
+ if or(D1<>0) then error('Delay mismatch'),end
+ H=mlist(['rd','H','iodelay'],H1.H+H2,D1)
+endfunction
+function H=%r_a_rd(H1,H2)
+// Author: Serge Steer, Copyright INRIA
+ D2=H2.iodelay
+ if or(D2<>0) then error('Delay mismatch'),end
+ H=mlist(['rd','H','iodelay'],H1+H2.H,D2)
+endfunction
+
+//Substraction overloading functions
+function H=%rd_a_rd(H1,H2)
+// Author: Serge Steer, Copyright INRIA
+ D1=H1.iodelay
+ D2=H2.iodelay
+ if or(D1<>D2) then error('Delay mismatch'),end
+ H=mlist(['rd','H','iodelay'],H1.H-H2.H,D1)
+endfunction
+function H=%rd_a_r(H1,H2)
+// Author: Serge Steer, Copyright INRIA
+ D1=H1.iodelay
+ if or(D1<>0) then error('Delay mismatch'),end
+ H=mlist(['rd','H','iodelay'],H1.H-H2,D1)
+endfunction
+function H=%rd_a_s(H1,H2)
+// Author: Serge Steer, Copyright INRIA
+ D1=H1.iodelay
+ if or(D1<>0) then error('Delay mismatch'),end
+ H=mlist(['rd','H','iodelay'],H1.H-H2,D1)
+endfunction
+function H=%r_a_rd(H1,H2)
+// Author: Serge Steer, Copyright INRIA
+ D2=H2.iodelay
+ if or(D2<>0) then error('Delay mismatch'),end
+ H=mlist(['rd','H','iodelay'],H1-H2.H,D2)
+endfunction
+
+//extraction overloading functions
+function H=%rd_e(varargin)
+// Author: Serge Steer, Copyright INRIA
+ H=varargin($)
+ Hh= H.H;Hh.num=Hh.num(varargin(1:$-1));Hh.den=Hh.den(varargin(1:$-1));
+ H.H=Hh
+ H.iodelay=H.iodelay(varargin(1:$-1))
+endfunction
+
+//insertion overloading functions
+function H=%rd_i_rd(varargin)
+// Author: Serge Steer, Copyright INRIA
+ H=varargin($)
+ Hfrom=varargin($-1)
+
+ Hh= H.H
+ Hh.num(varargin(1:$-2))=Hfrom.H.num
+ Hh.den(varargin(1:$-2))=Hfrom.H.den
+ H.H=Hh
+ H.iodelay(varargin(1:$-2))=Hfrom.iodelay
+endfunction
+function H=%r_i_rd(varargin)
+// Author: Serge Steer, Copyright INRIA
+ H=varargin($)
+ Hfrom=varargin($-1)
+
+ Hh= H.H
+ Hh.num(varargin(1:$-2))=Hfrom.num
+ Hh.den(varargin(1:$-2))=Hfrom.den
+ H.H=Hh
+ H.iodelay(varargin(1:$-2))=zeros(Hfrom.num)
+endfunction
+function H=%s_i_rd(varargin)
+// Author: Serge Steer, Copyright INRIA
+ H=varargin($)
+ Hfrom=varargin($-1)
+
+ Hh= H.H
+ Hh.num(varargin(1:$-2))=Hfrom
+ Hh.den(varargin(1:$-2))=ones(Hfrom)
+ H.H=Hh
+ H.iodelay(varargin(1:$-2))=zeros(Hfrom)
+endfunction
+function H=%rd_i_r(varargin)
+// Author: Serge Steer, Copyright INRIA
+ Hh=varargin($)
+ Hfrom=varargin($-1)
+ Hh.num(varargin(1:$-2))=Hfrom.H.num
+ Hh.den(varargin(1:$-2))=Hfrom.H.den
+ D=zeros(Hh.num)
+ D(varargin(1:$-2))=Hfrom.iodelay
+ H=mlist(['rd','H','iodelay'],Hh,D)
+endfunction
+function H=%rd_i_s(varargin)
+// Author: Serge Steer, Copyright INRIA
+ num=varargin($);den=ones(num)
+ Hfrom=varargin($-1)
+
+ num(varargin(1:$-2))=Hfrom.H.num
+ den(varargin(1:$-2))=Hfrom.H.den
+ D=zeros(num)
+ D(varargin(1:$-2))=Hfrom.iodelay
+ H=mlist(['rd','H','iodelay'],rlist(num,den,'c'),D)
+endfunction
+
+//transpostion overloading function
+function H=%rd_t(H)
+// Author: Serge Steer, Copyright INRIA
+ H.H=H.H'
+ H.iodelay=H.iodelay'
+endfunction
+
+//string overloading function
+function txt=%rd_string(H)
+// Author: Serge Steer, Copyright INRIA
+ r=H.H
+ N=string(r.num)
+ D=string(r.den)
+ ln=max(matrix(length(N),2,-1),'r')
+ ld=max(matrix(length(D),2,-1),'r')
+ l=max(ln,ld)
+ [m,n]=size(r.num);
+ S=emptystr(m,n)
+
+ kz=find(H.iodelay==0);//zero delay entries
+ w='exp('+string(-H.iodelay)+'*s)*';w(kz)='';
+
+ for i=1:m*n
+ s=2*i-1:2*i
+ pw=part(' ',1:length(w(i)))
+ N(s)=pw+part(' ',1:(l(i)-ln(i))/2)+N(s)
+ D(s)=pw+part(' ',1:(l(i)-ld(i))/2)+D(s)
+ S(i) =w(i)+part('-',ones(1,l(i)))
+ end
+ txt=emptystr(5*m,n);
+ txt(1:5:$,:)=N(1:2:$,:)
+ txt(2:5:$,:)=N(2:2:$,:)
+ txt(3:5:$,:)=S(1:$,:)
+ txt(4:5:$,:)=D(1:2:$,:)
+ txt(5:5:$,:)=D(2:2:$,:)
+endfunction
+
+
+//Display overloading function
+function %rd_p(H)
+// Author: Serge Steer, Copyright INRIA
+//used to display rational fraction with complex coefficients
+//The real case is hard coded
+// Copyright INRIA
+ T=string(H)
+ l=max(length(T),'r')
+ Te=emptystr(size(T,1),1)
+ for k=1:size(T,2)
+ Te=Te+part(T(:,k),1:l(k)+1)+' '
+ end
+ mprintf("%s\n",Te)
+endfunction
+
+
+//frequency response computation overloading
+function [frq,repf,splitf]=repfreq(varargin)
+// Author: Serge Steer, Copyright INRIA
+ H=varargin(1)
+ autolib.repfreq
+ if typeof(H)=='rd' then
+ D=H.iodelay
+ varargin(1)=H.H
+ [frq,repf,splitf]=repfreq(varargin(:))
+
+ w=-2*%i*%pi*frq;
+ for k=1:size(D,'*')
+ repf(k,:)=repf(k,:).*exp(D(k)*w)
+ end
+ else
+ [frq,repf,splitf]=repfreq(varargin(:))
+ end
+endfunction
+
+//bode plot computation overloading
+function bode(varargin)
+// Author: Serge Steer, Copyright INRIA
+ H=varargin(1)
+ xdesslib.bode
+ if typeof(H)=='rd' then
+ if type(varargin($))==10 then //a comment
+ cm=varargin($);varagin($)=null()
+ else
+ cm=''
+ end
+ D=H.iodelay
+ varargin(1)=H.H
+ [frq,repf,splitf]=repfreq(varargin(1:$))
+ w=-2*%i*%pi*frq;
+ for k=1:size(D,'*')
+ repf(k,:)=repf(k,:).*exp(D(k)*w)
+ end
+ bode(frq,repf)
+ else
+ bode(varargin(:))
+ end
+endfunction
+
+//nyquist plot computation overloading
+function nyquist(varargin)
+// Author: Serge Steer, Copyright INRIA
+ H=varargin(1)
+ xdesslib.nyquist
+ if typeof(H)=='rd' then
+ if type(varargin($))==10 then //a comment
+ cm=varargin($);varagin($)=null()
+ else
+ cm=''
+ end
+ D=H.iodelay
+ varargin(1)=H.H
+ [frq,repf,splitf]=repfreq(varargin(1:$))
+ w=-2*%i*%pi*frq;
+ for k=1:size(D,'*')
+ repf(k,:)=repf(k,:).*exp(D(k)*w)
+ end
+ nyquist(frq,repf)
+ else
+ nyquist(varargin(:))
+ end
+endfunction
+
diff --git a/2048/DEPENDENCIES/l2r.sci b/2048/DEPENDENCIES/l2r.sci new file mode 100755 index 000000000..26e61e0c8 --- /dev/null +++ b/2048/DEPENDENCIES/l2r.sci @@ -0,0 +1,13 @@ +// function [Rnum,Rnumdeg,Rden,Rdendeg] = l2r(N,degN,D,degD)
+// given Numerator and Denominator polynomial matrices in left form,
+// not necessarily coprime, finds right coprime factorisation.
+
+function [Rnum,Rnumdeg,Rden,Rdendeg] = l2r(N,degN,D,degD)
+
+[N,degN] = transp(N,degN);
+[D,degD] = transp(D,degD);
+
+[Rnum,Rnumdeg,Rden,Rdendeg] = left_prm(N,degN,D,degD);
+[Rnum,Rnumdeg] = transp(Rnum,Rnumdeg);
+[Rden,Rdendeg] = transp(Rden,Rdendeg);
+endfunction;
diff --git a/2048/DEPENDENCIES/label.sci b/2048/DEPENDENCIES/label.sci new file mode 100755 index 000000000..fd7639958 --- /dev/null +++ b/2048/DEPENDENCIES/label.sci @@ -0,0 +1,23 @@ +// Updated (9-12-06)
+// Input arguments:title, xlabel, ylabel and their font sizes
+
+function label(tname,tfont,labelx,labely,xyfont)
+a = get("current_axes")
+xtitle(tname,labelx,labely)
+xgrid
+t = a.title;
+t.font_size = tfont; //Title font size
+t.font_style = 2; //Title font style
+t.text = tname;
+
+u = a.x_label;
+u.font_size = xyfont; //Label font size
+u.font_style = 2; //Label font style
+
+v = a.y_label;
+v.font_size = xyfont; //Label font size
+v.font_style = 2; //Label font style
+
+// a.label_font_size = 3;
+
+endfunction;
diff --git a/2048/DEPENDENCIES/left_prm.sci b/2048/DEPENDENCIES/left_prm.sci new file mode 100755 index 000000000..cd0ffcbf7 --- /dev/null +++ b/2048/DEPENDENCIES/left_prm.sci @@ -0,0 +1,89 @@ +// function [B,degB,A,degA,Y,degY,X,degX] = ...
+// left_prm(N,degN,D,degD,job,gap)
+//
+// does three different things according to integers that 'job' takes
+// job = 1.
+// this is the default. It is always done for all jobs.
+// -1 -1 -1
+// Given ND , returns coprime B and A where ND = A B
+// It is enough if one sends the first four input arguments
+// If gap is required to be sent, then one can send either 1 or a null
+// entry for job
+// job = 2.
+// first solve for job = 1 and then solve XA + YB = I
+// job = 3.
+// used in solving XD + YN = C
+// after finding coprime factorization, data are returned
+//
+// convention: the variable with prefix deg stand for degrees
+// of the corresponding polynomial matrices
+//
+// input:
+// N: right fraction numerator polynomial matrix
+// D: right fraction denominator polynomial matrix
+// N and D are not neccessarily coprime
+// gap: variable used to zero entries; default value is 1.0e+8
+//
+// output
+// b and A are left coprime num. and den. polynomial matrices
+// X and Y are solutions to Aryabhatta identity, only for job = 2
+
+function [B,degB,A,degA,Y,degY,X,degX] = left_prm(N,degN,D,degD,job,gap)
+if argn(2) == 4 | argn(2) == 5
+ gap = 1.0e8 ;
+end
+// pause
+if argn(2) == 4,
+job = 1; end
+[F,degF] = rowjoin(D,degD,N,degN);
+[Frows,Fbcols] = polsize(F,degF); // Fbcols = block columns
+Fcols = Fbcols * (degF+1) ; // actual columns of F
+T1 = [];pr =[];degT1 = 0; T1rows = 0;shft = 0;
+S=F; sel = ones(Frows,1); T1bcols =1;
+abar = (Fbcols + 1):Frows; // a_super_bar of B-C.Chang
+while isempty(T1) | T1rows < Frows - Fbcols
+ Srows = Frows*T1bcols; // max actual columns of result
+ [T1,T1rows,sel,pr] = ...
+ t1calc(S,Srows,T1,T1rows,sel,pr,Frows,Fbcols,abar,gap);
+ [T1rows,T1cols] = size(T1);
+ if T1rows < Frows - Fbcols
+ T1 = [T1 zeros(T1rows,Frows)];
+ T1bcols = T1bcols + 1; // max. block columns of result
+ degT1 = degT1 + 1; // degree of result
+ shft = shft +Fbcols;
+ S = seshft(S,F,shft);
+ sel = [sel;sel(Srows-Frows+1:Srows)];
+ rowvec = (T1bcols-1)*Frows+(Fbcols+1):T1bcols * Frows;
+ abar = [abar rowvec]; // A_super_bar of B-C.chang
+ end
+end
+
+[B,degB,A,degA] = colsplit(T1,degT1,Fbcols,Frows-Fbcols);
+[B,degB] = clcoef(B,degB);
+B = -B;
+[A,degA] = clcoef(A,degA);
+// pause
+if job == 2
+ S = S(mtlb_logical(sel),:); // columns
+ [redSrows,Scols] = size(S);
+ C = [eye(Fbcols,Fbcols) zeros(Fbcols,Scols-Fbcols)]; // append with zeros
+ T2 = C/S;
+ T2 = makezero(T2,gap);
+ T2 = move_sci(T2,find(sel),Srows);
+ [X,degX,Y,degY] = colsplit(T2,degT1,Fbcols,Frows - Fbcols);
+ [X,degX] = clcoef(X,degX);
+ [Y,degY] = clcoef(Y,degY);
+elseif job == 3
+ Y = S;
+ degY = sel;
+ X = degT1;
+ degX = Fbcols;
+else
+ if job ~= 1
+ error('Message from left_prm:no legal job number specified')
+ end
+end
+endfunction
+
+
+
diff --git a/2048/DEPENDENCIES/makezero.sci b/2048/DEPENDENCIES/makezero.sci new file mode 100755 index 000000000..82644a6ba --- /dev/null +++ b/2048/DEPENDENCIES/makezero.sci @@ -0,0 +1,19 @@ +// function B = makezero(B,gap)
+// where B is a vector and gap acts as a tolerance
+
+function B = makezero(B,gap)
+
+if argn(2) == 1
+ gap = 1.0e8;
+end
+temp = B(find(B)); // non zero entries of B
+temp = -gsort(-abs(temp),'g','d'); // absolute values sorted in descending order
+len = length(temp);
+ratio = temp(1:len-1) ./temp(2:len); // each ratio >1
+min_ind = min(find(ratio>gap));
+if ~isempty(min_ind)
+ our_eps = temp(min_ind+1);
+ zeroind = find(abs(B)<=our_eps);
+ B(zeroind) = zeros(1,length(zeroind));
+end
+endfunction
diff --git a/2048/DEPENDENCIES/move_sci.sci b/2048/DEPENDENCIES/move_sci.sci new file mode 100755 index 000000000..6bd822942 --- /dev/null +++ b/2048/DEPENDENCIES/move_sci.sci @@ -0,0 +1,13 @@ +// function result = move_sci(b,nonred,max_sci)
+// Moves matrix b to matrix result with the information on where to move,
+// decided by the indices of nonred.
+// The matrix result will have as many rows as b has and max number of columns.
+// b is augumented with zeros to have nonred number of columns;
+// The columns of b put into those of result as decided by nonred.
+
+function result = move_sci(b,nonred,max_sci)
+[brows,bcols] = size(b);
+b = [b zeros(brows,length(nonred)-bcols)];
+result = zeros(brows,max_sci);
+result(:,nonred') = b;
+endfunction
diff --git a/2048/DEPENDENCIES/oe.sci b/2048/DEPENDENCIES/oe.sci new file mode 100755 index 000000000..767bd39d0 --- /dev/null +++ b/2048/DEPENDENCIES/oe.sci @@ -0,0 +1,146 @@ +
+// OE Model parameter estimation
+
+///////////////////////////////////////
+/////////// ARX Model /////////////////
+
+
+function [thetaN_arx,covt_arx,nvar,res] = arxc(data,na,nb,nk)
+az = max(na,nb+nk-1);
+zer = zeros(az,1);
+zd = data;
+// Zeros appended
+zd1(:,1) = [zer; zd(:,1)];
+zd1(:,2) = [zer; zd(:,2)];
+[r,c] = size(zd1);
+t = az+1:r;
+yt = zd1(:,1); ut = zd1(:,2);
+yt1 = yt'; ut1 = ut'; // row vector
+len1 = length(yt1);
+yt2 = zeros(1,len1-az); ut2 = zeros(1,len1-az);
+
+// arx(Data,[na nb nk])
+ for i=1:na
+ yt2 = [yt2; -yt1(t-i)];
+ end;
+ for i=nk:nb+nk-1
+ ut2 = [ut2; ut1(t-i)];
+ end;
+[r1,c1] = size(yt2); [r2,c2] = size(ut2);
+phit = [yt2(2:r1,:); ut2(2:r2,:)];
+m1 = phit*phit';
+[qm,rm] = qr(m1);
+m2 = phit*zd(:,1);
+thetaN_arx = inv(rm)*qm'*m2;
+// thetaN_arx = inv(m1)*m2;
+// thetaN_arx = m1\m2;
+
+[r11,c11] = size(thetaN_arx);
+a = thetaN_arx(1:na); b = thetaN_arx(na+1:r11);
+
+// Sum of squared residuals
+
+yhat = phit'*thetaN_arx;
+res = zd(:,1) - yhat;
+N = length(res);
+q = rank(phit);
+ssr = res'*res;
+sig2 = ssr/(N-q);
+nvar = sqrt(sig2);
+cov_arx = inv(m1);
+covt_arx = diag(cov_arx);
+endfunction;
+///////////////////////////////
+//////////////////////////////
+
+///////////////////////////////////////
+////////// Model Display /////////////
+
+function disp_mod(N1,covN1)
+len = length(covN1);
+B1 = pol2str(N1);
+ind = strindex(B1,['+','-']);
+ind = ind - 1;
+B2 = strsplit(B1,ind);
+covB = string(covN1);
+
+ if ascii(B2(1)) == 32
+ B2 = B2(2:len+1);
+ end;
+
+ B3(1) = ' ';
+ for i=1:len
+ B3(i) = strsubst(B2(i),'*x','(+-' + covB(i) + ')*x');
+ end;
+
+ B4 = B3(1);
+
+ for i=2:len
+ B4 = B4 + ' ' + B3(i);
+ end;
+
+disp(B4);
+endfunction;
+///////////////////////////////////////
+///////////////////////////////////////
+
+function [thetaN_oe,covN_oe,nvar,resid] = oe(zd,nb,nf,nk)
+
+[thetaN,covfN,nvar,res] = arxc(zd,nf,nb,nk);
+[r1,c1] = size(thetaN);
+yt = zd(:,1);
+m=50;
+ if nf==0
+ thetaN_oe = thetaN;
+ covN_oe = covfN;
+ else
+ for k=1:m
+ a = thetaN(1:nf);
+ b = thetaN(nf+1:r1);
+ A = [1 a']; // Filter
+ y = yt(1:length(u))';
+ yf = deconvol(y,A);
+ uf = deconvol(u,A);
+ zf = [yf(1:length(uf))' uf'];
+ zdf = detrend(zf,'constant');
+ [thetaNf,covf_a,nvar,resid] = arxc(zdf,nf,nb,nk)
+ thetaN = thetaNf;
+ a1 = thetaN(1:nf);
+ b = (norm(a-a1))/norm(a1);
+ if b<0.005
+ break;
+ end;
+ end;
+ thetaN_oe = thetaN;
+ covN_oe = covf_a;
+ end;
+[rt,ct] = size(thetaN_oe);
+f_oe = [1 thetaN_oe(1:nf)']; b1 = zeros(1,nk);
+b_oe = [b1 thetaN_oe(nf+1:rt)'];
+cov_f = covN_oe(1:nf); cov_b = covN_oe(nf+1:rt);
+
+x = poly(0,'x');
+ if nf ==0
+ disp('Discrete time model: y(t) = B(x)u(t) + e(t)');
+ else
+ disp('Discrete time model: y(t) = [B(x)/F(x)]u(t) + e(t)');
+ end;
+
+F = poly(f_oe,'x','coeff');
+cov_f1 = [0 cov_f'];
+B = poly( b_oe,'x','coeff');
+cov_b1 = cov_b';
+
+ if nb==0
+ error('All B parameters are zero');
+ else
+ disp('B(x) = ');
+ disp_mod(B,cov_b);
+ end;
+
+ if nf~=0
+ disp('F(x) = ');
+ disp_mod(F,cov_f1);
+ end;
+endfunction;
+
diff --git a/2048/DEPENDENCIES/plotacf.sci b/2048/DEPENDENCIES/plotacf.sci new file mode 100755 index 000000000..a36f8caff --- /dev/null +++ b/2048/DEPENDENCIES/plotacf.sci @@ -0,0 +1,39 @@ +// Procedure to plot the ACF, as discussed in Sec. 6.4.3. An example usage is given in 6.5.
+// 6.6
+
+// PLOTACF: Plots normalized autocorrelation function
+//
+// USAGE:: [acf]=plotacf(x,errlim,len,print_code)
+//
+// WHERE:: acf = autocorrelation values
+// x = time series data
+// errlim > 0; error limit = 2/sqrt(data_len)
+// len = length of acf that need to to be plotted
+// NOTE: if len=0 then len=data_length/2;
+// print_code = 0 ==> does not plot OR ELSE plots
+//
+// Pranob Banerjee
+
+function [x]=plotacf(y,errlim,len,code)
+exec('normacf.sci',-1);
+exec('label.sci',-1)
+x = normacf(y);
+l = length(y);
+
+r=l:2*(l-1); lim=2/sqrt(l); rl=1:length(r) ;
+N=length(rl); x=x(r);
+if len>0 & len<N, rl=1:len; x=x(rl); N=len; end;
+if(code > 0 )
+ if(errlim > 0 )
+ rl=rl-1;
+ plot(rl,x,rl,x,'o' , rl,lim*ones(N,1),'--', ...
+ rl,-lim*ones(N,1),'--')
+ xgrid
+ else
+ plot(rl,x)
+ end
+end;
+a = gca();
+a.data_bounds = [0 min(min(x),-lim-0.1); len-1 1.1];
+label(' ',4,'Lag','ACF',4)
+endfunction;
diff --git a/2048/DEPENDENCIES/pol2cart.sci b/2048/DEPENDENCIES/pol2cart.sci new file mode 100755 index 000000000..b12dbc32e --- /dev/null +++ b/2048/DEPENDENCIES/pol2cart.sci @@ -0,0 +1,6 @@ +// Polar to cartesian co ordinate conversion
+
+function [x,y] = pol2cart(theta,r)
+x = r*cos(theta);
+y = r*sin(theta);
+endfunction;
diff --git a/2048/DEPENDENCIES/poladd.sci b/2048/DEPENDENCIES/poladd.sci new file mode 100755 index 000000000..79fcdabed --- /dev/null +++ b/2048/DEPENDENCIES/poladd.sci @@ -0,0 +1,16 @@ +function [C,degC] = poladd(A,degA,B,degB)
+[rA,cA] = polsize(A,degA);
+[rB,cB] = polsize(B,degB);
+if cA ~= cB | rA ~= rB
+ error('poladd: Inconsistent dimensions');
+end
+
+degC = max(degA,degB);
+if degC >= degA
+ A = [A zeros(rA,(degC-degA)*cA)];
+end
+if degC >= degB
+ B = [B zeros(rB,(degC-degB)*cB)];
+end
+C = A+B;
+endfunction;
diff --git a/2048/DEPENDENCIES/polmul.sci b/2048/DEPENDENCIES/polmul.sci new file mode 100755 index 000000000..67e4adb9a --- /dev/null +++ b/2048/DEPENDENCIES/polmul.sci @@ -0,0 +1,34 @@ +// polmul
+// The command
+// [C,degA] = polmul(A,degA,B,degB)
+// produces the polynomial matrix C that equals the product A*B of the
+// polynomial matrices A and B.
+//
+// H. Kwakernaak, July, 1990
+
+
+function [C,degC] = polmul(A,degA,B,degB)
+[rA,cA] = polsize(A,degA);
+[rB,cB] = polsize(B,degB);
+if cA ~= rB
+ error('polmul: Inconsistent dimensions of input matrices');
+end
+
+degC = degA+degB;
+C = [];
+for k = 0:degA+degB
+ mi = 0;
+ if k-degB > mi
+ mi = k-degB;
+ end
+ ma = degA;
+ if k < ma
+ ma = k;
+ end
+ Ck = zeros(rA,cB);
+ for i = mi:ma
+ Ck = Ck + A(:,i*cA+1:(i+1)*cA)*B(:,(k-i)*cB+1:(k-i+1)*cB);
+ end
+ C = [C Ck];
+end
+endfunction
diff --git a/2048/DEPENDENCIES/polsize.sci b/2048/DEPENDENCIES/polsize.sci new file mode 100755 index 000000000..7c0c8a463 --- /dev/null +++ b/2048/DEPENDENCIES/polsize.sci @@ -0,0 +1,15 @@ +// function [rQ,cQ] = polsize(Q,degQ)
+// FUNCTION polsize TO DETERMINE THE DIMENSIONS
+// OF A POLYNOMIAL MATRIX
+//
+// H. Kwakernaak, August, 1990
+
+function [rQ,cQ] = polsize(Q,degQ)
+
+[rQ,cQ] = size(Q); cQ = cQ/(degQ+1);
+if abs(round(cQ)-cQ) > 1e-6
+ error('polsize: Degree of input inconsistent with number of columns');
+else
+ cQ = round(cQ);
+end
+endfunction
diff --git a/2048/DEPENDENCIES/polyno.sci b/2048/DEPENDENCIES/polyno.sci new file mode 100755 index 000000000..e54960cdb --- /dev/null +++ b/2048/DEPENDENCIES/polyno.sci @@ -0,0 +1,30 @@ +// Updated(1-8-07)
+// Operations:
+// Polynomial definition
+// Flipping of coefficients
+// Variables ------- passed as input argument (either 's' or 'z')
+// Both num and den are used mostly used in scicos files,
+// to get rid of negative powers of z
+
+// Polynomials with powers of s need to
+// be flipped only
+
+function [polynu,polyde] = polyno(zc,a)
+zc = clean(zc);
+polynu = poly(zc(length(zc):-1:1),a,'coeff');
+ if a == 'z'
+ polyde = %z^(length(zc) - 1);
+ else
+ polyde = 1;
+ end
+
+// Scicos(4.1) Filter block shouldn't have constant/constant
+ if type(polynu)==1 & type(polyde)==1
+ if a == 'z'
+ polynu = %z; polyde = %z;
+ else
+ polynu = %s; polyde = %s;
+ end;
+ end;
+
+endfunction
diff --git a/2048/DEPENDENCIES/putin.sci b/2048/DEPENDENCIES/putin.sci new file mode 100755 index 000000000..706898856 --- /dev/null +++ b/2048/DEPENDENCIES/putin.sci @@ -0,0 +1,22 @@ +// function [A,degA] = putin(A,degA,B,degB,i,j)
+// PUTS IN THE POLYNOMIAL B INTO THE MATRIX A AT (i,j)TH POSITION
+// Modified by Kannan Moudgalya on 8 November 1992
+
+function [A,degA] = putin(A,degA,B,degB,i,j)
+
+[rA,cA] = polsize(A,degA);
+if degB > degA
+ A = [A zeros(rA,(degB-degA)*cA)];
+ degA = degB;
+end
+
+for k = 0:degB
+ A(i,(k*cA)+j) = B(1,k+1);
+end
+if degA > degB
+ for k = degB+1:degA
+ A(i,(k*cA)+j) = 0;
+ end
+ [A,degA] = clcoef(A,degA);
+end
+endfunction;
diff --git a/2048/DEPENDENCIES/respol.sci b/2048/DEPENDENCIES/respol.sci new file mode 100755 index 000000000..c62c7e05e --- /dev/null +++ b/2048/DEPENDENCIES/respol.sci @@ -0,0 +1,112 @@ +// Computation of residues
+// 4.5
+// Numerator and denominator coefficients
+// are passed in decreasing powers of z(say)
+
+function [res,pol,q] = respol(num,den)
+len = length(num);
+if num(len) == 0
+ num = num(1:len-1);
+end
+
+[resi,q] = pfe(num,den);
+res = resi(:,2);
+res = int(res) + (clean(res - int(res),1.d-04));
+pol = resi(:,1);
+pol = int(pol) + (clean(pol - int(pol),1.d-04));
+endfunction;
+
+///////////////////////////////////////////////////
+// Partial fraction expansion
+
+function [resid1,q] = pfe(num,den)
+x = poly(0,'x');
+s = %s;
+
+exec flip.sci;
+
+num2 = flip(num);
+den2 = flip(den);
+num = poly(num2,'s','coeff');
+den = poly(den2,'s','coeff');
+[fac,g] = factors(den);
+polf = polfact(den);
+n = 1;
+
+r = clean(real(roots(den)),1.d-5);
+//disp(r);
+l = length(r);
+r = gsort(r,'g','i');
+r = [r; 0];
+j = 1;
+t1 = 1; q = [];
+rr = 0;
+rr1 = [0 0];
+m = 1;
+
+ for i = j:l
+ if abs(r(i)- r(i+1)) < 0.01 then
+ r(i);
+ r(i+1);
+ n = n+1;
+ m = n;
+ //pause
+ t1 = i;
+ //disp('Repeated roots')
+ else
+ m = n;
+ //pause
+ n = 1;
+ end
+ i;
+ if n == 1 then
+ rr1 = [rr1; r(i) m];
+ //pause
+ end;
+ j = t1 + 1;
+ end;
+rr2 = rr1(2:$,:);
+[r1,c1] = size(rr2);
+den1 = 1;
+
+for i = 1:r1
+ den1 = den1 * ((s-rr2(i,1))^(rr2(i,2)));
+end;
+[rem,quo] = pdiv(num,den);
+q = quo;
+if quo ~= 0
+ num = rem;
+end
+
+tf = num/den1;
+res1 = 0;
+res3 = [s 0];
+res5 = [0 0];
+for i = 1:r1
+ j = rr2(i,2) + 1;
+ tf1 = tf; //strictly proper
+ k = rr2(i,2);
+ tf2 = ((s-rr2(i,1))^k)*tf1;
+ rr2(i,1);
+ res1 = horner(tf2,rr2(i,1));
+ res2 = [(s - rr2(i,1))^(rr2(i,2)) res1];
+ res4 = [rr2(i,1) res1];
+ res3 = [res3; res2];
+ res5 = [res5; res4];
+ res = res1;
+ for m = 2:j-1
+ k;
+ rr2(i,1);
+ tf1 = derivat(tf2)/factorial(m-1); //ith derivative
+ res = horner(tf1,rr2(i,1));
+ res2 = [(s - rr2(i,1))^(j-m) res];
+ res4 = [rr2(i,1) res];
+ res5 = [res5; res4];
+ res3 = [res3; res2];
+ tf2 = tf1;
+ end;
+end;
+resid = res3(2:$,:); //with s terms
+resid1 = res5(2:$,:); //only poles(in decreasing no. of repetitions)
+endfunction;
+////////////////////////////////////////////////////////////
diff --git a/2048/DEPENDENCIES/rowjoin.sci b/2048/DEPENDENCIES/rowjoin.sci new file mode 100755 index 000000000..c305183cf --- /dev/null +++ b/2048/DEPENDENCIES/rowjoin.sci @@ -0,0 +1,31 @@ +// function [P,degP] = rowjoin(P1,degP1,P2,degP2)
+// MATLAB FUNCTION rowjoin TO SUPERPOSE TWO POLYNOMIAL
+// MATRICES
+
+// H. Kwakernaak, July, 1990
+
+function [P,degP] = rowjoin(P1,degP1,P2,degP2)
+
+[rP1,cP1] = polsize(P1,degP1);
+[rP2,cP2] = polsize(P2,degP2);
+if cP1 ~= cP2
+ error('rowjoin: Inconsistent numbers of columns');
+end
+
+rP = rP1+rP2; cP = cP1;
+if degP1 >= degP2
+ degP = degP1;
+else
+ degP = degP2;
+end
+
+if isempty(P1)
+ P = P2;
+elseif isempty(P2)
+ P = P1;
+else
+ P = zeros(rP,(degP+1)*cP);
+ P(1:rP1,1:(degP1+1)*cP1) = P1;
+ P(rP1+1:rP,1:(degP2+1)*cP2) = P2;
+end
+endfunction
diff --git a/2048/DEPENDENCIES/seshft.sci b/2048/DEPENDENCIES/seshft.sci new file mode 100755 index 000000000..e69cdf5b3 --- /dev/null +++ b/2048/DEPENDENCIES/seshft.sci @@ -0,0 +1,22 @@ +// function C = seshft(A,B,N)
+//given A and B matrices, returns C = [<-A-> 0
+// 0 <-B->] with B shifted east by N cols
+
+function C = seshft(A,B,N)
+[Arows,Acols] = size(A);
+[Brows,Bcols] = size(B);
+if N >= 0
+ B = [zeros(Brows,N) B];
+ Bcols = Bcols + N;
+elseif N < 0
+ A = [zeros(Arows,abs(N)) A];
+ Acols = Acols +abs(N);
+end
+ if Acols < Bcols
+ A = [A zeros(Arows,Bcols-Acols)];
+ elseif Acols > Bcols
+ B = [B zeros(Brows,Acols-Bcols)];
+ end
+ C = [A
+ B];
+endfunction
diff --git a/2048/DEPENDENCIES/stem.sci b/2048/DEPENDENCIES/stem.sci new file mode 100755 index 000000000..3042267fb --- /dev/null +++ b/2048/DEPENDENCIES/stem.sci @@ -0,0 +1,20 @@ +// Stem plot
+// Updated (1-12-06)
+function stem(x,y,xy)
+if argn(2) == 2
+xy = 5;
+end;
+set("figure_style","new")
+plot2d3(x,y,style=2);
+//plot2d4(x,y,style=-9); // default mark foreground colour: black
+// Can be used instead of xpoly
+// But default values cannot be changed
+xpoly(x,y,"marks");
+p=get("hdl");
+p.mark_size_unit = 'point';
+p.mark_size = xy;
+p.mark_style = 9;
+p.mark_foreground = 2;
+p.mark_background = 2;
+endfunction;
+
diff --git a/2048/DEPENDENCIES/t1calc.sci b/2048/DEPENDENCIES/t1calc.sci new file mode 100755 index 000000000..a9a265ba9 --- /dev/null +++ b/2048/DEPENDENCIES/t1calc.sci @@ -0,0 +1,35 @@ +// function [T1,T1rows,sel,pr] = ...
+// t1calc(S,Srows,T1,T1rows,sel,pr,Frows,Fbcols,abar,gap)
+// calculates the coefficient matrix T1
+// redundant row information is kept in sel: redundant rows are marked
+// with zeros. The undeleted rows are marked with ones.
+
+function [T1,T1rows,sel,pr] = t1calc(S,Srows,T1,T1rows,sel,pr,Frows,Fbcols,abar,gap)
+b = 1; // vector of primary red.rows
+
+while (T1rows < Frows - Fbcols) & or(sel==1) & ~isempty(b)
+ S = clean(S);
+ b = indep(S(mtlb_logical(sel),:),gap); // send selected rows of S
+ if ~isempty(b)
+ b = clean(b);
+ b = move_sci(b,find(sel),Srows);
+ j = length(b);
+ while ~(b(j) & or(abar==j)) // pick largest nonzero entry
+ j = j-1; // of coeff. belonging to abar
+ if ~j
+ fprintf('\nMessage from t1calc, called from left_prm\n\n')
+ error('Denominator is noninvertible')
+ end
+ end
+ if ~or(j<pr & pmodulo(pr,Frows) == pmodulo(j,Frows)) // pr(2),pr(1)
+ T1 = [T1; b]; // condition is not violated
+ T1rows = T1rows +1; // accept this vector
+ end // else don't accept
+ pr = [pr; j]; // update prime red row info
+ while j <= Srows
+ sel(j) = 0;
+ j = j + Frows;
+ end
+ end
+end
+endfunction
diff --git a/2048/DEPENDENCIES/tf.sci b/2048/DEPENDENCIES/tf.sci new file mode 100755 index 000000000..bba921d14 --- /dev/null +++ b/2048/DEPENDENCIES/tf.sci @@ -0,0 +1,33 @@ +//User defined function
+//Forms a transfer function
+//Scilab: Co efficients are given in increasing power of variable
+//Matlab: Co efficients are given in decreasing power of variable
+//Hence co efficients are flipped here
+
+//Input arguments: (1) Numerator co efficients(decreasing order)
+//(2) Denominator co efficients
+//(3) Variable to specify domain
+
+// Updated (30-11-06)
+// System is continuous => a is not passed
+// System is discrete => a = -1
+// System is discretized (sampled system) => a = Ts
+// Uses syslin
+
+function trfu = tf(num,den,a)
+ if argn(2) == 2
+ d = 'c';
+ elseif a == -1
+ d = 'd';
+ else
+ d = a
+ end;
+num = clean(num);
+den = clean(den);
+num1 = poly(num(length(num):-1:1),'x','coeff');
+den1 = poly(den(length(den):-1:1),'x','coeff');
+trfu = syslin(d,num1,den1);
+endfunction;
+
+
+
diff --git a/2048/DEPENDENCIES/transp.sci b/2048/DEPENDENCIES/transp.sci new file mode 100755 index 000000000..39854c9ce --- /dev/null +++ b/2048/DEPENDENCIES/transp.sci @@ -0,0 +1,17 @@ +// function [P,degP] = transp(Q,degQ)
+// MATLAB FUNCTION transp TO TRANSPOSE
+// A POLYNOMIAL MATRIX
+
+// H. Kwakernaak, July, 1990
+
+function [P,degP] = transp(Q,degQ)
+
+[rQ,cQ] = polsize(Q,degQ);
+
+rP = cQ; cP = rQ; degP = degQ;
+P = zeros(rP,(degP+1)*cP);
+for i = 1:degP+1
+ P(:,(i-1)*cP+1:i*cP) = Q(:,(i-1)*cQ+1:i*cQ)';
+end
+
+endfunction;
diff --git a/2048/DEPENDENCIES/xdync.sci b/2048/DEPENDENCIES/xdync.sci new file mode 100755 index 000000000..e082e587f --- /dev/null +++ b/2048/DEPENDENCIES/xdync.sci @@ -0,0 +1,37 @@ +// function [Y,degY,X,degX,B,degB,A,degA] = xdync(N,degN,D,degD,C,degC,gap)
+// given coefficient matrix in T1, primary redundant row information sel,
+// solves XD + YN = C
+
+// calling order changed on 16 April 2005. Old order:
+// function [B,degB,A,degA,Y,degY,X,degX] = xdync(N,degN,D,degD,C,degC,gap)
+
+function [Y,degY,X,degX,B,degB,A,degA] = xdync(N,degN,D,degD,C,degC,gap)
+if argn(2) == 6
+ gap = 1.0e+8;
+end
+
+[F,degF] = rowjoin(D,degD,N,degN);
+
+[Frows,Fbcols] = polsize(F,degF); //Fbcols = block columns
+
+[B,degB,A,degA,S,sel,degT1,Fbcols] = left_prm(N,degN,D,degD,3,gap);
+//if issoln(D,degD,C,degC,B,degB,A,degA)
+ [Crows,Ccols] = size(C);
+ [Srows,Scols] = size(S);
+ S = clean(S);
+ S = S(mtlb_logical(sel),:);
+ T2 =[];
+
+ for i = 1:Crows,
+ Saug = seshft(S,C(i,:),0);
+ b = cindep(Saug);
+ b = move_sci(b,find(sel),Srows);
+ T2 =[T2; b];
+ end
+
+ [X,degX,Y,degY] = colsplit(T2,degT1,Fbcols,Frows-Fbcols);
+
+ [X,degX] = clcoef(X,degX);
+ [Y,degY] = clcoef(Y,degY);
+ Y = clean(Y); X = clean(X);
+endfunction
|