summaryrefslogtreecommitdiff
diff options
context:
space:
mode:
Diffstat
-rwxr-xr-xGet.sci113
-rw-r--r--ar.sci102
-rw-r--r--armaX.sci51
-rw-r--r--arx.sci110
-rw-r--r--bj.sci122
-rw-r--r--calModelPara.sci16
-rw-r--r--compare.sci145
-rwxr-xr-xcompareBJ.sci56
-rw-r--r--dataSlice.sci67
-rw-r--r--detrend.sci100
-rw-r--r--estpoly.sci136
-rw-r--r--etfe.sci29
-rw-r--r--etfeeTest.sce13
-rw-r--r--fitch.sci42
-rw-r--r--frd.sci84
-rw-r--r--frdplot.sci11
-rwxr-xr-xiddata.sci68
-rw-r--r--iddataplot.sci69
-rw-r--r--identTime.sci32
-rwxr-xr-xidpoly.sci142
-rw-r--r--impulseest.sci78
-rw-r--r--iv.sci98
-rw-r--r--iv4.sci75
-rw-r--r--misData.matbin0 -> 113297 bytes
-rw-r--r--misdata.sci41
-rw-r--r--nInputData.sci15
-rw-r--r--oe.sci118
-rw-r--r--pe.sci112
-rw-r--r--predict.sci192
-rw-r--r--rARX.matbin0 -> 20743 bytes
-rw-r--r--rData.matbin0 -> 2319 bytes
-rw-r--r--rarx.sci63
-rw-r--r--sim.sci65
-rw-r--r--spa.sci73
-rw-r--r--stem.sci19
-rw-r--r--stepest.sci71
36 files changed, 2528 insertions, 0 deletions
diff --git a/Get.sci b/Get.sci
new file mode 100755
index 0000000..f14851a
--- /dev/null
+++ b/Get.sci
@@ -0,0 +1,113 @@
+function []=Get(varargin)
+ //disp('Get')
+ tempData = varargin($)
+ f = fieldnames(tempData)
+// disp(size(f,'*'))
+ if ~size(f,'*') then
+ if size(tempData,'*') == 0 then
+ printString = 'empty'
+ elseif isnan(tempData) then
+ printString = 'Nan'
+ elseif isinf(tempData) then
+ printString = 'inf'
+ else
+ printString = typeof(tempData)
+ end
+ mprintf('The type of ""Input"" is : %s \n',printString)
+ error(msprintf(gettext("%s: There is no objects in the ""Input"".\n"),"Get"))
+ else
+ maxLength = []
+ for ii = 1:size(f,'*')
+ maxLength = [maxLength length(f(ii))]
+ end
+ maxLength = max(maxLength)
+ for ii = 1:size(f,'*')
+ blanckSpace = ' '
+ for jj = 1:maxLength-length(f(ii))
+ blanckSpace = blanckSpace + ' '
+ end
+ mprintf('\t%s%s : ',blanckSpace,f(ii))
+ objectData = tempData(f(ii))
+ if typeof(objectData) == 'ce' then
+ //mprintf('cell')
+ sizeData = size(objectData)
+ prodData = prod(sizeData)
+ //pause
+ if prodData == 1 then
+ objectData = (cell2mat(objectData))
+ objectData = getString(objectData)
+ tempString = '{' + objectData + '}'
+ mprintf('%s',tempString)
+ else
+ objectData = getString(sizeData,[])
+ objectData = strsubst(objectData,']',' cell ]')
+ mprintf('%s',objectData)
+ end
+ elseif typeof(objectData) == 'string' then
+ if length(objectData) == 0 then
+ mprintf("{''''}")
+ else
+ mprintf('%s ',objectData)
+ end
+
+ elseif typeof(objectData) == 'constant' then
+ //disp('in constant')
+ //
+ if ~size(objectData,'*') then
+ mprintf('%s%s','[',']')
+ //pause
+ //mprintf('%d',objectData)
+ elseif size(objectData,'*') == 1 then
+ if round(objectData)-objectData == 0 then
+ mprintf('%d',objectData)
+ else
+ mprintf('%.2f',objectData)
+ end
+ else
+ sizeData = size(objectData)
+ objectData = getString(sizeData,[])
+ objectData = strsubst(objectData,']',' double ]')
+ mprintf('%s',objectData)
+ end
+ elseif typeof(objectData) == 'polynomial' then
+ if size(objectData,'*') == 1 then
+ polyData = pol2str(objectData)
+ mprintf('%s',polyData)
+ else
+ objectData = size(objectData)
+ objectData = getString(objectData,[])
+ objectData = strsubst(objectData,']',' polynomial ]')
+ mprintf('%s',objectData)
+ end
+ elseif typeof(objectData) == 'hypermat' then
+ objectData = size(objectData)
+ objectData = getString(objectData,[])
+ objectData = strsubst(objectData,']',' hypermat ]')
+ mprintf('%s',objectData)
+ end
+ mprintf('\n')
+ end
+ end
+endfunction
+
+
+// converting from matrix to string matrix
+function Output = getString(varargin)
+ [lhs,rhs] = argn(0)
+ if rhs == 1 then
+ inBetween = ' '
+ elseif rhs == 2 then
+ inBetween = 'x'
+ end
+ stringData = '['
+ matData = varargin(1)
+ for ii = 1 : length(matData)
+ stringData = stringData + string(matData(ii)) //+ inBetween
+ if ii == length(matData) then
+ else
+ stringData = stringData + inBetween
+ end
+ end
+ stringData = stringData + ']'
+ Output = stringData
+endfunction
diff --git a/ar.sci b/ar.sci
new file mode 100644
index 0000000..4239834
--- /dev/null
+++ b/ar.sci
@@ -0,0 +1,102 @@
+// Estimates Discrete time AR model
+// A(q)y(t) = e(t)
+// Current version uses random initial guess
+
+// Authors: Ashutosh,Harpreet,Inderpreet
+// Updated(12-6-16)
+function sys = ar(varargin)
+//
+ [lhs , rhs] = argn();
+ if ( rhs < 2 ) then
+ errmsg = msprintf(gettext("%s: Unexpected number of input arguments : %d provided while should be 2"), "ar", rhs);
+ error(errmsg)
+ end
+
+ z = varargin(1)
+ if typeof(z) == 'iddata' then
+ Ts = z.Ts;unit = z.TimeUnit
+ z = [z.OutputData z.InputData]
+ elseif typeof(z) == 'constant' then
+ Ts = 1;unit = 'seconds'
+ end
+ if ~iscolumn(z) then
+ errmsg = msprintf(gettext("%s: time series output data only"), "ar");
+ error(errmsg);
+ end
+
+ if (~isreal(z)) then
+ errmsg = msprintf(gettext("%s: input and output data matrix should be a real matrix"), "ar");
+ error(errmsg);
+ end
+
+ n = varargin(2)
+ if (size(n,"*") ~=1 )then
+ errmsg = msprintf(gettext("%s: order should be nonnegative integer number "), "ar");
+ error(errmsg);
+ end
+
+ if (size(find(n<0),"*") | size(find(((n-floor(n))<%eps)== %f))) then
+ errmsg = msprintf(gettext("%s: values of order and delay matrix [na] should be nonnegative integer number "), "ar");
+ error(errmsg);
+ end
+
+ na = n; nb = 0; nk = 0;
+ // storing U(k) , y(k) and n data in UDATA,YDATA and NDATA respectively
+ YDATA = z(:,1);
+ UDATA = zeros(size(z,1),1)
+ NDATA = size(UDATA,"*");
+ function e = G(p,m)
+ e = YDATA - _objfun(UDATA,YDATA,p,na,nb,nk);
+ endfunction
+ tempSum = na+nb
+ p0 = linspace(0.1,0.9,tempSum)';
+ [var,errl] = lsqrsolve(p0,G,size(UDATA,"*"));
+ err = (norm(errl)^2);
+ opt_err = err;
+ resid = G(var,[]);
+ a = 1-poly([var(nb+1:nb+na)]',"q","coeff");
+ b = poly([repmat(0,nk,1);var(1:nb)]',"q","coeff");
+ a = (poly([1,-coeff(a)],'q','coeff'))
+ t = idpoly(coeff(a),1,1,1,1,Ts)
+
+ // estimating the other parameters
+ [temp1,temp2,temp3] = predict([YDATA UDATA],t)
+ [temp11,temp22,temp33] = pe([YDATA UDATA],t)
+
+ estData = calModelPara(temp1,temp1,n(1))
+ //pause
+ t.Report.Fit.MSE = estData.MSE
+ t.Report.Fit.FPE = estData.FPE
+ t.Report.Fit.FitPer = estData.FitPer
+ t.Report.Fit.AIC = estData.AIC
+ t.Report.Fit.AICc = estData.AICc
+ t.Report.Fit.nAIC = estData.nAIC
+ t.Report.Fit.BIC = estData.BIC
+ t.TimeUnit = unit
+ sys = t
+ //sys = idpoly(coeff(a),1,1,1,1,Ts)
+// sys.TimeUnit = unit
+endfunction
+
+function yhat = _objfun(UDATA,YDATA,x,na,nb,nk)
+ x=x(:)
+ q = poly(0,'q')
+ tempSum = nb+na
+ // making polynomials
+ b = poly([repmat(0,nk,1);x(1:nb)]',"q","coeff");
+ a = 1 - poly([x(nb+1:nb+na)]',"q","coeff")
+ aSize = coeff(a);bSize = coeff(b)
+ maxDelay = max([length(aSize) length(bSize)])
+ yhat = [YDATA(1:maxDelay)]
+ for k=maxDelay+1:size(UDATA,"*")
+ tempB = 0
+ for ii = 1:size(bSize,'*')
+ tempB = tempB + bSize(ii)*UDATA(k-ii+1)
+ end
+ tempA = 0
+ for ii = 1:size(aSize,"*")
+ tempA = tempA + aSize(ii)*YDATA(k-ii)
+ end
+ yhat = [yhat; [ tempA+tempB ]];
+ end
+endfunction
diff --git a/armaX.sci b/armaX.sci
new file mode 100644
index 0000000..93376e2
--- /dev/null
+++ b/armaX.sci
@@ -0,0 +1,51 @@
+
+function sys = armaX(varargin)
+ [lhs , rhs] = argn();
+ if ( rhs < 2 ) then
+ errmsg = msprintf(gettext("%s: Unexpected number of input arguments : %d provided while should be 2"), "armaX", rhs);
+ error(errmsg)
+ end
+
+ z = varargin(1)
+ if typeof(z) == 'iddata' then
+ Ts = z.Ts;unit = z.TimeUnit
+ z = [z.OutputData z.InputData]
+ elseif typeof(z) == 'constant' then
+ Ts = 1;unit = 'seconds'
+ end
+ if ((~size(z,2)==2) & (~size(z,1)==2)) then
+ errmsg = msprintf(gettext("%s: input and output data matrix should be of size (number of data)*2"), "armaX");
+ error(errmsg);
+ end
+
+ if (~isreal(z)) then
+ errmsg = msprintf(gettext("%s: input and output data matrix should be a real matrix"), "armaX");
+ error(errmsg);
+ end
+
+ n = varargin(2)
+ if (size(n,"*")<3| size(n,"*")>4) then
+ errmsg = msprintf(gettext("%s: The order and delay matrix [na nb nc nk] should be of size [3 or 4]"), "armaX");
+ error(errmsg);
+ end
+
+ if (size(find(n<0),"*") | size(find(((n-floor(n))<%eps)== %f))) then
+ errmsg = msprintf(gettext("%s: values of order and delay matrix [na nb nc nk] should be nonnegative integer number "), "armaX");
+ error(errmsg);
+ end
+
+ na = n(1); nb = n(2); nc = n(3); //nd = n(4);nf = n(5);
+
+ if (size(n,"*") == 3) then
+ nk = 1
+ else
+ nk = n(4);
+ end
+
+ // storing U(k) , y(k) and n data in UDATA,YDATA and NDATA respectively
+ YDATA = z(:,1);
+ UDATA = z(:,2);
+
+ sys = estpoly(z,[na,nb,nc,0,0,nk])
+
+endfunction
diff --git a/arx.sci b/arx.sci
new file mode 100644
index 0000000..0289cc8
--- /dev/null
+++ b/arx.sci
@@ -0,0 +1,110 @@
+
+// Estimates Discrete time ARX model
+// A(q)y(t) = B(q)u(t) + e(t)
+// Current version uses random initial guess
+//
+
+// Authors: Ashutosh,Harpreet,Inderpreet
+// Updated(12-6-16)
+function sys = arx(varargin)
+ [lhs , rhs] = argn();
+ if ( rhs < 2 ) then
+ errmsg = msprintf(gettext("%s: Unexpected number of input arguments : %d provided while should be 2"), "arx", rhs);
+ error(errmsg)
+ end
+
+ z = varargin(1)
+ if typeof(z) == 'iddata' then
+ Ts = z.Ts;unit = z.TimeUnit
+ z = [z.OutputData z.InputData]
+ elseif typeof(z) == 'constant' then
+ Ts = 1;unit = 'seconds'
+ end
+ if ((~size(z,2)==2) & (~size(z,1)==2)) then
+ errmsg = msprintf(gettext("%s: input and output data matrix should be of size (number of data)*2"), "arx");
+ error(errmsg);
+ end
+
+ if (~isreal(z)) then
+ errmsg = msprintf(gettext("%s: input and output data matrix should be a real matrix"), "arx");
+ error(errmsg);
+ end
+
+ n = varargin(2)
+ if (size(n,"*")<2| size(n,"*")>3) then
+ errmsg = msprintf(gettext("%s: The order and delay matrix [na nb nk] should be of size [2 3]"), "arx");
+ error(errmsg);
+ end
+
+ if (size(find(n<0),"*") | size(find(((n-floor(n))<%eps)== %f))) then
+ errmsg = msprintf(gettext("%s: values of order and delay matrix [na nb nk] should be nonnegative integer number "), "arx");
+ error(errmsg);
+ end
+
+ na = n(1); nb = n(2); //nk = n(3); //nf = n(4);
+//
+ if (size(n,"*") == 2) then
+ nk = 1
+ else
+ nk = n(3);
+ end
+
+ // storing U(k) , y(k) and n data in UDATA,YDATA and NDATA respectively
+ YDATA = z(:,1);
+ UDATA = z(:,2);
+ NDATA = size(UDATA,"*");
+ function e = G(p,m)
+ e = YDATA - _objfunarx(UDATA,YDATA,p,na,nb,nk);
+ endfunction
+ tempSum = na+nb
+ p0 = linspace(0.1,0.9,tempSum)';
+ [var,errl] = lsqrsolve(p0,G,size(UDATA,"*"));
+ err = (norm(errl)^2);
+ opt_err = err;
+ resid = G(var,[]);
+
+ a = 1-poly([var(nb+1:nb+na)]',"q","coeff");
+ b = poly([repmat(0,nk,1);var(1:nb)]',"q","coeff");
+ a = (poly([1,-coeff(a)],'q','coeff'))
+ t = idpoly(coeff(a),coeff(b),1,1,1,Ts)
+
+ // estimating the other parameters
+ [temp1,temp2,temp3] = predict(z,t)
+ [temp11,temp22,temp33] = pe(z,t)
+
+ estData = calModelPara(temp1,temp1,n(1)+n(2))
+ //pause
+ t.Report.Fit.MSE = estData.MSE
+ t.Report.Fit.FPE = estData.FPE
+ t.Report.Fit.FitPer = estData.FitPer
+ t.Report.Fit.AIC = estData.AIC
+ t.Report.Fit.AICc = estData.AICc
+ t.Report.Fit.nAIC = estData.nAIC
+ t.Report.Fit.BIC = estData.BIC
+ t.TimeUnit = unit
+ sys = t
+
+endfunction
+
+function yhat = _objfunarx(UDATA,YDATA,x,na,nb,nk)
+ x=x(:)
+ q = poly(0,'q')
+ tempSum = nb+na
+ // making polynomials
+ b = poly([repmat(0,nk,1);x(1:nb)]',"q","coeff");
+ a = 1 - poly([x(nb+1:nb+na)]',"q","coeff")
+ aSize = coeff(a);bSize = coeff(b)
+ maxDelay = max([length(aSize) length(bSize)])
+ yhat = [YDATA(1:maxDelay)]
+ for k=maxDelay+1:size(UDATA,"*")
+ tempB = 0
+ for ii = 1:size(bSize,'*')
+ tempB = tempB + bSize(ii)*UDATA(k-ii+1)
+ end
+ tempA = 0
+ for ii = 1:size(aSize,"*")
+ tempA = tempA + aSize(ii)*YDATA(k-ii)
+ end
+ yhat = [yhat; [ tempA+tempB ]];
+ end
+endfunction
diff --git a/bj.sci b/bj.sci
new file mode 100644
index 0000000..f41ce65
--- /dev/null
+++ b/bj.sci
@@ -0,0 +1,122 @@
+// Estimates Discrete time BJ model
+// y(t) = [B(q)/F(q)]u(t) + [C(q)/D(q)]e(t)
+// Current version uses random initial guess
+// Need to get appropriate guess from OE and noise models
+
+// Authors: Ashutosh,Harpreet,Inderpreet
+// Updated(12-6-16)
+
+//function [theta_bj,opt_err,resid] = bj(varargin)
+function sys = bj(varargin)
+
+ [lhs , rhs] = argn();
+ if ( rhs < 2 ) then
+ errmsg = msprintf(gettext("%s: Unexpected number of input arguments : %d provided while should be 2"), "bj", rhs);
+ error(errmsg)
+ end
+
+ z = varargin(1)
+ if typeof(z) == 'iddata' then
+ Ts = z.Ts;unit = z.TimeUnit
+ z = [z.OutputData z.InputData]
+ elseif typeof(z) == 'constant' then
+ Ts = 1;unit = 'seconds'
+ end
+ if ((~size(z,2)==2) & (~size(z,1)==2)) then
+ errmsg = msprintf(gettext("%s: input and output data matrix should be of size (number of data)*2"), "bj");
+ error(errmsg);
+ end
+
+ if (~isreal(z)) then
+ errmsg = msprintf(gettext("%s: input and output data matrix should be a real matrix"), "bj");
+ error(errmsg);
+ end
+
+ n = varargin(2)
+ if (size(n,"*")<4| size(n,"*")>5) then
+ errmsg = msprintf(gettext("%s: The order and delay matrix [nb nc nd nf nk] should be of size [4 5]"), "bj");
+ error(errmsg);
+ end
+
+ if (size(find(n<0),"*") | size(find(((n-floor(n))<%eps)== %f))) then
+ errmsg = msprintf(gettext("%s: values of order and delay matrix [nb nc nd nf nk] should be nonnegative integer number "), "bj");
+ error(errmsg);
+ end
+
+ nb = n(1); nc = n(2); nd = n(3); nf = n(4);
+
+ if (size(n,"*") == 4) then
+ nk = 1
+ else
+ nk = n(5);
+ end
+
+ // storing U(k) , y(k) and n data in UDATA,YDATA and NDATA respectively
+ YDATA = z(:,1);
+ UDATA = z(:,2);
+ NDATA = size(UDATA,"*");
+ function e = G(p,m)
+ e = YDATA - _objfun(UDATA,p,nd,nc,nf,nb,nk);
+ endfunction
+ tempSum = nb+nc+nd+nf
+ p0 = linspace(0.5,0.9,tempSum)';
+ [var,errl] = lsqrsolve(p0,G,size(UDATA,"*"));
+
+ err = (norm(errl)^2);
+ opt_err = err;
+ resid = G(var,[]);
+ b = poly([repmat(0,nk,1);var(1:nb)]',"q","coeff");
+ c = poly([1; var(nb+1:nb+nc)]',"q","coeff");
+ d = poly([1; var(nb+nc+1:nb+nc+nd)]',"q","coeff");
+ f = poly([1; var(nb+nd+nc+1:nd+nc+nf+nb)]',"q","coeff");
+ t = idpoly(1,coeff(b),coeff(c),coeff(d),coeff(f),Ts)
+
+ // estimating the other parameters
+ [temp1,temp2,temp3] = predict(z,t)
+ [temp11,temp22,temp33] = pe(z,t)
+
+ estData = calModelPara(temp1,temp11,sum(n(1:4)))
+ //pause
+ t.Report.Fit.MSE = estData.MSE
+ t.Report.Fit.FPE = estData.FPE
+ t.Report.Fit.FitPer = estData.FitPer
+ t.Report.Fit.AIC = estData.AIC
+ t.Report.Fit.AICc = estData.AICc
+ t.Report.Fit.nAIC = estData.nAIC
+ t.Report.Fit.BIC = estData.BIC
+ t.TimeUnit = unit
+ sys = t
+endfunction
+
+function yhat = _objfun(UDATA,x,nd,nc,nf,nb,nk)
+ x=x(:)
+ q = poly(0,'q')
+ tempSum = nb+nc+nd+nf
+ // making polynomials
+ b = poly([repmat(0,nk,1);x(1:nb)]',"q","coeff");
+ c = poly([1; x(nb+1:nb+nc)]',"q","coeff");
+ d = poly([1; x(nb+nc+1:nb+nc+nd)]',"q","coeff");
+ f = poly([1; x(nb+nd+nc+1:nd+nc+nf+nb)]',"q","coeff");
+ bd = coeff(b*d); cf = coeff(c*f); fc_d = coeff(f*(c-d));
+ if size(bd,"*") == 1 then
+ bd = repmat(0,nb+nd+1,1)
+ end
+ maxDelay = max([length(bd) length(cf) length(fc_d)])
+ yhat = [YDATA(1:maxDelay)]
+ for k=maxDelay+1:size(UDATA,"*")
+ bdadd = 0
+ for i = 1:size(bd,"*")
+ bdadd = bdadd + bd(i)*UDATA(k-i+1)
+ end
+
+ fc_dadd = 0
+ for i = 1:size(fc_d,"*")
+ fc_dadd = fc_dadd + fc_d(i)*YDATA(k-i+1)
+ end
+ cfadd = 0
+ for i = 2:size(cf,"*")
+ cfadd = cfadd + cf(i)*yhat(k-i+1)
+ end
+ yhat = [yhat; [ bdadd + fc_dadd - cfadd ]];
+ end
+endfunction
diff --git a/calModelPara.sci b/calModelPara.sci
new file mode 100644
index 0000000..e8fbc02
--- /dev/null
+++ b/calModelPara.sci
@@ -0,0 +1,16 @@
+function varargout = calModelPara(varargin)
+ y = varargin(1)+varargin(2)
+ N = size(y,'r')
+ ek = varargin(2)
+ np = varargin(3)
+ //pause
+ mse = sum(ek.^2)/N;//disp(mse)
+ fpe = mse * (1 + np/N)/(1 - np/N);//disp(fpe)
+ nrmse = 1 - sqrt(sum(ek^2))/sqrt(sum((y - mean(y))^2));//disp(nrmse)
+ AIC = N * log(mse) + 2 * np + N * size(y,'c') * (log(2 * %pi) + 1);//disp(AIC)
+ AICc = AIC * 2 * np * (np + 1)/(N - np - 1);//disp(AICc)
+ nAIC = log(mse) + 2 * np/N;//disp(nAIC)
+ BIC = N * log(mse) + N * size(y,'c') * (log(2 * %pi) + 1) + np * log(N);//disp(BIC)
+ //pause
+ varargout(1) = struct('MSE',mse,'FPE',fpe,'FitPer',nrmse*100,'AIC',AIC,'AICc',AICc,'nAIC',nAIC,'BIC',BIC)
+endfunction
diff --git a/compare.sci b/compare.sci
new file mode 100644
index 0000000..4723130
--- /dev/null
+++ b/compare.sci
@@ -0,0 +1,145 @@
+function varargout = compare(varargin)
+ [lhs,rhs] = argn(0)
+//------------------------------------------------------------------------------
+// checking the number of inputs
+ if rhs < 2 || rhs > 3 then
+ error(msprintf(gettext("%s:Wrong number of input arguments.\n"),"compare"))
+ end
+//------------------------------------------------------------------------------
+ data = varargin(1)
+ model = varargin(2)
+ if rhs == 3 then
+ kStep = varargin(3)
+ elseif rhs == 2 then
+ kStep = %inf
+ end
+
+//------------------------------------------------------------------------------
+
+// k step analysis
+ if typeof(kStep) <> 'constant' || isnan(kStep) then
+ error(msprintf(gettext("%s:Prediction horizon(k) must be a non-negative integer number or inf.\n"),"compare"))
+ end
+// if given k step is infinity or []
+// if isinf(kStep) || ~size(kStep,'*') then
+// kStep = 1
+// end
+// checking the dimensions
+// if size(kStep,'*') <> 1 || (ceil(kStep)-kStep) then
+// error(msprintf(gettext("%s:Prediction horizon(k) must be a non-negative integer number or inf.\n"),"predict"))
+// end
+//------------------------------------------------------------------------------
+// checking the plant model
+ if typeof(model) ~= 'idpoly' then
+ error(msprintf(gettext("%s:Plant model must be ""idpoly"" type.\n"),"compare"))
+ end
+ modelSampleTime = model.Ts
+ modelTimeUnit = model.TimeUnit
+//------------------------------------------------------------------------------
+//checking the data type
+ if typeof(data) <> 'iddata' && typeof(data) <> 'constant' then
+ error(msprintf(gettext("%s:Sample data must be ""iddata"" type or ""n x 2"" matrix type.\n"),"compare"))
+ end
+// checking the plant data
+ if typeof(data) == 'iddata' then
+ if ~size(data.OutputData,'*') || ~size(data.InputData,'*') then
+ error(msprintf(gettext("%s:Number of sample data in input and output must be equal.\n"),"compare"))
+ end
+ plantSampleTime = data.Ts
+ plantTimeUnit = data.TimeUnit
+ data = [data.OutputData data.InputData]
+ //disp('iddata')
+ elseif typeof(data) == 'constant' then
+ if size(data,'c') ~= 2 then
+ error(msprintf(gettext("%s:Number of sample data in input and output must be equal.\n"),"compare"))
+ end
+ plantSampleTime = model.Ts
+ plantTimeUnit = model.TimeUnit
+ end
+//------------------------------------------------------------------------------
+// comparing the sampling time
+ if modelSampleTime-plantSampleTime <> 0 then
+ error(msprintf(gettext("%s:The sample time of the model and plant data must be equal.\n"),"compare"))
+ end
+// Comparing the time units
+ if ~strcmp(modelTimeUnit,plantTimeUnit) then
+ else
+ error(msprintf(gettext("%s:Time unit of the model and plant data must be equal.\n"),"compare"))
+ end
+//------------------------------------------------------------------------------
+// ckecking the k step size. if it greater than number of sample size then the
+// k step will become 1
+// if kStep >= size(data,'r') then
+// kStep = 1
+// end
+//------------------------------------------------------------------------------
+ //storing the plant data
+ // B(z) C(z)
+ // y(n) = ---------- u(n) + ---------- e(n)
+ // A(z)*F(z) A(z)*D(z)
+ aPoly = poly(model.a,'q','coeff');
+ bPoly = poly(model.b,'q','coeff');
+ cPoly = poly(model.c,'q','coeff');
+ dPoly = poly(model.d,'q','coeff');
+ fPoly = poly(model.f,'q','coeff');
+ Gq = bPoly/(aPoly*fPoly)
+ Hq = cPoly/(aPoly*dPoly)
+ //disp(kStep)
+ if kStep == %inf then
+ //disp('in inf')
+ outData = sim(data(:,2),model)
+ else
+ if kStep == 1 then
+ Wkq = Hq^-1
+ elseif kStep > 1 then
+ adCoeff = coeff(aPoly*dPoly);adCoeff = adCoeff(length(adCoeff):-1:1);
+ adPoly = poly(adCoeff,'q','coeff')
+ cCoeff = model.c;cCoeff = cCoeff(length(cCoeff):-1:1);
+ cPoly = poly(cCoeff,'q','coeff')
+ hBar = clean((ldiv(cPoly,adPoly,kStep))',0.00001)
+ hBarPoly = poly(hBar,'q','coeff')
+ Wkq = hBarPoly*Hq^-1//*bPoly/(dPoly*fPoly)
+ end
+ // pause
+ WkqGq = Wkq * Gq
+ tempWkqGq = coeff(WkqGq.den)
+ if tempWkqGq(1) <> 1 then
+ WkqGq.num = WkqGq.num/tempWkqGq(1)
+ WkqGq.den = WkqGq.den/tempWkqGq(1)
+ end
+ Wkq1 = 1-Wkq
+ tempWkq1 = coeff(Wkq1.den)
+ if tempWkq1(1) == 1 then
+ Wkq1.num = Wkq1.num/tempWkq1(1)
+ Wkq1.den = Wkq1.den/tempWkq1(1)
+ end
+ z = poly(0,'z');
+ WkqGqNum = horner(WkqGq.num,1/z);WkqGqDen = horner(WkqGq.den,1/z);
+ uPoly = WkqGqNum/WkqGqDen;
+ uPoly = syslin('d',uPoly);uData = flts(data(:,2)',uPoly);
+
+ Wkq1Num = horner(Wkq1.num,1/z);Wkq1Den = horner(Wkq1.den,1/z);
+ yPoly = Wkq1Num/Wkq1Den;
+ yPoly = syslin('d',yPoly);yData = flts(data(:,1)',yPoly);
+ outData = (uData+yData)'
+ end
+
+ tData = (1:size(data,'r'))'*plantSampleTime
+ fitData = fitValue(data(:,1),outData)
+ if lhs == 1 then
+ varargout(1) = []
+ plot(tData,data(:,1),'m')
+ plot(tData,outData,'b')
+ legend('Plant Data','Model Data : '+string(fitData)+'%')
+ xtitle('Comparison of Time Response','Time ('+ plantTimeUnit+')','Amplitude')
+ xgrid()
+ elseif lhs == 2 then
+ varargout(1) = fitData
+ varargout(2) = outData
+ elseif lhs == 3 then
+ varargout(1) = outData
+ varargout(2) = tData
+ varargout(3) = fitData
+ end
+
+endfunction
diff --git a/compareBJ.sci b/compareBJ.sci
new file mode 100755
index 0000000..06d2fb8
--- /dev/null
+++ b/compareBJ.sci
@@ -0,0 +1,56 @@
+function varargout = compareBJ(varargin)
+ //varargin(1) -> idpoly data about oe
+ //varargin(2) -> [y u] matrix of "nx2" dimension
+ //disp('compareBj')
+ bjData = varargin(1)
+ //disp(typeof(bjData))
+ plantData = varargin(2)
+ //disp(typeof(plantData))
+ yData = plantData(:,1)
+ uData = plantData(:,2)
+ uData = [0;uData]
+ //storing the data in A,B,C,D,F matrix
+ // B(z) C(z)
+ // y(n)=---------- u(n)+ ----------- e(n)
+ // A(z)F(z) A(z)D(z)
+ polyA = poly(bjData.a,'x','coeff')
+ polyB = poly(bjData.b,'x','coeff')
+ polyC = poly(bjData.c,'x','coeff')
+ polyD = poly(bjData.d,'x','coeff')
+ polyF = poly(bjData.f,'x','coeff')
+ adf = polyA*polyD*polyF
+ bd = polyB*polyD
+ cf = polyC*polyF
+ delay = max(size(coeff(adf),'*'),size(coeff(bd),'*'),size(coeff(cf),'*'))
+ yHat = [0]
+ bdCoeff = coeff(bd)
+ adfCoeff = coeff(adf)
+ adfCoeff = -adfCoeff(2:length(adfCoeff))
+ for ii = 1:length(uData)
+ uSum = 0;ySum = 0;
+ for jj = 1:length(bdCoeff)
+ if ii-jj <= 0 then
+ uSum = uSum + 0
+ else
+ uSum = uSum + uData(ii-jj+1)*bdCoeff(jj)
+ end
+ end
+ for jj = 1:length(adfCoeff)
+ if ii-jj <= 0 then
+ ySum = ySum + 0
+ else
+ ySum = ySum + yHat(ii-jj+1)*adfCoeff(jj)
+ end
+ end
+ yHat = [yHat; uSum+ySum]
+ end
+ tempStart = 1
+ if size(yHat,'r')- size(yData,'r') > 0 then
+ tempStart = size(yHat,'r')-size(yData,'r')+1
+ end
+ varargout(1) = yHat(tempStart:length(yHat))
+// plot(yHat(tempStart:length(yHat)),'m')
+// plot(yData)
+// xgrid()
+// pause
+endfunction
diff --git a/dataSlice.sci b/dataSlice.sci
new file mode 100644
index 0000000..6b3c160
--- /dev/null
+++ b/dataSlice.sci
@@ -0,0 +1,67 @@
+function sys = dataSlice(data,Start,End,Freq)
+ [lhs,rhs] = argn()
+ // storing the model data
+ modelData = data
+ // storing the statrting point
+ try
+ startData = Start
+ catch
+ startData = 1
+ end
+ //storing the end point
+ try
+ endData = End
+ catch
+ endData = LastIndex(data)
+ end
+ //Storing the frequency
+ try
+ freqData = Freq
+ catch
+ freqData = 1
+ end
+ // error message generate
+ if startData > endData then
+ error(msprintf(gettext("%s:Start index can not greater than End index.\n"),"dataSlice"))
+ end
+ if size(startData,'*') ~= 1 then
+ error(msprintf(gettext("%s:Start index must be non negative scalar integer number.\n"),"dataSlice"))
+ end
+ if size(endData,'*') ~= 1 then
+ error(msprintf(gettext("%s:End index must be non negative scalar integer number.\n"),"dataSlice"))
+ end
+ if ~freqData || size(freqData,'*') ~= 1 then
+ error(msprintf(gettext("%s:Frequency must be non negative scalar number.\n"),"dataSlice"))
+ end
+ //--------------------------------------------------------------------------
+ if typeof(modelData) == 'constant' then
+ Ts = 1
+ elseif typeof(modelData) == 'iddata' then
+ Ts = modelData.Ts
+ end
+ //--------------------------------------------------------------------------
+ if freqData> Ts || modulo(Ts,freqData) then
+ warning(msprintf(gettext("%s:inconsistent frequency.\n"),"dataSlice"))
+ freqData = Ts
+ end
+ if typeof(modelData) == 'constant' then
+ temp = modelData(startData:Ts/freqData:endData,:)
+ elseif typeof(modelData) == 'iddata' then
+ tempY = modelData.OutputData;tempU = modelData.InputData
+ tempY = tempY(startData:Ts/freqData:endData,:);tempU = tempU(startData:Ts/freqData:endData,:)
+ temp = iddata(tempY,tempU,Ts/freqData)
+ temp.TimeUnit = modelData.TimeUnit
+ end
+ sys = temp
+endfunction
+
+function varargout = LastIndex(modelData)
+ //finding the sample size
+ if typeof(modelData) == "constant" then
+ varargout(1) = length(modelData(:,1))
+
+ elseif typeof(modelData) == "iddata" then
+ temp = modelData.OutputData
+ varargout(1) = length(temp(:,1))
+ end
+endfunction
diff --git a/detrend.sci b/detrend.sci
new file mode 100644
index 0000000..993857e
--- /dev/null
+++ b/detrend.sci
@@ -0,0 +1,100 @@
+// Scilab ( http://www.scilab.org/ ) - This file is part of Scilab
+// Copyright (C) Bruno Pincon
+//
+// Copyright (C) 2012 - 2016 - Scilab Enterprises
+//
+// This file is hereby licensed under the terms of the GNU GPL v2.0,
+// pursuant to article 5.3.4 of the CeCILL v.2.1.
+// This file was originally licensed under the terms of the CeCILL v2.1,
+// and continues to be available under such terms.
+// For more information, see the COPYING file which you should have received
+// along with this program.
+//
+function [y] = detrend(x, flag, bp)
+ //
+ // this function removes the constant or linear or
+ // piecewise linear trend from a vector x. If x is
+ // matrix this function removes the trend of each
+ // column of x.
+ //
+ // flag = "constant" or "c" to removes the constant trend
+ // (simply the mean of the signal)
+ // flag = "linear" or "l" to remove the linear or piecewise
+ // linear trend.
+ //
+ // To define the piecewise linear trend the function needs the
+ // breakpoints and these must be given as the third argument bp.
+ //
+ // The "instants" of the signal x are in fact from 0 to m-1
+ // (m = length(x) if x is a vector and m = size(x,1) in case
+ // x is a matrix). So bp must be reals in [0 m-1].
+ //
+
+ rhs = argn(2)
+ if rhs < 1 | rhs > 3 then
+ error(msprintf(gettext("%s: Wrong number of input arguments: %d to %d expected.\n"),"detrend",1,3));
+ elseif rhs == 1
+ flag = "linear"; bp = []
+ elseif rhs == 2
+ bp = []
+ end
+
+ if type(x)~=1 then
+ error(msprintf(gettext("%s: Wrong type for input argument #%d: Array of floating point numbers expected.\n"),"detrend",1));
+ end
+ if type(flag)~=10 then
+ error(msprintf(gettext("%s: Wrong type for input argument #%d: Array of floating point numbers expected.\n"),"detrend",2));
+ end
+ if ~(type(bp)==1 & isreal(bp)) then
+ error(msprintf(gettext("%s: Wrong type for input argument #%d: Array of floating point numbers expected.\n"),"detrend",3));
+ end
+
+ [mx,nx] = size(x)
+ if mx==1 | nx==1 then
+ x_is_vector = %t; x = x(:); m = mx*nx; n = 1
+ elseif mx*nx == 0 then
+ y = []; return
+ else
+ x_is_vector = %f; m = mx; n = nx
+ end
+
+
+ if flag == "constant" | flag == "c" then
+ y = x - ones(m,1)*mean(x,"r")
+ elseif flag == "linear" | flag == "l"
+ bp = unique([0 ; bp(:) ; m-1])
+ // delete all the breakpoints outside [0,m-1]
+ while bp(1) < 0, bp(1)=[], end
+ while bp($) > m-1, bp($)=[], end
+ // breakpoints are 0-based so add one to
+ // compare them with signal vector indices (1-based)
+ bp = bp + 1;
+ nbp = length(bp);
+ // build the least square matrix with hat functions
+ // (as a basis for continuous piecewise linear functions)
+ A = zeros(m, nbp)
+ k1 = 1
+ delta_bp = diff(bp)
+ for j = 2:nbp-1
+ k2 = ceil(bp(j))-1
+ ind = (k1:k2)'
+ A(ind,j-1) = (bp(j) - ind)/delta_bp(j-1)
+ A(ind,j) = (ind - bp(j-1))/delta_bp(j-1)
+ k1 = k2+1
+ end
+ ind = (k1:m)'
+ A(ind,nbp-1) = (m - ind)/delta_bp(nbp-1)
+ A(ind,nbp) = (ind - bp(nbp-1))/delta_bp(nbp-1)
+ // solve the least square pb and retrieve the fitted
+ // piecewise linear func off the signal
+ y = x - A*(A\x)
+ else
+ error(msprintf(gettext("%s: Wrong value for input argument #%d: Must be in the set {%s}.\n") ,..
+ "detrend",2,"''constant'',''c'',''linear'',''l''"));
+ end
+
+ if x_is_vector then
+ y = matrix(y,mx,nx)
+ end
+
+endfunction
diff --git a/estpoly.sci b/estpoly.sci
new file mode 100644
index 0000000..f1ea48d
--- /dev/null
+++ b/estpoly.sci
@@ -0,0 +1,136 @@
+// Estimates Discrete time estpoly model
+// y(t) = [B(q)/F(q)]u(t) + [C(q)/D(q)]e(t)
+// Current version uses random initial guess
+// Need to get appropriate guess from OE and noise models
+
+// Authors: Ashutosh,Harpreet,Inderpreet
+// Updated(12-6-16)
+
+//function [theta_estpoly,opt_err,resid] = estpoly(varargin)
+function sys = estpoly(varargin)
+
+ [lhs , rhs] = argn();
+ if ( rhs < 2 ) then
+ errmsg = msprintf(gettext("%s: Unexpected number of input arguments : %d provided while should be 2"), "estpoly", rhs);
+ error(errmsg)
+ end
+
+ z = varargin(1)
+ if typeof(z) == 'iddata' then
+ Ts = z.Ts;unit = z.TimeUnit
+ z = [z.OutputData z.InputData]
+ elseif typeof(z) == 'constant' then
+ Ts = 1;unit = 'seconds'
+ end
+ if ((~size(z,2)==2) & (~size(z,1)==2)) then
+ errmsg = msprintf(gettext("%s: input and output data matrix should be of size (number of data)*2"), "estpoly");
+ error(errmsg);
+ end
+
+ if (~isreal(z)) then
+ errmsg = msprintf(gettext("%s: input and output data matrix should be a real matrix"), "estpoly");
+ error(errmsg);
+ end
+
+ n = varargin(2)
+ if (size(n,"*")<5| size(n,"*")>6) then
+ errmsg = msprintf(gettext("%s: The order and delay matrix [na nb nc nd nf nk] should be of size [5 6]"), "estpoly");
+ error(errmsg);
+ end
+
+ if (size(find(n<0),"*") | size(find(((n-floor(n))<%eps)== %f))) then
+ errmsg = msprintf(gettext("%s: values of order and delay matrix [nb nc nd nf nk] should be nonnegative integer number "), "estpoly");
+ error(errmsg);
+ end
+
+ na = n(1); nb = n(2); nc = n(3); nd = n(4);nf = n(5);
+
+ if (size(n,"*") == 5) then
+ nk = 1
+ else
+ nk = n(6);
+ end
+
+ // storing U(k) , y(k) and n data in UDATA,YDATA and NDATA respectively
+ YDATA = z(:,1);
+ UDATA = z(:,2);
+ NDATA = size(UDATA,"*");
+ function e = G(p,m)
+ e = YDATA - _oestpolyfun(UDATA,p,na,nb,nc,nd,nf,nk);//_oestpolyfun(UDATA,p,nd,nc,nf,nb,nk);
+ endfunction
+ tempSum = na+nb+nc+nd+nf
+ p0 = linspace(0.0001,0.001,tempSum)';
+ [var,errl] = lsqrsolve(p0,G,size(UDATA,"*"));
+ //disp(errl)
+ err = (norm(errl)^2);
+ //disp(err)
+ opt_err = err;
+ resid = G(var,[]);
+ x = var
+// b = poly([repmat(0,nk,1);var(1:nb)]',"q","coeff");
+// c = poly([1; var(nb+1:nb+nc)]',"q","coeff");
+// d = poly([1; var(nb+nc+1:nb+nc+nd)]',"q","coeff");
+// f = poly([1; var(nb+nd+nc+1:nd+nc+nf+nb)]',"q","coeff");
+ a = poly([1; x(1:na)]',"q","coeff");
+ b = poly([repmat(0,nk,1);x(na+1:na+nb)]',"q","coeff");
+ c = poly([1; x(na+nb+1:na+nb+nc)]',"q","coeff");
+ d = poly([1; x(na+nb+nc+1:na+nb+nc+nd)]',"q","coeff");
+ f = poly([1; x(na+nb+nd+nc+1:na+nd+nc+nf+nb)]',"q","coeff");
+ t = idpoly(coeff(a),coeff(b),coeff(c),coeff(d),coeff(f),Ts)
+
+ //t = sys;//idpoly(1,coeff(b),coeff(c),coeff(d),coeff(f),Ts)
+
+ // estimating the other parameters
+ [temp1,temp2,temp3] = predict(z,t)
+ [temp11,temp22,temp33] = pe(z,t)
+ //pause
+ estData = calModelPara(temp1,temp11,na+nb+nc+nd+nf)
+ //pause
+ t.Report.Fit.MSE = estData.MSE
+ t.Report.Fit.FPE = estData.FPE
+ t.Report.Fit.FitPer = estData.FitPer
+ t.Report.Fit.AIC = estData.AIC
+ t.Report.Fit.AICc = estData.AICc
+ t.Report.Fit.nAIC = estData.nAIC
+ t.Report.Fit.BIC = estData.BIC
+ t.TimeUnit = unit
+ sys = t
+
+ //sys.TimeUnit = unit
+endfunction
+
+function yhat = _oestpolyfun(UDATA,x,na,nb,nc,nd,nf,nk)//(UDATA,x,nd,nc,nf,nb,nk)
+ x=x(:)
+ q = poly(0,'q')
+ tempSum = na+nb+nc+nd+nf
+ // making polynomials
+ a = poly([1; x(1:na)]',"q","coeff");
+ b = poly([repmat(0,nk,1);x(na+1:na+nb)]',"q","coeff");
+ c = poly([1; x(na+nb+1:na+nb+nc)]',"q","coeff");
+ d = poly([1; x(na+nb+nc+1:na+nb+nc+nd)]',"q","coeff");
+ f = poly([1; x(na+nb+nd+nc+1:na+nd+nc+nf+nb)]',"q","coeff");
+ bd = coeff(b*d); cf = coeff(c*f); fc_d = coeff(f*(c-a*d));
+ if size(bd,"*") == 1 then
+ bd = repmat(0,nb+nd+1,1)
+ end
+ //pause
+ maxDelay = max([length(bd) length(cf) length(fc_d)])
+ yhat = [YDATA(1:maxDelay)]
+ for k=maxDelay+1:size(UDATA,"*")
+ bdadd = 0
+ for i = 1:size(bd,"*")
+ bdadd = bdadd + bd(i)*UDATA(k-i+1)
+ end
+
+ fc_dadd = 0
+ for i = 1:size(fc_d,"*")
+ fc_dadd = fc_dadd + fc_d(i)*YDATA(k-i+1)
+ end
+ cfadd = 0
+ for i = 2:size(cf,"*")
+ cfadd = cfadd + cf(i)*yhat(k-i+1)
+ end
+ //pause
+ yhat = [yhat; [ bdadd + fc_dadd - cfadd ]];
+ end
+endfunction
diff --git a/etfe.sci b/etfe.sci
new file mode 100644
index 0000000..b09b7c4
--- /dev/null
+++ b/etfe.sci
@@ -0,0 +1,29 @@
+function varargout = etfe(varargin)
+ [lhs,rhs] = argn()
+ data = varargin(1)
+ if rhs == 1 then
+ n = 128
+ elseif rhs == 2 then
+ n = varargin(2)
+ end
+ y = data.OutputData;
+ u = data.InputData
+ N = size(y,'r')
+ v = linspace(1,N,n)
+ y1 = y(v);u1 = u(v)
+ //y1 = y((1:(N-1)/(n-1):N));u1 = u((1:(N-1)/(n-1):N))
+ //y1($) = y(N);u1($) = u(N)
+ data12 = [y1,u1]
+ z = [fft(y1),fft(u1)]
+ z = z/size(z,'r')
+ magData1 = abs(z(:,1));magData2 = abs(z(:,2))
+ argData1 = phasemag(z(:,1),'m');argData2 = phasemag(z(:,2),'m')
+ magData = magData1./magData2;argData = argData1-argData2
+ argData = [cosd(argData) sind(argData)]
+ data = [magData.*argData(:,1) magData.*argData(:,2)]
+ output = data(:,1)+%i*data(:,2)
+ resp = output(1:ceil(length(output)/2))
+ frq = (1: ceil(n/2)) * %pi/floor(n/2)
+ output = frd(frq,resp,1)
+ varargout(1)= output
+endfunction
diff --git a/etfeeTest.sce b/etfeeTest.sce
new file mode 100644
index 0000000..46f250c
--- /dev/null
+++ b/etfeeTest.sce
@@ -0,0 +1,13 @@
+loadmatfile('rData.mat')
+z = [fft(data1(:,1)),fft(data1(:,2))]
+z = z/size(z,'r')
+magData1 = abs(z(:,1));magData2 = abs(z(:,2))
+argData1 = phasemag(z(:,1),'m');argData2 = phasemag(z(:,2),'m')
+magData = magData1./magData2;argData = argData1-argData2
+argData = [cosd(argData) sind(argData)]
+data = [magData.*argData(:,1) magData.*argData(:,2)]
+output = data(:,1)+%i*data(:,2)
+resp = output(1:ceil(length(output)/2))
+frq = (1: ceil(256/2)) * %pi/floor(256/2)
+output = frd(frq,resp,1)
+disp(output)
diff --git a/fitch.sci b/fitch.sci
new file mode 100644
index 0000000..b72b472
--- /dev/null
+++ b/fitch.sci
@@ -0,0 +1,42 @@
+function varargout = fitch(varargin)
+ [lhs , rhs] = argn()
+ if ( rhs <> 1 ) then
+ errmsg = msprintf(gettext("%s: Wrong number of input arguments"), "fitch");
+ error(errmsg)
+ elseif typeof(varargin(1)) <> "idpoly" then
+ error(msprintf(gettext("%s:Input model must be ""idpoly"" type.\n"),"fitch"))
+ end
+ model = varargin(1)
+ MSE = model.Report.Fit.MSE
+ FPE = model.Report.Fit.FPE
+ FitPer = model.Report.Fit.FitPer
+ AIC = model.Report.Fit.AIC
+ AICc = model.Report.Fit.AICc
+ nAIC = model.Report.Fit.nAIC
+ BIC = model.Report.Fit.BIC
+ t = tlist(['fitch','MSE','FPE','FitPer','AIC','AICc','nAIC','BIC'],MSE,FPE,FitPer,AIC,AICc,nAIC,BIC)
+ varargout(1) = t
+endfunction
+
+function %fitch_p(mytlist)
+ f = fieldnames(mytlist)
+ maxLength = []
+ for ii = 1:size(f,'*')
+ maxLength = [maxLength length(f(ii))]
+ end
+ maxLength = max(maxLength)
+ for ii = 1:size(f,'*')
+ blanckSpace = ' '
+ for jj = 1:maxLength-length(f(ii))
+ blanckSpace = blanckSpace + ' '
+ end
+ mprintf('\t%s%s : ',blanckSpace,f(ii))
+ objectData = mytlist(f(ii))
+ if ceil(objectData)-objectData then
+ mprintf("%.4f",objectData)
+ else
+ mprintf("%d",objectData)
+ end
+ mprintf("\n")
+ end
+endfunction
diff --git a/frd.sci b/frd.sci
new file mode 100644
index 0000000..c3838b4
--- /dev/null
+++ b/frd.sci
@@ -0,0 +1,84 @@
+function varargout = frd(varargin)
+ [lhs,rhs] = argn(0)
+ if rhs < 2 || rhs > 4 then
+ errmsg = msprintf(gettext("%s: Wrong numbers of input arguments."), "frd");
+ error(errmsg)
+ end
+ frequency = varargin(2)
+ freqUnit = 'rad/TimeUnit'
+ if ~iscolumn(frequency) then
+ errmsg = msprintf(gettext("%s: frequency must be a finite row vector."), "frd");
+ error(errmsg)
+ end
+ respData = varargin(1)
+// pause
+ if size(frequency,'r') <> size(respData,'r') then
+ errmsg = msprintf(gettext("%s: input output matrix dimension must be equal."), "frd");
+ error(errmsg)
+ end
+ if rhs == 2 then
+ Ts = 0
+ elseif rhs >2 then
+ Ts = varargin(3)
+ end
+ if Ts < 0 || size(Ts,'*') <> 1 || typeof(Ts) <> 'constant' then
+ errmsg = msprintf(gettext("%s: Sampling time must be a scalar non negative real number."), "frd");
+ error(errmsg)
+ end
+ // saving the spectrum value
+ if rhs == 4 then
+ spect = varargin(4)
+ else
+ spect = []
+ end
+ /// matching its dimensions
+ if ~size(spect) then
+ elseif size(frequency,'r') <> size(spect,'r') then
+ errmsg = msprintf(gettext("%s: Numbers of power spectra must be equal to the numbers of frequency."), "frd");
+ error(errmsg)
+ end
+ TUnit = 'seconds'
+ t = tlist(['frd','Frequency','FrequencyUnit','ResponseData','Ts','TimeUnit','Spect'],frequency,freqUnit,respData,Ts,TUnit,spect)
+ varargout(1) = t
+endfunction
+
+//overloading
+function %frd_p(varargin)
+ myTlist = varargin(1)
+ f = fieldnames(myTlist)
+ freqData = myTlist.Frequency
+ tempRespData= myTlist.ResponseData
+ for jj = 1:size(tempRespData,'c')
+ respData = tempRespData(:,jj)
+ mprintf("\t -------------------------")
+ mprintf("\n")
+ mprintf("\t Frequency \t Response")
+ mprintf("\n")
+ mprintf("\t -------------------------")
+ mprintf("\n")
+ for ii = 1:length(myTlist.Frequency)
+ temp = ''
+ if real(respData(ii))>=0 then
+ temp = temp + ' '
+ end
+ temp = temp + string(real(respData(ii)))
+ // temp = string(real(respData(ii)))
+ if imag(respData(ii)) > 0 then
+ temp = temp +"+"
+ end
+ if ~imag(respData(ii)) then
+ else
+ temp = temp + string(imag(respData(ii))) +"i"
+ end
+ // temp = temp + string(imag(respData(ii))) + " i"
+ mprintf("\n\t %f \t %s",freqData(ii),temp)//real(respData(ii)),imag(respData(ii)))
+ end
+ mprintf("\n\n")
+ end
+ if ~myTlist.Ts then
+ mprintf("\n Continuous Domain frequency response.")
+ else
+ mprintf("\n Sampling Time = "+string(myTlist.Ts)+" "+myTlist.TimeUnit)
+ mprintf("\n Discrete Domain frequency response.")
+ end
+endfunction
diff --git a/frdplot.sci b/frdplot.sci
new file mode 100644
index 0000000..c923429
--- /dev/null
+++ b/frdplot.sci
@@ -0,0 +1,11 @@
+function frdplot(varargin)
+ [lhs rhs] = argn(0)
+ if rhs <> 1 then
+ error(msprintf(gettext("%s: Wrong number of input arguments."),"frdplot"))
+ end
+ frdData = varargin(1)
+ if typeof(frdData) <> 'frd' then
+ error(msprintf(gettext("%s:Wrong type for input argument %d: ""frd"" expected."),"frdplot",1))
+ end
+ bode((frdData.Frequency)',(frdData.ResponseData)')
+endfunction
diff --git a/iddata.sci b/iddata.sci
new file mode 100755
index 0000000..e748029
--- /dev/null
+++ b/iddata.sci
@@ -0,0 +1,68 @@
+function varargout = iddata(varargin)
+ [lhs,rhs] = argn(0)
+ if rhs == 2 | rhs == 1 then
+ Ts = 1
+ elseif rhs == 3 then
+ Ts = varargin(3)
+ else
+ error(msprintf(gettext("%s:Incorrect number of input arguments.\n"),"iddata"))
+ end
+ if Ts <= 0 || typeof(Ts) <> 'constant' || ~size(Ts,'*') then
+ error(msprintf(gettext("%s:Inconsist sampling time ""Ts"", ""non negative real number"" expected.\n"),"iddata"))
+ end
+ if rhs == 1 then
+ OutputData = varargin(1)
+ InputData = []
+ elseif rhs == 2 | rhs == 3 then
+ OutputData = varargin(1)
+ InputData = varargin(2)
+ if isrow(InputData) then
+ InputData = InputData'
+ end
+ end
+ if size(OutputData,'*') & size(InputData,'*') then
+ if size(OutputData,'r') ~= size(InputData,'r') then
+ error(msprintf(gettext("%s:The numbers of the plant out datas must be equal to the numbers of the plant input datas.\n"),"iddata"))
+ end
+ end
+ t = tlist(['iddata','OutputData','InputData','Ts','TimeUnit'],OutputData,InputData,Ts,'seconds')
+ varargout(1) = t
+endfunction
+
+function %iddata_p(mytlist)
+ f = fieldnames(mytlist)
+ if ~size(mytlist(f(1)),'*') & ~size(mytlist(f(2)),'*') then
+ mprintf(' Empty sample data.\n')
+ else
+ outputSize = size(mytlist(f(1)))
+ inputSize = size(mytlist(f(2)))
+ if prod(outputSize) then
+ sampleSize = max(outputSize)
+ elseif prod(inputSize) then
+ sampleSize = max(inputSize)
+ end
+ mprintf(' Time domain sample data having %d samples.',sampleSize)
+ if round(mytlist(f(3)))-mytlist(f(3)) == 0 then
+ mprintf('\n Sampling Time = %d',mytlist(f(3)))
+ else
+ mprintf('\n Sampling Time = %f',mytlist(f(3)))
+ end
+ mprintf(' %s',mytlist(f(4)))
+ mprintf('\n')
+ if prod(outputSize) then
+ mprintf('\n Output channel \n')
+ for ii = 1:min(outputSize)
+ yString = 'y' + string(ii)
+ mprintf(' %s\n',yString)
+ end
+ end
+ if prod(inputSize) then
+ mprintf('\n Input channel \n')
+ for ii = 1:min(inputSize)
+ uString = 'u' + string(ii)
+ mprintf(' %s\n',uString)
+ end
+ end
+ end
+ mprintf('\n')
+endfunction
diff --git a/iddataplot.sci b/iddataplot.sci
new file mode 100644
index 0000000..4b1818d
--- /dev/null
+++ b/iddataplot.sci
@@ -0,0 +1,69 @@
+function iddataplot(varargin)
+ [lhs rhs] = argn(0)
+ if rhs <> 1 then
+ error(msprintf(gettext("%s: Wrong number of input arguments."),"iddataplot"))
+ end
+ iddataData = varargin(1)
+ if typeof(iddataData) <> 'iddata' then
+ error(msprintf(gettext("%s:Wrong type for input argument %d: ""iddata"" expected."),"iddataplot",1))
+ end
+// figure()
+// xtitle('Plant i/o Data','Time('+iddataData.TimeUnit+')','Amplitude')
+ uData = iddataData.InputData
+ yData = iddataData.OutputData
+ timeLength = max(size(uData,'r'),size(yData,'r'))
+ timeData = ((1:1:timeLength)*iddataData.Ts)'
+ //timeData = timeData(1:length(iddataData)-1)
+ if size(uData,'*') && size(yData,'*') then
+ firstIndex = 2
+ else
+ firstIndex = 1
+ end
+ // ploting y data
+ h = gcf()
+ if h.figure_name == "Plant Input Output Data"; then
+ clf()
+ end
+ if size(yData,'*') then
+ secondIndex = size(yData,'c')
+ if secondIndex == 1 then
+ outputString = 'y'
+ else
+ outputString = []
+ for ii = 1:secondIndex
+ outputString = [outputString 'y'+string(ii)]
+ end
+ end
+ //disp(outputString)
+ for ii = 1:secondIndex
+ subplot(firstIndex,secondIndex,ii);plot(timeData,yData(:,ii))
+ xtitle(outputString(ii))
+ end
+ end
+ if size(uData,'*') then
+ secondIndex = size(uData,'c')
+ // disp(secondIndex)
+ if secondIndex == 1 then
+ outputString = 'u'
+ else
+ outputString = []
+ for ii = 1:secondIndex
+ outputString = [outputString 'u'+string(ii)]
+ end
+ end
+ //disp(outputString)
+ for ii = 1:secondIndex
+ if size(yData,'*') then
+ temp = 1
+ else
+ temp = 0
+ end
+ subplot(firstIndex,secondIndex,ii+secondIndex*temp);plot2d(timeData,uData(:,ii),2)
+ //pause
+ xtitle(outputString(ii))
+ end
+ end
+ h = gcf()
+ //disp(h)
+ h.figure_name= "Plant Input Output Data";
+endfunction
diff --git a/identTime.sci b/identTime.sci
new file mode 100644
index 0000000..c4c47c8
--- /dev/null
+++ b/identTime.sci
@@ -0,0 +1,32 @@
+function varargout = identTime(varargin)
+ [lhs , rhs] = argn()
+ if ( rhs <> 1 ) then
+ errmsg = msprintf(gettext("%s: Wrong number of input arguments"), "identTime");
+ error(errmsg)
+ elseif typeof(varargin(1)) <> "iddata" then
+ error(msprintf(gettext("%s:Input model must be ""iddata"" type.\n"),"identTime"))
+ end
+ plantData = varargin(1)
+ //disp('yolo')
+ inputData = plantData.InputData;inputData = size(inputData,'r')
+ outputData = plantData.OutputData;outputData = size(outputData,'r')
+ sampleNumb = max(inputData,outputData)
+ timeData = (0:sampleNumb-1)*plantData.Ts
+ t = tlist(['identTime','samples','start','end','Frequency','TimeSeries'],sampleNumb,0,timeData($),1/plantData.Ts,timeData)
+ varargout(1) = t
+endfunction
+
+function %identTime_p(mytlist)
+ f = fieldnames(mytlist)
+ mprintf("\t samples : %d\n",mytlist.samples)
+ mprintf("\t start : %d\n",mytlist.start)
+ mprintf("\t end : %d\n",mytlist.end)
+ if ceil(mytlist.Frequency)-mytlist.Frequency then
+ mprintf("\t Frequency : %.4f\n",mytlist.Frequency)
+ else
+ mprintf("\t Frequency : %d\n",mytlist.Frequency)
+ end
+ timeData = mytlist.TimeSeries
+ mprintf("\t TimeSeries : %.2f, %.2f, . ,%.2f",timeData(1),timeData(2),timeData($))
+
+endfunction
diff --git a/idpoly.sci b/idpoly.sci
new file mode 100755
index 0000000..0be3754
--- /dev/null
+++ b/idpoly.sci
@@ -0,0 +1,142 @@
+function sys = idpoly(varargin)
+ [lhs,rhs] = argn(0)
+ tempCell = cell(6,1)
+ tempTs = 0
+ for ii = 1:rhs
+ if typeof(varargin(ii)) == 'string' & varargin(ii) == 'Ts' then
+ tempCell{6,1} = varargin(ii+1)
+ break
+ elseif ~size(varargin(ii),'*') then
+ tempCell{ii,1} = 1
+ elseif size(varargin(ii),'*') then
+ tempCell{ii,1} = varargin(ii)
+ end
+ end
+ //disp(tempCell)
+ for ii = 1:6
+ if ~size(cell2mat(tempCell(ii,1)),'*') & ii ~= 6 then
+ tempCell{ii,1} = 1
+ elseif ~size(cell2mat(tempCell(ii,1)),'*') & ii == 6 then
+ tempCell{ii,1} = -1
+ end
+ end
+ //storing the data in A,B,C,D,F matrix
+ // B(z) D(z)
+ // y(n)=---------- u(n)+ ----------- e(n)
+ // A(z)F(z) A(z)D(z)
+ A = cell2mat(tempCell(1,1)); B = cell2mat(tempCell(2,1));
+ C = cell2mat(tempCell(3,1)); D = cell2mat(tempCell(4,1));
+ F = cell2mat(tempCell(5,1)); Ts = cell2mat(tempCell(6,1));
+ if A(1,1) ~= 1 then
+ error(msprintf(gettext("%s: The first coefficient of ""A(z)"" polynomial must be 1.\n"),"idpoly"))
+ elseif C(1,1) ~= 1 then
+ error(msprintf(gettext("%s: The first coefficient of ""C(z)"" polynomial must be 1.\n"),"idpoly"))
+ elseif D(1,1) ~= 1 then
+ error(msprintf(gettext("%s: The first coefficient of ""D(z)"" polynomial must be 1.\n"),"idpoly"))
+ elseif F(1,1) ~= 1 then
+ error(msprintf(gettext("%s: The first coefficient of ""F(z)"" polynomial must be 1.\n"),"idpoly"))
+ end
+ report = struct('MSE',0,'FPE',0,'FitPer',0,'AIC',0,'AICc',0,'nAIC',0,'BIC',0)
+ errors = [0 0 0]
+ report = struct('Fit',report,'Uncertainty',errors)
+ t = tlist(['idpoly','a','b','c','d','f','Variable','TimeUnit','Ts','Report'],A,B,C,D,F,'z^-1','seconds',Ts,report)
+ //t = tlist(['idpoly','a','b','c','d','f','Variable','TimeUnit','Ts'],A,B,C,D,F,'z^-1','seconds',Ts)
+
+ sys = t
+endfunction
+
+
+function %idpoly_p(mytlist)
+ f = fieldnames(mytlist)
+ //A polynomial
+ if mytlist(f(1)) == 1 && size(mytlist(f(1)),'*') == 1 then
+ else
+ mprintf('\n A(z) =')
+ temp = poly2str(mytlist(f(1)))
+ mprintf('%s\n',temp)
+ end
+
+ //B polynomial
+ if mytlist(f(2)) == 1 then
+ else
+ mprintf('\n B(z) =')
+ temp = poly2str(mytlist(f(2)))
+ mprintf('%s\n',temp)
+ end
+
+ //C polynomial
+ if mytlist(f(3)) == 1 && size(mytlist(f(3)),'*') == 1 then
+ else
+ mprintf('\n C(z) =')
+ temp = poly2str(mytlist(f(3)))
+ mprintf('%s\n',temp)
+ end
+ //D polynomial
+ if mytlist(f(4)) == 1 && size(mytlist(f(4)),'*') == 1 then
+ elseif size(mytlist(f(4)),'*') > 1 then
+ mprintf('\n D(z) =')
+ temp = poly2str(mytlist(f(4)))
+ mprintf('%s\n',temp)
+ end
+
+ //F polynomial
+ if mytlist(f(5)) == 1 && size(mytlist(f(5)),'*') == 1 then
+ else
+ mprintf('\n F(z) =')
+ temp = poly2str(mytlist(f(5)))
+ mprintf('%s\n',temp)
+ end
+
+ mprintf('\n Sampling Time = ')
+
+ if mytlist.Ts == -1 then
+ mprintf('undefined')
+ else
+ if (ceil(mytlist.Ts)-mytlist.Ts) == 0 then
+ mprintf('%d %s',mytlist.Ts,mytlist.TimeUnit)
+ else
+ mprintf('%0.4f %s',mytlist.Ts,mytlist.TimeUnit)
+ end
+ end
+ //disp(mytlist.Ts)
+ mprintf('\n')
+ if mytlist.Report.Fit.MSE then
+ temp = ['MSE','FPE','FitPer','AIC','AICc','nAIC','BIC']
+ spaces = ' '
+ for ii = 1:size(temp,'c')
+ digiLength = length(string(round(mytlist.Report.Fit(temp(ii)))))
+ digiLength = digiLength + 5-length(temp(ii))
+ blank = ''
+ for jj = 1:digiLength+1
+ blank = blank + " "
+ end
+ spaces = spaces+blank+temp(ii)+' '
+ end
+ mprintf('\n')
+ mprintf(spaces)
+ mprintf("\n %.4f %.4f %.4f %.4f %.4f %.4f %.4f",mytlist.Report.Fit.MSE,mytlist.Report.Fit.FPE,mytlist.Report.Fit.FitPer,mytlist.Report.Fit.AIC,mytlist.Report.Fit.AICc,mytlist.Report.Fit.nAIC,mytlist.Report.Fit.BIC)
+ end
+
+
+endfunction
+
+
+function strout = poly2str(h)
+ temp = poly(h,'x','coeff')
+ temp = pol2str(temp)
+ temp = strsubst(temp,'-',' - ')
+ temp = strsubst(temp,'x^',' z^-')
+ temp = strsubst(temp,'x',' z^-1')
+ temp = strsubst(temp,'*','')
+ temp = strsubst(temp,'+',' + ')
+ [ind which]=strindex(temp,'-')
+ //disp(ind)
+// disp(which)
+ if ind(1,1) ~= 2 then
+ temp = ' ' + temp
+ elseif ind(1,1) == 2 then
+ temp = part(temp,4:length(temp))
+ temp = ' -' + temp
+ end
+ strout = temp
+endfunction
diff --git a/impulseest.sci b/impulseest.sci
new file mode 100644
index 0000000..8817bcc
--- /dev/null
+++ b/impulseest.sci
@@ -0,0 +1,78 @@
+function varargout = impulseest(varargin)
+ [lhs,rhs] = argn(0)
+ //checking the number of inputs
+ if rhs > 2 then
+ error(msprintf(gettext("%s: Unexpected number of input arguments "),"impulseest"))
+ end
+ modelData = varargin(1)
+ if typeof(modelData) <> "idpoly" then
+ error(msprintf(gettext("%s: Plant model must be ""idpoly"" type. "),"impulseest"))
+ end
+ //adding noise
+ if rhs == 2 then
+ noiseFlag = varargin(2)
+ if typeof(noiseFlag) <> 'boolean' then
+ error(msprintf(gettext("%s: Last input data must be ""boolean"".type "),"impulseest"))
+ end
+ else
+ noiseFlag = %F
+ end
+ z = poly(0,'z')
+ aPoly = poly(modelData.a(length(modelData.a):-1:1),'z','coeff')
+ bPoly = poly(modelData.b,'z','coeff')
+ fPoly = poly(modelData.f(length(modelData.f):-1:1),'z','coeff')
+ afPoly = aPoly*fPoly
+
+ bCoeff = modelData.b
+ extra = 1
+ if ~bCoeff(1,1) then
+ afCoeff = coeff(afPoly)
+ bLength = length(bCoeff);afLength = length(afCoeff)
+ if bLength == afLength then
+ extra = 1
+ else
+ extra = z^-(bLength-afLength)
+ end
+ end
+
+ bPoly = poly(modelData.b(length(modelData.b):-1:1),'z','coeff')
+ sys = syslin('d',bPoly,afPoly)*extra
+ if size(find(gsort(abs(roots(sys.den)))>1),'*') then
+ tempRoot = clean(roots(sys.den))
+ [sorted index] = gsort(abs(real(tempRoot)))
+ tempRoot = tempRoot(index)
+ n = 26/log10(abs(real(tempRoot(1))))
+
+ else
+ finalValue = horner(sys,1)*1
+ n = 10
+ tempValue = sum(ldiv(sys.num,sys.den,n))
+ if finalValue >= 0 then
+ while finalValue*0.99 >= tempValue || n >=1000
+ n = n+5
+ tempValue = sum(ldiv(sys.num,sys.den,n))
+ end
+ elseif finalValue < 0 then
+ while finalValue*0.99 <tempValue || n >=1000
+ n = n+5
+ tempValue = sum(ldiv(sys.num,sys.den,n))
+ end
+ if n > 100 then
+ n = 100
+ end
+ end
+ end
+ uData = [1 zeros(1,n)]
+ yData = flts(uData,sys)
+ timeData = (0:(n))*modelData.Ts
+// pause
+ if noiseFlag then
+ varargout(1) = yData'
+ else
+ stem(timeData,yData)
+ h = gcf()
+
+ h.figure_name= "Plant Impulse Response"
+ varargout(1) = []
+ end
+endfunction
diff --git a/iv.sci b/iv.sci
new file mode 100644
index 0000000..a819a4e
--- /dev/null
+++ b/iv.sci
@@ -0,0 +1,98 @@
+function varargout = iv(varargin)
+
+ [lhs , rhs] = argn(0);
+ if ( rhs < 2 || rhs > 3) then
+ errmsg = msprintf(gettext("%s: Unexpected number of input arguments : %d provided while should be 2 or 3"), "iv", rhs);
+ error(errmsg)
+ end
+ plantData = varargin(1)
+ if typeof(plantData) == 'iddata' then
+ Ts = plantData.Ts;unit = plantData.TimeUnit
+ plantData = [plantData.OutputData plantData.InputData]
+ elseif typeof(plantData) == 'constant' then
+ Ts = 1;unit = 'seconds'
+ end
+ if ((~size(plantData,2)==2) & (~size(plantData,1)==2)) then
+ errmsg = msprintf(gettext("%s: input and output data matrix should be of size (number of data)*2"), "iv");
+ error(errmsg);
+ end
+
+ if (~isreal(plantData)) then
+ errmsg = msprintf(gettext("%s: input and output data matrix should be a real matrix"), "arx");
+ error(errmsg);
+ end
+
+ n = varargin(2)
+ if (size(n,"*")<2| size(n,"*")>3) then
+ errmsg = msprintf(gettext("%s: The order and delay matrix [na nb nk] should be of size [2 or 3]"), "iv");
+ error(errmsg);
+ end
+
+ if (size(find(n<0),"*") | size(find(((n-floor(n))<%eps)== %f))) then
+ errmsg = msprintf(gettext("%s: values of order and delay matrix [na nb nk] should be nonnegative integer number "), "iv");
+ error(errmsg);
+ end
+ na = n(1);nb = n(2)
+ if size(n,'*') == 2 then
+ nk = 1
+ elseif size(n,'*') == 3 then
+ nk = n(3)
+ end
+ yData = plantData(:,1)
+ uData = plantData(:,2)
+ noOfSample = size(plantData,'r')
+ nb1 = nb + nk - 1
+ n = max(na, nb1)
+ if rhs == 3 then
+ if typeof(varargin(3)) <> 'constant' then
+ errmsg = msprintf(gettext("%s: Incompatible last input argument. "), "iv");
+ error(errmsg)
+ elseif size(varargin(3),'r') <> size(uData,'r') then
+ errmsg = msprintf(gettext("%s: number of samples of output must be equal to the dimensions of plant data "), "iv");
+ error(errmsg);
+ end
+ x = varargin(3)
+ elseif rhs == 2
+ arxModel = arx(plantData,[na nb nk])
+ x = sim(uData,arxModel)
+ end
+ phif = zeros(noOfSample,na+nb)
+ psif = zeros(noOfSample,na+nb)
+ // arranging samples of y matrix
+ for ii = 1:na
+ phif(ii+1:ii+noOfSample,ii) = yData
+ psif(ii+1:ii+noOfSample,ii) = x
+ end
+ // arranging samples of u matrix
+ for ii = 1:nb
+ phif(ii+nk:ii+noOfSample+nk-1,ii+na) = uData
+ psif(ii+nk:ii+noOfSample+nk-1,ii+na) = uData
+ end
+ lhs = psif'*phif
+ lhsinv = pinv(lhs)
+ //pause
+ theta = lhsinv * (psif)' * [yData;zeros(n,1)]
+ ypred = (phif * theta)
+ ypred = ypred(1:size(yData,'r'))
+ e = yData - ypred
+ sigma2 = norm(e)^2/(size(yData,'r') - na - nb)
+ vcov = sigma2 * pinv((phif)' * phif)
+ t = idpoly([1; -theta(1:na)],[zeros(nk,1); theta(na+1:$)],1,1,1,1)
+
+ // estimating the other parameters
+ [temp1,temp2,temp3] = predict(z,t)
+ [temp11,temp22,temp33] = pe(z,t)
+
+ estData = calModelPara(temp1,temp11,na+nb)
+ //pause
+ t.Report.Fit.MSE = estData.MSE
+ t.Report.Fit.FPE = estData.FPE
+ t.Report.Fit.FitPer = estData.FitPer
+ t.Report.Fit.AIC = estData.AIC
+ t.Report.Fit.AICc = estData.AICc
+ t.Report.Fit.nAIC = estData.nAIC
+ t.Report.Fit.BIC = estData.BIC
+ t.TimeUnit = unit
+ //sys = t
+ varargout(1) = t
+endfunction
diff --git a/iv4.sci b/iv4.sci
new file mode 100644
index 0000000..32fd778
--- /dev/null
+++ b/iv4.sci
@@ -0,0 +1,75 @@
+function varargout = iv4(varargin)
+ [lhs, rhs] = argn(0)
+
+ plantData = varargin(1)
+ orderData = varargin(2)
+ na = orderData(1);nb = orderData(2)
+ // arranging na ,nb,nk
+ if size(orderData,"*") == 2 then
+ nk = 1
+ elseif size(orderData,'*') == 3 then
+ nk = orderData(3)
+ end
+ nb1 = nb + nk - 1
+ n = max(na, nb1)
+ // arranging the plant data
+ if typeof(plantData) == 'constant' then
+ Ts = 1;unitData = 'second'
+ elseif typeof(plantData) == 'iddata' then
+ Ts = plantData.Ts;unitData = plantData.TimeUnit
+ plantData = [plantData.OutputData plantData.InputData]
+ end
+ noOfSample = size(plantData,'r')
+ // finding the iv model
+ ivTest = iv(plantData,[na nb nk]);
+ // residual
+ [aTemp,bTemp,cTemp] = pe(plantData,ivTest);
+ Lhat = ar(aTemp,na+nb);
+ x = sim(plantData(:,2),ivTest);
+ yData = plantData(:,1);uData = plantData(:,2)
+ Yf = filter(Lhat.a,Lhat.b,[plantData(:,1);zeros(n,1)]);
+ phif = zeros(noOfSample,na+nb)
+ psif = zeros(noOfSample,na+nb)
+ // arranging samples of y matrix
+ for ii = 1:na
+ phif(ii+1:ii+noOfSample,ii) = -yData
+ psif(ii+1:ii+noOfSample,ii) = -x
+ end
+ // arranging samples of u matrix
+ for ii = 1:nb
+ phif(ii+nk:ii+noOfSample+nk-1,ii+na) = uData
+ psif(ii+nk:ii+noOfSample+nk-1,ii+na) = uData
+ end
+ // passing it through the filters
+ for ii = 1:na+nb
+ phif(:,ii) = filter(Lhat.a,Lhat.b,phif(:,ii));
+ psif(:,ii) = filter(Lhat.a,Lhat.b,psif(:,ii));
+ end
+ lhs = psif'*phif
+ lhsinv = pinv(lhs)
+ theta = lhsinv * (psif)' * Yf
+ ypred = (phif * theta)
+ ypred = ypred(1:size(yData,'r'))
+ e = yData - ypred
+ sigma2 = norm(e)^2/(size(yData,'r') - na - nb)
+ vcov = sigma2 * pinv((phif)' * phif)
+ t = idpoly([1; theta(1:na)],[zeros(nk,1); theta(na+1:$)],1,1,1,Ts)
+
+ // estimating the other parameters
+ [temp1,temp2,temp3] = predict(z,t)
+ [temp11,temp22,temp33] = pe(z,t)
+
+ estData = calModelPara(temp1,temp11,na+nb)
+ //pause
+ t.Report.Fit.MSE = estData.MSE
+ t.Report.Fit.FPE = estData.FPE
+ t.Report.Fit.FitPer = estData.FitPer
+ t.Report.Fit.AIC = estData.AIC
+ t.Report.Fit.AICc = estData.AICc
+ t.Report.Fit.nAIC = estData.nAIC
+ t.Report.Fit.BIC = estData.BIC
+ t.TimeUnit = unitData
+ //sys = t
+ varargout(1) = t
+ //varargout(1) = idpoly([1; -theta(1:na)],[zeros(nk,1); theta(na+1:$)],1,1,1,Ts)
+endfunction
diff --git a/misData.mat b/misData.mat
new file mode 100644
index 0000000..6d26b91
--- /dev/null
+++ b/misData.mat
Binary files differ
diff --git a/misdata.sci b/misdata.sci
new file mode 100644
index 0000000..c1ced50
--- /dev/null
+++ b/misdata.sci
@@ -0,0 +1,41 @@
+function varargout = misdata(varargin)
+ [lhs,rhs] = argn(0)
+//------------------------------------------------------------------------------
+// checking the number of inputs
+ if rhs <> 1 then
+ error(msprintf(gettext("%s:Wrong number of input arguments.\n"),"misdata"))
+ end
+//------------------------------------------------------------------------------
+ ioData = varargin(1)
+ if typeof(ioData) <> "iddata" then
+ error(msprintf(gettext("%s: Plant input data must be ""iddata"" type. "),"misdata"))
+ end
+ inputMat = ioData.InputData;inputMat = linearINTRP(inputMat,abs(ioData.Ts));ioData.InputData = inputMat;
+ outputMat = ioData.OutputData;outputMat = linearINTRP(outputMat,abs(ioData.Ts));ioData.OutputData = outputMat;
+ varargout(1) = ioData
+endfunction
+
+function varargout = linearINTRP(matData,Ts)
+ // looking for overall nan values
+ nanData = isnan(matData);nanIndex = find(nanData == %T)
+ if ~size(nanIndex,'*') then
+ varargout(1) = matData
+ else
+ tempMat = []
+ matSize = size(matData,'r')
+ // looking for nan in each column
+ for ii = 1:size(matData,'c')
+ nanData = isnan(matData(:,ii));nanIndex = find(nanData == %T);
+ if ~size(nanData,'*') then
+ tempMat = [tempMat matData(,ii)]
+ else
+ timeData = (linspace(1*Ts,matSize*Ts , matSize))';
+ nanMat = isnan(matData(:,ii));
+ data = matData(:,ii)
+ data(nanMat) = interp1(timeData(~nanMat), data(~nanMat), timeData(nanMat));
+ tempMat = [tempMat data]
+ end
+ end
+ varargout(1) = tempMat
+ end
+endfunction
diff --git a/nInputData.sci b/nInputData.sci
new file mode 100644
index 0000000..59a3b64
--- /dev/null
+++ b/nInputData.sci
@@ -0,0 +1,15 @@
+function varargout = nInputSeries(varargin)
+ [lhs rhs] = argn(0)
+ if rhs <> 1 then
+ error(msprintf(gettext("%s: Wrong number of input arguments."),"iddataplot"))
+ end
+ iddataData = varargin(1)
+ if typeof(iddataData) <> 'iddata' then
+ error(msprintf(gettext("%s:Wrong type for input argument %d: ""iddata"" expected."),"nInputSeries",1))
+ end
+ if ~size(iddataData.InputData,'*') then
+ varargout(1) = 1
+ else
+ varargout(1) = size(iddataData.InputData,'c')
+ end
+endfunction
diff --git a/oe.sci b/oe.sci
new file mode 100644
index 0000000..8e249ec
--- /dev/null
+++ b/oe.sci
@@ -0,0 +1,118 @@
+// Copyright (C) 2015 - IIT Bombay - FOSSEE
+//
+// 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
+// Authors: Harpreet, Ashutosh
+// Organization: FOSSEE, IIT Bombay
+
+function sys = oe(varargin)
+
+// Estimates Discrete time BJ model
+// y(t) = [B(q)/F(q)]u(t) + [C(q)/D(q)]e(t)
+// Current version uses random initial guess
+// Need to get appropriate guess from OE and noise models
+
+ [lhs , rhs] = argn();
+ if ( rhs < 2 ) then
+ errmsg = msprintf(gettext("%s: Unexpected number of input arguments : %d provided while should be 2"), "oe_2", rhs);
+ error(errmsg)
+ end
+
+
+ z = varargin(1)
+ if typeof(z) == 'iddata' then
+ Ts = z.Ts;unit = z.TimeUnit
+ z = [z.OutputData z.InputData]
+ elseif typeof(z) == 'constant' then
+ Ts = 1;unit = 'seconds'
+ end
+ if ((~size(z,2)==2) & (~size(z,1)==2)) then
+ errmsg = msprintf(gettext("%s: input and output data matrix should be of size (number of data)*2"), "oe_2");
+ error(errmsg);
+ end
+
+ if (~isreal(z)) then
+ errmsg = msprintf(gettext("%s: input and output data matrix should be a real matrix"), "oe_2");
+ error(errmsg);
+ end
+//
+ n = varargin(2)
+ if (size(n,"*")<2| size(n,"*")>3) then
+ errmsg = msprintf(gettext("%s: The order and delay matrix [nb nf nk] should be of size [2 4]"), "oe_2");
+ error(errmsg);
+ end
+
+ if (size(find(n<0),"*") | size(find(((n-floor(n))<%eps)== %f))) then
+ errmsg = msprintf(gettext("%s: values of order and delay matrix [nb nf nk] should be nonnegative integer number "), "oe_2");
+ error(errmsg);
+ end
+//
+ nb= n(1); nf = n(2);
+//
+ if (size(n,"*") == 2) then
+ nk = 1
+ else
+ nk = n(3);
+ end
+
+ // storing U(k) , y(k) and n data in UDATA,YDATA and NDATA respectively
+ YDATA = z(:,1);
+ UDATA = z(:,2);
+ NDATA = size(UDATA,"*");
+ function e = G(p,m)
+ e = YDATA - _objoefun(UDATA,YDATA,p,nf,nb,nk);
+ endfunction
+ tempSum = nf+nb
+ p0 = linspace(0.04,0.041,tempSum)';
+ [var,errl] = lsqrsolve(p0,G,size(UDATA,"*"));
+ err = (norm(errl)^2);
+ opt_err = err;
+ resid = G(var,[]);
+ f = poly([1; var(nb+1:nb+nf)],"q","coeff");
+ b = poly([repmat(0,nk,1);var(1:nb)]',"q","coeff");
+ t = idpoly(1,coeff(b),1,1,coeff(f),Ts)
+
+ // estimating the other parameters
+ [temp1,temp2,temp3] = predict(z,t)
+ [temp11,temp22,temp33] = pe(z,t)
+
+ estData = calModelPara(temp1,temp11,n(1)+n(2))
+ //pause
+ t.Report.Fit.MSE = estData.MSE
+ t.Report.Fit.FPE = estData.FPE
+ t.Report.Fit.FitPer = estData.FitPer
+ t.Report.Fit.AIC = estData.AIC
+ t.Report.Fit.AICc = estData.AICc
+ t.Report.Fit.nAIC = estData.nAIC
+ t.Report.Fit.BIC = estData.BIC
+ t.TimeUnit = unit
+ sys = t
+ //sys = t
+ //sys.TimeUnit = unit
+endfunction
+
+function yhat = _objoefun(UDATA,YDATA,x,nf,nb,nk)
+ x=x(:)
+ q = poly(0,'q')
+ tempSum = nb+nf
+ // making polynomials
+ b = poly([repmat(0,nk,1);x(1:nb)]',"q","coeff");
+ f = -1*poly([x(nb+1:nb+nf)]',"q","coeff")
+ fSize = coeff(f);bSize = coeff(b)
+ maxDelay = max([length(fSize) length(bSize)])
+ yhat = [YDATA(1:maxDelay)]
+ for k=maxDelay+1:size(UDATA,"*")
+ tempB = 0
+ for ii = 1:size(bSize,'*')
+ tempB = tempB + bSize(ii)*UDATA(k-ii+1)
+ end
+ tempF = 0
+ for ii = 1:size(fSize,"*")
+ tempF = tempF + fSize(ii)*yhat(k-ii)
+ end
+ yhat = [yhat; [ tempB + tempF ]];
+ end
+endfunction
diff --git a/pe.sci b/pe.sci
new file mode 100644
index 0000000..3a097a8
--- /dev/null
+++ b/pe.sci
@@ -0,0 +1,112 @@
+function varargout = pe(varargin)
+ [lhs,rhs] = argn(0)
+
+
+//------------------------------------------------------------------------------
+ data = varargin(1)
+ model = varargin(2)
+ if rhs == 3 then
+ kStep = varargin(3)
+ elseif rhs == 2 then
+ kStep = 1
+ end
+//------------------------------------------------------------------------------
+// k step analysis
+ if typeof(kStep) <> 'constant' || isnan(kStep) then
+ error(msprintf(gettext("%s:Prediction horizon(k) must be a non-negative integer number or inf.\n"),"pe"))
+ end
+// if given k step is infinity or []
+ if isinf(kStep) || ~size(kStep,'*') then
+ kStep = 1
+ end
+// checking the dimensions
+ if size(kStep,'*') <> 1 || (ceil(kStep)-kStep) then
+ error(msprintf(gettext("%s:Prediction horizon(k) must be a non-negative integer number or inf.\n"),"pe"))
+ end
+//------------------------------------------------------------------------------
+// checking the plant model
+ if typeof(model) ~= 'idpoly' then
+ error(msprintf(gettext("%s:Plant model must be ""idpoly"" type.\n"),"pe"))
+ end
+ modelSampleTime = model.Ts
+ modelTimeUnit = model.TimeUnit
+//------------------------------------------------------------------------------
+//checking the data type
+ if typeof(data) <> 'iddata' && typeof(data) <> 'constant' then
+ error(msprintf(gettext("%s:Sample data must be ""iddata"" type or ""n x 2"" matrix type.\n"),"pe"))
+ end
+// checking the plant data
+ if typeof(data) == 'iddata' then
+ if ~size(data.OutputData,'*') || ~size(data.InputData,'*') then
+ error(msprintf(gettext("%s:Number of sample data in input and output must be equal.\n"),"pe"))
+ end
+ plantSampleTime = data.Ts
+ plantTimeUnit = data.TimeUnit
+ data = [data.OutputData data.InputData]
+ //disp('iddata')
+ elseif typeof(data) == 'constant' then
+ if size(data,'c') ~= 2 then
+ error(msprintf(gettext("%s:Number of sample data in input and output must be equal.\n"),"pe"))
+ end
+ plantSampleTime = model.Ts
+ plantTimeUnit = model.TimeUnit
+ end
+//------------------------------------------------------------------------------
+// comparing the sampling time
+ if modelSampleTime-plantSampleTime <> 0 then
+ error(msprintf(gettext("%s:The sample time of the model and plant data must be equal.\n"),"pe"))
+ end
+// Comparing the time units
+ if ~strcmp(modelTimeUnit,plantTimeUnit) then
+ else
+ error(msprintf(gettext("%s:Time unit of the model and plant data must be equal.\n"),"pe"))
+ end
+//------------------------------------------------------------------------------
+ sampleLength = size(data,'r')
+ // one step ahead prediction
+ [y1 ,x0] = predict(data,model)
+ errorData1 = data(:,1)-y1
+ if kStep == 1 then
+ eCapData = errorData1
+ elseif kStep > 1 then
+ aPoly = poly(model.a,'q','coeff');
+ dPoly = poly(model.d,'q','coeff');
+ adCoeff = coeff(aPoly*dPoly);adCoeff = adCoeff(length(adCoeff):-1:1);
+ adPoly = poly(adCoeff,'q','coeff')
+ cCoeff = model.c;cCoeff = cCoeff(length(cCoeff):-1:1);
+ cPoly = poly(cCoeff,'q','coeff')
+ hBar = clean((ldiv(cPoly,adPoly,kStep))',0.00001)
+ hBar = hBar(length(hBar):-1:1)
+ hBarLength = length(hBar)
+ errorData1 = [zeros(hBarLength,1);errorData1]
+ eCapData = []
+
+ for ii = 1:sampleLength+1
+ eCapData = [eCapData; hBar*errorData1(ii:ii+hBarLength-1)]
+ end
+ end
+
+ timeData = (modelSampleTime:modelSampleTime:(sampleLength)*modelSampleTime)'
+ pseudoData = size(eCapData,'r')
+ eCapData = eCapData(abs(pseudoData-sampleLength)+1:$)
+ //pause
+ if lhs == 1 then
+ clf()
+ plot(timeData,eCapData)
+ axisData = gca()
+
+ tempTimeUnit = 'Time('+modelTimeUnit+')'
+ xtitle('Predicted Response',tempTimeUnit,'y')
+ xgrid
+ //pause
+ varargout(1) = 0
+
+ elseif lhs == 2 then
+ varargout(1) = eCapData
+ varargout(2) = 0
+ elseif lhs == 3 then
+ varargout(1) = eCapData
+ varargout(2) = timeData
+ varargout(3) = 0
+ end
+endfunction
diff --git a/predict.sci b/predict.sci
new file mode 100644
index 0000000..fcd4b23
--- /dev/null
+++ b/predict.sci
@@ -0,0 +1,192 @@
+//[y,x0] = predict(data,idpoly,k)
+//References
+//Digital Control(11.1.2) by Kanna M.Moudgalya
+//System Identification Theory for User Second Edition (3.2) by Lennart Ljung
+// Code Aurthor - Ashutosh Kumar Bhargava
+
+function varargout = predict(varargin)
+ [lhs,rhs] = argn(0)
+//------------------------------------------------------------------------------
+// checking the number of inputs
+ if rhs < 2 || rhs > 3 then
+ error(msprintf(gettext("%s:Wrong number of input arguments.\n"),"predict"))
+ end
+//------------------------------------------------------------------------------
+ data = varargin(1)
+ model = varargin(2)
+ if rhs == 3 then
+ kStep = varargin(3)
+ elseif rhs == 2 then
+ kStep = 1
+ end
+
+//------------------------------------------------------------------------------
+
+// k step analysis
+ if typeof(kStep) <> 'constant' || isnan(kStep) then
+ error(msprintf(gettext("%s:Prediction horizon(k) must be a non-negative integer number or inf.\n"),"predict"))
+ end
+// if given k step is infinity or []
+ if isinf(kStep) || ~size(kStep,'*') then
+ kStep = 1
+ end
+// checking the dimensions
+ if size(kStep,'*') <> 1 || (ceil(kStep)-kStep) then
+ error(msprintf(gettext("%s:Prediction horizon(k) must be a non-negative integer number or inf.\n"),"predict"))
+ end
+//------------------------------------------------------------------------------
+// checking the plant model
+ if typeof(model) ~= 'idpoly' then
+ error(msprintf(gettext("%s:Plant model must be ""idpoly"" type.\n"),"predict"))
+ end
+ modelSampleTime = model.Ts
+ modelTimeUnit = model.TimeUnit
+//------------------------------------------------------------------------------
+//checking the data type
+ if typeof(data) <> 'iddata' && typeof(data) <> 'constant' then
+ error(msprintf(gettext("%s:Sample data must be ""iddata"" type or ""n x 2"" matrix type.\n"),"predict"))
+ end
+// checking the plant data
+ if typeof(data) == 'iddata' then
+ if ~size(data.OutputData,'*') || ~size(data.InputData,'*') then
+ error(msprintf(gettext("%s:Number of sample data in input and output must be equal.\n"),"predict"))
+ end
+ plantSampleTime = data.Ts
+ plantTimeUnit = data.TimeUnit
+ data = [data.OutputData data.InputData]
+ //disp('iddata')
+ elseif typeof(data) == 'constant' then
+ if size(data,'c') ~= 2 then
+ error(msprintf(gettext("%s:Number of sample data in input and output must be equal.\n"),"predict"))
+ end
+ plantSampleTime = model.Ts
+ plantTimeUnit = model.TimeUnit
+ end
+//------------------------------------------------------------------------------
+// comparing the sampling time
+ if modelSampleTime-plantSampleTime <> 0 then
+ error(msprintf(gettext("%s:The sample time of the model and plant data must be equal.\n"),"predict"))
+ end
+// Comparing the time units
+ if ~strcmp(modelTimeUnit,plantTimeUnit) then
+ else
+ error(msprintf(gettext("%s:Time unit of the model and plant data must be equal.\n"),"predict"))
+ end
+//------------------------------------------------------------------------------
+// ckecking the k step size. if it greater than number of sample size then the
+// k step will become 1
+ if kStep >= size(data,'r') then
+ kStep = 1
+ end
+//------------------------------------------------------------------------------
+ //storing the plant data
+ // B(z) C(z)
+ // y(n) = ---------- u(n) + ---------- e(n)
+ // A(z)*F(z) A(z)*D(z)
+ aPoly = poly(model.a,'q','coeff');
+ bPoly = poly(model.b,'q','coeff');
+ cPoly = poly(model.c,'q','coeff');
+ dPoly = poly(model.d,'q','coeff');
+ fPoly = poly(model.f,'q','coeff');
+ Gq = bPoly/(aPoly*fPoly)
+ Hq = cPoly/(aPoly*dPoly)
+
+ if kStep == 1 then
+ Wkq = Hq^-1
+ elseif kStep > 1 then
+ adCoeff = coeff(aPoly*dPoly);adCoeff = adCoeff(length(adCoeff):-1:1);
+ adPoly = poly(adCoeff,'q','coeff')
+ cCoeff = model.c;cCoeff = cCoeff(length(cCoeff):-1:1);
+ cPoly = poly(cCoeff,'q','coeff')
+ hBar = clean((ldiv(cPoly,adPoly,kStep))',0.00001)
+ hBarPoly = poly(hBar,'q','coeff')
+ Wkq = hBarPoly*Hq^-1
+ end
+ WkqGq = Wkq * Gq
+ tempWkqGq = coeff(WkqGq.den)
+ if tempWkqGq(1) <> 1 then
+ WkqGq.num = WkqGq.num/tempWkqGq(1)
+ WkqGq.den = WkqGq.den/tempWkqGq(1)
+ end
+ Wkq1 = 1-Wkq
+ tempWkq1 = coeff(Wkq1.den)
+ if tempWkq1(1) == 1 then
+ Wkq1.num = Wkq1.num/tempWkq1(1)
+ Wkq1.den = Wkq1.den/tempWkq1(1)
+ end
+ //pause
+//------------------------------------------------------------------------------
+ // storing the plant data
+ uCoeff = coeff(WkqGq.num*Wkq1.den)
+ yCoeff = coeff(WkqGq.den*Wkq1.num)
+ yCapCoeff = coeff(WkqGq.den*Wkq1.den)
+ //pause
+ lengthuCoeff = length(uCoeff)
+ lengthyCoeff = length(yCoeff)
+ lengthyCapCoeff = length(yCapCoeff)
+//------------------------------------------------------------------------------
+ // keeping initial conditions equal to zero
+ uData = zeros(lengthuCoeff,1)
+ yData = zeros(lengthyCoeff,1)
+ yCapData = zeros(lengthyCapCoeff-1,1)
+ uData = [uData;data(:,2)]
+ yData = [yData;data(:,1)]
+ sampleData = size(data,'r')
+ //pause
+ // reversing the coefficients
+ if ~size(uCoeff,'*') then
+ uCoeff = 0
+ else
+ uCoeff = uCoeff(lengthuCoeff:-1:1)
+ end
+ if ~size(yCoeff) then
+ yCoeff = 0
+ else
+ yCoeff = yCoeff(lengthyCoeff:-1:1)
+ end
+ if ~size(yCapCoeff,'*') then
+ yCapCoeff = 0
+ else
+ yCapCoeff = -yCapCoeff(lengthyCapCoeff:-1:2)
+ end
+ //pause
+ for ii = 1:sampleData+1
+ //pause
+ if ~size(uData(ii:ii+lengthuCoeff-1),'*') then
+ tempu = 0
+ else
+ tempu = uCoeff*uData(ii:ii+lengthuCoeff-1);
+ end
+ if ~size(yData(ii:ii+lengthyCoeff-1),'*')
+ tempy = 0
+ else
+ tempy = yCoeff*yData(ii:ii+lengthyCoeff-1);
+ end
+ if ~size(yCapData(ii:ii+lengthyCapCoeff-2),'*') then
+ tempyCap = 0
+ else
+ tempyCap = yCapCoeff*yCapData(ii:ii+lengthyCapCoeff-2);
+ end
+ yCapData = [yCapData;tempu+tempy+tempyCap];
+ end
+ // pause
+ extraSample = abs(size(yCapData,'r')-sampleData)
+ yCapData = yCapData(extraSample+1:$)
+ timeData = ((modelSampleTime:modelSampleTime:(size(yCapData,'r')*modelSampleTime))');
+ if lhs == 1 then
+ clf()
+ plot(timeData,yCapData)
+ axisData = gca()
+ tempTimeUnit = 'Time('+modelTimeUnit+')'
+ xtitle('Predicted Response',tempTimeUnit,'y')
+ xgrid
+ varargout(1) = 0
+ elseif lhs == 2 then
+ varargout(1) = yCapData
+ varargout(2) = 0
+ elseif lhs == 3 then
+ varargout(1) = yCapData
+ varargout(2) = timeData
+ varargout(3) = 0
+ end
+endfunction
diff --git a/rARX.mat b/rARX.mat
new file mode 100644
index 0000000..22171f2
--- /dev/null
+++ b/rARX.mat
Binary files differ
diff --git a/rData.mat b/rData.mat
new file mode 100644
index 0000000..605635f
--- /dev/null
+++ b/rData.mat
Binary files differ
diff --git a/rarx.sci b/rarx.sci
new file mode 100644
index 0000000..e2c291c
--- /dev/null
+++ b/rarx.sci
@@ -0,0 +1,63 @@
+function varargout = rarx(varargin)
+ [lhs,rhs] = argn(0)
+//
+ plantData = varargin(1)
+ orderData = varargin(2)
+ na = orderData(1);nb = orderData(2)
+ // arranging na ,nb,nk
+ if size(orderData,"*") == 2 then
+ nk = 1
+ elseif size(orderData,'*') == 3 then
+ nk = orderData(3)
+ end
+ // storing the lambda value
+ if rhs == 3 then
+ lambda = varargin(3)
+ else
+ lambda = 0.95
+ end
+
+ nb1 = nb + nk - 1
+ n = max(na, nb1)
+ // arranging the plant data
+ if typeof(plantData) == 'constant' then
+ Ts = 1;unitData = 'second'
+ elseif typeof(plantData) == 'iddata' then
+ Ts = plantData.Ts;unitData = plantData.TimeUnit
+ plantData = [plantData.OutputData plantData.InputData]
+ end
+ N = size(plantData,'r')
+ uIndex = nk:nb1
+ yIndex = []
+ if na <> 0 then
+ yIndex = 1:na
+ end
+ df = N - na - nb
+ Plast = 10^4 * (eye(na+nb,na + nb))
+ theta = zeros(N + 1, na + nb)
+ yHat = plantData(:,1);yData = plantData(:,1)
+ tempData = zeros(N,na+nb)
+ for ii = 1:na
+ tempData(ii+1:ii+N,ii) = -plantData(:,1)
+ end
+ // arranging samples of u matrix
+ for ii = 1:nb
+ tempData(ii+nk:ii+N+nk-1,ii+na) = plantData(:,2)
+ end
+ //tempData = [zeros(1,na+nb);tempData]
+ tempData = tempData(1:N+1,:)
+
+ for ii = 1:N
+ temp = tempData(ii,:)
+ yHat(ii) = temp*theta(ii,:)'
+ eps_i = yData(ii)-yHat(ii)
+ kappa_i = Plast * temp'/(lambda + temp * Plast * temp')
+ theta(ii+1,:) = ((theta(ii,:))' + eps_i * kappa_i)'
+ Plast = (eye(na + nb,na + nb) - kappa_i * (temp)) * Plast/lambda
+
+ end
+ theta = theta(1:N,:)
+ yHat = yHat(1:N)
+
+ varargout(1) = struct('theta',theta,'yHat',yHat)
+endfunction
diff --git a/sim.sci b/sim.sci
new file mode 100644
index 0000000..9cd03f1
--- /dev/null
+++ b/sim.sci
@@ -0,0 +1,65 @@
+function varargout = sim(varargin)
+ [lhs,rhs] = argn(0)
+ //checking the number of inputs
+ if rhs < 2 || rhs > 3 then
+ error(msprintf(gettext("%s: Unexpected number of input arguments "),"sim"))
+ end
+ modelData = varargin(2)
+ inputData = varargin(1)
+ //checking the first input
+ if typeof(inputData) <> "constant" && typeof(inputData) <> "iddata" then
+ error(msprintf(gettext("%s: Plant input data must be ""iddata"" type or ""double vector"". "),"sim"))
+ end
+ if typeof(inputData) == "iddata" then
+ inputData = inputData.InputData
+ end
+ if ~iscolumn(inputData) then
+ error(msprintf(gettext("%s: Plant input data must be ""double column vector"". "),"sim"))
+ end
+ //checking the plant model type
+ if typeof(modelData) <> "idpoly" then
+ error(msprintf(gettext("%s: Plant model must be ""idpoly"" type. "),"sim"))
+ end
+ //adding noise
+ if rhs == 3 then
+ noiseFlag = varargin(3)
+ if typeof(noiseFlag) <> 'boolean' then
+ error(msprintf(gettext("%s: Last input data must be ""boolean"".type "),"sim"))
+ end
+ else
+ noiseFlag = %F
+ end
+ //adding noise of zero mean and 1 as standard deviation
+ if noiseFlag then
+ numberOfSamples = size(inputData,'r')
+ R = = grand(numberOfSamples,1,"nor",0,1)
+ inputData = inputData+R
+ end
+ z = poly(0,'z')
+ aPoly = poly(modelData.a(length(modelData.a):-1:1),'z','coeff')
+ bPoly = poly(modelData.b,'z','coeff')
+ fPoly = poly(modelData.f(length(modelData.f):-1:1),'z','coeff')
+ afPoly = aPoly*fPoly
+
+ bCoeff = modelData.b
+ extra = 1
+ if ~bCoeff(1,1) then
+ afCoeff = coeff(afPoly)
+ bLength = length(bCoeff);afLength = length(afCoeff)
+ if bLength == afLength then
+ extra = 1
+ else
+ extra = z^-(bLength-afLength)
+ end
+ end
+
+ bPoly = poly(modelData.b(length(modelData.b):-1:1),'z','coeff')
+ sys = syslin('d',bPoly,afPoly)*extra
+ outputData = (flts(inputData',sys))'
+
+ if typeof(varargin(1)) == "iddata" then
+ outputData = iddata(outputData,[],varargin(1).Ts)
+ outputData.TimeUnit = varargin(1).TimeUnit
+ end
+ varargout(1) = outputData
+endfunction
diff --git a/spa.sci b/spa.sci
new file mode 100644
index 0000000..4b91cc9
--- /dev/null
+++ b/spa.sci
@@ -0,0 +1,73 @@
+function varargout = spa(varargin)
+ [lhs , rhs] = argn();
+ if ( rhs < 1 || rhs > 3 ) then
+ errmsg = msprintf(gettext("%s: Wrong number of input arguments" ),"spa");
+ error(errmsg)
+ elseif typeof(varargin(1)) <> "iddata" then
+ error(msprintf(gettext("%s:Plant data must be ""iddata"" type.\n"),"spa"))
+ end
+
+ plantData = varargin(1)
+ windowSize = %F
+ inputFreq = %F
+//------------------------------------------------------------------------------
+// arranging the plant data
+ inputData = plantData.InputData;
+ if ~size(inputData,"*") then
+ error(msprintf(gettext("%s:Input data must be non-zero vector. \n"),"spa"))
+ end
+ outputData = plantData.OutputData
+ if ~size(outputData,"*") then
+ error(msprintf(gettext("%s:Output data must be non-zero vector. \n"),"spa"))
+ end
+//------------------------------------------------------------------------------
+ N = size(inputData,'r')
+ nout = size(outputData,'c')
+ nin = size(inputData,'c')
+
+ if ~windowSize then
+ windowSize = min(N/10, 30)
+ end
+
+ if inputFreq then
+ else
+ inputFreq = (1:128)/128 * %pi/plantData.Ts
+ end
+
+ M = windowSize
+ Ryu = xcov(outputData,inputData,M,'biased')
+ Ruu = xcov(inputData,inputData,M,'biased')
+ Ryy = xcov(outputData,outputData,M,'biased')
+ G = [];spect = [];
+ for ii = 1:nout
+ phiY = spaCalculation(inputFreq,Ryy,M)
+ temp = phiY
+ for jj = 1:nin
+ phiYU = spaCalculation(inputFreq,Ryu,M)//sapply(freq, cov2spec, Ryu[i, j, ], M)
+ phiU = spaCalculation(inputFreq,Ruu,M)
+ G = [G phiYU./phiU]
+ //pause
+ temp = temp - phiYU .* conj(phiYU)./phiU
+ end
+ spect = [spect temp]
+ end
+// disp(size(spect))
+// disp(spect)
+ frdData = frd(G',inputFreq',plantData.Ts,spect')
+ varargout(1) = frdData
+endfunction
+
+function varargout = spaCalculation(varargin)
+ inputFreq = varargin(1)
+ xcovData = varargin(2)
+ Mnumber = varargin(3)
+ temp = []
+ win_l=window('hn',2*Mnumber+1)
+ for ii = 1:size(inputFreq,'c')
+ seq1 = exp(-(%i) * (-Mnumber:Mnumber) * inputFreq(ii))
+ data = (seq1.*win_l)
+ data2 = sum(data.*xcovData')
+ temp = [temp data2]
+ end
+ varargout(1) = temp
+endfunction
diff --git a/stem.sci b/stem.sci
new file mode 100644
index 0000000..2bf5b0a
--- /dev/null
+++ b/stem.sci
@@ -0,0 +1,19 @@
+// 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/stepest.sci b/stepest.sci
new file mode 100644
index 0000000..6310971
--- /dev/null
+++ b/stepest.sci
@@ -0,0 +1,71 @@
+function varargout = stepest(varargin)
+ [lhs,rhs] = argn(0)
+ //checking the number of inputs
+ if rhs > 2 then
+ error(msprintf(gettext("%s: Unexpected number of input arguments "),"stepest"))
+ end
+ modelData = varargin(1)
+ if typeof(modelData) <> "idpoly" then
+ error(msprintf(gettext("%s: Plant model must be ""idpoly"" type. "),"stepest"))
+ end
+ //adding noise
+ if rhs == 2 then
+ noiseFlag = varargin(2)
+ if typeof(noiseFlag) <> 'boolean' then
+ error(msprintf(gettext("%s: Last input data must be ""boolean"".type "),"stepest"))
+ end
+ else
+ noiseFlag = %F
+ end
+ z = poly(0,'z')
+ aPoly = poly(modelData.a(length(modelData.a):-1:1),'z','coeff')
+ bPoly = poly(modelData.b,'z','coeff')
+ fPoly = poly(modelData.f(length(modelData.f):-1:1),'z','coeff')
+ afPoly = aPoly*fPoly
+
+ bCoeff = modelData.b
+ extra = 1
+ if ~bCoeff(1,1) then
+ afCoeff = coeff(afPoly)
+ bLength = length(bCoeff);afLength = length(afCoeff)
+ if bLength == afLength then
+ extra = 1
+ else
+ extra = z^-(bLength-afLength)
+ end
+ end
+
+ bPoly = poly(modelData.b(length(modelData.b):-1:1),'z','coeff')
+ sys = syslin('d',bPoly,afPoly)*extra
+ if size(find(gsort(abs(roots(sys.den)))>1),'*') then
+ tempRoot = clean(roots(sys.den))
+ [sorted index] = gsort(abs(real(tempRoot)))
+ tempRoot = tempRoot(index)
+ n = 26/log10(abs(real(tempRoot(1))))
+
+ else
+ finalValue = horner(sys,1)*1
+ n = 10
+ tempValue = sum(ldiv(sys.num,sys.den,n))
+ if finalValue >= 0 then
+ while finalValue*0.99 >= tempValue || n >=1000
+ n = n+10
+ tempValue = sum(ldiv(sys.num,sys.den,n))
+ end
+ elseif finalValue < 0 then
+ while finalValue*0.99 <tempValue || n >=1000
+ n = n+10
+ tempValue = sum(ldiv(sys.num,sys.den,n))
+ end
+ end
+ end
+ uData = ones(1,n)
+ yData = flts(uData,sys)
+ timeData = (0:(n-1))*modelData.Ts
+ if noiseFlag then
+ varargout(1) = yData'
+ else
+ plot2d2(timeData,yData,2)
+ varargout(1) = []
+ end
+endfunction
generated by cgit (git 2.34.1) at 2024-12-19 09:16:47 +0000